diff --git a/src/core/cell_system/CellStructure.cpp b/src/core/cell_system/CellStructure.cpp index 1c700736c3..91c5bee4a7 100644 --- a/src/core/cell_system/CellStructure.cpp +++ b/src/core/cell_system/CellStructure.cpp @@ -44,6 +44,7 @@ #include #include #include +#include #include #include #include @@ -249,8 +250,7 @@ void CellStructure::set_atom_decomposition() { } void CellStructure::set_regular_decomposition( - double range, - boost::optional> fully_connected_boundary) { + double range, std::optional> fully_connected_boundary) { auto &system = get_system(); auto &local_geo = *system.local_geo; auto const &box_geo = *system.box_geo; diff --git a/src/core/cell_system/CellStructure.hpp b/src/core/cell_system/CellStructure.hpp index 1cf449dfe6..9993fffc3d 100644 --- a/src/core/cell_system/CellStructure.hpp +++ b/src/core/cell_system/CellStructure.hpp @@ -46,6 +46,7 @@ #include #include #include +#include #include #include #include @@ -559,7 +560,7 @@ struct CellStructure : public System::Leaf { */ void set_regular_decomposition( double range, - boost::optional> fully_connected_boundary); + std::optional> fully_connected_boundary); /** * @brief Set the particle decomposition to @ref HybridDecomposition. diff --git a/src/core/cell_system/HybridDecomposition.cpp b/src/core/cell_system/HybridDecomposition.cpp index 4f3adcffcd..4c21fe6475 100644 --- a/src/core/cell_system/HybridDecomposition.cpp +++ b/src/core/cell_system/HybridDecomposition.cpp @@ -37,6 +37,7 @@ #include #include #include +#include #include #include @@ -48,7 +49,7 @@ HybridDecomposition::HybridDecomposition(boost::mpi::communicator comm, std::set n_square_types) : m_comm(std::move(comm)), m_box(box_geo), m_cutoff_regular(cutoff_regular), m_regular_decomposition(RegularDecomposition( - m_comm, cutoff_regular + skin, m_box, local_box, {})), + m_comm, cutoff_regular + skin, m_box, local_box, std::nullopt)), m_n_square(AtomDecomposition(m_comm, m_box)), m_n_square_types(std::move(n_square_types)), m_get_global_ghost_flags(std::move(get_ghost_flags)) { diff --git a/src/core/cell_system/RegularDecomposition.cpp b/src/core/cell_system/RegularDecomposition.cpp index dc66337d1c..df48efd35c 100644 --- a/src/core/cell_system/RegularDecomposition.cpp +++ b/src/core/cell_system/RegularDecomposition.cpp @@ -397,30 +397,29 @@ void RegularDecomposition::init_cell_interactions() { auto const &node_grid = ::communicator.node_grid; auto const global_halo_offset = hadamard_product(node_pos, cell_grid) - halo; auto const global_size = hadamard_product(node_grid, cell_grid); - auto at_boundary = [&global_size](int coord, Utils::Vector3i cell_idx) { + auto const at_boundary = [&global_size](int coord, Utils::Vector3i cell_idx) { return (cell_idx[coord] == 0 or cell_idx[coord] == global_size[coord]); }; // For the fully connected feature (cells that don't share at least a corner) // only apply if one cell is a ghost cell (i.e. connections across the - // periodic boudnary. - auto fcb_is_inner_connection = [&global_size, this](Utils::Vector3i a, - Utils::Vector3i b) { - if (not fully_connected_boundary()) - return false; - - const int fc_normal = (*fully_connected_boundary()).first; - const int fc_dir = (*fully_connected_boundary()).second; - - bool involves_ghost_cell = - (a[fc_normal] == -1 or a[fc_normal] == global_size[fc_normal] or - b[fc_normal] == -1 or b[fc_normal] == global_size[fc_normal]); - if (involves_ghost_cell) - return false; - return (abs((a - b)[fc_dir]) > 1); // cells do not share at least a corner + // periodic boundary. + auto const fcb_is_inner_connection = [&global_size, this](Utils::Vector3i a, + Utils::Vector3i b) { + if (fully_connected_boundary()) { + auto const [fc_normal, fc_dir] = *fully_connected_boundary(); + auto const involves_ghost_cell = + (a[fc_normal] == -1 or a[fc_normal] == global_size[fc_normal] or + b[fc_normal] == -1 or b[fc_normal] == global_size[fc_normal]); + if (not involves_ghost_cell) { + // check if cells do not share at least a corner + return std::abs((a - b)[fc_dir]) > 1; + } + } + return false; }; - /* Tanslate a node local index (relative to the origin of the local grid) + /* Translate a node local index (relative to the origin of the local grid) * to a global index. */ auto global_index = [&](Utils::Vector3i const &local_index) -> Utils::Vector3i { @@ -439,18 +438,19 @@ void RegularDecomposition::init_cell_interactions() { [&](Utils::Vector3i const &global_index) -> Utils::Vector3i { return (global_index - global_halo_offset); }; - // Validation + + // sanity checks if (fully_connected_boundary()) { - if ((*fully_connected_boundary()).first == - (*fully_connected_boundary()).second) { - throw std::runtime_error("fully_connectd_boundary normal and connection " - "coordinates need to differ."); + auto const [fc_normal, fc_dir] = *fully_connected_boundary(); + if (fc_normal == fc_dir) { + throw std::domain_error("fully_connected_boundary normal and connection " + "coordinates need to differ."); } - if (node_grid[(*fully_connected_boundary()).second] != 1) { + if (node_grid[fc_dir] != 1) { throw std::runtime_error( "The MPI nodegrid must be 1 in the fully connected direction."); - }; - }; + } + } /* We only consider local cells (e.g. not halo cells), which * span the range [(1,1,1), cell_grid) in local coordinates. */ @@ -469,8 +469,7 @@ void RegularDecomposition::init_cell_interactions() { * in the direction as neighbors, not only the nearest ones. // */ if (fully_connected_boundary()) { - int fc_boundary = (*fully_connected_boundary()).first; - int fc_direction = (*fully_connected_boundary()).second; + auto const [fc_boundary, fc_direction] = *fully_connected_boundary(); // Fully connected is only needed at the box surface if (at_boundary(fc_boundary, {m, n, o})) { @@ -671,7 +670,7 @@ GhostCommunicator RegularDecomposition::prepare_comm() { RegularDecomposition::RegularDecomposition( boost::mpi::communicator comm, double range, BoxGeometry const &box_geo, LocalBox const &local_geo, - boost::optional> fully_connected) + std::optional> fully_connected) : m_comm(std::move(comm)), m_box(box_geo), m_local_box(local_geo), m_fully_connected_boundary(std::move(fully_connected)) { diff --git a/src/core/cell_system/RegularDecomposition.hpp b/src/core/cell_system/RegularDecomposition.hpp index 9fe2a87479..f79a71ff10 100644 --- a/src/core/cell_system/RegularDecomposition.hpp +++ b/src/core/cell_system/RegularDecomposition.hpp @@ -79,7 +79,7 @@ struct RegularDecomposition : public ParticleDecomposition { boost::mpi::communicator m_comm; BoxGeometry const &m_box; LocalBox m_local_box; - boost::optional> m_fully_connected_boundary = {}; + std::optional> m_fully_connected_boundary = {}; std::vector cells; std::vector m_local_cells; std::vector m_ghost_cells; @@ -89,11 +89,8 @@ struct RegularDecomposition : public ParticleDecomposition { public: RegularDecomposition(boost::mpi::communicator comm, double range, BoxGeometry const &box_geo, LocalBox const &local_geo, - boost::optional> fully_connected); + std::optional> fully_connected); - boost::optional> fully_connected_boundary() const { - return m_fully_connected_boundary; - }; GhostCommunicator const &exchange_ghosts_comm() const override { return m_exchange_ghosts_comm; } @@ -120,6 +117,8 @@ struct RegularDecomposition : public ParticleDecomposition { Utils::Vector3d max_cutoff() const override; Utils::Vector3d max_range() const override; + auto fully_connected_boundary() const { return m_fully_connected_boundary; } + std::optional minimum_image_distance() const override { return {m_box}; } diff --git a/src/core/forces.cpp b/src/core/forces.cpp index 3580607fbc..bc317eb119 100644 --- a/src/core/forces.cpp +++ b/src/core/forces.cpp @@ -221,7 +221,6 @@ void System::System::calculate_forces() { lb_couple_particles(); } - particles = cell_structure->local_particles(); #ifdef CUDA #ifdef CALIPER CALI_MARK_BEGIN("copy_forces_from_GPU"); diff --git a/src/python/espressomd/cell_system.py b/src/python/espressomd/cell_system.py index 2befdf16b0..17a5364337 100644 --- a/src/python/espressomd/cell_system.py +++ b/src/python/espressomd/cell_system.py @@ -115,14 +115,14 @@ def set_regular_decomposition(self, **kwargs): use_verlet_lists : :obj:`bool`, optional Activates or deactivates the usage of Verlet lists. Defaults to ``True``. - fully_connected_boundary : dict, optional - If set, connects all cells on a given boundary along the given direciton. - Example: {"boundary":"z", "direction":"x"} connects all cells on the boundary normal to - the z-direcion along the x axis. This corresponds to - z as shear plane normal and x as shear direction in Lees-Edwards - boundary conditions. - """ + fully_connected_boundary : :obj:`dict`, optional + If set, connects all cells on a given boundary along the given direction. + Example: ``{"boundary": "z", "direction": "x"}`` connects all + cells on the boundary normal to the z-direction along the x-axis. + This corresponds to z-axis as shear plane normal and x-axis as + shear direction in Lees-Edwards boundary conditions. + """ self.call_method("initialize", name="regular_decomposition", **kwargs) def set_n_square(self, **kwargs): diff --git a/src/script_interface/cell_system/CellSystem.cpp b/src/script_interface/cell_system/CellSystem.cpp index 4af8d63c3d..1e30b0e7bc 100644 --- a/src/script_interface/cell_system/CellSystem.cpp +++ b/src/script_interface/cell_system/CellSystem.cpp @@ -40,7 +40,9 @@ #include #include +#include #include +#include #include #include #include @@ -49,29 +51,25 @@ #include #include -namespace { - -int coord(const std::string &s) { +static int coord(std::string const &s) { if (s == "x") return 0; if (s == "y") return 1; if (s == "z") return 2; - throw std::runtime_error("Invalid argument for Cartesian coordinate: " + s); + throw std::invalid_argument("Invalid Cartesian coordinate: '" + s + "'"); } -std::string coord_letter(int c) { +static std::string coord_letter(int c) { if (c == 0) return "x"; if (c == 1) return "y"; - if (c == 2) - return "z"; - throw std::runtime_error("Invalid numeric coordinate"); + assert(c == 2); + return "z"; } -} // namespace namespace ScriptInterface { namespace CellSystem { @@ -299,19 +297,14 @@ void CellSystem::initialize(CellStructureType const &cs_type, auto n_square_types = std::set{ns_types.begin(), ns_types.end()}; m_cell_structure->set_hybrid_decomposition(cutoff_regular, n_square_types); } else if (cs_type == CellStructureType::REGULAR) { - boost::optional> fcb_pair = {}; - if (params.count("fully_connected_boundary") == 0 || - is_type(params.at("fully_connected_boundary"))) { - // pass - } else { - auto const fully_connected_boundary = - get_value>( - params, "fully_connected_boundary"); - const int boundary = coord( - boost::get(fully_connected_boundary.at("boundary"))); - const int direction = coord( - boost::get(fully_connected_boundary.at("direction"))); - fcb_pair = {{boundary, direction}}; + std::optional> fcb_pair = std::nullopt; + if (params.contains("fully_connected_boundary") and + not is_none(params.at("fully_connected_boundary"))) { + auto const variant = get_value(params, "fully_connected_boundary"); + context()->parallel_try_catch([&fcb_pair, &variant]() { + fcb_pair = {{coord(boost::get(variant.at("boundary"))), + coord(boost::get(variant.at("direction")))}}; + }); } context()->parallel_try_catch([this, &fcb_pair]() { m_cell_structure->set_regular_decomposition( diff --git a/testsuite/python/lees_edwards.py b/testsuite/python/lees_edwards.py index 384c24dd36..4ebdb2dafb 100644 --- a/testsuite/python/lees_edwards.py +++ b/testsuite/python/lees_edwards.py @@ -57,26 +57,29 @@ def axis(coord): class LeesEdwards(ut.TestCase): box_l = [5, 5, 5] system = espressomd.System(box_l=box_l) + node_grid = np.copy(system.cell_system.node_grid) + n_nodes = np.prod(node_grid) direction_permutations = list(itertools.permutations(["x", "y", "z"], 2)) def setUp(self): system = self.system - system.cell_system.skin = 0.0 + system.box_l = self.box_l + system.cell_system.skin = 0. system.cell_system.set_n_square(use_verlet_lists=True) system.time = 0.0 system.time_step = 0.5 + system.min_global_cut = 0. + system.cell_system.node_grid = self.node_grid + + def tearDown(self): + system = self.system system.part.clear() - system.box_l = self.box_l system.bonded_inter.clear() - system.min_global_cut = 0 - system.cell_system.skin = 0 + system.lees_edwards.protocol = None if espressomd.has_features("COLLISION_DETECTION"): system.collision_detection.set_params(mode="off") - def tearDown(self): - self.system.lees_edwards.protocol = None - def test_00_is_none_by_default(self): system = self.system @@ -196,6 +199,21 @@ def test_protocols(self): shear_direction=valid, shear_plane_normal=valid, protocol=lin_protocol) + with self.assertRaisesRegex(ValueError, "fully_connected_boundary normal and connection coordinates need to differ"): + system.cell_system.set_regular_decomposition( + fully_connected_boundary={"boundary": "z", "direction": "z"}) + self.assertEqual(system.cell_system.decomposition_type, "n_square") + with self.assertRaisesRegex(ValueError, "Invalid Cartesian coordinate: 't'"): + system.cell_system.set_regular_decomposition( + fully_connected_boundary={"boundary": "z", "direction": "t"}) + self.assertEqual(system.cell_system.decomposition_type, "n_square") + if self.n_nodes > 1: + with self.assertRaisesRegex(RuntimeError, "The MPI nodegrid must be 1 in the fully connected direction"): + system.cell_system.node_grid = [1, self.n_nodes, 1] + system.cell_system.set_regular_decomposition( + fully_connected_boundary={"boundary": "z", "direction": "y"}) + self.assertEqual(system.cell_system.decomposition_type, "n_square") + def test_boundary_crossing_lin(self): """ A particle crosses the upper and lower boundary with linear shear. @@ -544,9 +562,11 @@ def test_virt_sites_interaction(self): @utx.skipIfMissingFeatures( ["EXTERNAL_FORCES", "VIRTUAL_SITES_RELATIVE"]) def test__virt_sites_rotation(self): - """A particle with virtual sites is plces on the boundary. We check if - the forces yield the correct torque and if a rotation frequency is - transmitted backl to the virtual sites.""" + """ + A particle with virtual sites is placed on the boundary. We check if + the forces yield the correct torque and if a rotation frequency is + transmitted back to the virtual sites. + """ system = self.system system.part.clear() @@ -742,7 +762,7 @@ def run_lj_pair_visibility(self, shear_direction, shear_plane_normal): """ Simulate LJ particles coming into contact under linear shear and verify forces. This is to make sure that no pairs get lost or are outdated - in the short range loop. + in the short range loop. """ shear_axis, normal_axis = axis( shear_direction), axis(shear_plane_normal) @@ -779,32 +799,37 @@ def run_lj_pair_visibility(self, shear_direction, shear_plane_normal): assert have_interacted def test_zz_lj_pair_visibility(self): - # check that regular decomposition without fully connected doesn't catch the particle + # check that regular decomposition without fully connected doesn't + # catch the particle system = self.system - system.box_l = 10, 10, 10 + system.box_l = [10, 10, 10] with self.assertRaises(AssertionError): system.cell_system.set_regular_decomposition( fully_connected_boundary=None) - n_nodes = system.cell_system.get_state()["n_nodes"] - system.cell_system.node_grid = [1, n_nodes, 1] + self.assertIsNone(system.cell_system.fully_connected_boundary) + system.cell_system.node_grid = [1, self.n_nodes, 1] self.run_lj_pair_visibility("x", "y") - for verlet in False, True: + for verlet in (False, True): for shear_direction, shear_plane_normal in self.direction_permutations: system.cell_system.set_n_square(use_verlet_lists=verlet) self.run_lj_pair_visibility( shear_direction, shear_plane_normal) - for verlet in False, True: + for verlet in (False, True): for shear_direction, shear_plane_normal in self.direction_permutations: system.cell_system.set_regular_decomposition( fully_connected_boundary=None) - n_nodes = system.cell_system.get_state()["n_nodes"] normal_axis = axis(shear_plane_normal) system.cell_system.node_grid = [ - n_nodes if normal_axis[i] == 1 else 1 for i in range(3)] + self.n_nodes if normal_axis[i] == 1 else 1 for i in range(3)] + fully_connected_boundary = {"boundary": shear_plane_normal, + "direction": shear_direction} system.cell_system.set_regular_decomposition( - fully_connected_boundary={"boundary": shear_plane_normal, "direction": shear_direction}, use_verlet_lists=verlet) + use_verlet_lists=verlet, + fully_connected_boundary=fully_connected_boundary) + self.assertEqual(system.cell_system.fully_connected_boundary, + fully_connected_boundary) self.run_lj_pair_visibility( shear_direction, shear_plane_normal) diff --git a/testsuite/python/regular_decomposition.py b/testsuite/python/regular_decomposition.py index a15f8694aa..d94720f11e 100644 --- a/testsuite/python/regular_decomposition.py +++ b/testsuite/python/regular_decomposition.py @@ -115,7 +115,7 @@ def test_fully_connected_boundary(self): ng = system.cell_system.node_grid system.cell_system.node_grid = [ng[0], 1, ng[2] * ng[1]] system.periodic = [True] * 3 - # Check fthat it's initially disabled + # Check that it's initially disabled self.assertEqual(system.cell_system.get_params()[ "fully_connected_boundary"], None) @@ -124,15 +124,14 @@ def test_fully_connected_boundary(self): fully_connected_boundary=dict(direction="y", boundary="z")) self.assertEqual(system.cell_system.get_params()[ "fully_connected_boundary"], dict(direction="y", boundary="z")) - # Check that the setting survies cell system re-intis + # Check that the setting survives cell system re-initialization system.cell_system.min_global_cut = system.box_l / 4.1 self.assertEqual(system.cell_system.get_params()[ "fully_connected_boundary"], dict(direction="y", boundary="z")) - # Check that particle visibility + # Check particle visibility. # Place particles on a cubic lattice and use the - # non_bonded_loop_trace() to check that all pairs are seen - # as expected + # non_bonded_loop_trace() to check that all pairs are seen as expected fc_normal = np.array((0, 0, 1)) # z fc_normal_coord = 2 # z fc_dir = np.array((0, 1, 0)) # y @@ -144,6 +143,7 @@ def test_fully_connected_boundary(self): def id_for_idx(idx): return ( idx[0] % N) * N * N + (idx[1] % N) * N + idx[2] % N + ids = [id_for_idx(idx) for idx in indices] dx = system.box_l / N positions = [idx * dx for idx in indices] @@ -165,7 +165,7 @@ def distance(id1, id2): # next neighbors must_find_nn = [i for i, d in distances.items() if d <= max_range] - # Fully connected neighbos + # Fully connected neighbors indices_lower_boundary = [ idx for idx in indices if idx[fc_normal_coord] == 0] must_find_fc = [tuple(sorted((id_for_idx(idx), id_for_idx(idx + i * fc_dir - fc_normal)))) @@ -197,9 +197,9 @@ def assert_can_find(pair): cs_pairs = system.cell_system.non_bonded_loop_trace() found = [] for id1, id2, _rest1, _rest2, _rest3, _rest4 in cs_pairs: - p = tuple(sorted((id1, id2))) # Make teh pair unique - found.append(p) # to check for double countin - if (p in must_find): + p = tuple(sorted((id1, id2))) # Make the pair unique + found.append(p) # to check for double counting + if p in must_find: must_find.remove(p) else: assert_can_find(p) # close enough so that cells share a corner