LB stokes sphere test case update

src/core/grid_based_algorithms/lbboundaries.cpp

@@ -352,11 +352,14 @@ int lbboundary_get_force(void *lbb, double *f) {
     mpi_gather_stats(8, forces.data(), nullptr, nullptr, nullptr);
 
     f[0] = forces[3 * no + 0] * lbpar.agrid /
-           lbpar.tau; // lbpar.tau; TODO this makes the units wrong and
+           (lbpar.tau * lbpar.tau *
+            lbpar.rho); // lbpar.tau; TODO this makes the units wrong and
     f[1] = forces[3 * no + 1] * lbpar.agrid /
-           lbpar.tau; // lbpar.tau; the result correct. But it's 3.13AM
+           (lbpar.tau * lbpar.tau *
+            lbpar.rho); // lbpar.tau; the result correct. But it's 3.13AM
     f[2] = forces[3 * no + 2] * lbpar.agrid /
-           lbpar.tau; // lbpar.tau; on a Saturday at the ICP. Someone fix.
+           (lbpar.tau * lbpar.tau *
+            lbpar.rho); // lbpar.tau; on a Saturday at the ICP. Someone fix.
 #else
     return ES_ERROR;
 #endif

src/script_interface/lbboundaries/LBBoundary.hpp

@@ -21,6 +21,7 @@ along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
 #include "ScriptInterface.hpp"
 #include "auto_parameters/AutoParameters.hpp"
+#include "core/communication.hpp"
 #include "core/grid_based_algorithms/lbboundaries/LBBoundary.hpp"
 #include "shapes/Shape.hpp"
 
@@ -65,7 +66,11 @@ class LBBoundary : public AutoParameters<LBBoundary> {
 
   Variant call_method(const std::string &name, const VariantMap &) override {
     if (name == "get_force") {
-      return m_lbboundary->get_force();
+      // The get force method uses mpi callbacks on lb cpu
+      if (this_node == 0)
+        return m_lbboundary->get_force();
+      else
+        return none;
     }
     return none;
   }

testsuite/python/CMakeLists.txt

@@ -93,7 +93,7 @@ python_test(FILE reaction_ensemble.py MAX_NUM_PROC 4)
 python_test(FILE constant_pH.py MAX_NUM_PROC 4)
 python_test(FILE swimmer_reaction.py MAX_NUM_PROC 1)
 python_test(FILE writevtf.py MAX_NUM_PROC 4)
-python_test(FILE lb_stokes_sphere_gpu.py MAX_NUM_PROC 1)
+python_test(FILE lb_stokes_sphere.py MAX_NUM_PROC 1)
 python_test(FILE ek_eof_one_species_x.py MAX_NUM_PROC 1)
 python_test(FILE ek_eof_one_species_y.py MAX_NUM_PROC 1)
 python_test(FILE ek_eof_one_species_z.py MAX_NUM_PROC 1)

testsuite/python/lb_stokes_sphere_gpu.py → testsuite/python/lb_stokes_sphere.py

@@ -32,39 +32,43 @@
 import numpy as np
 import sys
 
-
-@ut.skipIf(not has_features(["LB_GPU", "LB_BOUNDARIES_GPU"]),
-           "Features not available, skipping test!")
-class Stokes(ut.TestCase):
+# Define the LB Parameters
+TIME_STEP = 0.4
+AGRID = 1.0
+KVISC = 5.0
+DENS = 1.0
+LB_PARAMS = {'agrid': AGRID,
+             'dens': DENS,
+             'visc': KVISC,
+             'fric': 1.0,
+             'tau': TIME_STEP}
+# System setup
+radius = 5.4
+box_width = 44
+real_width = box_width + 2 * AGRID
+box_length = 50
+v = [0, 0, 0.01]  # The boundary slip
+
+
+class Stokes(object):
+    lbf = None
+    system = espressomd.System(box_l=[real_width, real_width, box_length])
+    system.box_l = [real_width, real_width, box_length]
+    system.time_step = TIME_STEP
+    system.cell_system.skin = 0.4
+
+    # The temperature is zero.
+    system.thermostat.set_lb(kT=0)
 
     def test_stokes(self):
-        # System setup
-        agrid = 1
-        radius = 5.5
-        box_width = 54
-        real_width = box_width + 2 * agrid
-        box_length = 54
-        system = espressomd.System(box_l=[real_width, real_width, box_length])
-        system.box_l = [real_width, real_width, box_length]
-        system.time_step = 0.4
-        system.cell_system.skin = 0.4
-
-        # The temperature is zero.
-        system.thermostat.set_lb(kT=0)
-
-        # LB Parameters
-        v = [0, 0, 0.01]  # The boundary slip
-        kinematic_visc = 5.0
-
-        # Invoke LB fluid
-        lbf = lb.LBFluidGPU(visc=kinematic_visc, dens=1,
-                            agrid=agrid, tau=system.time_step, fric=1)
-        system.actors.add(lbf)
+        self.system.actors.clear()
+        self.system.lbboundaries.clear()
+        self.system.actors.add(self.lbf)
 
         # Setup walls
         walls = [None] * 4
         walls[0] = lbboundaries.LBBoundary(shape=shapes.Wall(
-            normal=[-1, 0, 0], dist=-(1 + box_width)), velocity=v)
+                                           normal=[-1, 0, 0], dist=-(1 + box_width)), velocity=v)
         walls[1] = lbboundaries.LBBoundary(
             shape=shapes.Wall(
                 normal=[
@@ -85,32 +89,52 @@ def test_stokes(self):
             velocity=v)
 
         for wall in walls:
-            system.lbboundaries.add(wall)
+            self.system.lbboundaries.add(wall)
 
         # setup sphere without slip in the middle
         sphere = lbboundaries.LBBoundary(shape=shapes.Sphere(
             radius=radius, center=[real_width / 2] * 2 + [box_length / 2], direction=1))
 
-        system.lbboundaries.add(sphere)
+        self.system.lbboundaries.add(sphere)
 
         def size(vector):
             tmp = 0
             for k in vector:
                 tmp += k * k
             return np.sqrt(tmp)
 
-        system.integrator.run(800)
+        self.system.integrator.run(800)
 
-        stokes_force = 6 * np.pi * kinematic_visc * radius * size(v)
+        stokes_force = 6 * np.pi * KVISC * radius * size(v)
         print("Stokes' Law says: f=%f" % stokes_force)
-        
+
         # get force that is exerted on the sphere
-        for i in range(5):
-            system.integrator.run(200)
+        for i in range(4):
+            self.system.integrator.run(200)
             force = sphere.get_force()
             print("Measured force: f=%f" % size(force))
             self.assertLess(abs(1.0 - size(force) / stokes_force), 0.06)
 
+##Invoke the GPU LB
+
+
+@ut.skipIf(not espressomd.has_features(
+    ['LB_GPU', 'LB_BOUNDARIES_GPU', 'EXTERNAL_FORCES']), "Skipping test due to missing features.")
+class LBGPUStokes(ut.TestCase, Stokes):
+
+    def setUp(self):
+        self.lbf = espressomd.lb.LBFluidGPU(**LB_PARAMS)
+
+#Invoke the CPU LB
+
+
+@ut.skipIf(not espressomd.has_features(
+    ['LB', 'LB_BOUNDARIES', 'EXTERNAL_FORCES']), "Skipping test due to missing features.")
+class LBCPUStokes(ut.TestCase, Stokes):
+
+    def setUp(self):
+        self.lbf = espressomd.lb.LBFluid(**LB_PARAMS)
+
 
 if __name__ == "__main__":
     ut.main()