Branch data 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 : 732 : FieldLayout<Dim>::FieldLayout(const mpi::Communicator& communicator)
38 : 732 : : comm(communicator)
39 [ + - ]: 732 : , dLocalDomains_m("local domains (device)", 0)
40 [ + - ]: 1464 : , hLocalDomains_m(Kokkos::create_mirror_view(dLocalDomains_m)) {
41 [ + + ]: 3258 : for (unsigned int d = 0; d < Dim; ++d) {
42 : 2526 : minWidth_m[d] = 0;
43 : : }
44 : 732 : }
45 : :
46 : : template <unsigned Dim>
47 : 452 : FieldLayout<Dim>::FieldLayout(mpi::Communicator communicator, const NDIndex<Dim>& domain,
48 : : std::array<bool, Dim> isParallel, bool isAllPeriodic)
49 : 452 : : FieldLayout(communicator) {
50 [ + - ]: 452 : initialize(domain, isParallel, isAllPeriodic);
51 : 452 : }
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 : 452 : void FieldLayout<Dim>::initialize(const NDIndex<Dim>& domain, std::array<bool, Dim> isParallel,
72 : : bool isAllPeriodic) {
73 : 452 : int nRanks = comm.size();
74 : :
75 : 452 : gDomain_m = domain;
76 : :
77 : 452 : isAllPeriodic_m = isAllPeriodic;
78 : :
79 : 452 : isParallelDim_m = isParallel;
80 : :
81 [ + - ]: 452 : if (nRanks < 2) {
82 [ + - ]: 452 : Kokkos::resize(dLocalDomains_m, nRanks);
83 [ + - ]: 452 : Kokkos::resize(hLocalDomains_m, nRanks);
84 : 452 : hLocalDomains_m(0) = domain;
85 [ + - ]: 452 : Kokkos::deep_copy(dLocalDomains_m, hLocalDomains_m);
86 : 452 : 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 : 954 : const typename FieldLayout<Dim>::NDIndex_t& FieldLayout<Dim>::getLocalNDIndex(int rank) const {
117 [ + - ]: 1908 : 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
|
