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 : #ifndef IPPL_FIELD_LAYOUT_H
9 : #define IPPL_FIELD_LAYOUT_H
10 :
11 : #include <array>
12 : #include <iostream>
13 : #include <map>
14 : #include <vector>
15 :
16 : #include "Types/ViewTypes.h"
17 :
18 : #include "Communicate/Communicator.h"
19 : #include "Index/NDIndex.h"
20 : #include "Partition/Partitioner.h"
21 :
22 : namespace ippl {
23 :
24 : template <unsigned Dim>
25 : class FieldLayout;
26 :
27 : template <unsigned Dim>
28 : std::ostream& operator<<(std::ostream&, const FieldLayout<Dim>&);
29 :
30 : // enumeration used to describe a hypercube's relation to
31 : // a particular axis in a given bounded domain
32 : enum e_cube_tag {
33 : UPPER = 0,
34 : LOWER = 1,
35 : IS_PARALLEL = 2
36 : };
37 :
38 : namespace detail {
39 : /*!
40 : * Counts the hypercubes in a given dimension
41 : * @param dim the dimension
42 : * @return 3^n
43 : */
44 4368 : constexpr unsigned int countHypercubes(unsigned int dim) {
45 4368 : unsigned int ret = 1;
46 28428 : for (unsigned int d = 0; d < dim; d++) {
47 24060 : ret *= 3;
48 : }
49 4368 : return ret;
50 : }
51 :
52 : /*!
53 : * Compile-time evaluation of x!
54 : * @param x input value
55 : * @return x factorial
56 : */
57 : constexpr unsigned int factorial(unsigned x) {
58 : return x == 0 ? 1 : x * factorial(x - 1);
59 : }
60 :
61 : /*!
62 : * Compile-time evaluation of binomial coefficients, aka
63 : * elements of Pascal's triangle, aka choices from combinatorics, etc
64 : * @param a number of options
65 : * @param b number of choices
66 : * @return a choose b
67 : */
68 : constexpr unsigned int binomialCoefficient(unsigned a, unsigned b) {
69 : return factorial(a) / (factorial(b) * factorial(a - b));
70 : }
71 :
72 : /*!
73 : * Compile-time evaluation of the number of hypercubes of dimension m in a hypercube
74 : * of dimension Dim
75 : * @tparam Dim parent hypercube dimension
76 : * @param m sub-hypercube dimension
77 : * @return The number of m-cubes in an n-cube
78 : */
79 : template <unsigned Dim>
80 : constexpr unsigned int countCubes(unsigned m) {
81 : return (1 << (Dim - m)) * binomialCoefficient(Dim, m);
82 : }
83 :
84 : /*!
85 : * Determines whether a facet is on the upper boundary
86 : * of its domain. For lower dimension hypercubes, determines
87 : * whether the component is on the upper boundary of
88 : * the domain along at least one axis
89 : * @param face the hypercube's index
90 : * @return Whether it touches any upper boundary
91 : */
92 : bool isUpper(unsigned int face);
93 :
94 : /*!
95 : * Determine the axis perpendicular to a given facet
96 : * (throws an exception if the index does not correspond
97 : * to a facet)
98 : * @param face the facet's index
99 : * @return The index of the axis perpendicular to that facet
100 : */
101 : unsigned int getFaceDim(unsigned int face);
102 :
103 : /*!
104 : * Converts between ternary encoding and face set indexing
105 : * @tparam Dim the number of dimensions
106 : * @param index the ternary-encoded index of a facet in [0, 3^Dim)
107 : * @return The index of that facet in a set of faces in [0, 2*Dim)
108 : */
109 : template <unsigned Dim>
110 : unsigned int indexToFace(unsigned int index) {
111 : // facets are group low/high by axis
112 : unsigned int axis = index / 2;
113 : // the digit to subtract is determined by whether the index describes an upper
114 : // face (even index) or lower face (odd index) and that digit's position is
115 : // determined by the axis of the face
116 : unsigned int toRemove = (2 - index % 2) * countHypercubes(axis);
117 : // start with all 2s (in base 3) and change the correct digit to get the encoded face
118 : return countHypercubes(Dim) - 1 - toRemove;
119 : }
120 :
121 : /*!
122 : * Computes the ternary-encoded index of a hypercube
123 : * @tparam Dim the number of dimensions in the full hypercube
124 : * @tparam CubeTags... variadic argument list, must be all e_cube_tag
125 : * @param tag(s...) the tags describing the hypercube of interest
126 : * @return The index of the desired hypercube
127 : */
128 : template <
129 : unsigned Dim, typename... CubeTags,
130 : typename = std::enable_if_t<sizeof...(CubeTags) == Dim - 1>,
131 : typename = std::enable_if_t<std::conjunction_v<std::is_same<e_cube_tag, CubeTags>...>>>
132 0 : unsigned int getCube(e_cube_tag tag, CubeTags... tags) {
133 : if constexpr (Dim == 1) {
134 0 : return tag;
135 : } else {
136 0 : return tag + 3 * getCube<Dim - 1>(tags...);
137 : }
138 : }
139 :
140 : /*!
141 : * Utility function for getFace
142 : */
143 : template <size_t... Idx>
144 : unsigned int getFace_impl(const std::array<e_cube_tag, sizeof...(Idx)>& args,
145 : const std::index_sequence<Idx...>&) {
146 : return getCube<sizeof...(Idx)>(args[Idx]...);
147 : }
148 :
149 : /*!
150 : * Convenience alias for getCube for getting facets
151 : * @tparam Dim the number of dimensions in the parent hypercube
152 : * @param axis the axis perpendicular to the facet
153 : * @param side whether the facet is an upper or lower facet
154 : * @return The index of the facet
155 : */
156 : template <unsigned Dim>
157 : unsigned int getFace(unsigned int axis, e_cube_tag side) {
158 : std::array<e_cube_tag, Dim> args;
159 : args.fill(IS_PARALLEL);
160 : args[axis] = side;
161 : return getFace_impl(args, std::make_index_sequence<Dim>{});
162 : }
163 : } // namespace detail
164 :
165 : template <unsigned Dim>
166 : class FieldLayout {
167 : public:
168 : using NDIndex_t = NDIndex<Dim>;
169 : using view_type = typename detail::ViewType<NDIndex_t, 1>::view_type;
170 : using host_mirror_type = typename view_type::host_mirror_type;
171 :
172 : struct bound_type {
173 : // lower bounds (ordering: x, y, z, ...)
174 : std::array<long, Dim> lo;
175 : // upper bounds (ordering: x, y, z, ...)
176 : std::array<long, Dim> hi;
177 :
178 : /*!
179 : * Compute the size of the region described by the bounds
180 : * @return Product of the axial dimensions of the region
181 : */
182 0 : long size() const {
183 0 : long total = 1;
184 0 : for (unsigned d = 0; d < Dim; d++) {
185 0 : total *= hi[d] - lo[d];
186 : }
187 0 : return total;
188 : }
189 : };
190 :
191 : using rank_list = std::vector<int>;
192 : using bounds_list = std::vector<bound_type>;
193 :
194 : using neighbor_list = std::array<rank_list, detail::countHypercubes(Dim) - 1>;
195 : using neighbor_range_list = std::array<bounds_list, detail::countHypercubes(Dim) - 1>;
196 :
197 : /*!
198 : * Default constructor, which should only be used if you are going to
199 : * call 'initialize' soon after (before using in any context)
200 : */
201 : FieldLayout(const mpi::Communicator& = MPI_COMM_WORLD);
202 :
203 : FieldLayout(mpi::Communicator, const NDIndex<Dim>& domain, std::array<bool, Dim> decomp,
204 : bool isAllPeriodic = false);
205 :
206 : // Destructor: Everything deletes itself automatically ... the base
207 : // class destructors inform all the FieldLayoutUser's we're going away.
208 996 : virtual ~FieldLayout() = default;
209 :
210 : // Initialization functions, only to be called by the user of FieldLayout
211 : // objects when the FieldLayout was created using the default constructor;
212 : // otherwise these are only called internally by the various non-default
213 : // FieldLayout constructors:
214 :
215 : void initialize(const NDIndex<Dim>& domain, std::array<bool, Dim> decomp,
216 : bool isAllPeriodic = false);
217 :
218 : // Return the domain.
219 1722 : const NDIndex<Dim>& getDomain() const { return gDomain_m; }
220 :
221 : // Compare FieldLayouts to see if they represent the same domain; if
222 : // dimensionalities are different, the NDIndex operator==() will return
223 : // false:
224 : template <unsigned Dim2>
225 : bool operator==(const FieldLayout<Dim2>& x) const {
226 : return gDomain_m == x.getDomain();
227 : }
228 :
229 : bool operator==(const FieldLayout<Dim>& x) const {
230 : for (unsigned int i = 0; i < Dim; ++i) {
231 : if (hLocalDomains_m(comm.rank())[i] != x.getLocalNDIndex()[i]) {
232 : return false;
233 : }
234 : }
235 : return true;
236 : }
237 :
238 : // for the requested dimension, report if the distribution is
239 : // SERIAL or PARALLEL
240 : bool getDistribution(unsigned int d) const {
241 : return minWidth_m[d] == (unsigned int)gDomain_m[d].length();
242 : }
243 :
244 : // for the requested dimension, report if the distribution was requested to
245 : // be SERIAL or PARALLEL
246 0 : std::array<bool, Dim> isParallel() const { return isParallelDim_m; }
247 :
248 : const NDIndex_t& getLocalNDIndex(int rank = -1) const;
249 :
250 : const host_mirror_type getHostLocalDomains() const;
251 :
252 : const view_type getDeviceLocalDomains() const;
253 :
254 : /*!
255 : * Get a list of all the neighbors, arranged by ternary encoding
256 : * of the hypercubes
257 : * @return List of list of neighbor ranks touching each boundary component
258 : */
259 : const neighbor_list& getNeighbors() const;
260 :
261 : /*!
262 : * Get the domain ranges corresponding to regions that should be sent
263 : * to neighbor ranks
264 : * @return Ranges to send
265 : */
266 : const neighbor_range_list& getNeighborsSendRange() const;
267 :
268 : /*!
269 : * Get the domain ranges corresponding to regions that should be received
270 : * from neighbor ranks
271 : * @return Ranges to receive
272 : */
273 : const neighbor_range_list& getNeighborsRecvRange() const;
274 :
275 : /*!
276 : * Given the index of a hypercube, find the index of the opposite hypercube,
277 : * i.e. the component with the same codimension belonging to a neighboring domain
278 : * that touches the hypercube with the given index, as determined by the
279 : * ternary encoding for hypercubes.
280 : *
281 : * For neighbor communication, the opposite component is the one that receives
282 : * sent data or sends us data to receive for a given component.
283 : *
284 : * The matching index is given by swapping alls 1s for 0s and vice versa in
285 : * the ternary encoding, while keeping the 2s unchanged. This can be understood
286 : * from the fact that if the local component is on the upper boundary of the local
287 : * domain, the neighbor component must be on the lower boundary of its local domain,
288 : * and vice versa. The 2s are unchanged because both the local component and the
289 : * neighbor component must be parallel to the same axes, otherwise their intersection
290 : * would have lower or higher dimension than the components themselves.
291 : * @param index index of the known component
292 : * @return Index of the matching component
293 : */
294 : static int getMatchingIndex(int index);
295 :
296 : /*!
297 : * Recursively finds neighbor ranks for layouts with all periodic boundary
298 : * conditions
299 : * @param nghost number of ghost cells
300 : * @param localDomain the rank's local domain
301 : * @param grown the local domain, grown by the number of ghost cells
302 : * @param neighborDomain a candidate neighbor rank's domain
303 : * @param rank the candidate neighbor's rank
304 : * @param offsets a dictionary containing offsets along different dimensions
305 : * @param d0 the dimension from which to start checking (default 0)
306 : * @param codim the codimension of overlapping regions to check (default 0)
307 : */
308 : void findPeriodicNeighbors(const int nghost, const NDIndex<Dim>& localDomain,
309 : NDIndex<Dim>& grown, NDIndex<Dim>& neighborDomain,
310 : const int rank, std::map<unsigned int, int>& offsets,
311 : unsigned d0 = 0, unsigned codim = 0);
312 :
313 : /*!
314 : * Finds all neighboring ranks based on the field layout
315 : * @param nghost number of ghost cells (default 1)
316 : */
317 : void findNeighbors(int nghost = 1);
318 :
319 : /*!
320 : * Adds a neighbor to the neighbor list
321 : * @param gnd the local domain, including ghost cells
322 : * @param nd the local domain
323 : * @param ndNeighbor the neighbor rank's domain
324 : * @param intersect the intersection of the domains
325 : * @param nghost number of ghost cells
326 : * @param rank the neighbor's rank
327 : */
328 : void addNeighbors(const NDIndex_t& gnd, const NDIndex_t& nd, const NDIndex_t& ndNeighbor,
329 : const NDIndex_t& intersect, int nghost, int rank);
330 :
331 : void write(std::ostream& = std::cout) const;
332 :
333 : void updateLayout(const std::vector<NDIndex_t>& domains);
334 :
335 : bool isAllPeriodic_m;
336 :
337 : mpi::Communicator comm;
338 :
339 : private:
340 : /*!
341 : * Obtain the bounds to send / receive. The second domain, i.e.,
342 : * nd2, is grown by nghost cells in each dimension in order to
343 : * figure out the intersecting cells.
344 : * @param nd1 either remote or owned domain
345 : * @param nd2 either remote or owned domain
346 : * @param offset to map global to local grid point
347 : * @param nghost number of ghost cells per dimension
348 : */
349 : bound_type getBounds(const NDIndex_t& nd1, const NDIndex_t& nd2, const NDIndex_t& offset,
350 : int nghost);
351 :
352 : int getPeriodicOffset(const NDIndex_t& nd, const unsigned int d, const int k);
353 :
354 : private:
355 : //! Global domain
356 : NDIndex_t gDomain_m;
357 :
358 : //! Local domains (device view)
359 : view_type dLocalDomains_m;
360 :
361 : //! Local domains (host mirror view)
362 : host_mirror_type hLocalDomains_m;
363 :
364 : std::array<bool, Dim> isParallelDim_m;
365 :
366 : unsigned int minWidth_m[Dim];
367 :
368 : neighbor_list neighbors_m;
369 : neighbor_range_list neighborsSendRange_m, neighborsRecvRange_m;
370 :
371 : void calcWidths();
372 : };
373 :
374 : template <unsigned Dim>
375 0 : inline std::ostream& operator<<(std::ostream& out, const FieldLayout<Dim>& f) {
376 0 : f.write(out);
377 0 : return out;
378 : }
379 : } // namespace ippl
380 :
381 : #include "FieldLayout/FieldLayout.hpp"
382 :
383 : #endif
|
