From 42730857ebaf626c00b53157483344358c23ce35 Mon Sep 17 00:00:00 2001 From: Patrick Kreissl Date: Wed, 20 Oct 2021 12:29:35 +0200 Subject: [PATCH 1/3] Implement HybridDecomposition. --- src/core/AtomDecomposition.hpp | 4 + src/core/CMakeLists.txt | 1 + src/core/CellStructure.cpp | 13 +- src/core/CellStructure.hpp | 24 ++ src/core/CellStructureType.hpp | 4 +- src/core/HybridDecomposition.cpp | 164 +++++++++++ src/core/HybridDecomposition.hpp | 125 +++++++++ src/core/RegularDecomposition.hpp | 4 + src/core/cells.cpp | 71 +++++ src/core/cells.hpp | 48 ++++ src/core/event.cpp | 7 +- src/core/grid_based_algorithms/lb.cpp | 4 +- src/python/espressomd/cellsystem.pxd | 15 +- src/python/espressomd/cellsystem.pyx | 46 ++++ testsuite/python/CMakeLists.txt | 10 + testsuite/python/cellsystem.py | 19 +- testsuite/python/engine_lb.py | 9 +- testsuite/python/engine_lb_hybrid.py | 77 ++++++ testsuite/python/engine_lb_nsquare.py | 2 - testsuite/python/hybrid_decomposition.py | 258 ++++++++++++++++++ .../python/lb_momentum_conservation_hybrid.py | 72 +++++ .../lb_momentum_conservation_nsquare.py | 2 - testsuite/python/lb_stats.py | 8 +- testsuite/python/lb_stats_hybrid.py | 47 ++++ testsuite/python/lb_stats_nsquare.py | 2 - testsuite/python/regular_decomposition.py | 4 +- 26 files changed, 1017 insertions(+), 23 deletions(-) create mode 100644 src/core/HybridDecomposition.cpp create mode 100644 src/core/HybridDecomposition.hpp create mode 100644 testsuite/python/engine_lb_hybrid.py create mode 100644 testsuite/python/hybrid_decomposition.py create mode 100644 testsuite/python/lb_momentum_conservation_hybrid.py create mode 100644 testsuite/python/lb_stats_hybrid.py diff --git a/src/core/AtomDecomposition.hpp b/src/core/AtomDecomposition.hpp index 006292f9b57..0b07ba1fcd5 100644 --- a/src/core/AtomDecomposition.hpp +++ b/src/core/AtomDecomposition.hpp @@ -80,6 +80,10 @@ class AtomDecomposition : public ParticleDecomposition { return Utils::make_span(m_ghost_cells); } + /* Getter needed for HybridDecomposition */ + std::vector get_local_cells() const { return m_local_cells; } + std::vector get_ghost_cells() const { return m_ghost_cells; } + /** * @brief Determine which cell a particle id belongs to. * diff --git a/src/core/CMakeLists.txt b/src/core/CMakeLists.txt index 3dc075f4d7b..92a74af5d99 100644 --- a/src/core/CMakeLists.txt +++ b/src/core/CMakeLists.txt @@ -15,6 +15,7 @@ set(EspressoCore_SRC galilei.cpp ghosts.cpp grid.cpp + HybridDecomposition.cpp immersed_boundaries.cpp interactions.cpp event.cpp diff --git a/src/core/CellStructure.cpp b/src/core/CellStructure.cpp index 0df13eab4f4..71c2a825ded 100644 --- a/src/core/CellStructure.cpp +++ b/src/core/CellStructure.cpp @@ -23,6 +23,7 @@ #include "AtomDecomposition.hpp" #include "CellStructureType.hpp" +#include "HybridDecomposition.hpp" #include "RegularDecomposition.hpp" #include "grid.hpp" #include "lees_edwards/lees_edwards.hpp" @@ -154,7 +155,7 @@ void CellStructure::remove_all_particles() { /* Map the data parts flags from cells to those used internally * by the ghost communication */ -static unsigned map_data_parts(unsigned data_parts) { +unsigned map_data_parts(unsigned data_parts) { using namespace Cells; /* clang-format off */ @@ -260,3 +261,13 @@ void CellStructure::set_regular_decomposition( m_type = CellStructureType::CELL_STRUCTURE_REGULAR; local_geo.set_cell_structure_type(m_type); } + +void CellStructure::set_hybrid_decomposition( + boost::mpi::communicator const &comm, double cutoff_regular, + BoxGeometry const &box, LocalBox &local_geo, + std::set n_square_types) { + set_particle_decomposition(std::make_unique( + comm, cutoff_regular, box, local_geo, n_square_types)); + m_type = CellStructureType::CELL_STRUCTURE_HYBRID; + local_geo.set_cell_structure_type(m_type); +} diff --git a/src/core/CellStructure.hpp b/src/core/CellStructure.hpp index fd7d77fcd09..5d464237e17 100644 --- a/src/core/CellStructure.hpp +++ b/src/core/CellStructure.hpp @@ -26,6 +26,7 @@ #include "BoxGeometry.hpp" #include "Cell.hpp" #include "CellStructureType.hpp" +#include "HybridDecomposition.hpp" #include "LocalBox.hpp" #include "Particle.hpp" #include "ParticleDecomposition.hpp" @@ -74,6 +75,15 @@ enum DataPart : unsigned { }; } // namespace Cells +/** + * @brief Map the data parts flags from cells to those + * used internally by the ghost communication. + * + * @param data_parts data parts flags + * @return ghost communication flags + */ +unsigned map_data_parts(unsigned data_parts); + namespace Cells { inline ParticleRange particles(Utils::Span cells) { /* Find first non-empty cell */ @@ -522,6 +532,20 @@ struct CellStructure { double range, BoxGeometry const &box, LocalBox &local_geo); + /** + * @brief Set the particle decomposition to HybridDecomposition. + * + * @param comm Communicator to use. + * @param cutoff_regular Interaction cutoff_regular. + * @param box Box geometry. + * @param local_geo Geometry of the local box. + * @param n_square_types Particle types to put into n_square decomposition. + */ + void set_hybrid_decomposition(boost::mpi::communicator const &comm, + double cutoff_regular, BoxGeometry const &box, + LocalBox &local_geo, + std::set n_square_types); + public: template void bond_loop(BondKernel const &bond_kernel) { for (auto &p : local_particles()) { diff --git a/src/core/CellStructureType.hpp b/src/core/CellStructureType.hpp index b3d7ffa56e5..17436585d58 100644 --- a/src/core/CellStructureType.hpp +++ b/src/core/CellStructureType.hpp @@ -27,7 +27,9 @@ enum class CellStructureType : int { /** cell structure regular decomposition */ CELL_STRUCTURE_REGULAR = 1, /** cell structure n square */ - CELL_STRUCTURE_NSQUARE = 2 + CELL_STRUCTURE_NSQUARE = 2, + /** cell structure hybrid */ + CELL_STRUCTURE_HYBRID = 3 }; #endif // ESPRESSO_CELLSTRUCTURETYPE_HPP diff --git a/src/core/HybridDecomposition.cpp b/src/core/HybridDecomposition.cpp new file mode 100644 index 00000000000..e5c21e4b8c1 --- /dev/null +++ b/src/core/HybridDecomposition.cpp @@ -0,0 +1,164 @@ +/* + * Copyright (C) 2010-2020 The ESPResSo project + * Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010 + * Max-Planck-Institute for Polymer Research, Theory Group + * + * This file is part of ESPResSo. + * + * ESPResSo is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * ESPResSo is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#include "HybridDecomposition.hpp" + +#include "CellStructure.hpp" +#include "event.hpp" + +#include +#include +#include + +HybridDecomposition::HybridDecomposition(boost::mpi::communicator comm, + double cutoff_regular, + const BoxGeometry &box_geo, + const LocalBox &local_box, + 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_n_square(AtomDecomposition(m_comm, m_box)), + m_n_square_types(std::move(n_square_types)) { + + /* Vector containing cells of both child decompositions */ + m_local_cells = m_regular_decomposition.get_local_cells(); + auto local_cells_n_square = m_n_square.get_local_cells(); + std::copy(local_cells_n_square.begin(), local_cells_n_square.end(), + std::back_inserter(m_local_cells)); + + /* Vector containing ghost cells of both child decompositions */ + m_ghost_cells = m_regular_decomposition.get_ghost_cells(); + auto ghost_cells_n_square = m_n_square.get_ghost_cells(); + std::copy(ghost_cells_n_square.begin(), ghost_cells_n_square.end(), + std::back_inserter(m_ghost_cells)); + + /* Communicators that contain communications of both child decompositions */ + m_exchange_ghosts_comm = m_regular_decomposition.exchange_ghosts_comm(); + auto exchange_ghosts_comm_n_square = m_n_square.exchange_ghosts_comm(); + std::copy(exchange_ghosts_comm_n_square.communications.begin(), + exchange_ghosts_comm_n_square.communications.end(), + std::back_inserter(m_exchange_ghosts_comm.communications)); + + m_collect_ghost_force_comm = + m_regular_decomposition.collect_ghost_force_comm(); + auto collect_ghost_force_comm_n_square = + m_n_square.collect_ghost_force_comm(); + std::copy(collect_ghost_force_comm_n_square.communications.begin(), + collect_ghost_force_comm_n_square.communications.end(), + std::back_inserter(m_collect_ghost_force_comm.communications)); + + /* coupling between the child decompositions via neighborship relation */ + std::vector additional_reds = m_n_square.get_local_cells(); + std::copy(ghost_cells_n_square.begin(), ghost_cells_n_square.end(), + std::back_inserter(additional_reds)); + for (auto &local_cell : m_regular_decomposition.local_cells()) { + std::vector red_neighbors(local_cell->m_neighbors.red().begin(), + local_cell->m_neighbors.red().end()); + std::vector black_neighbors(local_cell->m_neighbors.black().begin(), + local_cell->m_neighbors.black().end()); + std::copy(additional_reds.begin(), additional_reds.end(), + std::back_inserter(red_neighbors)); + local_cell->m_neighbors = Neighbors(red_neighbors, black_neighbors); + } +} + +void HybridDecomposition::resort(bool global, + std::vector &diff) { + ParticleList displaced_parts; + + /* Check for n_square type particles in regular decomposition */ + for (auto &c : m_regular_decomposition.local_cells()) { + for (auto it = c->particles().begin(); it != c->particles().end();) { + /* Particle is in the right decomposition, i.e. has no n_square type */ + if (not is_n_square_type(it->p.type)) { + std::advance(it, 1); + continue; + } + + /* else remove from current cell ... */ + auto p = std::move(*it); + it = c->particles().erase(it); + diff.emplace_back(ModifiedList{c->particles()}); + diff.emplace_back(RemovedParticle{p.identity()}); + + /* ... and insert into a n_square cell */ + auto const first_local_cell = m_n_square.get_local_cells()[0]; + first_local_cell->particles().insert(std::move(p)); + diff.emplace_back(ModifiedList{first_local_cell->particles()}); + } + + /* Now check for regular decomposition type particles in n_square */ + for (auto &c : m_n_square.local_cells()) { + for (auto it = c->particles().begin(); it != c->particles().end();) { + /* Particle is of n_square type */ + if (is_n_square_type(it->p.type)) { + std::advance(it, 1); + continue; + } + + /* else remove from current cell ... */ + auto p = std::move(*it); + it = c->particles().erase(it); + diff.emplace_back(ModifiedList{c->particles()}); + diff.emplace_back(RemovedParticle{p.identity()}); + + /* ... and insert in regular decomposition */ + auto const target_cell = particle_to_cell(p); + /* if particle belongs to this node insert it into correct cell */ + if (target_cell != nullptr) { + target_cell->particles().insert(std::move(p)); + diff.emplace_back(ModifiedList{target_cell->particles()}); + } + /* otherwise just put into regular decomposition */ + else { + auto first_local_cell = m_regular_decomposition.get_local_cells()[0]; + first_local_cell->particles().insert(std::move(p)); + diff.emplace_back(ModifiedList{first_local_cell->particles()}); + } + } + } + } + + /* now resort into correct cells within the respective decompositions */ + m_regular_decomposition.resort(global, diff); + m_n_square.resort(global, diff); + + /* basically do CellStructure::ghost_count() */ + ghost_communicator(exchange_ghosts_comm(), GHOSTTRANS_PARTNUM); + + /* basically do CellStructure::ghost_update(unsigned data_parts) */ + ghost_communicator(exchange_ghosts_comm(), + map_data_parts(global_ghost_flags())); +} + +Utils::Vector +HybridDecomposition::parts_per_decomposition_local() const { + std::size_t parts_in_regular = 0; + std::size_t parts_in_n_square = 0; + for (auto &c : m_regular_decomposition.get_local_cells()) { + parts_in_regular += c->particles().size(); + } + for (auto &c : m_n_square.get_local_cells()) { + parts_in_n_square += c->particles().size(); + } + return {parts_in_regular, parts_in_n_square}; +} diff --git a/src/core/HybridDecomposition.hpp b/src/core/HybridDecomposition.hpp new file mode 100644 index 00000000000..53eff6acdbd --- /dev/null +++ b/src/core/HybridDecomposition.hpp @@ -0,0 +1,125 @@ +/* + * Copyright (C) 2010-2020 The ESPResSo project + * Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010 + * Max-Planck-Institute for Polymer Research, Theory Group + * + * This file is part of ESPResSo. + * + * ESPResSo is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * ESPResSo is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#ifndef ESPRESSO_HYBRID_DECOMPOSITION_HPP +#define ESPRESSO_HYBRID_DECOMPOSITION_HPP + +#include "AtomDecomposition.hpp" +#include "Cell.hpp" +#include "Particle.hpp" +#include "ParticleDecomposition.hpp" +#include "RegularDecomposition.hpp" +#include "ghosts.hpp" +#include "integrate.hpp" + +#include +#include + +#include +#include + +#include +#include +#include + +namespace Utils { +using Vector2i = Vector; +} + +class HybridDecomposition : public ParticleDecomposition { + boost::mpi::communicator m_comm; + BoxGeometry m_box; + double m_cutoff_regular; + std::vector m_local_cells; + std::vector m_ghost_cells; + + GhostCommunicator m_exchange_ghosts_comm; + GhostCommunicator m_collect_ghost_force_comm; + + /** RegularDecomposition to hold the small particles */ + RegularDecomposition m_regular_decomposition; + /** N-Square Decomposition to hold large particles */ + AtomDecomposition m_n_square; + /** Set containing the types that should be handled using n_square */ + std::set const m_n_square_types; + + bool is_n_square_type(int type_id) { + return (m_n_square_types.find(type_id) != m_n_square_types.end()); + } + +public: + HybridDecomposition(boost::mpi::communicator comm, double cutoff_regular, + const BoxGeometry &box_geo, + const LocalBox &local_box, + std::set n_square_types); + + Utils::Vector3i get_cell_grid() const { + return m_regular_decomposition.cell_grid; + } + + Utils::Vector3d get_cell_size() const { + return m_regular_decomposition.cell_size; + } + + std::set get_n_square_types() const { return m_n_square_types; } + + void resort(bool global, std::vector &diff) override; + + double get_cutoff_regular() const { return m_cutoff_regular; } + + GhostCommunicator const &exchange_ghosts_comm() const override { + return m_exchange_ghosts_comm; + } + + GhostCommunicator const &collect_ghost_force_comm() const override { + return m_collect_ghost_force_comm; + } + + Utils::Span local_cells() override { + return Utils::make_span(m_local_cells); + } + + Utils::Span ghost_cells() override { + return Utils::make_span(m_ghost_cells); + } + + Cell *particle_to_cell(Particle const &p) override { + if (is_n_square_type(p.p.type)) { + return m_n_square.particle_to_cell(p); + } + return m_regular_decomposition.particle_to_cell(p); + } + + Utils::Vector3d max_cutoff() const override { + return m_n_square.max_cutoff(); + } + Utils::Vector3d max_range() const override { return m_n_square.max_range(); } + + boost::optional minimum_image_distance() const override { + return m_box; + } + + BoxGeometry box() const override { return m_box; } + + Utils::Vector parts_per_decomposition_local() const; +}; + +#endif // ESPRESSO_HYBRID_DECOMPOSITION_HPP diff --git a/src/core/RegularDecomposition.hpp b/src/core/RegularDecomposition.hpp index 8713f58788f..e607637f3a5 100644 --- a/src/core/RegularDecomposition.hpp +++ b/src/core/RegularDecomposition.hpp @@ -105,6 +105,10 @@ struct RegularDecomposition : public ParticleDecomposition { return Utils::make_span(m_ghost_cells); } + /* Getter needed for HybridDecomposition */ + std::vector get_local_cells() const { return m_local_cells; } + std::vector get_ghost_cells() const { return m_ghost_cells; } + Cell *particle_to_cell(Particle const &p) override { return position_to_cell(p.pos()); } diff --git a/src/core/cells.cpp b/src/core/cells.cpp index cc76b46f406..f3160de3082 100644 --- a/src/core/cells.cpp +++ b/src/core/cells.cpp @@ -43,6 +43,7 @@ #include #include +#include #include #include #include @@ -195,6 +196,61 @@ std::vector mpi_resort_particles(int global_flag) { return n_parts; } +void set_regular_decomposition() { + cell_structure.set_regular_decomposition(comm_cart, interaction_range(), + box_geo, local_geo); + on_cell_structure_change(); +} + +void set_atom_decomposition() { + cell_structure.set_atom_decomposition(comm_cart, box_geo, local_geo); + on_cell_structure_change(); +} + +void set_hybrid_decomposition(std::set n_square_types, + double cutoff_regular) { + cell_structure.set_hybrid_decomposition(comm_cart, cutoff_regular, box_geo, + local_geo, n_square_types); + on_cell_structure_change(); +} + +REGISTER_CALLBACK(set_regular_decomposition) +REGISTER_CALLBACK(set_atom_decomposition) +REGISTER_CALLBACK(set_hybrid_decomposition) + +void mpi_set_regular_decomposition() { + mpi_call_all(set_regular_decomposition); +} +void mpi_set_atom_decomposition() { mpi_call_all(set_atom_decomposition); } +void mpi_set_hybrid_decomposition(std::set n_square_types, + double cutoff_regular) { + mpi_call_all(set_hybrid_decomposition, n_square_types, cutoff_regular); +} + +static Utils::Vector hybrid_parts_per_decomposition_local() { + assert(cell_structure.decomposition_type() == + CellStructureType::CELL_STRUCTURE_HYBRID); + /* Get current HybridDecomposition to extract n_square_types */ + auto ¤t_hybrid_decomposition = + dynamic_cast( + Utils::as_const(cell_structure).decomposition()); + + Utils::Vector parts_per_decomposition; + boost::mpi::reduce( + comm_cart, current_hybrid_decomposition.parts_per_decomposition_local(), + parts_per_decomposition, std::plus<>{}, 0); + + return parts_per_decomposition; +} + +REGISTER_CALLBACK_MAIN_RANK(hybrid_parts_per_decomposition_local) + +std::pair hybrid_parts_per_decomposition() { + auto const parts_per_decomposition = mpi_call( + Communication::Result::main_rank, hybrid_parts_per_decomposition_local); + return std::make_pair(parts_per_decomposition[0], parts_per_decomposition[1]); +} + void cells_re_init(CellStructureType new_cs) { switch (new_cs) { case CellStructureType::CELL_STRUCTURE_REGULAR: @@ -204,6 +260,16 @@ void cells_re_init(CellStructureType new_cs) { case CellStructureType::CELL_STRUCTURE_NSQUARE: cell_structure.set_atom_decomposition(comm_cart, box_geo, local_geo); break; + case CellStructureType::CELL_STRUCTURE_HYBRID: { + /* Get current HybridDecomposition to extract n_square_types */ + auto ¤t_hybrid_decomposition = + dynamic_cast( + Utils::as_const(cell_structure).decomposition()); + cell_structure.set_hybrid_decomposition( + comm_cart, current_hybrid_decomposition.get_cutoff_regular(), box_geo, + local_geo, current_hybrid_decomposition.get_n_square_types()); + break; + } default: throw std::runtime_error("Unknown cell system type"); } @@ -270,6 +336,11 @@ const RegularDecomposition *get_regular_decomposition() { Utils::as_const(cell_structure).decomposition()); } +const HybridDecomposition *get_hybrid_decomposition() { + return &dynamic_cast( + Utils::as_const(cell_structure).decomposition()); +} + void mpi_set_use_verlet_lists_local(bool use_verlet_lists) { cell_structure.use_verlet_list = use_verlet_lists; } diff --git a/src/core/cells.hpp b/src/core/cells.hpp index 52b708c3a6a..132049b1834 100644 --- a/src/core/cells.hpp +++ b/src/core/cells.hpp @@ -36,6 +36,14 @@ * regardless their spatial position (see \ref AtomDecomposition.hpp). * This is suitable for long range interactions that cannot be treated by a * special method like P3M. + * - hybrid decomposition: Initializes both regular decomposition and nsquare + * at the same time and and has is given a set of particle types + * n_square_types (see \ref HybridDecomposition.hpp). By default, particles + * will be distributed using the regular decomposition. For particles of the + * types defined as n_square_types the nsquare method is used. This is + * suitable for systems containing lots of small particles with short + * interaction range mixed with a few large particles with long interaction + * range. There, the large particles should be treated using nsquare. */ #include "Cell.hpp" @@ -61,6 +69,39 @@ enum { /** Type of cell structure in use. */ extern CellStructure cell_structure; +/** Initialize cell sturcture DomainDecomposition */ +void set_regular_decomposition(); + +/** Initialize cell structure AtomDecomposition */ +void set_atom_decomposition(); + +/** Initialize cell structure HybridDecomposition + * @param n_square_types Set of particle types that will use n_square. + */ +void set_hybrid_decomposition(std::set n_square_types, + double cutoff_regular); + +/** Change cell structure to DomainDecomposition on all nodes. */ +void mpi_set_regular_decomposition(); + +/** Change cell structure to AtomDecomposition on all nodes. */ +void mpi_set_atom_decomposition(); + +/** Change cell structure to HybridDecomposition on all nodes. + * @param n_square_types Set of particle types that will use n_square. + */ +void mpi_set_hybrid_decomposition(std::set n_square_types, + double cutoff_regular); + +/** + * @brief Number of parts in child decompositions of a HybridDecomposition. + * + * Throws if called when active decomposition type is not hybrid. + * + * @return pair(parts_in_domain_dec, parts_in_n_square) + */ +std::pair hybrid_parts_per_decomposition(); + /** Reinitialize the cell structures. * @param new_cs The new topology to use afterwards. */ @@ -128,6 +169,13 @@ Cell *find_current_cell(const Particle &p); */ const RegularDecomposition *get_regular_decomposition(); +/** + * @brief Return a pointer to the global HybridDecomposition. + * + * @return Pointer to the decomposition if it is set, nullptr otherwise. + */ +const HybridDecomposition *get_hybrid_decomposition(); + class PairInfo { public: PairInfo() = default; diff --git a/src/core/event.cpp b/src/core/event.cpp index 9acc8934515..94a276e87a0 100644 --- a/src/core/event.cpp +++ b/src/core/event.cpp @@ -193,7 +193,12 @@ void on_particle_charge_change() { } void on_particle_change() { - cell_structure.set_resort_particles(Cells::RESORT_LOCAL); + if (cell_structure.decomposition_type() == + CellStructureType::CELL_STRUCTURE_HYBRID) { + cell_structure.set_resort_particles(Cells::RESORT_GLOBAL); + } else { + cell_structure.set_resort_particles(Cells::RESORT_LOCAL); + } reinit_electrostatics = true; reinit_magnetostatics = true; recalc_forces = true; diff --git a/src/core/grid_based_algorithms/lb.cpp b/src/core/grid_based_algorithms/lb.cpp index 3401746e2d0..c3388e57b9b 100644 --- a/src/core/grid_based_algorithms/lb.cpp +++ b/src/core/grid_based_algorithms/lb.cpp @@ -614,8 +614,8 @@ void lb_sanity_checks(const LB_Parameters &lb_parameters) { if (lb_parameters.viscosity <= 0.0) { runtimeErrorMsg() << "Lattice Boltzmann fluid viscosity not set"; } - if (local_geo.cell_structure_type() == - CellStructureType::CELL_STRUCTURE_NSQUARE) { + if (local_geo.cell_structure_type() != + CellStructureType::CELL_STRUCTURE_REGULAR) { if (n_nodes > 1) { runtimeErrorMsg() << "LB only works with regular decomposition, " "if using more than one MPI node"; diff --git a/src/python/espressomd/cellsystem.pxd b/src/python/espressomd/cellsystem.pxd index c734e55d1dd..3ff9080dc86 100644 --- a/src/python/espressomd/cellsystem.pxd +++ b/src/python/espressomd/cellsystem.pxd @@ -18,8 +18,9 @@ # from libcpp cimport bool -from libcpp.vector cimport vector +from libcpp.set cimport set as cpp_set from libcpp.pair cimport pair, tuple +from libcpp.vector cimport vector from .utils cimport Vector3i, Vector3d cdef extern from "cells.hpp": @@ -38,6 +39,7 @@ cdef extern from "CellStructureType.hpp": ctypedef enum CellStructureType: CELL_STRUCTURE_REGULAR "CellStructureType::CELL_STRUCTURE_REGULAR" CELL_STRUCTURE_NSQUARE "CellStructureType::CELL_STRUCTURE_NSQUARE" + CELL_STRUCTURE_HYBRID "CellStructureType::CELL_STRUCTURE_HYBRID" cdef extern from "cells.hpp": ctypedef struct CellStructure: @@ -47,14 +49,18 @@ cdef extern from "cells.hpp": CellStructure cell_structure const RegularDecomposition * get_regular_decomposition() + const HybridDecomposition * get_hybrid_decomposition() vector[pair[int, int]] mpi_get_pairs(double distance) except + vector[pair[int, int]] mpi_get_pairs_of_types(double distance, vector[int] types) except + vector[PairInfo] mpi_non_bonded_loop_trace() vector[int] mpi_resort_particles(int global_flag) void mpi_bcast_cell_structure(int cs) + void mpi_set_hybrid_decomposition(cpp_set[int] n_square_types, double cutoff_regular) void mpi_set_use_verlet_lists(bool use_verlet_lists) + pair[size_t, size_t] hybrid_parts_per_decomposition() + cdef extern from "tuning.hpp": cdef void c_tune_skin "tune_skin" (double min_skin, double max_skin, double tol, int int_steps, bool adjust_max_skin) @@ -68,6 +74,13 @@ cdef extern from "RegularDecomposition.hpp": Vector3i cell_grid double cell_size[3] +cdef extern from "HybridDecomposition.hpp": + cppclass HybridDecomposition: + Vector3i get_cell_grid() + Vector3d get_cell_size() + cpp_set[int] get_n_square_types() + double get_cutoff_regular() + cdef extern from "grid.hpp": void mpi_set_node_grid(const Vector3i & node_grid) diff --git a/src/python/espressomd/cellsystem.pyx b/src/python/espressomd/cellsystem.pyx index 8a76ea2f759..e668a6daeb9 100644 --- a/src/python/espressomd/cellsystem.pyx +++ b/src/python/espressomd/cellsystem.pyx @@ -63,6 +63,31 @@ cdef class CellSystem: return True + def set_hybrid_decomposition( + self, n_square_types=None, + cutoff_regular=None, use_verlet_lists=True): + """ + Activates the hybrid domain decomposition. + + Parameters + ---------- + n_square_types : list of :obj:`int` + Particles types that should be handled in nsquare cell system. + cutoff_regular : :obj:`float` + Maximum cutoff to consider in regular decomposition. + Should be as low as the system permits. + use_verlet_lists : :obj:`bool`, optional + Activates or deactivates the usage of the Verlet + lists for this algorithm. + + """ + mpi_set_use_verlet_lists(use_verlet_lists) + if n_square_types is None: + n_square_types = set() + mpi_set_hybrid_decomposition(n_square_types, cutoff_regular) + + return True + def get_state(self): s = self.__getstate__() @@ -72,6 +97,22 @@ cdef class CellSystem: [rd.cell_grid[0], rd.cell_grid[1], rd.cell_grid[2]]) s["cell_size"] = np.array( [rd.cell_size[0], rd.cell_size[1], rd.cell_size[2]]) + elif cell_structure.decomposition_type() == CellStructureType.CELL_STRUCTURE_HYBRID: + hd = get_hybrid_decomposition() + cell_grid = hd.get_cell_grid() + cell_size = hd.get_cell_size() + s["cell_grid"] = np.array([cell_grid[d] for d in range(3)]) + s["cell_size"] = np.array([cell_size[d] for d in range(3)]) + s["n_square_types"] = hd.get_n_square_types() + s["cutoff_regular"] = hd.get_cutoff_regular() + # manually trigger resort to make sure that particles end up + # in their respective child decomposition before querying + self.resort() + n_regular, n_n_square = hybrid_parts_per_decomposition() + s["parts_per_decomposition"] = { + "regular": n_regular, + "n_square": n_n_square + } s["verlet_reuse"] = get_verlet_reuse() s["n_nodes"] = n_nodes @@ -85,6 +126,8 @@ cdef class CellSystem: s["type"] = "regular_decomposition" if cell_structure.decomposition_type() == CellStructureType.CELL_STRUCTURE_NSQUARE: s["type"] = "nsquare" + if cell_structure.decomposition_type() == CellStructureType.CELL_STRUCTURE_HYBRID: + s["type"] = "hybrid_decomposition" s["skin"] = skin s["node_grid"] = np.array([node_grid[0], node_grid[1], node_grid[2]]) @@ -99,6 +142,9 @@ cdef class CellSystem: use_verlet_lists=d['use_verlet_list']) elif d['type'] == "nsquare": self.set_n_square(use_verlet_lists=d['use_verlet_list']) + elif d['type'] == "hybrid_decomposition": + self.set_hybrid_decomposition( + use_verlet_lists=d['use_verlet_list']) def get_pairs(self, distance, types='all'): """ diff --git a/testsuite/python/CMakeLists.txt b/testsuite/python/CMakeLists.txt index e70cae5c15c..4e75b73ff56 100644 --- a/testsuite/python/CMakeLists.txt +++ b/testsuite/python/CMakeLists.txt @@ -121,6 +121,8 @@ python_test(FILE dds-and-bh-gpu.py MAX_NUM_PROC 4 LABELS gpu) python_test(FILE electrostaticInteractions.py MAX_NUM_PROC 2) python_test(FILE engine_langevin.py MAX_NUM_PROC 4) python_test(FILE engine_lb.py MAX_NUM_PROC 2 LABELS gpu) +python_test(FILE engine_lb_hybrid.py MAX_NUM_PROC 2 LABELS gpu) +python_test(FILE engine_lb_hybrid.py MAX_NUM_PROC 1 LABELS gpu SUFFIX 1_core) python_test(FILE engine_lb_nsquare.py MAX_NUM_PROC 1 LABELS gpu) python_test(FILE icc.py MAX_NUM_PROC 4) python_test(FILE mass-and-rinertia_per_particle.py MAX_NUM_PROC 2 LABELS long) @@ -175,6 +177,8 @@ python_test(FILE virtual_sites_tracers.py MAX_NUM_PROC 2 DEPENDENCIES python_test(FILE virtual_sites_tracers_gpu.py MAX_NUM_PROC 2 LABELS gpu DEPENDENCIES virtual_sites_tracers_common.py) python_test(FILE regular_decomposition.py MAX_NUM_PROC 4) +python_test(FILE hybrid_decomposition.py MAX_NUM_PROC 1 SUFFIX 1_core) +python_test(FILE hybrid_decomposition.py MAX_NUM_PROC 4) python_test(FILE integrator_npt.py MAX_NUM_PROC 4) python_test(FILE integrator_npt_stats.py MAX_NUM_PROC 4 LABELS long) python_test(FILE integrator_steepest_descent.py MAX_NUM_PROC 4) @@ -186,6 +190,9 @@ python_test(FILE dipolar_interface.py MAX_NUM_PROC 1 LABELS gpu) python_test(FILE dipolar_mpi_exceptions.py MAX_NUM_PROC 2) python_test(FILE lb.py MAX_NUM_PROC 2 LABELS gpu) python_test(FILE lb_stats.py MAX_NUM_PROC 2 LABELS gpu long) +python_test(FILE lb_stats_hybrid.py MAX_NUM_PROC 2 LABELS gpu long) +python_test(FILE lb_stats_hybrid.py MAX_NUM_PROC 1 LABELS gpu long SUFFIX + 1_core) python_test(FILE lb_stats_nsquare.py MAX_NUM_PROC 1 LABELS gpu long) python_test(FILE lb_vtk.py MAX_NUM_PROC 2 LABELS gpu) python_test(FILE force_cap.py MAX_NUM_PROC 2) @@ -246,6 +253,9 @@ python_test(FILE lb_shear.py MAX_NUM_PROC 2 LABELS gpu) python_test(FILE lb_thermostat.py MAX_NUM_PROC 2 LABELS gpu) python_test(FILE lb_buoyancy_force.py MAX_NUM_PROC 4 LABELS gpu) python_test(FILE lb_momentum_conservation.py MAX_NUM_PROC 4 LABELS gpu) +python_test(FILE lb_momentum_conservation_hybrid.py MAX_NUM_PROC 4 LABELS gpu) +python_test(FILE lb_momentum_conservation_hybrid.py MAX_NUM_PROC 1 LABELS gpu + SUFFIX 1_core) python_test(FILE lb_momentum_conservation_nsquare.py MAX_NUM_PROC 1 LABELS gpu) python_test(FILE p3m_electrostatic_pressure.py MAX_NUM_PROC 2) python_test(FILE sigint.py DEPENDENCIES sigint_child.py MAX_NUM_PROC 1) diff --git a/testsuite/python/cellsystem.py b/testsuite/python/cellsystem.py index c7f4b72070a..c0868f2710e 100644 --- a/testsuite/python/cellsystem.py +++ b/testsuite/python/cellsystem.py @@ -35,10 +35,16 @@ def test_cell_system(self): s = self.system.cell_system.get_state() self.assertEqual( [s['use_verlet_list'], s['type']], [1, "regular_decomposition"]) + self.system.cell_system.set_hybrid_decomposition( + use_verlet_lists=False, n_square_types={1, 3, 5}, cutoff_regular=1.27) + s = self.system.cell_system.get_state() + self.assertEqual( + [s['use_verlet_list'], s['type'], + s['n_square_types'], s['cutoff_regular']], + [0, "hybrid_decomposition", {1, 3, 5}, 1.27]) @ut.skipIf(n_nodes == 1, "Skipping test: only runs for n_nodes >= 2") - def test_node_grid(self): - self.system.cell_system.set_regular_decomposition() + def check_node_grid(self): for i in range(3): node_grid_ref = [1, 1, 1] node_grid_ref[i] = self.n_nodes @@ -46,6 +52,15 @@ def test_node_grid(self): node_grid = self.system.cell_system.get_state()['node_grid'] np.testing.assert_array_equal(node_grid, node_grid_ref) + def test_node_grid_regular(self): + self.system.cell_system.set_regular_decomposition() + self.check_node_grid() + + def test_node_grid_hybrid(self): + self.system.cell_system.set_hybrid_decomposition( + n_square_types={1}, cutoff_regular=0) + self.check_node_grid() + if __name__ == "__main__": ut.main() diff --git a/testsuite/python/engine_lb.py b/testsuite/python/engine_lb.py index 24acf55aa5e..4a6ca91b369 100644 --- a/testsuite/python/engine_lb.py +++ b/testsuite/python/engine_lb.py @@ -53,22 +53,23 @@ def add_all_types_of_swimmers( pos2 = [2, 3, 5.99] if put_in_corners else [2.9, 2.5, 3] pos3 = [1.5, 5.99, 1] if put_in_corners else [2, 2, 2.5] + # particle type is only relevant for engine_lb_hybrid test system.part.add(pos=pos0, quat=minus_y, fix=3 * [fix], mass=0.9, rinertia=3 * [7], rotation=3 * [rotation], swimming={"mode": "pusher", "f_swim": 0.10, - "dipole_length": 0.5}) + "dipole_length": 0.5}, type=1) system.part.add(pos=pos1, quat=plus_x, fix=3 * [fix], mass=1.9, rinertia=3 * [8], rotation=3 * [rotation], swimming={"mode": "pusher", "v_swim": 0.02, - "dipole_length": 0.6}) + "dipole_length": 0.6}, type=0) system.part.add(pos=pos2, quat=plus_z, fix=3 * [fix], mass=2.9, rinertia=3 * [9], rotation=3 * [rotation], swimming={"mode": "puller", "f_swim": 0.08, - "dipole_length": 0.7}) + "dipole_length": 0.7}, type=1) system.part.add(pos=pos3, quat=plus_y, fix=3 * [fix], mass=3.9, rinertia=3 * [10], rotation=3 * [rotation], swimming={"mode": "puller", "v_swim": 0.05, - "dipole_length": 0.8}) + "dipole_length": 0.8}, type=0) def setUp(self): self.set_cellsystem() diff --git a/testsuite/python/engine_lb_hybrid.py b/testsuite/python/engine_lb_hybrid.py new file mode 100644 index 00000000000..605b60bcfb8 --- /dev/null +++ b/testsuite/python/engine_lb_hybrid.py @@ -0,0 +1,77 @@ +# Copyright (C) 2010-2022 The ESPResSo project +# +# This file is part of ESPResSo. +# +# ESPResSo is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# ESPResSo is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . + +import espressomd +import espressomd.lb +import unittest as ut +import unittest_decorators as utx +from engine_lb import SwimmerTest + + +@utx.skipIfMissingFeatures(["ENGINE", "ROTATIONAL_INERTIA", "MASS"]) +@ut.skipIf(SwimmerTest.n_nodes > 1, + "LB with N-square only works on 1 MPI rank") +class SwimmerTestHybrid0CPU(SwimmerTest, ut.TestCase): + + lb_class = espressomd.lb.LBFluid + tol = 1e-10 + + def set_cellsystem(self): + self.system.cell_system.set_hybrid_decomposition( + n_square_types={0}, cutoff_regular=1) + + +@utx.skipIfMissingFeatures(["ENGINE", "ROTATIONAL_INERTIA", "MASS"]) +@ut.skipIf(SwimmerTest.n_nodes > 1, + "LB with N-square only works on 1 MPI rank") +class SwimmerTestHybrid1CPU(SwimmerTest, ut.TestCase): + + lb_class = espressomd.lb.LBFluid + tol = 1e-10 + + def set_cellsystem(self): + self.system.cell_system.set_hybrid_decomposition( + n_square_types={1}, cutoff_regular=1) + + +@utx.skipIfMissingGPU() +@utx.skipIfMissingFeatures(["ENGINE", "ROTATIONAL_INERTIA", "MASS"]) +class SwimmerTestHybrid0GPU(SwimmerTest, ut.TestCase): + + lb_class = espressomd.lb.LBFluidGPU + tol = 1e-5 + + def set_cellsystem(self): + self.system.cell_system.set_hybrid_decomposition( + n_square_types={0}, cutoff_regular=1) + + +@utx.skipIfMissingGPU() +@utx.skipIfMissingFeatures(["ENGINE", "ROTATIONAL_INERTIA", "MASS"]) +class SwimmerTestHybrid1GPU(SwimmerTest, ut.TestCase): + + lb_class = espressomd.lb.LBFluidGPU + tol = 1e-5 + + def set_cellsystem(self): + self.system.cell_system.set_hybrid_decomposition( + n_square_types={1}, cutoff_regular=1) + print(self.system.cell_system.get_state()) + + +if __name__ == "__main__": + ut.main() diff --git a/testsuite/python/engine_lb_nsquare.py b/testsuite/python/engine_lb_nsquare.py index 0d5415d6c8d..44a3e44cf7e 100644 --- a/testsuite/python/engine_lb_nsquare.py +++ b/testsuite/python/engine_lb_nsquare.py @@ -36,8 +36,6 @@ def set_cellsystem(self): @utx.skipIfMissingGPU() @utx.skipIfMissingFeatures(["ENGINE", "ROTATIONAL_INERTIA", "MASS"]) -@ut.skipIf(SwimmerTest.n_nodes > 1, - "LB with N-square only works on 1 MPI rank") class SwimmerTestNSquareGPU(SwimmerTest, ut.TestCase): lb_class = espressomd.lb.LBFluidGPU diff --git a/testsuite/python/hybrid_decomposition.py b/testsuite/python/hybrid_decomposition.py new file mode 100644 index 00000000000..56448c2d17f --- /dev/null +++ b/testsuite/python/hybrid_decomposition.py @@ -0,0 +1,258 @@ +# +# Copyright (C) 2013-2019 The ESPResSo project +# +# This file is part of ESPResSo. +# +# ESPResSo is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# ESPResSo is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# +import unittest as ut +import unittest_decorators as utx +import espressomd +import itertools +import numpy as np + + +class HybridDecomposition(ut.TestCase): + system = espressomd.System(box_l=3 * [50.0]) + original_node_grid = tuple(system.cell_system.node_grid) + + def setUp(self): + self.system.cell_system.set_hybrid_decomposition( + use_verlet_lists=False, n_square_types={1}, cutoff_regular=2.5) + self.system.cell_system.node_grid = self.original_node_grid + self.system.time_step = 1e-3 + + def tearDown(self): + self.system.part.clear() + + def check_resort(self): + n_part = 2352 + + # number of particles that should end up in the respective child + # decomposition + n_n_square = n_part // 4 + n_regular = n_part - n_n_square + + # Add the particles on node 0, so that they have to be resorted + particles = self.system.part.add( + pos=n_part * [(0, 0, 0)], type=n_n_square * [0, 0, 0, 1]) + parts_per_decomposition = self.system.cell_system.get_state()[ + 'parts_per_decomposition'] + self.assertEqual(parts_per_decomposition['n_square'], n_n_square) + self.assertEqual(parts_per_decomposition['regular'], n_regular) + + # And now change their positions + particles.pos = self.system.box_l * \ + np.random.random((n_part, 3)) + + # Add an interacting particle in a corner of the box + self.system.part.add(pos=(0.01, 0.01, 0.01), type=0) + if espressomd.has_features(['LENNARD_JONES']): + self.system.non_bonded_inter[0, 1].lennard_jones.set_params( + epsilon=1.0, sigma=3.0, cutoff=6.0, shift=0.1) + ref_energy = self.system.analysis.energy()['total'] + assert ref_energy > 10. + + # Distribute the particles on the nodes + part_dist = self.system.cell_system.resort() + + # Check that we did not lose particles + self.assertEqual(sum(part_dist), n_part + 1) + + # Check that we can still access all the particles + # This basically checks if part_node and local_particles + # is still in a valid state after the particle exchange + self.assertEqual(len(self.system.part.all()), n_part + 1) + self.assertEqual(sum(self.system.part.all().type), n_n_square) + + # Check that the system is still valid + if espressomd.has_features(['LENNARD_JONES']): + # energy calculation + new_energy = self.system.analysis.energy()['total'] + self.assertEqual(new_energy, ref_energy) + # force calculation + self.system.integrator.run(0, recalc_forces=True) + + def prepare_hybrid_setup(self, n_part_small=0, n_part_large=0): + """Setup system with small and large particles, minimize + energy and setup thermostat. Particles have random + initial velocities; cutoff of small particles is 2.5. + + """ + box_l = self.system.box_l[0] + if n_part_small > 0: + self.system.part.add( + type=[0] * n_part_small, + pos=np.random.random( + (n_part_small, + 3)) * box_l, + v=np.random.randn( + n_part_small, + 3)) + if n_part_large > 0: + self.system.part.add( + type=[1] * n_part_large, + pos=np.random.random( + (n_part_large, + 3)) * box_l, + v=np.random.randn( + n_part_large, + 3)) + self.assertEqual(len(self.system.part.all()), + n_part_small + n_part_large) + + # setup interactions + self.system.non_bonded_inter[0, 0].lennard_jones.set_params( + epsilon=1, sigma=1, cutoff=2.5, shift="auto") + # mixing rule: sigma = 0.5 * (s_large + s_small) + self.system.non_bonded_inter[0, 1].lennard_jones.set_params( + epsilon=1, sigma=2.5, cutoff=2**(1 / 6) * 2.5, shift="auto") + self.system.non_bonded_inter[1, 1].lennard_jones.set_params( + epsilon=1, sigma=4, cutoff=10, shift="auto") + + # remove overlap + self.system.integrator.set_steepest_descent( + f_max=0, gamma=30, max_displacement=0.01) + self.system.integrator.run(0) + old_force = np.max(np.linalg.norm(self.system.part.all().f, axis=1)) + while old_force > n_part_small + n_part_large: + self.system.integrator.run(20) + force = np.max(np.linalg.norm(self.system.part.all().f, axis=1)) + old_force = force + + self.system.integrator.set_vv() + + def test_resort(self): + self.check_resort() + + @ut.skipIf(system.cell_system.get_state()["n_nodes"] != 4, + "Skipping test: only runs for n_nodes >= 4") + def test_resort_alternating(self): + # check particle resorting when the left and right cells are different + self.system.cell_system.node_grid = [4, 1, 1] + self.check_resort() + + def test_position_rounding(self): + """This places a particle on the box boundary, + with parameters that could cause problems with + rounding.""" + self.system.box_l = [50.0, 50.0, 50.0] + self.system.cell_system.skin = 0.4 + self.system.min_global_cut = 12.0 / 4.25 + self.system.part.add(pos=[25, 25, 0]) + self.assertEqual(1, len(self.system.part)) + + @utx.skipIfMissingFeatures(["LENNARD_JONES"]) + def test_non_bonded_loop_trace(self): + """Validates that the distances used by the non-bonded loop + match with the minimum image distance accessible by Python, + checks that no pairs are lost or double-counted. + + """ + self.prepare_hybrid_setup(n_part_small=50, n_part_large=50) + # regular cutoff used in preparation + cutoff_regular = 2.5 + + cs_pairs = self.system.cell_system.non_bonded_loop_trace() + distance_vec = self.system.distance_vec + + # Distance for all pairs of particles obtained by Python + py_distances = {} + for p1, p2 in self.system.part.pairs(): + py_distances[p1.id, p2.id] = np.copy(distance_vec(p1, p2)) + + # Go through pairs found by the non-bonded loop and check distance + for p in cs_pairs: + # p is a tuple with (id1, id2, pos1, pos2, vec21, mpi_node) + # Note that system.distance_vec uses the opposite sign convention + # as the minimum image distance in the core + if (p[0], p[1]) in py_distances: + np.testing.assert_allclose( + np.copy(p[4]), -py_distances[p[0], p[1]]) + del py_distances[p[0], p[1]] + elif (p[1], p[0]) in py_distances: + np.testing.assert_allclose( + np.copy(p[4]), py_distances[p[1], p[0]]) + del py_distances[p[1], p[0]] + else: + raise Exception("Extra pair from core", p) + + # test for pairs from regular child decomposition with more than one + # cell + for ids, dist in py_distances.items(): + if np.linalg.norm(dist) < cutoff_regular: + raise Exception("Pair not found by the core", ids) + + @utx.skipIfMissingFeatures(["LENNARD_JONES"]) + def test_against_nsquare(self): + self.prepare_hybrid_setup(n_part_small=150, n_part_large=50) + + steps_per_round = 20 + for _ in range(40): + # integrate using hybrid and calculate energy and forces + self.system.cell_system.set_hybrid_decomposition( + n_square_types={1}, cutoff_regular=2.5) + self.system.integrator.run(steps_per_round) + energy = self.system.analysis.energy() + forces = self.system.part.all().f + + # compare to n_square for consistency + self.system.cell_system.set_n_square() + self.system.integrator.run(0) + + energy_n_square = self.system.analysis.energy() + forces_n_square = self.system.part.all().f + self.assertAlmostEqual( + energy_n_square["non_bonded"], + energy["non_bonded"]) + self.assertAlmostEqual( + energy_n_square[("non_bonded", 0, 0)], energy[("non_bonded", 0, 0)]) + self.assertAlmostEqual( + energy_n_square[("non_bonded", 0, 1)], energy[("non_bonded", 0, 1)]) + self.assertAlmostEqual( + energy_n_square[("non_bonded", 1, 1)], energy[("non_bonded", 1, 1)]) + self.assertTrue(np.allclose(forces_n_square, forces)) + + def test_sort_into_child_decs(self): + """Assert that particles end up in the respective child + decomposition, depending on their type. Also, check that + changing particle type or n_square_types results in the + expected resort. + + """ + n_parts = 3 + parts = self.system.part.add( + pos=np.random.random( + (n_parts, + 3)) * self.system.box_l[0], + type=np.random.randint( + 2, + size=n_parts)) + for ndx, types in enumerate(itertools.product([0, 1], repeat=n_parts)): + parts.type = types + n_square_type = ndx % 2 + n_n_square = n_parts - \ + np.sum(types) if n_square_type == 0 else np.sum(types) + n_regular = n_parts - n_n_square + + self.system.cell_system.set_hybrid_decomposition( + n_square_types={n_square_type}, cutoff_regular=0) + parts_per_decomposition = self.system.cell_system.get_state()[ + 'parts_per_decomposition'] + self.assertEqual(parts_per_decomposition['n_square'], n_n_square) + self.assertEqual(parts_per_decomposition['regular'], n_regular) + + +if __name__ == "__main__": + ut.main() diff --git a/testsuite/python/lb_momentum_conservation_hybrid.py b/testsuite/python/lb_momentum_conservation_hybrid.py new file mode 100644 index 00000000000..6c2bf2cb7e3 --- /dev/null +++ b/testsuite/python/lb_momentum_conservation_hybrid.py @@ -0,0 +1,72 @@ +# Copyright (C) 2010-2022 The ESPResSo project +# +# This file is part of ESPResSo. +# +# ESPResSo is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# ESPResSo is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . + +import espressomd +import espressomd.lb +import unittest as ut +import unittest_decorators as utx +from lb_momentum_conservation import Momentum + + +@utx.skipIfMissingGPU() +@utx.skipIfMissingFeatures(['EXTERNAL_FORCES']) +class TestLBGPUMomentumConservationNSquare(Momentum, ut.TestCase): + + lb_class = espressomd.lb.LBFluidGPU + + def set_cellsystem(self): + self.system.cell_system.set_hybrid_decomposition( + n_square_types={0}, cutoff_regular=1) + + +@utx.skipIfMissingGPU() +@utx.skipIfMissingFeatures(['EXTERNAL_FORCES']) +class TestLBGPUMomentumConservationRegular(Momentum, ut.TestCase): + + lb_class = espressomd.lb.LBFluidGPU + + def set_cellsystem(self): + self.system.cell_system.set_hybrid_decomposition( + n_square_types={1}, cutoff_regular=1) + + +@utx.skipIfMissingFeatures(['EXTERNAL_FORCES']) +@ut.skipIf(Momentum.n_nodes > 1, + "LB with N-square only works on 1 MPI rank") +class TestLBCPUMomentumConservationNSquare(Momentum, ut.TestCase): + + lb_class = espressomd.lb.LBFluid + + def set_cellsystem(self): + self.system.cell_system.set_hybrid_decomposition( + n_square_types={0}, cutoff_regular=1) + + +@utx.skipIfMissingFeatures(['EXTERNAL_FORCES']) +@ut.skipIf(Momentum.n_nodes > 1, + "LB with N-square only works on 1 MPI rank") +class TestLBCPUMomentumConservationRegular(Momentum, ut.TestCase): + + lb_class = espressomd.lb.LBFluid + + def set_cellsystem(self): + self.system.cell_system.set_hybrid_decomposition( + n_square_types={1}, cutoff_regular=1) + + +if __name__ == "__main__": + ut.main() diff --git a/testsuite/python/lb_momentum_conservation_nsquare.py b/testsuite/python/lb_momentum_conservation_nsquare.py index 4e87f011e53..453e4e1d773 100644 --- a/testsuite/python/lb_momentum_conservation_nsquare.py +++ b/testsuite/python/lb_momentum_conservation_nsquare.py @@ -24,8 +24,6 @@ @utx.skipIfMissingGPU() @utx.skipIfMissingFeatures(['EXTERNAL_FORCES']) -@ut.skipIf(Momentum.n_nodes > 1, - "LB with N-square only works on 1 MPI rank") class TestLBGPUMomentumConservation(Momentum, ut.TestCase): lb_class = espressomd.lb.LBFluidGPU diff --git a/testsuite/python/lb_stats.py b/testsuite/python/lb_stats.py index 46f50094825..b77a45d8049 100644 --- a/testsuite/python/lb_stats.py +++ b/testsuite/python/lb_stats.py @@ -51,10 +51,12 @@ def tearDown(self): def test_mass_momentum_thermostat(self): self.n_col_part = 100 - partcls = self.system.part.add(pos=np.random.random( - (self.n_col_part, 3)) * self.system.box_l[0]) + # different types needed for lb_stats_hybrid test + particles = self.system.part.add( + type=self.n_col_part // 2 * [0, 1], pos=np.random.random( + (self.n_col_part, 3)) * self.system.box_l[0]) if espressomd.has_features("MASS"): - partcls.mass = 0.1 + np.random.random( + particles.mass = 0.1 + np.random.random( len(self.system.part)) self.system.thermostat.turn_off() diff --git a/testsuite/python/lb_stats_hybrid.py b/testsuite/python/lb_stats_hybrid.py new file mode 100644 index 00000000000..f2324dc24cd --- /dev/null +++ b/testsuite/python/lb_stats_hybrid.py @@ -0,0 +1,47 @@ +# Copyright (C) 2010-2022 The ESPResSo project +# +# This file is part of ESPResSo. +# +# ESPResSo is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# ESPResSo is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +import unittest as ut +import unittest_decorators as utx +from lb_stats import TestLB + +import espressomd +import espressomd.lb + + +@ut.skipIf(TestLB.n_nodes > 1, + "LB with N-square only works on 1 MPI rank") +class TestLBCPU(TestLB, ut.TestCase): + + def setUp(self): + self.system.cell_system.set_hybrid_decomposition( + n_square_types={0}, cutoff_regular=0) + self.lb_class = espressomd.lb.LBFluid + self.params.update({"mom_prec": 1E-9, "mass_prec_per_node": 5E-8}) + + +@utx.skipIfMissingGPU() +class TestLBGPU(TestLB, ut.TestCase): + + def setUp(self): + self.system.cell_system.set_hybrid_decomposition( + n_square_types={1}, cutoff_regular=0) + self.lb_class = espressomd.lb.LBFluidGPU + self.params.update({"mom_prec": 1E-3, "mass_prec_per_node": 1E-5}) + + +if __name__ == "__main__": + ut.main() diff --git a/testsuite/python/lb_stats_nsquare.py b/testsuite/python/lb_stats_nsquare.py index 9bf6e992952..f6fcc11db1c 100644 --- a/testsuite/python/lb_stats_nsquare.py +++ b/testsuite/python/lb_stats_nsquare.py @@ -33,8 +33,6 @@ def setUp(self): @utx.skipIfMissingGPU() -@ut.skipIf(TestLB.n_nodes > 1, - "LB with N-square only works on 1 MPI rank") class TestLBGPU(TestLB, ut.TestCase): def setUp(self): diff --git a/testsuite/python/regular_decomposition.py b/testsuite/python/regular_decomposition.py index 55b2ee62f29..fc667f36681 100644 --- a/testsuite/python/regular_decomposition.py +++ b/testsuite/python/regular_decomposition.py @@ -40,11 +40,11 @@ def check_resort(self): n_part = 2351 # Add the particles on node 0, so that they have to be resorted - partcls = self.system.part.add( + particles = self.system.part.add( pos=n_part * [(0, 0, 0)], type=n_part * [1]) # And now change their positions - partcls.pos = self.system.box_l * \ + particles.pos = self.system.box_l * \ np.random.random((n_part, 3)) # Add an interacting particle in a corner of the box From f0d0533f677e8a80129a3d42bb09ae0714eb6798 Mon Sep 17 00:00:00 2001 From: Patrick Kreissl Date: Mon, 25 Apr 2022 15:54:21 +0200 Subject: [PATCH 2/3] sphinx doc: Hybrid decomposition. --- doc/sphinx/system_setup.rst | 48 +++++++++++++++++++++++++++++++++++++ 1 file changed, 48 insertions(+) diff --git a/doc/sphinx/system_setup.rst b/doc/sphinx/system_setup.rst index 62dc61fe2a7..f17d66dcb4f 100644 --- a/doc/sphinx/system_setup.rst +++ b/doc/sphinx/system_setup.rst @@ -309,6 +309,54 @@ this node is twice as high. For 3 processors, the interactions are 0-0, Therefore it is highly recommended that you use N-squared only with an odd number of nodes, if with multiple processors at all. +.. _Hybrid: + +Hybrid decomposition +~~~~~~~~~~~~~~~~~~~~ + +If for a simulation setup the interaction range is much smaller than the +system size, use of a :ref:`Regular decomposition` leads to efficient +scaling behavior (order :math:`N` instead of order :math:`N^2`). +Consider a system with many small particles, e.g. a polymer solution. +There, already the addition of one single large particle increases the maximum +interaction range and thus the minimum cell size of the decomposition. +Due to this larger cell size, throughout the simulation box a large number +of non-interacting pairs of small particles is visited during the short +range calculation. This can considerably increase the computational cost of +the simulation. + +For such simulation setups, i.e. systems with a few large particles and much +more small particles, the hybrid decomposition can be used. This hybrid +decomposition is backed by two coupled particle decompositions which can +be used to efficiently deal with the differently sized particles. +Specifically that means putting the small particles into a +:ref:`Regular decomposition`. There, the minimum cell size is limited only +by the maximum interaction range of all particles within this decomposition. +The few large particles are put into a :ref:`N-squared` cellsystem. Particles +within this decomposition interact both, amongst each other and with all small +particles in the :ref:`Regular decomposition`. The hybrid decomposition can therefore +effectively recover the computational efficiency of the regular decomposition, +given that only a few large particles have been added. + +Invoking :py:meth:`~espressomd.cellsystem.CellSystem.set_hybrid_decomposition` +selects the hybrid decomposition. :: + + system = espressomd.System(box_l=[10, 10, 10]) + system.cell_system.set_hybrid_decomposition(n_square_types={1, 3}, cutoff_regular=1.2) + +Here, ``n_square_types`` is a python set containing the types of particles to +put into the :ref:`N-squared` cellsystem, i.e. the particle types of the +large particles. Particles with other types will by default be put into the +:ref:`Regular decomposition`. Note that for now it is also necessary to manually set +the maximum cutoff to consider for interactions within the +:ref:`Regular decomposition`, i.e. the maximum interaction range among all +small particle types. Set this via the ``cutoff_regular`` parameter. + +.. note:: + + The hybrid particle decomposition has been added to |es| only recently and + for now should be considered an experimental feature. If you notice some unexpected + behavior please let us know via github or the mailing list. .. _CUDA: From c11bb01ab5dba669b55ba5a6b07992e21909d44c Mon Sep 17 00:00:00 2001 From: Patrick Kreissl Date: Tue, 10 May 2022 16:34:05 +0200 Subject: [PATCH 3/3] Increase codecov. --- src/core/cells.cpp | 17 ----------------- src/core/cells.hpp | 12 ------------ testsuite/python/hybrid_decomposition.py | 12 ++++++++---- 3 files changed, 8 insertions(+), 33 deletions(-) diff --git a/src/core/cells.cpp b/src/core/cells.cpp index f3160de3082..f0535eb3eaf 100644 --- a/src/core/cells.cpp +++ b/src/core/cells.cpp @@ -196,17 +196,6 @@ std::vector mpi_resort_particles(int global_flag) { return n_parts; } -void set_regular_decomposition() { - cell_structure.set_regular_decomposition(comm_cart, interaction_range(), - box_geo, local_geo); - on_cell_structure_change(); -} - -void set_atom_decomposition() { - cell_structure.set_atom_decomposition(comm_cart, box_geo, local_geo); - on_cell_structure_change(); -} - void set_hybrid_decomposition(std::set n_square_types, double cutoff_regular) { cell_structure.set_hybrid_decomposition(comm_cart, cutoff_regular, box_geo, @@ -214,14 +203,8 @@ void set_hybrid_decomposition(std::set n_square_types, on_cell_structure_change(); } -REGISTER_CALLBACK(set_regular_decomposition) -REGISTER_CALLBACK(set_atom_decomposition) REGISTER_CALLBACK(set_hybrid_decomposition) -void mpi_set_regular_decomposition() { - mpi_call_all(set_regular_decomposition); -} -void mpi_set_atom_decomposition() { mpi_call_all(set_atom_decomposition); } void mpi_set_hybrid_decomposition(std::set n_square_types, double cutoff_regular) { mpi_call_all(set_hybrid_decomposition, n_square_types, cutoff_regular); diff --git a/src/core/cells.hpp b/src/core/cells.hpp index 132049b1834..ad1be861238 100644 --- a/src/core/cells.hpp +++ b/src/core/cells.hpp @@ -69,24 +69,12 @@ enum { /** Type of cell structure in use. */ extern CellStructure cell_structure; -/** Initialize cell sturcture DomainDecomposition */ -void set_regular_decomposition(); - -/** Initialize cell structure AtomDecomposition */ -void set_atom_decomposition(); - /** Initialize cell structure HybridDecomposition * @param n_square_types Set of particle types that will use n_square. */ void set_hybrid_decomposition(std::set n_square_types, double cutoff_regular); -/** Change cell structure to DomainDecomposition on all nodes. */ -void mpi_set_regular_decomposition(); - -/** Change cell structure to AtomDecomposition on all nodes. */ -void mpi_set_atom_decomposition(); - /** Change cell structure to HybridDecomposition on all nodes. * @param n_square_types Set of particle types that will use n_square. */ diff --git a/testsuite/python/hybrid_decomposition.py b/testsuite/python/hybrid_decomposition.py index 56448c2d17f..dd07d831674 100644 --- a/testsuite/python/hybrid_decomposition.py +++ b/testsuite/python/hybrid_decomposition.py @@ -177,6 +177,10 @@ def test_non_bonded_loop_trace(self): # p is a tuple with (id1, id2, pos1, pos2, vec21, mpi_node) # Note that system.distance_vec uses the opposite sign convention # as the minimum image distance in the core + self.assertTrue( + (p[0], p[1]) in py_distances or + (p[1], p[0]) in py_distances, + msg=f"Extra pair from core {p}") if (p[0], p[1]) in py_distances: np.testing.assert_allclose( np.copy(p[4]), -py_distances[p[0], p[1]]) @@ -185,14 +189,14 @@ def test_non_bonded_loop_trace(self): np.testing.assert_allclose( np.copy(p[4]), py_distances[p[1], p[0]]) del py_distances[p[1], p[0]] - else: - raise Exception("Extra pair from core", p) # test for pairs from regular child decomposition with more than one # cell for ids, dist in py_distances.items(): - if np.linalg.norm(dist) < cutoff_regular: - raise Exception("Pair not found by the core", ids) + self.assertGreaterEqual( + np.linalg.norm(dist), + cutoff_regular, + msg=f"Pair not found by the core {ids}") @utx.skipIfMissingFeatures(["LENNARD_JONES"]) def test_against_nsquare(self):