LB+particles: guard against missed coupling due to round-off error (#4827)

Fixes #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.

Author: kodiakhq[bot]

src/core/lb/particle_coupling.cpp

@@ -249,11 +249,12 @@ void ParticleCoupling::kernel(Particle &p) {
 
   // Calculate coupling force
   Utils::Vector3d force_on_particle = {};
+  auto folded_pos = box_geo.folded_position(p.pos());
 
 #ifdef ENGINE
   if (not p.swimming().is_engine_force_on_fluid)
 #endif
-    for (auto pos : positions_in_halo(p.pos(), box_geo, agrid)) {
+    for (auto pos : positions_in_halo(folded_pos, box_geo, agrid)) {
       if (in_local_halo(pos, agrid)) {
         auto const vel_offset = lb_drift_velocity_offset(p);
         auto const drag_force = lb_drag_force(m_lb, p, pos, vel_offset);
@@ -272,7 +273,7 @@ void ParticleCoupling::kernel(Particle &p) {
 
   // couple positions including shifts by one box length to add
   // forces to ghost layers
-  for (auto pos : positions_in_halo(p.pos(), box_geo, agrid)) {
+  for (auto pos : positions_in_halo(folded_pos, box_geo, agrid)) {
     if (in_local_domain(pos)) {
       /* Particle is in our LB volume, so this node
        * is responsible to adding its force */

testsuite/python/lb.py

@@ -566,6 +566,16 @@ def test_viscous_coupling_pairs(self):
             np.testing.assert_allclose(
                 np.sum(applied_forces, axis=0), [0, 0, 0])
 
+    def test_viscous_coupling_rounding(self):
+        lbf = self.lb_class(**self.params, **self.lb_params)
+        self.system.lb = lbf
+        self.system.thermostat.set_lb(LB_fluid=lbf, seed=3, gamma=self.gamma)
+        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))
+
     def test_thermalization_force_balance(self):
         system = self.system