Fix inertialess tracers coupling (#4903)

Fixes #4902

Description of changes:
- couple inertialess tracers to the correct LB grid points

src/core/lb/particle_coupling.cpp

@@ -125,8 +125,9 @@ std::vector<Utils::Vector3d> positions_in_halo(Utils::Vector3d const &pos,
   auto const halo_vec = Utils::Vector3d::broadcast(halo);
   auto const fully_inside_lower = local_box.my_left() + 2. * halo_vec;
   auto const fully_inside_upper = local_box.my_right() - 2. * halo_vec;
-  if (in_box(pos, fully_inside_lower, fully_inside_upper)) {
-    return {pos};
+  auto const pos_folded = box_geo.folded_position(pos);
+  if (in_box(pos_folded, fully_inside_lower, fully_inside_upper)) {
+    return {pos_folded};
   }
   auto const halo_lower_corner = local_box.my_left() - halo_vec;
   auto const halo_upper_corner = local_box.my_right() + halo_vec;
@@ -137,11 +138,11 @@ std::vector<Utils::Vector3d> positions_in_halo(Utils::Vector3d const &pos,
       for (int k : {-1, 0, 1}) {
         Utils::Vector3d shift{{double(i), double(j), double(k)}};
         Utils::Vector3d pos_shifted =
-            pos + Utils::hadamard_product(box_geo.length(), shift);
+            pos_folded + Utils::hadamard_product(box_geo.length(), shift);
 
         if (box_geo.type() == BoxType::LEES_EDWARDS) {
           auto le = box_geo.lees_edwards_bc();
-          auto normal_shift = (pos_shifted - pos)[le.shear_plane_normal];
+          auto normal_shift = (pos_shifted - pos_folded)[le.shear_plane_normal];
           if (normal_shift > std::numeric_limits<double>::epsilon())
             pos_shifted[le.shear_direction] += le.pos_offset;
           if (normal_shift < -std::numeric_limits<double>::epsilon())

src/core/virtual_sites/lb_tracers.cpp

@@ -44,7 +44,6 @@ void lb_tracers_add_particle_force_to_fluid(CellStructure &cell_structure,
     return;
   }
   auto const agrid = lb.get_agrid();
-  auto const to_lb_units = 1. / agrid;
 
   // Distribute summed-up forces from physical particles to ghosts
   init_forces_ghosts(cell_structure.ghost_particles());
@@ -62,7 +61,7 @@ void lb_tracers_add_particle_force_to_fluid(CellStructure &cell_structure,
       if (bookkeeping.should_be_coupled(p)) {
         for (auto const &pos :
              positions_in_halo(p.pos(), box_geo, local_box, agrid)) {
-          add_md_force(lb, pos * to_lb_units, p.force(), time_step);
+          add_md_force(lb, pos, p.force(), time_step);
         }
       }
     }

testsuite/python/lb.py

@@ -588,6 +588,25 @@ def test_viscous_coupling_rounding(self):
             self.system.integrator.run(1)
             self.assertTrue(np.all(p.f != 0.0))
 
+    def test_tracers_coupling_rounding(self):
+        import espressomd.propagation
+        lbf = self.lb_class(**self.params, **self.lb_params)
+        self.system.lb = lbf
+        ext_f = np.array([0.01, 0.02, 0.03])
+        lbf.ext_force_density = ext_f
+        self.system.thermostat.set_lb(LB_fluid=lbf, seed=3, gamma=self.gamma)
+        rtol = self.rtol
+        if lbf.single_precision:
+            rtol *= 100.
+        mode_tracer = espressomd.propagation.Propagation.TRANS_LB_TRACER
+        self.system.time = 0.
+        p = self.system.part.add(pos=[-1E-30] * 3, propagation=mode_tracer)
+        for _ in range(10):
+            self.system.integrator.run(1)
+            vel = np.copy(p.v) * self.params["density"]
+            vel_ref = (self.system.time + lbf.tau * 0.5) * ext_f
+            np.testing.assert_allclose(vel, vel_ref, rtol=rtol, atol=0.)
+
     def test_thermalization_force_balance(self):
         system = self.system
 

testsuite/python/virtual_sites_tracers_common.py

@@ -29,8 +29,9 @@
 
 
 class VirtualSitesTracersCommon:
-    box_height = 10.
-    box_lw = 8.
+    agrid = 0.5
+    box_height = 10. * agrid
+    box_lw = 8. * agrid
     system = espressomd.System(box_l=(box_lw, box_lw, box_height))
     system.time_step = 0.08
     system.cell_system.skin = 0.1
@@ -46,19 +47,20 @@ def tearDown(self):
     def set_lb(self, ext_force_density=(0, 0, 0), dir_walls=2):
         self.system.lb = None
         self.lbf = self.LBClass(
-            kT=0.0, agrid=1., density=1., kinematic_viscosity=1.8,
+            kT=0.0, agrid=self.agrid, density=1., kinematic_viscosity=1.8,
             tau=self.system.time_step, ext_force_density=ext_force_density)
         self.system.lb = self.lbf
         self.system.thermostat.set_lb(LB_fluid=self.lbf, gamma=1.)
 
         # Setup boundaries
         normal = [0, 0, 0]
         normal[dir_walls] = 1
-        wall_shape = espressomd.shapes.Wall(normal=normal, dist=0.5)
+        wall_shape = espressomd.shapes.Wall(
+            normal=normal, dist=0.5 * self.agrid)
         self.lbf.add_boundary_from_shape(wall_shape)
         normal[dir_walls] = -1
         wall_shape = espressomd.shapes.Wall(
-            normal=normal, dist=-(self.system.box_l[dir_walls] - 0.5))
+            normal=normal, dist=-(self.system.box_l[dir_walls] - 0.5 * self.agrid))
         self.lbf.add_boundary_from_shape(wall_shape)
 
         espressomd.utils.handle_errors("setup")
@@ -72,7 +74,7 @@ def test_ab_single_step(self):
         self.lbf[:, :, :].velocity = np.random.random((*self.lbf.shape, 3))
         force = [1, -2, 3]
         # Test several particle positions
-        for pos in [[3, 2, 1], [0, 0, 0],
+        for pos in [[3 * self.agrid, 2 * self.agrid, 1 * self.agrid], [0, 0, 0],
                     self.system.box_l * 0.49,
                     self.system.box_l,
                     self.system.box_l * 0.99]:
@@ -133,8 +135,8 @@ def test_advection(self):
             system.integrator.run(400)
 
             # Add tracer in the fluid domain
-            pos_initial = [3.5, 3.5, 3.5]
-            pos_initial[direction] = 0.5
+            pos_initial = 3 * [3.5 * self.agrid]
+            pos_initial[direction] = 0.5 * self.agrid
             p = system.part.add(
                 pos=pos_initial,
                 propagation=espressomd.propagation.Propagation.TRANS_LB_TRACER)