MAIA bb96820c
Multiphysics at AIA
Loading...
Searching...
No Matches
structuredpartition.cpp
Go to the documentation of this file.
1// Copyright (C) 2024 The m-AIA AUTHORS
2//
3// This file is part of m-AIA (https://git.rwth-aachen.de/aia/m-AIA/m-AIA)
4//
5// SPDX-License-Identifier: LGPL-3.0-only
6
8#include <cstring>
9#include "COMM/mpioverride.h"
10#include "IO/parallelio.h"
11#include "globals.h"
12
13using namespace std;
14
15
25template <MInt nDim>
27 const MString fileName_,
28 MPI_Comm comm_,
29 const MInt noDomains_)
30 : m_noBlocks(noBlocks_),
31 m_mpiComm(comm_),
32 m_noDomains(noDomains_),
33 m_noPartitions(1),
34 m_eps(std::numeric_limits<MFloat>::epsilon()) {
35 for(MInt j = 0; j < m_noBlocks; j++) {
36 m_blockInfo.push_back(make_unique<PartitionInfo<nDim>>());
37 }
38
39 MIntScratchSpace dataSetSize(m_noBlocks * nDim, AT_, "dataSetSize");
40 MFloatScratchSpace weight(m_noBlocks, AT_, "weights");
41
42 dataSetSize.fill(0);
43
44 if(globalDomainId() == 0) {
45 ParallelIoHdf5 pio(fileName_, maia::parallel_io::PIO_APPEND, MPI_COMM_SELF);
46 for(MInt i = 0; i < m_noBlocks; i++) {
47 MString sBlockName = "/block" + std::to_string(i);
48 std::vector<ParallelIo::size_type> tmp(nDim, 0);
49 pio.getArraySize("x", sBlockName, &tmp[0]);
50 std::copy(tmp.begin(), tmp.end(), &dataSetSize[i * nDim]);
51 weight[i] = 1.0;
52 if(pio.hasAttribute("weight", sBlockName)) { // get the weighting factor
53 pio.getAttribute<MFloat>(&weight[i], "weight", sBlockName);
54 }
55 }
56
57 MPI_Bcast(&dataSetSize[0], m_noBlocks * nDim, MPI_INT, 0, m_mpiComm, AT_, "dataSetSize.getPointer()");
58 MPI_Bcast(weight.getPointer(), m_noBlocks, MPI_DOUBLE, 0, m_mpiComm, AT_, "weight.getPointer()");
59 } else {
60 MPI_Bcast(&dataSetSize[0], m_noBlocks * nDim, MPI_INT, 0, m_mpiComm, AT_, "dataSetSize.getPointer()");
61 MPI_Bcast(weight.getPointer(), m_noBlocks, MPI_DOUBLE, 0, m_mpiComm, AT_, "weight.getPointer()");
62 }
63
64 for(MInt i = 0; i < m_noBlocks; i++) {
65 m_blockInfo[i]->blockId = i;
66 for(MInt j = 0; j < nDim; j++) {
67 m_blockInfo[i]->size[j] = dataSetSize[i * nDim + j] - 1;
68 m_blockInfo[i]->weight = weight[i];
69 m_blockInfo[i]->offset[j] = 0;
70 }
71 }
72}
73
74
84template <MInt nDim>
86 ParallelIoHdf5* pio,
87 MPI_Comm comm_,
88 const MInt noDomains_)
89 : m_noBlocks(noBlocks_),
90 m_mpiComm(comm_),
91 m_noDomains(noDomains_),
92 m_noPartitions(1),
93 m_eps(std::numeric_limits<MFloat>::epsilon()) {
94 for(MInt j = 0; j < m_noBlocks; j++) {
95 m_blockInfo.push_back(make_unique<PartitionInfo<nDim>>());
96 }
97
98 ParallelIo::size_type ijk_max[3] = {0};
99 for(MInt i = 0; i < m_noBlocks; i++) {
100 m_blockInfo[i]->blockId = i;
101 MString sBlockName = "/block" + std::to_string(i);
102 pio->getArraySize("x", sBlockName, &ijk_max[0]);
103
104 // assign ijk_max to m_blockInfo and also assign the offset
105 for(MInt j = 0; j < nDim; j++) {
106 m_blockInfo[i]->size[j] = ijk_max[j] - 1;
107 m_blockInfo[i]->offset[j] = 0;
108 }
109 }
110}
111
126template <MInt nDim>
128 const MInt blockId,
129 const MInt (&level2dimension1D)[nDim],
130 MInt& partitionCounter,
131 MInt (&beginPartitionHelper)[nDim],
132 MInt (&endPartitionHelper)[nDim],
133 MInt (&divisor)[nDim]) {
134 treeNode* myTreePointer;
135 MInt beginBlock[nDim]{0};
136 MInt endBlock[nDim]{0};
137
138 divisor[level2dimension1D[treePointer->level]] = treePointer->noLeaves;
139
140 if(treePointer->level < (nDim - 1)) {
141 beginPartitionHelper[level2dimension1D[treePointer->level]] = 0;
142 for(MInt i = 0; i < treePointer->noChilds; i++) {
143 if(i > 0) {
144 beginPartitionHelper[level2dimension1D[treePointer->level]] =
145 beginPartitionHelper[level2dimension1D[treePointer->level]] + treePointer->childs[i - 1]->noLeaves;
146 }
147 endPartitionHelper[level2dimension1D[treePointer->level]] =
148 beginPartitionHelper[level2dimension1D[treePointer->level]] + treePointer->childs[i]->noLeaves;
149 myTreePointer = treePointer->childs[i];
150 setPartitionInfoHelper(myTreePointer, blockId, level2dimension1D, partitionCounter, beginPartitionHelper,
151 endPartitionHelper, divisor);
152 }
153 } else {
154 for(MInt i = 0; i < treePointer->noChilds; i++) {
155 beginPartitionHelper[level2dimension1D[treePointer->level]] = i;
156 endPartitionHelper[level2dimension1D[treePointer->level]] = i + 1;
157 for(MInt j = 0; j < nDim; j++) {
158 beginBlock[j] = rounding(((MFloat)m_blockInfo[blockId]->size[j]) * ((MFloat)beginPartitionHelper[j])
159 / ((MFloat)divisor[j]));
160 endBlock[j] =
161 rounding(((MFloat)m_blockInfo[blockId]->size[j]) * ((MFloat)endPartitionHelper[j]) / ((MFloat)divisor[j]));
162 }
163 m_partitionInfo[partitionCounter]->weight = m_blockInfo[blockId]->weight;
164 m_partitionInfo[partitionCounter]->totalSize = 1;
165
166 for(MInt j = 0; j < nDim; j++) {
167 m_partitionInfo[partitionCounter]->size[j] = endBlock[j] - beginBlock[j] + 1;
168 m_partitionInfo[partitionCounter]->totalSize *= (m_partitionInfo[partitionCounter]->size[j]);
169 m_partitionInfo[partitionCounter]->offset[j] = beginBlock[j];
170 }
171
172 m_partitionInfo[partitionCounter]->blockId = blockId;
173 m_partitionInfo[partitionCounter]->partitionId = partitionCounter;
174 partitionCounter++;
175 }
176 }
177}
178
184template <MInt nDim>
186 treeNode* myTreeRootPointer{};
187 MInt beginPartitionHelper[nDim]{0};
188 MInt endPartitionHelper[nDim]{0};
189 MInt divisor[nDim]{0};
190
191 m_partitionInfo.clear();
192 for(MInt j = 0; j < m_noPartitions; j++) {
193 m_partitionInfo.push_back(make_unique<PartitionInfo<nDim>>());
194 }
195
196 MInt partitionCounter = 0;
197 for(MInt i = 0; i < m_noBlocks; i++) {
198 myTreeRootPointer = treeRoot[i];
199 MInt level2dimension1D[nDim]{0};
200
201 for(MInt j = 0; j < nDim; j++) {
202 level2dimension1D[j] = level2dimension(j, i);
203 }
204
205 setPartitionInfoHelper(myTreeRootPointer, i, level2dimension1D, partitionCounter, beginPartitionHelper,
206 endPartitionHelper, divisor);
207 }
208}
209
215template <MInt nDim>
217 const MInt blockId,
218 const MInt (&level2dimension1D)[nDim],
219 treeNode*& insertTreePointer) {
220 treeNode* myTreePointer{};
221
222 // compare visited nodes with current candidate insertion node
223 // by their normalized number of children
224 if(((((MFloat)treePointer->noChilds)
225 * ((MFloat)m_blockInfo[blockId]->size[level2dimension1D[insertTreePointer->level]] + 1.0))
226 < (((MFloat)insertTreePointer->noChilds)
227 * ((MFloat)m_blockInfo[blockId]->size[level2dimension1D[treePointer->level]] + 1.0)))
228 || (approx(((MFloat)(treePointer->noChilds)
229 * ((MFloat)(m_blockInfo[blockId]->size[level2dimension1D[insertTreePointer->level]]) + 1.0)),
230 ((MFloat)(insertTreePointer->noChilds)
231 * ((MFloat)(m_blockInfo[blockId]->size[level2dimension1D[treePointer->level]]) + 1.0)),
232 m_eps)
233 && (treePointer->level < insertTreePointer->level))) {
234 insertTreePointer = treePointer;
235 }
236
237 if(treePointer->level < (nDim - 1)) {
238 // look for leftmost child with least leaves
239 myTreePointer = treePointer->childs[0];
240 for(MInt i = 1; i < treePointer->noChilds; i++) {
241 if(treePointer->childs[i]->noLeaves < myTreePointer->noLeaves) {
242 myTreePointer = treePointer->childs[i];
243 }
244 }
245
246 // now traverse in the leftmost child with least leaves
247 traverseForInsertionNode(myTreePointer, blockId, level2dimension1D, insertTreePointer);
248 }
249}
250
251
256template <MInt nDim>
258 if(treePointer->level < (nDim - 1)) {
259 treePointer->childs = new treeNode*[treePointer->noChilds];
260 for(MInt i = 0; i < treePointer->noChilds; i++) {
261 treePointer->childs[i] = new treeNode;
262 }
263 for(MInt countChild = 0; countChild < treePointer->noChilds; countChild++) {
264 treePointer->childs[countChild]->level = (treePointer->level) + 1;
265 treeNode* myTreePointer = treePointer->childs[countChild];
266 initializeChilds(myTreePointer);
267 }
268 }
269}
270
275template <MInt nDim>
277 if(treePointer->level < (nDim - 1)) {
278 for(MInt i = 0; i < treePointer->noChilds; i++) {
279 treeNode* myTreePointer = treePointer->childs[i];
280 destroyChilds(myTreePointer);
281 }
282 delete[] treePointer->childs;
283 }
284}
285
291template <MInt nDim>
293 treeNode* myTreePointer{};
294 MInt sumLeaf{};
295
296 if(treePointer->level < (nDim - 1)) {
297 sumLeaf = 0;
298 for(MInt i = 0; i < treePointer->noChilds; i++) {
299 myTreePointer = treePointer->childs[i];
300 sumLeaf += sumLeaves(myTreePointer);
301 }
302 } else {
303 sumLeaf = treePointer->noChilds;
304 }
305 treePointer->noLeaves = sumLeaf;
306 return sumLeaf;
307}
308
315template <MInt nDim>
317 const MInt blockId,
318 const MInt (&level2dimension1D)[nDim]) {
319 destroyChilds(insertTreePointer);
320 insertTreePointer->noChilds++;
321 MInt myNoLeaves = insertTreePointer->noLeaves + 1;
322 initializeChilds(insertTreePointer);
323
324 if(insertTreePointer->level < (nDim - 1)) {
325 for(MInt i = insertTreePointer->noChilds + 1; i <= myNoLeaves; i++) {
326 treeNode* myInsertTreePointer = insertTreePointer;
327 traverseForInsertionNode(insertTreePointer, blockId, level2dimension1D, myInsertTreePointer);
328 insertChildAtNode(myInsertTreePointer, blockId, level2dimension1D);
329 insertTreePointer->noLeaves = sumLeaves(insertTreePointer);
330 }
331 }
332}
333
339template <MInt nDim>
341 // select the block id which currently contains the largest partition
342 MInt blockId = 0;
343 MLongFloat maxWeightedSize = 0;
344 for(MInt i = 0; i < m_noPartitions; i++) {
345 const MLongFloat weightedSize = m_partitionInfo[i]->totalSize * m_partitionInfo[i]->weight;
346 if(weightedSize > maxWeightedSize) {
347 blockId = m_partitionInfo[i]->blockId;
348 maxWeightedSize = weightedSize;
349 }
350 }
351
352 // set this tree root for further processing
353 treeNode* myTreeRootPointer = treeRoot[blockId];
354 treeNode* myInsertTreePointer = treeRoot[blockId];
355
356 // pick corresponding level->dimension map
357 MInt level2dimension1D[nDim];
358 for(MInt i = 0; i < nDim; i++) {
359 level2dimension1D[i] = level2dimension(i, blockId);
360 }
361
362 // now find insertion point
363 traverseForInsertionNode(myTreeRootPointer, blockId, level2dimension1D, myInsertTreePointer);
364 insertChildAtNode(myInsertTreePointer, blockId, level2dimension1D);
365 myTreeRootPointer->noLeaves = sumLeaves(myTreeRootPointer);
366
367 // update partition information
368 m_noPartitions++;
369 setPartitionInfo(treeRoot, level2dimension);
370}
371
375template <MInt nDim>
377 treeNode* myTreeRootPointer;
378 treeNode** myTreeRoot;
379 MInt countBlock;
380 m_noPartitions = m_noBlocks;
381 myTreeRoot = new treeNode*[m_noBlocks];
382
383 // Create one tree root for each block
384 for(MInt i = 0; i < m_noBlocks; i++) {
385 myTreeRoot[i] = new treeNode;
386 }
387
388 for(countBlock = 0; countBlock < m_noBlocks; countBlock++) {
389 myTreeRootPointer = myTreeRoot[countBlock];
390 initializeChilds(myTreeRootPointer);
391 }
392
393 // prepare by filling level2dimension
394 MIntScratchSpace level2dimensionA(nDim, m_noBlocks, AT_, "level2dimensionA");
395 level2dimensionA.fill(0);
396
397 for(MInt k = 0; k < m_noBlocks; k++) {
398 std::vector<std::pair<MInt, MInt>> dummy{};
399 for(MInt i = 0; i < nDim; i++) {
400 dummy.push_back({m_blockInfo[k]->size[i], i});
401 }
402 // sort pairs by first value in descending order
403 std::sort(dummy.begin(), dummy.end(), [](auto& left, auto& right) { return left.first > right.first; });
404 for(MInt i = 0; i < nDim; i++) {
405 level2dimensionA(i, k) = dummy[i].second;
406 }
407 }
408
409 setPartitionInfo(myTreeRoot, level2dimensionA);
410
411 // now add leaves successively
412 while(m_noPartitions < m_noDomains) {
413 addLeaf(myTreeRoot, level2dimensionA);
414 }
415
416 // Clean up tree
417 for(MInt i = 0; i < m_noBlocks; i++) {
418 myTreeRootPointer = myTreeRoot[i];
419 destroyChilds(myTreeRootPointer);
420 }
421
422 for(MInt j = 0; j < m_noBlocks; j++) {
423 delete myTreeRoot[j];
424 }
425
426 delete[] myTreeRoot;
427}
428
433template <MInt nDim>
435 // temporary scratch to send/receive the partition Info
436 MIntScratchSpace filePartitionInfo((3 + 3 + 4) * m_noDomains, "filePartitionInfo", AT_);
437 filePartitionInfo.fill(0);
438 MInt noCPUs = 0;
439
440 if(globalDomainId() == 0) { // only root reads in the information
441 ParallelIoHdf5 pio("partitionFile.hdf5", maia::parallel_io::PIO_APPEND, MPI_COMM_SELF);
442 pio.getAttribute<MInt>(&noCPUs, "noCPUs", "");
443 }
444
445 // broadcast the number of cpus
446 MPI_Bcast(&noCPUs, 1, MPI_INT, 0, m_mpiComm, AT_, "noCPUs");
447 if(noCPUs != m_noDomains) {
448 m_log << "WARNING no CPUs is not equal to ones given in partitionFile" << endl;
449 return false;
450 }
451
452 if(globalDomainId() == 0) {
453 // read the information out of the file
454 ParallelIoHdf5 pio("partitionFile.hdf5", maia::parallel_io::PIO_APPEND, MPI_COMM_SELF);
455 ParallelIo::size_type asize[2] = {noCPUs, 10};
456 pio.readArray(filePartitionInfo.getPointer(), "", "partitionData", 2, asize);
457 }
458
459 MPI_Bcast(filePartitionInfo.getPointer(), noCPUs * (10), MPI_INT, 0, m_mpiComm, AT_,
460 "filePartitionInfo.getPointer()");
461
462 m_noPartitions = noCPUs;
463 for(MInt j = 0; j < m_noPartitions; j++) {
464 m_partitionInfo.push_back(make_unique<PartitionInfo<nDim>>());
465 }
466
467 // put the information into the right place
468 for(MInt j = 0; j < m_noPartitions; j++) {
469 for(MInt dim = 0; dim < nDim; dim++) {
470 m_partitionInfo[j]->size[dim] = filePartitionInfo[j * 10 + dim];
471 m_partitionInfo[j]->offset[dim] = filePartitionInfo[j * 10 + 3 + dim];
472 }
473
474 m_partitionInfo[j]->blockId = filePartitionInfo[j * 10 + 6];
475 m_partitionInfo[j]->partitionId = filePartitionInfo[j * 10 + 7];
476 m_partitionInfo[j]->cpu = filePartitionInfo[j * 10 + 8];
477 m_partitionInfo[j]->totalSize = filePartitionInfo[j * 10 + 9];
478 }
479
480 return true;
481}
482
483template <MInt nDim>
485 MFloat a;
486 ((x - floor(x)) >= 0.5) ? a = ceil(x) : a = floor(x);
487 return (MInt)a;
488}
489
490// Explicit instantiations for 2D and 3D
491template class StructuredDecomposition<2>;
492template class StructuredDecomposition<3>;
void getAttribute(T *value, const MString &name, const MString &datasetName="")
Retrieve a file or dataset attribute.
Definition: parallelio.h:1788
MBool hasAttribute(const MString &name, const MString &path="")
Check if a given attribute exists in the file.
Definition: parallelio.h:714
size_type getArraySize(const MString &name, const size_type dimensionId=0)
Get the length of an array in the file.
Definition: parallelio.h:667
void readArray(T *array, const MString &name, size_type memoryStride=-1, size_type diskStride=-1)
Read array data from file. [MPI]
Definition: parallelio.h:1523
This class is a ScratchSpace.
Definition: scratch.h:758
T * getPointer() const
Deprecated: use begin() instead!
Definition: scratch.h:316
void fill(T val)
fill the scratch with a given value
Definition: scratch.h:311
Class for the decomposition (partition) of structured grids.
MBool readFromFile()
Read a precomputed partition info from an external file.
MInt sumLeaves(treeNode *&treePointer)
Compute the number of leaves of a tree nodes.
void addLeaf(treeNode **&treeRoot, MIntScratchSpace &level2dimension)
Add a leaf into the tree.
void decompose()
Decompose the grid into partitions for the set number of domains.
std::vector< std::unique_ptr< PartitionInfo< nDim > > > m_blockInfo
StructuredDecomposition(const MInt, const MString, MPI_Comm, const MInt)
Constructor for HDF5 grid file name Reads the size of the coordinate datasets for all blocks.
void destroyChilds(treeNode *&treePointer)
Recursively destroy all childs of a tree node.
void setPartitionInfo(treeNode **&, MIntScratchSpace &)
Computes partition information from given BCT.
void insertChildAtNode(treeNode *&, const MInt, const MInt(&)[nDim])
Inserts new child node below given tree node.
void traverseForInsertionNode(treeNode *&treePointer, const MInt, const MInt(&)[nDim], treeNode *&)
Recursively traverses the given tree to find the next insertion position.
void initializeChilds(treeNode *&treePointer)
Initializes the childs of the given tree node.
void setPartitionInfoHelper(treeNode *&, const MInt, const MInt(&)[nDim], MInt &, MInt(&)[nDim], MInt(&)[nDim], MInt(&)[nDim])
Helper function for setPartionInfo.
treeNode ** childs
MBool approx(const T &, const U &, const T)
Definition: functions.h:272
MInt globalDomainId()
Return global domain id.
InfoOutFile m_log
int32_t MInt
Definition: maiatypes.h:62
std::basic_string< char > MString
Definition: maiatypes.h:55
double MFloat
Definition: maiatypes.h:52
long double MLongFloat
Definition: maiatypes.h:53
bool MBool
Definition: maiatypes.h:58
int MPI_Bcast(void *buffer, int count, MPI_Datatype datatype, int root, MPI_Comm comm, const MString &name, const MString &varname)
same as MPI_Bcast
const MInt PIO_APPEND
Definition: parallelio.h:38
Definition: contexttypes.h:19