Branch data 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 : 756 : 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 : 36 : void HaloCells<T, Dim, ViewArgs...>::applyPeriodicSerialDim(view_type& view,
235 : : const Layout_t* layout,
236 : : const int nghost) {
237 : 36 : int myRank = layout->comm.rank();
238 [ + - ]: 36 : const auto& lDomains = layout->getHostLocalDomains();
239 : 36 : 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 [ + + ]: 162 : for (size_t i = 0; i < Dim; ++i) {
247 : 126 : ext[i] = view.extent(i);
248 : 126 : begin[i] = 0;
249 : : }
250 : :
251 : : Op op;
252 : :
253 [ + + ]: 162 : for (unsigned d = 0; d < Dim; ++d) {
254 : 126 : end = ext;
255 [ + - ]: 126 : end[d] = nghost;
256 : :
257 [ + - ]: 126 : if (lDomains[myRank][d].length() == domain[d].length()) {
258 : 126 : 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 [ + - + - ]: 126 : ippl::parallel_for(
264 : 126 : "applyPeriodicSerialDim", createRangePolicy<Dim, exec_space>(begin, end),
265 [ + - + - : 7234056 : 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 : 3616902 : coords[d] += nghost;
273 : 3616902 : auto&& left = apply(view, coords);
274 : :
275 : : // N - nghost - i
276 : 3616902 : coords[d] = N - coords[d];
277 : 3616902 : auto&& right = apply(view, coords);
278 : :
279 : : // nghost - 1 - i
280 : 3616902 : coords[d] += 2 * nghost - 1 - N;
281 : 3616902 : op(apply(view, coords), right);
282 : :
283 : : // N - (nghost - 1 - i) = N - (nghost - 1) + i
284 : 3616902 : coords[d] = N - coords[d];
285 : 3616902 : op(apply(view, coords), left);
286 : : });
287 : :
288 [ + - + - ]: 252 : Kokkos::fence();
289 : : }
290 : : }
291 : 36 : }
292 : : } // namespace detail
293 : : } // namespace ippl
|
