From 24a419747a8bec481b7a2bb09f413967e8f102a9 Mon Sep 17 00:00:00 2001 From: Rudolf Weeber Date: Fri, 24 Nov 2023 15:43:14 +0100 Subject: [PATCH 1/2] 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 From f4181475b8b4bf572cb6d95602278033840e8444 Mon Sep 17 00:00:00 2001 From: Rudolf Weeber Date: Mon, 27 Nov 2023 17:22:14 +0100 Subject: [PATCH 2/2] Formatting --- src/core/lb/particle_coupling.cpp | 1 - testsuite/python/lb.py | 7 +++---- 2 files changed, 3 insertions(+), 5 deletions(-) diff --git a/src/core/lb/particle_coupling.cpp b/src/core/lb/particle_coupling.cpp index f8b2faf1c8..416522eae5 100644 --- a/src/core/lb/particle_coupling.cpp +++ b/src/core/lb/particle_coupling.cpp @@ -251,7 +251,6 @@ void ParticleCoupling::kernel(Particle &p) { 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 diff --git a/testsuite/python/lb.py b/testsuite/python/lb.py index 3cc7536ac6..f0d0cd4d52 100644 --- a/testsuite/python/lb.py +++ b/testsuite/python/lb.py @@ -565,18 +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]) + p = self.system.part.add(pos=[-1E-30] * 3, v=[-1, 0, 0]) self.system.integrator.run(1) - for i in range(20): + 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