1264 {
1265 TRACE();
1266
1267 using namespace std;
1268
1275
1277 return;
1278 }
1279
1280
1282 isPartitionCell.fill(0);
1285 isPartitionCell(cellId) = 1;
1286 }
1287
1288
1289 MLong firstMinLevelCell = std::numeric_limits<MLong>::max();
1290 for(
MInt i = 0; i < noLocalMinLevelCells; i++) {
1291 firstMinLevelCell =
mMin(firstMinLevelCell,
a_globalId(localMinLevelCells[i]));
1292 }
1293
1294
1296 MPI_Allgather(&firstMinLevelCell, 1, MPI_LONG, minLevelCellGlobalIdOffsets.data(), 1, MPI_LONG,
mpiComm(), AT_,
1297 "firstMinLevelCell", "minLevelCellGlobalIdOffsets.data()");
1299
1302 minLevelCellGlobalIdOffsets[cpu] = minLevelCellGlobalIdOffsets[cpu + 1];
1303 }
1304 }
1305
1306 std::vector<std::tuple<MLong, MInt, MInt>> partitionLevelAncestorsOnMinLevel;
1307 std::vector<std::tuple<MLong, MInt, MInt>> partitionLevelAncestorsOther;
1309 isPartitionLevelAncestor.fill(1);
1310
1311 MInt noPartitionLevelAncestorsGlobal = 0;
1314 "Error: tree should only contain internal cells at this moment.");
1315
1316
1319 noPartitionLevelAncestorsLocal -=
1321 &isPartitionLevelAncestor[0]);
1322 }
1323
1325
1326 MPI_Allreduce(&noPartitionLevelAncestorsLocal, &noPartitionLevelAncestorsGlobal, 1, MPI_INT, MPI_SUM,
mpiComm(),
1327 AT_, "noPartitionLevelAncestorsLocal", "noPartitionLevelAncestorsGlobal");
1329
1331 a_hasProperty(i, Cell::IsPartLvlAncestor) = isPartitionLevelAncestor(i);
1332 }
1333
1336 if(localMinLevelId[cellId] < 0) {
1337 isPartitionLevelAncestor[
cellId] = 1;
1338 }
1339 }
1340
1343 if(isPartitionLevelAncestor(cellId)) {
1344 if(localMinLevelId[cellId] > -1) {
1345
1346 partitionLevelAncestorsOnMinLevel.push_back(std::make_tuple(globalId, cellId,
domainId()));
1347 } else {
1348
1349
1352 if(globalId >= minLevelCellGlobalIdOffsets[cpu] && globalId < minLevelCellGlobalIdOffsets[cpu + 1]) {
1353 ndom = cpu;
1354 break;
1355 }
1356 }
1357 TERMM_IF_COND(ndom < 0 || ndom >=
noDomains(),
"domain for globalId not found " + std::to_string(ndom));
1358 partitionLevelAncestorsOther.push_back(std::make_tuple(globalId, cellId, ndom));
1359 }
1360 }
1361 }
1362
1365 noQuery.fill(0);
1366
1367
1368 for(
MUint i = 0; i < partitionLevelAncestorsOther.size(); i++) {
1369 noQuery[get<2>(partitionLevelAncestorsOther[i])]++;
1370 }
1371 queryOffsets(0) = 0;
1373 queryOffsets(c + 1) = queryOffsets(c) + noQuery[c];
1374 }
1375
1379 noQuery.fill(0);
1380
1381
1382 for(
MUint i = 0; i < partitionLevelAncestorsOther.size(); i++) {
1383 MLong globalId = get<0>(partitionLevelAncestorsOther[i]);
1384 MInt cellId = get<1>(partitionLevelAncestorsOther[i]);
1385 MInt cpu = get<2>(partitionLevelAncestorsOther[i]);
1386 queryGlobalId(queryOffsets(cpu) + noQuery(cpu)) = globalId;
1387 queryCellInfo(queryOffsets(cpu) + noQuery(cpu)) = cellInfo[
cellId];
1388 queryIsPartitionCell(queryOffsets(cpu) + noQuery(cpu)) = isPartitionCell[
cellId];
1389
1390 noQuery[cpu]++;
1391 }
1392
1395
1396
1397 MPI_Alltoall(&noQuery[0], 1, MPI_INT, &noCollect[0], 1, MPI_INT,
mpiComm(), AT_,
"noQuery[0]",
"noCollect[0]");
1398
1399 collectOffsets(0) = 0;
1401 collectOffsets(c + 1) = collectOffsets(c) + noCollect[c];
1402 }
1403
1404
1409 queryReq.fill(MPI_REQUEST_NULL);
1410
1411
1414 if(noQuery[c] == 0) continue;
1415 MPI_Issend(&queryGlobalId[queryOffsets[c]], noQuery[c], MPI_LONG, c, 456,
mpiComm(), &queryReq[queryCnt++], AT_,
1416 "queryGlobalId[queryOffsets[c]]");
1417 MPI_Issend(&queryCellInfo[queryOffsets[c]], noQuery[c], MPI_UNSIGNED_CHAR, c, 457,
mpiComm(),
1418 &queryReq[queryCnt++], AT_, "queryCellInfo[queryOffsets[c]]");
1419 MPI_Issend(&queryIsPartitionCell[queryOffsets[c]], noQuery[c], MPI_INT, c, 458,
mpiComm(), &queryReq[queryCnt++],
1420 AT_, "queryIsPartitionCell[queryOffsets[c]]");
1421 }
1422
1423
1424 MInt collectCnt = 0;
1426 if(noCollect[c] == 0) continue;
1427 MPI_Recv(&collectGlobalId[collectOffsets[c]], noCollect[c], MPI_LONG, c, 456,
mpiComm(), MPI_STATUS_IGNORE, AT_,
1428 "collectGlobalId[collectOffsets[c]]");
1429 MPI_Recv(&collectCellInfo[collectOffsets[c]], noCollect[c], MPI_UNSIGNED_CHAR, c, 457,
mpiComm(),
1430 MPI_STATUS_IGNORE, AT_, "collectCellInfo[collectOffsets[c]]");
1431 MPI_Recv(&collectIsPartitionCell[collectOffsets[c]], noCollect[c], MPI_INT, c, 458,
mpiComm(), MPI_STATUS_IGNORE,
1432 AT_, "collectIsPartitionCell[collectOffsets[c]]");
1433 collectCnt++;
1434 }
1435
1436 if(queryCnt > 0)
MPI_Waitall(queryCnt, &queryReq[0], MPI_STATUSES_IGNORE, AT_);
1437
1438
1439 std::vector<std::tuple<MLong, MUchar, MInt>> collectData;
1440 for(
MUint i = 0; i < partitionLevelAncestorsOnMinLevel.size(); i++) {
1441 MLong globalId = get<0>(partitionLevelAncestorsOnMinLevel[i]);
1442 MInt localId = get<1>(partitionLevelAncestorsOnMinLevel[i]);
1443 ASSERT(
domainId() == get<2>(partitionLevelAncestorsOnMinLevel[i]),
"");
1444 collectData.push_back(std::make_tuple(globalId, cellInfo[localId], isPartitionCell[localId]));
1445 }
1446
1447
1449 for(
MInt d = 0; d < noCollect[c]; d++) {
1450 MInt id = collectOffsets[c] + d;
1451 collectData.push_back(std::make_tuple(collectGlobalId[id], collectCellInfo[id], collectIsPartitionCell[id]));
1452 }
1453 }
1454
1455 sort(collectData.begin(), collectData.end());
1456
1457
1459 MLongScratchSpace collectMinLevelTreeId(collectData.size(), AT_,
"collectMinLevelTreeId");
1462 MIntScratchSpace collectNoChildren(collectData.size(), AT_,
"collectNoChildren");
1463 MIntScratchSpace collectNoOffspring(collectData.size(), AT_,
"collectNoOffspring");
1466 std::vector<MInt> newCells;
1468 std::map<MLong, MInt> collectMap;
1469 collectMinLevelTreeId.fill(-1);
1470 collectParentId.fill(-1);
1471 collectLevel.fill(-1);
1472 collectNoChildren.fill(-1);
1473 collectNoOffspring.fill(-1);
1474 collectChildIds.fill(-1);
1475 collectNghbrIds.fill(-1);
1476
1477
1478 for(
MUint i = 0; i < collectData.size(); i++) {
1479 collectMap[get<0>(collectData[i])] = i;
1480 }
1481
1482
1483 for(
MUint i = 0; i < partitionLevelAncestorsOnMinLevel.size(); i++) {
1484 MLong globalId = get<0>(partitionLevelAncestorsOnMinLevel[i]);
1485 MInt localId = get<1>(partitionLevelAncestorsOnMinLevel[i]);
1486 MInt id = collectMap[globalId];
1488 collectParentId(id) = -1;
1491 }
1492 MInt minId = localMinLevelId[localId];
1493 ASSERT(minId > -1, "");
1494 collectMinLevelTreeId(id) = minLevelCellsTreeId[minId];
1495 collectNoOffspring(id) = (minId == noLocalMinLevelCells - 1)
1496 ? minLevelCellGlobalIdOffsets[
domainId() + 1] - globalId
1497 :
a_globalId(localMinLevelCells[minId + 1]) - globalId;
1498 }
1499
1500
1501
1502 for(
MUint i = 0; i < partitionLevelAncestorsOnMinLevel.size(); i++) {
1503 MInt counter = collectMap[get<0>(partitionLevelAncestorsOnMinLevel[i])];
1505 &collectNoChildren[0], &collectChildIds[0], &collectNoOffspring[0]);
1506 }
1507
1508
1509 for(
MUint i = 0; i < partitionLevelAncestorsOnMinLevel.size(); i++) {
1510 MLong globalId = get<0>(partitionLevelAncestorsOnMinLevel[i]);
1511 MInt localId = get<1>(partitionLevelAncestorsOnMinLevel[i]);
1512 MInt id = collectMap[globalId];
1515 a_childId(localId, j) = collectChildIds(
id, j);
1516 }
1517 }
1518
1519
1521 std::vector<MLong> returnData;
1522
1523 if(collectCnt > 0) {
1524 returnOffsets(0) = 0;
1525 collectCnt = 0;
1527 if(noCollect[c] > 0) {
1528
1529 std::map<MLong, MInt> returnChain;
1530
1531 for(
MInt d = 0; d < noCollect[c]; d++) {
1532 const MInt id0 = collectOffsets[c] + d;
1533 MLong globalId = collectGlobalId[id0];
1534 MInt id1 = collectMap[globalId];
1535 if(collectIsPartitionCell[id0]) {
1536 id1 = collectMap[collectParentId[id1]];
1537 globalId = get<0>(collectData[id1]);
1538 }
1539 while(globalId > -1) {
1540 if(returnChain.count(globalId) == 0) {
1541 returnChain[globalId] = id1;
1542 if(collectParentId[id1] > -1) {
1543 id1 = collectMap[collectParentId[id1]];
1544 globalId = get<0>(collectData[id1]);
1545 } else {
1546 globalId = -1;
1547 }
1548 } else {
1549 globalId = -1;
1550 }
1551 }
1552 }
1553
1554 for(auto it = returnChain.begin(); it != returnChain.end(); it++) {
1555 MLong globalId = it->first;
1556 MInt id1 = it->second;
1557
1558
1563
1570 } else {
1571
1572 localId = it1->second;
1574 ASSERT(
a_hasProperty(localId, Cell::IsPartLvlAncestor),
"");
1575 }
1578 a_level(localId) = collectLevel(id1);
1579
1581 a_childId(localId, i) = collectChildIds(id1, i);
1583
1586 }
1587 }
1588
1589 newCells.push_back(localId);
1590 newCellsMaxLevel =
mMax(newCellsMaxLevel,
a_level(localId));
1591 } else {
1592
1593 MInt offs = returnData.size();
1594 returnData.resize(offs + noData);
1595 returnData[offs++] = globalId;
1596 returnData[offs++] = collectMinLevelTreeId(id1);
1597 returnData[offs++] = collectLevel(id1);
1598 returnData[offs++] = collectParentId(id1);
1599 returnData[offs++] = collectNoOffspring(id1);
1601 returnData[offs++] = collectChildIds(id1, i);
1602 }
1604 returnData[offs++] = collectNghbrIds(id1, i);
1605 }
1606 }
1607 }
1609 returnOffsets(collectCnt + 1) = returnData.size();
1610 collectCnt++;
1611 }
1612 }
1613 }
1614 }
1615
1618 returnReq.fill(MPI_REQUEST_NULL);
1620
1621
1622 if(collectCnt > 0) {
1625 if(noCollect[c] == 0) continue;
1626
1627 noReturnDataBuffer[returnCnt] = returnOffsets(returnCnt + 1) - returnOffsets(returnCnt);
1628 MPI_Issend(&noReturnDataBuffer[returnCnt], 1, MPI_INT, c, 459,
mpiComm(), &returnReq[returnCnt], AT_,
1629 "&noReturnDataBuffer[returnCnt]");
1630 returnCnt++;
1631 }
1632 }
1633
1634 MIntScratchSpace noReceiveData(std::max(queryCnt, 1), AT_,
"noReceiveData");
1636 queryCnt = 0;
1637 receiveOffs(0) = 0;
1638
1641 if(noQuery[c] == 0) continue;
1642 MPI_Recv(&noReceiveData[queryCnt], 1, MPI_INT, c, 459,
mpiComm(), MPI_STATUS_IGNORE, AT_,
1643 "noReceiveData[queryCnt]");
1644 receiveOffs(queryCnt + 1) = receiveOffs(queryCnt) + noReceiveData[queryCnt];
1645 queryCnt++;
1646 }
1647
1648 if(returnCnt > 0)
MPI_Waitall(returnCnt, &returnReq[0], MPI_STATUSES_IGNORE, AT_);
1649
1650 MLongScratchSpace receiveData(std::max(receiveOffs(queryCnt), 1), AT_,
"receiveData");
1651
1652 returnReq.fill(MPI_REQUEST_NULL);
1653 returnCnt = 0;
1654
1655 if(collectCnt > 0) {
1658 if(noCollect[c] == 0) continue;
1659 MInt noReturnData = returnOffsets(returnCnt + 1) - returnOffsets(returnCnt);
1660 MInt offs = returnOffsets(returnCnt);
1661 MPI_Issend(&returnData[offs], noReturnData, MPI_LONG, c, 460,
mpiComm(), &returnReq[returnCnt++], AT_,
1662 "returnData[offs]");
1663 }
1664 }
1665 queryCnt = 0;
1668 if(noQuery[c] == 0) continue;
1669 MPI_Recv(&receiveData(receiveOffs(queryCnt)), noReceiveData(queryCnt), MPI_LONG, c, 460,
mpiComm(),
1670 MPI_STATUS_IGNORE, AT_, "receiveData(receiveOffs(queryCnt))");
1671 queryCnt++;
1672 }
1673
1674 if(returnCnt > 0)
MPI_Waitall(returnCnt, &returnReq[0], MPI_STATUSES_IGNORE, AT_);
1675
1676 queryCnt = 0;
1679 if(noQuery[c] == 0) continue;
1680 MInt offs = receiveOffs(queryCnt);
1681
1682 while(offs < receiveOffs(queryCnt + 1)) {
1683 MLong globalId = receiveData(offs++);
1687
1694 } else {
1695
1696 localId = it->second;
1698 ASSERT(
a_hasProperty(localId, Cell::IsPartLvlAncestor),
"");
1699 }
1700 MLong treeId = receiveData(offs++);
1701 a_level(localId) = receiveData(offs++);
1704
1705 newCells.push_back(localId);
1706 newCellsMaxLevel =
mMax(newCellsMaxLevel,
a_level(localId));
1707
1709 a_childId(localId, i) = receiveData(offs++);
1711
1714 }
1715 }
1718 }
1719
1721 maia::grid::hilbert::treeIdToCoordinates<nDim>(
1723 &
grid_.m_targetGridCenterOfGravity[0],
grid_.m_targetGridLengthLevel0);
1724 }
1725 }
1726 queryCnt++;
1727 }
1728
1730
1734 if(localMinLevelId[cellId] < 0) {
1737 }
1738 }
1739
1740 static const MFloat childStencil[8][3] = {{-F1, -F1, -F1}, {F1, -F1, -F1}, {-F1, F1, -F1}, {F1, F1, -F1},
1741 {-F1, -F1, F1}, {F1, -F1, F1}, {-F1, F1, F1}, {F1, F1, F1}};
1742
1744 for(
MUint i = 0; i < newCells.size(); i++) {
1745 if(
a_level(newCells[i]) == lvl) {
1747 if(
a_childId(newCells[i], child) > -1) {
1753 }
1754 }
1755 }
1756 }
1757 }
1758 }
1759 }
1760
1763 }
1764
1765 for(
MUint i = 0; i < newCells.size(); i++) {
1768
1772 }
1773
1779 for(
MInt d = firstDom; d <= lastDom; d++) {
1782 if(idx < 0) {
1787 }
1788 }
1789 } else {
1792 if(idx < 0) {
1797 }
1798 }
1799 }
1800
1803 .end());
1804
1805 for(
MUint i = 0; i < newCells.size(); i++) {
1813 for(
MInt d = firstDom; d <= lastDom; d++) {
1816 ASSERT(idx > -1, "");
1818 }
1819 } else {
1822 ASSERT(idx > -1, "");
1824 }
1825 }
1826
1827
1830 mTerm(1, AT_,
"Invalid sorting of neighbor domains.");
1831 }
1832 }
1836 mTerm(1, AT_,
"Invalid sorting of halo cells.");
1837 }
1838 }
1841 mTerm(1, AT_,
"Invalid sorting of window cells.");
1842 }
1843 }
1844 }
1845 }
static constexpr MInt m_noDirs
MLong & a_neighborId(const MInt cellId, const MInt dir)
MInt traverseAndFlagSubtree(const MUint &cellId, const MUchar *const cellInfo, const MInt N, MInt *const flag)
Recursively flag subtree cells from partition level to the leaf level.
static constexpr MInt m_maxNoChilds
MFloat & a_coordinate(const MInt cellId, const MInt dir)
MFloat cellLengthAtLevel(const MInt level) const
MInt findNeighborDomainId(MLong globalId)
MInt findIndex(ITERATOR first, ITERATOR last, const U &val)
MLong & a_parentId(const MInt cellId)
void resetCell(const MInt &cellId)
void traverseAndSetupTruncatedSubtree(MInt &counter, const std::vector< std::tuple< MLong, MUchar, MInt > > &collectData, MLong *const collectParentId, MInt *const collectLevel, MInt *const collectNoChildren, MLong *const collectChildIds, MInt *const collectNoOffspring)
Recursively setup truncated subtree between minLevel and partitionLevel (used to setup partition leve...
void mTerm(const MInt errorCode, const MString &location, const MString &message)
int MPI_Issend(const void *buf, int count, MPI_Datatype datatype, int dest, int tag, MPI_Comm comm, MPI_Request *request, const MString &name, const MString &varname)
same as MPI_Issend
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_Alltoall(const void *sendbuf, int sendcount, MPI_Datatype sendtype, void *recvbuf, int recvcount, MPI_Datatype recvtype, MPI_Comm comm, const MString &name, const MString &sndvarname, const MString &rcvvarname)
same as MPI_Alltoall
int MPI_Allgather(const void *sendbuf, int sendcount, MPI_Datatype sendtype, void *recvbuf, int recvcount, MPI_Datatype recvtype, MPI_Comm comm, const MString &name, const MString &sndvarname, const MString &rcvvarname)
same as MPI_Allgather
1.9.5