From 4f4040a26c0e29bfbbc1883d16b13a51885a09d2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Wed, 6 Apr 2022 17:45:52 +0200 Subject: [PATCH 1/7] core: Separate particle storage from particle data Move particle creation, deletion, storage, cache and tracking logic to a dedicated source file particle_node.cpp and leave particle property setters and getters in particle_data.cpp. --- src/core/CMakeLists.txt | 1 + src/core/particle_data.cpp | 487 +--------------------- src/core/particle_data.hpp | 126 +----- src/core/particle_node.cpp | 511 ++++++++++++++++++++++++ src/core/particle_node.hpp | 137 +++++++ src/python/espressomd/particle_data.pxd | 13 +- 6 files changed, 671 insertions(+), 604 deletions(-) create mode 100644 src/core/particle_node.cpp create mode 100644 src/core/particle_node.hpp diff --git a/src/core/CMakeLists.txt b/src/core/CMakeLists.txt index 1b632c873d6..3dc075f4d7b 100644 --- a/src/core/CMakeLists.txt +++ b/src/core/CMakeLists.txt @@ -22,6 +22,7 @@ set(EspressoCore_SRC npt.cpp partCfg_global.cpp particle_data.cpp + particle_node.cpp polymer.cpp pressure.cpp rattle.cpp diff --git a/src/core/particle_data.cpp b/src/core/particle_data.cpp index 17d8864a258..c515e91f9d8 100644 --- a/src/core/particle_data.cpp +++ b/src/core/particle_data.cpp @@ -18,55 +18,36 @@ * You should have received a copy of the GNU General Public License * along with this program. If not, see . */ -/** \file - * Particles and particle lists. - * - * The corresponding header file is particle_data.hpp. - */ + #include "particle_data.hpp" #include "Particle.hpp" #include "cells.hpp" #include "communication.hpp" +#include "config.hpp" #include "event.hpp" -#include "grid.hpp" #include "nonbonded_interactions/nonbonded_interaction_data.hpp" #include "partCfg_global.hpp" #include "rotation.hpp" -#include -#include -#include -#include +#include +#include #include -#include -#include -#include -#include -#include -#include +#include +#include #include #include #include -#include -#include -#include -#include -#include -#include +#include + +#include +#include #include -#include -#include #include constexpr auto some_tag = 42; -static bool type_list_enable; -static std::unordered_map> particle_type_map; - -void remove_id_from_map(int part_id, int type); -void add_id_to_type_map(int part_id, int type); namespace { /** @@ -420,324 +401,12 @@ void mpi_update_particle_property(int id, const T &value) { mpi_update_particle(id, value); } -/** - * @brief id -> rank - */ -std::unordered_map particle_node; - void delete_exclusion(Particle *part, int part2); void add_exclusion(Particle *part, int part2); void auto_exclusion(int distance); -static void mpi_who_has_local() { - static std::vector sendbuf; - - auto local_particles = cell_structure.local_particles(); - auto const n_part = static_cast(local_particles.size()); - boost::mpi::gather(comm_cart, n_part, 0); - - if (n_part == 0) - return; - - sendbuf.resize(n_part); - - std::transform(local_particles.begin(), local_particles.end(), - sendbuf.begin(), [](Particle const &p) { return p.id(); }); - - MPI_Send(sendbuf.data(), n_part, MPI_INT, 0, some_tag, comm_cart); -} - -REGISTER_CALLBACK(mpi_who_has_local) - -void mpi_who_has() { - mpi_call(mpi_who_has_local); - - auto local_particles = cell_structure.local_particles(); - - static std::vector n_parts; - boost::mpi::gather(comm_cart, static_cast(local_particles.size()), - n_parts, 0); - - static std::vector pdata; - - /* then fetch particle locations */ - for (int pnode = 0; pnode < n_nodes; pnode++) { - if (pnode == this_node) { - for (auto const &p : local_particles) - particle_node[p.id()] = this_node; - - } else if (n_parts[pnode] > 0) { - pdata.resize(n_parts[pnode]); - MPI_Recv(pdata.data(), n_parts[pnode], MPI_INT, pnode, some_tag, - comm_cart, MPI_STATUS_IGNORE); - for (int i = 0; i < n_parts[pnode]; i++) - particle_node[pdata[i]] = pnode; - } - } -} - -/** - * @brief Rebuild the particle index. - */ -void build_particle_node() { mpi_who_has(); } - -/** - * @brief Get the mpi rank which owns the particle with id. - */ -int get_particle_node(int id) { - if (id < 0) { - throw std::domain_error("Invalid particle id: " + std::to_string(id)); - } - - if (particle_node.empty()) - build_particle_node(); - - auto const needle = particle_node.find(id); - - // Check if particle has a node, if not, we assume it does not exist. - if (needle == particle_node.end()) { - throw std::runtime_error("Particle node for id " + std::to_string(id) + - " not found!"); - } - return needle->second; -} - -void clear_particle_node() { particle_node.clear(); } - -namespace { -/* Limit cache to 100 MiB */ -std::size_t const max_cache_size = (100ul * 1048576ul) / sizeof(Particle); -Utils::Cache particle_fetch_cache(max_cache_size); -} // namespace - -void invalidate_fetch_cache() { particle_fetch_cache.invalidate(); } -std::size_t fetch_cache_max_size() { return particle_fetch_cache.max_size(); } - -boost::optional get_particle_data_local(int id) { - auto p = cell_structure.get_local_particle(id); - - if (p and (not p->is_ghost())) { - return *p; - } - - return {}; -} - -REGISTER_CALLBACK_ONE_RANK(get_particle_data_local) - -const Particle &get_particle_data(int part) { - auto const pnode = get_particle_node(part); - - if (pnode == this_node) { - assert(cell_structure.get_local_particle(part)); - return *cell_structure.get_local_particle(part); - } - - /* Query the cache */ - auto const p_ptr = particle_fetch_cache.get(part); - if (p_ptr) { - return *p_ptr; - } - - /* Cache miss, fetch the particle, - * put it into the cache and return a pointer into the cache. */ - auto const cache_ptr = particle_fetch_cache.put( - part, Communication::mpiCallbacks().call(Communication::Result::one_rank, - get_particle_data_local, part)); - return *cache_ptr; -} - -static void mpi_get_particles_local() { - std::vector ids; - boost::mpi::scatter(comm_cart, ids, 0); - - std::vector parts(ids.size()); - std::transform(ids.begin(), ids.end(), parts.begin(), [](int id) { - assert(cell_structure.get_local_particle(id)); - return *cell_structure.get_local_particle(id); - }); - - Utils::Mpi::gatherv(comm_cart, parts.data(), static_cast(parts.size()), - 0); -} - -REGISTER_CALLBACK(mpi_get_particles_local) - -/** - * @brief Get multiple particles at once. - * - * *WARNING* Particles are returned in an arbitrary order. - * - * @param ids The ids of the particles that should be returned. - * - * @returns The particle list. - */ -std::vector mpi_get_particles(Utils::Span ids) { - mpi_call(mpi_get_particles_local); - /* Return value */ - std::vector parts(ids.size()); - - /* Group ids per node */ - static std::vector> node_ids(comm_cart.size()); - for (auto &per_node : node_ids) { - per_node.clear(); - } - - for (auto const &id : ids) { - auto const pnode = get_particle_node(id); - - if (pnode != this_node) - node_ids[pnode].push_back(id); - } - - /* Distributed ids to the nodes */ - { - static std::vector ignore; - boost::mpi::scatter(comm_cart, node_ids, ignore, 0); - } - - /* Copy local particles */ - std::transform(node_ids[this_node].cbegin(), node_ids[this_node].cend(), - parts.begin(), [](int id) { - assert(cell_structure.get_local_particle(id)); - return *cell_structure.get_local_particle(id); - }); - - static std::vector node_sizes(comm_cart.size()); - std::transform( - node_ids.cbegin(), node_ids.cend(), node_sizes.begin(), - [](std::vector const &ids) { return static_cast(ids.size()); }); - - Utils::Mpi::gatherv(comm_cart, parts.data(), static_cast(parts.size()), - parts.data(), node_sizes.data(), 0); - - return parts; -} - -void prefetch_particle_data(Utils::Span in_ids) { - /* Nothing to do on a single node. */ - // NOLINTNEXTLINE(clang-analyzer-core.NonNullParamChecker) - if (comm_cart.size() == 1) - return; - - static std::vector ids; - ids.clear(); - - boost::algorithm::copy_if(in_ids, std::back_inserter(ids), [](int id) { - return (get_particle_node(id) != this_node) && particle_fetch_cache.has(id); - }); - - /* Don't prefetch more particles than fit the cache. */ - if (ids.size() > particle_fetch_cache.max_size()) - ids.resize(particle_fetch_cache.max_size()); - - /* Fetch the particles... */ - for (auto &p : mpi_get_particles(ids)) { - auto id = p.id(); - particle_fetch_cache.put(id, std::move(p)); - } -} - -/** Move a particle to a new position. If it does not exist, it is created. - * The position must be on the local node! - * - * @param id the identity of the particle to move - * @param pos its new position - * @param _new if true, the particle is allocated, else has to exists already - * - * @return Pointer to the particle. - */ -Particle *local_place_particle(int id, const Utils::Vector3d &pos, int _new) { - auto pp = Utils::Vector3d{pos[0], pos[1], pos[2]}; - auto i = Utils::Vector3i{}; - fold_position(pp, i, box_geo); - - if (_new) { - Particle new_part; - new_part.id() = id; - new_part.pos() = pp; - new_part.image_box() = i; - - return cell_structure.add_local_particle(std::move(new_part)); - } - - auto pt = cell_structure.get_local_particle(id); - pt->pos() = pp; - pt->image_box() = i; - - return pt; -} - -boost::optional mpi_place_new_particle_local(int p_id, - Utils::Vector3d const &pos) { - auto p = local_place_particle(p_id, pos, 1); - on_particle_change(); - if (p) { - return comm_cart.rank(); - } - return {}; -} - -REGISTER_CALLBACK_ONE_RANK(mpi_place_new_particle_local) - -/** Create particle at a position on a node. - * Also calls \ref on_particle_change. - * \param p_id the particle to create. - * \param pos the particles position. - */ -int mpi_place_new_particle(int p_id, const Utils::Vector3d &pos) { - return mpi_call(Communication::Result::one_rank, mpi_place_new_particle_local, - p_id, pos); -} - -void mpi_place_particle_local(int pnode, int p_id) { - if (pnode == this_node) { - Utils::Vector3d pos; - comm_cart.recv(0, some_tag, pos); - local_place_particle(p_id, pos, 0); - } - - cell_structure.set_resort_particles(Cells::RESORT_GLOBAL); - on_particle_change(); -} - -REGISTER_CALLBACK(mpi_place_particle_local) - -/** Move particle to a position on a node. - * Also calls \ref on_particle_change. - * \param node the node to attach it to. - * \param p_id the particle to move. - * \param pos the particles position. - */ -void mpi_place_particle(int node, int p_id, const Utils::Vector3d &pos) { - mpi_call(mpi_place_particle_local, node, p_id); - - if (node == this_node) - local_place_particle(p_id, pos, 0); - else { - comm_cart.send(node, some_tag, pos); - } - - cell_structure.set_resort_particles(Cells::RESORT_GLOBAL); - on_particle_change(); -} - -int place_particle(int p_id, Utils::Vector3d const &pos) { - if (p_id < 0) { - throw std::domain_error("Invalid particle id: " + std::to_string(p_id)); - } - if (particle_exists(p_id)) { - mpi_place_particle(get_particle_node(p_id), p_id, pos); - - return ES_PART_OK; - } - particle_node[p_id] = mpi_place_new_particle(p_id, pos); - - return ES_PART_CREATED; -} - void set_particle_v(int part, Utils::Vector3d const &v) { mpi_update_particle(part, v); @@ -861,19 +530,7 @@ void set_particle_mu_E(int part, Utils::Vector3d const &mu_E) { void set_particle_type(int p_id, int type) { make_particle_type_exist(type); - - if (type_list_enable) { - // check if the particle exists already and the type is changed, then remove - // it from the list which contains it - auto const &cur_par = get_particle_data(p_id); - int prev_type = cur_par.type(); - if (prev_type != type) { - // particle existed before so delete it from the list - remove_id_from_map(p_id, prev_type); - } - add_id_to_type_map(p_id, type); - } - + on_particle_type_change(p_id, type); mpi_update_particle_property(p_id, type); } @@ -991,45 +648,6 @@ const std::vector &get_particle_bonds(int part) { return ret; } -static void mpi_remove_particle_local(int p_id) { - if (p_id != -1) { - cell_structure.remove_particle(p_id); - } else { - cell_structure.remove_all_particles(); - } - on_particle_change(); -} - -REGISTER_CALLBACK(mpi_remove_particle_local) - -/** Remove a particle. - * Also calls \ref on_particle_change. - * \param p_id the particle to remove, use -1 to remove all particles. - */ -void mpi_remove_particle(int p_id) { - mpi_call_all(mpi_remove_particle_local, p_id); -} - -void remove_all_particles() { - mpi_remove_particle(-1); - clear_particle_node(); -} - -int remove_particle(int p_id) { - auto const &cur_par = get_particle_data(p_id); - if (type_list_enable) { - // remove particle from its current type_list - int type = cur_par.type(); - remove_id_from_map(p_id, type); - } - - particle_node[p_id] = -1; - mpi_remove_particle(p_id); - particle_node.erase(p_id); - - return ES_OK; -} - /** Locally rescale all particles on current node. * @param dir direction to scale (0/1/2 = x/y/z, 3 = x+y+z isotropically) * @param scale factor by which to rescale (>1: stretch, <1: contract) @@ -1187,86 +805,3 @@ void auto_exclusions(int distance) { } } #endif // EXCLUSIONS - -void init_type_map(int type) { - type_list_enable = true; - if (type < 0) - throw std::runtime_error("Types may not be negative"); - - auto &map_for_type = particle_type_map[type]; - map_for_type.clear(); - for (auto const &p : partCfg()) { - if (p.type() == type) - map_for_type.insert(p.id()); - } -} - -void remove_id_from_map(int part_id, int type) { - auto it = particle_type_map.find(type); - if (it != particle_type_map.end()) - it->second.erase(part_id); -} - -int get_random_p_id(int type, int random_index_in_type_map) { - auto it = particle_type_map.find(type); - if (it == particle_type_map.end()) { - throw std::runtime_error("The provided particle type " + - std::to_string(type) + - " is currently not tracked by the system."); - } - - if (random_index_in_type_map + 1 > it->second.size()) - throw std::runtime_error("The provided index exceeds the number of " - "particle types listed in the particle_type_map"); - return *std::next(it->second.begin(), random_index_in_type_map); -} - -void add_id_to_type_map(int part_id, int type) { - auto it = particle_type_map.find(type); - if (it != particle_type_map.end()) - it->second.insert(part_id); -} - -int number_of_particles_with_type(int type) { - auto it = particle_type_map.find(type); - if (it == particle_type_map.end()) { - throw std::runtime_error("The provided particle type " + - std::to_string(type) + - " is currently not tracked by the system."); - } - - return static_cast(it->second.size()); -} - -bool particle_exists(int part_id) { - if (particle_node.empty()) - build_particle_node(); - return particle_node.count(part_id); -} - -std::vector get_particle_ids() { - if (particle_node.empty()) - build_particle_node(); - - auto ids = Utils::keys(particle_node); - boost::sort(ids); - - return ids; -} - -int get_maximal_particle_id() { - if (particle_node.empty()) - build_particle_node(); - - return boost::accumulate(particle_node, -1, - [](int max, const std::pair &kv) { - return std::max(max, kv.first); - }); -} - -int get_n_part() { - if (particle_node.empty()) - build_particle_node(); - - return static_cast(particle_node.size()); -} diff --git a/src/core/particle_data.hpp b/src/core/particle_data.hpp index 7b39c92f154..ed7e6c86d74 100644 --- a/src/core/particle_data.hpp +++ b/src/core/particle_data.hpp @@ -21,95 +21,27 @@ #ifndef CORE_PARTICLE_DATA_HPP #define CORE_PARTICLE_DATA_HPP /** \file - * Particles and particle lists. + * Particles property access. * - * This file contains everything related to particle storage. If you want to + * This file contains everything related to particle properties. If you want to * add a new property to the particles, it is probably a good idea to modify * Particle to give scripts access to that property. You always have to modify * two positions: first the print section, where you should add your new * data at the end, and second the read section where you have to find a nice - * and short name for your property to appear in the Python code. Then you - * just parse your part out of argc and argv. - * - * Implementation in particle_data.cpp. + * and short name for your property to appear in the Python code. */ #include "config.hpp" #include "Particle.hpp" -#include "ParticleList.hpp" +#include "particle_node.hpp" #include #include #include -#include #include -/************************************************ - * defines - ************************************************/ - -enum { - /// ok code for \ref place_particle - ES_PART_OK = 0, - /// error code for \ref place_particle - ES_PART_ERROR = -1, - /// ok code for \ref place_particle, particle is new - ES_PART_CREATED = 1 -}; - -/************************************************ - * Functions - ************************************************/ - -/** Invalidate \ref particle_node. This has to be done - * at the beginning of the integration. - */ -void clear_particle_node(); - -/** - * @brief Get particle data. - * - * @param part the identity of the particle to fetch - * @return Pointer to copy of particle if it exists, - * nullptr otherwise; - */ -const Particle &get_particle_data(int part); - -/** - * @brief Fetch a range of particle into the fetch cache. - * - * - * If the range is larger than the cache size, only - * the particle that fit into the cache are fetched. - * - * The particles have to exist, an exception it throw - * if one of the the particles can not be found. - * - * @param ids Ids of the particles that should be fetched. - */ -void prefetch_particle_data(Utils::Span ids); - -/** @brief Invalidate the fetch cache for get_particle_data. */ -void invalidate_fetch_cache(); - -/** @brief Return the maximal number of particles that are - * kept in the fetch cache. - */ -std::size_t fetch_cache_max_size(); - -/** Call only on the head node. - * Move a particle to a new position. - * If it does not exist, it is created. - * @param part the identity of the particle to move - * @param p its new position - * @retval ES_PART_OK if particle existed - * @retval ES_PART_CREATED if created - * @retval ES_PART_ERROR if id is illegal - */ -int place_particle(int part, Utils::Vector3d const &p); - /** Call only on the head node: set particle velocity. * @param part the particle. * @param v its new velocity. @@ -339,17 +271,6 @@ const std::vector &get_particle_bonds(int part); int change_exclusion(int part, int part2, int _delete); #endif -/** Remove particle with a given identity. Also removes all bonds to the - * particle. - * @param part identity of the particle to remove - * @retval ES_OK on success - * @retval ES_ERROR if particle does not exist - */ -int remove_particle(int part); - -/** Remove all particles. */ -void remove_all_particles(); - /** Rescale all particle positions in direction @p dir by a factor @p scale. */ void mpi_rescale_particles(int dir, double scale); @@ -362,12 +283,6 @@ void mpi_rescale_particles(int dir, double scale); */ void auto_exclusions(int distance); -void init_type_map(int type); - -/** Find a particle of given type and return its id */ -int get_random_p_id(int type, int random_index_in_type_map); -int number_of_particles_with_type(int type); - // The following functions are used by the python interface to obtain // properties of a particle, which are only compiled in in some configurations // This is needed, because cython does not support conditional compilation @@ -421,37 +336,4 @@ inline Utils::Vector3i get_particle_rotation(Particle const *p) { } #endif // ROTATION -/** - * @brief Check if particle exists. - * - * @param part Id of the particle - * @return True iff the particle exists. - */ -bool particle_exists(int part); - -/** - * @brief Get the mpi rank which owns the particle with id. - * - * @param id Id of the particle - * @return The MPI rank the particle is on. - */ -int get_particle_node(int id); - -/** - * @brief Get all particle ids. - * - * @return Sorted ids of all existing particles. - */ -std::vector get_particle_ids(); - -/** - * @brief Get maximal particle id. - */ -int get_maximal_particle_id(); - -/** - * @brief Get number of particles. - */ -int get_n_part(); - #endif diff --git a/src/core/particle_node.cpp b/src/core/particle_node.cpp new file mode 100644 index 00000000000..dfdf2b14539 --- /dev/null +++ b/src/core/particle_node.cpp @@ -0,0 +1,511 @@ +/* + * Copyright (C) 2010-2019 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 "particle_node.hpp" + +#include "Particle.hpp" +#include "cells.hpp" +#include "communication.hpp" +#include "event.hpp" +#include "grid.hpp" +#include "partCfg_global.hpp" +#include "particle_data.hpp" + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +enum { + /// ok code for \ref place_particle + ES_PART_OK = 0, + /// error code for \ref place_particle + ES_PART_ERROR = -1, + /// ok code for \ref place_particle, particle is new + ES_PART_CREATED = 1 +}; + +constexpr auto some_tag = 42; +static bool type_list_enable; +static std::unordered_map> particle_type_map; +std::unordered_map particle_node; + +void init_type_map(int type) { + type_list_enable = true; + if (type < 0) + throw std::runtime_error("Types may not be negative"); + + auto &map_for_type = particle_type_map[type]; + map_for_type.clear(); + for (auto const &p : partCfg()) { + if (p.type() == type) + map_for_type.insert(p.id()); + } +} + +void remove_id_from_map(int part_id, int type) { + auto it = particle_type_map.find(type); + if (it != particle_type_map.end()) + it->second.erase(part_id); +} + +void add_id_to_type_map(int part_id, int type) { + auto it = particle_type_map.find(type); + if (it != particle_type_map.end()) + it->second.insert(part_id); +} + +void on_particle_type_change(int p_id, int type) { + if (type_list_enable) { + // check if the particle exists already and the type is changed, then remove + // it from the list which contains it + auto const &cur_par = get_particle_data(p_id); + int prev_type = cur_par.type(); + if (prev_type != type) { + // particle existed before so delete it from the list + remove_id_from_map(p_id, prev_type); + } + add_id_to_type_map(p_id, type); + } +} + +namespace { +/* Limit cache to 100 MiB */ +std::size_t const max_cache_size = (100ul * 1048576ul) / sizeof(Particle); +Utils::Cache particle_fetch_cache(max_cache_size); +} // namespace + +void invalidate_fetch_cache() { particle_fetch_cache.invalidate(); } +std::size_t fetch_cache_max_size() { return particle_fetch_cache.max_size(); } + +boost::optional get_particle_data_local(int id) { + auto p = cell_structure.get_local_particle(id); + + if (p and (not p->is_ghost())) { + return *p; + } + + return {}; +} + +REGISTER_CALLBACK_ONE_RANK(get_particle_data_local) + +const Particle &get_particle_data(int part) { + auto const pnode = get_particle_node(part); + + if (pnode == this_node) { + assert(cell_structure.get_local_particle(part)); + return *cell_structure.get_local_particle(part); + } + + /* Query the cache */ + auto const p_ptr = particle_fetch_cache.get(part); + if (p_ptr) { + return *p_ptr; + } + + /* Cache miss, fetch the particle, + * put it into the cache and return a pointer into the cache. */ + auto const cache_ptr = particle_fetch_cache.put( + part, Communication::mpiCallbacks().call(Communication::Result::one_rank, + get_particle_data_local, part)); + return *cache_ptr; +} + +static void mpi_get_particles_local() { + std::vector ids; + boost::mpi::scatter(comm_cart, ids, 0); + + std::vector parts(ids.size()); + std::transform(ids.begin(), ids.end(), parts.begin(), [](int id) { + assert(cell_structure.get_local_particle(id)); + return *cell_structure.get_local_particle(id); + }); + + Utils::Mpi::gatherv(comm_cart, parts.data(), static_cast(parts.size()), + 0); +} + +REGISTER_CALLBACK(mpi_get_particles_local) + +/** + * @brief Get multiple particles at once. + * + * *WARNING* Particles are returned in an arbitrary order. + * + * @param ids The ids of the particles that should be returned. + * + * @returns The particle list. + */ +std::vector mpi_get_particles(Utils::Span ids) { + mpi_call(mpi_get_particles_local); + /* Return value */ + std::vector parts(ids.size()); + + /* Group ids per node */ + static std::vector> node_ids(comm_cart.size()); + for (auto &per_node : node_ids) { + per_node.clear(); + } + + for (auto const &id : ids) { + auto const pnode = get_particle_node(id); + + if (pnode != this_node) + node_ids[pnode].push_back(id); + } + + /* Distributed ids to the nodes */ + { + static std::vector ignore; + boost::mpi::scatter(comm_cart, node_ids, ignore, 0); + } + + /* Copy local particles */ + std::transform(node_ids[this_node].cbegin(), node_ids[this_node].cend(), + parts.begin(), [](int id) { + assert(cell_structure.get_local_particle(id)); + return *cell_structure.get_local_particle(id); + }); + + static std::vector node_sizes(comm_cart.size()); + std::transform( + node_ids.cbegin(), node_ids.cend(), node_sizes.begin(), + [](std::vector const &ids) { return static_cast(ids.size()); }); + + Utils::Mpi::gatherv(comm_cart, parts.data(), static_cast(parts.size()), + parts.data(), node_sizes.data(), 0); + + return parts; +} + +void prefetch_particle_data(Utils::Span in_ids) { + /* Nothing to do on a single node. */ + // NOLINTNEXTLINE(clang-analyzer-core.NonNullParamChecker) + if (comm_cart.size() == 1) + return; + + static std::vector ids; + ids.clear(); + + boost::algorithm::copy_if(in_ids, std::back_inserter(ids), [](int id) { + return (get_particle_node(id) != this_node) && particle_fetch_cache.has(id); + }); + + /* Don't prefetch more particles than fit the cache. */ + if (ids.size() > particle_fetch_cache.max_size()) + ids.resize(particle_fetch_cache.max_size()); + + /* Fetch the particles... */ + for (auto &p : mpi_get_particles(ids)) { + auto id = p.id(); + particle_fetch_cache.put(id, std::move(p)); + } +} + +static void mpi_who_has_local() { + static std::vector sendbuf; + + auto local_particles = cell_structure.local_particles(); + auto const n_part = static_cast(local_particles.size()); + boost::mpi::gather(comm_cart, n_part, 0); + + if (n_part == 0) + return; + + sendbuf.resize(n_part); + + std::transform(local_particles.begin(), local_particles.end(), + sendbuf.begin(), [](Particle const &p) { return p.id(); }); + + MPI_Send(sendbuf.data(), n_part, MPI_INT, 0, some_tag, comm_cart); +} + +REGISTER_CALLBACK(mpi_who_has_local) + +void mpi_who_has() { + mpi_call(mpi_who_has_local); + + auto local_particles = cell_structure.local_particles(); + + static std::vector n_parts; + boost::mpi::gather(comm_cart, static_cast(local_particles.size()), + n_parts, 0); + + static std::vector pdata; + + /* then fetch particle locations */ + for (int pnode = 0; pnode < n_nodes; pnode++) { + if (pnode == this_node) { + for (auto const &p : local_particles) + particle_node[p.id()] = this_node; + + } else if (n_parts[pnode] > 0) { + pdata.resize(n_parts[pnode]); + MPI_Recv(pdata.data(), n_parts[pnode], MPI_INT, pnode, some_tag, + comm_cart, MPI_STATUS_IGNORE); + for (int i = 0; i < n_parts[pnode]; i++) + particle_node[pdata[i]] = pnode; + } + } +} + +/** + * @brief Rebuild the particle index. + */ +void build_particle_node() { mpi_who_has(); } + +/** + * @brief Get the mpi rank which owns the particle with id. + */ +int get_particle_node(int id) { + if (id < 0) { + throw std::domain_error("Invalid particle id: " + std::to_string(id)); + } + + if (particle_node.empty()) + build_particle_node(); + + auto const needle = particle_node.find(id); + + // Check if particle has a node, if not, we assume it does not exist. + if (needle == particle_node.end()) { + throw std::runtime_error("Particle node for id " + std::to_string(id) + + " not found!"); + } + return needle->second; +} + +void clear_particle_node() { particle_node.clear(); } + +/** Move a particle to a new position. If it does not exist, it is created. + * The position must be on the local node! + * + * @param id the identity of the particle to move + * @param pos its new position + * @param _new if true, the particle is allocated, else has to exists already + * + * @return Pointer to the particle. + */ +Particle *local_place_particle(int id, const Utils::Vector3d &pos, int _new) { + auto pp = Utils::Vector3d{pos[0], pos[1], pos[2]}; + auto i = Utils::Vector3i{}; + fold_position(pp, i, box_geo); + + if (_new) { + Particle new_part; + new_part.id() = id; + new_part.pos() = pp; + new_part.image_box() = i; + + return cell_structure.add_local_particle(std::move(new_part)); + } + + auto pt = cell_structure.get_local_particle(id); + pt->pos() = pp; + pt->image_box() = i; + + return pt; +} + +boost::optional mpi_place_new_particle_local(int p_id, + Utils::Vector3d const &pos) { + auto p = local_place_particle(p_id, pos, 1); + on_particle_change(); + if (p) { + return comm_cart.rank(); + } + return {}; +} + +REGISTER_CALLBACK_ONE_RANK(mpi_place_new_particle_local) + +/** Create particle at a position on a node. + * Also calls \ref on_particle_change. + * \param p_id the particle to create. + * \param pos the particles position. + */ +int mpi_place_new_particle(int p_id, const Utils::Vector3d &pos) { + return mpi_call(Communication::Result::one_rank, mpi_place_new_particle_local, + p_id, pos); +} + +void mpi_place_particle_local(int pnode, int p_id) { + if (pnode == this_node) { + Utils::Vector3d pos; + comm_cart.recv(0, some_tag, pos); + local_place_particle(p_id, pos, 0); + } + + cell_structure.set_resort_particles(Cells::RESORT_GLOBAL); + on_particle_change(); +} + +REGISTER_CALLBACK(mpi_place_particle_local) + +/** Move particle to a position on a node. + * Also calls \ref on_particle_change. + * \param node the node to attach it to. + * \param p_id the particle to move. + * \param pos the particles position. + */ +void mpi_place_particle(int node, int p_id, const Utils::Vector3d &pos) { + mpi_call(mpi_place_particle_local, node, p_id); + + if (node == this_node) + local_place_particle(p_id, pos, 0); + else { + comm_cart.send(node, some_tag, pos); + } + + cell_structure.set_resort_particles(Cells::RESORT_GLOBAL); + on_particle_change(); +} + +int place_particle(int p_id, Utils::Vector3d const &pos) { + if (p_id < 0) { + throw std::domain_error("Invalid particle id: " + std::to_string(p_id)); + } + if (particle_exists(p_id)) { + mpi_place_particle(get_particle_node(p_id), p_id, pos); + + return ES_PART_OK; + } + particle_node[p_id] = mpi_place_new_particle(p_id, pos); + + return ES_PART_CREATED; +} + +static void mpi_remove_particle_local(int p_id) { + if (p_id != -1) { + cell_structure.remove_particle(p_id); + } else { + cell_structure.remove_all_particles(); + } + on_particle_change(); +} + +REGISTER_CALLBACK(mpi_remove_particle_local) + +/** Remove a particle. + * Also calls \ref on_particle_change. + * \param p_id the particle to remove, use -1 to remove all particles. + */ +void mpi_remove_particle(int p_id) { + mpi_call_all(mpi_remove_particle_local, p_id); +} + +void remove_all_particles() { + mpi_remove_particle(-1); + clear_particle_node(); +} + +int remove_particle(int p_id) { + auto const &cur_par = get_particle_data(p_id); + if (type_list_enable) { + // remove particle from its current type_list + int type = cur_par.type(); + remove_id_from_map(p_id, type); + } + + particle_node[p_id] = -1; + mpi_remove_particle(p_id); + particle_node.erase(p_id); + + return ES_OK; +} + +int get_random_p_id(int type, int random_index_in_type_map) { + auto it = particle_type_map.find(type); + if (it == particle_type_map.end()) { + throw std::runtime_error("The provided particle type " + + std::to_string(type) + + " is currently not tracked by the system."); + } + + if (random_index_in_type_map + 1 > it->second.size()) + throw std::runtime_error("The provided index exceeds the number of " + "particle types listed in the particle_type_map"); + return *std::next(it->second.begin(), random_index_in_type_map); +} + +int number_of_particles_with_type(int type) { + auto it = particle_type_map.find(type); + if (it == particle_type_map.end()) { + throw std::runtime_error("The provided particle type " + + std::to_string(type) + + " is currently not tracked by the system."); + } + + return static_cast(it->second.size()); +} + +bool particle_exists(int part_id) { + if (particle_node.empty()) + build_particle_node(); + return particle_node.count(part_id); +} + +std::vector get_particle_ids() { + if (particle_node.empty()) + build_particle_node(); + + auto ids = Utils::keys(particle_node); + boost::sort(ids); + + return ids; +} + +int get_maximal_particle_id() { + if (particle_node.empty()) + build_particle_node(); + + return boost::accumulate(particle_node, -1, + [](int max, const std::pair &kv) { + return std::max(max, kv.first); + }); +} + +int get_n_part() { + if (particle_node.empty()) + build_particle_node(); + + return static_cast(particle_node.size()); +} diff --git a/src/core/particle_node.hpp b/src/core/particle_node.hpp new file mode 100644 index 00000000000..733a74cdfbe --- /dev/null +++ b/src/core/particle_node.hpp @@ -0,0 +1,137 @@ +/* + * Copyright (C) 2010-2019 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 CORE_PARTICLE_NODE_HPP +#define CORE_PARTICLE_NODE_HPP +/** \file + * Particles creation and deletion. + * + * This file contains everything related to particle storage and tracking. + * + * Implementation in particle_node.cpp. + */ + +#include "Particle.hpp" + +#include +#include + +#include +#include + +/** + * @brief Get particle data. + * + * @param part the identity of the particle to fetch + * @return Pointer to copy of particle if it exists, + * nullptr otherwise; + */ +const Particle &get_particle_data(int part); + +/** + * @brief Fetch a range of particle into the fetch cache. + * + * + * If the range is larger than the cache size, only + * the particle that fit into the cache are fetched. + * + * The particles have to exist, an exception it throw + * if one of the the particles can not be found. + * + * @param ids Ids of the particles that should be fetched. + */ +void prefetch_particle_data(Utils::Span ids); + +/** @brief Invalidate the fetch cache for get_particle_data. */ +void invalidate_fetch_cache(); + +/** @brief Return the maximal number of particles that are + * kept in the fetch cache. + */ +std::size_t fetch_cache_max_size(); + +/** Invalidate \ref particle_node. This has to be done + * at the beginning of the integration. + */ +void clear_particle_node(); + +/** Call only on the head node. + * Move a particle to a new position. + * If it does not exist, it is created. + * @param part the identity of the particle to move + * @param p its new position + * @retval ES_PART_OK if particle existed + * @retval ES_PART_CREATED if created + * @retval ES_PART_ERROR if id is illegal + */ +int place_particle(int part, Utils::Vector3d const &p); + +/** Remove particle with a given identity. Also removes all bonds to the + * particle. + * @param part identity of the particle to remove + * @retval ES_OK on success + * @retval ES_ERROR if particle does not exist + */ +int remove_particle(int part); + +/** Remove all particles. */ +void remove_all_particles(); + +void init_type_map(int type); +void on_particle_type_change(int p_id, int type); + +/** Find a particle of given type and return its id */ +int get_random_p_id(int type, int random_index_in_type_map); +int number_of_particles_with_type(int type); + +/** + * @brief Check if particle exists. + * + * @param part Id of the particle + * @return True iff the particle exists. + */ +bool particle_exists(int part); + +/** + * @brief Get the mpi rank which owns the particle with id. + * + * @param id Id of the particle + * @return The MPI rank the particle is on. + */ +int get_particle_node(int id); + +/** + * @brief Get all particle ids. + * + * @return Sorted ids of all existing particles. + */ +std::vector get_particle_ids(); + +/** + * @brief Get maximal particle id. + */ +int get_maximal_particle_id(); + +/** + * @brief Get number of particles. + */ +int get_n_part(); + +#endif diff --git a/src/python/espressomd/particle_data.pxd b/src/python/espressomd/particle_data.pxd index f1683a478b8..a3baa3910c1 100644 --- a/src/python/espressomd/particle_data.pxd +++ b/src/python/espressomd/particle_data.pxd @@ -26,8 +26,7 @@ include "myconfig.pxi" from .utils cimport Span -# Import particle data structures and setter functions from particle_data.hpp -cdef extern from "particle_data.hpp": +cdef extern from "Particle.hpp": cppclass BondView: int bond_id() Span[const int] partner_ids() @@ -71,11 +70,10 @@ cdef extern from "particle_data.hpp": bool has_exclusion(int pid) except + particle_parameters_swimming swimming() +cdef extern from "particle_data.hpp": # Setter/getter/modifier functions functions void prefetch_particle_data(vector[int] ids) - int place_particle(int part, const Vector3d & p) except + - void set_particle_v(int part, const Vector3d & v) void set_particle_f(int part, const Vector3d & f) @@ -156,12 +154,15 @@ cdef extern from "particle_data.hpp": IF ENGINE: void set_particle_swimming(int part, particle_parameters_swimming swim) + void remove_all_bonds_to(int part) + +cdef extern from "particle_node.hpp": + int place_particle(int part, const Vector3d & p) except + + int remove_particle(int part) except + void remove_all_particles() except + - void remove_all_bonds_to(int part) - bool particle_exists(int part) int get_particle_node(int pid) except + From 7af7b9c3e1388c330b42a267fc5baa48f1a50647 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Wed, 6 Apr 2022 18:06:09 +0200 Subject: [PATCH 2/7] core: Narrow interface of particle management Narrow includes of particle_data.hpp and particle_node.hpp. Make particle creation/deletion functions return void, since internal details of the particle storage are irrelevant. Make utility functions private (internal linkage). Replace MPI primitives by their equivalent boost::mpi functions. --- src/core/PartCfg.cpp | 2 +- src/core/cells.cpp | 2 +- src/core/cluster_analysis/Cluster.cpp | 2 +- .../cluster_analysis/ClusterStructure.cpp | 1 + .../mdlc_correction.cpp | 1 + src/core/event.cpp | 2 +- src/core/observables/fetch_particles.hpp | 2 +- src/core/pair_criteria/pair_criteria.hpp | 2 +- src/core/particle_data.cpp | 17 +++-- src/core/particle_data.hpp | 6 +- src/core/particle_node.cpp | 62 +++++++------------ src/core/particle_node.hpp | 15 ++--- .../reaction_methods/ReactionAlgorithm.cpp | 1 + .../tests/ConstantpHEnsemble_test.cpp | 1 - .../tests/ReactionAlgorithm_test.cpp | 1 + .../tests/ReactionEnsemble_test.cpp | 1 + src/core/statistics_chain.cpp | 2 +- .../EspressoSystemInterface_test.cpp | 1 + .../EspressoSystemStandAlone_test.cpp | 1 + src/core/unit_tests/ParticleFactory.hpp | 1 + src/core/unit_tests/Verlet_list_test.cpp | 2 +- src/core/virtual_sites.cpp | 1 + src/python/espressomd/particle_data.pxd | 4 +- src/python/espressomd/particle_data.pyx | 10 +-- src/python/espressomd/system.pxd | 2 +- testsuite/python/particle.py | 4 +- 26 files changed, 66 insertions(+), 80 deletions(-) diff --git a/src/core/PartCfg.cpp b/src/core/PartCfg.cpp index a2437195526..62d7d3494f7 100644 --- a/src/core/PartCfg.cpp +++ b/src/core/PartCfg.cpp @@ -20,7 +20,7 @@ #include "PartCfg.hpp" #include "grid.hpp" -#include "particle_data.hpp" +#include "particle_node.hpp" #include diff --git a/src/core/cells.cpp b/src/core/cells.cpp index 0c7b95ff83b..cc76b46f406 100644 --- a/src/core/cells.cpp +++ b/src/core/cells.cpp @@ -32,7 +32,7 @@ #include "event.hpp" #include "grid.hpp" #include "integrate.hpp" -#include "particle_data.hpp" +#include "particle_node.hpp" #include "ParticleDecomposition.hpp" #include "RegularDecomposition.hpp" diff --git a/src/core/cluster_analysis/Cluster.cpp b/src/core/cluster_analysis/Cluster.cpp index 0841aaccefd..dd579daa64d 100644 --- a/src/core/cluster_analysis/Cluster.cpp +++ b/src/core/cluster_analysis/Cluster.cpp @@ -29,7 +29,7 @@ #include "Particle.hpp" #include "errorhandling.hpp" #include "grid.hpp" -#include "particle_data.hpp" +#include "particle_node.hpp" #include diff --git a/src/core/cluster_analysis/ClusterStructure.cpp b/src/core/cluster_analysis/ClusterStructure.cpp index 2a092519dfd..005e36db3bb 100644 --- a/src/core/cluster_analysis/ClusterStructure.cpp +++ b/src/core/cluster_analysis/ClusterStructure.cpp @@ -22,6 +22,7 @@ #include "PartCfg.hpp" #include "errorhandling.hpp" #include "partCfg_global.hpp" +#include "particle_node.hpp" #include diff --git a/src/core/electrostatics_magnetostatics/mdlc_correction.cpp b/src/core/electrostatics_magnetostatics/mdlc_correction.cpp index f280ef1aa4f..5885741a92f 100644 --- a/src/core/electrostatics_magnetostatics/mdlc_correction.cpp +++ b/src/core/electrostatics_magnetostatics/mdlc_correction.cpp @@ -33,6 +33,7 @@ #include "errorhandling.hpp" #include "grid.hpp" #include "particle_data.hpp" +#include "particle_node.hpp" #include #include diff --git a/src/core/event.cpp b/src/core/event.cpp index 5324e5af7e8..9acc8934515 100644 --- a/src/core/event.cpp +++ b/src/core/event.cpp @@ -47,7 +47,7 @@ #include "nonbonded_interactions/nonbonded_interaction_data.hpp" #include "npt.hpp" #include "partCfg_global.hpp" -#include "particle_data.hpp" +#include "particle_node.hpp" #include "thermostat.hpp" #include "virtual_sites.hpp" diff --git a/src/core/observables/fetch_particles.hpp b/src/core/observables/fetch_particles.hpp index cf5200166dd..cd15435c855 100644 --- a/src/core/observables/fetch_particles.hpp +++ b/src/core/observables/fetch_particles.hpp @@ -21,7 +21,7 @@ #define FETCH_PARTICLES_HPP #include "grid.hpp" -#include "particle_data.hpp" +#include "particle_node.hpp" #include diff --git a/src/core/pair_criteria/pair_criteria.hpp b/src/core/pair_criteria/pair_criteria.hpp index 7f249b95a7a..d8b8a22db42 100644 --- a/src/core/pair_criteria/pair_criteria.hpp +++ b/src/core/pair_criteria/pair_criteria.hpp @@ -22,7 +22,7 @@ #include "Particle.hpp" #include "energy_inline.hpp" #include "grid.hpp" -#include "particle_data.hpp" +#include "particle_node.hpp" namespace PairCriteria { /** @brief Criterion which provides a true/false for a pair of particles */ diff --git a/src/core/particle_data.cpp b/src/core/particle_data.cpp index c515e91f9d8..5b38e348779 100644 --- a/src/core/particle_data.cpp +++ b/src/core/particle_data.cpp @@ -28,6 +28,7 @@ #include "event.hpp" #include "nonbonded_interactions/nonbonded_interaction_data.hpp" #include "partCfg_global.hpp" +#include "particle_node.hpp" #include "rotation.hpp" #include @@ -40,8 +41,6 @@ #include #include -#include - #include #include #include @@ -342,7 +341,7 @@ void local_remove_pair_bonds_to(Particle &p, int other_pid) { RemovePairBondsTo{other_pid}(p); } -void mpi_send_update_message_local(int node, int id) { +static void mpi_send_update_message_local(int node, int id) { if (node == comm_cart.rank()) { UpdateMessage msg{}; comm_cart.recv(0, some_tag, msg); @@ -372,7 +371,7 @@ REGISTER_CALLBACK(mpi_send_update_message_local) * @param id Id of the particle to update * @param msg The message */ -void mpi_send_update_message(int id, const UpdateMessage &msg) { +static void mpi_send_update_message(int id, const UpdateMessage &msg) { auto const pnode = get_particle_node(id); mpi_call(mpi_send_update_message_local, pnode, id); @@ -664,7 +663,7 @@ void local_rescale_particles(int dir, double scale) { static void mpi_rescale_particles_local(int dir) { double scale = 0.0; - MPI_Recv(&scale, 1, MPI_DOUBLE, 0, some_tag, comm_cart, MPI_STATUS_IGNORE); + comm_cart.recv(0, some_tag, scale); local_rescale_particles(dir, scale); on_particle_change(); } @@ -677,7 +676,7 @@ void mpi_rescale_particles(int dir, double scale) { if (pnode == this_node) { local_rescale_particles(dir, scale); } else { - MPI_Send(&scale, 1, MPI_DOUBLE, pnode, some_tag, comm_cart); + comm_cart.send(pnode, some_tag, scale); } } on_particle_change(); @@ -689,7 +688,7 @@ void mpi_rescale_particles(int dir, double scale) { * @param part2 the identity of the second exclusion partner * @param _delete if true, delete the exclusion instead of add */ -void local_change_exclusion(int part1, int part2, int _delete) { +static void local_change_exclusion(int part1, int part2, int _delete) { /* part1, if here */ auto part = cell_structure.get_local_particle(part1); if (part) { @@ -724,7 +723,7 @@ void add_partner(std::vector &il, int i, int j, int distance) { } } // namespace -void mpi_send_exclusion_local(int part1, int part2, int _delete) { +static void mpi_send_exclusion_local(int part1, int part2, int _delete) { local_change_exclusion(part1, part2, _delete); on_particle_change(); } @@ -737,7 +736,7 @@ REGISTER_CALLBACK(mpi_send_exclusion_local) * \param part2 identity of second particle of the exclusion. * \param _delete if true, do not add the exclusion, rather delete it if found */ -void mpi_send_exclusion(int part1, int part2, int _delete) { +static void mpi_send_exclusion(int part1, int part2, int _delete) { mpi_call_all(mpi_send_exclusion_local, part1, part2, _delete); } diff --git a/src/core/particle_data.hpp b/src/core/particle_data.hpp index ed7e6c86d74..c65eafbb8a2 100644 --- a/src/core/particle_data.hpp +++ b/src/core/particle_data.hpp @@ -34,7 +34,6 @@ #include "config.hpp" #include "Particle.hpp" -#include "particle_node.hpp" #include #include @@ -237,7 +236,8 @@ void delete_particle_bond(int part, Utils::Span bond); */ void delete_particle_bonds(int part); -/** @brief Removs the specified bond from the particle +/** @brief Removes the specified bond from the particle + * @param p The particle to update * @param bond The bond in the form * {bond_id, partner_1, partner_2, ...} */ @@ -245,6 +245,8 @@ void local_remove_bond(Particle &p, std::vector const &bond); /** @brief Removes all pair bonds on the particle which have the specified * particle id as partner + * @param p The particle to update + * @param other_pid The particle id to filter for */ void local_remove_pair_bonds_to(Particle &p, int other_pid); diff --git a/src/core/particle_node.cpp b/src/core/particle_node.cpp index dfdf2b14539..8808063d19d 100644 --- a/src/core/particle_node.cpp +++ b/src/core/particle_node.cpp @@ -42,8 +42,6 @@ #include #include -#include - #include #include #include @@ -53,19 +51,16 @@ #include #include -enum { - /// ok code for \ref place_particle - ES_PART_OK = 0, - /// error code for \ref place_particle - ES_PART_ERROR = -1, - /// ok code for \ref place_particle, particle is new - ES_PART_CREATED = 1 -}; - constexpr auto some_tag = 42; + +/** @brief Enable particle type tracking in @ref particle_type_map. */ static bool type_list_enable; + +/** @brief Mapping particle types to lists of particle ids. */ static std::unordered_map> particle_type_map; -std::unordered_map particle_node; + +/** @brief Mapping particle ids to MPI ranks. */ +static std::unordered_map particle_node; void init_type_map(int type) { type_list_enable = true; @@ -80,13 +75,13 @@ void init_type_map(int type) { } } -void remove_id_from_map(int part_id, int type) { +static void remove_id_from_map(int part_id, int type) { auto it = particle_type_map.find(type); if (it != particle_type_map.end()) it->second.erase(part_id); } -void add_id_to_type_map(int part_id, int type) { +static void add_id_to_type_map(int part_id, int type) { auto it = particle_type_map.find(type); if (it != particle_type_map.end()) it->second.insert(part_id); @@ -115,7 +110,7 @@ Utils::Cache particle_fetch_cache(max_cache_size); void invalidate_fetch_cache() { particle_fetch_cache.invalidate(); } std::size_t fetch_cache_max_size() { return particle_fetch_cache.max_size(); } -boost::optional get_particle_data_local(int id) { +static boost::optional get_particle_data_local(int id) { auto p = cell_structure.get_local_particle(id); if (p and (not p->is_ghost())) { @@ -174,7 +169,7 @@ REGISTER_CALLBACK(mpi_get_particles_local) * * @returns The particle list. */ -std::vector mpi_get_particles(Utils::Span ids) { +static std::vector mpi_get_particles(Utils::Span ids) { mpi_call(mpi_get_particles_local); /* Return value */ std::vector parts(ids.size()); @@ -255,12 +250,12 @@ static void mpi_who_has_local() { std::transform(local_particles.begin(), local_particles.end(), sendbuf.begin(), [](Particle const &p) { return p.id(); }); - MPI_Send(sendbuf.data(), n_part, MPI_INT, 0, some_tag, comm_cart); + comm_cart.send(0, some_tag, sendbuf); } REGISTER_CALLBACK(mpi_who_has_local) -void mpi_who_has() { +static void mpi_who_has() { mpi_call(mpi_who_has_local); auto local_particles = cell_structure.local_particles(); @@ -279,8 +274,7 @@ void mpi_who_has() { } else if (n_parts[pnode] > 0) { pdata.resize(n_parts[pnode]); - MPI_Recv(pdata.data(), n_parts[pnode], MPI_INT, pnode, some_tag, - comm_cart, MPI_STATUS_IGNORE); + comm_cart.recv(pnode, some_tag, pdata); for (int i = 0; i < n_parts[pnode]; i++) particle_node[pdata[i]] = pnode; } @@ -290,11 +284,8 @@ void mpi_who_has() { /** * @brief Rebuild the particle index. */ -void build_particle_node() { mpi_who_has(); } +static void build_particle_node() { mpi_who_has(); } -/** - * @brief Get the mpi rank which owns the particle with id. - */ int get_particle_node(int id) { if (id < 0) { throw std::domain_error("Invalid particle id: " + std::to_string(id)); @@ -345,8 +336,8 @@ Particle *local_place_particle(int id, const Utils::Vector3d &pos, int _new) { return pt; } -boost::optional mpi_place_new_particle_local(int p_id, - Utils::Vector3d const &pos) { +static boost::optional mpi_place_new_particle_local( + int p_id, Utils::Vector3d const &pos) { auto p = local_place_particle(p_id, pos, 1); on_particle_change(); if (p) { @@ -362,7 +353,7 @@ REGISTER_CALLBACK_ONE_RANK(mpi_place_new_particle_local) * \param p_id the particle to create. * \param pos the particles position. */ -int mpi_place_new_particle(int p_id, const Utils::Vector3d &pos) { +static int mpi_place_new_particle(int p_id, const Utils::Vector3d &pos) { return mpi_call(Communication::Result::one_rank, mpi_place_new_particle_local, p_id, pos); } @@ -386,7 +377,7 @@ REGISTER_CALLBACK(mpi_place_particle_local) * \param p_id the particle to move. * \param pos the particles position. */ -void mpi_place_particle(int node, int p_id, const Utils::Vector3d &pos) { +static void mpi_place_particle(int node, int p_id, const Utils::Vector3d &pos) { mpi_call(mpi_place_particle_local, node, p_id); if (node == this_node) @@ -399,18 +390,15 @@ void mpi_place_particle(int node, int p_id, const Utils::Vector3d &pos) { on_particle_change(); } -int place_particle(int p_id, Utils::Vector3d const &pos) { +void place_particle(int p_id, Utils::Vector3d const &pos) { if (p_id < 0) { throw std::domain_error("Invalid particle id: " + std::to_string(p_id)); } if (particle_exists(p_id)) { mpi_place_particle(get_particle_node(p_id), p_id, pos); - - return ES_PART_OK; + } else { + particle_node[p_id] = mpi_place_new_particle(p_id, pos); } - particle_node[p_id] = mpi_place_new_particle(p_id, pos); - - return ES_PART_CREATED; } static void mpi_remove_particle_local(int p_id) { @@ -428,7 +416,7 @@ REGISTER_CALLBACK(mpi_remove_particle_local) * Also calls \ref on_particle_change. * \param p_id the particle to remove, use -1 to remove all particles. */ -void mpi_remove_particle(int p_id) { +static void mpi_remove_particle(int p_id) { mpi_call_all(mpi_remove_particle_local, p_id); } @@ -437,7 +425,7 @@ void remove_all_particles() { clear_particle_node(); } -int remove_particle(int p_id) { +void remove_particle(int p_id) { auto const &cur_par = get_particle_data(p_id); if (type_list_enable) { // remove particle from its current type_list @@ -448,8 +436,6 @@ int remove_particle(int p_id) { particle_node[p_id] = -1; mpi_remove_particle(p_id); particle_node.erase(p_id); - - return ES_OK; } int get_random_p_id(int type, int random_index_in_type_map) { diff --git a/src/core/particle_node.hpp b/src/core/particle_node.hpp index 733a74cdfbe..b87b7550f2a 100644 --- a/src/core/particle_node.hpp +++ b/src/core/particle_node.hpp @@ -40,8 +40,6 @@ * @brief Get particle data. * * @param part the identity of the particle to fetch - * @return Pointer to copy of particle if it exists, - * nullptr otherwise; */ const Particle &get_particle_data(int part); @@ -75,21 +73,16 @@ void clear_particle_node(); /** Call only on the head node. * Move a particle to a new position. * If it does not exist, it is created. - * @param part the identity of the particle to move - * @param p its new position - * @retval ES_PART_OK if particle existed - * @retval ES_PART_CREATED if created - * @retval ES_PART_ERROR if id is illegal + * @param p_id identity of the particle to move (or create) + * @param pos position */ -int place_particle(int part, Utils::Vector3d const &p); +void place_particle(int p_id, Utils::Vector3d const &pos); /** Remove particle with a given identity. Also removes all bonds to the * particle. * @param part identity of the particle to remove - * @retval ES_OK on success - * @retval ES_ERROR if particle does not exist */ -int remove_particle(int part); +void remove_particle(int part); /** Remove all particles. */ void remove_all_particles(); diff --git a/src/core/reaction_methods/ReactionAlgorithm.cpp b/src/core/reaction_methods/ReactionAlgorithm.cpp index decf32655b4..c7839ce1f5e 100644 --- a/src/core/reaction_methods/ReactionAlgorithm.cpp +++ b/src/core/reaction_methods/ReactionAlgorithm.cpp @@ -25,6 +25,7 @@ #include "grid.hpp" #include "partCfg_global.hpp" #include "particle_data.hpp" +#include "particle_node.hpp" #include "statistics.hpp" #include diff --git a/src/core/reaction_methods/tests/ConstantpHEnsemble_test.cpp b/src/core/reaction_methods/tests/ConstantpHEnsemble_test.cpp index b483424276d..b5aaab65764 100644 --- a/src/core/reaction_methods/tests/ConstantpHEnsemble_test.cpp +++ b/src/core/reaction_methods/tests/ConstantpHEnsemble_test.cpp @@ -29,7 +29,6 @@ #include "reaction_methods/utils.hpp" #include "communication.hpp" -#include "particle_data.hpp" #include diff --git a/src/core/reaction_methods/tests/ReactionAlgorithm_test.cpp b/src/core/reaction_methods/tests/ReactionAlgorithm_test.cpp index 6ffab4d522e..7ff6e2ecb27 100644 --- a/src/core/reaction_methods/tests/ReactionAlgorithm_test.cpp +++ b/src/core/reaction_methods/tests/ReactionAlgorithm_test.cpp @@ -31,6 +31,7 @@ #include "Particle.hpp" #include "communication.hpp" #include "particle_data.hpp" +#include "particle_node.hpp" #include diff --git a/src/core/reaction_methods/tests/ReactionEnsemble_test.cpp b/src/core/reaction_methods/tests/ReactionEnsemble_test.cpp index a50457d57b2..804f6c36a9c 100644 --- a/src/core/reaction_methods/tests/ReactionEnsemble_test.cpp +++ b/src/core/reaction_methods/tests/ReactionEnsemble_test.cpp @@ -31,6 +31,7 @@ #include "EspressoSystemStandAlone.hpp" #include "communication.hpp" #include "particle_data.hpp" +#include "particle_node.hpp" #include diff --git a/src/core/statistics_chain.cpp b/src/core/statistics_chain.cpp index 3aca1593094..53d392065a4 100644 --- a/src/core/statistics_chain.cpp +++ b/src/core/statistics_chain.cpp @@ -25,7 +25,7 @@ #include "Particle.hpp" #include "grid.hpp" -#include "particle_data.hpp" +#include "particle_node.hpp" #include #include diff --git a/src/core/unit_tests/EspressoSystemInterface_test.cpp b/src/core/unit_tests/EspressoSystemInterface_test.cpp index 61eace5d331..57039b073f5 100644 --- a/src/core/unit_tests/EspressoSystemInterface_test.cpp +++ b/src/core/unit_tests/EspressoSystemInterface_test.cpp @@ -27,6 +27,7 @@ #include "communication.hpp" #include "config.hpp" #include "particle_data.hpp" +#include "particle_node.hpp" #include "virtual_sites.hpp" #include "virtual_sites/VirtualSitesOff.hpp" diff --git a/src/core/unit_tests/EspressoSystemStandAlone_test.cpp b/src/core/unit_tests/EspressoSystemStandAlone_test.cpp index 4eb46a4761a..e1d34113c9c 100644 --- a/src/core/unit_tests/EspressoSystemStandAlone_test.cpp +++ b/src/core/unit_tests/EspressoSystemStandAlone_test.cpp @@ -42,6 +42,7 @@ namespace utf = boost::unit_test; #include "nonbonded_interactions/lj.hpp" #include "observables/ParticleVelocities.hpp" #include "particle_data.hpp" +#include "particle_node.hpp" #include #include diff --git a/src/core/unit_tests/ParticleFactory.hpp b/src/core/unit_tests/ParticleFactory.hpp index f1a66c00abb..072ed43484c 100644 --- a/src/core/unit_tests/ParticleFactory.hpp +++ b/src/core/unit_tests/ParticleFactory.hpp @@ -20,6 +20,7 @@ #define REACTION_ENSEMBLE_TESTS_PARTICLE_FACTORY_HPP #include "particle_data.hpp" +#include "particle_node.hpp" #include diff --git a/src/core/unit_tests/Verlet_list_test.cpp b/src/core/unit_tests/Verlet_list_test.cpp index de340f6f593..2c1431f4459 100644 --- a/src/core/unit_tests/Verlet_list_test.cpp +++ b/src/core/unit_tests/Verlet_list_test.cpp @@ -37,7 +37,7 @@ namespace bdata = boost::unit_test::data; #include "ParticleFactory.hpp" #include "integrate.hpp" #include "nonbonded_interactions/lj.hpp" -#include "particle_data.hpp" +#include "particle_node.hpp" #include "thermostat.hpp" #include diff --git a/src/core/virtual_sites.cpp b/src/core/virtual_sites.cpp index 629468a9b00..c78d2b318a5 100644 --- a/src/core/virtual_sites.cpp +++ b/src/core/virtual_sites.cpp @@ -32,6 +32,7 @@ #include "integrate.hpp" #include "nonbonded_interactions/nonbonded_interaction_data.hpp" #include "particle_data.hpp" +#include "particle_node.hpp" #include #include diff --git a/src/python/espressomd/particle_data.pxd b/src/python/espressomd/particle_data.pxd index a3baa3910c1..cbb714866d2 100644 --- a/src/python/espressomd/particle_data.pxd +++ b/src/python/espressomd/particle_data.pxd @@ -157,9 +157,9 @@ cdef extern from "particle_data.hpp": void remove_all_bonds_to(int part) cdef extern from "particle_node.hpp": - int place_particle(int part, const Vector3d & p) except + + void place_particle(int p_id, const Vector3d & pos) except + - int remove_particle(int part) except + + void remove_particle(int part) except + void remove_all_particles() except + diff --git a/src/python/espressomd/particle_data.pyx b/src/python/espressomd/particle_data.pyx index 02032b3b037..5ca0b14e7ad 100644 --- a/src/python/espressomd/particle_data.pyx +++ b/src/python/espressomd/particle_data.pyx @@ -170,8 +170,7 @@ cdef class ParticleHandle: raise ValueError("invalid particle position") check_type_or_throw_except( _pos, 3, float, "Position must be 3 floats") - if place_particle(self._id, make_Vector3d(_pos)) == -1: - raise Exception("particle could not be set") + place_particle(self._id, make_Vector3d(_pos)) def __get__(self): self.update_particle_data() @@ -1188,8 +1187,7 @@ cdef class ParticleHandle: espressomd.particle_data.ParticleList.clear """ - if remove_particle(self._id): - raise Exception("Could not remove particle.") + remove_particle(self._id) del self def add_verified_bond(self, bond): @@ -1795,9 +1793,7 @@ Set quat and scalar dipole moment (dipm) instead.") # done here. check_type_or_throw_except( p_dict["pos"], 3, float, "Position must be 3 floats.") - error_code = place_particle(p_dict["id"], make_Vector3d(p_dict["pos"])) - if error_code == -1: - raise Exception("particle could not be set.") + place_particle(p_dict["id"], make_Vector3d(p_dict["pos"])) # position is taken care of del p_dict["pos"] diff --git a/src/python/espressomd/system.pxd b/src/python/espressomd/system.pxd index b02f9410ed6..b587c67160f 100644 --- a/src/python/espressomd/system.pxd +++ b/src/python/espressomd/system.pxd @@ -33,7 +33,7 @@ IF EXCLUSIONS: cdef extern from "particle_data.hpp": void auto_exclusions(int distance) -cdef extern from "particle_data.hpp": +cdef extern from "particle_node.hpp": int init_type_map(int type) except + int number_of_particles_with_type(int type) except + diff --git a/testsuite/python/particle.py b/testsuite/python/particle.py index e812848bcac..b217ed8fa95 100644 --- a/testsuite/python/particle.py +++ b/testsuite/python/particle.py @@ -284,9 +284,11 @@ def test_invalid_particle_ids_exceptions(self): p._id = 42 p.node for i in range(1, 10): + p._id = -i with self.assertRaisesRegex(ValueError, f"Invalid particle id: {-i}"): - p._id = -i p.node + with self.assertRaisesRegex(ValueError, f"Invalid particle id: {-i}"): + p.remove() with self.assertRaisesRegex(ValueError, f"Invalid particle id: {-i}"): self.system.part.add(pos=[0., 0., 0.], id=-i) From b388ce5d941afaf5ac9bc80ebe57a6f364b3afd8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Wed, 6 Apr 2022 18:09:02 +0200 Subject: [PATCH 3/7] core: Use consistent argument names --- src/core/particle_node.cpp | 62 ++++++++++++------------- src/core/particle_node.hpp | 16 +++---- src/python/espressomd/particle_data.pxd | 8 ++-- 3 files changed, 43 insertions(+), 43 deletions(-) diff --git a/src/core/particle_node.cpp b/src/core/particle_node.cpp index 8808063d19d..d4bd02b8006 100644 --- a/src/core/particle_node.cpp +++ b/src/core/particle_node.cpp @@ -75,16 +75,16 @@ void init_type_map(int type) { } } -static void remove_id_from_map(int part_id, int type) { +static void remove_id_from_map(int p_id, int type) { auto it = particle_type_map.find(type); if (it != particle_type_map.end()) - it->second.erase(part_id); + it->second.erase(p_id); } -static void add_id_to_type_map(int part_id, int type) { +static void add_id_to_type_map(int p_id, int type) { auto it = particle_type_map.find(type); if (it != particle_type_map.end()) - it->second.insert(part_id); + it->second.insert(p_id); } void on_particle_type_change(int p_id, int type) { @@ -110,8 +110,8 @@ Utils::Cache particle_fetch_cache(max_cache_size); void invalidate_fetch_cache() { particle_fetch_cache.invalidate(); } std::size_t fetch_cache_max_size() { return particle_fetch_cache.max_size(); } -static boost::optional get_particle_data_local(int id) { - auto p = cell_structure.get_local_particle(id); +static boost::optional get_particle_data_local(int p_id) { + auto p = cell_structure.get_local_particle(p_id); if (p and (not p->is_ghost())) { return *p; @@ -122,16 +122,16 @@ static boost::optional get_particle_data_local(int id) { REGISTER_CALLBACK_ONE_RANK(get_particle_data_local) -const Particle &get_particle_data(int part) { - auto const pnode = get_particle_node(part); +const Particle &get_particle_data(int p_id) { + auto const pnode = get_particle_node(p_id); if (pnode == this_node) { - assert(cell_structure.get_local_particle(part)); - return *cell_structure.get_local_particle(part); + assert(cell_structure.get_local_particle(p_id)); + return *cell_structure.get_local_particle(p_id); } /* Query the cache */ - auto const p_ptr = particle_fetch_cache.get(part); + auto const p_ptr = particle_fetch_cache.get(p_id); if (p_ptr) { return *p_ptr; } @@ -139,8 +139,8 @@ const Particle &get_particle_data(int part) { /* Cache miss, fetch the particle, * put it into the cache and return a pointer into the cache. */ auto const cache_ptr = particle_fetch_cache.put( - part, Communication::mpiCallbacks().call(Communication::Result::one_rank, - get_particle_data_local, part)); + p_id, Communication::mpiCallbacks().call(Communication::Result::one_rank, + get_particle_data_local, p_id)); return *cache_ptr; } @@ -180,11 +180,11 @@ static std::vector mpi_get_particles(Utils::Span ids) { per_node.clear(); } - for (auto const &id : ids) { - auto const pnode = get_particle_node(id); + for (auto const &p_id : ids) { + auto const p_node = get_particle_node(p_id); - if (pnode != this_node) - node_ids[pnode].push_back(id); + if (p_node != this_node) + node_ids[p_node].push_back(p_id); } /* Distributed ids to the nodes */ @@ -195,9 +195,9 @@ static std::vector mpi_get_particles(Utils::Span ids) { /* Copy local particles */ std::transform(node_ids[this_node].cbegin(), node_ids[this_node].cend(), - parts.begin(), [](int id) { - assert(cell_structure.get_local_particle(id)); - return *cell_structure.get_local_particle(id); + parts.begin(), [](int p_id) { + assert(cell_structure.get_local_particle(p_id)); + return *cell_structure.get_local_particle(p_id); }); static std::vector node_sizes(comm_cart.size()); @@ -286,19 +286,19 @@ static void mpi_who_has() { */ static void build_particle_node() { mpi_who_has(); } -int get_particle_node(int id) { - if (id < 0) { - throw std::domain_error("Invalid particle id: " + std::to_string(id)); +int get_particle_node(int p_id) { + if (p_id < 0) { + throw std::domain_error("Invalid particle id: " + std::to_string(p_id)); } if (particle_node.empty()) build_particle_node(); - auto const needle = particle_node.find(id); + auto const needle = particle_node.find(p_id); // Check if particle has a node, if not, we assume it does not exist. if (needle == particle_node.end()) { - throw std::runtime_error("Particle node for id " + std::to_string(id) + + throw std::runtime_error("Particle node for id " + std::to_string(p_id) + " not found!"); } return needle->second; @@ -309,27 +309,27 @@ void clear_particle_node() { particle_node.clear(); } /** Move a particle to a new position. If it does not exist, it is created. * The position must be on the local node! * - * @param id the identity of the particle to move + * @param p_id the identity of the particle to move * @param pos its new position * @param _new if true, the particle is allocated, else has to exists already * * @return Pointer to the particle. */ -Particle *local_place_particle(int id, const Utils::Vector3d &pos, int _new) { +Particle *local_place_particle(int p_id, const Utils::Vector3d &pos, int _new) { auto pp = Utils::Vector3d{pos[0], pos[1], pos[2]}; auto i = Utils::Vector3i{}; fold_position(pp, i, box_geo); if (_new) { Particle new_part; - new_part.id() = id; + new_part.id() = p_id; new_part.pos() = pp; new_part.image_box() = i; return cell_structure.add_local_particle(std::move(new_part)); } - auto pt = cell_structure.get_local_particle(id); + auto pt = cell_structure.get_local_particle(p_id); pt->pos() = pp; pt->image_box() = i; @@ -463,10 +463,10 @@ int number_of_particles_with_type(int type) { return static_cast(it->second.size()); } -bool particle_exists(int part_id) { +bool particle_exists(int p_id) { if (particle_node.empty()) build_particle_node(); - return particle_node.count(part_id); + return particle_node.count(p_id); } std::vector get_particle_ids() { diff --git a/src/core/particle_node.hpp b/src/core/particle_node.hpp index b87b7550f2a..0ee27247450 100644 --- a/src/core/particle_node.hpp +++ b/src/core/particle_node.hpp @@ -41,7 +41,7 @@ * * @param part the identity of the particle to fetch */ -const Particle &get_particle_data(int part); +const Particle &get_particle_data(int p_id); /** * @brief Fetch a range of particle into the fetch cache. @@ -80,9 +80,9 @@ void place_particle(int p_id, Utils::Vector3d const &pos); /** Remove particle with a given identity. Also removes all bonds to the * particle. - * @param part identity of the particle to remove + * @param p_id identity of the particle to remove */ -void remove_particle(int part); +void remove_particle(int p_id); /** Remove all particles. */ void remove_all_particles(); @@ -97,18 +97,18 @@ int number_of_particles_with_type(int type); /** * @brief Check if particle exists. * - * @param part Id of the particle + * @param p_id identity of the particle * @return True iff the particle exists. */ -bool particle_exists(int part); +bool particle_exists(int p_id); /** - * @brief Get the mpi rank which owns the particle with id. + * @brief Get the MPI rank which owns the a specific particle. * - * @param id Id of the particle + * @param p_id identity of the particle * @return The MPI rank the particle is on. */ -int get_particle_node(int id); +int get_particle_node(int p_id); /** * @brief Get all particle ids. diff --git a/src/python/espressomd/particle_data.pxd b/src/python/espressomd/particle_data.pxd index cbb714866d2..5ee1ee45b4f 100644 --- a/src/python/espressomd/particle_data.pxd +++ b/src/python/espressomd/particle_data.pxd @@ -159,15 +159,15 @@ cdef extern from "particle_data.hpp": cdef extern from "particle_node.hpp": void place_particle(int p_id, const Vector3d & pos) except + - void remove_particle(int part) except + + void remove_particle(int p_id) except + void remove_all_particles() except + - bool particle_exists(int part) + bool particle_exists(int p_id) - int get_particle_node(int pid) except + + int get_particle_node(int p_id) except + - const particle & get_particle_data(int pid) except + + const particle & get_particle_data(int p_id) except + vector[int] get_particle_ids() except + From 5a4c1c42d3fb80d41775d22b6587ad2dd47742e8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Wed, 6 Apr 2022 19:58:54 +0200 Subject: [PATCH 4/7] core: Simplify particle callbacks Split callbacks that mixed multiple responsibilities via a flag into multiple callbacks that have only one responsibility each. --- src/core/exclusions.cpp | 17 +++-- src/core/exclusions.hpp | 8 +- src/core/particle_data.cpp | 99 ++++++++++++++----------- src/core/particle_data.hpp | 28 +++---- src/core/particle_node.cpp | 86 +++++++++++---------- src/python/espressomd/particle_data.pxd | 3 +- src/python/espressomd/particle_data.pyx | 6 +- 7 files changed, 135 insertions(+), 112 deletions(-) diff --git a/src/core/exclusions.cpp b/src/core/exclusions.cpp index 30626bb3ea1..51effe844ee 100644 --- a/src/core/exclusions.cpp +++ b/src/core/exclusions.cpp @@ -25,20 +25,23 @@ #include "exclusions.hpp" +#include "Particle.hpp" + #include #include -void add_exclusion(Particle *part, int part2) { - if (Utils::contains(part->exclusions(), part2)) +void add_exclusion(Particle &p, int p_id) { + if (Utils::contains(p.exclusions(), p_id)) return; - part->exclusions().push_back(part2); + p.exclusions().push_back(p_id); } -void delete_exclusion(Particle *part, int part2) { - auto &el = part->exclusions(); +void delete_exclusion(Particle &p, int p_id) { + auto &el = p.exclusions(); - el.erase(std::remove(el.begin(), el.end(), part2), el.end()); + el.erase(std::remove(el.begin(), el.end(), p_id), el.end()); } -#endif + +#endif // EXCLUSIONS diff --git a/src/core/exclusions.hpp b/src/core/exclusions.hpp index 67b4f7439b0..a199d6f960a 100644 --- a/src/core/exclusions.hpp +++ b/src/core/exclusions.hpp @@ -28,6 +28,7 @@ #include #ifdef EXCLUSIONS + /** Determine if the non-bonded interactions between @p p1 and @p p2 should be * calculated. */ @@ -39,9 +40,10 @@ inline bool do_nonbonded(Particle const &p1, Particle const &p2) { } /** Remove exclusion from particle if possible */ -void delete_exclusion(Particle *part, int part2); +void delete_exclusion(Particle &p, int p_id); /** Insert an exclusion if not already set */ -void add_exclusion(Particle *part, int part2); -#endif +void add_exclusion(Particle &p, int p_id); + +#endif // EXCLUSIONS #endif // ESPRESSO_EXCLUSIONS_HPP diff --git a/src/core/particle_data.cpp b/src/core/particle_data.cpp index 5b38e348779..97d4df0f6ce 100644 --- a/src/core/particle_data.cpp +++ b/src/core/particle_data.cpp @@ -26,6 +26,7 @@ #include "communication.hpp" #include "config.hpp" #include "event.hpp" +#include "exclusions.hpp" #include "nonbonded_interactions/nonbonded_interaction_data.hpp" #include "partCfg_global.hpp" #include "particle_node.hpp" @@ -42,6 +43,7 @@ #include #include +#include #include #include #include @@ -400,12 +402,6 @@ void mpi_update_particle_property(int id, const T &value) { mpi_update_particle(id, value); } -void delete_exclusion(Particle *part, int part2); - -void add_exclusion(Particle *part, int part2); - -void auto_exclusion(int distance); - void set_particle_v(int part, Utils::Vector3d const &v) { mpi_update_particle(part, v); @@ -683,35 +679,41 @@ void mpi_rescale_particles(int dir, double scale) { } #ifdef EXCLUSIONS -/** Locally add an exclusion to a particle. - * @param part1 the identity of the first exclusion partner - * @param part2 the identity of the second exclusion partner - * @param _delete if true, delete the exclusion instead of add +/** + * @brief Locally remove an exclusion to a particle. + * @param part1 the identity of the first exclusion partner + * @param part2 the identity of the second exclusion partner */ -static void local_change_exclusion(int part1, int part2, int _delete) { - /* part1, if here */ - auto part = cell_structure.get_local_particle(part1); - if (part) { - if (_delete) - delete_exclusion(part, part2); - else - add_exclusion(part, part2); +static void local_remove_exclusion(int part1, int part2) { + auto *p1 = cell_structure.get_local_particle(part1); + if (p1) { + delete_exclusion(*p1, part2); + } + auto *p2 = cell_structure.get_local_particle(part2); + if (p2) { + delete_exclusion(*p2, part1); } +} - /* part2, if here */ - part = cell_structure.get_local_particle(part2); - if (part) { - if (_delete) - delete_exclusion(part, part1); - else - add_exclusion(part, part1); +/** + * @brief Locally add an exclusion to a particle. + * @param part1 the identity of the first exclusion partner + * @param part2 the identity of the second exclusion partner + */ +static void local_add_exclusion(int part1, int part2) { + auto *p1 = cell_structure.get_local_particle(part1); + if (p1) { + add_exclusion(*p1, part2); + } + auto *p2 = cell_structure.get_local_particle(part2); + if (p2) { + add_exclusion(*p2, part1); } } -namespace { /* keep a unique list for particle i. Particle j is only added if it is not i and not already in the list. */ -void add_partner(std::vector &il, int i, int j, int distance) { +static void add_partner(std::vector &il, int i, int j, int distance) { if (j == i) return; for (int k = 0; k < il.size(); k += 2) @@ -721,31 +723,38 @@ void add_partner(std::vector &il, int i, int j, int distance) { il.push_back(j); il.push_back(distance); } -} // namespace -static void mpi_send_exclusion_local(int part1, int part2, int _delete) { - local_change_exclusion(part1, part2, _delete); +static void mpi_remove_exclusion_local(int part1, int part2) { + local_remove_exclusion(part1, part2); on_particle_change(); } -REGISTER_CALLBACK(mpi_send_exclusion_local) +REGISTER_CALLBACK(mpi_remove_exclusion_local) -/** Send exclusions. - * Also calls \ref on_particle_change. - * \param part1 identity of first particle of the exclusion. - * \param part2 identity of second particle of the exclusion. - * \param _delete if true, do not add the exclusion, rather delete it if found - */ -static void mpi_send_exclusion(int part1, int part2, int _delete) { - mpi_call_all(mpi_send_exclusion_local, part1, part2, _delete); +static void mpi_add_exclusion_local(int part1, int part2) { + local_add_exclusion(part1, part2); + on_particle_change(); } -int change_exclusion(int part1, int part2, int _delete) { - if (particle_exists(part1) && particle_exists(part2)) { - mpi_send_exclusion(part1, part2, _delete); - return ES_OK; +REGISTER_CALLBACK(mpi_add_exclusion_local) + +static void check_particle_exists(int p_id) { + if (not particle_exists(p_id)) { + throw std::runtime_error("Particle with id " + std::to_string(p_id) + + " not found"); } - return ES_ERROR; +} + +void remove_particle_exclusion(int part1, int part2) { + check_particle_exists(part1); + check_particle_exists(part2); + mpi_call_all(mpi_remove_exclusion_local, part1, part2); +} + +void add_particle_exclusion(int part1, int part2) { + check_particle_exists(part1); + check_particle_exists(part2); + mpi_call_all(mpi_add_exclusion_local, part1, part2); } void auto_exclusions(int distance) { @@ -800,7 +809,7 @@ void auto_exclusions(int distance) { auto const partner_list = kv.second; for (int j : partner_list) if (id < j) - change_exclusion(id, j, 0); + add_particle_exclusion(id, j); } } #endif // EXCLUSIONS diff --git a/src/core/particle_data.hpp b/src/core/particle_data.hpp index c65eafbb8a2..748cef0cdc3 100644 --- a/src/core/particle_data.hpp +++ b/src/core/particle_data.hpp @@ -260,21 +260,19 @@ void add_particle_bond(int part, Utils::Span bond); const std::vector &get_particle_bonds(int part); #ifdef EXCLUSIONS -/** Call only on the head node: change particle constraints. - * @param part identity of particle for which the exclusion is set. +/** @brief Remove particle exclusion. + * Call only on the head node. + * @param part1 identity of particle for which the exclusion is set. * @param part2 identity of particle for which the exclusion is set. - * If -1, delete all exclusions. - * @param _delete if true, do not add the exclusion, rather delete it if - * found - * @retval ES_OK on success - * @retval ES_ERROR on failure (e.g. particles do not exist / did not have - * exclusion set) - */ -int change_exclusion(int part, int part2, int _delete); -#endif + */ +void remove_particle_exclusion(int part1, int part2); -/** Rescale all particle positions in direction @p dir by a factor @p scale. */ -void mpi_rescale_particles(int dir, double scale); +/** @brief Add particle exclusion. + * Call only on the head node. + * @param part1 identity of particle for which the exclusion is set. + * @param part2 identity of particle for which the exclusion is set. + */ +void add_particle_exclusion(int part1, int part2); /** Automatically add the next \ neighbors in each molecule to the * exclusion list. @@ -284,6 +282,10 @@ void mpi_rescale_particles(int dir, double scale); * particles, you should avoid this function and setup exclusions manually. */ void auto_exclusions(int distance); +#endif + +/** Rescale all particle positions in direction @p dir by a factor @p scale. */ +void mpi_rescale_particles(int dir, double scale); // The following functions are used by the python interface to obtain // properties of a particle, which are only compiled in in some configurations diff --git a/src/core/particle_node.cpp b/src/core/particle_node.cpp index d4bd02b8006..a3d54cd495d 100644 --- a/src/core/particle_node.cpp +++ b/src/core/particle_node.cpp @@ -306,39 +306,52 @@ int get_particle_node(int p_id) { void clear_particle_node() { particle_node.clear(); } -/** Move a particle to a new position. If it does not exist, it is created. - * The position must be on the local node! +/** + * @brief Create a new particle. + * The position must be on the local node! * - * @param p_id the identity of the particle to move - * @param pos its new position - * @param _new if true, the particle is allocated, else has to exists already + * @param p_id identity of the particle to create + * @param pos position * - * @return Pointer to the particle. + * @return Pointer to the particle. */ -Particle *local_place_particle(int p_id, const Utils::Vector3d &pos, int _new) { - auto pp = Utils::Vector3d{pos[0], pos[1], pos[2]}; - auto i = Utils::Vector3i{}; - fold_position(pp, i, box_geo); - - if (_new) { - Particle new_part; - new_part.id() = p_id; - new_part.pos() = pp; - new_part.image_box() = i; - - return cell_structure.add_local_particle(std::move(new_part)); - } +static Particle *local_insert_particle(int p_id, const Utils::Vector3d &pos) { + auto folded_pos = pos; + auto image_box = Utils::Vector3i{}; + fold_position(folded_pos, image_box, box_geo); + + Particle new_part; + new_part.id() = p_id; + new_part.pos() = folded_pos; + new_part.image_box() = image_box; + + return cell_structure.add_local_particle(std::move(new_part)); +} + +/** + * @brief Move a particle to a new position. + * The position must be on the local node! + * + * @param p_id identity of the particle to move + * @param pos new position + * + * @return Pointer to the particle. + */ +static Particle *local_move_particle(int p_id, const Utils::Vector3d &pos) { + auto folded_pos = pos; + auto image_box = Utils::Vector3i{}; + fold_position(folded_pos, image_box, box_geo); auto pt = cell_structure.get_local_particle(p_id); - pt->pos() = pp; - pt->image_box() = i; + pt->pos() = folded_pos; + pt->image_box() = image_box; return pt; } -static boost::optional mpi_place_new_particle_local( - int p_id, Utils::Vector3d const &pos) { - auto p = local_place_particle(p_id, pos, 1); +static boost::optional +mpi_place_new_particle_local(int p_id, Utils::Vector3d const &pos) { + auto p = local_insert_particle(p_id, pos); on_particle_change(); if (p) { return comm_cart.rank(); @@ -362,7 +375,7 @@ void mpi_place_particle_local(int pnode, int p_id) { if (pnode == this_node) { Utils::Vector3d pos; comm_cart.recv(0, some_tag, pos); - local_place_particle(p_id, pos, 0); + local_move_particle(p_id, pos); } cell_structure.set_resort_particles(Cells::RESORT_GLOBAL); @@ -381,7 +394,7 @@ static void mpi_place_particle(int node, int p_id, const Utils::Vector3d &pos) { mpi_call(mpi_place_particle_local, node, p_id); if (node == this_node) - local_place_particle(p_id, pos, 0); + local_move_particle(p_id, pos); else { comm_cart.send(node, some_tag, pos); } @@ -402,26 +415,21 @@ void place_particle(int p_id, Utils::Vector3d const &pos) { } static void mpi_remove_particle_local(int p_id) { - if (p_id != -1) { - cell_structure.remove_particle(p_id); - } else { - cell_structure.remove_all_particles(); - } + cell_structure.remove_particle(p_id); on_particle_change(); } REGISTER_CALLBACK(mpi_remove_particle_local) -/** Remove a particle. - * Also calls \ref on_particle_change. - * \param p_id the particle to remove, use -1 to remove all particles. - */ -static void mpi_remove_particle(int p_id) { - mpi_call_all(mpi_remove_particle_local, p_id); +static void mpi_remove_all_particles_local() { + cell_structure.remove_all_particles(); + on_particle_change(); } +REGISTER_CALLBACK(mpi_remove_all_particles_local) + void remove_all_particles() { - mpi_remove_particle(-1); + mpi_call_all(mpi_remove_all_particles_local); clear_particle_node(); } @@ -434,7 +442,7 @@ void remove_particle(int p_id) { } particle_node[p_id] = -1; - mpi_remove_particle(p_id); + mpi_call_all(mpi_remove_particle_local, p_id); particle_node.erase(p_id); } diff --git a/src/python/espressomd/particle_data.pxd b/src/python/espressomd/particle_data.pxd index 5ee1ee45b4f..6efeff9019e 100644 --- a/src/python/espressomd/particle_data.pxd +++ b/src/python/espressomd/particle_data.pxd @@ -149,7 +149,8 @@ cdef extern from "particle_data.hpp": const vector[BondView] & get_particle_bonds(int part) IF EXCLUSIONS: - int change_exclusion(int part, int part2, int _delete) + void remove_particle_exclusion(int part1, int part2) + void add_particle_exclusion(int part1, int part2) IF ENGINE: void set_particle_swimming(int part, particle_parameters_swimming swim) diff --git a/src/python/espressomd/particle_data.pyx b/src/python/espressomd/particle_data.pyx index 5ca0b14e7ad..bcfc26166e1 100644 --- a/src/python/espressomd/particle_data.pyx +++ b/src/python/espressomd/particle_data.pyx @@ -1036,8 +1036,7 @@ cdef class ParticleHandle: raise Exception( "Cannot exclude of a particle with itself!\n" f"->particle id {self._id}, partner {_partner}.") - if change_exclusion(self._id, _partner, 0) == 1: - raise Exception(f"Particle with id {_partner} does not exist.") + add_particle_exclusion(self._id, _partner) def delete_exclusion(self, _partner): check_type_or_throw_except( @@ -1046,8 +1045,7 @@ cdef class ParticleHandle: if not self.particle_data.has_exclusion(_partner): raise Exception( f"Particle with id {_partner} is not in exclusion list.") - if change_exclusion(self._id, _partner, 1) == 1: - raise Exception(f"Particle with id {_partner} does not exist.") + remove_particle_exclusion(self._id, _partner) IF ENGINE: property swimming: From 255d0c82aa69f1951ec45a35ecc2e6573d302ec7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Wed, 6 Apr 2022 23:58:57 +0200 Subject: [PATCH 5/7] core: Rewrite particle exclusion exception mechanism --- src/core/particle_data.cpp | 13 +++-- src/python/espressomd/particle_data.pxd | 4 +- src/python/espressomd/particle_data.pyx | 69 ++++++++++++++++--------- testsuite/python/particle.py | 18 +++++++ 4 files changed, 75 insertions(+), 29 deletions(-) diff --git a/src/core/particle_data.cpp b/src/core/particle_data.cpp index 97d4df0f6ce..f3247dbca91 100644 --- a/src/core/particle_data.cpp +++ b/src/core/particle_data.cpp @@ -745,15 +745,22 @@ static void check_particle_exists(int p_id) { } } -void remove_particle_exclusion(int part1, int part2) { +static void particle_exclusion_sanity_checks(int part1, int part2) { + if (part1 == part2) { + throw std::runtime_error("Particles cannot exclude themselves (id " + + std::to_string(part1) + ")"); + } check_particle_exists(part1); check_particle_exists(part2); +} + +void remove_particle_exclusion(int part1, int part2) { + particle_exclusion_sanity_checks(part1, part2); mpi_call_all(mpi_remove_exclusion_local, part1, part2); } void add_particle_exclusion(int part1, int part2) { - check_particle_exists(part1); - check_particle_exists(part2); + particle_exclusion_sanity_checks(part1, part2); mpi_call_all(mpi_add_exclusion_local, part1, part2); } diff --git a/src/python/espressomd/particle_data.pxd b/src/python/espressomd/particle_data.pxd index 6efeff9019e..7da8dd4abd9 100644 --- a/src/python/espressomd/particle_data.pxd +++ b/src/python/espressomd/particle_data.pxd @@ -149,8 +149,8 @@ cdef extern from "particle_data.hpp": const vector[BondView] & get_particle_bonds(int part) IF EXCLUSIONS: - void remove_particle_exclusion(int part1, int part2) - void add_particle_exclusion(int part1, int part2) + void remove_particle_exclusion(int part1, int part2) except + + void add_particle_exclusion(int part1, int part2) except + IF ENGINE: void set_particle_swimming(int part, particle_parameters_swimming swim) diff --git a/src/python/espressomd/particle_data.pyx b/src/python/espressomd/particle_data.pyx index bcfc26166e1..ee75e5edc79 100644 --- a/src/python/espressomd/particle_data.pyx +++ b/src/python/espressomd/particle_data.pyx @@ -989,24 +989,26 @@ cdef class ParticleHandle: IF EXCLUSIONS: property exclusions: """ - The exclusion list of particles where nonbonded interactions are ignored. + The exclusion list of particles where non-bonded interactions are ignored. .. note:: This needs the feature ``EXCLUSIONS``. + exclusions : (N,) array_like of :obj:`int` + """ - def __set__(self, _partners): + def __set__(self, partners): # Delete all for e in self.exclusions: self.delete_exclusion(e) - nlvl = nesting_level(_partners) + nlvl = nesting_level(partners) if nlvl == 0: # Single item - self.add_exclusion(_partners) + self.add_exclusion(partners) elif nlvl == 1: # List of items - for partner in _partners: + for partner in partners: self.add_exclusion(partner) else: raise ValueError( @@ -1016,36 +1018,55 @@ cdef class ParticleHandle: self.update_particle_data() return array_locked(self.particle_data.exclusions_as_vector()) - def add_exclusion(self, _partner): + def add_exclusion(self, partner): """ - Excluding interaction with the given partner. + Exclude non-bonded interactions with the given partner. + + .. note:: + This needs the feature ``EXCLUSIONS``. Parameters ----------- - _partner : :obj:`int` - partner + partner : :class:`~espressomd.particle_data.ParticleHandle` or :obj:`int` + Particle to exclude. """ + if isinstance(partner, ParticleHandle): + p_id = partner.id + else: + p_id = partner check_type_or_throw_except( - _partner, 1, int, "PID of partner has to be an int.") + p_id, 1, int, "Argument 'partner' has to be a ParticleHandle or int.") self.update_particle_data() - if self.particle_data.has_exclusion(_partner): - raise Exception( - f"Exclusion id {_partner} already in exclusion list of particle {self._id}") - if self._id == _partner: - raise Exception( - "Cannot exclude of a particle with itself!\n" - f"->particle id {self._id}, partner {_partner}.") - add_particle_exclusion(self._id, _partner) + if self.particle_data.has_exclusion(p_id): + raise RuntimeError( + f"Particle with id {p_id} is already in exclusion list of particle with id {self._id}") + add_particle_exclusion(self._id, p_id) - def delete_exclusion(self, _partner): + def delete_exclusion(self, partner): + """ + Remove exclusion of non-bonded interactions with the given partner. + + .. note:: + This needs the feature ``EXCLUSIONS``. + + Parameters + ----------- + partner : :class:`~espressomd.particle_data.ParticleHandle` or :obj:`int` + Particle to remove from exclusions. + + """ + if isinstance(partner, ParticleHandle): + p_id = partner.id + else: + p_id = partner check_type_or_throw_except( - _partner, 1, int, "PID of partner has to be an int.") + p_id, 1, int, "Argument 'partner' has to be a ParticleHandle or int.") self.update_particle_data() - if not self.particle_data.has_exclusion(_partner): - raise Exception( - f"Particle with id {_partner} is not in exclusion list.") - remove_particle_exclusion(self._id, _partner) + if not self.particle_data.has_exclusion(p_id): + raise RuntimeError( + f"Particle with id {p_id} is not in exclusion list of particle with id {self._id}") + remove_particle_exclusion(self._id, p_id) IF ENGINE: property swimming: diff --git a/testsuite/python/particle.py b/testsuite/python/particle.py index b217ed8fa95..9f5a198f56c 100644 --- a/testsuite/python/particle.py +++ b/testsuite/python/particle.py @@ -210,6 +210,24 @@ def test_vs_relative(self): with self.assertRaisesRegex(ValueError, "quaternion is zero"): p1.vs_quat = [0, 0, 0, 0] + @utx.skipIfMissingFeatures(["EXCLUSIONS"]) + def test_invalid_exclusions(self): + pid1 = self.pid + pid2 = self.pid + 1 + + p1 = self.partcl + with self.assertRaisesRegex(RuntimeError, rf"Particles cannot exclude themselves \(id {self.pid}\)"): + p1.add_exclusion(pid1) + with self.assertRaisesRegex(RuntimeError, f"Particle with id {pid2} not found"): + p1.add_exclusion(pid2) + + self.system.part.add(id=pid2, pos=(0, 0, 0)) + with self.assertRaisesRegex(RuntimeError, f"Particle with id {pid2} is not in exclusion list of particle with id {pid1}"): + p1.delete_exclusion(pid2) + with self.assertRaisesRegex(RuntimeError, f"Particle with id {pid2} is already in exclusion list of particle with id {pid1}"): + p1.add_exclusion(pid2) + p1.add_exclusion(pid2) + @utx.skipIfMissingFeatures("DIPOLES") def test_contradicting_properties_dip_dipm(self): with self.assertRaises(ValueError): From cf061d778e3f46d68f099fd3528d3c4c501bba59 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Thu, 7 Apr 2022 19:01:48 +0200 Subject: [PATCH 6/7] core: Make particle insertion run in constant time --- src/core/particle_node.cpp | 49 ++++++++++++++++++++++++++++++------ src/core/particle_node.hpp | 2 +- testsuite/python/particle.py | 38 +++++++++++++++++++++++++++- 3 files changed, 80 insertions(+), 9 deletions(-) diff --git a/src/core/particle_node.cpp b/src/core/particle_node.cpp index a3d54cd495d..380a0a18307 100644 --- a/src/core/particle_node.cpp +++ b/src/core/particle_node.cpp @@ -62,6 +62,15 @@ static std::unordered_map> particle_type_map; /** @brief Mapping particle ids to MPI ranks. */ static std::unordered_map particle_node; +/** + * @brief Keep track of the largest particle id. + * This book-keeping variable is necessary to make particle insertion run + * in constant time. Traversing the @ref particle_node to find the largest + * particle id scales with O(N) and traversing the local cells in parallel + * followed by a reduction scales with O(N^2). + */ +static int max_seen_pid = -1; + void init_type_map(int type) { type_list_enable = true; if (type < 0) @@ -265,18 +274,22 @@ static void mpi_who_has() { n_parts, 0); static std::vector pdata; + max_seen_pid = -1; /* then fetch particle locations */ for (int pnode = 0; pnode < n_nodes; pnode++) { if (pnode == this_node) { - for (auto const &p : local_particles) + for (auto const &p : local_particles) { particle_node[p.id()] = this_node; - + max_seen_pid = std::max(max_seen_pid, p.id()); + } } else if (n_parts[pnode] > 0) { pdata.resize(n_parts[pnode]); comm_cart.recv(pnode, some_tag, pdata); - for (int i = 0; i < n_parts[pnode]; i++) + for (int i = 0; i < n_parts[pnode]; i++) { particle_node[pdata[i]] = pnode; + max_seen_pid = std::max(max_seen_pid, pdata[i]); + } } } } @@ -306,6 +319,21 @@ int get_particle_node(int p_id) { void clear_particle_node() { particle_node.clear(); } +/** + * @brief Calculate the largest particle id. + * Traversing the @ref particle_node to find the largest particle id + * scales with O(N). Consider using the cached value in @ref max_seen_pid + * if possible. This function is only necessary when the cached value is + * invalidated, for example when removing the particle which has the + * largest id. + */ +static int calculate_max_seen_id() { + return boost::accumulate(particle_node, -1, + [](int max, const std::pair &kv) { + return std::max(max, kv.first); + }); +} + /** * @brief Create a new particle. * The position must be on the local node! @@ -411,6 +439,7 @@ void place_particle(int p_id, Utils::Vector3d const &pos) { mpi_place_particle(get_particle_node(p_id), p_id, pos); } else { particle_node[p_id] = mpi_place_new_particle(p_id, pos); + max_seen_pid = std::max(max_seen_pid, p_id); } } @@ -444,6 +473,15 @@ void remove_particle(int p_id) { particle_node[p_id] = -1; mpi_call_all(mpi_remove_particle_local, p_id); particle_node.erase(p_id); + if (p_id == max_seen_pid) { + --max_seen_pid; + // if there is a gap (i.e. there is no particle with id max_seen_pid - 1, + // then the cached value is invalidated and has to be recomputed (slow) + if (particle_node.count(max_seen_pid) == 0 or + particle_node[max_seen_pid] == -1) { + max_seen_pid = calculate_max_seen_id(); + } + } } int get_random_p_id(int type, int random_index_in_type_map) { @@ -491,10 +529,7 @@ int get_maximal_particle_id() { if (particle_node.empty()) build_particle_node(); - return boost::accumulate(particle_node, -1, - [](int max, const std::pair &kv) { - return std::max(max, kv.first); - }); + return max_seen_pid; } int get_n_part() { diff --git a/src/core/particle_node.hpp b/src/core/particle_node.hpp index 0ee27247450..b951714bb63 100644 --- a/src/core/particle_node.hpp +++ b/src/core/particle_node.hpp @@ -39,7 +39,7 @@ /** * @brief Get particle data. * - * @param part the identity of the particle to fetch + * @param p_id the identity of the particle to fetch */ const Particle &get_particle_data(int p_id); diff --git a/testsuite/python/particle.py b/testsuite/python/particle.py index 9f5a198f56c..dca864992f3 100644 --- a/testsuite/python/particle.py +++ b/testsuite/python/particle.py @@ -292,6 +292,42 @@ def test_image_box(self): np.testing.assert_equal(np.copy(p.image_box), [1, 1, 1]) + def test_particle_numbering(self): + """ + Instantiate particles on different nodes and check they are + numbered sequentially. + """ + system = self.system + offset = system.box_l / 100. + # clear the cached value for largest particle id + system.part.clear() + # update the cached value to 20 + p_id_start = self.pid + 3 + system.part.add(pos=offset, id=p_id_start) + # check the cached value is updated regardless of the MPI node + # where particles are inserted + for i in range(3): + pos = [0., 0., 0.] + pos[i] = -offset[i] # guaranteed to hit different MPI domains + p = self.system.part.add(pos=pos) + self.assertEqual(p.id, p_id_start + i + 1) + # removing the particle with highest id should update the cached value + p.remove() + self.assertEqual(system.part.highest_particle_id, p_id_start + 2) + + # update the cached value to 32 + p_id_start += 10 + p0 = system.part.add(pos=offset, id=p_id_start + 0) + p1 = system.part.add(pos=offset, id=p_id_start + 1) + p2 = system.part.add(pos=offset, id=p_id_start + 2) + # invalidate the cache by introducing a gap at 31 and then removing 32 + p1.remove() + p2.remove() + # the cache should now be 30 + self.assertEqual(system.part.highest_particle_id, p0.id) + p3 = system.part.add(pos=offset) + self.assertEqual(p3.id, p0.id + 1) + def test_invalid_particle_ids_exceptions(self): self.system.part.clear() handle_to_non_existing_particle = self.system.part.by_id(42) @@ -472,7 +508,7 @@ def test_to_dict(self): def test_update(self): self.system.part.clear() p = self.system.part.add(pos=0.5 * self.system.box_l) - # cannot change id + # cannot change id (to avoid corrupting caches in the core) with self.assertRaisesRegex(Exception, "Cannot change particle id."): p.update({'id': 1}) # check value change From 40547141df689b39fdd207f8df9a67261c611bcd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Thu, 7 Apr 2022 21:56:26 +0200 Subject: [PATCH 7/7] python: Fix broken particle creation logic The np.shape() call was raising a numpy deprecation warning on the virtual sites relative property. The dictionary type check was not correct and would fall through for any non-dict object. --- src/python/espressomd/particle_data.pyx | 8 ++++---- testsuite/python/particle.py | 18 +++++++++++++++++- 2 files changed, 21 insertions(+), 5 deletions(-) diff --git a/src/python/espressomd/particle_data.pyx b/src/python/espressomd/particle_data.pyx index ee75e5edc79..1c65ea3295d 100644 --- a/src/python/espressomd/particle_data.pyx +++ b/src/python/espressomd/particle_data.pyx @@ -1767,9 +1767,9 @@ cdef class ParticleList: """ # Did we get a dictionary - if len(args) == 1: - if hasattr(args[0], "__getitem__"): - particles_dict = args[0] + if len(args) == 1 and isinstance( + args[0], (dict, collections.OrderedDict)): + particles_dict = args[0] else: if len(args) == 0 and len(kwargs) != 0: particles_dict = kwargs @@ -1826,7 +1826,7 @@ Set quat and scalar dipole moment (dipm) instead.") def _place_new_particles(self, p_list_dict): # Check if all entries have the same length n_parts = len(p_list_dict["pos"]) - if not all(np.shape(v) and len(v) == + if not all(np.array(v, dtype=object).shape and len(v) == n_parts for v in p_list_dict.values()): raise ValueError( "When adding several particles at once, all lists of attributes have to have the same size") diff --git a/testsuite/python/particle.py b/testsuite/python/particle.py index dca864992f3..7f9c87eecaa 100644 --- a/testsuite/python/particle.py +++ b/testsuite/python/particle.py @@ -21,6 +21,7 @@ import espressomd import espressomd.interactions import numpy as np +import collections class ParticleProperties(ut.TestCase): @@ -346,6 +347,17 @@ def test_invalid_particle_ids_exceptions(self): with self.assertRaisesRegex(ValueError, f"Invalid particle id: {-i}"): self.system.part.add(pos=[0., 0., 0.], id=-i) + def test_invalid_particle_creation(self): + err_msg = r"add\(\) takes either a dictionary or a bunch of keyword args" + with self.assertRaisesRegex(ValueError, err_msg): + self.system.part.add(1) + with self.assertRaisesRegex(ValueError, err_msg): + self.system.part.add(1, 2) + with self.assertRaisesRegex(ValueError, err_msg): + self.system.part.add(1, id=2) + with self.assertRaisesRegex(ValueError, err_msg): + self.system.part.add([1, 2]) + def test_parallel_property_setters(self): s = self.system s.part.clear() @@ -502,8 +514,12 @@ def test_to_dict(self): pp = str(p) pdict = p.to_dict() p.remove() - self.system.part.add(pdict) + p = self.system.part.add(pdict) + self.assertEqual(str(self.system.part.select()), pp) + p.remove() + p = self.system.part.add(collections.OrderedDict(pdict)) self.assertEqual(str(self.system.part.select()), pp) + p.remove() def test_update(self): self.system.part.clear()