walberla: Restore all LB samples

samples/object_in_fluid/motivation.py

@@ -26,6 +26,7 @@
 espressomd.assert_features(required_features)
 
 import os
+import tqdm
 import argparse
 import warnings
 
@@ -46,8 +47,8 @@
 
 boxX = 22.0
 boxY = 14.0
-boxZ = 15.0
-time_step = 0.1
+boxZ = 6.0
+time_step = 0.05
 
 system = espressomd.System(box_l=(boxX, boxY, boxZ))
 system.time_step = time_step
@@ -63,18 +64,15 @@
 # creating the RBCs
 cell0 = oif.OifCell(cell_type=cell_type,
                     particle_type=0, origin=[5.0, 5.0, 3.0])
-cell1 = oif.OifCell(cell_type=cell_type,
-                    particle_type=1, origin=[5.0, 5.0, 7.0])
 
 # cell-wall interactions
-system.non_bonded_inter[0, 10].soft_sphere.set_params(
-    a=0.0001, n=1.2, cutoff=0.1, offset=0.0)
-system.non_bonded_inter[1, 10].soft_sphere.set_params(
+system.non_bonded_inter[cell0.particle_type, 10].soft_sphere.set_params(
     a=0.0001, n=1.2, cutoff=0.1, offset=0.0)
 
 # fluid
-lbf = espressomd.lb.LBFluidWalberla(agrid=1., density=1., kinematic_viscosity=1.5,
-                                    tau=0.1, ext_force_density=[0.002, 0., 0.])
+lbf = espressomd.lb.LBFluidWalberla(
+    agrid=1., density=1., kinematic_viscosity=1.5, tau=system.time_step,
+    ext_force_density=[0.025, 0., 0.], single_precision=True)
 system.actors.add(lbf)
 system.thermostat.set_lb(LB_fluid=lbf, gamma=1.5)
 
@@ -140,15 +138,15 @@
     system.constraints.add(shape=shape, particle_type=10)
 
 
-maxCycle = 50
+def write_cells_vtk(i):
+    filepath = os.path.join(output_path, "cell{cell_id}_{index}.vtk")
+    cell0.output_vtk_pos_folded(file_name=filepath.format(cell_id=0, index=i))
+
+
+maxCycle = 100
 # main integration loop
-cell0.output_vtk_pos_folded(file_name=os.path.join(output_path, "cell0_0.vtk"))
-cell1.output_vtk_pos_folded(file_name=os.path.join(output_path, "cell1_0.vtk"))
-for i in range(1, maxCycle):
-    system.integrator.run(steps=500)
-    cell0.output_vtk_pos_folded(
-        file_name=os.path.join(output_path, f"cell0_{i}.vtk"))
-    cell1.output_vtk_pos_folded(
-        file_name=os.path.join(output_path, f"cell1_{i}.vtk"))
-    print(f"time: {i * time_step:.1f}")
+for i in tqdm.tqdm(range(maxCycle)):
+    write_cells_vtk(i)
+    system.integrator.run(steps=100)
+write_cells_vtk(maxCycle)
 print("Simulation completed.")

testsuite/python/ek_diffusion.py

@@ -95,8 +95,7 @@ def detail_test_diffusion(self, single_precision: bool):
         # check that the density is conserved
         np.testing.assert_almost_equal(
             np.sum(simulated_density), self.DENSITY, decimal_precision)
-        if np.any(simulated_density < 0.):
-            self.fail("EK density array contains negative densities!")
+        assert np.all(simulated_density >= 0.), "EK has negative densities"
 
         # check that the maximum is in the right place
         peak = np.unravel_index(
@@ -124,7 +123,7 @@ def detail_test_diffusion(self, single_precision: bool):
             calc_density,
             simulated_density,
             atol=1e-5,
-            rtol=0)
+            rtol=0.)
 
 
 if __name__ == "__main__":

testsuite/python/tests_common.py

@@ -18,6 +18,7 @@
 #
 
 import pathlib
+import itertools
 import numpy as np
 
 
@@ -356,24 +357,16 @@ def fold_index(idx, shape):
 
 
 def get_lb_nodes_around_pos(pos, lbf):
