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 0 : 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 0 : int match = 0;
29 0 : for (unsigned d = 1; d < detail::countHypercubes(Dim); d *= 3) {
30 0 : match += digit_swap[index % 3] * d;
31 0 : index /= 3;
32 : }
33 0 : return match;
34 : }
35 :
36 : template <unsigned Dim>
37 934 : FieldLayout<Dim>::FieldLayout(const mpi::Communicator& communicator)
38 934 : : comm(communicator)
39 934 : , dLocalDomains_m("local domains (device)", 0)
40 1868 : , hLocalDomains_m(Kokkos::create_mirror_view(dLocalDomains_m)) {
41 3878 : for (unsigned int d = 0; d < Dim; ++d) {
42 2944 : minWidth_m[d] = 0;
43 : }
44 934 : }
45 :
46 : template <unsigned Dim>
47 654 : FieldLayout<Dim>::FieldLayout(mpi::Communicator communicator, const NDIndex<Dim>& domain,
48 : std::array<bool, Dim> isParallel, bool isAllPeriodic)
49 654 : : FieldLayout(communicator) {
50 654 : initialize(domain, isParallel, isAllPeriodic);
51 654 : }
52 :
53 : template <unsigned Dim>
54 24 : void FieldLayout<Dim>::updateLayout(const std::vector<NDIndex<Dim>>& domains) {
55 24 : if (domains.empty()) {
56 0 : return;
57 : }
58 :
59 48 : for (unsigned int i = 0; i < domains.size(); i++) {
60 48 : hLocalDomains_m(i) = domains[i];
61 : }
62 :
63 24 : findNeighbors();
64 :
65 24 : Kokkos::deep_copy(dLocalDomains_m, hLocalDomains_m);
66 :
67 24 : calcWidths();
68 : }
69 :
70 : template <unsigned Dim>
71 654 : void FieldLayout<Dim>::initialize(const NDIndex<Dim>& domain, std::array<bool, Dim> isParallel,
72 : bool isAllPeriodic) {
73 654 : int nRanks = comm.size();
74 :
75 654 : gDomain_m = domain;
76 :
77 654 : isAllPeriodic_m = isAllPeriodic;
78 :
79 654 : isParallelDim_m = isParallel;
80 :
81 654 : if (nRanks < 2) {
82 654 : Kokkos::resize(dLocalDomains_m, nRanks);
83 654 : Kokkos::resize(hLocalDomains_m, nRanks);
84 654 : hLocalDomains_m(0) = domain;
85 654 : Kokkos::deep_copy(dLocalDomains_m, hLocalDomains_m);
86 654 : return;
87 : }
88 :
89 : /* Check to see if we have too few elements to partition. If so, reduce
90 : * the number of ranks (if necessary) to just the number of elements along
91 : * parallel dims.
92 : */
93 0 : long totparelems = 1;
94 0 : for (unsigned d = 0; d < Dim; ++d) {
95 0 : totparelems *= domain[d].length();
96 : }
97 :
98 0 : if (totparelems < nRanks) {
99 0 : nRanks = totparelems;
100 : }
101 :
102 0 : Kokkos::resize(dLocalDomains_m, nRanks);
103 0 : Kokkos::resize(hLocalDomains_m, nRanks);
104 :
105 : detail::Partitioner<Dim> partitioner;
106 0 : partitioner.split(domain, hLocalDomains_m, isParallel, nRanks);
107 :
108 0 : findNeighbors();
109 :
110 0 : Kokkos::deep_copy(dLocalDomains_m, hLocalDomains_m);
111 :
112 0 : calcWidths();
113 : }
114 :
115 : template <unsigned Dim>
116 1254 : const typename FieldLayout<Dim>::NDIndex_t& FieldLayout<Dim>::getLocalNDIndex(int rank) const {
117 2508 : return hLocalDomains_m(rank > 0 ? rank : comm.rank());
118 : }
119 :
120 : template <unsigned Dim>
121 1368 : const typename FieldLayout<Dim>::host_mirror_type FieldLayout<Dim>::getHostLocalDomains()
122 : const {
123 1368 : return hLocalDomains_m;
124 : }
125 :
126 : template <unsigned Dim>
127 : const typename FieldLayout<Dim>::view_type FieldLayout<Dim>::getDeviceLocalDomains() const {
128 : return dLocalDomains_m;
129 : }
130 :
131 : template <unsigned Dim>
132 24 : const typename FieldLayout<Dim>::neighbor_list& FieldLayout<Dim>::getNeighbors() const {
133 24 : return neighbors_m;
134 : }
135 :
136 : template <unsigned Dim>
137 0 : const typename FieldLayout<Dim>::neighbor_range_list& FieldLayout<Dim>::getNeighborsSendRange()
138 : const {
139 0 : return neighborsSendRange_m;
140 : }
141 :
142 : template <unsigned Dim>
143 0 : const typename FieldLayout<Dim>::neighbor_range_list& FieldLayout<Dim>::getNeighborsRecvRange()
144 : const {
145 0 : return neighborsRecvRange_m;
146 : }
147 :
148 : template <unsigned Dim>
149 0 : void FieldLayout<Dim>::write(std::ostream& out) const {
150 0 : if (comm.rank() > 0) {
151 0 : return;
152 : }
153 :
154 0 : out << "Domain = " << gDomain_m << "\n"
155 0 : << "Total number of boxes = " << hLocalDomains_m.size() << "\n";
156 :
157 : using size_type = typename host_mirror_type::size_type;
158 0 : for (size_type i = 0; i < hLocalDomains_m.size(); ++i) {
159 0 : out << " Box " << i << " " << hLocalDomains_m(i) << "\n";
160 : }
161 : }
162 :
163 : template <unsigned Dim>
164 24 : void FieldLayout<Dim>::calcWidths() {
165 : // initialize widths first
166 108 : for (unsigned int d = 0; d < Dim; ++d) {
167 84 : minWidth_m[d] = gDomain_m[d].length();
168 : }
169 :
170 : using size_type = typename host_mirror_type::size_type;
171 48 : for (size_type i = 0; i < hLocalDomains_m.size(); ++i) {
172 24 : const NDIndex_t& dom = hLocalDomains_m(i);
173 108 : for (unsigned int d = 0; d < Dim; ++d) {
174 84 : if ((unsigned int)dom[d].length() < minWidth_m[d]) {
175 0 : minWidth_m[d] = dom[d].length();
176 : }
177 : }
178 : }
179 24 : }
180 :
181 : template <unsigned Dim>
182 0 : void FieldLayout<Dim>::findPeriodicNeighbors(const int nghost, const NDIndex<Dim>& localDomain,
183 : NDIndex<Dim>& grown, NDIndex<Dim>& neighborDomain,
184 : const int rank,
185 : std::map<unsigned int, int>& offsets, unsigned d0,
186 : unsigned codim) {
187 0 : for (unsigned int d = d0; d < Dim; ++d) {
188 : // 0 - check upper boundary
189 : // 1 - check lower boundary
190 0 : for (int k = 0; k < 2; ++k) {
191 0 : auto offset = offsets[d] = getPeriodicOffset(localDomain, d, k);
192 0 : if (offset == 0) {
193 0 : continue;
194 : }
195 :
196 0 : grown[d] += offset;
197 :
198 0 : if (grown.touches(neighborDomain)) {
199 0 : auto intersect = grown.intersect(neighborDomain);
200 0 : for (auto& [d, offset] : offsets) {
201 0 : neighborDomain[d] -= offset;
202 : }
203 0 : addNeighbors(grown, localDomain, neighborDomain, intersect, nghost, rank);
204 0 : for (auto& [d, offset] : offsets) {
205 0 : neighborDomain[d] += offset;
206 : }
207 : }
208 0 : if (codim + 1 < Dim) {
209 0 : findPeriodicNeighbors(nghost, localDomain, grown, neighborDomain, rank, offsets,
210 : d + 1, codim + 1);
211 : }
212 :
213 0 : grown[d] -= offset;
214 0 : offsets.erase(d);
215 : }
216 : }
217 0 : }
218 :
219 : template <unsigned Dim>
220 24 : void FieldLayout<Dim>::findNeighbors(int nghost) {
221 : /* We need to reset the neighbor list
222 : * and its ranges because of the repartitioner.
223 : */
224 4368 : for (size_t i = 0; i < detail::countHypercubes(Dim) - 1; i++) {
225 4344 : neighbors_m[i].clear();
226 4344 : neighborsSendRange_m[i].clear();
227 4344 : neighborsRecvRange_m[i].clear();
228 : }
229 :
230 24 : int myRank = comm.rank();
231 :
232 : // get my local box
233 24 : auto& nd = hLocalDomains_m[myRank];
234 :
235 : // grow the box by nghost cells in each dimension
236 24 : auto gnd = nd.grow(nghost);
237 :
238 24 : static IpplTimings::TimerRef findInternalNeighborsTimer =
239 24 : IpplTimings::getTimer("findInternal");
240 24 : static IpplTimings::TimerRef findPeriodicNeighborsTimer =
241 24 : IpplTimings::getTimer("findPeriodic");
242 48 : for (int rank = 0; rank < comm.size(); ++rank) {
243 24 : if (rank == myRank) {
244 : // do not compare with my domain
245 24 : continue;
246 : }
247 :
248 0 : auto& ndNeighbor = hLocalDomains_m[rank];
249 0 : IpplTimings::startTimer(findInternalNeighborsTimer);
250 : // For inter-processor neighbors
251 0 : if (gnd.touches(ndNeighbor)) {
252 0 : auto intersect = gnd.intersect(ndNeighbor);
253 0 : addNeighbors(gnd, nd, ndNeighbor, intersect, nghost, rank);
254 : }
255 0 : IpplTimings::stopTimer(findInternalNeighborsTimer);
256 :
257 0 : IpplTimings::startTimer(findPeriodicNeighborsTimer);
258 0 : if (isAllPeriodic_m) {
259 0 : std::map<unsigned int, int> offsets;
260 0 : findPeriodicNeighbors(nghost, nd, gnd, ndNeighbor, rank, offsets);
261 0 : }
262 0 : IpplTimings::stopTimer(findPeriodicNeighborsTimer);
263 : }
264 24 : }
265 :
266 : template <unsigned Dim>
267 0 : void FieldLayout<Dim>::addNeighbors(const NDIndex_t& gnd, const NDIndex_t& nd,
268 : const NDIndex_t& ndNeighbor, const NDIndex_t& intersect,
269 : int nghost, int rank) {
270 : bound_type rangeSend, rangeRecv;
271 0 : rangeSend = getBounds(nd, ndNeighbor, nd, nghost);
272 :
273 0 : rangeRecv = getBounds(ndNeighbor, nd, nd, nghost);
274 :
275 0 : int index = 0;
276 0 : for (unsigned d = 0, digit = 1; d < Dim; d++, digit *= 3) {
277 : // For each dimension, check what kind of intersection we have and construct
278 : // an index for the component based on its base-3 representation with the
279 : // following properties for each digit:
280 : // 0 - touching the lower axis value
281 : // 1 - touching the upper axis value
282 : // 2 - parallel to the axis
283 0 : if (intersect[d].length() == 1) {
284 0 : if (gnd[d].first() != intersect[d].first()) {
285 0 : index += digit;
286 : }
287 : } else {
288 0 : index += 2 * digit;
289 : }
290 : }
291 0 : neighbors_m[index].push_back(rank);
292 0 : neighborsSendRange_m[index].push_back(rangeSend);
293 0 : neighborsRecvRange_m[index].push_back(rangeRecv);
294 0 : }
295 :
296 : template <unsigned Dim>
297 0 : typename FieldLayout<Dim>::bound_type FieldLayout<Dim>::getBounds(const NDIndex_t& nd1,
298 : const NDIndex_t& nd2,
299 : const NDIndex_t& offset,
300 : int nghost) {
301 0 : NDIndex<Dim> gnd = nd2.grow(nghost);
302 :
303 0 : NDIndex<Dim> overlap = gnd.intersect(nd1);
304 :
305 : bound_type intersect;
306 :
307 : /* Obtain the intersection bounds with local ranges of the view.
308 : * Add "+1" to the upper bound since Kokkos loops always to "< extent".
309 : */
310 0 : for (size_t i = 0; i < Dim; ++i) {
311 0 : intersect.lo[i] = overlap[i].first() - offset[i].first() /*offset*/ + nghost;
312 0 : intersect.hi[i] = overlap[i].last() - offset[i].first() /*offset*/ + nghost + 1;
313 : }
314 :
315 0 : return intersect;
316 : }
317 :
318 : template <unsigned Dim>
319 0 : int FieldLayout<Dim>::getPeriodicOffset(const NDIndex_t& nd, const unsigned int d,
320 : const int k) {
321 0 : switch (k) {
322 0 : case 0:
323 0 : if (nd[d].max() == gDomain_m[d].max()) {
324 0 : return -gDomain_m[d].length();
325 : }
326 0 : break;
327 0 : case 1:
328 0 : if (nd[d].min() == gDomain_m[d].min()) {
329 0 : return gDomain_m[d].length();
330 : }
331 0 : break;
332 0 : default:
333 0 : throw IpplException("FieldLayout:getPeriodicOffset", "k has to be either 0 or 1");
334 : }
335 0 : return 0;
336 : }
337 :
338 : } // namespace ippl
|
