BUG: Fixes bug in plicRDF interface reconstruction across cyclic boundaries (fixes #3279)

Description:
The plicRDF interface reconstruction scheme calculates a reconstructed distance
function (RDF) in all interface cells and their point neighbours.
In point neighbours to interface cells, the RDF is calculated as a weighted
average of the distances to all adjacent interface cell centres with the weight
being the inverse distance to the cell centre.
By using the zoneDistribute class written by Henning Scheufler, the required
stencil data is communicated efficiently for stencil cells living on different
sides of one or more processor patches.
Some of the data required for the RDF reconstruction are cell centre and
interface centre positions. When a stencil extends across a cyclic patch these
positions have so far not been properly transformed in OpenFOAM. This issue is
fixed by the current contribution.
The fix is done by modifying the zoneDistribute class to hold the required
information about zone cells adjacent to cyclic patches. Positions are then
communicated with a new getPosition function which replaces getValue for
position data in the reconstructedDistanceFunction and plicRDF classes.
The implementation does not change the behaviour for cells not on
a cyclic patch and should have insignificant effect on efficiency for these.
The implementation can probably be optimised in terms of efficiency for zone
cells on cyclic patches, but we note that there are generally only very few of
these (interface cells and their point neighbours on cyclic patches) and so
the potential for speedup is expected to be limited.

Current limitations:

- In parallel, the user must constrain the decomposition to preserve cyclic
patches on the same processor, for the implementation to work properly.
  - See an example here: tutorials/discInConstantFlowCyclicBCs/system/decomposeParDict
- In the case of parallel rotational cyclics that are split by the decomposition
the current bugfix does not work and therefore throws an error. This is ongoing
work and should be reported and fixed by a future patch.

For further details, please see the modified files and the comments therein:
- $FOAM_SRC/finiteVolume/fvMesh/zoneDistribute/zoneDistribute.H
- $FOAM_SRC/finiteVolume/fvMesh/zoneDistribute/zoneDistribute.C
- $FOAM_SRC/finiteVolume/fvMesh/zoneDistribute/zoneDistributeI.H
- $FOAM_SRC/transportModels/geometricVoF/reconstructedDistanceFunction/reconstructedDistanceFunction.C
- $FOAM_SRC/src/transportModels/geometricVoF/reconstructionSchemes/plicSchemes/plicRDF/plicRDF.C

Co-authored-by: David Müller <> KIT
Co-authored-by: Konstantinos Missios <> Roskilde University
Co-authored-by: Johan Roenby <> Roskilde University and STROMNING
This commit is contained in:
Johan Roenby 2025-01-10 15:32:29 +00:00 committed by Kutalmis Bercin
parent cb1aa273fd
commit 4fe3f55e4d
10 changed files with 424 additions and 18 deletions

View File

@ -27,6 +27,7 @@ License
\*---------------------------------------------------------------------------*/
#include "zoneDistribute.H"
#include "processorPolyPatch.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -44,10 +45,48 @@ Foam::zoneDistribute::zoneDistribute(const fvMesh& mesh)
stencil_(zoneCPCStencil::New(mesh)),
globalNumbering_(stencil_.globalNumbering()),
send_(UPstream::nProcs()),
pBufs_(UPstream::commsTypes::nonBlocking)
pBufs_(UPstream::commsTypes::nonBlocking),
cyclicBoundaryCells_(mesh.nCells(), false)
{
// Don't clear storage on persistent buffer
pBufs_.allowClearRecv(false);
// Loop over boundary patches and store cells with a face on a cyclic patch
bool hasCyclicPatches = false;
forAll(mesh.boundaryMesh(), patchi)
{
const cyclicPolyPatch* cpp =
isA<cyclicPolyPatch>(mesh.boundaryMesh()[patchi]);
if (cpp)
{
cyclicBoundaryCells_.set(cpp->faceCells());
hasCyclicPatches = true;
}
}
// Populate cyclicCentres_
if(hasCyclicPatches)
{
// Make a boolList from the bitSet
boolList isCyclicCell(mesh.nCells(), false);
forAll(cyclicBoundaryCells_, celli)
{
if (cyclicBoundaryCells_.test(celli))
{
isCyclicCell[celli] = true;
}
}
// Use getFields to get map of cell centres across processor boundaries
setUpCommforZone(isCyclicCell, true);
cyclicCentres_.reset
(
new Map<vectorField>(getFields(isCyclicCell, mesh_.C()))
);
}
}
@ -141,5 +180,226 @@ void Foam::zoneDistribute::setUpCommforZone
}
}
Foam::List<Foam::label> Foam::zoneDistribute::getCyclicPatches
(
const label celli,
const label globalIdx
) const
{
// Initialise cyclic patch label list
List<label> patches(0);
// If celli is not on a cyclic boundary, return the empty list
if (!cyclicBoundaryCells_.test(celli))
{
return patches;
}
const polyBoundaryMesh& bMesh = mesh_.boundaryMesh();
// So celli belongs to at least one cyclic patch.
// Let us figure out which.
if (globalNumbering_.isLocal(globalIdx)) // celli and globalIdx on same proc
{
// Making list of cyclic patches to which celli belongs
List<label> celliCyclicPatches;
forAll(bMesh, patchi)
{
if (isA<cyclicPolyPatch>(bMesh[patchi]))
{
// Note: Probably not efficient due to use of found(celli) but
// typically only used for very cells (interface cells and their
// point neighbours on cyclic boundaries).
if (bMesh[patchi].faceCells().found(celli))
{
celliCyclicPatches.append(patchi);
}
}
}
// Get all local point neighbor cells of celli, i.e. all point
// neighbours that are not on the other side of a cyclic patch.
List<label> localPointNeiCells(0);
const labelList& cellPoints = mesh_.cellPoints()[celli];
for (const label cellPoint : cellPoints)
{
const labelList& pointKCells = mesh_.pointCells()[cellPoint];
for (const label pointKCell : pointKCells)
{
if (!localPointNeiCells.found(pointKCell))
{
localPointNeiCells.append(pointKCell);
}
}
}
// Since globalIdx is a global cell index obtained from the point
// neighbour list, stencil[celli], all cells in this that are not in
// localPointNeiCells must be cyclic neighbour cells.
const label localIdx = globalNumbering_.toLocal(globalIdx);
if (!localPointNeiCells.found(localIdx))
{
for (const label patchi : celliCyclicPatches)
{
// find the corresponding cyclic neighbor patch ID
const cyclicPolyPatch& cpp =
static_cast<const cyclicPolyPatch&>(bMesh[patchi]);
const label neiPatch = cpp.neighbPatchID();
// Check if the cell globalIdx is on neiPatch.
// If it is, append neiPatch to list of
if (bMesh[neiPatch].faceCells().found(localIdx))
{
patches.append(neiPatch);
// Here it may be possible to append patchi and do:
//
// cpp.transformPosition()
//
// instead of
//
// cpp.neighbPatch().transformPosition() in getPosition()
}
}
}
}
else // celli and globalIdx on differet processors
{
// Note: The following is needed if a celli is located at the interface
// (plicRDF), on a cyclic patch and a processor patch. In this case
// globalIdx may be on a different processor, but requires the
// transformation from a cyclic patch on the processor of celli.
const List<label>& faces = mesh_.cells()[celli];
// Loop over all faces of celli and find cyclic patches
for (const label facei : faces)
{
if (mesh_.isInternalFace(facei)) continue;
const label patchi = bMesh.whichPatch(facei);
const cyclicPolyPatch* cpp = isA<cyclicPolyPatch>(bMesh[patchi]);
if (cpp)
{
// Get the neighbor cell across the cyclic face
const label cycNeiPatch = cpp->neighbPatchID();
const label cycNeiCell =
bMesh[cycNeiPatch].faceCells()[facei - cpp->start()];
const List<label>& cycNeiCellFaces = mesh_.cells()[cycNeiCell];
// Loop over all the faces of the neighbor cell
for (const label cycNeiCellFace : cycNeiCellFaces)
{
if (mesh_.isInternalFace(cycNeiCellFace))
{
continue;
}
// Check if the neighbor cell has processor patches
const label neiPatch = bMesh.whichPatch(cycNeiCellFace);
const processorPolyPatch* ppp =
isA<processorPolyPatch>(bMesh[neiPatch]);
if (ppp)
{
// Avoid duplicate entries
if (patches.found(cycNeiPatch))
{
continue;
}
// Since we can not access any information on globalIdx
// we use the cell centre map from the stencil to
// identify the cell.
const label localFaceID = cycNeiCellFace - ppp->start();
const vector neiCentre =
ppp->neighbFaceCellCentres()[localFaceID];
forAll(cyclicCentres_()[celli], k)
{
if
(
(
mag(cyclicCentres_()[celli][k] - neiCentre)
< 100*SMALL
)
&& (stencil_[celli][k] == globalIdx)
)
{
patches.append(cycNeiPatch);
// Here an alternative might be to append patchi
// and do:
// cpp.transformPosition()
// instead of
// cpp.neighbPatch().transformPosition()
// in getValue()
}
}
}
}
}
}
}
return patches;
}
Foam::vector Foam::zoneDistribute::getPosition
(
const VolumeField<vector>& positions,
const Map<vector>& valuesFromOtherProc,
const label gblIdx,
const List<label> cyclicPatchID
) const
{
// Position vector, possibly from other processor, to be returned
vector position(getValue(positions, valuesFromOtherProc, gblIdx));
// Dealing with position transformation across cyclic patches.
// If no transformation is required (most cases), cyclicPatchID is empty
forAll(cyclicPatchID, i)
{
const label patchi = cyclicPatchID[i];
const cyclicPolyPatch& cpp =
static_cast<const cyclicPolyPatch&>
(
positions.mesh().boundaryMesh()[patchi]
);
if (cpp.transform() != coupledPolyPatch::transformType::ROTATIONAL)
{
cpp.neighbPatch().transformPosition(position, 0);
}
else if (globalNumbering_.isLocal(gblIdx))
{
const label localIdx = globalNumbering_.toLocal(gblIdx);
for (const label facei : mesh_.cells()[localIdx])
{
if (mesh_.boundaryMesh().whichPatch(facei) == cyclicPatchID[i])
{
cpp.neighbPatch().transformPosition(position, facei);
continue;
}
}
}
else
{
FatalErrorInFunction
<< "Rotational cyclic patches are not supported in parallel.\n"
<< "Try to decompose the domain so that the rotational cyclic"
<< "patch is not split in between processors."
<< exit(FatalError);
}
}
return position;
}
// ************************************************************************* //

