26#include "communication.hpp"
29#include "system/System.hpp"
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>
51int RegularDecomposition::position_to_cell_index(
55 for (
auto i = 0u; i < 3u; i++) {
56 cpos[i] =
static_cast<int>(std::floor(pos[i] *
inv_cell_size[i])) + 1 -
82void RegularDecomposition::move_if_local(
84 std::vector<ParticleChange> &modified_cells) {
85 for (
auto &part : src) {
86 auto target_cell = position_to_cell(part.pos());
89 target_cell->particles().insert(std::move(part));
90 modified_cells.emplace_back(
ModifiedList{target_cell->particles()});
92 rest.
insert(std::move(part));
99void RegularDecomposition::move_left_or_right(
ParticleList &src,
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;
109 for (
auto it = src.
begin(); it != src.
end();) {
110 auto const pos = it->pos()[dir];
112 right.
insert(std::move(*it));
115 left.
insert(std::move(*it));
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;
128 for (
int dir = 0; dir < 3; dir++) {
130 if (Utils::Mpi::cart_get<3>(
m_comm).dims[dir] == 1) {
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);
139 node_neighbors[2 * dir], 0, recv_buf_l);
143 using boost::mpi::request;
146 move_left_or_right(pl, send_buf_l, send_buf_r, dir);
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);
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());
160 move_if_local(recv_buf_l, pl, modified_cells);
161 move_if_local(recv_buf_r, pl, modified_cells);
175 std::vector<ParticleChange> &diff) {
179 for (
auto it = c->particles().begin(); it != c->particles().end();) {
185 if (target_cell == c) {
190 auto p = std::move(*it);
191 it = c->particles().erase(it);
195 if (target_cell ==
nullptr) {
197 displaced_parts.
insert(std::move(p));
200 else if (target_cell != c) {
201 target_cell->particles().insert(std::move(p));
202 diff.emplace_back(
ModifiedList{target_cell->particles()});
208 auto const grid = Utils::Mpi::cart_get<3>(
m_comm).dims;
212 int rounds_left = grid[0] + grid[1] + grid[2] - 3;
213 for (; rounds_left > 0; rounds_left--) {
214 exchange_neighbors(displaced_parts, diff);
216 auto left_over = boost::mpi::all_reduce(
m_comm, displaced_parts.
size(),
217 std::plus<std::size_t>());
219 if (left_over == 0) {
224 exchange_neighbors(displaced_parts, diff);
227 if (not displaced_parts.
empty()) {
230 for (
auto &part : displaced_parts) {
232 <<
"than one local box length in one timestep";
233 sort_cell->particles().insert(std::move(part));
235 diff.emplace_back(
ModifiedList{sort_cell->particles()});
240void RegularDecomposition::mark_cells() {
256void RegularDecomposition::fill_comm_cell_lists(
ParticleList **part_lists,
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++) {
264 *part_lists++ = &(
cells.at(i).particles());
269 auto dir_max_range = [
this](
unsigned int i) {
273 return {dir_max_range(0u), dir_max_range(1u), dir_max_range(2u)};
277int RegularDecomposition::calc_processor_min_num_cells()
const {
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;
289void RegularDecomposition::create_cell_grid(
double range) {
290 auto const cart_info = Utils::Mpi::cart_get<3>(
m_comm);
294 auto const min_num_cells = calc_processor_min_num_cells();
298 auto const cells_per_dir =
299 static_cast<int>(std::ceil(std::cbrt(min_num_cells)));
307 auto const scale = std::cbrt(RegularDecomposition::max_num_cells / volume);
309 for (
auto i = 0u; i < 3u; i++) {
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]);
314 if (cell_range[i] < range) {
316 cell_grid[i] =
static_cast<int>(std::floor(local_box_l[i] / range));
319 <<
"interaction range " << range <<
" in direction " << i
320 <<
" is larger than the local box size " << local_box_l[i];
323 cell_range[i] = local_box_l[i] /
static_cast<double>(
cell_grid[i]);
335 if (n_local_cells <= RegularDecomposition::max_num_cells)
340 auto min_size = cell_range[0];
342 for (
auto i = 1u; i < 3u; ++i) {
343 if (
cell_grid[i] > 1 and cell_range[i] < min_size) {
345 min_size = cell_range[i];
354 if (n_local_cells < min_num_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";
364 if (n_local_cells > RegularDecomposition::max_num_cells) {
368 auto const node_pos = cart_info.coords;
372 for (
auto i = 0u; i < 3u; i++) {
382 cells.resize(
static_cast<unsigned int>(new_cells));
388 return boost::container::flat_set<K, std::remove_reference_t<Comparator>>(
389 std::forward<Comparator>(comp));
392void RegularDecomposition::init_cell_interactions() {
395 auto const cart_info = Utils::Mpi::cart_get<3>(
m_comm);
396 auto const &node_pos = cart_info.coords;
401 return (cell_idx[
coord] == 0 or cell_idx[
coord] == global_size[
coord]);
407 auto const fcb_is_inner_connection = [&global_size,
this](
Utils::Vector3i a,
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) {
416 return std::abs((a - b)[fc_dir]) > 1;
426 return (global_halo_offset + local_index);
431 auto const folded_index = (global_index + global_size) % global_size;
439 return (global_index - global_halo_offset);
445 if (fc_normal == fc_dir) {
446 throw std::domain_error(
"fully_connected_boundary normal and connection "
447 "coordinates need to differ.");
449 if (node_grid[fc_dir] != 1) {
450 throw std::runtime_error(
451 "The MPI nodegrid must be 1 in the fully connected direction.");
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++) {
475 if (at_boundary(fc_boundary, {m, n, o})) {
476 lower_index[fc_direction] = -1;
477 upper_index[fc_direction] = global_size[fc_boundary];
483 for (
auto i = 0u; i < 3u; i++) {
485 lower_index[i] = std::max(0, lower_index[i]);
486 upper_index[i] = std::min(global_size[i] - 1, upper_index[i]);
492 auto neighbors = make_flat_set<Utils::Vector3i>(
494 return folded_linear_index(a) < folded_linear_index(b);
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++) {
504 if (fcb_is_inner_connection({m, n, o}, {r, q, p}))
511 auto const ind1 = folded_linear_index({m, n, o});
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);
521 auto cell = &
cells.at(
524 red_neighbors.push_back(cell);
526 black_neighbors.push_back(cell);
551 boost::reverse(c.part_lists);
562 auto next = std::next(it);
573 int dir, lr, i, cnt, n_comm_cells[3];
576 auto const comm_info = Utils::Mpi::cart_get<3>(
m_comm);
577 auto const node_neighbors = Utils::Mpi::cart_neighbors<3>(
m_comm);
581 for (dir = 0; dir < 3; dir++) {
582 for (lr = 0; lr < 2; lr++) {
584 if (comm_info.dims[dir] == 1)
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];
610 for (lr = 0; lr < 2; lr++) {
611 if (comm_info.dims[dir] == 1) {
613 ghost_comm.communications[cnt].type =
GHOST_LOCL;
614 ghost_comm.communications[cnt].node =
m_comm.rank();
617 ghost_comm.communications[cnt].part_lists.resize(2 * n_comm_cells[dir]);
620 lc[dir] = hc[dir] = 1 + lr * (
cell_grid[dir] - 1);
622 fill_comm_cell_lists(ghost_comm.communications[cnt].part_lists.data(),
626 lc[dir] = hc[dir] = 0 + (1 - lr) * (
cell_grid[dir] + 1);
629 fill_comm_cell_lists(
630 &ghost_comm.communications[cnt].part_lists[n_comm_cells[dir]], lc,
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]);
643 lc[dir] = hc[dir] = 1 + lr * (
cell_grid[dir] - 1);
645 fill_comm_cell_lists(
646 ghost_comm.communications[cnt].part_lists.data(), lc, hc);
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]);
655 lc[dir] = hc[dir] = (1 - lr) * (
cell_grid[dir] + 1);
657 fill_comm_cell_lists(
658 ghost_comm.communications[cnt].part_lists.data(), lc, hc);
671 boost::mpi::communicator comm,
double range,
BoxGeometry const &box_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)) {
678 create_cell_grid(range);
681 init_cell_interactions();
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.
auto const & boundary() const
Boundary information for the local box.
auto const & my_left() const
Left (bottom, front) corner of this nodes local box.
auto const & length() const
Dimensions of the box a single node is responsible for.
void clear()
Remove all elements form container.
bool empty() const
Is the container empty?
T & insert(T const &v)
Insert an element into the container.
iterator erase(iterator it)
Remove element from the list.
std::size_t size() const
Number of elements in the container.
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.
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
#define GHOST_PSTSTORE
additional flag for poststoring
#define GHOST_LOCL
transfer data from cell to cell on this node
#define GHOST_SEND
send to a single node
#define GHOST_PREFETCH
additional flag for prefetching
std::array< mpi::request, 2 > isendrecv(mpi::communicator const &comm, int dest, int stag, const T &sval, int src, int rtag, T &rval)
mpi::status sendrecv(mpi::communicator const &comm, int dest, int stag, const T &sval, int src, int rtag, T &rval)
T product(Vector< T, N > const &v)
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.
auto hadamard_product(Vector< T, N > const &a, Vector< U, N > const &b)
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.
std::vector< GhostCommunication > communications
List of ghost communications.
Struct holding all information for one particle.
auto & pos_at_last_verlet_update()
auto const & image_box() const
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
auto fully_connected_boundary() const
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