This method is structured in 4 parts: Step 1: Identify cells in slice and sort them by their sliceHilbertId. First the data is saved a collector and after identification sorted. Then the data is distrubted into single arrays. Step 2: Create MPI communicator only for domains with slice cells Step 3: Gather all data for slice grid file. First the hilbertId information is exchanged globally. Then the global offsets for the local data can be determined on every domain. It is important to understand, that the number of hilbertIds is different on every domain and that they can be distributed in in any order, since a slice of a 3D hilbert curve is performed. After that the neighbor domains are determined to exchange cellIds to determine the cell mapping to slice globalId completly. At the end the data arrays are determined. Step 4: Write grid file with slice data. Since the data is sorted by slice HilbertId, the data writing is done in data chunks with the same hilbertId. If a domain has less hilbertIds than others, then dummy calls a performed.
11415 {
11416 TRACE();
11417
11418
11421
11422 std::cerr << "Warning: CartesianGrid::createGridSlice minLevel and maxLevel grid are equal! "
11423 << "(hang-up might occur)" << std::endl;
11424 }
11425 }
11426
11427
11428 IF_CONSTEXPR(nDim != 3) { TERMM(1, "Can not create a 2D slice from a 2D grid."); }
11429
11430 m_log <<
"Creating 2D slice from grid... ";
11431
11432
11433 const MInt nDimSlice = 2;
11434
11435
11437
11438 array<MInt, nDimSlice> coordArray{};
11439 if(axis == "x") {
11440 axisNum = 0;
11441 coordArray = {{1, 2}};
11442 } else if(axis == "y") {
11443 axisNum = 1;
11444 coordArray = {{0, 2}};
11445 } else if(axis == "z") {
11446 axisNum = 2;
11447 coordArray = {{0, 1}};
11448 } else {
11449 TERMM(1, "Unknown axis! Check your property file.");
11450 }
11451
11452 MBool isAtCenterOfGravity =
false;
11453
11455 isAtCenterOfGravity = true;
11456 }
11457
11459
11461
11462
11466
11467 std::array<MFloat, nDimSlice> targetGridCenterOfGravity{};
11470
11471
11472 MInt noSlicePartitionCells = 0;
11473 MInt noSliceMinLevelCells = 0;
11474 MInt noSliceCells = 0;
11475 MLong noLeafCells = 0;
11476
11477
11478
11481 sliceCellCollector[i].fill(0);
11482 }
11483
11485
11486
11487 MIntScratchSpace localPartitionCellGlobalIds(noLocalPartitionCells + 1, AT_,
"localPartitionCellGlobalIds");
11489
11490
11492
11493
11494 std::map<MInt, MInt> hilbertIdCount;
11495 std::map<MInt, MInt> hilbertIdCountSolver;
11496
11497 for(
MInt i = 0, cellId = 0; i < noLocalPartitionCells; i++) {
11498
11500 MInt isPartitionCell = 0;
11501
11503 if(isAtCenterOfGravity) {
11504 eps = MFloatEps;
11505 }
11507 isPartitionCell = 1;
11508 noSlicePartitionCells++;
11509 } else {
11510
11512 }
11513
11514
11516 while(cellId < nextPartitionCell) {
11517
11520 continue;
11521 }
11522
11523
11524 if(cellId == partitionCellId) {
11525 sliceCellCollector[noSliceCells][2] = isPartitionCell;
11526 isPartitionCell = 0;
11527 }
11528
11529
11531
11532 array<MFloat, 2> normedSliceCoord{};
11533 normedSliceCoord[0] =
11536 normedSliceCoord[1] =
11539 const MLong sliceHilbertId =
11541
11542 hilbertIdCount[sliceHilbertId]++;
11543 sliceCellCollector[noSliceCells][1] = sliceHilbertId;
11544
11545
11546 if(solverId > -1 &&
a_solver(cellId, solverId)) hilbertIdCountSolver[sliceHilbertId]++;
11547
11548
11550 sliceCellCollector[noSliceCells][3] = 1;
11551 noSliceMinLevelCells++;
11552 }
11553
11554
11556 noLeafCells++;
11557 }
11558
11559 sliceCellCollector[noSliceCells][0] =
cellId;
11560 noSliceCells++;
11561
11562
11564 }
11565 }
11566
11568
11570
11571
11572 MIntScratchSpace slicePartitionCells(noSlicePartitionCells + 1, AT_,
"slicePartitionCells");
11573
11574 MIntScratchSpace sliceMinLevelCells(noSliceMinLevelCells, AT_,
"sliceMinLevelCells");
11575
11577
11579
11580 std::map<MInt, MInt> hilbertIdMinLevelCellsCount;
11581 std::map<MInt, MInt> hilbertIdPartitionCellsCount;
11582
11583 if(noSliceCells > 0) {
11584
11585 MInt noHilbertIds = 0;
11586 for(auto const& ent1 : hilbertIdCount) {
11587 hIdKeys[noHilbertIds] = ent1.first;
11588 noHilbertIds++;
11589 }
11590
11591
11592 std::stable_sort(sliceCellCollector.begin(), sliceCellCollector.begin() + noSliceCells,
11593 [](
const array<MInt, 4>&
a,
const array<MInt, 4>&
b) { return a[1] < b[1]; });
11594
11595
11596 for(
MInt i = 0, j = 0, k = 0; i < noSliceCells; i++) {
11597 cellIdsInSlice[i] = sliceCellCollector[i][0];
11598
11599 if(sliceCellCollector[i][2] == 1) {
11600 slicePartitionCells[j] = sliceCellCollector[i][0];
11601 hilbertIdPartitionCellsCount[sliceCellCollector[i][1]]++;
11602 j++;
11603 }
11604
11605 if(sliceCellCollector[i][3] == 1) {
11606 sliceMinLevelCells[k] = sliceCellCollector[i][0];
11607 hilbertIdMinLevelCellsCount[sliceCellCollector[i][1]]++;
11608 k++;
11609 }
11610 }
11611
11612 slicePartitionCells[noSlicePartitionCells] = cellIdsInSlice[noSliceCells - 1] + 1;
11613 }
11614
11615 if(solverId < 0) {
11616
11617 if(noSliceCellIds != nullptr) {
11618 *noSliceCellIds = noSliceCells;
11619 }
11620 if(sliceCellIds != nullptr && noSliceCells > 0) {
11621 std::copy_n(&cellIdsInSlice[0], noSliceCells, sliceCellIds);
11622 }
11623 } else {
11624
11625 MIntScratchSpace cellIdsInSliceSolver(noSliceCells, AT_,
"cellIdsInSliceSolver");
11626 MInt noSliceCellsSolver = 0;
11627 for(
MInt i = 0; i < noSliceCells; i++) {
11628 if(
a_solver(cellIdsInSlice[i], solverId)) {
11629 cellIdsInSliceSolver[noSliceCellsSolver] = cellIdsInSlice[i];
11630 noSliceCellsSolver++;
11631 }
11632 }
11633
11634
11635 if(noSliceCellIds != nullptr) {
11636 *noSliceCellIds = noSliceCellsSolver;
11637 }
11638 if(sliceCellIds != nullptr && noSliceCellsSolver > 0) {
11639 std::copy_n(&cellIdsInSliceSolver[0], noSliceCellsSolver, sliceCellIds);
11640 }
11641 }
11642
11644
11646
11647
11648 MPI_Comm mpiCommSlice = MPI_COMM_NULL;
11650
11651
11653 if(noSliceCells > 0) {
11654 MPI_Comm_rank(
mpiComm(), &rank);
11655 }
11656
11658 "rank", "sliceRanks[0]");
11659
11660
11661 const MInt noRelDomains = count_if(sliceRanks.begin(), sliceRanks.end(), [](
const MInt a) { return a != -1; });
11663
11664 map<MInt, MInt> domainMap;
11666 for(auto&& slice : sliceRanks) {
11667 if(slice != -1) {
11668 relDomains[position] = slice;
11669 domainMap[slice] = position;
11670 position++;
11671 }
11672 }
11673
11674 MPI_Group globalGroup, localGroup;
11676 MPI_Group_incl(globalGroup, noRelDomains, &relDomains[0], &localGroup, AT_);
11677
11678
11680
11683
11684
11685 if(noSliceCells < 1) {
11686 return;
11687 }
11688
11689
11691
11693
11694
11695
11696
11697
11698
11699 MInt noSliceDomains = -1;
11700 MPI_Comm_size(mpiCommSlice, &noSliceDomains);
11701
11702 MInt sliceDomain = -1;
11703 MPI_Comm_rank(mpiCommSlice, &sliceDomain);
11704
11705
11706 MIntScratchSpace noLocalHilbertIds(noSliceDomains, AT_,
"noLocalHilbertIds");
11707 noLocalHilbertIds[sliceDomain] = hilbertIdCount.size();
11708
11709
11710 MPI_Allgather(MPI_IN_PLACE, 1, MPI_INT, &noLocalHilbertIds[0], 1, MPI_INT, mpiCommSlice, AT_,
"MPI_IN_PLACE",
11711 "noLocalHilbertIds[0]");
11712
11713
11715 offsetsRecvData[0] = 0;
11716 for(
MInt i = 1; i < noSliceDomains; i++) {
11717 offsetsRecvData[i] = offsetsRecvData[i - 1] + noLocalHilbertIds[i - 1] * 6;
11718 }
11719 MInt noTotalRecvData = offsetsRecvData[noSliceDomains - 1] + noLocalHilbertIds[noSliceDomains - 1] * 6;
11720
11721
11722
11723
11724 MIntScratchSpace sendHilbertData(noLocalHilbertIds[sliceDomain] * 6, AT_,
"sendHilbertData");
11725 for(
MUlong i = 0; i < hIdKeys.size(); i++) {
11726 sendHilbertData[i * 6] = hIdKeys[i];
11727 sendHilbertData[i * 6 + 1] =
domainId();
11728 sendHilbertData[i * 6 + 2] = hilbertIdCount[hIdKeys[i]];
11729 sendHilbertData[i * 6 + 3] = hilbertIdPartitionCellsCount[hIdKeys[i]];
11730 sendHilbertData[i * 6 + 4] = hilbertIdMinLevelCellsCount[hIdKeys[i]];
11731 sendHilbertData[i * 6 + 5] = hilbertIdCountSolver[hIdKeys[i]];
11732 }
11733
11734
11735 MIntScratchSpace recvHilbertData(noTotalRecvData, AT_,
"recvHilbertData");
11736 MInt noSendData = sendHilbertData.size();
11738 for(
MInt i = 0; i < noSliceDomains; i++) {
11739 noRecvData[i] = noLocalHilbertIds[i] * 6;
11740 }
11741 MPI_Allgatherv(&sendHilbertData[0], noSendData, MPI_INT, &recvHilbertData[0], &noRecvData[0], &offsetsRecvData[0],
11742 MPI_INT, mpiCommSlice, AT_, "sendHilbertData[0]", "recvHilbertData[0]");
11743
11744
11746
11747 for(
MInt i = 0; i < noTotalRecvData / 6; i++) {
11748 sliceGlobalHilbertInfo[i][0] = recvHilbertData[i * 6];
11749 sliceGlobalHilbertInfo[i][1] = recvHilbertData[i * 6 + 1];
11750 sliceGlobalHilbertInfo[i][2] = recvHilbertData[i * 6 + 2];
11751 sliceGlobalHilbertInfo[i][3] = recvHilbertData[i * 6 + 3];
11752 sliceGlobalHilbertInfo[i][4] = recvHilbertData[i * 6 + 4];
11753 sliceGlobalHilbertInfo[i][5] = recvHilbertData[i * 6 + 5];
11754 }
11755
11756
11757 std::stable_sort(sliceGlobalHilbertInfo.begin(), sliceGlobalHilbertInfo.end(),
11758 [](
const array<MInt, 6>&
a,
const array<MInt, 6>&
b) { return a[1] < b[1]; });
11759
11760 std::stable_sort(sliceGlobalHilbertInfo.begin(), sliceGlobalHilbertInfo.end(),
11761 [](
const array<MInt, 6>&
a,
const array<MInt, 6>&
b) { return a[0] < b[0]; });
11762
11763
11764
11765
11766
11767 MInt maxNoHilbertIds = *std::max_element(noLocalHilbertIds.begin(), noLocalHilbertIds.end());
11768 MIntScratchSpace hilbertDomainOffset(maxNoHilbertIds, AT_,
"hilbertDomainOffset");
11769 hilbertDomainOffset.fill(0);
11770 MIntScratchSpace hilbertPartitionDomainOffset(maxNoHilbertIds, AT_,
"hilbertDomainOffset");
11771 hilbertPartitionDomainOffset.fill(0);
11772 MIntScratchSpace hilbertMinLevelDomainOffset(maxNoHilbertIds, AT_,
"hilbertDomainOffset");
11773 hilbertMinLevelDomainOffset.fill(0);
11774
11775
11776 MIntScratchSpace hilbertDomainOffsetSolver(maxNoHilbertIds, AT_,
"hilbertDomainOffsetSolver");
11777 hilbertDomainOffsetSolver.fill(0);
11778
11779
11780
11781
11782
11783
11784
11785
11786
11787
11788
11789
11790
11791
11792
11793 for(
MInt i = 0; i < noLocalHilbertIds[sliceDomain]; i++) {
11794 MInt offset = 0, offsetPart = 0, offsetMin = 0, offsetSolver = 0, j = 0;
11795
11796 while(sliceGlobalHilbertInfo[j][0] < hIdKeys[i]) {
11797 offset += sliceGlobalHilbertInfo[j][2];
11798 offsetPart += sliceGlobalHilbertInfo[j][3];
11799 offsetMin += sliceGlobalHilbertInfo[j][4];
11800 offsetSolver += sliceGlobalHilbertInfo[j][5];
11801 j++;
11802 }
11803
11804 while(sliceGlobalHilbertInfo[j][1] <
domainId()) {
11805 offset += sliceGlobalHilbertInfo[j][2];
11806 offsetPart += sliceGlobalHilbertInfo[j][3];
11807 offsetMin += sliceGlobalHilbertInfo[j][4];
11808 offsetSolver += sliceGlobalHilbertInfo[j][5];
11809 j++;
11810 }
11811 hilbertDomainOffset[i] = offset;
11812 hilbertPartitionDomainOffset[i] = offsetPart;
11813 hilbertMinLevelDomainOffset[i] = offsetMin;
11814
11815 if(solverId > -1) {
11816 hilbertDomainOffsetSolver[i] = offsetSolver;
11818 }
11819 }
11820
11821
11823
11825 mpiCommSlice, AT_, "noSliceCells", "domainOffset[0]");
11826
11827
11828 for(
MInt i = 0, offset = 0, tmp = 0; i < (noSliceDomains + 1); i++) {
11831 offset += tmp;
11832 }
11833
11834 MInt noLocalHilbertIdsSolver = 0;
11835
11836 if(solverId < 0) {
11837
11838 if(noSliceHilbertIds != nullptr) {
11839 *noSliceHilbertIds = noLocalHilbertIds[sliceDomain];
11840 }
11841
11842 if(sliceHilbertInfo != nullptr) {
11843 MIntScratchSpace hilbertInfo(noLocalHilbertIds[sliceDomain] * 3, AT_,
"hilbertInfo");
11844 for(
MInt i = 0; i < noLocalHilbertIds[sliceDomain]; i++) {
11845 hilbertInfo[i * 3] = hIdKeys[i];
11846 hilbertInfo[i * 3 + 1] = hilbertIdCount[hIdKeys[i]];
11847 hilbertInfo[i * 3 + 2] = hilbertDomainOffset[i];
11848 }
11849 std::copy_n(&hilbertInfo[0], noLocalHilbertIds[sliceDomain] * 3, sliceHilbertInfo);
11850 }
11851 } else {
11852 TERMM_IF_COND(noSliceHilbertIds == nullptr || sliceHilbertInfo == nullptr,
11853 "Error: solverId given, pointers need to be != nullptr.");
11854
11855 MIntScratchSpace hilbertInfoSolver(noLocalHilbertIds[sliceDomain] * 3, AT_,
"hilbertInfo");
11856 MInt noHIdsSolver = 0;
11857 for(
MInt i = 0; i < noLocalHilbertIds[sliceDomain]; i++) {
11858 if(hilbertIdCountSolver[hIdKeys[i]] > 0) {
11859 hilbertInfoSolver[noHIdsSolver * 3] = hIdKeys[i];
11860 hilbertInfoSolver[noHIdsSolver * 3 + 1] = hilbertIdCountSolver[hIdKeys[i]];
11861 hilbertInfoSolver[noHIdsSolver * 3 + 2] = hilbertDomainOffsetSolver[i];
11862 noHIdsSolver++;
11863 }
11864 }
11865
11866 noLocalHilbertIdsSolver = noHIdsSolver;
11867
11868 *noSliceHilbertIds = noHIdsSolver;
11869
11870 std::copy_n(&hilbertInfoSolver[0], noHIdsSolver * 3, sliceHilbertInfo);
11871 }
11872
11873
11875
11877
11878
11879
11880
11881
11882 array<MInt, nDimSlice * 2> nghbrArray{};
11883 if(axis == "x") {
11884 nghbrArray = {{2, 3, 4, 5}};
11885 } else if(axis == "y") {
11886 nghbrArray = {{0, 1, 4, 5}};
11887 } else if(axis == "z") {
11888 nghbrArray = {{0, 1, 2, 3}};
11889 }
11890
11891
11892 std::set<MInt> relatedDomains;
11893
11894
11895 for(
MInt i = 0; i < noSliceMinLevelCells; i++) {
11896
11897 for(
MInt j = 0; j < (nDimSlice * 2); j++) {
11898
11900 continue;
11901 }
11902
11906 relatedDomains.insert(nghbrDomainId);
11907 }
11908 }
11909 }
11910
11911
11912 MInt partitionLevelShift = 0;
11913 for(
MInt i = 0; i < noSlicePartitionCells; ++i) {
11915 partitionLevelShift =
mMax(levelDiff, partitionLevelShift);
11916 }
11917
11918
11919 set<MLong> partitionLevelAncestorIds;
11920 for(
MInt i = 0; i < noSliceCells; i++) {
11923 partitionLevelAncestorIds.insert(cellId);
11924 }
11925 }
11926
11927
11929
11931
11932
11933 map<MInt, MInt> idMap;
11934
11935 idMap[-1] = -1;
11936
11937 for(
MInt i = 0, j = 0, h = 0; i < noSliceCells; i++, j++) {
11938
11939 if(j == hilbertIdCount[hIdKeys[h]]) {
11940 h++;
11941 j = 0;
11942 }
11944 }
11945
11946
11947 if(!relatedDomains.empty()) {
11948
11949 std::vector<MInt> sliceRelatedDomains(relatedDomains.begin(), relatedDomains.end());
11950
11951
11952
11953 MInt noRecvCells = 0;
11954 MIntScratchSpace relDomainOffsets(sliceRelatedDomains.size() + 1, AT_,
"relDomainOffsets");
11955 for(
MUlong i = 0; i < sliceRelatedDomains.size(); i++) {
11956 relDomainOffsets[i] = noRecvCells;
11957 noRecvCells +=
11959 }
11960
11961 relDomainOffsets[sliceRelatedDomains.size()] = noRecvCells;
11962
11963
11965 recvIds.fill(-1);
11966
11967
11969 fill(sendRequests.begin(), sendRequests.end(), MPI_REQUEST_NULL);
11971 fill(recvRequests.begin(), recvRequests.end(), MPI_REQUEST_NULL);
11972
11973
11974 for(
MUlong i = 0; i < sliceRelatedDomains.size(); i++) {
11975 MInt d = sliceRelatedDomains[i];
11978 "recvIds[relDomainOffsets[i]]");
11979 }
11980
11981
11982 for(
MUlong i = 0; i < sliceRelatedDomains.size(); i++) {
11984 domainMap[
domainId()], mpiCommSlice, &sendRequests[i], AT_,
"cellIdsInSlice[0]");
11985 }
11986
11987
11988 MPI_Waitall(sliceRelatedDomains.size(), &recvRequests[0], MPI_STATUSES_IGNORE, AT_);
11989
11990
11991 MPI_Waitall(sliceRelatedDomains.size(), &sendRequests[0], MPI_STATUSES_IGNORE, AT_);
11992
11993
11994 for(
MUlong j = 0; j < sliceRelatedDomains.size(); j++) {
11995 MInt d = sliceRelatedDomains[j];
11996 std::vector<MInt> hId(noLocalHilbertIds[domainMap[d]]);
11997 std::vector<MInt> hCount(noLocalHilbertIds[domainMap[d]]);
11998 std::vector<MInt> hOffsets(noLocalHilbertIds[domainMap[d]]);
11999
12000
12001 for(
MUlong k = 0, i = 0; i < hId.size(); k++) {
12002 if(sliceGlobalHilbertInfo[k][1] == d) {
12003 hId[i] = sliceGlobalHilbertInfo[k][0];
12004 hCount[i] = sliceGlobalHilbertInfo[k][2];
12005
12008 while(sliceGlobalHilbertInfo[m][0] < hId[i]) {
12009 offset += sliceGlobalHilbertInfo[m][2];
12010 m++;
12011 }
12012 while(sliceGlobalHilbertInfo[m][1] < d) {
12013 offset += sliceGlobalHilbertInfo[m][2];
12014 m++;
12015 }
12016 hOffsets[i] = offset;
12017 i++;
12018 }
12019 }
12020
12021
12022 for(
MInt i = 0, n = 0, h = 0; i < (relDomainOffsets[j + 1] - relDomainOffsets[j]); i++, n++) {
12023 if(n == hCount[h]) {
12024 h++;
12025 n = 0;
12026 }
12027 idMap[recvIds[i + relDomainOffsets[j]] +
m_domainOffsets[d]] = n + hOffsets[h];
12028 }
12029 }
12030 }
12031
12032
12033 MBool onlySliceCells =
false;
12034 if(noSlicePartitionCells < 1) {
12035 onlySliceCells = true;
12036 }
12037
12038
12039
12040 if(noSlicePartitionCells > 0) {
12041 std::map<MInt, MInt>::iterator it;
12042 MInt lastKeyId = hIdKeys.size() - 1;
12043 MInt val = hilbertIdCount[hIdKeys[lastKeyId]] + hilbertDomainOffset[lastKeyId];
12044
12046
12047 auto result = std::find_if(idMap.begin(), idMap.end(), [val](const auto& mo) { return mo.second == val; });
12048
12049 if(it == idMap.end() && result == idMap.end()) {
12051 hilbertIdCount[hIdKeys[lastKeyId]] + hilbertDomainOffset[lastKeyId];
12052 } else if(result != idMap.end()) {
12053 MInt foundKey = result->first;
12055 }
12056 }
12057
12059
12061
12062
12063 const MUint noChildInRef = pow(2, nDimSlice);
12064
12065
12066 MIntScratchSpace partitionCellsId(noSlicePartitionCells, AT_,
"partitionCellsId");
12067 MFloatScratchSpace partitionCellsWorkload(noSlicePartitionCells, AT_,
"partitionCellsWorkLoad");
12068 MInt maxNoOffsprings = 0;
12069 MFloat maxWorkload = 0.0;
12070 MFloat totalWorkload = 0.0;
12071
12072 for(
MInt i = 0, c = 0, partitionCellNoOffsprings = 0; i < noSlicePartitionCells; i++) {
12074
12075 partitionCellNoOffsprings = 1;
12076
12077
12079
12080 if(
a_level(cellIdsInSlice[c]) >
a_level(slicePartitionCells[i])) {
12081 partitionCellNoOffsprings++;
12082 }
12083 c++;
12084 if(c == noSliceCells) {
12085 break;
12086 }
12087 }
12088
12089
12090 partitionCellsWorkload[i] = partitionCellNoOffsprings;
12091
12092 totalWorkload += partitionCellsWorkload[i];
12093
12094 maxNoOffsprings =
mMax(partitionCellNoOffsprings, maxNoOffsprings);
12095
12096 maxWorkload =
mMax(partitionCellsWorkload[i], maxWorkload);
12097 }
12098
12099
12100 MLong totalNoPartitionLevelAncestors = 0;
12101 MLong localPartitionLevelAncestorCount = (signed)partitionLevelAncestorIds.size();
12102 MPI_Allreduce(&localPartitionLevelAncestorCount, &totalNoPartitionLevelAncestors, 1, MPI_LONG, MPI_SUM, mpiCommSlice,
12103 AT_, "localPartitionLevelAncestorCount", "totalNoPartitionLevelAncestors");
12104
12107 "partitionLevelShift", "maxPartitionLevelShift");
12108
12109
12110 MPI_Allreduce(MPI_IN_PLACE, &maxNoOffsprings, 1, MPI_INT, MPI_MAX, mpiCommSlice, AT_,
"MPI_IN_PLACE",
12111 "maxNoOffsprings");
12112
12113 MPI_Allreduce(MPI_IN_PLACE, &totalWorkload, 1, MPI_DOUBLE, MPI_SUM, mpiCommSlice, AT_,
"MPI_IN_PLACE",
12114 "totalWorkload");
12115 MPI_Allreduce(MPI_IN_PLACE, &maxWorkload, 1, MPI_DOUBLE, MPI_MAX, mpiCommSlice, AT_,
"MPI_IN_PLACE",
"maxWorkload");
12116 MFloat maxNoCPUs = totalWorkload / maxWorkload;
12117
12118
12120
12122
12123
12124 array<MInt, 8> childArray{};
12125 if(axis == "x") {
12126 childArray = {{0, 0, 1, 1, 2, 2, 3, 3}};
12127 } else if(axis == "y") {
12128 childArray = {{0, 1, 0, 1, 2, 3, 2, 3}};
12129 } else if(axis == "z") {
12130 childArray = {{0, 1, 2, 3, 0, 1, 2, 3}};
12131 }
12132
12133
12137 const array<MFloat, nDimSlice* 2> targetGridBoundingBox = {
12140 array<MFloat, nDimSlice * 2> geometryExtents{};
12141 MInt decisiveDirection = 0;
12142 for(
MInt dir = 0; dir < nDimSlice; dir++) {
12143 geometryExtents[dir] = targetGridBoundingBox[dir + nDimSlice] - targetGridBoundingBox[dir];
12144 decisiveDirection = geometryExtents[dir] > geometryExtents[decisiveDirection] ? dir : decisiveDirection;
12145 }
12146
12147
12148 MLongScratchSpace sliceMinLevelCellsTreeId(noSliceMinLevelCells, AT_,
"sliceMinLevelCellsTreeId");
12149 for(
MInt i = 0; i < noSliceMinLevelCells; ++i) {
12150 std::array<MFloat, nDimSlice> coord = {
12152 maia::grid::hilbert::coordinatesToTreeId<nDimSlice>(sliceMinLevelCellsTreeId[i], &coord[0],
12155 }
12156
12157
12158 MLongScratchSpace sliceMinLevelCellsNghbrIds(noSliceMinLevelCells, 2 * nDimSlice, AT_,
"sliceMinLevelCellsNghbrIds");
12159 for(
MInt i = 0; i < noSliceMinLevelCells; ++i) {
12160 for(
MInt j = 0; j < (nDimSlice * 2); j++) {
12162 if(nghbrId == -1) {
12163 sliceMinLevelCellsNghbrIds(i, j) = -1;
12164 } else {
12165 sliceMinLevelCellsNghbrIds(i, j) = idMap[
a_globalId(nghbrId)];
12166 }
12167 }
12168 }
12169
12170
12172
12174
12177
12178 for(
MInt i = 0; i < noSliceCells; ++i) {
12180 MUint noChilds = 0;
12181 MUint curChildPos = 0;
12182
12183 for(
MInt j = 0; j < pow(2, nDim) && curChildPos < noChildInRef; j++) {
12184 if(
a_childId(cellIdsInSlice[i], j) != -1) {
12185
12188 continue;
12189 }
12190 noChilds += 1;
12191 curChildPos++;
12192 }
12193 }
12194
12200 pos = childArray[j];
12201 }
12202 }
12203 }
12204
12205 MUint isMinLvl = 0;
12207 isMinLvl = (
MUint)1;
12208 }
12209 MUint tmpBit = noChilds | (pos << 4) | (isMinLvl << 7);
12210 sliceCellInfo[i] =
static_cast<MUchar>(tmpBit);
12211
12212
12215 }
12216 }
12217
12218
12220
12222
12223
12225
12226
12227 const MString gridFileName = fileName + ParallelIo::fileExt();
12228 ParallelIo file(gridFileName, PIO_REPLACE, mpiCommSlice);
12229
12230
12231 ParallelIo::size_type sliceCellsOffset, noTotalCells;
12232 ParallelIo::calcOffset(noSliceCells, &sliceCellsOffset, &noTotalCells, mpiCommSlice);
12233
12234
12235 if(onlySliceCells) {
12236 noSlicePartitionCells = 0;
12237 }
12238
12240 ParallelIo::calcOffset(noSlicePartitionCells, &partitionOffset, &
noPartitionCells, mpiCommSlice);
12241
12242 ParallelIo::size_type minLevelCellsOffset, noTotalMinLevelCells;
12243 ParallelIo::calcOffset(noSliceMinLevelCells, &minLevelCellsOffset, &noTotalMinLevelCells, mpiCommSlice);
12244
12247
12248
12251 file.defineArray(PIO_LONG, "minLevelCellsTreeId", noTotalMinLevelCells);
12252 file.defineArray(PIO_LONG, "minLevelCellsNghbrIds", 2 * nDimSlice * noTotalMinLevelCells);
12253 file.defineArray(PIO_UCHAR, "cellInfo", noTotalCells);
12254
12258 if(writeSolver) {
12259 file.defineArray(PIO_UCHAR, "solver", noTotalCells);
12260
12261 for(
MInt i = 0; i < noSliceCells; ++i) {
12266 tmpBit |= (1 << solver);
12267 }
12268 }
12269 solverBits[i] =
static_cast<MUchar>(tmpBit);
12270 }
12271 }
12272
12273
12275 MPI_Allreduce(MPI_IN_PLACE, &
maxLevel, 1, MPI_INT, MPI_MAX, mpiCommSlice, AT_,
"MPI_IN_PLACE",
"maxLevel");
12276
12277 MPI_Allreduce(MPI_IN_PLACE, &noLeafCells, 1, MPI_LONG, MPI_SUM, mpiCommSlice, AT_,
"MPI_IN_PLACE",
"noLeafCells");
12278
12279
12282 MLong solverCount = 0;
12283 for(
MInt i = 0; i < noSliceCells; i++) {
12286 solverCount++;
12287 }
12288 }
12289 MPI_Allreduce(MPI_IN_PLACE, &solverCount, 1, MPI_LONG, MPI_SUM, mpiCommSlice, AT_,
"MPI_IN_PLACE",
"solverCount");
12290 file.setAttributes(&solverCount,
"noCells_" + std::to_string(
b), 1);
12291 }
12292 }
12293
12294
12295 file.setAttributes(&nDimSlice, "nDim", 1);
12296 file.setAttributes(&
noSolvers,
"noSolvers", 1);
12297 file.setAttributes(&tstep, "globalTimeStep", 1);
12298 file.setAttributes(&noTotalCells, "noCells", 1);
12299 file.setAttributes(&noLeafCells, "noLeafCells", 1);
12300 file.setAttributes(&noTotalMinLevelCells, "noMinLevelCells", 1);
12302 file.setAttributes(&totalNoPartitionLevelAncestors, "noPartitionLevelAncestors", 1);
12303 file.setAttributes(&
m_minLevel,
"minLevel", 1);
12304 file.setAttributes(&
maxLevel,
"maxLevel", 1);
12308 file.setAttributes(&
centerOfGravity[0],
"centerOfGravity", nDimSlice);
12309 file.setAttributes(&
boundingBox[0],
"boundingBox", 2 * nDimSlice);
12310
12311
12315 file.setAttributes(&targetGridCenterOfGravity[0], "multiSolverCenterOfGravity", nDimSlice);
12316 file.setAttributes(&targetGridBoundingBox[0], "multiSolverBoundingBox", 2 * nDimSlice);
12317 }
12318
12320 file.setAttributes(&decisiveDirection, "decisiveDirection", 1);
12321 file.setAttributes(&totalWorkload, "totalWorkload", 1);
12322 file.setAttributes(&maxWorkload, "partitionCellMaxWorkload", 1);
12323 file.setAttributes(&avgWorkload, "partitionCellAverageWorkload", 1);
12324 file.setAttributes(&maxNoOffsprings, "partitionCellMaxNoOffspring", 1);
12325 file.setAttributes(&avgOffspring, "partitionCellAverageNoOffspring", 1);
12328 file.setAttributes(&maxNoCPUs, "maxNoBalancedCPUs", 1);
12329
12330
12331 file.setAttribute(axis, "sliceAxis");
12332 file.setAttribute(intercept, "sliceIntercept");
12333
12334 MBool optimizedSliceIo =
true;
12335 optimizedSliceIo = Context::getBasicProperty<MBool>("optimizedSliceIo", AT_, &optimizedSliceIo);
12336
12337
12338 if(optimizedSliceIo) {
12339 std::vector<MInt> contHilbertIdsPartitionCounts{};
12340 std::vector<MInt> contHilbertIdsPartitionOffset{};
12341
12342 std::vector<MInt> contHilbertIdCount{};
12343 std::vector<MInt> contHilbertIdOffset{};
12344
12345 std::vector<MInt> contHilbertIdMinCellCount{};
12346 std::vector<MInt> contHilbertIdMinCellOffset{};
12347
12348 {
12349 MInt contHIdMinCellCount = -1;
12350 std::vector<MInt> minCellCount{};
12351 std::vector<MInt> minCellOffset{};
12352
12353 for(
MInt i = 0, localMinOffset = 0; i < noLocalHilbertIds[sliceDomain]; i++) {
12354 const MInt count = hilbertIdMinLevelCellsCount[hIdKeys[i]];
12355 const MInt offset = hilbertMinLevelDomainOffset[i];
12356 if(localMinOffset < noSliceMinLevelCells) {
12357 if(count > 0) {
12358 minCellCount.push_back(count);
12359 minCellOffset.push_back(offset);
12360 }
12361 }
12362 localMinOffset += count;
12363 }
12364
12365 contHIdMinCellCount = (noSliceMinLevelCells > 0) ? minCellCount[0] : -1;
12366 if(noSliceMinLevelCells == 1) {
12367 contHilbertIdMinCellCount.push_back(contHIdMinCellCount);
12368 contHilbertIdMinCellOffset.push_back(minCellOffset[0]);
12369 }
12370
12371
12372 for(
MInt h = 1, i = 0; h < noSliceMinLevelCells; h++) {
12373 if(minCellCount[h - 1] + minCellOffset[h - 1] == minCellOffset[h]) {
12374 contHIdMinCellCount += minCellCount[h];
12375 } else {
12376 contHilbertIdMinCellCount.push_back(contHIdMinCellCount);
12377 contHilbertIdMinCellOffset.push_back(minCellOffset[i]);
12378
12379 i = h;
12380 contHIdMinCellCount = minCellCount[h];
12381 }
12382
12383
12384 if(h == noSliceMinLevelCells - 1) {
12385 contHilbertIdMinCellCount.push_back(contHIdMinCellCount);
12386 contHilbertIdMinCellOffset.push_back(minCellOffset[i]);
12387 }
12388 }
12389 }
12390
12391 {
12392 MInt contHIdPartitionCellCount = hilbertIdPartitionCellsCount[hIdKeys[0]];
12393 MInt contHIdCount = hilbertIdCount[hIdKeys[0]];
12394
12395 if(noLocalHilbertIds[sliceDomain] == 1) {
12396 contHilbertIdsPartitionCounts.push_back(contHIdPartitionCellCount);
12397 contHilbertIdsPartitionOffset.push_back(hilbertPartitionDomainOffset[0]);
12398
12399 contHilbertIdCount.push_back(contHIdCount);
12400 contHilbertIdOffset.push_back(hilbertDomainOffset[0]);
12401 }
12402
12403
12404 for(
MInt h = 1, i = 0; h < noLocalHilbertIds[sliceDomain]; h++) {
12405
12406 if(hilbertIdPartitionCellsCount[hIdKeys[h - 1]] + hilbertPartitionDomainOffset[h - 1]
12407 == hilbertPartitionDomainOffset[h]) {
12408 TERMM_IF_NOT_COND(hilbertIdCount[hIdKeys[h - 1]] + hilbertDomainOffset[h - 1] == hilbertDomainOffset[h],
12409 "Error: cells not contiguous");
12410 contHIdPartitionCellCount += hilbertIdPartitionCellsCount[hIdKeys[h]];
12411 contHIdCount += hilbertIdCount[hIdKeys[h]];
12412 } else {
12413
12414 contHilbertIdsPartitionCounts.push_back(contHIdPartitionCellCount);
12415 contHilbertIdsPartitionOffset.push_back(hilbertPartitionDomainOffset[i]);
12416
12417 contHilbertIdCount.push_back(contHIdCount);
12418 contHilbertIdOffset.push_back(hilbertDomainOffset[i]);
12419
12420 i = h;
12421 contHIdPartitionCellCount = hilbertIdPartitionCellsCount[hIdKeys[h]];
12422 contHIdCount = hilbertIdCount[hIdKeys[h]];
12423 }
12424
12425
12426 if(h == noLocalHilbertIds[sliceDomain] - 1) {
12427 contHilbertIdsPartitionCounts.push_back(contHIdPartitionCellCount);
12428 contHilbertIdsPartitionOffset.push_back(hilbertPartitionDomainOffset[i]);
12429
12430 contHilbertIdCount.push_back(contHIdCount);
12431 contHilbertIdOffset.push_back(hilbertDomainOffset[i]);
12432 }
12433 }
12434 }
12435 const MInt noContHilbertIds = contHilbertIdsPartitionCounts.size();
12436
12437 if(solverId > -1 && noSliceContHilbertIds != nullptr) {
12438 TERMM_IF_COND(sliceContiguousHilbertInfo == nullptr,
12439 "Error: sliceContiguousHilbertInfo is a nullptr but noSliceContHilbertIds is not.");
12440
12441 std::vector<MInt> hIdCountSolver;
12442 std::vector<MInt> hDomainOffsetSolver;
12443 for(
MInt i = 0; i < noLocalHilbertIds[sliceDomain]; i++) {
12444 if(hilbertIdCountSolver[hIdKeys[i]] > 0) {
12445 hIdCountSolver.push_back(hilbertIdCountSolver[hIdKeys[i]]);
12446 hDomainOffsetSolver.push_back(hilbertDomainOffsetSolver[i]);
12447 }
12448 }
12449
12450 std::vector<MInt> contHilbertIdCountSolver{};
12451 std::vector<MInt> contHilbertIdOffsetSolver{};
12452
12453
12454 MInt contHIdCountSolver = (hIdCountSolver.size() > 0) ? hIdCountSolver[0] : 0;
12455
12456 if(noLocalHilbertIdsSolver == 1) {
12457 contHilbertIdCountSolver.push_back(contHIdCountSolver);
12458 contHilbertIdOffsetSolver.push_back(hDomainOffsetSolver[0]);
12459 }
12460
12461
12462 for(
MInt h = 1, i = 0; h < noLocalHilbertIdsSolver; h++) {
12463
12464 if(hIdCountSolver[h - 1] + hDomainOffsetSolver[h - 1] == hDomainOffsetSolver[h]) {
12465 contHIdCountSolver += hIdCountSolver[h];
12466 } else {
12467
12468 contHilbertIdCountSolver.push_back(contHIdCountSolver);
12469 contHilbertIdOffsetSolver.push_back(hDomainOffsetSolver[i]);
12470
12471 i = h;
12472 contHIdCountSolver = hIdCountSolver[h];
12473 }
12474
12475
12476 if(h == noLocalHilbertIdsSolver - 1) {
12477 contHilbertIdCountSolver.push_back(contHIdCountSolver);
12478 contHilbertIdOffsetSolver.push_back(hDomainOffsetSolver[i]);
12479 }
12480 }
12481
12482 const MInt noContHilbertIdsSolver = contHilbertIdCountSolver.size();
12483 *noSliceContHilbertIds = noContHilbertIdsSolver;
12484
12485 if(noContHilbertIdsSolver > 0) {
12486 MIntScratchSpace contHilbertInfo(noContHilbertIdsSolver * 3, AT_,
"contHilbertInfo");
12487 for(
MInt i = 0; i < noContHilbertIdsSolver; i++) {
12488 contHilbertInfo[i * 3] = -1;
12489 contHilbertInfo[i * 3 + 1] = contHilbertIdCountSolver[i];
12490 contHilbertInfo[i * 3 + 2] = contHilbertIdOffsetSolver[i];
12491 }
12492 std::copy_n(&contHilbertInfo[0], noContHilbertIdsSolver * 3, sliceContiguousHilbertInfo);
12493 }
12494 } else if(solverId < 0) {
12495
12496 if(noSliceContHilbertIds != nullptr) {
12497 *noSliceContHilbertIds = noContHilbertIds;
12498 }
12499
12500 if(sliceContiguousHilbertInfo != nullptr) {
12501 MIntScratchSpace contHilbertInfo(noContHilbertIds * 3, AT_,
"contHilbertInfo");
12502 for(
MInt i = 0; i < noContHilbertIds; i++) {
12503 contHilbertInfo[i * 3] = -1;
12504 contHilbertInfo[i * 3 + 1] = contHilbertIdCount[i];
12505 contHilbertInfo[i * 3 + 2] = contHilbertIdOffset[i];
12506 }
12507 std::copy_n(&contHilbertInfo[0], noContHilbertIds * 3, sliceContiguousHilbertInfo);
12508 }
12509 }
12510
12512 for(
MInt i = 0, localOffset = 0, localPartOffset = 0; i < hilbertDomainOffset.size0(); i++) {
12513
12514 if(i < noContHilbertIds) {
12515
12516 file.setOffset(contHilbertIdsPartitionCounts[i], contHilbertIdsPartitionOffset[i]);
12517 file.writeArray(&partitionCellsId[0] + localPartOffset, "partitionCellsGlobalId");
12518 file.writeArray(&partitionCellsWorkload[0] + localPartOffset, "partitionCellsWorkload");
12519 localPartOffset += contHilbertIdsPartitionCounts[i];
12520
12521
12522 file.setOffset(contHilbertIdCount[i], contHilbertIdOffset[i]);
12523 file.writeArray(&sliceCellInfo[0] + localOffset, "cellInfo");
12524
12525 if(writeSolver) {
12526 file.writeArray(&solverBits[0] + localOffset, "solver");
12527 }
12528 localOffset += contHilbertIdCount[i];
12529 } else {
12530
12531 file.setOffset(0, 0);
12532 file.writeArray(&partitionCellsId[0], "partitionCellsGlobalId");
12533 file.writeArray(&partitionCellsWorkload[0], "partitionCellsWorkload");
12534
12535 file.writeArray(&sliceCellInfo[0], "cellInfo");
12536 if(writeSolver) {
12537 file.writeArray(&solverBits[0], "solver");
12538 }
12539 }
12540 }
12541
12542 const MInt noContMinCells = contHilbertIdMinCellCount.size();
12543 for(
MInt i = 0, localMinOffset = 0; i < hilbertDomainOffset.size0(); i++) {
12544
12545 if(i < noContMinCells) {
12546
12547 file.setOffset(contHilbertIdMinCellCount[i], contHilbertIdMinCellOffset[i]);
12548 file.writeArray(&sliceMinLevelCellsTreeId[0] + localMinOffset, "minLevelCellsTreeId");
12549
12550 file.setOffset(contHilbertIdMinCellCount[i] * 2 * nDimSlice, contHilbertIdMinCellOffset[i] * 2 * nDimSlice);
12551 file.writeArray(&sliceMinLevelCellsNghbrIds[0] + localMinOffset * 2 * nDimSlice, "minLevelCellsNghbrIds");
12552 localMinOffset += contHilbertIdMinCellCount[i];
12553 } else {
12554
12555 file.setOffset(0, 0);
12556 file.writeArray(&partitionCellsId[0], "minLevelCellsTreeId");
12557 file.writeArray(&partitionCellsId[0], "minLevelCellsNghbrIds");
12558 }
12559 }
12560
12562 if(sliceDomain == 0) std::cerr << "Slice grid " << gridFileName << " write time: " << writeTimeTotal << std::endl;
12563 } else {
12564
12565
12566
12568
12569
12570 if(noSliceContHilbertIds != nullptr || sliceContiguousHilbertInfo != nullptr) {
12571 TERMM(1, "Error: using non-optimized slice IO but contiguous Hilbert info requested by passing non-nullptr as "
12572 "arguments.");
12573 }
12574
12575
12576
12577 for(
MInt i = 0, localOffset = 0, localPartOffset = 0, localMinOffset = 0; i < hilbertDomainOffset.size0(); i++) {
12578 if(i < noLocalHilbertIds[sliceDomain]) {
12579
12580
12581
12582 file.setOffset(hilbertIdPartitionCellsCount[hIdKeys[i]], hilbertPartitionDomainOffset[i]);
12583 file.writeArray(&partitionCellsId[0] + localPartOffset, "partitionCellsGlobalId");
12584 file.writeArray(&partitionCellsWorkload[0] + localPartOffset, "partitionCellsWorkload");
12585 localPartOffset += hilbertIdPartitionCellsCount[hIdKeys[i]];
12586
12587
12588
12589
12590
12591 file.setOffset(hilbertIdCount[hIdKeys[i]], hilbertDomainOffset[i]);
12592 file.writeArray(&sliceCellInfo[0] + localOffset, "cellInfo");
12593
12594
12595 if(writeSolver) {
12596 file.writeArray(&solverBits[0] + localOffset, "solver");
12597 }
12598
12599
12600 if(localMinOffset < noSliceMinLevelCells) {
12601
12602 file.setOffset(hilbertIdMinLevelCellsCount[hIdKeys[i]], hilbertMinLevelDomainOffset[i]);
12603 file.writeArray(&sliceMinLevelCellsTreeId[0] + localMinOffset, "minLevelCellsTreeId");
12604
12605 file.setOffset(hilbertIdMinLevelCellsCount[hIdKeys[i]] * 2 * nDimSlice,
12606 hilbertMinLevelDomainOffset[i] * 2 * nDimSlice);
12607 file.writeArray(&sliceMinLevelCellsNghbrIds[0] + localMinOffset * 2 * nDimSlice, "minLevelCellsNghbrIds");
12608 } else {
12609
12610 file.setOffset(0, 0);
12611 file.writeArray(&partitionCellsId[0], "minLevelCellsTreeId");
12612 file.writeArray(&partitionCellsId[0], "minLevelCellsNghbrIds");
12613 }
12614 localMinOffset += hilbertIdMinLevelCellsCount[hIdKeys[i]];
12615
12616 localOffset += hilbertIdCount[hIdKeys[i]];
12617 } else {
12618
12619
12620 file.setOffset(0, 0);
12621 file.writeArray(&partitionCellsId[0], "partitionCellsGlobalId");
12622 file.writeArray(&partitionCellsWorkload[0], "partitionCellsWorkload");
12623
12624 file.writeArray(&partitionCellsId[0], "minLevelCellsTreeId");
12625 file.writeArray(&partitionCellsId[0], "minLevelCellsNghbrIds");
12626
12627 file.writeArray(&sliceCellInfo[0], "cellInfo");
12628 if(writeSolver) {
12629 file.writeArray(&solverBits[0], "solver");
12630 }
12631 }
12632 }
12634 if(sliceDomain == 0) std::cerr << "Slice grid " << gridFileName << " write time: " << writeTimeTotal << std::endl;
12635 }
12636
12637
12638 file.close();
12639
12640
12642
12643 m_log <<
"done!" << endl;
12644}
MFloat m_partitionCellWorkloadThreshold
MInt maxPartitionLevelShift() const
MInt m_partitionCellOffspringThreshold
MBool m_hasMultiSolverBoundingBox
maia::grid::tree::SolverBitsetType::reference a_solver(const MInt cellId, const MInt solverId)
Returns if cell cellId is used by solver solverId.
MInt findNeighborDomainId(const MLong globalId)
void centerOfGravity(MFloat *const center) const
void boundingBox(MFloat *const box) const
const MLong & domainOffset(const MInt id) const
Return domain offset.
MFloat m_targetGridBoundingBox[6]
Multisolver grid information.
std::basic_string< char > MString
int MPI_Comm_free(MPI_Comm *comm, const MString &name, const MString &varname)
same as MPI_Comm_free, but updates the number of MPI communicators
int MPI_Comm_create(MPI_Comm comm, MPI_Group group, MPI_Comm *newcomm, const MString &name, const MString &varname)
same as MPI_Comm_create, but updates the number of MPI communicators
int MPI_Group_incl(MPI_Group group, int n, const int ranks[], MPI_Group *newgroup, const MString &name)
same as MPI_Group_incl
int MPI_Comm_group(MPI_Comm comm, MPI_Group *group, const MString &name, const MString &varname)
same as MPI_Comm_group
int MPI_Allreduce(const void *sendbuf, void *recvbuf, int count, MPI_Datatype datatype, MPI_Op op, MPI_Comm comm, const MString &name, const MString &sndvarname, const MString &rcvvarname)
same as MPI_Allreduce
int MPI_Group_free(MPI_Group *group, const MString &name)
same as MPI_Group_free
1.9.5