View File

@ -115,6 +115,15 @@ class zoneDistribute
//- Persistent set of exchange buffers
PstreamBuffers pBufs_;
//- Cell labels for all cells with a face on a cyclic boundary
bitSet cyclicBoundaryCells_;
//- Holds for each cell on cyclic patch the centres of the cells in its
// point neighbour stencil.
// Used in getCyclicPatches to identify neigbhour patch ID of point
// neighbours on other side of a processorPolyPatch.
autoPtr<Map<vectorField>> cyclicCentres_;
// Private Member Functions
@ -175,6 +184,19 @@ public:
return globalNumbering_;
}
//- Finds and returns list of all cyclic patch labels to which celli's
// point neighbour cell, globalIdx, belongs. celli and globalIdx touch
// in at least one point on these patces. globalIdx typically belongs
// to stencil_[celli]. The returned label list is used to transform
// positions across cyclic boundaries e.g. to be able to calculate
// distances between cell centres and interface centres in plicRDF
// across such boundaries.
List<label> getCyclicPatches
(
const label celli,
const label globalIdx
) const;
//- Gives patchNumber and patchFaceNumber for a given
//- Geometric volume field
template<typename Type>
@ -185,6 +207,14 @@ public:
const label gblIdx
) const;
vector getPosition
(
const VolumeField<vector>& positions,
const Map<vector>& valuesFromOtherProc,
const label gblIdx,
const List<label> cyclicPatchID = List<label>()
) const;
//- Returns stencil and provides a Map with globalNumbering
template<typename Type>
Map<Field<Type>> getFields
@ -193,6 +223,17 @@ public:
const VolumeField<Type>& phi
);
//- Returns stencil and provides a Map with globalNumbering
// Note: Not used currently (v2412) but needed for future surface
// tension modelling.
template<typename Type>
Map<Field<Type>> getPositionFields
(
const boolList& zone,
const VolumeField<Type>& phi,
const bool& checkTransformation = false
);
//- Returns stencil and provides a Map with globalNumbering
template<typename Type>
Map<Type> getDatafromOtherProc

