Commit 80d28942 authored by mattijs's avatar mattijs Committed by Andrew Heather
Browse files

ENH: overset: GAMG normalisation of hole cells

parent 9cd9d879
......@@ -357,6 +357,71 @@ bool Foam::dynamicOversetFvMesh::updateAddressing() const
}
Foam::scalar Foam::dynamicOversetFvMesh::cellAverage
(
const labelList& types,
const labelList& nbrTypes,
const scalarField& norm,
const scalarField& nbrNorm,
const label celli,
bitSet& isFront
) const
{
const labelList& own = faceOwner();
const labelList& nei = faceNeighbour();
const cell& cFaces = cells()[celli];
scalar avg = 0.0;
label n = 0;
label nFront = 0;
for (const label facei : cFaces)
{
if (isInternalFace(facei))
{
label nbrCelli = (own[facei] == celli ? nei[facei] : own[facei]);
if (norm[nbrCelli] == -GREAT)
{
// Invalid neighbour. Add to front
if (isFront.set(facei))
{
nFront++;
}
}
else
{
// Valid neighbour. Add to average
avg += norm[nbrCelli];
n++;
}
}
else
{
if (nbrNorm[facei-nInternalFaces()] == -GREAT)
{
if (isFront.set(facei))
{
nFront++;
}
}
else
{
avg += nbrNorm[facei-nInternalFaces()];
n++;
}
}
}
if (n > 0)
{
return avg/n;
}
else
{
return norm[celli];
}
}
void Foam::dynamicOversetFvMesh::writeAgglomeration
(
const GAMGAgglomeration& agglom
......@@ -375,7 +440,8 @@ void Foam::dynamicOversetFvMesh::writeAgglomeration
this->time().timeName(),
*this,
IOobject::NO_READ,
IOobject::AUTO_WRITE
IOobject::NO_WRITE,
false
),
*this,
dimensionedScalar(dimless, Zero)
......@@ -430,7 +496,8 @@ void Foam::dynamicOversetFvMesh::writeAgglomeration
this->time().timeName(),
*this,
IOobject::NO_READ,
IOobject::AUTO_WRITE
IOobject::NO_WRITE,
false
),
*this,
dimensionedScalar(dimless, Zero)
......
......@@ -146,6 +146,17 @@ protected:
template<class GeoField>
static void correctCoupledBoundaryConditions(GeoField& fld);
//- Average norm of valid neighbours
scalar cellAverage
(
const labelList& types,
const labelList& nbrTypes,
const scalarField& norm,
const scalarField& nbrNorm,
const label celli,
bitSet& isFront
) const;
//- Debug: dump agglomeration
void writeAgglomeration
(
......
......@@ -30,6 +30,7 @@ License
#include "calculatedProcessorFvPatchField.H"
#include "lduInterfaceFieldPtrsList.H"
#include "processorFvPatch.H"
#include "syncTools.H"
// * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * * //
......@@ -171,19 +172,152 @@ Foam::tmp<Foam::scalarField> Foam::dynamicOversetFvMesh::normalisation
}
}
forAll(norm, celli)
// Count number of problematic cells
label nZeroDiag = 0;
for (const scalar n : norm)
{
scalar& n = norm[celli];
if (mag(n) < SMALL)
if (magSqr(n) < sqr(SMALL))
{
n = 1.0; //?
nZeroDiag++;
}
else
}
reduce(nZeroDiag, sumOp<label>());
if (debug)
{
Pout<< "For field " << m.psi().name() << " have zero diagonals for "
<< nZeroDiag << " cells" << endl;
}
if (nZeroDiag > 0)
{
// Walk out the norm across hole cells
const labelList& own = faceOwner();
const labelList& nei = faceNeighbour();
const cellCellStencilObject& overlap = Stencil::New(*this);
const labelUList& types = overlap.cellTypes();
label nHoles = 0;
scalarField extrapolatedNorm(norm);
forAll(types, celli)
{
if (types[celli] == cellCellStencil::HOLE)
{
extrapolatedNorm[celli] = -GREAT;
nHoles++;
}
}
bitSet isFront(nFaces());
for (label facei = 0; facei < nInternalFaces(); facei++)
{
label ownType = types[own[facei]];
label neiType = types[nei[facei]];
if
(
(ownType == cellCellStencil::HOLE)
!= (neiType == cellCellStencil::HOLE)
)
{
isFront.set(facei);
}
}
labelList nbrTypes;
syncTools::swapBoundaryCellList(*this, types, nbrTypes);
for (label facei = nInternalFaces(); facei < nFaces(); facei++)
{
label ownType = types[own[facei]];
label neiType = nbrTypes[facei-nInternalFaces()];
if
(
(ownType == cellCellStencil::HOLE)
!= (neiType == cellCellStencil::HOLE)
)
{
isFront.set(facei);
}
}
while (true)
{
scalarField nbrNorm;
syncTools::swapBoundaryCellList(*this, extrapolatedNorm, nbrNorm);
bitSet newIsFront(nFaces());
scalarField newNorm(extrapolatedNorm);
label nChanged = 0;
for (const label facei : isFront)
{
if (extrapolatedNorm[own[facei]] == -GREAT)
{
// Average owner cell, add faces to newFront
newNorm[own[facei]] = cellAverage
(
types,
nbrTypes,
extrapolatedNorm,
nbrNorm,
own[facei],
newIsFront
);
nChanged++;
}
if
(
isInternalFace(facei)
&& extrapolatedNorm[nei[facei]] == -GREAT
)
{
// Average nei cell, add faces to newFront
newNorm[nei[facei]] = cellAverage
(
types,
nbrTypes,
extrapolatedNorm,
nbrNorm,
nei[facei],
newIsFront
);
nChanged++;
}
}
reduce(nChanged, sumOp<label>());
if (nChanged == 0)
{
break;
}
// Transfer new front
extrapolatedNorm.transfer(newNorm);
isFront.transfer(newIsFront);
syncTools::syncFaceList(*this, isFront, maxEqOp<unsigned int>());
}
forAll(norm, celli)
{
// Restore original diagonal
n = m.diag()[celli];
scalar& n = norm[celli];
if (mag(n) < SMALL)
{
n = extrapolatedNorm[celli];
}
else
{
// Use original diagonal
n = m.diag()[celli];
}
}
}
else
{
// Use original diagonal
norm = m.diag();
}
return tnorm;
}
......@@ -545,6 +679,38 @@ Foam::SolverPerformance<Type> Foam::dynamicOversetFvMesh::solve
// Calculate stabilised diagonal as normalisation for interpolation
const scalarField norm(normalisation(m));
if (debug)
{
volScalarField scale
(
IOobject
(
m.psi().name() + "_normalisation",
this->time().timeName(),
*this,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
*this,
dimensionedScalar(dimless, Zero)
);
scale.ref().field() = norm;
correctBoundaryConditions
<
volScalarField,
oversetFvPatchField<scalar>
>(scale.boundaryFieldRef(), false);
scale.write();
if (debug)
{
Pout<< "dynamicOversetFvMesh::solve() :"
<< " writing matrix normalisation for field " << m.psi().name()
<< " to " << scale.name() << endl;
}
}
// Switch to extended addressing (requires mesh::update() having been
// called)
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment