Fix inertialess tracers

doc/sphinx/particles.rst

@@ -386,6 +386,7 @@ For correct results, the LB thermostat has to be deactivated for virtual sites::
    system.thermostat.set_lb(kT=0, act_on_virtual=False)
 
 Please note that the velocity attribute of the virtual particles does not carry valid information for this virtual sites scheme.
+With the LB GPU implementation, inertialess tracers only work on 1 MPI rank.
 
 .. _Interacting with groups of particles:
 

src/core/virtual_sites/VirtualSitesInertialessTracers.cpp

@@ -30,6 +30,14 @@
 
 #include <algorithm>
 
+static void check_no_vs_exist(char const *const message) {
+  if (std::any_of(cell_structure.local_particles().begin(),
+                  cell_structure.local_particles().end(),
+                  [](Particle const &p) { return p.is_virtual(); })) {
+    runtimeErrorMsg() << "Inertialess Tracers: " << message;
+  }
+}
+
 void VirtualSitesInertialessTracers::after_force_calc() {
   // Now the forces are computed and need to go into the LB fluid
   if (lattice_switch == ActiveLB::CPU) {
@@ -39,16 +47,15 @@ void VirtualSitesInertialessTracers::after_force_calc() {
 #ifdef CUDA
   if (lattice_switch == ActiveLB::GPU) {
     IBM_ForcesIntoFluid_GPU(cell_structure.local_particles(), this_node);
+    if (comm_cart.size() != 1 and this_node != 0) {
+      check_no_vs_exist("The LB GPU method cannot integrate virtual sites when "
+                        "more than 1 MPI ranks are used. The particles on MPI "
+                        "rank >= 2 are now in an undeterminate state.");
+    }
     return;
   }
 #endif
-  if (std::any_of(cell_structure.local_particles().begin(),
-                  cell_structure.local_particles().end(),
-                  [](Particle &p) { return p.is_virtual(); })) {
-    runtimeErrorMsg() << "Inertialess Tracers: No LB method was active but "
-                         "virtual sites present.";
-    return;
-  }
+  check_no_vs_exist("No LB method was active but virtual sites present.");
 }
 
 void VirtualSitesInertialessTracers::after_lb_propagation(double time_step) {

src/core/virtual_sites/lb_inertialess_tracers.cpp

@@ -92,7 +92,7 @@ void IBM_UpdateParticlePositions(ParticleRange const &particles,
   // Euler integrator
   for (auto &p : particles) {
     if (p.is_virtual()) {
-      for (int axis = 0; axis < 2; axis++) {
+      for (int axis = 0; axis < 3; axis++) {
 #ifdef EXTERNAL_FORCES
         if (not p.is_fixed_along(axis))
 #endif

testsuite/python/virtual_sites_tracers.py

@@ -25,13 +25,9 @@
 
 @utx.skipIfMissingFeatures(
     ['VIRTUAL_SITES_INERTIALESS_TRACERS', 'LB_BOUNDARIES'])
-class VirtualSitesTracers(ut.TestCase, VirtualSitesTracersCommon):
+class VirtualSitesTracers(VirtualSitesTracersCommon, ut.TestCase):
 
-    def setUp(self):
-        self.LBClass = espressomd.lb.LBFluid
-
-    def tearDown(self):
-        VirtualSitesTracersCommon.tearDown(self)
+    LBClass = espressomd.lb.LBFluid
 
 
 if __name__ == "__main__":

testsuite/python/virtual_sites_tracers_common.py

@@ -17,6 +17,7 @@
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 #
 import espressomd
+import espressomd.lb
 import espressomd.shapes
 import espressomd.lbboundaries
 import espressomd.virtual_sites
@@ -30,13 +31,18 @@ class VirtualSitesTracersCommon:
     system.time_step = 0.05
     system.cell_system.skin = 0.1
 
+    def setUp(self):
+        self.system.box_l = (self.box_lw, self.box_lw, self.box_height)
+
     def tearDown(self):
-        self.system.part.clear()
-        self.system.actors.clear()
         self.system.thermostat.turn_off()
+        self.system.lbboundaries.clear()
+        self.system.actors.clear()
+        self.system.part.clear()
 
-    def reset_lb(self, ext_force_density=(0, 0, 0)):
+    def reset_lb(self, ext_force_density=(0, 0, 0), dir_walls=2):
         self.system.lbboundaries.clear()
+        self.system.actors.clear()
         self.lbf = self.LBClass(
             kT=0.0, agrid=1, dens=1, visc=1.8,
             tau=self.system.time_step, ext_force_density=ext_force_density)
@@ -47,11 +53,14 @@ def reset_lb(self, ext_force_density=(0, 0, 0)):
             gamma=1)
 
         # Setup boundaries
+        normal = [0, 0, 0]
+        normal[dir_walls] = 1
         walls = [espressomd.lbboundaries.LBBoundary() for k in range(2)]
         walls[0].set_params(shape=espressomd.shapes.Wall(
-            normal=[0, 0, 1], dist=0.5))
+            normal=normal, dist=0.5))
+        normal[dir_walls] = -1
         walls[1].set_params(shape=espressomd.shapes.Wall(
-            normal=[0, 0, -1], dist=-self.box_height - 0.5))
+            normal=normal, dist=-(self.system.box_l[dir_walls] - 0.5)))
 
         for wall in walls:
             self.system.lbboundaries.add(wall)
@@ -70,27 +79,59 @@ def test_aa_method_switching(self):
             self.system.virtual_sites, espressomd.virtual_sites.VirtualSitesInertialessTracers)
 
     def test_advection(self):
-        self.reset_lb(ext_force_density=[0.1, 0, 0])
-        # System setup
-        system = self.system
-
-        system.virtual_sites = espressomd.virtual_sites.VirtualSitesInertialessTracers()
-
-        # Establish steady state flow field
-        p = system.part.add(pos=(0, 5.5, 5.5), virtual=True)
-        system.integrator.run(400)
-
-        p.pos = (0, 5.5, 5.5)
-        system.time = 0
-
-        # Perform integration
-        for _ in range(2):
-            system.integrator.run(100)
-            # compute expected position
-            dist = self.lbf.get_interpolated_velocity(p.pos)[0] * system.time
-            self.assertAlmostEqual(p.pos[0] / dist, 1, delta=0.005)
+        node_grid = self.system.cell_system.node_grid
+        lb_on_gpu = self.LBClass is espressomd.lb.LBFluidGPU
+        for direction in [0, 1, 2]:
+            # System setup
+            system = self.system
+            system.virtual_sites = espressomd.virtual_sites.VirtualSitesInertialessTracers()
+
+            # LB setup with walls
+            ext_force = [0., 0., 0.]
+            ext_force[direction] = 0.1
+            dir_walls = (direction + 2) % 3
+            box_l = 3 * [self.box_lw]
+            box_l[dir_walls] = self.box_height
+            system.box_l = box_l
+            self.reset_lb(ext_force_density=ext_force, dir_walls=dir_walls)
+
+            # Establish steady state flow field
+            system.integrator.run(400)
+
+            # Add tracer in the fluid domain
+            pos_initial = [3.5, 3.5, 3.5]
+            pos_initial[direction] = 0.5
+            p = system.part.add(pos=pos_initial, virtual=True)
+
+            # Perform integration
+            n_steps = 100
+            if node_grid[direction] != 1 and lb_on_gpu:
+                n_steps = 50  # with GPU, tracers must stay on MPI rank 0
+            system.time = 0
+            for _ in range(2):
+                system.integrator.run(n_steps)
+                # compute expected position
+                lb_vel = self.lbf.get_interpolated_velocity(p.pos)
+                ref_dist = lb_vel[direction] * system.time
+                tracer_dist = p.pos[direction] - pos_initial[direction]
+                self.assertAlmostEqual(tracer_dist / ref_dist, 1., delta=0.01)
+
+            self.tearDown()
+
+    def test_zz_exceptions_with_lb(self):
+        node_grid = self.system.cell_system.node_grid
+        lb_on_gpu = self.LBClass is espressomd.lb.LBFluidGPU
+        if lb_on_gpu and sum(node_grid) != 3:
+            self.reset_lb()
+            system = self.system
+            system.virtual_sites = espressomd.virtual_sites.VirtualSitesInertialessTracers()
+            p = system.part.add(pos=(0, 0, 0), virtual=True)
+            system.integrator.run(1)
+            p.pos = (-0.5, -0.5, -0.5)
+            with self.assertRaisesRegex(Exception, "The LB GPU method cannot integrate virtual sites when more than 1 MPI ranks are used"):
+                system.integrator.run(1)
 
-    def test_zz_without_lb(self):
+    def test_zz_exceptions_without_lb(self):
         """Check behaviour without lb. Ignore non-virtual particles, complain on
         virtual ones.
 
@@ -103,5 +144,5 @@ def test_zz_without_lb(self):
         p = system.part.add(pos=(0, 0, 0))
         system.integrator.run(1)
         p.virtual = True
-        with self.assertRaises(Exception):
+        with self.assertRaisesRegex(Exception, "No LB method was active but virtual sites present"):
             system.integrator.run(1)

testsuite/python/virtual_sites_tracers_gpu.py

@@ -26,13 +26,9 @@
 @utx.skipIfMissingGPU()
 @utx.skipIfMissingFeatures(
     ['VIRTUAL_SITES_INERTIALESS_TRACERS', 'LB_BOUNDARIES'])
-class VirtualSitesTracers(ut.TestCase, VirtualSitesTracersCommon):
+class VirtualSitesTracers(VirtualSitesTracersCommon, ut.TestCase):
 
-    def setUp(self):
-        self.LBClass = espressomd.lb.LBFluidGPU
-
-    def tearDown(self):
-        VirtualSitesTracersCommon.tearDown(self)
+    LBClass = espressomd.lb.LBFluidGPU
 
 
 if __name__ == "__main__":