Line data Source code
1 : //
2 : // Class HaloCells
3 : // The guard / ghost cells of BareField.
4 : //
5 :
6 : #include <memory>
7 : #include <vector>
8 :
9 : #include "Utility/IpplException.h"
10 :
11 : #include "Communicate/Communicator.h"
12 :
13 : namespace ippl {
14 : namespace detail {
15 : template <typename T, unsigned Dim, class... ViewArgs>
16 726 : HaloCells<T, Dim, ViewArgs...>::HaloCells() {}
17 :
18 : template <typename T, unsigned Dim, class... ViewArgs>
19 0 : void HaloCells<T, Dim, ViewArgs...>::accumulateHalo(view_type& view, Layout_t* layout) {
20 0 : exchangeBoundaries<lhs_plus_assign>(view, layout, HALO_TO_INTERNAL);
21 0 : }
22 :
23 : template <typename T, unsigned Dim, class... ViewArgs>
24 0 : void HaloCells<T, Dim, ViewArgs...>::accumulateHalo_noghost(view_type& view, Layout_t* layout, int nghost) {
25 0 : exchangeBoundaries<lhs_plus_assign>(view, layout, HALO_TO_INTERNAL_NOGHOST, nghost);
26 0 : }
27 : template <typename T, unsigned Dim, class... ViewArgs>
28 0 : void HaloCells<T, Dim, ViewArgs...>::fillHalo(view_type& view, Layout_t* layout) {
29 0 : exchangeBoundaries<assign>(view, layout, INTERNAL_TO_HALO);
30 0 : }
31 :
32 : template <typename T, unsigned Dim, class... ViewArgs>
33 : template <class Op>
34 0 : void HaloCells<T, Dim, ViewArgs...>::exchangeBoundaries(view_type& view, Layout_t* layout,
35 : SendOrder order, int nghost) {
36 : using neighbor_list = typename Layout_t::neighbor_list;
37 : using range_list = typename Layout_t::neighbor_range_list;
38 :
39 0 : auto& comm = layout->comm;
40 :
41 0 : const neighbor_list& neighbors = layout->getNeighbors();
42 0 : const range_list &sendRanges = layout->getNeighborsSendRange(),
43 0 : &recvRanges = layout->getNeighborsRecvRange();
44 :
45 0 : auto ldom = layout->getLocalNDIndex();
46 0 : for (const auto& axis : ldom) {
47 0 : if ((axis.length() == 1) && (Dim != 1)) {
48 0 : throw std::runtime_error(
49 : "HaloCells: Cannot do neighbour exchange when domain decomposition "
50 : "contains planes!");
51 : }
52 : }
53 :
54 : // needed for the NOGHOST approach - we want to remove the ghost
55 : // cells on the boundaries of the global domain from the halo
56 : // exchange when we set HALO_TO_INTERNAL_NOGHOST
57 0 : const auto domain = layout->getDomain();
58 0 : const auto& ldomains = layout->getHostLocalDomains();
59 :
60 0 : size_t totalRequests = 0;
61 0 : for (const auto& componentNeighbors : neighbors) {
62 0 : totalRequests += componentNeighbors.size();
63 : }
64 :
65 0 : int me=Comm->rank();
66 :
67 : using memory_space = typename view_type::memory_space;
68 : using buffer_type = mpi::Communicator::buffer_type<memory_space>;
69 0 : std::vector<MPI_Request> requests(totalRequests);
70 : // sending loop
71 0 : constexpr size_t cubeCount = detail::countHypercubes(Dim) - 1;
72 0 : size_t requestIndex = 0;
73 0 : for (size_t index = 0; index < cubeCount; index++) {
74 0 : int tag = mpi::tag::HALO + index;
75 0 : const auto& componentNeighbors = neighbors[index];
76 0 : for (size_t i = 0; i < componentNeighbors.size(); i++) {
77 0 : int targetRank = componentNeighbors[i];
78 :
79 : bound_type range;
80 0 : if (order == INTERNAL_TO_HALO) {
81 : /*We store only the sending and receiving ranges
82 : * of INTERNAL_TO_HALO and use the fact that the
83 : * sending range of HALO_TO_INTERNAL is the receiving
84 : * range of INTERNAL_TO_HALO and vice versa
85 : */
86 0 : range = sendRanges[index][i];
87 0 : } else if (order == HALO_TO_INTERNAL_NOGHOST) {
88 0 : range = recvRanges[index][i];
89 :
90 0 : for (size_t j = 0; j < Dim; ++j) {
91 0 : bool isLower = ((range.lo[j] + ldomains[me][j].first()
92 0 : - nghost) == domain[j].min());
93 0 : bool isUpper = ((range.hi[j] - 1 +
94 0 : ldomains[me][j].first() - nghost)
95 0 : == domain[j].max());
96 0 : range.lo[j] += isLower * (nghost);
97 0 : range.hi[j] -= isUpper * (nghost);
98 : }
99 : } else {
100 0 : range = recvRanges[index][i];
101 : }
102 :
103 : size_type nsends;
104 0 : pack(range, view, haloData_m, nsends);
105 :
106 0 : buffer_type buf = comm.template getBuffer<memory_space, T>(nsends);
107 :
108 0 : comm.isend(targetRank, tag, haloData_m, *buf, requests[requestIndex++], nsends);
109 0 : buf->resetWritePos();
110 : }
111 : }
112 :
113 : // receiving loop
114 0 : for (size_t index = 0; index < cubeCount; index++) {
115 0 : int tag = mpi::tag::HALO + Layout_t::getMatchingIndex(index);
116 0 : const auto& componentNeighbors = neighbors[index];
117 0 : for (size_t i = 0; i < componentNeighbors.size(); i++) {
118 0 : int sourceRank = componentNeighbors[i];
119 :
120 : bound_type range;
121 0 : if (order == INTERNAL_TO_HALO) {
122 0 : range = recvRanges[index][i];
123 0 : } else if (order == HALO_TO_INTERNAL_NOGHOST) {
124 0 : range = sendRanges[index][i];
125 :
126 0 : for (size_t j = 0; j < Dim; ++j) {
127 0 : bool isLower = ((range.lo[j] + ldomains[me][j].first()
128 0 : - nghost) == domain[j].min());
129 0 : bool isUpper = ((range.hi[j] - 1 +
130 0 : ldomains[me][j].first() - nghost)
131 0 : == domain[j].max());
132 0 : range.lo[j] += isLower * (nghost);
133 0 : range.hi[j] -= isUpper * (nghost);
134 : }
135 : } else {
136 0 : range = sendRanges[index][i];
137 : }
138 :
139 0 : size_type nrecvs = range.size();
140 :
141 0 : buffer_type buf = comm.template getBuffer<memory_space, T>(nrecvs);
142 :
143 0 : comm.recv(sourceRank, tag, haloData_m, *buf, nrecvs * sizeof(T), nrecvs);
144 0 : buf->resetReadPos();
145 :
146 0 : unpack<Op>(range, view, haloData_m);
147 : }
148 : }
149 :
150 0 : if (totalRequests > 0) {
151 0 : MPI_Waitall(totalRequests, requests.data(), MPI_STATUSES_IGNORE);
152 : }
153 :
154 0 : comm.freeAllBuffers();
155 0 : }
156 :
157 : template <typename T, unsigned Dim, class... ViewArgs>
158 0 : void HaloCells<T, Dim, ViewArgs...>::pack(const bound_type& range, const view_type& view,
159 : databuffer_type& fd, size_type& nsends) {
160 0 : auto subview = makeSubview(view, range);
161 :
162 0 : auto& buffer = fd.buffer;
163 :
164 0 : size_t size = subview.size();
165 0 : nsends = size;
166 0 : if (buffer.size() < size) {
167 0 : int overalloc = Comm->getDefaultOverallocation();
168 0 : Kokkos::realloc(buffer, size * overalloc);
169 : }
170 :
171 : using index_array_type =
172 : typename RangePolicy<Dim, typename view_type::execution_space>::index_array_type;
173 0 : ippl::parallel_for(
174 0 : "HaloCells::pack()", getRangePolicy(subview),
175 0 : KOKKOS_LAMBDA(const index_array_type& args) {
176 0 : int l = 0;
177 :
178 0 : for (unsigned d1 = 0; d1 < Dim; d1++) {
179 0 : int next = args[d1];
180 0 : for (unsigned d2 = 0; d2 < d1; d2++) {
181 0 : next *= subview.extent(d2);
182 : }
183 0 : l += next;
184 : }
185 :
186 0 : buffer(l) = apply(subview, args);
187 : });
188 0 : Kokkos::fence();
189 0 : }
190 :
191 : template <typename T, unsigned Dim, class... ViewArgs>
192 : template <typename Op>
193 0 : void HaloCells<T, Dim, ViewArgs...>::unpack(const bound_type& range, const view_type& view,
194 : databuffer_type& fd) {
195 0 : auto subview = makeSubview(view, range);
196 0 : auto buffer = fd.buffer;
197 :
198 : // 29. November 2020
199 : // https://stackoverflow.com/questions/3735398/operator-as-template-parameter
200 : Op op;
201 :
202 : using index_array_type =
203 : typename RangePolicy<Dim, typename view_type::execution_space>::index_array_type;
204 0 : ippl::parallel_for(
205 0 : "HaloCells::unpack()", getRangePolicy(subview),
206 0 : KOKKOS_LAMBDA(const index_array_type& args) {
207 0 : int l = 0;
208 :
209 0 : for (unsigned d1 = 0; d1 < Dim; d1++) {
210 0 : int next = args[d1];
211 0 : for (unsigned d2 = 0; d2 < d1; d2++) {
212 0 : next *= subview.extent(d2);
213 : }
214 0 : l += next;
215 : }
216 :
217 0 : op(apply(subview, args), buffer(l));
218 : });
219 0 : Kokkos::fence();
220 0 : }
221 :
222 : template <typename T, unsigned Dim, class... ViewArgs>
223 0 : auto HaloCells<T, Dim, ViewArgs...>::makeSubview(const view_type& view,
224 : const bound_type& intersect) {
225 0 : auto makeSub = [&]<size_t... Idx>(const std::index_sequence<Idx...>&) {
226 : return Kokkos::subview(view,
227 0 : Kokkos::make_pair(intersect.lo[Idx], intersect.hi[Idx])...);
228 : };
229 0 : return makeSub(std::make_index_sequence<Dim>{});
230 : }
231 :
232 : template <typename T, unsigned Dim, class... ViewArgs>
233 : template <typename Op>
234 9 : void HaloCells<T, Dim, ViewArgs...>::applyPeriodicSerialDim(view_type& view,
235 : const Layout_t* layout,
236 : const int nghost) {
237 9 : int myRank = layout->comm.rank();
238 9 : const auto& lDomains = layout->getHostLocalDomains();
239 9 : const auto& domain = layout->getDomain();
240 :
241 : using exec_space = typename view_type::execution_space;
242 : using index_type = typename RangePolicy<Dim, exec_space>::index_type;
243 :
244 : Kokkos::Array<index_type, Dim> ext, begin, end;
245 :
246 39 : for (size_t i = 0; i < Dim; ++i) {
247 30 : ext[i] = view.extent(i);
248 30 : begin[i] = 0;
249 : }
250 :
251 : Op op;
252 :
253 39 : for (unsigned d = 0; d < Dim; ++d) {
254 30 : end = ext;
255 30 : end[d] = nghost;
256 :
257 30 : if (lDomains[myRank][d].length() == domain[d].length()) {
258 30 : int N = view.extent(d) - 1;
259 :
260 : using index_array_type =
261 : typename RangePolicy<Dim,
262 : typename view_type::execution_space>::index_array_type;
263 30 : ippl::parallel_for(
264 30 : "applyPeriodicSerialDim", createRangePolicy<Dim, exec_space>(begin, end),
265 3002796 : KOKKOS_LAMBDA(index_array_type & coords) {
266 : // The ghosts are filled starting from the inside
267 : // of the domain proceeding outwards for both lower
268 : // and upper faces. The extra brackets and explicit
269 : // mention
270 :
271 : // nghost + i
272 0 : coords[d] += nghost;
273 0 : auto&& left = apply(view, coords);
274 :
275 : // N - nghost - i
276 0 : coords[d] = N - coords[d];
277 0 : auto&& right = apply(view, coords);
278 :
279 : // nghost - 1 - i
280 0 : coords[d] += 2 * nghost - 1 - N;
281 0 : op(apply(view, coords), right);
282 :
283 : // N - (nghost - 1 - i) = N - (nghost - 1) + i
284 0 : coords[d] = N - coords[d];
285 0 : op(apply(view, coords), left);
286 : });
287 :
288 60 : Kokkos::fence();
289 : }
290 : }
291 9 : }
292 : } // namespace detail
293 : } // namespace ippl
|