View File

@ -128,7 +128,59 @@ Foam::Map<Foam::Field<Type>> Foam::zoneDistribute::getFields
for (const label gblIdx : stencil_[celli])
{
tmpField.append(getValue(phi,neiValues,gblIdx));
tmpField.append(getValue(phi, neiValues, gblIdx));
}
stencilWithValues.emplace(celli, tmpField);
}
}
return stencilWithValues;
}
template<typename Type>
Foam::Map<Foam::Field<Type>> Foam::zoneDistribute::getPositionFields
(
const boolList& zone,
const VolumeField<Type>& phi,
const bool& checkTransformation
)
{
if (zone.size() != phi.size())
{
FatalErrorInFunction
<< "size of zone: " << zone.size()
<< "size of phi:" << phi.size()
<< "do not match. Did the mesh change?"
<< exit(FatalError);
}
// Get values from other proc
Map<Type> neiValues = getDatafromOtherProc(zone, phi);
Map<Field<Type>> stencilWithValues;
DynamicField<Type> tmpField(128);
forAll(zone, celli)
{
if (zone[celli])
{
tmpField.clear();
for (const label gblIdx : stencil_[celli])
{
List<label> cyclicPatches(0);
if(checkTransformation)
{
cyclicPatches = getCyclicPatches(celli, gblIdx);
}
tmpField.append
(
getPosition(phi, neiValues, gblIdx, cyclicPatches)
);
}
stencilWithValues.emplace(celli, tmpField);

View File

@ -282,9 +282,18 @@ const Foam::volScalarField& Foam::reconstructedDistanceFunction::constructRDF
if (mag(n) != 0)
{
n /= mag(n);
vector c = distribute.getValue(centre,mapCentres,gblIdx);
vector distanceToIntSeg = (c - p);
scalar distToSurf = distanceToIntSeg & (n);
vector c
(
distribute.getPosition
(
centre,
mapCentres,
gblIdx,
distribute.getCyclicPatches(celli, gblIdx)
)
);
vector distanceToIntSeg(c - p);
scalar distToSurf = distanceToIntSeg & n;
scalar weight = 0;
if (mag(distanceToIntSeg) != 0)
@ -336,10 +345,18 @@ const Foam::volScalarField& Foam::reconstructedDistanceFunction::constructRDF
if (mag(n) != 0)
{
n /= mag(n);
vector c =
distribute.getValue(centre, mapCentres, gblIdx);
vector distanceToIntSeg = (c - p);
scalar distToSurf = distanceToIntSeg & (n);
vector c
(
distribute.getPosition
(
centre,
mapCentres,
gblIdx,
distribute.getCyclicPatches(pCellI, gblIdx)
)
);
vector distanceToIntSeg(c - p);
scalar distToSurf = distanceToIntSeg & n;
scalar weight = 0;
if (mag(distanceToIntSeg) != 0)

View File

@ -96,8 +96,16 @@ void Foam::reconstruction::plicRDF::interpolateNormal()
if (mag(n) != 0)
{
n /= mag(n);
vector centre =
exchangeFields.getValue(centre_, mapCentre, gblIdx);
vector centre
(
exchangeFields.getPosition
(
centre_,
mapCentre,
gblIdx,
exchangeFields.getCyclicPatches(celli, gblIdx)
)
);
vector distanceToIntSeg = (tensor::I- n*n) & (p - centre);
estimatedNormal += n /max(mag(distanceToIntSeg), SMALL);
weight += 1/max(mag(distanceToIntSeg), SMALL);
@ -146,7 +154,11 @@ void Foam::reconstruction::plicRDF::interpolateNormal()
const label gblIdx = stencil[celli][i];
cellCentre.append
(
exchangeFields.getValue(mesh_.C(), mapCC, gblIdx)
exchangeFields.getPosition
(
mesh_.C(), mapCC, gblIdx,
exchangeFields.getCyclicPatches(celli, gblIdx)
)
);
alphaValues.append
(
@ -193,7 +205,11 @@ void Foam::reconstruction::plicRDF::gradSurf(const volScalarField& phi)
{
cellCentre.append
(
exchangeFields.getValue(mesh_.C(), mapCC, gblIdx)
exchangeFields.getPosition
(
mesh_.C(), mapCC, gblIdx,
exchangeFields.getCyclicPatches(celli, gblIdx)
)
);
phiValues.append
(

View File

@ -5,4 +5,6 @@ cd "${0%/*}" || exit # Run from this directory
cleanCase0
rm -f Phase1
#------------------------------------------------------------------------------

View File

@ -5,10 +5,20 @@ cd "${0%/*}" || exit # Run from this directory
restore0Dir
touch case.foam
runApplication blockMesh
runApplication setAlphaField
runApplication $(getApplication)
runApplication decomposePar
runParallel $(getApplication)
grep "Phase-1" log.$(getApplication) | cut -d' ' -f5 > Phase1
echo "Volume conservation:"
echo "Initial Phase-1: " $(head -1 Phase1)
echo "Final Phase-1 : " $(tail -1 Phase1)
#------------------------------------------------------------------------------

View File

@ -35,7 +35,7 @@ purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writePrecision 12;
writeCompression off;

View File

@ -24,5 +24,13 @@ coeffs
n ( 2 1 2 );
}
constraints
{
preservePatches
{
type preservePatches;
patches (left right top bottom);
}
}
// ************************************************************************* //

View File

@ -22,10 +22,10 @@ solvers
{
isoFaceTol 1e-6;
surfCellTol 1e-6;
snapTol 1e-10;
snapTol 0;
nAlphaBounds 3;
clip true;
reconstructionScheme isoAlpha; // isoAlpha
clip false;
reconstructionScheme plicRDF;
nAlphaSubCycles 1;
cAlpha 1; // Read but not used by interIsoFoam