ENH: regionSplit: improve algorithm order

The old version of regionSplit would hand out regions one by one. This
is a big problem when there are lots of regions - the extreme being
in the decompositionMethods, where it is used to cluster cells and most clusters
being only one cell. This rewrite uses a mesh wave to determine disconnected
regions in one go. This produced non-compact numbering which is then compacted
in a second phase.
On a 14M cell case with cyclic constraints this reduced decompose
time from 40 mins down to 5.
This commit is contained in:
mattijs
2015-11-26 16:52:18 +00:00
parent 710f189456
commit 8d0154ba32
2 changed files with 197 additions and 497 deletions

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -28,6 +28,8 @@ License
#include "processorPolyPatch.H" #include "processorPolyPatch.H"
#include "globalIndex.H" #include "globalIndex.H"
#include "syncTools.H" #include "syncTools.H"
#include "FaceCellWave.H"
#include "minData.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -40,340 +42,93 @@ defineTypeNameAndDebug(regionSplit, 0);
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
// Handle (non-processor) coupled faces. void Foam::regionSplit::calcNonCompactRegionSplit
void Foam::regionSplit::transferCoupledFaceRegion
(
const label faceI,
const label otherFaceI,
labelList& faceRegion,
DynamicList<label>& newChangedFaces
) const
{
if (faceRegion[faceI] >= 0)
{
if (faceRegion[otherFaceI] == -1)
{
faceRegion[otherFaceI] = faceRegion[faceI];
newChangedFaces.append(otherFaceI);
}
else if (faceRegion[otherFaceI] == -2)
{
// otherFaceI blocked but faceI is not. Is illegal for coupled
// faces, not for explicit connections.
}
else if (faceRegion[otherFaceI] != faceRegion[faceI])
{
FatalErrorIn
(
"regionSplit::transferCoupledFaceRegion"
"(const label, const label, labelList&, labelList&) const"
) << "Problem : coupled face " << faceI
<< " on patch " << mesh().boundaryMesh().whichPatch(faceI)
<< " has region " << faceRegion[faceI]
<< " but coupled face " << otherFaceI
<< " has region " << faceRegion[otherFaceI]
<< endl
<< "Is your blocked faces specification"
<< " synchronized across coupled boundaries?"
<< abort(FatalError);
}
}
else if (faceRegion[faceI] == -1)
{
if (faceRegion[otherFaceI] >= 0)
{
faceRegion[faceI] = faceRegion[otherFaceI];
newChangedFaces.append(faceI);
}
else if (faceRegion[otherFaceI] == -2)
{
// otherFaceI blocked but faceI is not. Is illegal for coupled
// faces, not for explicit connections.
}
}
}
void Foam::regionSplit::fillSeedMask
(
const List<labelPair>& explicitConnections,
labelList& cellRegion,
labelList& faceRegion,
const label seedCellID,
const label markValue
) const
{
// Do seed cell
cellRegion[seedCellID] = markValue;
// Collect faces on seed cell
const cell& cFaces = mesh().cells()[seedCellID];
label nFaces = 0;
labelList changedFaces(cFaces.size());
forAll(cFaces, i)
{
label faceI = cFaces[i];
if (faceRegion[faceI] == -1)
{
faceRegion[faceI] = markValue;
changedFaces[nFaces++] = faceI;
}
}
changedFaces.setSize(nFaces);
// Loop over changed faces. FaceCellWave in small.
while (changedFaces.size())
{
//if (debug)
//{
// Pout<< "regionSplit::fillSeedMask : changedFaces:"
// << changedFaces.size() << endl;
//}
DynamicList<label> changedCells(changedFaces.size());
forAll(changedFaces, i)
{
label faceI = changedFaces[i];
label own = mesh().faceOwner()[faceI];
if (cellRegion[own] == -1)
{
cellRegion[own] = markValue;
changedCells.append(own);
}
if (mesh().isInternalFace(faceI))
{
label nei = mesh().faceNeighbour()[faceI];
if (cellRegion[nei] == -1)
{
cellRegion[nei] = markValue;
changedCells.append(nei);
}
}
}
//if (debug)
//{
// Pout<< "regionSplit::fillSeedMask : changedCells:"
// << changedCells.size() << endl;
//}
// Loop over changedCells and collect faces
DynamicList<label> newChangedFaces(changedCells.size());
forAll(changedCells, i)
{
label cellI = changedCells[i];
const cell& cFaces = mesh().cells()[cellI];
forAll(cFaces, cFaceI)
{
label faceI = cFaces[cFaceI];
if (faceRegion[faceI] == -1)
{
faceRegion[faceI] = markValue;
newChangedFaces.append(faceI);
}
}
}
//if (debug)
//{
// Pout<< "regionSplit::fillSeedMask : changedFaces before sync:"
// << changedFaces.size() << endl;
//}
// Check for changes to any locally coupled face.
// Global connections are done later.
const polyBoundaryMesh& patches = mesh().boundaryMesh();
forAll(patches, patchI)
{
const polyPatch& pp = patches[patchI];
if
(
isA<cyclicPolyPatch>(pp)
&& refCast<const cyclicPolyPatch>(pp).owner()
)
{
// Transfer from neighbourPatch to here or vice versa.
const cyclicPolyPatch& cycPatch =
refCast<const cyclicPolyPatch>(pp);
label faceI = cycPatch.start();
forAll(cycPatch, i)
{
label otherFaceI = cycPatch.transformGlobalFace(faceI);
transferCoupledFaceRegion
(
faceI,
otherFaceI,
faceRegion,
newChangedFaces
);
faceI++;
}
}
}
forAll(explicitConnections, i)
{
transferCoupledFaceRegion
(
explicitConnections[i][0],
explicitConnections[i][1],
faceRegion,
newChangedFaces
);
}
//if (debug)
//{
// Pout<< "regionSplit::fillSeedMask : changedFaces after sync:"
// << newChangedFaces.size() << endl;
//}
changedFaces.transfer(newChangedFaces);
}
}
Foam::label Foam::regionSplit::calcLocalRegionSplit
( (
const globalIndex& globalFaces,
const boolList& blockedFace, const boolList& blockedFace,
const List<labelPair>& explicitConnections, const List<labelPair>& explicitConnections,
labelList& cellRegion labelList& cellRegion
) const ) const
{ {
if (debug) // Field on cells and faces.
{ List<minData> cellData(mesh().nCells());
if (blockedFace.size()) List<minData> faceData(mesh().nFaces());
{
// Check that blockedFace is synced.
boolList syncBlockedFace(blockedFace);
syncTools::swapFaceList(mesh(), syncBlockedFace);
forAll(syncBlockedFace, faceI) // Take over blockedFaces by seeding a negative number
{ // (so is always less than the decomposition)
if (syncBlockedFace[faceI] != blockedFace[faceI]) label nUnblocked = 0;
{ forAll(faceData, faceI)
FatalErrorIn {
( if (blockedFace.size() && blockedFace[faceI])
"regionSplit::calcLocalRegionSplit(..)" {
) << "Face " << faceI << " not synchronised. My value:" faceData[faceI] = minData(-2);
<< blockedFace[faceI] << " coupled value:" }
<< syncBlockedFace[faceI] else
<< abort(FatalError); {
} nUnblocked++;
}
} }
} }
// Region per face. // Seed unblocked faces
// -1 unassigned labelList seedFaces(nUnblocked);
// -2 blocked List<minData> seedData(nUnblocked);
labelList faceRegion(mesh().nFaces(), -1); nUnblocked = 0;
if (blockedFace.size())
forAll(faceData, faceI)
{ {
forAll(blockedFace, faceI) if (blockedFace.empty() || !blockedFace[faceI])
{ {
if (blockedFace[faceI]) seedFaces[nUnblocked] = faceI;
{ // Seed face with globally unique number
faceRegion[faceI] = -2; seedData[nUnblocked] = minData(globalFaces.toGlobal(faceI));
} nUnblocked++;
} }
} }
// Assign local regions // Propagate information inwards
// ~~~~~~~~~~~~~~~~~~~~ FaceCellWave<minData> deltaCalc
(
mesh(),
explicitConnections,
false, // disable walking through cyclicAMI for backwards compatibility
seedFaces,
seedData,
faceData,
cellData,
mesh().globalData().nTotalCells()+1
);
// Start with region 0
label nLocalRegions = 0;
label unsetCellI = 0; // And extract
cellRegion.setSize(mesh().nCells());
do forAll(cellRegion, cellI)
{ {
// Find first unset cell if (cellData[cellI].valid(deltaCalc.data()))
for (; unsetCellI < mesh().nCells(); unsetCellI++)
{ {
if (cellRegion[unsetCellI] == -1) cellRegion[cellI] = cellData[cellI].data();
{
break;
}
} }
else
if (unsetCellI >= mesh().nCells())
{ {
break; // Unvisited cell -> only possible if surrounded by blocked faces.
} // If so make up region from any of the faces
const cell& cFaces = mesh().cells()[cellI];
label faceI = cFaces[0];
fillSeedMask if (blockedFace.size() && !blockedFace[faceI])
(
explicitConnections,
cellRegion,
faceRegion,
unsetCellI,
nLocalRegions
);
// Current unsetCell has now been handled. Go to next region.
nLocalRegions++;
unsetCellI++;
}
while (true);
if (debug)
{
forAll(cellRegion, cellI)
{
if (cellRegion[cellI] < 0)
{ {
FatalErrorIn("regionSplit::calcLocalRegionSplit(..)") FatalErrorIn("regionSplit::calcNonCompactRegionSplit(..)")
<< "cell:" << cellI << " region:" << cellRegion[cellI] << "Problem: unblocked face " << faceI
<< abort(FatalError); << " at " << mesh().faceCentres()[faceI]
} << " on unassigned cell " << cellI
} << mesh().cellCentres()[faceI]
<< exit(FatalError);
forAll(faceRegion, faceI)
{
if (faceRegion[faceI] == -1)
{
FatalErrorIn("regionSplit::calcLocalRegionSplit(..)")
<< "face:" << faceI << " region:" << faceRegion[faceI]
<< abort(FatalError);
} }
cellRegion[cellI] = globalFaces.toGlobal(faceI);
} }
} }
return nLocalRegions;
} }
@ -388,176 +143,142 @@ Foam::autoPtr<Foam::globalIndex> Foam::regionSplit::calcRegionSplit
{ {
// See header in regionSplit.H // See header in regionSplit.H
// 1. Do local analysis
label nLocalRegions = calcLocalRegionSplit if (!doGlobalRegions)
{
// Block all parallel faces to avoid comms across
boolList coupledOrBlockedFace(blockedFace);
const polyBoundaryMesh& pbm = mesh().boundaryMesh();
if (coupledOrBlockedFace.size())
{
forAll(pbm, patchI)
{
const polyPatch& pp = pbm[patchI];
if (isA<processorPolyPatch>(pp))
{
label faceI = pp.start();
forAll(pp, i)
{
coupledOrBlockedFace[faceI++] = true;
}
}
}
}
// Create dummy (local only) globalIndex
labelList offsets(Pstream::nProcs()+1, 0);
for (label i = Pstream::myProcNo()+1; i < offsets.size(); i++)
{
offsets[i] = mesh().nFaces();
}
const globalIndex globalRegions(offsets.xfer());
// Minimise regions across connected cells
// Note: still uses global decisions so all processors are running
// in lock-step, i.e. slowest determines overall time.
// To avoid this we could switch off Pstream::parRun.
calcNonCompactRegionSplit
(
globalRegions,
coupledOrBlockedFace,
explicitConnections,
cellRegion
);
// Compact
Map<label> globalToCompact(mesh().nCells()/8);
forAll(cellRegion, cellI)
{
label region = cellRegion[cellI];
label globalRegion;
Map<label>::const_iterator fnd = globalToCompact.find(region);
if (fnd == globalToCompact.end())
{
globalRegion = globalRegions.toGlobal(globalToCompact.size());
globalToCompact.insert(region, globalRegion);
}
else
{
globalRegion = fnd();
}
cellRegion[cellI] = globalRegion;
}
// Return globalIndex with size = localSize and all regions local
labelList compactOffsets(Pstream::nProcs()+1, 0);
for (label i = Pstream::myProcNo()+1; i < compactOffsets.size(); i++)
{
compactOffsets[i] = globalToCompact.size();
}
return autoPtr<globalIndex>(new globalIndex(compactOffsets.xfer()));
}
// Initial global region numbers
const globalIndex globalRegions(mesh().nFaces());
// Minimise regions across connected cells (including parallel)
calcNonCompactRegionSplit
( (
globalRegions,
blockedFace, blockedFace,
explicitConnections, explicitConnections,
cellRegion cellRegion
); );
if (!doGlobalRegions)
// Now our cellRegion will have
// - non-local regions (i.e. originating from other processors)
// - non-compact locally originating regions
// so we'll need to compact
// 4a: count per originating processor the number of regions
labelList nOriginating(Pstream::nProcs(), 0);
{ {
return autoPtr<globalIndex>(new globalIndex(nLocalRegions)); labelHashSet haveRegion(mesh().nCells()/8);
}
// 2. Assign global regions
// ~~~~~~~~~~~~~~~~~~~~~~~~
// Offset local regions to create unique global regions.
globalIndex globalRegions(nLocalRegions);
// Convert regions to global ones
forAll(cellRegion, cellI)
{
cellRegion[cellI] = globalRegions.toGlobal(cellRegion[cellI]);
}
// 3. Merge global regions
// ~~~~~~~~~~~~~~~~~~~~~~~
// Regions across non-blocked proc patches get merged.
// This will set merged global regions to be the min of both.
// (this will create gaps in the global region list so they will get
// merged later on)
while (true)
{
if (debug)
{
Pout<< nl << "-- Starting Iteration --" << endl;
}
const polyBoundaryMesh& patches = mesh().boundaryMesh();
labelList nbrRegion(mesh().nFaces()-mesh().nInternalFaces(), -1);
forAll(patches, patchI)
{
const polyPatch& pp = patches[patchI];
if (pp.coupled())
{
const labelList& patchCells = pp.faceCells();
SubList<label> patchNbrRegion
(
nbrRegion,
pp.size(),
pp.start()-mesh().nInternalFaces()
);
forAll(patchCells, i)
{
label faceI = pp.start()+i;
if (!blockedFace.size() || !blockedFace[faceI])
{
patchNbrRegion[i] = cellRegion[patchCells[i]];
}
}
}
}
syncTools::swapBoundaryFaceList(mesh(), nbrRegion);
Map<label> globalToMerged(mesh().nFaces()-mesh().nInternalFaces());
forAll(patches, patchI)
{
const polyPatch& pp = patches[patchI];
if (pp.coupled())
{
const labelList& patchCells = pp.faceCells();
SubList<label> patchNbrRegion
(
nbrRegion,
pp.size(),
pp.start()-mesh().nInternalFaces()
);
forAll(patchCells, i)
{
label faceI = pp.start()+i;
if (!blockedFace.size() || !blockedFace[faceI])
{
if (patchNbrRegion[i] < cellRegion[patchCells[i]])
{
//Pout<< "on patch:" << pp.name()
// << " cell:" << patchCells[i]
// << " at:"
// << mesh().cellCentres()[patchCells[i]]
// << " was:" << cellRegion[patchCells[i]]
// << " nbr:" << patchNbrRegion[i]
// << endl;
globalToMerged.insert
(
cellRegion[patchCells[i]],
patchNbrRegion[i]
);
}
}
}
}
}
label nMerged = returnReduce(globalToMerged.size(), sumOp<label>());
if (debug)
{
Pout<< "nMerged:" << nMerged << endl;
}
if (nMerged == 0)
{
break;
}
// Renumber the regions according to the globalToMerged
forAll(cellRegion, cellI) forAll(cellRegion, cellI)
{ {
label regionI = cellRegion[cellI]; label region = cellRegion[cellI];
Map<label>::const_iterator iter = globalToMerged.find(regionI);
if (iter != globalToMerged.end()) // Count originating processor. Use isLocal as efficiency since
// most cells are locally originating.
if (globalRegions.isLocal(region))
{ {
cellRegion[cellI] = iter(); if (haveRegion.insert(region))
{
nOriginating[Pstream::myProcNo()]++;
}
} }
} else
}
// Now our cellRegion will have non-local elements in it. So compact
// it.
// 4a: count. Use a labelHashSet to count regions only once.
label nCompact = 0;
{
labelHashSet localRegion(mesh().nFaces()-mesh().nInternalFaces());
forAll(cellRegion, cellI)
{
if
(
globalRegions.isLocal(cellRegion[cellI])
&& localRegion.insert(cellRegion[cellI])
)
{ {
nCompact++; label procI = globalRegions.whichProcID(region);
if (haveRegion.insert(region))
{
nOriginating[procI]++;
}
} }
} }
} }
if (debug) if (debug)
{ {
Pout<< "Compacted from " << nLocalRegions Pout<< "Counted " << nOriginating[Pstream::myProcNo()]
<< " down to " << nCompact << " local regions." << endl; << " local regions." << endl;
} }
// Global numbering for compacted local regions // Global numbering for compacted local regions
autoPtr<globalIndex> globalCompactPtr(new globalIndex(nCompact)); autoPtr<globalIndex> globalCompactPtr
(
new globalIndex(nOriginating[Pstream::myProcNo()])
);
const globalIndex& globalCompact = globalCompactPtr(); const globalIndex& globalCompact = globalCompactPtr();
@ -568,12 +289,15 @@ Foam::autoPtr<Foam::globalIndex> Foam::regionSplit::calcRegionSplit
// labelList. // labelList.
// Local compaction map // Local compaction map
Map<label> globalToCompact(2*nCompact); Map<label> globalToCompact(2*nOriginating[Pstream::myProcNo()]);
// Remote regions we want the compact number for // Remote regions we want the compact number for
List<labelHashSet> nonLocal(Pstream::nProcs()); List<labelHashSet> nonLocal(Pstream::nProcs());
forAll(nonLocal, procI) forAll(nonLocal, procI)
{ {
nonLocal[procI].resize((nLocalRegions-nCompact)/Pstream::nProcs()); if (procI != Pstream::myProcNo())
{
nonLocal[procI].resize(2*nOriginating[procI]);
}
} }
forAll(cellRegion, cellI) forAll(cellRegion, cellI)
@ -581,15 +305,12 @@ Foam::autoPtr<Foam::globalIndex> Foam::regionSplit::calcRegionSplit
label region = cellRegion[cellI]; label region = cellRegion[cellI];
if (globalRegions.isLocal(region)) if (globalRegions.isLocal(region))
{ {
Map<label>::const_iterator iter = globalToCompact.find(region); // Insert new compact region (if not yet present)
if (iter == globalToCompact.end()) globalToCompact.insert
{ (
label compactRegion = globalCompact.toGlobal region,
( globalCompact.toGlobal(globalToCompact.size())
globalToCompact.size() );
);
globalToCompact.insert(region, compactRegion);
}
} }
else else
{ {
@ -603,22 +324,17 @@ Foam::autoPtr<Foam::globalIndex> Foam::regionSplit::calcRegionSplit
// Convert the nonLocal (labelHashSets) to labelLists. // Convert the nonLocal (labelHashSets) to labelLists.
labelListList sendNonLocal(Pstream::nProcs()); labelListList sendNonLocal(Pstream::nProcs());
labelList nNonLocal(Pstream::nProcs(), 0);
forAll(sendNonLocal, procI) forAll(sendNonLocal, procI)
{ {
sendNonLocal[procI].setSize(nonLocal[procI].size()); sendNonLocal[procI] = nonLocal[procI].toc();
forAllConstIter(labelHashSet, nonLocal[procI], iter)
{
sendNonLocal[procI][nNonLocal[procI]++] = iter.key();
}
} }
if (debug) if (debug)
{ {
forAll(nNonLocal, procI) forAll(sendNonLocal, procI)
{ {
Pout<< " from processor " << procI Pout<< " from processor " << procI
<< " want " << nNonLocal[procI] << " want " << sendNonLocal[procI].size()
<< " region numbers." << " region numbers."
<< endl; << endl;
} }
@ -627,7 +343,7 @@ Foam::autoPtr<Foam::globalIndex> Foam::regionSplit::calcRegionSplit
// Get the wanted region labels into recvNonLocal // Get the wanted region labels into recvNonLocal
labelListList recvNonLocal; labelListList recvNonLocal(Pstream::nProcs());
labelListList sizes; labelListList sizes;
Pstream::exchange<labelList, label> Pstream::exchange<labelList, label>
( (
@ -655,6 +371,7 @@ Foam::autoPtr<Foam::globalIndex> Foam::regionSplit::calcRegionSplit
// Send back (into recvNonLocal) // Send back (into recvNonLocal)
recvNonLocal.clear(); recvNonLocal.clear();
recvNonLocal.setSize(sendWantedLocal.size());
sizes.clear(); sizes.clear();
Pstream::exchange<labelList, label> Pstream::exchange<labelList, label>
( (

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -89,6 +89,9 @@ Description
Can optionally keep all regions local to the processor. Can optionally keep all regions local to the processor.
Note: does not walk across cyclicAMI/cyclicACMI - since not 'coupled()'
at the patch level.
SourceFiles SourceFiles
regionSplit.C regionSplit.C
@ -126,30 +129,10 @@ class regionSplit
// Private Member Functions // Private Member Functions
//- Transfer faceRegion data from one face to the other (or vice versa) //- Calculate region split in non-compact (global) numbering.
void transferCoupledFaceRegion void calcNonCompactRegionSplit
(
const label faceI,
const label otherFaceI,
labelList& faceRegion,
DynamicList<label>& newChangedFaces
) const;
//- Given a seed cell label, fill cellRegion/faceRegion with markValue
// for contiguous region around it
void fillSeedMask
(
const List<labelPair>& explicitConnections,
labelList& cellRegion,
labelList& faceRegion,
const label seedCellID,
const label markValue
) const;
//- Calculate local region split. Return number of regions.
label calcLocalRegionSplit
( (
const globalIndex& globalFaces,
const boolList& blockedFace, const boolList& blockedFace,
const List<labelPair>& explicitConnections, const List<labelPair>& explicitConnections,
labelList& cellRegion labelList& cellRegion