Line data Source code
1 : //
2 : // Class FieldLayout
3 : // FieldLayout describes how a given index space (represented by an NDIndex
4 : // object) is distributed among MPI ranks. It performs the initial
5 : // partitioning. The user may request that a particular dimension not be
6 : // partitioned by flagging that axis as 'SERIAL' (instead of 'PARALLEL').
7 : //
8 : #include "Ippl.h"
9 :
10 : #include <cstdlib>
11 : #include <limits>
12 :
13 : #include "Utility/IpplException.h"
14 : #include "Utility/IpplTimings.h"
15 : #include "Utility/PAssert.h"
16 :
17 : #include "FieldLayout/FieldLayout.h"
18 :
19 : namespace ippl {
20 :
21 : template <unsigned Dim>
22 109800 : int FieldLayout<Dim>::getMatchingIndex(int index) {
23 : // If we are touching another boundary component, that component must be parallel to us,
24 : // so any 2s are unchanged. All other digits must be swapped because otherwise we would
25 : // share a higher dimension boundary region with that other component and the digit
26 : // would be 2
27 : static const int digit_swap[3] = {1, 0, 2};
28 109800 : int match = 0;
29 712576 : for (unsigned d = 1; d < detail::countHypercubes(Dim); d *= 3) {
30 602776 : match += digit_swap[index % 3] * d;
31 602776 : index /= 3;
32 : }
33 109800 : return match;
34 : }
35 :
36 : template <unsigned Dim>
37 2652 : FieldLayout<Dim>::FieldLayout(const mpi::Communicator& communicator)
38 2652 : : comm(communicator)
39 2652 : , dLocalDomains_m("local domains (device)", 0)
40 5304 : , hLocalDomains_m(Kokkos::create_mirror_view(dLocalDomains_m)) {
41 11268 : for (unsigned int d = 0; d < Dim; ++d) {
42 8616 : minWidth_m[d] = 0;
43 : }
44 2652 : }
45 :
46 : template <unsigned Dim>
47 1324 : FieldLayout<Dim>::FieldLayout(mpi::Communicator communicator, const NDIndex<Dim>& domain,
48 : std::array<bool, Dim> isParallel, bool isAllPeriodic)
49 1324 : : FieldLayout(communicator) {
50 1324 : initialize(domain, isParallel, isAllPeriodic);
51 1324 : }
52 :
53 : template <unsigned Dim>
54 48 : void FieldLayout<Dim>::updateLayout(const std::vector<NDIndex<Dim>>& domains) {
55 48 : if (domains.empty()) {
56 0 : return;
57 : }
58 :
59 144 : for (unsigned int i = 0; i < domains.size(); i++) {
60 192 : hLocalDomains_m(i) = domains[i];
61 : }
62 :
63 48 : findNeighbors();
64 :
65 48 : Kokkos::deep_copy(dLocalDomains_m, hLocalDomains_m);
66 :
67 48 : calcWidths();
68 : }
69 :
70 : template <unsigned Dim>
71 1876 : void FieldLayout<Dim>::initialize(const NDIndex<Dim>& domain, std::array<bool, Dim> isParallel,
72 : bool isAllPeriodic) {
73 1876 : int nRanks = comm.size();
74 :
75 1876 : gDomain_m = domain;
76 :
77 1876 : isAllPeriodic_m = isAllPeriodic;
78 :
79 1876 : isParallelDim_m = isParallel;
80 :
81 1876 : if (nRanks < 2) {
82 0 : Kokkos::resize(dLocalDomains_m, nRanks);
83 0 : Kokkos::resize(hLocalDomains_m, nRanks);
84 0 : hLocalDomains_m(0) = domain;
85 0 : Kokkos::deep_copy(dLocalDomains_m, hLocalDomains_m);
86 0 : return;
87 : }
88 :
89 : /* Check in to how many local domains we can partition the given domain.
90 : If there are more ranks than local domains, we throw an error, ippl does not allow this.
91 : */
92 1876 : long totparelems = 1;
93 7808 : for (unsigned d = 0; d < Dim; ++d) {
94 5932 : totparelems *= isParallelDim_m[d] ? domain[d].length() : 1;
95 : }
96 :
97 1876 : if (totparelems < nRanks) {
98 0 : throw std::runtime_error("FieldLayout:initialize: domain can only be partitioned in to "
99 0 : + std::to_string(totparelems) + " local domains, but there are " + std::to_string(nRanks) + " ranks, decrease the number of ranks or increase the domain.");
100 : }
101 :
102 1876 : Kokkos::resize(dLocalDomains_m, nRanks);
103 1876 : Kokkos::resize(hLocalDomains_m, nRanks);
104 :
105 : detail::Partitioner<Dim> partitioner;
106 1876 : partitioner.split(domain, hLocalDomains_m, isParallel, nRanks);
107 :
108 1876 : findNeighbors();
109 :
110 1876 : Kokkos::deep_copy(dLocalDomains_m, hLocalDomains_m);
111 :
112 1876 : calcWidths();
113 : }
114 :
115 : template <unsigned Dim>
116 4164 : const typename FieldLayout<Dim>::NDIndex_t& FieldLayout<Dim>::getLocalNDIndex() const {
117 8328 : return hLocalDomains_m(comm.rank());
118 : }
119 :
120 : template <unsigned Dim>
121 : const typename FieldLayout<Dim>::NDIndex_t& FieldLayout<Dim>::getLocalNDIndex(int rank) const {
122 : // Rank must be valid for the communicator
123 : PAssert(rank >= 0 && rank < comm.size());
124 :
125 : return hLocalDomains_m(rank);
126 : }
127 :
128 : template <unsigned Dim>
129 7260 : const typename FieldLayout<Dim>::host_mirror_type FieldLayout<Dim>::getHostLocalDomains()
130 : const {
131 7260 : return hLocalDomains_m;
132 : }
133 :
134 : template <unsigned Dim>
135 : const typename FieldLayout<Dim>::view_type FieldLayout<Dim>::getDeviceLocalDomains() const {
136 : return dLocalDomains_m;
137 : }
138 :
139 : template <unsigned Dim>
140 996 : const typename FieldLayout<Dim>::neighbor_list& FieldLayout<Dim>::getNeighbors() const {
141 996 : return neighbors_m;
142 : }
143 :
144 : template <unsigned Dim>
145 684 : const typename FieldLayout<Dim>::neighbor_range_list& FieldLayout<Dim>::getNeighborsSendRange()
146 : const {
147 684 : return neighborsSendRange_m;
148 : }
149 :
150 : template <unsigned Dim>
151 684 : const typename FieldLayout<Dim>::neighbor_range_list& FieldLayout<Dim>::getNeighborsRecvRange()
152 : const {
153 684 : return neighborsRecvRange_m;
154 : }
155 :
156 : template <unsigned Dim>
157 0 : void FieldLayout<Dim>::write(std::ostream& out) const {
158 0 : if (comm.rank() > 0) {
159 0 : return;
160 : }
161 :
162 0 : out << "Domain = " << gDomain_m << "\n"
163 0 : << "Total number of boxes = " << hLocalDomains_m.size() << "\n";
164 :
165 : using size_type = typename host_mirror_type::size_type;
166 0 : for (size_type i = 0; i < hLocalDomains_m.size(); ++i) {
167 0 : out << " Box " << i << " " << hLocalDomains_m(i) << "\n";
168 : }
169 : }
170 :
171 : template <unsigned Dim>
172 2476 : void FieldLayout<Dim>::calcWidths() {
173 : // initialize widths first
174 10508 : for (unsigned int d = 0; d < Dim; ++d) {
175 8032 : minWidth_m[d] = gDomain_m[d].length();
176 : }
177 :
178 : using size_type = typename host_mirror_type::size_type;
179 7428 : for (size_type i = 0; i < hLocalDomains_m.size(); ++i) {
180 4952 : const NDIndex_t& dom = hLocalDomains_m(i);
181 21016 : for (unsigned int d = 0; d < Dim; ++d) {
182 16064 : if ((unsigned int)dom[d].length() < minWidth_m[d]) {
183 2484 : minWidth_m[d] = dom[d].length();
184 : }
185 : }
186 : }
187 2476 : }
188 :
189 : template <unsigned Dim>
190 10640 : void FieldLayout<Dim>::findPeriodicNeighbors(const int nghost, const NDIndex<Dim>& localDomain,
191 : NDIndex<Dim>& grown, NDIndex<Dim>& neighborDomain,
192 : const int rank,
193 : std::map<unsigned int, int>& offsets, unsigned d0,
194 : unsigned codim) {
195 16464 : for (unsigned int d = d0; d < Dim; ++d) {
196 : // 0 - check upper boundary
197 : // 1 - check lower boundary
198 17472 : for (int k = 0; k < 2; ++k) {
199 11648 : auto offset = offsets[d] = getPeriodicOffset(localDomain, d, k);
200 11648 : if (offset == 0) {
201 96 : continue;
202 : }
203 :
204 11552 : grown[d] += offset;
205 :
206 11552 : if (grown.touches(neighborDomain)) {
207 11552 : auto intersect = grown.intersect(neighborDomain);
208 55248 : for (auto& [d, offset] : offsets) {
209 43696 : neighborDomain[d] -= offset;
210 : }
211 11552 : addNeighbors(grown, localDomain, neighborDomain, intersect, nghost, rank);
212 55248 : for (auto& [d, offset] : offsets) {
213 43696 : neighborDomain[d] += offset;
214 : }
215 : }
216 11552 : if (codim + 1 < Dim) {
217 10544 : findPeriodicNeighbors(nghost, localDomain, grown, neighborDomain, rank, offsets,
218 : d + 1, codim + 1);
219 : }
220 :
221 11552 : grown[d] -= offset;
222 11552 : offsets.erase(d);
223 : }
224 : }
225 10640 : }
226 :
227 : template <unsigned Dim>
228 2476 : void FieldLayout<Dim>::findNeighbors(int nghost) {
229 : /* We need to reset the neighbor list
230 : * and its ranges because of the repartitioner.
231 : */
232 373452 : for (size_t i = 0; i < detail::countHypercubes(Dim) - 1; i++) {
233 370976 : neighbors_m[i].clear();
234 370976 : neighborsSendRange_m[i].clear();
235 370976 : neighborsRecvRange_m[i].clear();
236 : }
237 :
238 2476 : int myRank = comm.rank();
239 :
240 : // get my local box
241 2476 : auto& nd = hLocalDomains_m[myRank];
242 :
243 : // grow the box by nghost cells in each dimension
244 2476 : auto gnd = nd.grow(nghost);
245 :
246 2476 : static IpplTimings::TimerRef findInternalNeighborsTimer =
247 134 : IpplTimings::getTimer("findInternal");
248 2476 : static IpplTimings::TimerRef findPeriodicNeighborsTimer =
249 134 : IpplTimings::getTimer("findPeriodic");
250 7428 : for (int rank = 0; rank < comm.size(); ++rank) {
251 4952 : if (rank == myRank) {
252 : // do not compare with my domain
253 2476 : continue;
254 : }
255 :
256 2476 : auto& ndNeighbor = hLocalDomains_m[rank];
257 2476 : IpplTimings::startTimer(findInternalNeighborsTimer);
258 : // For inter-processor neighbors
259 2476 : if (gnd.touches(ndNeighbor)) {
260 2476 : auto intersect = gnd.intersect(ndNeighbor);
261 2476 : addNeighbors(gnd, nd, ndNeighbor, intersect, nghost, rank);
262 : }
263 2476 : IpplTimings::stopTimer(findInternalNeighborsTimer);
264 :
265 2476 : IpplTimings::startTimer(findPeriodicNeighborsTimer);
266 2476 : if (isAllPeriodic_m) {
267 96 : std::map<unsigned int, int> offsets;
268 96 : findPeriodicNeighbors(nghost, nd, gnd, ndNeighbor, rank, offsets);
269 96 : }
270 2476 : IpplTimings::stopTimer(findPeriodicNeighborsTimer);
271 : }
272 2476 : }
273 :
274 : template <unsigned Dim>
275 14028 : void FieldLayout<Dim>::addNeighbors(const NDIndex_t& gnd, const NDIndex_t& nd,
276 : const NDIndex_t& ndNeighbor, const NDIndex_t& intersect,
277 : int nghost, int rank) {
278 : bound_type rangeSend, rangeRecv;
279 14028 : rangeSend = getBounds(nd, ndNeighbor, nd, nghost);
280 :
281 14028 : rangeRecv = getBounds(ndNeighbor, nd, nd, nghost);
282 :
283 14028 : int index = 0;
284 85884 : for (unsigned d = 0, digit = 1; d < Dim; d++, digit *= 3) {
285 : // For each dimension, check what kind of intersection we have and construct
286 : // an index for the component based on its base-3 representation with the
287 : // following properties for each digit:
288 : // 0 - touching the lower axis value
289 : // 1 - touching the upper axis value
290 : // 2 - parallel to the axis
291 71856 : if (intersect[d].length() == 1) {
292 49036 : if (gnd[d].first() != intersect[d].first()) {
293 24518 : index += digit;
294 : }
295 : } else {
296 22820 : index += 2 * digit;
297 : }
298 : }
299 14028 : neighbors_m[index].push_back(rank);
300 14028 : neighborsSendRange_m[index].push_back(rangeSend);
301 14028 : neighborsRecvRange_m[index].push_back(rangeRecv);
302 14028 : }
303 :
304 : template <unsigned Dim>
305 28056 : typename FieldLayout<Dim>::bound_type FieldLayout<Dim>::getBounds(const NDIndex_t& nd1,
306 : const NDIndex_t& nd2,
307 : const NDIndex_t& offset,
308 : int nghost) {
309 28056 : NDIndex<Dim> gnd = nd2.grow(nghost);
310 :
311 28056 : NDIndex<Dim> overlap = gnd.intersect(nd1);
312 :
313 : bound_type intersect;
314 :
315 : /* Obtain the intersection bounds with local ranges of the view.
316 : * Add "+1" to the upper bound since Kokkos loops always to "< extent".
317 : */
318 171768 : for (size_t i = 0; i < Dim; ++i) {
319 143712 : intersect.lo[i] = overlap[i].first() - offset[i].first() /*offset*/ + nghost;
320 143712 : intersect.hi[i] = overlap[i].last() - offset[i].first() /*offset*/ + nghost + 1;
321 : }
322 :
323 55176 : return intersect;
324 : }
325 :
326 : template <unsigned Dim>
327 11648 : int FieldLayout<Dim>::getPeriodicOffset(const NDIndex_t& nd, const unsigned int d,
328 : const int k) {
329 11648 : switch (k) {
330 5824 : case 0:
331 5824 : if (nd[d].max() == gDomain_m[d].max()) {
332 5776 : return -gDomain_m[d].length();
333 : }
334 48 : break;
335 5824 : case 1:
336 5824 : if (nd[d].min() == gDomain_m[d].min()) {
337 5776 : return gDomain_m[d].length();
338 : }
339 48 : break;
340 0 : default:
341 0 : throw IpplException("FieldLayout:getPeriodicOffset", "k has to be either 0 or 1");
342 : }
343 96 : return 0;
344 : }
345 :
346 : } // namespace ippl
|
