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 : : #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 : 792 : 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
|
