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 66256 : 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 66256 : int match = 0;
29 428960 : for (unsigned d = 1; d < detail::countHypercubes(Dim); d *= 3) {
30 362704 : match += digit_swap[index % 3] * d;
31 362704 : index /= 3;
32 : }
33 66256 : return match;
34 : }
35 :
36 : template <unsigned Dim>
37 1884 : FieldLayout<Dim>::FieldLayout(const mpi::Communicator& communicator)
38 1884 : : comm(communicator)
39 1884 : , dLocalDomains_m("local domains (device)", 0)
40 3768 : , hLocalDomains_m(Kokkos::create_mirror_view(dLocalDomains_m)) {
41 7812 : for (unsigned int d = 0; d < Dim; ++d) {
42 5928 : minWidth_m[d] = 0;
43 : }
44 1884 : }
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 1324 : void FieldLayout<Dim>::initialize(const NDIndex<Dim>& domain, std::array<bool, Dim> isParallel,
72 : bool isAllPeriodic) {
73 1324 : int nRanks = comm.size();
74 :
75 1324 : gDomain_m = domain;
76 :
77 1324 : isAllPeriodic_m = isAllPeriodic;
78 :
79 1324 : isParallelDim_m = isParallel;
80 :
81 1324 : 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 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 1324 : long totparelems = 1;
94 5324 : for (unsigned d = 0; d < Dim; ++d) {
95 4000 : totparelems *= domain[d].length();
96 : }
97 :
98 1324 : if (totparelems < nRanks) {
99 0 : nRanks = totparelems;
100 : }
101 :
102 1324 : Kokkos::resize(dLocalDomains_m, nRanks);
103 1324 : Kokkos::resize(hLocalDomains_m, nRanks);
104 :
105 : detail::Partitioner<Dim> partitioner;
106 1324 : partitioner.split(domain, hLocalDomains_m, isParallel, nRanks);
107 :
108 1324 : findNeighbors();
109 :
110 1324 : Kokkos::deep_copy(dLocalDomains_m, hLocalDomains_m);
111 :
112 1324 : calcWidths();
113 : }
114 :
115 : template <unsigned Dim>
116 2900 : const typename FieldLayout<Dim>::NDIndex_t& FieldLayout<Dim>::getLocalNDIndex(int rank) const {
117 5800 : return hLocalDomains_m(rank > 0 ? rank : comm.rank());
118 : }
119 :
120 : template <unsigned Dim>
121 3848 : const typename FieldLayout<Dim>::host_mirror_type FieldLayout<Dim>::getHostLocalDomains()
122 : const {
123 3848 : 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 680 : const typename FieldLayout<Dim>::neighbor_list& FieldLayout<Dim>::getNeighbors() const {
133 680 : return neighbors_m;
134 : }
135 :
136 : template <unsigned Dim>
137 440 : const typename FieldLayout<Dim>::neighbor_range_list& FieldLayout<Dim>::getNeighborsSendRange()
138 : const {
139 440 : return neighborsSendRange_m;
140 : }
141 :
142 : template <unsigned Dim>
143 440 : const typename FieldLayout<Dim>::neighbor_range_list& FieldLayout<Dim>::getNeighborsRecvRange()
144 : const {
145 440 : 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 1372 : void FieldLayout<Dim>::calcWidths() {
165 : // initialize widths first
166 5540 : for (unsigned int d = 0; d < Dim; ++d) {
167 4168 : minWidth_m[d] = gDomain_m[d].length();
168 : }
169 :
170 : using size_type = typename host_mirror_type::size_type;
171 4116 : for (size_type i = 0; i < hLocalDomains_m.size(); ++i) {
172 2744 : const NDIndex_t& dom = hLocalDomains_m(i);
173 11080 : for (unsigned int d = 0; d < Dim; ++d) {
174 8336 : if ((unsigned int)dom[d].length() < minWidth_m[d]) {
175 1380 : minWidth_m[d] = dom[d].length();
176 : }
177 : }
178 : }
179 1372 : }
180 :
181 : template <unsigned Dim>
182 10640 : 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 16464 : for (unsigned int d = d0; d < Dim; ++d) {
188 : // 0 - check upper boundary
189 : // 1 - check lower boundary
190 17472 : for (int k = 0; k < 2; ++k) {
191 11648 : auto offset = offsets[d] = getPeriodicOffset(localDomain, d, k);
192 11648 : if (offset == 0) {
193 96 : continue;
194 : }
195 :
196 11552 : grown[d] += offset;
197 :
198 11552 : if (grown.touches(neighborDomain)) {
199 11552 : auto intersect = grown.intersect(neighborDomain);
200 55248 : for (auto& [d, offset] : offsets) {
201 43696 : neighborDomain[d] -= offset;
202 : }
203 11552 : addNeighbors(grown, localDomain, neighborDomain, intersect, nghost, rank);
204 55248 : for (auto& [d, offset] : offsets) {
205 43696 : neighborDomain[d] += offset;
206 : }
207 : }
208 11552 : if (codim + 1 < Dim) {
209 10544 : findPeriodicNeighbors(nghost, localDomain, grown, neighborDomain, rank, offsets,
210 : d + 1, codim + 1);
211 : }
212 :
213 11552 : grown[d] -= offset;
214 11552 : offsets.erase(d);
215 : }
216 : }
217 10640 : }
218 :
219 : template <unsigned Dim>
220 1372 : void FieldLayout<Dim>::findNeighbors(int nghost) {
221 : /* We need to reset the neighbor list
222 : * and its ranges because of the repartitioner.
223 : */
224 172524 : for (size_t i = 0; i < detail::countHypercubes(Dim) - 1; i++) {
225 171152 : neighbors_m[i].clear();
226 171152 : neighborsSendRange_m[i].clear();
227 171152 : neighborsRecvRange_m[i].clear();
228 : }
229 :
230 1372 : int myRank = comm.rank();
231 :
232 : // get my local box
233 1372 : auto& nd = hLocalDomains_m[myRank];
234 :
235 : // grow the box by nghost cells in each dimension
236 1372 : auto gnd = nd.grow(nghost);
237 :
238 1372 : static IpplTimings::TimerRef findInternalNeighborsTimer =
239 98 : IpplTimings::getTimer("findInternal");
240 1372 : static IpplTimings::TimerRef findPeriodicNeighborsTimer =
241 98 : IpplTimings::getTimer("findPeriodic");
242 4116 : for (int rank = 0; rank < comm.size(); ++rank) {
243 2744 : if (rank == myRank) {
244 : // do not compare with my domain
245 1372 : continue;
246 : }
247 :
248 1372 : auto& ndNeighbor = hLocalDomains_m[rank];
249 1372 : IpplTimings::startTimer(findInternalNeighborsTimer);
250 : // For inter-processor neighbors
251 1372 : if (gnd.touches(ndNeighbor)) {
252 1372 : auto intersect = gnd.intersect(ndNeighbor);
253 1372 : addNeighbors(gnd, nd, ndNeighbor, intersect, nghost, rank);
254 : }
255 1372 : IpplTimings::stopTimer(findInternalNeighborsTimer);
256 :
257 1372 : IpplTimings::startTimer(findPeriodicNeighborsTimer);
258 1372 : if (isAllPeriodic_m) {
259 96 : std::map<unsigned int, int> offsets;
260 96 : findPeriodicNeighbors(nghost, nd, gnd, ndNeighbor, rank, offsets);
261 96 : }
262 1372 : IpplTimings::stopTimer(findPeriodicNeighborsTimer);
263 : }
264 1372 : }
265 :
266 : template <unsigned Dim>
267 12924 : 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 12924 : rangeSend = getBounds(nd, ndNeighbor, nd, nghost);
272 :
273 12924 : rangeRecv = getBounds(ndNeighbor, nd, nd, nghost);
274 :
275 12924 : int index = 0;
276 80916 : 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 67992 : if (intersect[d].length() == 1) {
284 47932 : if (gnd[d].first() != intersect[d].first()) {
285 23966 : index += digit;
286 : }
287 : } else {
288 20060 : index += 2 * digit;
289 : }
290 : }
291 12924 : neighbors_m[index].push_back(rank);
292 12924 : neighborsSendRange_m[index].push_back(rangeSend);
293 12924 : neighborsRecvRange_m[index].push_back(rangeRecv);
294 12924 : }
295 :
296 : template <unsigned Dim>
297 25848 : 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 25848 : NDIndex<Dim> gnd = nd2.grow(nghost);
302 :
303 25848 : 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 161832 : for (size_t i = 0; i < Dim; ++i) {
311 135984 : intersect.lo[i] = overlap[i].first() - offset[i].first() /*offset*/ + nghost;
312 135984 : intersect.hi[i] = overlap[i].last() - offset[i].first() /*offset*/ + nghost + 1;
313 : }
314 :
315 51128 : return intersect;
316 : }
317 :
318 : template <unsigned Dim>
319 11648 : int FieldLayout<Dim>::getPeriodicOffset(const NDIndex_t& nd, const unsigned int d,
320 : const int k) {
321 11648 : switch (k) {
322 5824 : case 0:
323 5824 : if (nd[d].max() == gDomain_m[d].max()) {
324 5776 : return -gDomain_m[d].length();
325 : }
326 48 : break;
327 5824 : case 1:
328 5824 : if (nd[d].min() == gDomain_m[d].min()) {
329 5776 : return gDomain_m[d].length();
330 : }
331 48 : break;
332 0 : default:
333 0 : throw IpplException("FieldLayout:getPeriodicOffset", "k has to be either 0 or 1");
334 : }
335 96 : return 0;
336 : }
337 :
338 : } // namespace ippl
|