-    """Returns lb node(s) relevant for interpolation around the given position"""
-
-    pos_lb_units = pos / lbf.agrid - .5  # Rel to node centers
+    """Returns LB node(s) relevant for interpolation around the given position"""
 
+    pos_lb_units = pos / lbf.agrid - 0.5  # relative to node centers
     lower_left_index = np.array(np.floor(pos_lb_units), dtype=int)
 
-    # Position exactly on a node? Return just that node
-    if np.all(lower_left_index == pos_lb_units):
-        return [lbf[fold_index(lower_left_index, lbf.shape)]]
-
-    # Otherwise return 8 surrounding nodes
-    res = []
-    for i in (0, 1):
-        for j in (0, 1):
-            for k in (0, 1):
-                res.append(lbf[fold_index(lower_left_index +
-                                          np.array((i, j, k), dtype=int), lbf.shape)])
-    return res
+    nodes = []
+    for i, j, k in itertools.product([0, 1], repeat=3):
+        index = lower_left_index + np.array((i, j, k), dtype=int)
+        nodes.append(lbf[fold_index(index, lbf.shape)])
+    return nodes
 
 
 def random_dipoles(n_particles):
@@ -416,11 +409,11 @@ def check_non_bonded_loop_trace(ut_obj, system, cutoff=None):
             msg=f"Extra pair from core {p}")
         if (p[0], p[1]) in py_distances:
             np.testing.assert_allclose(
-                np.copy(p[4]), -np.copy(py_distances[p[0], p[1]]))
+                np.copy(p[4]), -py_distances[p[0], p[1]])
             del py_distances[p[0], p[1]]
         elif (p[1], p[0]) in py_distances:
             np.testing.assert_allclose(
-                np.copy(p[4]), np.copy(py_distances[p[1], p[0]]))
+                np.copy(p[4]), py_distances[p[1], p[0]])
             del py_distances[p[1], p[0]]
 
     for ids, dist in py_distances.items():

testsuite/scripts/samples/CMakeLists.txt

@@ -99,7 +99,7 @@ sample_test(FILE test_visualization_elc.py)
 sample_test(FILE test_visualization_npt.py)
 sample_test(FILE test_visualization_poiseuille.py)
 sample_test(FILE test_widom_insertion.py)
-# sample_test(FILE test_object_in_fluid__motivation.py) # TODO WALBERLA
+sample_test(FILE test_object_in_fluid__motivation.py)
 sample_test(FILE test_immersed_boundary.py)
 
 add_custom_target(

testsuite/scripts/samples/test_object_in_fluid__motivation.py

@@ -19,6 +19,7 @@
 import importlib_wrapper
 import os
 import pathlib
+import numpy as np
 
 os.chdir("@SAMPLES_DIR@/object_in_fluid")
 sample, skipIfMissingFeatures = importlib_wrapper.configure_and_import(
@@ -38,15 +39,25 @@ def test_file_generation(self):
             "wallBottom.vtk",
             "wallBack.vtk",
             "wallFront.vtk"]
-        for i in [0, 1]:
-            for j in range(sample.maxCycle):
-                basenames.append(f"cell{i}_{j}.vtk")
+        for j in range(sample.maxCycle + 1):
+            basenames.append(f"cell0_{j}.vtk")
 
         # test .vtk files exist
         path_vtk_root = pathlib.Path("output")
         for name in basenames:
             filepath = path_vtk_root / "sim0" / name
-            self.assertTrue(filepath.is_file(), f"File {filepath} not created")
+            self.assertTrue(
+                filepath.is_file(),
+                f"File '{filepath}' not created")
+
+        # make sure we are still in the LB linear regime
+        lb_density = np.copy(sample.lbf[:, :, :].density)
+        np.testing.assert_allclose(lb_density, 1., atol=2e-3)
+
+        # verify cell momentum
+        cell_vel = np.mean(self.system.part.all().v, axis=0)
+        np.testing.assert_allclose(
+            cell_vel, [1.54e-2, 1.8e-3, 0.], rtol=1e-2, atol=1e-6)
 
 
 if __name__ == "__main__":