ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
RegularDecomposition.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
23
24#include "cell_system/Cell.hpp"
25
26#include "communication.hpp"
28#include "errorhandling.hpp"
29#include "system/System.hpp"
30
31#include <utils/Vector.hpp>
32#include <utils/index.hpp>
35
36#include <boost/container/flat_set.hpp>
37#include <boost/mpi/collectives.hpp>
38#include <boost/mpi/communicator.hpp>
39#include <boost/range/algorithm/reverse.hpp>
40#include <boost/range/numeric.hpp>
41
42#include <algorithm>
43#include <array>
44#include <cmath>
45#include <cstddef>
46#include <functional>
47#include <iterator>
48#include <utility>
49#include <vector>
50
51int RegularDecomposition::position_to_cell_index(
52 Utils::Vector3d const &pos) const {
53 Utils::Vector3i cpos;
54
55 for (auto i = 0u; i < 3u; i++) {
56 cpos[i] = static_cast<int>(std::floor(pos[i] * inv_cell_size[i])) + 1 -
57 cell_offset[i];
58
59 /* particles outside our box. Still take them if
60 nonperiodic boundary. We also accept the particle if we are at
61 the box boundary, and the particle is within the box. In this case
62 the particle belongs here and could otherwise potentially be dismissed
63 due to rounding errors. */
64 if (cpos[i] < 1) {
65 if ((!m_box.periodic(i) or (pos[i] >= m_box.length()[i])) and
66 m_local_box.boundary()[2u * i])
67 cpos[i] = 1;
68 else
69 return -1;
70 } else if (cpos[i] > cell_grid[i]) {
71 if ((!m_box.periodic(i) or (pos[i] < m_box.length()[i])) and
72 m_local_box.boundary()[2u * i + 1u])
73 cpos[i] = cell_grid[i];
74 else
75 return -1;
76 }
77 }
78
80}
81
82void RegularDecomposition::move_if_local(
83 ParticleList &src, ParticleList &rest,
84 std::vector<ParticleChange> &modified_cells) {
85 for (auto &part : src) {
86 auto target_cell = position_to_cell(part.pos());
87
88 if (target_cell) {
89 target_cell->particles().insert(std::move(part));
90 modified_cells.emplace_back(ModifiedList{target_cell->particles()});
91 } else {
92 rest.insert(std::move(part));
93 }
94 }
95
96 src.clear();
97}
98
99void RegularDecomposition::move_left_or_right(ParticleList &src,
100 ParticleList &left,
101 ParticleList &right,
102 int dir) const {
103 auto const is_open_boundary_left = m_local_box.boundary()[2 * dir] != 0;
104 auto const is_open_boundary_right = m_local_box.boundary()[2 * dir + 1] != 0;
105 auto const can_move_left = m_box.periodic(dir) or not is_open_boundary_left;
106 auto const can_move_right = m_box.periodic(dir) or not is_open_boundary_right;
107 auto const my_left = m_local_box.my_left()[dir];
108 auto const my_right = m_local_box.my_right()[dir];
109 for (auto it = src.begin(); it != src.end();) {
110 auto const pos = it->pos()[dir];
111 if (m_box.get_mi_coord(pos, my_right, dir) >= 0. and can_move_right) {
112 right.insert(std::move(*it));
113 it = src.erase(it);
114 } else if (m_box.get_mi_coord(pos, my_left, dir) < 0. and can_move_left) {
115 left.insert(std::move(*it));
116 it = src.erase(it);
117 } else {
118 ++it;
119 }
120 }
121}
122
123void RegularDecomposition::exchange_neighbors(
124 ParticleList &pl, std::vector<ParticleChange> &modified_cells) {
125 auto const node_neighbors = Utils::Mpi::cart_neighbors<3>(m_comm);
126 static ParticleList send_buf_l, send_buf_r, recv_buf_l, recv_buf_r;
127
128 for (int dir = 0; dir < 3; dir++) {
129 /* Single node direction, no action needed. */
130 if (Utils::Mpi::cart_get<3>(m_comm).dims[dir] == 1) {
131 continue;
132 /* In this (common) case left and right neighbors are
133 the same, and we need only one communication */
134 }
135 if (Utils::Mpi::cart_get<3>(m_comm).dims[dir] == 2) {
136 move_left_or_right(pl, send_buf_l, send_buf_l, dir);
137
138 Utils::Mpi::sendrecv(m_comm, node_neighbors[2 * dir], 0, send_buf_l,
139 node_neighbors[2 * dir], 0, recv_buf_l);
140
141 send_buf_l.clear();
142 } else {
143 using boost::mpi::request;
145
146 move_left_or_right(pl, send_buf_l, send_buf_r, dir);
147
148 auto req_l = isendrecv(m_comm, node_neighbors[2 * dir], 0, send_buf_l,
149 node_neighbors[2 * dir], 0, recv_buf_l);
150 auto req_r = isendrecv(m_comm, node_neighbors[2 * dir + 1], 0, send_buf_r,
151 node_neighbors[2 * dir + 1], 0, recv_buf_r);
152
153 std::array<request, 4> reqs{{req_l[0], req_l[1], req_r[0], req_r[1]}};
154 boost::mpi::wait_all(reqs.begin(), reqs.end());
155
156 send_buf_l.clear();
157 send_buf_r.clear();
158 }
159
160 move_if_local(recv_buf_l, pl, modified_cells);
161 move_if_local(recv_buf_r, pl, modified_cells);
162 }
163}
164
165/**
166 * @brief Fold coordinates to box and reset the old position.
167 */
168static void fold_and_reset(Particle &p, BoxGeometry const &box_geo) {
169 box_geo.fold_position(p.pos(), p.image_box());
170
172}
173
175 std::vector<ParticleChange> &diff) {
176 ParticleList displaced_parts;
177
178 for (auto &c : local_cells()) {
179 for (auto it = c->particles().begin(); it != c->particles().end();) {
180 fold_and_reset(*it, m_box);
181
182 auto target_cell = particle_to_cell(*it);
183
184 /* Particle is in place */
185 if (target_cell == c) {
186 std::advance(it, 1);
187 continue;
188 }
189
190 auto p = std::move(*it);
191 it = c->particles().erase(it);
192 diff.emplace_back(ModifiedList{c->particles()});
193
194 /* Particle is not local */
195 if (target_cell == nullptr) {
196 diff.emplace_back(RemovedParticle{p.id()});
197 displaced_parts.insert(std::move(p));
198 }
199 /* Particle belongs on this node but is in the wrong cell. */
200 else if (target_cell != c) {
201 target_cell->particles().insert(std::move(p));
202 diff.emplace_back(ModifiedList{target_cell->particles()});
203 }
204 }
205 }
206
207 if (global) {
208 auto const grid = Utils::Mpi::cart_get<3>(m_comm).dims;
209 /* Worst case we need grid - 1 rounds per direction.
210 * This correctly implies that if there is only one node,
211 * no action should be taken. */
212 int rounds_left = grid[0] + grid[1] + grid[2] - 3;
213 for (; rounds_left > 0; rounds_left--) {
214 exchange_neighbors(displaced_parts, diff);
215
216 auto left_over = boost::mpi::all_reduce(m_comm, displaced_parts.size(),
217 std::plus<std::size_t>());
218
219 if (left_over == 0) {
220 break;
221 }
222 }
223 } else {
224 exchange_neighbors(displaced_parts, diff);
225 }
226
227 if (not displaced_parts.empty()) {
228 auto sort_cell = local_cells()[0];
229
230 for (auto &part : displaced_parts) {
231 runtimeErrorMsg() << "Particle " << part.id() << " moved more "
232 << "than one local box length in one timestep";
233 sort_cell->particles().insert(std::move(part));
234
235 diff.emplace_back(ModifiedList{sort_cell->particles()});
236 }
237 }
238}
239
240void RegularDecomposition::mark_cells() {
241 m_local_cells.clear();
242 m_ghost_cells.clear();
243
244 int cnt_c = 0;
245 for (int o = 0; o < ghost_cell_grid[2]; o++)
246 for (int n = 0; n < ghost_cell_grid[1]; n++)
247 for (int m = 0; m < ghost_cell_grid[0]; m++) {
248 if ((m > 0 && m < ghost_cell_grid[0] - 1 && n > 0 &&
249 n < ghost_cell_grid[1] - 1 && o > 0 && o < ghost_cell_grid[2] - 1))
250 m_local_cells.push_back(&cells.at(cnt_c++));
251 else
252 m_ghost_cells.push_back(&cells.at(cnt_c++));
253 }
254}
255
256void RegularDecomposition::fill_comm_cell_lists(ParticleList **part_lists,
257 Utils::Vector3i const &lc,
258 Utils::Vector3i const &hc) {
259 for (int o = lc[0]; o <= hc[0]; o++)
260 for (int n = lc[1]; n <= hc[1]; n++)
261 for (int m = lc[2]; m <= hc[2]; m++) {
262 auto const i = Utils::get_linear_index(o, n, m, ghost_cell_grid);
263
264 *part_lists++ = &(cells.at(i).particles());
265 }
266}
267
269 auto dir_max_range = [this](unsigned int i) {
270 return std::min(0.5 * m_box.length()[i], m_local_box.length()[i]);
271 };
272
273 return {dir_max_range(0u), dir_max_range(1u), dir_max_range(2u)};
274}
275
277int RegularDecomposition::calc_processor_min_num_cells() const {
278 /* the minimal number of cells can be lower if there are at least two nodes
279 serving a direction,
280 since this also ensures that the cell size is at most half the box
281 length. However, if there is only one processor for a direction, there
282 have to be at least two cells for this direction. */
283 return boost::accumulate(Utils::Mpi::cart_get<3>(m_comm).dims, 1,
284 [](int n_cells, int grid) {
285 return (grid == 1) ? 2 * n_cells : n_cells;
286 });
287}
288
289void RegularDecomposition::create_cell_grid(double range) {
290 auto const cart_info = Utils::Mpi::cart_get<3>(m_comm);
291
292 int n_local_cells;
293 auto cell_range = Utils::Vector3d::broadcast(range);
294 auto const min_num_cells = calc_processor_min_num_cells();
295
296 if (range <= 0.) {
297 /* this is the non-interacting case */
298 auto const cells_per_dir =
299 static_cast<int>(std::ceil(std::cbrt(min_num_cells)));
300
301 cell_grid = Utils::Vector3i::broadcast(cells_per_dir);
302 n_local_cells = Utils::product(cell_grid);
303 } else {
304 /* Calculate initial cell grid */
305 auto const &local_box_l = m_local_box.length();
306 auto const volume = Utils::product(local_box_l);
307 auto const scale = std::cbrt(RegularDecomposition::max_num_cells / volume);
308
309 for (auto i = 0u; i < 3u; i++) {
310 /* this is at least 1 */
311 cell_grid[i] = static_cast<int>(std::ceil(local_box_l[i] * scale));
312 cell_range[i] = local_box_l[i] / static_cast<double>(cell_grid[i]);
313
314 if (cell_range[i] < range) {
315 /* ok, too many cells for this direction, set to minimum */
316 cell_grid[i] = static_cast<int>(std::floor(local_box_l[i] / range));
317 if (cell_grid[i] < 1) {
319 << "interaction range " << range << " in direction " << i
320 << " is larger than the local box size " << local_box_l[i];
321 cell_grid[i] = 1;
322 }
323 cell_range[i] = local_box_l[i] / static_cast<double>(cell_grid[i]);
324 }
325 }
326
327 /* It may be necessary to asymmetrically assign the scaling to the
328 coordinates, which the above approach will not do.
329 For a symmetric box, it gives a symmetric result. Here we correct that.
330 */
331 for (;;) {
332 n_local_cells = Utils::product(cell_grid);
333
334 /* done */
335 if (n_local_cells <= RegularDecomposition::max_num_cells)
336 break;
337
338 /* find coordinate with the smallest cell range */
339 auto min_ind = 0u;
340 auto min_size = cell_range[0];
341
342 for (auto i = 1u; i < 3u; ++i) {
343 if (cell_grid[i] > 1 and cell_range[i] < min_size) {
344 min_ind = i;
345 min_size = cell_range[i];
346 }
347 }
348
349 cell_grid[min_ind]--;
350 cell_range[min_ind] = m_local_box.length()[min_ind] / cell_grid[min_ind];
351 }
352
353 /* sanity check */
354 if (n_local_cells < min_num_cells) {
355 runtimeErrorMsg() << "number of cells " << n_local_cells
356 << " is smaller than minimum " << min_num_cells
357 << ": either interaction range is too large for "
358 << "the current skin (range=" << range << ", "
359 << "half_local_box_l=[" << local_box_l / 2. << "]) "
360 << "or min_num_cells too large";
361 }
362 }
363
364 if (n_local_cells > RegularDecomposition::max_num_cells) {
365 runtimeErrorMsg() << "no suitable cell grid found";
366 }
367
368 auto const node_pos = cart_info.coords;
369
370 /* now set all dependent variables */
371 int new_cells = 1;
372 for (auto i = 0u; i < 3u; i++) {
373 ghost_cell_grid[i] = cell_grid[i] + 2;
374 new_cells *= ghost_cell_grid[i];
375 cell_size[i] = m_local_box.length()[i] / static_cast<double>(cell_grid[i]);
376 inv_cell_size[i] = 1.0 / cell_size[i];
377 cell_offset[i] = node_pos[i] * cell_grid[i];
378 }
379
380 /* allocate cell array and cell pointer arrays */
381 cells.clear();
382 cells.resize(static_cast<unsigned int>(new_cells));
383 m_local_cells.resize(n_local_cells);
384 m_ghost_cells.resize(new_cells - n_local_cells);
385}
386
387template <class K, class Comparator> auto make_flat_set(Comparator &&comp) {
388 return boost::container::flat_set<K, std::remove_reference_t<Comparator>>(
389 std::forward<Comparator>(comp));
390}
391
392void RegularDecomposition::init_cell_interactions() {
393
394 auto const halo = Utils::Vector3i{1, 1, 1};
395 auto const cart_info = Utils::Mpi::cart_get<3>(m_comm);
396 auto const &node_pos = cart_info.coords;
397 auto const &node_grid = ::communicator.node_grid;
398 auto const global_halo_offset = hadamard_product(node_pos, cell_grid) - halo;
399 auto const global_size = hadamard_product(node_grid, cell_grid);
400 auto const at_boundary = [&global_size](int coord, Utils::Vector3i cell_idx) {
401 return (cell_idx[coord] == 0 or cell_idx[coord] == global_size[coord]);
402 };
403
404 // For the fully connected feature (cells that don't share at least a corner)
405 // only apply if one cell is a ghost cell (i.e. connections across the
406 // periodic boundary.
407 auto const fcb_is_inner_connection = [&global_size, this](Utils::Vector3i a,
408 Utils::Vector3i b) {
410 auto const [fc_normal, fc_dir] = *fully_connected_boundary();
411 auto const involves_ghost_cell =
412 (a[fc_normal] == -1 or a[fc_normal] == global_size[fc_normal] or
413 b[fc_normal] == -1 or b[fc_normal] == global_size[fc_normal]);
414 if (not involves_ghost_cell) {
415 // check if cells do not share at least a corner
416 return std::abs((a - b)[fc_dir]) > 1;
417 }
418 }
419 return false;
420 };
421
422 /* Translate a node local index (relative to the origin of the local grid)
423 * to a global index. */
424 auto global_index =
425 [&](Utils::Vector3i const &local_index) -> Utils::Vector3i {
426 return (global_halo_offset + local_index);
427 };
428
429 /* Linear index in the global cell grid. */
430 auto folded_linear_index = [&](Utils::Vector3i const &global_index) {
431 auto const folded_index = (global_index + global_size) % global_size;
432
433 return get_linear_index(folded_index, global_size);
434 };
435
436 /* Translate a global index into a local one */
437 auto local_index =
438 [&](Utils::Vector3i const &global_index) -> Utils::Vector3i {
439 return (global_index - global_halo_offset);
440 };
441
442 // sanity checks
444 auto const [fc_normal, fc_dir] = *fully_connected_boundary();
445 if (fc_normal == fc_dir) {
446 throw std::domain_error("fully_connected_boundary normal and connection "
447 "coordinates need to differ.");
448 }
449 if (node_grid[fc_dir] != 1) {
450 throw std::runtime_error(
451 "The MPI nodegrid must be 1 in the fully connected direction.");
452 }
453 }
454
455 /* We only consider local cells (e.g. not halo cells), which
456 * span the range [(1,1,1), cell_grid) in local coordinates. */
457 auto const start = global_index(Utils::Vector3i{1, 1, 1});
458 auto const end = start + cell_grid;
459
460 /* loop all local cells */
461 for (int o = start[2]; o < end[2]; o++)
462 for (int n = start[1]; n < end[1]; n++)
463 for (int m = start[0]; m < end[0]; m++) {
464 /* next-nearest neighbors in every direction */
465 Utils::Vector3i lower_index = {m - 1, n - 1, o - 1};
466 Utils::Vector3i upper_index = {m + 1, n + 1, o + 1};
467
468 /* In the fully connected case, we consider all cells
469 * in the direction as neighbors, not only the nearest ones.
470 // */
472 auto const [fc_boundary, fc_direction] = *fully_connected_boundary();
473
474 // Fully connected is only needed at the box surface
475 if (at_boundary(fc_boundary, {m, n, o})) {
476 lower_index[fc_direction] = -1;
477 upper_index[fc_direction] = global_size[fc_boundary];
478 }
479 }
480
481 /* In non-periodic directions, the halo needs not
482 * be considered. */
483 for (auto i = 0u; i < 3u; i++) {
484 if (not m_box.periodic(i)) {
485 lower_index[i] = std::max(0, lower_index[i]);
486 upper_index[i] = std::min(global_size[i] - 1, upper_index[i]);
487 }
488 }
489
490 /* Unique set of neighbors, cells are compared by their linear
491 * index in the global cell grid. */
492 auto neighbors = make_flat_set<Utils::Vector3i>(
493 [&](Utils::Vector3i const &a, Utils::Vector3i const &b) {
494 return folded_linear_index(a) < folded_linear_index(b);
495 });
496
497 /* Collect neighbors */
498 for (int p = lower_index[2]; p <= upper_index[2]; p++)
499 for (int q = lower_index[1]; q <= upper_index[1]; q++)
500 for (int r = lower_index[0]; r <= upper_index[0]; r++) {
502 // Avoid fully connecting the boundary layer and the
503 // next INNER layer
504 if (fcb_is_inner_connection({m, n, o}, {r, q, p}))
505 continue;
506 }
507 neighbors.insert(Utils::Vector3i{r, q, p});
508 }
509
510 /* Red-black partition by global index. */
511 auto const ind1 = folded_linear_index({m, n, o});
512
513 std::vector<Cell *> red_neighbors;
514 std::vector<Cell *> black_neighbors;
515 for (auto const &neighbor : neighbors) {
516 auto const ind2 = folded_linear_index(neighbor);
517 /* Exclude cell itself */
518 if (ind1 == ind2)
519 continue;
520
521 auto cell = &cells.at(
522 get_linear_index(local_index(neighbor), ghost_cell_grid));
523 if (ind2 > ind1) {
524 red_neighbors.push_back(cell);
525 } else {
526 black_neighbors.push_back(cell);
527 }
528 }
529
530 cells[get_linear_index(local_index({m, n, o}), ghost_cell_grid)]
531 .m_neighbors = Neighbors<Cell *>(red_neighbors, black_neighbors);
532 }
533}
534
535namespace {
536/** Revert the order of a communicator: After calling this the
537 * communicator is working in reverted order with exchanged
538 * communication types GHOST_SEND <-> GHOST_RECV.
539 */
541 /* revert order */
542 boost::reverse(comm.communications);
543
544 /* exchange SEND/RECV */
545 for (auto &c : comm.communications) {
546 if (c.type == GHOST_SEND)
547 c.type = GHOST_RECV;
548 else if (c.type == GHOST_RECV)
549 c.type = GHOST_SEND;
550 else if (c.type == GHOST_LOCL) {
551 boost::reverse(c.part_lists);
552 }
553 }
554}
555
556/** Of every two communication rounds, set the first receivers to prefetch and
557 * poststore
558 */
560 for (auto it = comm.communications.begin(); it != comm.communications.end();
561 it += 2) {
562 auto next = std::next(it);
563 if (it->type == GHOST_RECV && next->type == GHOST_SEND) {
564 it->type |= GHOST_PREFETCH | GHOST_PSTSTORE;
565 next->type |= GHOST_PREFETCH | GHOST_PSTSTORE;
566 }
567 }
568}
569
570} // namespace
571
572GhostCommunicator RegularDecomposition::prepare_comm() {
573 int dir, lr, i, cnt, n_comm_cells[3];
574 Utils::Vector3i lc{}, hc{}, done{};
575
576 auto const comm_info = Utils::Mpi::cart_get<3>(m_comm);
577 auto const node_neighbors = Utils::Mpi::cart_neighbors<3>(m_comm);
578
579 /* calculate number of communications */
580 std::size_t num = 0;
581 for (dir = 0; dir < 3; dir++) {
582 for (lr = 0; lr < 2; lr++) {
583 /* No communication for border of non periodic direction */
584 if (comm_info.dims[dir] == 1)
585 num++;
586 else
587 num += 2;
588 }
589 }
590
591 /* prepare communicator */
592 auto ghost_comm = GhostCommunicator{m_comm, num};
593
594 /* number of cells to communicate in a direction */
595 n_comm_cells[0] = cell_grid[1] * cell_grid[2];
596 n_comm_cells[1] = cell_grid[2] * ghost_cell_grid[0];
597 n_comm_cells[2] = ghost_cell_grid[0] * ghost_cell_grid[1];
598
599 cnt = 0;
600 /* direction loop: x, y, z */
601 for (dir = 0; dir < 3; dir++) {
602 lc[(dir + 1) % 3] = 1 - done[(dir + 1) % 3];
603 lc[(dir + 2) % 3] = 1 - done[(dir + 2) % 3];
604 hc[(dir + 1) % 3] = cell_grid[(dir + 1) % 3] + done[(dir + 1) % 3];
605 hc[(dir + 2) % 3] = cell_grid[(dir + 2) % 3] + done[(dir + 2) % 3];
606 /* lr loop: left right */
607 /* here we could in principle build in a one sided ghost
608 communication, simply by taking the lr loop only over one
609 value */
610 for (lr = 0; lr < 2; lr++) {
611 if (comm_info.dims[dir] == 1) {
612 /* just copy cells on a single node */
613 ghost_comm.communications[cnt].type = GHOST_LOCL;
614 ghost_comm.communications[cnt].node = m_comm.rank();
615
616 /* Buffer has to contain Send and Recv cells -> factor 2 */
617 ghost_comm.communications[cnt].part_lists.resize(2 * n_comm_cells[dir]);
618
619 /* fill send ghost_comm cells */
620 lc[dir] = hc[dir] = 1 + lr * (cell_grid[dir] - 1);
621
622 fill_comm_cell_lists(ghost_comm.communications[cnt].part_lists.data(),
623 lc, hc);
624
625 /* fill recv ghost_comm cells */
626 lc[dir] = hc[dir] = 0 + (1 - lr) * (cell_grid[dir] + 1);
627
628 /* place receive cells after send cells */
629 fill_comm_cell_lists(
630 &ghost_comm.communications[cnt].part_lists[n_comm_cells[dir]], lc,
631 hc);
632
633 cnt++;
634 } else {
635 /* i: send/recv loop */
636 for (i = 0; i < 2; i++) {
637 if ((comm_info.coords[dir] + i) % 2 == 0) {
638 ghost_comm.communications[cnt].type = GHOST_SEND;
639 ghost_comm.communications[cnt].node = node_neighbors[2 * dir + lr];
640 ghost_comm.communications[cnt].part_lists.resize(n_comm_cells[dir]);
641 /* prepare folding of ghost positions */
642
643 lc[dir] = hc[dir] = 1 + lr * (cell_grid[dir] - 1);
644
645 fill_comm_cell_lists(
646 ghost_comm.communications[cnt].part_lists.data(), lc, hc);
647 cnt++;
648 }
649 if ((comm_info.coords[dir] + (1 - i)) % 2 == 0) {
650 ghost_comm.communications[cnt].type = GHOST_RECV;
651 ghost_comm.communications[cnt].node =
652 node_neighbors[2 * dir + (1 - lr)];
653 ghost_comm.communications[cnt].part_lists.resize(n_comm_cells[dir]);
654
655 lc[dir] = hc[dir] = (1 - lr) * (cell_grid[dir] + 1);
656
657 fill_comm_cell_lists(
658 ghost_comm.communications[cnt].part_lists.data(), lc, hc);
659 cnt++;
660 }
661 }
662 }
663 done[dir] = 1;
664 }
665 }
666
667 return ghost_comm;
668}
669
671 boost::mpi::communicator comm, double range, BoxGeometry const &box_geo,
672 LocalBox const &local_geo,
673 std::optional<std::pair<int, int>> fully_connected)
674 : m_comm(std::move(comm)), m_box(box_geo), m_local_box(local_geo),
675 m_fully_connected_boundary(std::move(fully_connected)) {
676
677 /* set up new regular decomposition cell structure */
678 create_cell_grid(range);
679
680 /* setup cell neighbors */
681 init_cell_interactions();
682
683 /* mark local and ghost cells */
684 mark_cells();
685
686 /* create communicators */
687 m_exchange_ghosts_comm = prepare_comm();
688 m_collect_ghost_force_comm = prepare_comm();
689
690 /* collect forces has to be done in reverted order! */
691 revert_comm_order(m_collect_ghost_force_comm);
692
693 assign_prefetches(m_exchange_ghosts_comm);
694 assign_prefetches(m_collect_ghost_force_comm);
695}
static int coord(std::string const &s)
auto make_flat_set(Comparator &&comp)
static void fold_and_reset(Particle &p, BoxGeometry const &box_geo)
Fold coordinates to box and reset the old position.
Vector implementation and trait types for boost qvm interoperability.
T get_mi_coord(T a, T b, unsigned coord) const
Get the minimum-image distance between two coordinates.
Utils::Vector3d const & length() const
Box length.
constexpr bool periodic(unsigned coord) const
Check periodicity in direction.
void fold_position(Utils::Vector3d &pos, Utils::Vector3i &image_box) const
Fold coordinates to primary simulation box in-place.
auto const & my_right() const
Right (top, back) corner of this nodes local box.
Definition LocalBox.hpp:47
auto const & boundary() const
Boundary information for the local box.
Definition LocalBox.hpp:59
auto const & my_left() const
Left (bottom, front) corner of this nodes local box.
Definition LocalBox.hpp:45
auto const & length() const
Dimensions of the box a single node is responsible for.
Definition LocalBox.hpp:49
iterator begin()
Definition Bag.hpp:82
void clear()
Remove all elements form container.
Definition Bag.hpp:134
bool empty() const
Is the container empty?
Definition Bag.hpp:96
T & insert(T const &v)
Insert an element into the container.
Definition Bag.hpp:147
iterator erase(iterator it)
Remove element from the list.
Definition Bag.hpp:165
iterator end()
Definition Bag.hpp:83
std::size_t size() const
Number of elements in the container.
Definition Bag.hpp:90
static DEVICE_QUALIFIER constexpr Vector< T, N > broadcast(typename Base::value_type const &value)
Create a vector that has all entries set to the same value.
Definition Vector.hpp:110
Communicator communicator
This file contains the errorhandling code for severe errors, like a broken bond or illegal parameter ...
#define runtimeErrorMsg()
#define GHOST_RECV
recv from a single node
Definition ghosts.hpp:106
#define GHOST_PSTSTORE
additional flag for poststoring
Definition ghosts.hpp:119
#define GHOST_LOCL
transfer data from cell to cell on this node
Definition ghosts.hpp:112
#define GHOST_SEND
send to a single node
Definition ghosts.hpp:104
#define GHOST_PREFETCH
additional flag for prefetching
Definition ghosts.hpp:117
std::array< mpi::request, 2 > isendrecv(mpi::communicator const &comm, int dest, int stag, const T &sval, int src, int rtag, T &rval)
Definition sendrecv.hpp:73
mpi::status sendrecv(mpi::communicator const &comm, int dest, int stag, const T &sval, int src, int rtag, T &rval)
Definition sendrecv.hpp:66
T product(Vector< T, N > const &v)
Definition Vector.hpp:374
int get_linear_index(int a, int b, int c, const Vector3i &adim, MemoryOrder memory_order=MemoryOrder::COLUMN_MAJOR)
get the linear index from the position (a,b,c) in a 3D grid of dimensions adim.
Definition index.hpp:43
auto hadamard_product(Vector< T, N > const &a, Vector< U, N > const &b)
Definition Vector.hpp:379
void assign_prefetches(GhostCommunicator &comm)
Of every two communication rounds, set the first receivers to prefetch and poststore.
void revert_comm_order(GhostCommunicator &comm)
Revert the order of a communicator: After calling this the communicator is working in reverted order ...
Utils::Vector3i node_grid
Properties for a ghost communication.
Definition ghosts.hpp:157
std::vector< GhostCommunication > communications
List of ghost communications.
Definition ghosts.hpp:166
Struct holding all information for one particle.
Definition Particle.hpp:395
auto & pos_at_last_verlet_update()
Definition Particle.hpp:442
auto const & image_box() const
Definition Particle.hpp:444
auto const & pos() const
Definition Particle.hpp:431
GhostCommunicator m_exchange_ghosts_comm
Utils::Vector3i ghost_cell_grid
linked cell grid with ghost frame.
Utils::Vector3d max_cutoff() const override
Utils::Vector3d inv_cell_size
inverse cell_size.
std::vector< Cell * > m_ghost_cells
BoxGeometry const & m_box
Cell * particle_to_cell(Particle const &p) override
void resort(bool global, std::vector< ParticleChange > &diff) override
std::vector< Cell > cells
RegularDecomposition(boost::mpi::communicator comm, double range, BoxGeometry const &box_geo, LocalBox const &local_geo, std::optional< std::pair< int, int > > fully_connected)
Utils::Vector3d cell_size
Cell size.
std::span< Cell *const > local_cells() const override
std::vector< Cell * > m_local_cells
Utils::Vector3d max_range() const override
GhostCommunicator m_collect_ghost_force_comm
Utils::Vector3i cell_grid
Grid dimensions per node.
Utils::Vector3i cell_offset
Offset in global grid.
boost::mpi::communicator m_comm