diff --git a/src/core/lb/particle_coupling.cpp b/src/core/lb/particle_coupling.cpp index 48b1afa9f1..aca44c4a87 100644 --- a/src/core/lb/particle_coupling.cpp +++ b/src/core/lb/particle_coupling.cpp @@ -125,8 +125,9 @@ std::vector positions_in_halo(Utils::Vector3d const &pos, auto const halo_vec = Utils::Vector3d::broadcast(halo); auto const fully_inside_lower = local_box.my_left() + 2. * halo_vec; auto const fully_inside_upper = local_box.my_right() - 2. * halo_vec; - if (in_box(pos, fully_inside_lower, fully_inside_upper)) { - return {pos}; + auto const pos_folded = box_geo.folded_position(pos); + if (in_box(pos_folded, fully_inside_lower, fully_inside_upper)) { + return {pos_folded}; } auto const halo_lower_corner = local_box.my_left() - halo_vec; auto const halo_upper_corner = local_box.my_right() + halo_vec; @@ -137,11 +138,11 @@ std::vector positions_in_halo(Utils::Vector3d const &pos, for (int k : {-1, 0, 1}) { Utils::Vector3d shift{{double(i), double(j), double(k)}}; Utils::Vector3d pos_shifted = - pos + Utils::hadamard_product(box_geo.length(), shift); + pos_folded + Utils::hadamard_product(box_geo.length(), shift); if (box_geo.type() == BoxType::LEES_EDWARDS) { auto le = box_geo.lees_edwards_bc(); - auto normal_shift = (pos_shifted - pos)[le.shear_plane_normal]; + auto normal_shift = (pos_shifted - pos_folded)[le.shear_plane_normal]; if (normal_shift > std::numeric_limits::epsilon()) pos_shifted[le.shear_direction] += le.pos_offset; if (normal_shift < -std::numeric_limits::epsilon()) diff --git a/src/core/virtual_sites/lb_tracers.cpp b/src/core/virtual_sites/lb_tracers.cpp index 9edea02a1d..2618337494 100644 --- a/src/core/virtual_sites/lb_tracers.cpp +++ b/src/core/virtual_sites/lb_tracers.cpp @@ -44,7 +44,6 @@ void lb_tracers_add_particle_force_to_fluid(CellStructure &cell_structure, return; } auto const agrid = lb.get_agrid(); - auto const to_lb_units = 1. / agrid; // Distribute summed-up forces from physical particles to ghosts init_forces_ghosts(cell_structure.ghost_particles()); @@ -62,7 +61,7 @@ void lb_tracers_add_particle_force_to_fluid(CellStructure &cell_structure, if (bookkeeping.should_be_coupled(p)) { for (auto const &pos : positions_in_halo(p.pos(), box_geo, local_box, agrid)) { - add_md_force(lb, pos * to_lb_units, p.force(), time_step); + add_md_force(lb, pos, p.force(), time_step); } } } diff --git a/testsuite/python/lb.py b/testsuite/python/lb.py index 5fdb5bef9d..68c2b095b1 100644 --- a/testsuite/python/lb.py +++ b/testsuite/python/lb.py @@ -588,6 +588,25 @@ def test_viscous_coupling_rounding(self): self.system.integrator.run(1) self.assertTrue(np.all(p.f != 0.0)) + def test_tracers_coupling_rounding(self): + import espressomd.propagation + lbf = self.lb_class(**self.params, **self.lb_params) + self.system.lb = lbf + ext_f = np.array([0.01, 0.02, 0.03]) + lbf.ext_force_density = ext_f + self.system.thermostat.set_lb(LB_fluid=lbf, seed=3, gamma=self.gamma) + rtol = self.rtol + if lbf.single_precision: + rtol *= 100. + mode_tracer = espressomd.propagation.Propagation.TRANS_LB_TRACER + self.system.time = 0. + p = self.system.part.add(pos=[-1E-30] * 3, propagation=mode_tracer) + for _ in range(10): + self.system.integrator.run(1) + vel = np.copy(p.v) * self.params["density"] + vel_ref = (self.system.time + lbf.tau * 0.5) * ext_f + np.testing.assert_allclose(vel, vel_ref, rtol=rtol, atol=0.) + def test_thermalization_force_balance(self): system = self.system diff --git a/testsuite/python/virtual_sites_tracers_common.py b/testsuite/python/virtual_sites_tracers_common.py index 6c393d7801..85edaa57ae 100644 --- a/testsuite/python/virtual_sites_tracers_common.py +++ b/testsuite/python/virtual_sites_tracers_common.py @@ -29,8 +29,9 @@ class VirtualSitesTracersCommon: - box_height = 10. - box_lw = 8. + agrid = 0.5 + box_height = 10. * agrid + box_lw = 8. * agrid system = espressomd.System(box_l=(box_lw, box_lw, box_height)) system.time_step = 0.08 system.cell_system.skin = 0.1 @@ -46,7 +47,7 @@ def tearDown(self): def set_lb(self, ext_force_density=(0, 0, 0), dir_walls=2): self.system.lb = None self.lbf = self.LBClass( - kT=0.0, agrid=1., density=1., kinematic_viscosity=1.8, + kT=0.0, agrid=self.agrid, density=1., kinematic_viscosity=1.8, tau=self.system.time_step, ext_force_density=ext_force_density) self.system.lb = self.lbf self.system.thermostat.set_lb(LB_fluid=self.lbf, gamma=1.) @@ -54,11 +55,12 @@ def set_lb(self, ext_force_density=(0, 0, 0), dir_walls=2): # Setup boundaries normal = [0, 0, 0] normal[dir_walls] = 1 - wall_shape = espressomd.shapes.Wall(normal=normal, dist=0.5) + wall_shape = espressomd.shapes.Wall( + normal=normal, dist=0.5 * self.agrid) self.lbf.add_boundary_from_shape(wall_shape) normal[dir_walls] = -1 wall_shape = espressomd.shapes.Wall( - normal=normal, dist=-(self.system.box_l[dir_walls] - 0.5)) + normal=normal, dist=-(self.system.box_l[dir_walls] - 0.5 * self.agrid)) self.lbf.add_boundary_from_shape(wall_shape) espressomd.utils.handle_errors("setup") @@ -72,7 +74,7 @@ def test_ab_single_step(self): self.lbf[:, :, :].velocity = np.random.random((*self.lbf.shape, 3)) force = [1, -2, 3] # Test several particle positions - for pos in [[3, 2, 1], [0, 0, 0], + for pos in [[3 * self.agrid, 2 * self.agrid, 1 * self.agrid], [0, 0, 0], self.system.box_l * 0.49, self.system.box_l, self.system.box_l * 0.99]: @@ -133,8 +135,8 @@ def test_advection(self): system.integrator.run(400) # Add tracer in the fluid domain - pos_initial = [3.5, 3.5, 3.5] - pos_initial[direction] = 0.5 + pos_initial = 3 * [3.5 * self.agrid] + pos_initial[direction] = 0.5 * self.agrid p = system.part.add( pos=pos_initial, propagation=espressomd.propagation.Propagation.TRANS_LB_TRACER)