ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
send_mesh.cpp
Go to the documentation of this file.
1/*
2 * Copyright (C) 2010-2024 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 "config/config.hpp"
23
24#if defined(P3M) or defined(DP3M)
25
26#include "fft/fft.hpp"
27#include "p3m/common.hpp"
28#include "p3m/packing.hpp"
29#include "p3m/send_mesh.hpp"
30
31#include <utils/Vector.hpp>
33
34#include <boost/mpi/communicator.hpp>
35#include <boost/mpi/datatype.hpp>
36
37#include <mpi.h>
38
39#include <algorithm>
40#include <cstddef>
41#include <span>
42#include <utility>
43
44template <typename T>
45static void mesh_sendrecv(T const *const sendbuf, int scount, int dest,
46 T *const recvbuf, int rcount, int source,
47 boost::mpi::communicator const &comm, int tag) {
48 auto const type = boost::mpi::get_mpi_datatype<T>(*sendbuf);
49 MPI_Sendrecv(reinterpret_cast<void const *>(sendbuf), scount, type, dest, tag,
50 reinterpret_cast<void *>(recvbuf), rcount, type, source, tag,
51 comm, MPI_STATUS_IGNORE);
52}
53
54/** Add values of a 3d-grid input block (size[3]) to values of 3d-grid
55 * output array with dimension dim[3] at start position start[3].
56 *
57 * \param in Pointer to first element of input block data.
58 * \param out Pointer to first element of output grid.
59 * \param start Start position of block in output grid.
60 * \param size Dimensions of the block
61 * \param dim Dimensions of the output grid.
62 */
63template <typename FloatType>
64static void p3m_add_block(FloatType const *in, FloatType *out,
65 int const start[3], int const size[3],
66 int const dim[3]) {
67 /* fast, mid and slow changing indices */
68 int f, m, s;
69 /* linear index of in grid, linear index of out grid */
70 int li_in = 0, li_out = 0;
71 /* offsets for indices in output grid */
72 int m_out_offset, s_out_offset;
73
74 li_out = start[2] + (dim[2] * (start[1] + (dim[1] * start[0])));
75 m_out_offset = dim[2] - size[2];
76 s_out_offset = (dim[2] * (dim[1] - size[1]));
77
78 for (s = 0; s < size[0]; s++) {
79 for (m = 0; m < size[1]; m++) {
80 for (f = 0; f < size[2]; f++) {
81 out[li_out++] += in[li_in++];
82 }
83 li_out += m_out_offset;
84 }
85 li_out += s_out_offset;
86 }
87}
88
89template <typename FloatType>
90void p3m_send_mesh<FloatType>::resize(boost::mpi::communicator const &comm,
91 P3MLocalMesh const &local_mesh) {
92 int done[3] = {0, 0, 0};
93 /* send grids */
94 for (int i = 0; i < 3; i++) {
95 for (int j = 0; j < 3; j++) {
96 /* left */
97 s_ld[i * 2][j] = 0 + done[j] * local_mesh.margin[j * 2];
98 if (j == i)
99 s_ur[i * 2][j] = local_mesh.margin[j * 2];
100 else
101 s_ur[i * 2][j] =
102 local_mesh.dim[j] - done[j] * local_mesh.margin[(j * 2) + 1];
103 /* right */
104 if (j == i)
105 s_ld[(i * 2) + 1][j] = local_mesh.in_ur[j];
106 else
107 s_ld[(i * 2) + 1][j] = 0 + done[j] * local_mesh.margin[j * 2];
108 s_ur[(i * 2) + 1][j] =
109 local_mesh.dim[j] - done[j] * local_mesh.margin[(j * 2) + 1];
110 }
111 done[i] = 1;
112 }
113 max = 0;
114 for (int i = 0; i < 6; i++) {
115 s_size[i] = 1;
116 for (int j = 0; j < 3; j++) {
117 s_dim[i][j] = s_ur[i][j] - s_ld[i][j];
118 s_size[i] *= s_dim[i][j];
119 }
120 max = std::max(max, s_size[i]);
121 }
122 /* communication */
123 auto const node_neighbors = Utils::Mpi::cart_neighbors<3>(comm);
124
125 int r_margin[6];
126 for (int i = 0; i < 6; i++) {
127 auto const j = (i % 2 == 0) ? i + 1 : i - 1;
128
129 if (node_neighbors[i] != comm.rank()) {
130 mesh_sendrecv(&(local_mesh.margin[i]), 1, node_neighbors[i],
131 &(r_margin[j]), 1, node_neighbors[j], comm, REQ_P3M_INIT);
132 } else {
133 r_margin[j] = local_mesh.margin[i];
134 }
135 }
136 /* recv grids */
137 for (int i = 0; i < 3; i++) {
138 for (int j = 0; j < 3; j++) {
139 if (j == i) {
140 r_ld[i * 2][j] = s_ld[i * 2][j] + local_mesh.margin[2 * j];
141 r_ur[i * 2][j] = s_ur[i * 2][j] + r_margin[2 * j];
142 r_ld[(i * 2) + 1][j] = s_ld[(i * 2) + 1][j] - r_margin[(2 * j) + 1];
143 r_ur[(i * 2) + 1][j] =
144 s_ur[(i * 2) + 1][j] - local_mesh.margin[(2 * j) + 1];
145 } else {
146 r_ld[i * 2][j] = s_ld[i * 2][j];
147 r_ur[i * 2][j] = s_ur[i * 2][j];
148 r_ld[(i * 2) + 1][j] = s_ld[(i * 2) + 1][j];
149 r_ur[(i * 2) + 1][j] = s_ur[(i * 2) + 1][j];
150 }
151 }
152 }
153 for (int i = 0; i < 6; i++) {
154 r_size[i] = 1;
155 for (int j = 0; j < 3; j++) {
156 r_dim[i][j] = r_ur[i][j] - r_ld[i][j];
157 r_size[i] *= r_dim[i][j];
158 }
159 max = std::max(max, r_size[i]);
160 }
161}
162
163template <typename FloatType>
164void p3m_send_mesh<FloatType>::gather_grid(boost::mpi::communicator const &comm,
165 std::span<FloatType *> meshes,
166 Utils::Vector3i const &dim) {
167 auto const node_neighbors = Utils::Mpi::cart_neighbors<3>(comm);
168 send_grid.resize(max * meshes.size());
169 recv_grid.resize(max * meshes.size());
170
171 /* direction loop */
172 for (int s_dir = 0; s_dir < 6; s_dir++) {
173 auto const r_dir = (s_dir % 2 == 0) ? s_dir + 1 : s_dir - 1;
174
175 /* pack send block */
176 if (s_size[s_dir] > 0) {
177 for (std::size_t i = 0; i < meshes.size(); i++) {
178 fft_pack_block(meshes[i], send_grid.data() + i * s_size[s_dir],
179 s_ld[s_dir], s_dim[s_dir], dim.data(), 1);
180 }
181 }
182
183 /* communication */
184 if (node_neighbors[s_dir] != comm.rank()) {
185 auto const send_size = static_cast<int>(meshes.size()) * s_size[s_dir];
186 auto const recv_size = static_cast<int>(meshes.size()) * r_size[r_dir];
187 mesh_sendrecv(send_grid.data(), send_size, node_neighbors[s_dir],
188 recv_grid.data(), recv_size, node_neighbors[r_dir], comm,
189 REQ_P3M_GATHER);
190 } else {
191 std::swap(send_grid, recv_grid);
192 }
193 /* add recv block */
194 if (r_size[r_dir] > 0) {
195 for (std::size_t i = 0; i < meshes.size(); i++) {
196 p3m_add_block(recv_grid.data() + i * r_size[r_dir], meshes[i],
197 r_ld[r_dir], r_dim[r_dir], dim.data());
198 }
199 }
200 }
201}
202
203template <typename FloatType>
204void p3m_send_mesh<FloatType>::spread_grid(boost::mpi::communicator const &comm,
205 std::span<FloatType *> meshes,
206 Utils::Vector3i const &dim) {
207 auto const node_neighbors = Utils::Mpi::cart_neighbors<3>(comm);
208 send_grid.resize(max * meshes.size());
209 recv_grid.resize(max * meshes.size());
210
211 /* direction loop */
212 for (int s_dir = 5; s_dir >= 0; s_dir--) {
213 auto const r_dir = (s_dir % 2 == 0) ? s_dir + 1 : s_dir - 1;
214
215 /* pack send block */
216 if (r_size[r_dir] > 0) {
217 for (std::size_t i = 0; i < meshes.size(); i++) {
218 fft_pack_block(meshes[i], send_grid.data() + i * r_size[r_dir],
219 r_ld[r_dir], r_dim[r_dir], dim.data(), 1);
220 }
221 }
222 /* communication */
223 if (node_neighbors[r_dir] != comm.rank()) {
224 auto const send_size = static_cast<int>(meshes.size()) * r_size[r_dir];
225 auto const recv_size = static_cast<int>(meshes.size()) * s_size[s_dir];
226 mesh_sendrecv(send_grid.data(), send_size, node_neighbors[r_dir],
227 recv_grid.data(), recv_size, node_neighbors[s_dir], comm,
228 REQ_P3M_SPREAD);
229 } else {
230 std::swap(send_grid, recv_grid);
231 }
232 /* unpack recv block */
233 if (s_size[s_dir] > 0) {
234 for (std::size_t i = 0; i < meshes.size(); i++) {
235 fft_unpack_block(recv_grid.data() + i * s_size[s_dir], meshes[i],
236 s_ld[s_dir], s_dim[s_dir], dim.data(), 1);
237 }
238 }
239 }
240}
241
242template class p3m_send_mesh<float>;
243template class p3m_send_mesh<double>;
244
245#endif // defined(P3M) or defined(DP3M)
Vector implementation and trait types for boost qvm interoperability.
P3M halo communicator.
Definition send_mesh.hpp:38
void gather_grid(boost::mpi::communicator const &comm, std::span< FloatType * > meshes, Utils::Vector3i const &dim)
void spread_grid(boost::mpi::communicator const &comm, std::span< FloatType * > meshes, Utils::Vector3i const &dim)
void resize(boost::mpi::communicator const &comm, P3MLocalMesh const &local_mesh)
Definition send_mesh.cpp:90
This file contains the defaults for ESPResSo.
Routines, row decomposition, data structures and communication for the 3D-FFT.
Common functions for dipolar and charge P3M.
void fft_pack_block(FloatType const *const in, FloatType *const out, int const *start, int const *size, int const *dim, int element)
Pack a 3D-block of size size starting at start of an input 3D-grid in with dimension dim into an outp...
Definition packing.hpp:44
void fft_unpack_block(FloatType const *const in, FloatType *const out, int const *start, int const *size, int const *dim, int element)
Unpack a 3D-block in of size size into an output 3D-grid out of size dim starting at position start.
Definition packing.hpp:83
static void p3m_add_block(FloatType const *in, FloatType *out, int const start[3], int const size[3], int const dim[3])
Add values of a 3d-grid input block (size[3]) to values of 3d-grid output array with dimension dim[3]...
Definition send_mesh.cpp:64
static void mesh_sendrecv(T const *const sendbuf, int scount, int dest, T *const recvbuf, int rcount, int source, boost::mpi::communicator const &comm, int tag)
Definition send_mesh.cpp:45
Properties of the local mesh.
Utils::Vector3i dim
dimension (size) of local mesh.
int in_ur[3]
inner up right grid point + (1,1,1)
int margin[6]
number of margin mesh points.
DEVICE_QUALIFIER constexpr pointer data() noexcept
Definition Array.hpp:120