diff --git a/doc/sphinx/particles.rst b/doc/sphinx/particles.rst index b1da59f0d94..7bbb7b46052 100644 --- a/doc/sphinx/particles.rst +++ b/doc/sphinx/particles.rst @@ -386,6 +386,7 @@ For correct results, the LB thermostat has to be deactivated for virtual sites:: system.thermostat.set_lb(kT=0, act_on_virtual=False) Please note that the velocity attribute of the virtual particles does not carry valid information for this virtual sites scheme. +With the LB GPU implementation, inertialess tracers only work on 1 MPI rank. .. _Interacting with groups of particles: diff --git a/src/core/virtual_sites/VirtualSitesInertialessTracers.cpp b/src/core/virtual_sites/VirtualSitesInertialessTracers.cpp index 30dd822599d..f22415d1499 100644 --- a/src/core/virtual_sites/VirtualSitesInertialessTracers.cpp +++ b/src/core/virtual_sites/VirtualSitesInertialessTracers.cpp @@ -30,6 +30,14 @@ #include +static void check_no_vs_exist(char const *const message) { + if (std::any_of(cell_structure.local_particles().begin(), + cell_structure.local_particles().end(), + [](Particle const &p) { return p.is_virtual(); })) { + runtimeErrorMsg() << "Inertialess Tracers: " << message; + } +} + void VirtualSitesInertialessTracers::after_force_calc() { // Now the forces are computed and need to go into the LB fluid if (lattice_switch == ActiveLB::CPU) { @@ -39,16 +47,15 @@ void VirtualSitesInertialessTracers::after_force_calc() { #ifdef CUDA if (lattice_switch == ActiveLB::GPU) { IBM_ForcesIntoFluid_GPU(cell_structure.local_particles(), this_node); + if (comm_cart.size() != 1 and this_node != 0) { + check_no_vs_exist("The LB GPU method cannot integrate virtual sites when " + "more than 1 MPI ranks are used. The particles on MPI " + "rank >= 2 are now in an undeterminate state."); + } return; } #endif - if (std::any_of(cell_structure.local_particles().begin(), - cell_structure.local_particles().end(), - [](Particle &p) { return p.is_virtual(); })) { - runtimeErrorMsg() << "Inertialess Tracers: No LB method was active but " - "virtual sites present."; - return; - } + check_no_vs_exist("No LB method was active but virtual sites present."); } void VirtualSitesInertialessTracers::after_lb_propagation(double time_step) { diff --git a/src/core/virtual_sites/lb_inertialess_tracers.cpp b/src/core/virtual_sites/lb_inertialess_tracers.cpp index 5aae5af66c5..b87b2d42327 100644 --- a/src/core/virtual_sites/lb_inertialess_tracers.cpp +++ b/src/core/virtual_sites/lb_inertialess_tracers.cpp @@ -92,7 +92,7 @@ void IBM_UpdateParticlePositions(ParticleRange const &particles, // Euler integrator for (auto &p : particles) { if (p.is_virtual()) { - for (int axis = 0; axis < 2; axis++) { + for (int axis = 0; axis < 3; axis++) { #ifdef EXTERNAL_FORCES if (not p.is_fixed_along(axis)) #endif diff --git a/testsuite/python/virtual_sites_tracers.py b/testsuite/python/virtual_sites_tracers.py index 0cf9e631733..811b51f00f1 100644 --- a/testsuite/python/virtual_sites_tracers.py +++ b/testsuite/python/virtual_sites_tracers.py @@ -25,13 +25,9 @@ @utx.skipIfMissingFeatures( ['VIRTUAL_SITES_INERTIALESS_TRACERS', 'LB_BOUNDARIES']) -class VirtualSitesTracers(ut.TestCase, VirtualSitesTracersCommon): +class VirtualSitesTracers(VirtualSitesTracersCommon, ut.TestCase): - def setUp(self): - self.LBClass = espressomd.lb.LBFluid - - def tearDown(self): - VirtualSitesTracersCommon.tearDown(self) + LBClass = espressomd.lb.LBFluid if __name__ == "__main__": diff --git a/testsuite/python/virtual_sites_tracers_common.py b/testsuite/python/virtual_sites_tracers_common.py index 054b0ece56f..eb61789db55 100644 --- a/testsuite/python/virtual_sites_tracers_common.py +++ b/testsuite/python/virtual_sites_tracers_common.py @@ -17,6 +17,7 @@ # along with this program. If not, see . # import espressomd +import espressomd.lb import espressomd.shapes import espressomd.lbboundaries import espressomd.virtual_sites @@ -30,13 +31,18 @@ class VirtualSitesTracersCommon: system.time_step = 0.05 system.cell_system.skin = 0.1 + def setUp(self): + self.system.box_l = (self.box_lw, self.box_lw, self.box_height) + def tearDown(self): - self.system.part.clear() - self.system.actors.clear() self.system.thermostat.turn_off() + self.system.lbboundaries.clear() + self.system.actors.clear() + self.system.part.clear() - def reset_lb(self, ext_force_density=(0, 0, 0)): + def reset_lb(self, ext_force_density=(0, 0, 0), dir_walls=2): self.system.lbboundaries.clear() + self.system.actors.clear() self.lbf = self.LBClass( kT=0.0, agrid=1, dens=1, visc=1.8, tau=self.system.time_step, ext_force_density=ext_force_density) @@ -47,11 +53,14 @@ def reset_lb(self, ext_force_density=(0, 0, 0)): gamma=1) # Setup boundaries + normal = [0, 0, 0] + normal[dir_walls] = 1 walls = [espressomd.lbboundaries.LBBoundary() for k in range(2)] walls[0].set_params(shape=espressomd.shapes.Wall( - normal=[0, 0, 1], dist=0.5)) + normal=normal, dist=0.5)) + normal[dir_walls] = -1 walls[1].set_params(shape=espressomd.shapes.Wall( - normal=[0, 0, -1], dist=-self.box_height - 0.5)) + normal=normal, dist=-(self.system.box_l[dir_walls] - 0.5))) for wall in walls: self.system.lbboundaries.add(wall) @@ -70,27 +79,59 @@ def test_aa_method_switching(self): self.system.virtual_sites, espressomd.virtual_sites.VirtualSitesInertialessTracers) def test_advection(self): - self.reset_lb(ext_force_density=[0.1, 0, 0]) - # System setup - system = self.system - - system.virtual_sites = espressomd.virtual_sites.VirtualSitesInertialessTracers() - - # Establish steady state flow field - p = system.part.add(pos=(0, 5.5, 5.5), virtual=True) - system.integrator.run(400) - - p.pos = (0, 5.5, 5.5) - system.time = 0 - - # Perform integration - for _ in range(2): - system.integrator.run(100) - # compute expected position - dist = self.lbf.get_interpolated_velocity(p.pos)[0] * system.time - self.assertAlmostEqual(p.pos[0] / dist, 1, delta=0.005) + node_grid = self.system.cell_system.node_grid + lb_on_gpu = self.LBClass is espressomd.lb.LBFluidGPU + for direction in [0, 1, 2]: + # System setup + system = self.system + system.virtual_sites = espressomd.virtual_sites.VirtualSitesInertialessTracers() + + # LB setup with walls + ext_force = [0., 0., 0.] + ext_force[direction] = 0.1 + dir_walls = (direction + 2) % 3 + box_l = 3 * [self.box_lw] + box_l[dir_walls] = self.box_height + system.box_l = box_l + self.reset_lb(ext_force_density=ext_force, dir_walls=dir_walls) + + # Establish steady state flow field + system.integrator.run(400) + + # Add tracer in the fluid domain + pos_initial = [3.5, 3.5, 3.5] + pos_initial[direction] = 0.5 + p = system.part.add(pos=pos_initial, virtual=True) + + # Perform integration + n_steps = 100 + if node_grid[direction] != 1 and lb_on_gpu: + n_steps = 50 # with GPU, tracers must stay on MPI rank 0 + system.time = 0 + for _ in range(2): + system.integrator.run(n_steps) + # compute expected position + lb_vel = self.lbf.get_interpolated_velocity(p.pos) + ref_dist = lb_vel[direction] * system.time + tracer_dist = p.pos[direction] - pos_initial[direction] + self.assertAlmostEqual(tracer_dist / ref_dist, 1., delta=0.01) + + self.tearDown() + + def test_zz_exceptions_with_lb(self): + node_grid = self.system.cell_system.node_grid + lb_on_gpu = self.LBClass is espressomd.lb.LBFluidGPU + if lb_on_gpu and sum(node_grid) != 3: + self.reset_lb() + system = self.system + system.virtual_sites = espressomd.virtual_sites.VirtualSitesInertialessTracers() + p = system.part.add(pos=(0, 0, 0), virtual=True) + system.integrator.run(1) + p.pos = (-0.5, -0.5, -0.5) + with self.assertRaisesRegex(Exception, "The LB GPU method cannot integrate virtual sites when more than 1 MPI ranks are used"): + system.integrator.run(1) - def test_zz_without_lb(self): + def test_zz_exceptions_without_lb(self): """Check behaviour without lb. Ignore non-virtual particles, complain on virtual ones. @@ -103,5 +144,5 @@ def test_zz_without_lb(self): p = system.part.add(pos=(0, 0, 0)) system.integrator.run(1) p.virtual = True - with self.assertRaises(Exception): + with self.assertRaisesRegex(Exception, "No LB method was active but virtual sites present"): system.integrator.run(1) diff --git a/testsuite/python/virtual_sites_tracers_gpu.py b/testsuite/python/virtual_sites_tracers_gpu.py index 12a2670509d..57d4c6873bc 100644 --- a/testsuite/python/virtual_sites_tracers_gpu.py +++ b/testsuite/python/virtual_sites_tracers_gpu.py @@ -26,13 +26,9 @@ @utx.skipIfMissingGPU() @utx.skipIfMissingFeatures( ['VIRTUAL_SITES_INERTIALESS_TRACERS', 'LB_BOUNDARIES']) -class VirtualSitesTracers(ut.TestCase, VirtualSitesTracersCommon): +class VirtualSitesTracers(VirtualSitesTracersCommon, ut.TestCase): - def setUp(self): - self.LBClass = espressomd.lb.LBFluidGPU - - def tearDown(self): - VirtualSitesTracersCommon.tearDown(self) + LBClass = espressomd.lb.LBFluidGPU if __name__ == "__main__":