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

src/core/lb/particle_coupling.cpp

@@ -249,11 +249,13 @@ 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 +274,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

@@ -565,6 +565,17 @@ def test_viscous_coupling_pairs(self):
                 [n.last_applied_force for n in all_coupling_nodes])
             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 i 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