From 24a419747a8bec481b7a2bb09f413967e8f102a9 Mon Sep 17 00:00:00 2001 From: Rudolf Weeber Date: Fri, 24 Nov 2023 15:43:14 +0100 Subject: [PATCH] LB+particles: guard against missed coupling due to round-off error --- src/core/lb/particle_coupling.cpp | 6 ++++-- testsuite/python/lb.py | 11 +++++++++++ 2 files changed, 15 insertions(+), 2 deletions(-) diff --git a/src/core/lb/particle_coupling.cpp b/src/core/lb/particle_coupling.cpp index 63fa308d26..f8b2faf1c8 100644 --- a/src/core/lb/particle_coupling.cpp +++ b/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 */ diff --git a/testsuite/python/lb.py b/testsuite/python/lb.py index 592d7f434e..3cc7536ac6 100644 --- a/testsuite/python/lb.py +++ b/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