From bcc87872cd4751278df4dfa35fe3e365b49ec06f Mon Sep 17 00:00:00 2001
From: "kodiakhq[bot]" <49736102+kodiakhq[bot]@users.noreply.github.com>
Date: Tue, 5 Dec 2023 18:26:00 +0000
Subject: [PATCH] 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.
---
 .../lb_particle_coupling.cpp                      |  6 ++++--
 testsuite/python/lb.py                            | 15 +++++++++++++++
 2 files changed, 19 insertions(+), 2 deletions(-)

diff --git a/src/core/grid_based_algorithms/lb_particle_coupling.cpp b/src/core/grid_based_algorithms/lb_particle_coupling.cpp
index 9fa0e16d57f..e6591a92c54 100644
--- a/src/core/grid_based_algorithms/lb_particle_coupling.cpp
+++ b/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 */
diff --git a/testsuite/python/lb.py b/testsuite/python/lb.py
index b15f9d090f6..8244d206ae6 100644
--- a/testsuite/python/lb.py
+++ b/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]