ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
particle_node.cpp
Go to the documentation of this file.
1/*
2 * Copyright (C) 2010-2022 The ESPResSo project
3 * Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010
4 * Max-Planck-Institute for Polymer Research, Theory Group
5 *
6 * This file is part of ESPResSo.
7 *
8 * ESPResSo is free software: you can redistribute it and/or modify
9 * it under the terms of the GNU General Public License as published by
10 * the Free Software Foundation, either version 3 of the License, or
11 * (at your option) any later version.
12 *
13 * ESPResSo is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 * GNU General Public License for more details.
17 *
18 * You should have received a copy of the GNU General Public License
19 * along with this program. If not, see <http://www.gnu.org/licenses/>.
20 */
21
22#include "particle_node.hpp"
23
24#include "BoxGeometry.hpp"
25#include "Particle.hpp"
27#include "cells.hpp"
28#include "communication.hpp"
30#include "system/System.hpp"
31
32#include <utils/Cache.hpp>
33#include <utils/Vector.hpp>
34#include <utils/keys.hpp>
35#include <utils/mpi/gatherv.hpp>
36
37#include <boost/mpi/collectives/all_gather.hpp>
38#include <boost/mpi/collectives/all_reduce.hpp>
39#include <boost/mpi/collectives/gather.hpp>
40#include <boost/mpi/collectives/reduce.hpp>
41#include <boost/mpi/collectives/scatter.hpp>
42
43#include <algorithm>
44#include <cmath>
45#include <functional>
46#include <iterator>
47#include <span>
48#include <stdexcept>
49#include <string>
50#include <unordered_map>
51#include <unordered_set>
52#include <utility>
53#include <vector>
54
55constexpr auto some_tag = 42;
56
57/** @brief Enable particle type tracking in @ref particle_type_map. */
58static bool type_list_enable;
59
60/** @brief Mapping particle types to lists of particle ids. */
61static std::unordered_map<int, std::unordered_set<int>> particle_type_map;
62
63/** @brief Mapping particle ids to MPI ranks. */
64static std::unordered_map<int, int> particle_node;
65
66static auto &get_cell_structure() {
68}
69
70/**
71 * @brief Keep track of the largest particle id.
72 * This book-keeping variable is necessary to make particle insertion run
73 * in constant time. Traversing the @ref particle_node to find the largest
74 * particle id scales with O(N) and traversing the local cells in parallel
75 * followed by a reduction scales with O(N^2).
76 */
77static int max_seen_pid = -1;
78
79static auto rebuild_needed() {
80 auto is_rebuild_needed = ::particle_node.empty();
81 boost::mpi::broadcast(::comm_cart, is_rebuild_needed, 0);
82 return is_rebuild_needed;
83}
84
86 boost::mpi::broadcast(::comm_cart, ::max_seen_pid, 0);
87}
88
89void init_type_map(int type) {
90 if (type < 0) {
91 throw std::runtime_error("Types may not be negative");
92 }
93 ::type_list_enable = true;
94 auto &nonbonded_ias = *System::get_system().nonbonded_ias;
95 nonbonded_ias.make_particle_type_exist(type);
96
97 std::vector<int> local_pids;
98 for (auto const &p : get_cell_structure().local_particles()) {
99 if (p.type() == type) {
100 local_pids.emplace_back(p.id());
101 }
102 }
103
104 std::vector<std::vector<int>> global_pids;
105 boost::mpi::all_gather(::comm_cart, local_pids, global_pids);
106 ::particle_type_map[type].clear();
107 for (auto const &vec : global_pids) {
108 for (auto const &p_id : vec) {
109 ::particle_type_map[type].insert(p_id);
110 }
111 }
112}
113
114static void remove_id_from_map(int p_id, int type) {
115 auto it = particle_type_map.find(type);
116 if (it != particle_type_map.end())
117 it->second.erase(p_id);
118}
119
120static void add_id_to_type_map(int p_id, int type) {
121 auto it = particle_type_map.find(type);
122 if (it != particle_type_map.end())
123 it->second.insert(p_id);
124}
125
126void on_particle_type_change(int p_id, int old_type, int new_type) {
127 if (::type_list_enable) {
128 if (old_type == type_tracking::any_type) {
129 for (auto &kv : ::particle_type_map) {
130 auto it = kv.second.find(p_id);
131 if (it != kv.second.end()) {
132 kv.second.erase(it);
133#ifndef NDEBUG
134 if (auto p = get_cell_structure().get_local_particle(p_id)) {
135 assert(p->type() == kv.first);
136 }
137#endif
138 break;
139 }
140 }
141 } else if (old_type != type_tracking::new_part) {
142 if (old_type != new_type) {
143 remove_id_from_map(p_id, old_type);
144 }
145 }
146 add_id_to_type_map(p_id, new_type);
147 }
148}
149
150namespace {
151/* Limit cache to 100 MiB */
152std::size_t const max_cache_size = (100ul * 1048576ul) / sizeof(Particle);
154} // namespace
155
157std::size_t fetch_cache_max_size() { return particle_fetch_cache.max_size(); }
158
159static void mpi_send_particle_data_local(int p_id) {
160 auto const p = get_cell_structure().get_local_particle(p_id);
161 auto const found = p and not p->is_ghost();
162 assert(1 == boost::mpi::all_reduce(::comm_cart, static_cast<int>(found),
163 std::plus<>()) &&
164 "Particle not found");
165 if (found) {
166 ::comm_cart.send(0, 42, *p);
167 }
168}
169
171
172const Particle &get_particle_data(int p_id) {
173 auto const pnode = get_particle_node(p_id);
174
175 if (pnode == this_node) {
176 auto const p = get_cell_structure().get_local_particle(p_id);
177 assert(p != nullptr);
178 return *p;
179 }
180
181 /* Query the cache */
182 auto const p_ptr = particle_fetch_cache.get(p_id);
183 if (p_ptr) {
184 return *p_ptr;
185 }
186
187 /* Cache miss, fetch the particle,
188 * put it into the cache and return a pointer into the cache. */
190 Particle result{};
191 ::comm_cart.recv(boost::mpi::any_source, boost::mpi::any_tag, result);
192 return *(particle_fetch_cache.put(p_id, std::move(result)));
193}
194
196 std::vector<int> ids;
197 boost::mpi::scatter(comm_cart, ids, 0);
198
199 std::vector<Particle> parts(ids.size());
200 std::transform(ids.begin(), ids.end(), parts.begin(), [](int p_id) {
201 auto const p = get_cell_structure().get_local_particle(p_id);
202 assert(p != nullptr);
203 return *p;
204 });
205
206 Utils::Mpi::gatherv(comm_cart, parts.data(), static_cast<int>(parts.size()),
207 0);
208}
209
211
212/**
213 * @brief Get multiple particles at once.
214 *
215 * *WARNING* Particles are returned in an arbitrary order.
216 *
217 * @param ids The ids of the particles that should be returned.
218 *
219 * @returns The particle list.
220 */
221static std::vector<Particle> mpi_get_particles(std::span<const int> ids) {
223 /* Return value */
224 std::vector<Particle> parts(ids.size());
225
226 /* Group ids per node */
227 static std::vector<std::vector<int>> node_ids(comm_cart.size());
228 for (auto &per_node : node_ids) {
229 per_node.clear();
230 }
231
232 for (auto const &p_id : ids) {
233 auto const p_node = get_particle_node(p_id);
234
235 if (p_node != this_node)
236 node_ids[p_node].push_back(p_id);
237 }
238
239 /* Distributed ids to the nodes */
240 {
241 static std::vector<int> ignore;
242 boost::mpi::scatter(comm_cart, node_ids, ignore, 0);
243 }
244
245 /* Copy local particles */
246 std::transform(node_ids[this_node].cbegin(), node_ids[this_node].cend(),
247 parts.begin(), [](int p_id) {
248 auto const p = get_cell_structure().get_local_particle(p_id);
249 assert(p != nullptr);
250 return *p;
251 });
252
253 static std::vector<int> node_sizes(comm_cart.size());
254 std::transform(
255 node_ids.cbegin(), node_ids.cend(), node_sizes.begin(),
256 [](std::vector<int> const &ids) { return static_cast<int>(ids.size()); });
257
258 Utils::Mpi::gatherv(comm_cart, parts.data(), static_cast<int>(parts.size()),
259 parts.data(), node_sizes.data(), 0);
260
261 return parts;
262}
263
264void prefetch_particle_data(std::span<const int> in_ids) {
265 /* Nothing to do on a single node. */
266 // NOLINTNEXTLINE(clang-analyzer-core.NonNullParamChecker)
267 if (comm_cart.size() == 1)
268 return;
269
270 static std::vector<int> ids;
271 ids.clear();
272 auto out_ids = std::back_inserter(ids);
273
274 std::copy_if(in_ids.begin(), in_ids.end(), out_ids, [](int id) {
275 return (get_particle_node(id) != this_node) && particle_fetch_cache.has(id);
276 });
277
278 /* Don't prefetch more particles than fit the cache. */
279 if (ids.size() > particle_fetch_cache.max_size())
280 ids.resize(particle_fetch_cache.max_size());
281
282 /* Fetch the particles... */
283 for (auto &p : mpi_get_particles(ids)) {
284 auto id = p.id();
285 particle_fetch_cache.put(id, std::move(p));
286 }
287}
288
289static void mpi_who_has_local() {
290 static std::vector<int> sendbuf;
291
292 auto local_particles = get_cell_structure().local_particles();
293 auto const n_part = static_cast<int>(local_particles.size());
294 boost::mpi::gather(comm_cart, n_part, 0);
295
296 if (n_part == 0) {
298 return;
299 }
300
301 sendbuf.resize(n_part);
302
303 std::transform(local_particles.begin(), local_particles.end(),
304 sendbuf.begin(), [](Particle const &p) { return p.id(); });
305
306 comm_cart.send(0, some_tag, sendbuf);
308}
309
311
312static void mpi_who_has_head() {
313 auto local_particles = get_cell_structure().local_particles();
314
315 static std::vector<int> n_parts;
316 boost::mpi::gather(comm_cart, static_cast<int>(local_particles.size()),
317 n_parts, 0);
318
319 static std::vector<int> pdata;
320 auto const n_nodes = ::comm_cart.size();
321 max_seen_pid = -1;
322
323 /* then fetch particle locations */
324 for (int pnode = 0; pnode < n_nodes; pnode++) {
325 if (pnode == this_node) {
326 for (auto const &p : local_particles) {
327 particle_node[p.id()] = this_node;
328 max_seen_pid = std::max(max_seen_pid, p.id());
329 }
330 } else if (n_parts[pnode] > 0) {
331 pdata.resize(n_parts[pnode]);
332 comm_cart.recv(pnode, some_tag, pdata);
333 for (int i = 0; i < n_parts[pnode]; i++) {
334 particle_node[pdata[i]] = pnode;
335 max_seen_pid = std::max(max_seen_pid, pdata[i]);
336 }
337 }
338 }
340}
341
342/**
343 * @brief Rebuild the particle index.
344 */
349
350/**
351 * @brief Rebuild the particle index.
352 */
354 if (this_node == 0) {
356 } else {
358 }
359}
360
361int get_particle_node(int p_id) {
362 if (p_id < 0) {
363 throw std::domain_error("Invalid particle id: " + std::to_string(p_id));
364 }
365
366 if (particle_node.empty())
368
369 auto const needle = particle_node.find(p_id);
370
371 // Check if particle has a node, if not, we assume it does not exist.
372 if (needle == particle_node.end()) {
373 throw std::runtime_error("Particle node for id " + std::to_string(p_id) +
374 " not found!");
375 }
376 return needle->second;
377}
378
380 if (p_id < 0) {
381 throw std::domain_error("Invalid particle id: " + std::to_string(p_id));
382 }
383
384 if (rebuild_needed()) {
386 }
387
388 if (this_node != 0) {
389 return -1;
390 }
391
392 auto const needle = particle_node.find(p_id);
393
394 // Check if particle has a node, if not, we assume it does not exist.
395 if (needle == particle_node.end()) {
396 throw std::runtime_error("Particle node for id " + std::to_string(p_id) +
397 " not found!");
398 }
399 return needle->second;
400}
401
403 ::max_seen_pid = -1;
404 particle_node.clear();
405}
406
408 for (auto &kv : ::particle_type_map) {
409 kv.second.clear();
410 }
411}
412
413/**
414 * @brief Calculate the largest particle id.
415 * Traversing the @ref particle_node to find the largest particle id
416 * scales with O(N). Consider using the cached value in @ref max_seen_pid
417 * if possible. This function is only necessary when the cached value is
418 * invalidated, for example when removing the particle which has the
419 * largest id.
420 */
422 return std::accumulate(
423 particle_node.begin(), particle_node.end(), -1,
424 [](int max, auto const &kv) { return std::max(max, kv.first); });
425}
426
427/**
428 * @brief Create a new particle and attach it to a cell.
429 * @param p_id The identity of the particle to create.
430 * @param pos The particle position.
431 * @return Whether the particle was created on that node.
432 */
433static bool maybe_insert_particle(int p_id, Utils::Vector3d const &pos) {
434 auto const &box_geo = *System::get_system().box_geo;
435 auto folded_pos = pos;
436 auto image_box = Utils::Vector3i{};
437 box_geo.fold_position(folded_pos, image_box);
438
439 Particle new_part;
440 new_part.id() = p_id;
441 new_part.pos() = folded_pos;
442 new_part.image_box() = image_box;
443
444 return get_cell_structure().add_local_particle(std::move(new_part)) !=
445 nullptr;
446}
447
448/**
449 * @brief Move particle to a new position.
450 * @param p_id The identity of the particle to move.
451 * @param pos The new particle position.
452 * @return Whether the particle was moved from that node.
453 */
454static bool maybe_move_particle(int p_id, Utils::Vector3d const &pos) {
455 auto const &system = System::get_system();
456 auto const &box_geo = *system.box_geo;
457 auto p = system.cell_structure->get_local_particle(p_id);
458 if (p == nullptr) {
459 return false;
460 }
461 auto folded_pos = pos;
462 auto image_box = Utils::Vector3i{};
463 box_geo.fold_position(folded_pos, image_box);
464 p->pos() = folded_pos;
465 p->image_box() = image_box;
466 return true;
467}
468
475
476void remove_particle(int p_id) {
477 if (::type_list_enable) {
478 auto p = get_cell_structure().get_local_particle(p_id);
479 auto p_type = -1;
480 if (p != nullptr and not p->is_ghost()) {
481 if (this_node == 0) {
482 p_type = p->type();
483 } else {
484 ::comm_cart.send(0, 42, p->type());
485 }
486 } else if (this_node == 0) {
487 ::comm_cart.recv(boost::mpi::any_source, 42, p_type);
488 }
489 assert(this_node != 0 or p_type != -1);
490 boost::mpi::broadcast(::comm_cart, p_type, 0);
491 remove_id_from_map(p_id, p_type);
492 }
493
494 if (this_node == 0) {
495 particle_node[p_id] = -1;
496 }
497 get_cell_structure().remove_particle(p_id);
500 if (this_node == 0) {
501 particle_node.erase(p_id);
502 if (p_id == ::max_seen_pid) {
504 // if there is a gap (i.e. there is no particle with id max_seen_pid - 1,
505 // then the cached value is invalidated and has to be recomputed (slow)
506 if (particle_node.count(::max_seen_pid) == 0 or
509 }
510 }
511 }
513}
514
515void make_new_particle(int p_id, Utils::Vector3d const &pos) {
516 if (rebuild_needed()) {
518 }
519 auto const has_created = maybe_insert_particle(p_id, pos);
521
522 auto node = -1;
523 auto const node_local = (has_created) ? ::comm_cart.rank() : 0;
524 boost::mpi::reduce(::comm_cart, node_local, node, std::plus<int>{}, 0);
525 if (::this_node == 0) {
526 particle_node[p_id] = node;
527 max_seen_pid = std::max(max_seen_pid, p_id);
528 assert(not has_created or node == 0);
529 }
531}
532
533void set_particle_pos(int p_id, Utils::Vector3d const &pos) {
534 auto const has_moved = maybe_move_particle(p_id, pos);
535 get_cell_structure().set_resort_particles(Cells::RESORT_GLOBAL);
537
538 auto success = false;
539 boost::mpi::reduce(::comm_cart, has_moved, success, std::plus<bool>{}, 0);
540 if (::this_node == 0 and !success) {
541 throw std::runtime_error("Particle node for id " + std::to_string(p_id) +
542 " not found!");
543 }
544}
545
546int get_random_p_id(int type, int random_index_in_type_map) {
547 auto it = particle_type_map.find(type);
548 if (it == particle_type_map.end()) {
549 throw std::runtime_error("The provided particle type " +
550 std::to_string(type) +
551 " is currently not tracked by the system.");
552 }
553
554 if (random_index_in_type_map + 1 > it->second.size())
555 throw std::runtime_error("The provided index exceeds the number of "
556 "particle types listed in the particle_type_map");
557 // there is no guarantee of order across MPI ranks
558 auto p_id = *std::next(it->second.begin(), random_index_in_type_map);
559 boost::mpi::broadcast(::comm_cart, p_id, 0);
560 return p_id;
561}
562
564 auto it = particle_type_map.find(type);
565 if (it == particle_type_map.end()) {
566 throw std::runtime_error("The provided particle type " +
567 std::to_string(type) +
568 " is currently not tracked by the system.");
569 }
570
571 return static_cast<int>(it->second.size());
572}
573
574bool particle_exists(int p_id) {
575 if (particle_node.empty())
577 return particle_node.contains(p_id);
578}
579
580std::vector<int> get_particle_ids() {
581 if (particle_node.empty())
583
584 auto ids = Utils::keys(particle_node);
585 std::ranges::sort(ids);
586
587 return ids;
588}
589
590std::vector<int> get_particle_ids_parallel() {
591 if (rebuild_needed()) {
593 }
594 auto pids = Utils::keys(particle_node);
595 boost::mpi::broadcast(::comm_cart, pids, 0);
596 return pids;
597}
598
600 if (rebuild_needed()) {
602 }
603
604 return max_seen_pid;
605}
606
608 if (particle_node.empty())
610
611 return static_cast<int>(particle_node.size());
612}
#define REGISTER_CALLBACK(cb)
Register a static callback without return value.
Vector implementation and trait types for boost qvm interoperability.
This file contains everything related to the global cell structure / cell system.
auto call_all(void(*fp)(Args...), ArgRef &&...args) const -> std::enable_if_t< std::is_void_v< decltype(fp(args...))> >
Call a callback on all nodes.
void on_particle_change()
Called every time a particle property changes.
std::shared_ptr< CellStructure > cell_structure
std::shared_ptr< BoxGeometry > box_geo
std::shared_ptr< InteractionsNonBonded > nonbonded_ias
boost::mpi::communicator comm_cart
The communicator.
int this_node
The number of this node.
void mpi_call(void(*fp)(Args...), ArgRef &&...args)
Call a local function.
MpiCallbacks & mpiCallbacks()
Returns a reference to the global callback class instance.
System & get_system()
void gatherv(const boost::mpi::communicator &comm, const T *in_values, int in_size, T *out_values, const int *sizes, const int *displs, int root)
Definition gatherv.hpp:87
auto keys(Map const &m) -> std::vector< typename Map::key_type >
Return the keys of a map type.
Definition keys.hpp:34
Utils::Cache< int, Particle > particle_fetch_cache(max_cache_size)
auto constexpr new_part
auto constexpr any_type
Various procedures concerning interactions between particles.
std::vector< int > get_particle_ids_parallel()
static int max_seen_pid
Keep track of the largest particle id.
static void add_id_to_type_map(int p_id, int type)
static bool maybe_move_particle(int p_id, Utils::Vector3d const &pos)
Move particle to a new position.
static void build_particle_node()
Rebuild the particle index.
void make_new_particle(int p_id, Utils::Vector3d const &pos)
Create a new particle and attach it to a cell.
const Particle & get_particle_data(int p_id)
Get particle data.
static std::unordered_map< int, int > particle_node
Mapping particle ids to MPI ranks.
int get_particle_node(int p_id)
Get the MPI rank which owns the a specific particle.
int number_of_particles_with_type(int type)
void set_particle_pos(int p_id, Utils::Vector3d const &pos)
Move particle to a new position.
void init_type_map(int type)
static void remove_id_from_map(int p_id, int type)
static void mpi_synchronize_max_seen_pid_local()
void remove_all_particles()
Remove all particles.
int get_random_p_id(int type, int random_index_in_type_map)
Find a particle of given type and return its id.
static int calculate_max_seen_id()
Calculate the largest particle id.
void remove_particle(int p_id)
Remove particle with a given identity.
static void mpi_get_particles_local()
static void mpi_who_has_head()
static std::vector< Particle > mpi_get_particles(std::span< const int > ids)
Get multiple particles at once.
int get_maximal_particle_id()
Get maximal particle id.
int get_particle_node_parallel(int p_id)
constexpr auto some_tag
std::vector< int > get_particle_ids()
Get all particle ids.
std::size_t fetch_cache_max_size()
Return the maximal number of particles that are kept in the fetch cache.
void prefetch_particle_data(std::span< const int > in_ids)
Fetch a range of particle into the fetch cache.
void invalidate_fetch_cache()
Invalidate the fetch cache for get_particle_data.
static auto rebuild_needed()
static auto & get_cell_structure()
static void build_particle_node_parallel()
Rebuild the particle index.
static bool maybe_insert_particle(int p_id, Utils::Vector3d const &pos)
Create a new particle and attach it to a cell.
void clear_particle_node()
Invalidate particle_node.
static bool type_list_enable
Enable particle type tracking in particle_type_map.
static void mpi_who_has_local()
void on_particle_type_change(int p_id, int old_type, int new_type)
static void mpi_send_particle_data_local(int p_id)
int get_n_part()
Get number of particles.
static void clear_particle_type_map()
static std::unordered_map< int, std::unordered_set< int > > particle_type_map
Mapping particle types to lists of particle ids.
bool particle_exists(int p_id)
Check if particle exists.
Particles creation and deletion.
Struct holding all information for one particle.
Definition Particle.hpp:395