LB+particles: guard against missed coupling due to round-off error (e…

…spressomd#4827)

Fixes espressomd#4825

Description of changes:
* Bugfix: particles outside the simulation box are now properly coupled using PBC. The coordinates are folded before considering shifted positions in the LB particle coupling code. Also, a test that fails without the fix is added. To my understanding, if the particle position is folded, and then all combination of `folded_pos[i]` +/- `box_l[i]` for all Cartesian directions are considered, both the particle in the primary box and all potential halo regions are caught.

src/core/grid_based_algorithms/lb_particle_coupling.cpp

@@ -312,9 +312,11 @@ void lb_lbcoupling_calc_particle_lattice_ia(bool couple_virtual,
           if (p.is_virtual() and !couple_virtual)
             return;
 
+          auto const folded_pos = folded_position(p.pos(), box_geo);
+
           // Calculate coupling force
           Utils::Vector3d force = {};
-          for (auto pos : positions_in_halo(p.pos(), box_geo)) {
+          for (auto pos : positions_in_halo(folded_pos, box_geo)) {
             if (in_local_halo(pos)) {
               force = lb_viscous_coupling(p, pos,
                                           noise_amplitude * f_random(p.id()));
@@ -324,7 +326,7 @@ void lb_lbcoupling_calc_particle_lattice_ia(bool couple_virtual,
 
           // couple positions including shifts by one box length to add
           // forces to ghost layers
-          for (auto pos : positions_in_halo(p.pos(), box_geo)) {
+          for (auto pos : positions_in_halo(folded_pos, box_geo)) {
             if (in_local_domain(pos)) {
               /* if the particle is in our LB volume, this node
                * is responsible to adding its force */

testsuite/python/lb.py

@@ -385,6 +385,21 @@ def test_viscous_coupling(self):
         np.testing.assert_allclose(
             np.copy(p.f), -self.params['friction'] * (v_part - v_fluid), atol=1E-6)
 
+    def test_viscous_coupling_rounding(self):
+        lbf = self.lb_class(
+            visc=self.params['viscosity'],
+            dens=self.params['dens'],
+            agrid=self.params['agrid'],
+            tau=self.params['time_step'],
+            kT=1., seed=1)
+        self.system.actors.add(lbf)
+        self.system.thermostat.set_lb(LB_fluid=lbf, gamma=0.1, seed=1)
+        p = self.system.part.add(pos=[-1E-30] * 3, v=[-1, 0, 0])
+        self.system.integrator.run(1)
+        for _ in range(20):
+            self.system.integrator.run(1)
+            self.assertTrue(np.all(p.f != 0.0))
+
     @utx.skipIfMissingFeatures("EXTERNAL_FORCES")
     def test_ext_force_density(self):
         ext_force_density = [2.3, 1.2, 0.1]