GGI conservation for partially covered faces

Changes necessary for correctly creating weights, deltas and deltaCoeffs for
bridged overlap.
This commit is contained in:
Vuko Vukcevic 2018-02-21 14:41:29 +01:00
parent b495da0d2c
commit 88baefbf7d
8 changed files with 204 additions and 91 deletions

View file

@ -188,10 +188,19 @@ tmp<Field<Type> > ggiFvPatchField<Type>::patchNeighbourField() const
this->patchInternalField()
);
// Note: bridging now takes into account fully uncovered and partially
// covered faces. VV, 18/Oct/2017.
// if (pTraits<Type>::rank == 0)
// {
// // Scale the field for scalars to ensure conservative and consistent
// // flux on both sides
// ggiPatch_.scaleForPartialCoverage(bridgeField, pnf);
// }
// else
{
// Bridge the field for higher order tensors to correctly take into
// account mirroring
ggiPatch_.bridge(bridgeField, pnf);
}
}
return tpnf;
}

View file

@ -62,8 +62,9 @@ void Foam::ggiFvPatch::makeWeights(scalarField& w) const
// HJ, 2/Aug/2007
if (ggiPolyPatch_.master())
{
// Master side. No need to bridge the overlap since reconFaceCellCentres
// already takes it into account. VV, 14/Feb/2018
// Master side. No need to scale partially uncovered or set fully
// uncovered faces since delta already takes it into account.
// VV, 25/Feb/2018.
const vectorField n = nf();
@ -79,8 +80,9 @@ void Foam::ggiFvPatch::makeWeights(scalarField& w) const
}
else
{
// Slave side. Interpolate the master side weights and scale them up for
// partially covered faces if the bridge overlap is switched on
// Slave side. Interpolate the master side weights, scale them for
// partially covered faces and set weights for fully uncovered faces if
// the bridge overlap is switched on. VV, 15/Feb/2018.
scalarField masterWeights(shadow().size());
shadow().makeWeights(masterWeights);
@ -90,13 +92,12 @@ void Foam::ggiFvPatch::makeWeights(scalarField& w) const
if (bridgeOverlap())
{
// Bridge weights for symmetry plane treatment
const scalarField bridgeField(w.size(), 0.5);
// Weights for fully uncovered faces
const scalarField uncoveredWeights(w.size(), 0.5);
bridge(bridgeField, w);
// Scale master weights for partially overlapping faces
// scaleForPartialCoverage(w);
// Scale partially overlapping faces and set uncovered weights
// for fully uncovered faces
scaleForPartialCoverage(uncoveredWeights, w);
}
// Finally construct these weights as 1 - master weights
@ -110,8 +111,9 @@ void Foam::ggiFvPatch::makeDeltaCoeffs(scalarField& dc) const
{
if (ggiPolyPatch_.master())
{
// Master side. No need to bridge the overlap since delta already takes
// it into account. VV, 18/Oct/2017.
// Master side. No need to scale partially uncovered or set fully
// uncovered faces since delta already takes it into account.
// VV, 25/Feb/2018.
// Stabilised form for bad meshes. HJ, 24/Aug/2011
const vectorField d = delta();
@ -120,8 +122,9 @@ void Foam::ggiFvPatch::makeDeltaCoeffs(scalarField& dc) const
}
else
{
// Slave side. Interpolate the master side and scale it up for partially
// covered faces if the bridge overlap is switched on
// Slave side. Interpolate the master side, scale it for partially
// covered faces and set deltaCoeffs for fully uncovered faces if the
// bridge overlap is switched on. VV, 15/Feb/2018.
scalarField masterDeltas(shadow().size());
shadow().makeDeltaCoeffs(masterDeltas);
@ -129,7 +132,15 @@ void Foam::ggiFvPatch::makeDeltaCoeffs(scalarField& dc) const
if (bridgeOverlap())
{
scaleForPartialCoverage(dc);
// Delta coeffs for fully uncovered faces obtained from deltas on
// this side
const vectorField d = delta();
const scalarField uncoveredDeltaCoeffs =
1.0/max(nf() & d, 0.05*mag(d));
// Scale partially overlapping faces and set uncovered deltaCoeffs
// for fully uncovered faces.
scaleForPartialCoverage(uncoveredDeltaCoeffs, dc);
}
}
}
@ -158,8 +169,9 @@ Foam::tmp<Foam::vectorField> Foam::ggiFvPatch::delta() const
{
if (ggiPolyPatch_.master())
{
// Master side. Note: bridging partially covered faces correctly taken
// into account in reconFaceCellCentres function
// Master side. Note: scaling partially covered faces and setting deltas
// to fully uncovered faces correctly taken into account in
// reconFaceCellCentres function. VV, 15/Feb/2018.
tmp<vectorField> tDelta = ggiPolyPatch_.reconFaceCellCentres() - Cn();
@ -167,8 +179,9 @@ Foam::tmp<Foam::vectorField> Foam::ggiFvPatch::delta() const
}
else
{
// Slave side. Interpolate the master side and scale it up for partially
// covered faces if the bridge overlap is switched on
// Slave side. Interpolate the master side, scale it for partially
// covered faces and set deltas for fully uncovered faces if the bridge
// overlap is switched on. VV, 15/Feb/2018.
tmp<vectorField> tDelta = interpolate
(
@ -177,7 +190,12 @@ Foam::tmp<Foam::vectorField> Foam::ggiFvPatch::delta() const
if (bridgeOverlap())
{
scaleForPartialCoverage(tDelta());
// Deltas for fully uncovered faces
const vectorField uncoveredDeltas(2.0*fvPatch::delta());
// Scale partially overlapping faces and set uncovered deltas for
// fully uncovered faces
scaleForPartialCoverage(uncoveredDeltas, tDelta());
}
return tDelta;

View file

@ -154,15 +154,22 @@ public:
return ggiPolyPatch_.bridge(bridgeField, ff);
}
//- Scale field for partially covered faces. Needed for correct
//- Scale field for partially covered faces and set
// uncoveredFaceField to fully uncovered faces. Needed for correct
// calculation of fvPatch data: delta, deltaCoeffs and weights
template<class Type>
void scaleForPartialCoverage
(
Field<Type>& fieldToScale
const Field<Type>& uncoveredFaceField,
Field<Type>& ff
) const
{
return ggiPolyPatch_.scaleForPartialCoverage(fieldToScale);
return
ggiPolyPatch_.scaleForPartialCoverage
(
uncoveredFaceField,
ff
);
}

View file

@ -207,16 +207,24 @@ template<class MasterPatch, class SlavePatch>
template<class Type>
void GGIInterpolation<MasterPatch, SlavePatch>::scale
(
Field<Type>& fieldToScale,
const Field<Type>& uncoveredFaceField,
Field<Type>& ff,
const labelList& addr,
const labelList& partiallyUncoveredAddr,
const scalarField& uncoveredFractions
)
{
// Loop through uncovered faces and set them
forAll (addr, ufI)
{
const label& faceI = addr[ufI];
ff[faceI] = uncoveredFaceField[faceI];
}
// Loop through partially covered faces and scale them up
forAll (partiallyUncoveredAddr, pcfI)
{
fieldToScale[partiallyUncoveredAddr[pcfI]] /=
ff[partiallyUncoveredAddr[pcfI]] /=
(1.0 - uncoveredFractions[pcfI]);
}
}
@ -226,7 +234,8 @@ template<class MasterPatch, class SlavePatch>
template<class Type>
void GGIInterpolation<MasterPatch, SlavePatch>::maskedScale
(
Field<Type>& fieldToScale,
const Field<Type>& uncoveredFaceField,
Field<Type>& ff,
const labelList& mask,
const labelList& uncoveredFaces,
const labelList& partiallyUncoveredAddr,
@ -241,6 +250,39 @@ void GGIInterpolation<MasterPatch, SlavePatch>::maskedScale
label maskAddrI = 0;
forAll (uncoveredFaces, uncoI)
{
// Pick the uncovered face
const label faceI = uncoveredFaces[uncoI];
// Search through the mask
for (; maskAddrI < mask.size(); ++maskAddrI)
{
if (faceI == mask[maskAddrI])
{
// Found masked bridged face
// Put the result into condensed list: masked faces only
ff[maskAddrI] = uncoveredFaceField[maskAddrI];
break;
}
else if (mask[maskAddrI] > faceI)
{
// Gone beyond my index: my face is not present in the mask
// Go one back and check for next uncovered face
if (maskAddrI > 0)
{
--maskAddrI;
}
break;
}
}
}
// Reset maskAddrI
maskAddrI = 0;
forAll (partiallyUncoveredAddr, pcfI)
{
// Pick partially covered face
@ -251,7 +293,7 @@ void GGIInterpolation<MasterPatch, SlavePatch>::maskedScale
if (faceI == mask[maskAddrI])
{
// Found masked partially covered face
fieldToScale[maskAddrI] /= (1.0 - uncoveredFractions[pcfI]);
ff[maskAddrI] /= (1.0 - uncoveredFractions[pcfI]);
break;
}
@ -624,15 +666,15 @@ void GGIInterpolation<MasterPatch, SlavePatch>::bridgeMaster
{
if
(
bridgeField.size() != masterPatch_.size()
|| ff.size() != masterPatch_.size()
(bridgeField.size() != masterPatch_.size())
|| (ff.size() != masterPatch_.size())
)
{
FatalErrorIn
(
"void GGIInterpolation<MasterPatch, SlavePatch>::bridgeMaster\n"
"(\n"
" const Field<Type>& bridgeField,\n"
" const Field<Type>& bridgeField\n,"
" Field<Type>& ff\n"
") const"
) << "given field does not correspond to patch. Patch size: "
@ -669,7 +711,7 @@ void GGIInterpolation<MasterPatch, SlavePatch>::maskedBridgeMaster
"void GGIInterpolation<MasterPatch, SlavePatch>::"
"maskedBridgeMaster\n"
"(\n"
" const Field<Type>& bridgeField,\n"
" const Field<Type>& bridgeField\n,"
" Field<Type>& ff,\n"
" const labelList& mask\n"
") const"
@ -703,8 +745,8 @@ void GGIInterpolation<MasterPatch, SlavePatch>::bridgeSlave
{
if
(
bridgeField.size() != slavePatch_.size()
|| ff.size() != slavePatch_.size()
(bridgeField.size() != slavePatch_.size())
|| (ff.size() != slavePatch_.size())
)
{
FatalErrorIn
@ -712,7 +754,7 @@ void GGIInterpolation<MasterPatch, SlavePatch>::bridgeSlave
"void GGIInterpolation<MasterPatch, SlavePatch>::"
"bridgeSlave\n"
"(\n"
" const Field<Type>& bridgeField,\n"
" const Field<Type>& bridgeField\n,"
" Field<Type>& ff"
") const"
) << "given field does not correspond to patch. Patch size: "
@ -749,7 +791,7 @@ void GGIInterpolation<MasterPatch, SlavePatch>::maskedBridgeSlave
"void GGIInterpolation<MasterPatch, SlavePatch>::"
"maskedBridgeSlave\n"
"(\n"
" const Field<Type>& bridgeField,\n"
" const Field<Type>& bridgeField\n,"
" Field<Type>& ff\n,"
" const labelList& mask\n"
") const"
@ -777,26 +819,34 @@ template<class MasterPatch, class SlavePatch>
template<class Type>
void GGIInterpolation<MasterPatch, SlavePatch>::scaleMaster
(
Field<Type>& fieldToScale
const Field<Type>& uncoveredFaceField,
Field<Type>& ff
) const
{
if (fieldToScale.size() != masterPatch_.size())
if
(
uncoveredFaceField.size() != masterPatch_.size()
|| ff.size() != masterPatch_.size()
)
{
FatalErrorIn
(
"void GGIInterpolation<MasterPatch, SlavePatch>::scaleMaster\n"
"(\n"
" Field<Type>& fieldToScale\n"
" const Field<Type>& uncoveredFaceField\n,"
" Field<Type>& ff\n"
") const"
) << "given field does not correspond to patch. Patch size: "
<< masterPatch_.size()
<< " scale field size: " << fieldToScale.size()
<< " uncovered face field size: " << uncoveredFaceField.size()
<< " scale field size: " << ff.size()
<< abort(FatalError);
}
scale
(
fieldToScale,
uncoveredFaceField,
ff,
uncoveredMasterFaces(),
partiallyUncoveredMasterFaces(),
masterFaceUncoveredFractions()
@ -808,30 +858,34 @@ template<class MasterPatch, class SlavePatch>
template<class Type>
void GGIInterpolation<MasterPatch, SlavePatch>::maskedScaleMaster
(
Field<Type>& fieldToScale,
const Field<Type>& uncoveredFaceField,
Field<Type>& ff,
const labelList& mask
) const
{
if (fieldToScale.size() != mask.size())
if (ff.size() != mask.size())
{
FatalErrorIn
(
"void GGIInterpolation<MasterPatch, SlavePatch>::"
"maskedScaleMaster\n"
"(\n"
" Field<Type>& fieldToScale,\n"
" const Field<Type>& uncoveredFaceField\n,"
" Field<Type>& ff,\n"
" const labelList& mask\n"
") const"
) << "given field does not correspond to patch. Patch (mask) size: "
<< masterPatch_.size()
<< " scale field size: " << fieldToScale.size()
<< " uncovered face field size: " << uncoveredFaceField.size()
<< " scale field size: " << ff.size()
<< " mask size: " << mask.size()
<< abort(FatalError);
}
maskedScale
(
fieldToScale,
uncoveredFaceField,
ff,
mask,
uncoveredMasterFaces(),
partiallyUncoveredMasterFaces(),
@ -844,27 +898,35 @@ template<class MasterPatch, class SlavePatch>
template<class Type>
void GGIInterpolation<MasterPatch, SlavePatch>::scaleSlave
(
Field<Type>& fieldToScale
const Field<Type>& uncoveredFaceField,
Field<Type>& ff
) const
{
if (fieldToScale.size() != slavePatch_.size())
if
(
uncoveredFaceField.size() != slavePatch_.size()
|| ff.size() != slavePatch_.size()
)
{
FatalErrorIn
(
"void GGIInterpolation<MasterPatch, SlavePatch>::"
"scaleSlave\n"
"(\n"
" Field<Type>& fieldToScale"
" const Field<Type>& uncoveredFaceField\n,"
" Field<Type>& ff\n"
") const"
) << "given field does not correspond to patch. Patch size: "
<< slavePatch_.size()
<< " scale field size: " << fieldToScale.size()
<< " uncovered face field size: " << uncoveredFaceField.size()
<< " scale field size: " << ff.size()
<< abort(FatalError);
}
scale
(
fieldToScale,
uncoveredFaceField,
ff,
uncoveredSlaveFaces(),
partiallyUncoveredSlaveFaces(),
slaveFaceUncoveredFractions()
@ -876,30 +938,34 @@ template<class MasterPatch, class SlavePatch>
template<class Type>
void GGIInterpolation<MasterPatch, SlavePatch>::maskedScaleSlave
(
Field<Type>& fieldToScale,
const Field<Type>& uncoveredFaceField,
Field<Type>& ff,
const labelList& mask
) const
{
if (fieldToScale.size() != mask.size())
if (ff.size() != mask.size())
{
FatalErrorIn
(
"void GGIInterpolation<MasterPatch, SlavePatch>::"
"maskedScaleSlave\n"
"(\n"
" Field<Type>& fieldToScale\n,"
" const Field<Type>& uncoveredFaceField\n,"
" Field<Type>& ff\n,"
" const labelList& mask\n"
") const"
) << "given field does not correspond to patch. Patch (mask) size: "
<< slavePatch_.size()
<< " scale field size: " << fieldToScale.size()
<< " uncovered face field size: " << uncoveredFaceField.size()
<< " scale field size: " << ff.size()
<< " mask size: " << mask.size()
<< abort(FatalError);
}
maskedScale
(
fieldToScale,
uncoveredFaceField,
ff,
mask,
uncoveredSlaveFaces(),
partiallyUncoveredSlaveFaces(),

View file

@ -507,22 +507,25 @@ class GGIInterpolation
// Scaling operations for partially covered faces
//- Scale partially covered faces given addressing
//- Scale partially covered faces and set fully uncovered faces
// given addressing
template<class Type>
static void scale
(
Field<Type>& fieldToScale,
const Field<Type>& uncoveredFaceField,
Field<Type>& ff,
const labelList& addr,
const labelList& partiallyUncoveredAddr,
const scalarField& uncoveredFractions
);
//- Scale partially covered faces given addressing
// for masked faces only
//- Scale partially covered faces and set fully uncovered faces
// given addressing for masked faces only
template<class Type>
static void maskedScale
(
Field<Type>& fieldToScale,
const Field<Type>& uncoveredFaceField,
Field<Type>& ff,
const labelList& mask,
const labelList& uncoveredFaces,
const labelList& partiallyUncoveredAddr,
@ -718,35 +721,41 @@ public:
// Scaling operations for correction of partially uncovered faces
//- Scale partially uncovered master patch field
//- Scale partially uncovered master patch field and set fully
// uncovered face field
template<class Type>
void scaleMaster
(
Field<Type>& fieldToScale
const Field<Type>& uncoveredFaceField,
Field<Type>& ff
) const;
//- Scale partially uncovered master patch field, only for marked
// master faces
//- Scale partially uncovered master patch field and set fully
// uncovered face field, only for marked master faces
template<class Type>
void maskedScaleMaster
(
Field<Type>& fieldToScale,
const Field<Type>& uncoveredFaceField,
Field<Type>& ff,
const labelList& mask
) const;
//- Scale partially uncovered slave patch field
//- Scale partially uncovered slave patch field and set fully
// uncovered face field
template<class Type>
void scaleSlave
(
Field<Type>& fieldToScale
const Field<Type>& uncoveredFaceField,
Field<Type>& ff
) const;
//- Scale partially uncovered slave patch field, only for marked
// slave faces
//- Scale partially uncovered slave patch field and set fully
// uncovered, only for marked slave faces
template<class Type>
void maskedScaleSlave
(
Field<Type>& fieldToScale,
const Field<Type>& uncoveredFaceField,
Field<Type>& ff,
const labelList& mask
) const;

View file

@ -294,8 +294,6 @@ void Foam::ggiPolyPatch::calcReconFaceCellCentres() const
- boundaryMesh()[shadowID].faceCentres()
);
// The field is now interpolated, but we need to bridge it as well for
// possible partially covered faces
if (bridgeOverlap_)
{
// Get necessary mesh data from polyPatch on this (master) side
@ -305,13 +303,12 @@ void Foam::ggiPolyPatch::calcReconFaceCellCentres() const
const vectorField ccf = faceCellCentres();
const vectorField nf = Sf/mag(Sf);
// Mirrored field since bridging assumes symmetry plane treatment.
// Mirror centre-to-face cell centre vectors to other side
const vectorField mirrorField =
transform(I - 2.0*sqr(nf), ccf - cf);
// Deltas for fully uncovered faces
const vectorField uncoveredDeltas(2.0*(cf - ccf));
// Take into account fully uncovered and partially covered faces
bridge(mirrorField, tdf());
// Scale partially overlapping faces and set uncovered deltas to
// fully uncovered faces
scaleForPartialCoverage(uncoveredDeltas, tdf());
}
// Calculate the reconstructed cell centres

View file

@ -378,12 +378,14 @@ public:
Field<Type>& ff
) const;
//- Scale field for partially covered faces. Needed for correct
//- Scale field for partially covered faces and set
// uncoveredFaceField to fully uncovered faces. Needed for correct
// calculation of fvPatchData: delta, deltaCoeffs and weights
template<class Type>
void scaleForPartialCoverage
(
Field<Type>& fieldToScale
const Field<Type>& uncoveredFaceField,
Field<Type>& ff
) const;

View file

@ -270,20 +270,23 @@ void Foam::ggiPolyPatch::bridge
template<class Type>
void Foam::ggiPolyPatch::scaleForPartialCoverage
(
Field<Type>& fieldToScale
const Field<Type>& uncoveredFaceField,
Field<Type>& ff
) const
{
// Check and expand the field from patch size to zone size
if (fieldToScale.size() != size())
if (ff.size() != size())
{
FatalErrorIn
(
"template<class Type> ggiPolyPatch::scaleForPartialCoverage\n"
"(\n"
" const Field<Type>& bridgeField,\n"
" const Field<Type>& uncoveredFaceField,\n"
" Field<Type>& ff,\n"
") const"
) << "Incorrect patch field size for scaling. Field size: "
<< fieldToScale.size() << " patch size: " << size()
<< ff.size() << " uncovered field size: "
<< uncoveredFaceField.size() << " patch size: " << size()
<< abort(FatalError);
}
@ -299,11 +302,11 @@ void Foam::ggiPolyPatch::scaleForPartialCoverage
{
if (master())
{
patchToPatch().scaleMaster(fieldToScale);
patchToPatch().scaleMaster(uncoveredFaceField, ff);
}
else
{
patchToPatch().scaleSlave(fieldToScale);
patchToPatch().scaleSlave(uncoveredFaceField, ff);
}
}
else
@ -313,7 +316,8 @@ void Foam::ggiPolyPatch::scaleForPartialCoverage
{
patchToPatch().maskedScaleMaster
(
fieldToScale,
uncoveredFaceField,
ff,
zoneAddressing()
);
}
@ -321,7 +325,8 @@ void Foam::ggiPolyPatch::scaleForPartialCoverage
{
patchToPatch().maskedScaleSlave
(
fieldToScale,
uncoveredFaceField,
ff,
zoneAddressing()
);
}