From a0f1e98d247419ee74af1969f22635da85b12940 Mon Sep 17 00:00:00 2001
From: Kutalmis Bercin <kutalmis.bercin@esi-group.com>
Date: Fri, 2 Dec 2022 15:49:20 +0000
Subject: [PATCH] ENH: finiteArea: improve robustness in code sections
vulnerable to math errors
It has been observed that the finite-area framework is prone to numerical
issues when zero-valued edge lenghts, edge/face normals and face areas exist.
To improve exception handling at identified code sections to gracefully
overcome math errors, the problematic entities are lower-bounded by SMALL.
---
.../liquidFilmFoam/capillaryCourantNo.H | 2 +-
.../faMesh/faMeshDemandDrivenData.C | 81 +++++++++++++++++++
.../faMesh/faPatches/faPatch/faPatch.C | 14 +++-
.../leastSquaresFaVectors.C | 14 ++++
.../edgeInterpolation/edgeInterpolation.C | 81 ++++++++++---------
.../schemes/NVDscheme/faNVDscheme.C | 12 +++
6 files changed, 166 insertions(+), 38 deletions(-)
diff --git a/applications/solvers/finiteArea/liquidFilmFoam/capillaryCourantNo.H b/applications/solvers/finiteArea/liquidFilmFoam/capillaryCourantNo.H
index 8e66fd171f2..f4e87630e02 100644
--- a/applications/solvers/finiteArea/liquidFilmFoam/capillaryCourantNo.H
+++ b/applications/solvers/finiteArea/liquidFilmFoam/capillaryCourantNo.H
@@ -4,7 +4,7 @@
sqrt
(
2*M_PI*sigma*sqr(aMesh.edgeInterpolation::deltaCoeffs())
- *aMesh.edgeInterpolation::deltaCoeffs()
+ *mag(aMesh.edgeInterpolation::deltaCoeffs())
/rhol
)
).value()*runTime.deltaT().value();
diff --git a/src/finiteArea/faMesh/faMeshDemandDrivenData.C b/src/finiteArea/faMesh/faMeshDemandDrivenData.C
index 1d98ea822d7..0a7620fc8c1 100644
--- a/src/finiteArea/faMesh/faMeshDemandDrivenData.C
+++ b/src/finiteArea/faMesh/faMeshDemandDrivenData.C
@@ -594,6 +594,12 @@ Foam::tmp<Foam::vectorField> Foam::faMesh::calcRawEdgeNormals(int order) const
edgeNormals[edgei].removeCollinear(edgeLine.unitVec());
edgeNormals[edgei].normalise();
+
+ // Do not allow any mag(val) < SMALL
+ if (mag(edgeNormals[edgei]) < SMALL)
+ {
+ edgeNormals[edgei] = vector::uniform(SMALL);
+ }
}
return tedgeNormals;
@@ -658,6 +664,12 @@ void Foam::faMesh::calcLe() const
edges_[edgei].line(localPoints),
edgeNormals[edgei]
);
+
+ // Do not allow any mag(val) < SMALL
+ if (mag(leVectors[edgei]) < SMALL)
+ {
+ leVectors[edgei] = vector::uniform(SMALL);
+ }
}
// Copy internal field
@@ -672,6 +684,15 @@ void Foam::faMesh::calcLe() const
{
const faPatch& fap = boundary()[patchi];
bfld[patchi] = fap.patchRawSlice(leVectors);
+
+ for (auto& patchEdge : bfld[patchi])
+ {
+ // Do not allow any mag(val) < SMALL
+ if (mag(patchEdge) < SMALL)
+ {
+ patchEdge = vector::uniform(SMALL);
+ }
+ }
}
}
else
@@ -692,6 +713,12 @@ void Foam::faMesh::calcLe() const
edges_[edgei].line(localPoints),
edgeNormals[edgei]
);
+
+ // Do not allow any mag(val) < SMALL
+ if (mag(fld[edgei]) < SMALL)
+ {
+ fld[edgei] = vector::uniform(SMALL);
+ }
}
}
@@ -715,6 +742,12 @@ void Foam::faMesh::calcLe() const
bndEdgeNormals[patchEdgei]
);
+ // Do not allow any mag(val) < SMALL
+ if (mag(pfld[patchEdgei]) < SMALL)
+ {
+ pfld[patchEdgei] = vector::uniform(SMALL);
+ }
+
++edgei;
}
}
@@ -759,6 +792,13 @@ void Foam::faMesh::calcMagLe() const
for (const edge& e : internalEdges())
{
*iter = e.mag(localPoints);
+
+ // Do not allow any mag(val) < SMALL
+ if (mag(*iter) < SMALL)
+ {
+ *iter = SMALL;
+ }
+
++iter;
}
}
@@ -774,6 +814,13 @@ void Foam::faMesh::calcMagLe() const
for (const edge& e : boundary()[patchi].patchSlice(edges_))
{
*iter = e.mag(localPoints);
+
+ // Do not allow any mag(val) < SMALL
+ if (mag(*iter) < SMALL)
+ {
+ *iter = SMALL;
+ }
+
++iter;
}
}
@@ -955,6 +1002,12 @@ void Foam::faMesh::calcS() const
forAll(fld, facei)
{
fld[facei] = Foam::mag(meshFaceAreas[facei]);
+
+ // Do not allow any mag(val) < SMALL
+ if (mag(fld[facei]) < SMALL)
+ {
+ fld[facei] = SMALL;
+ }
}
}
else
@@ -968,6 +1021,13 @@ void Foam::faMesh::calcS() const
for (const face& f : faces())
{
*iter = f.mag(localPoints);
+
+ // Do not allow any mag(val) < SMALL
+ if (mag(*iter) < SMALL)
+ {
+ *iter = SMALL;
+ }
+
++iter;
}
}
@@ -1028,6 +1088,15 @@ void Foam::faMesh::calcFaceAreaNormals() const
// Make unit normals
fld.normalise();
+
+ for (auto& f : fld)
+ {
+ // Do not allow any mag(val) < SMALL
+ if (mag(f) < SMALL)
+ {
+ f = vector::uniform(SMALL);
+ }
+ }
}
@@ -1122,6 +1191,12 @@ void Foam::faMesh::calcEdgeAreaNormals() const
fld[edgei].removeCollinear(edgeLine.unitVec());
fld[edgei].normalise();
+
+ // Do not allow any mag(val) < SMALL
+ if (mag(fld[edgei]) < SMALL)
+ {
+ fld[edgei] = vector::uniform(SMALL);
+ }
}
}
@@ -1148,6 +1223,12 @@ void Foam::faMesh::calcEdgeAreaNormals() const
pfld[patchEdgei].removeCollinear(edgeLine.unitVec());
pfld[patchEdgei].normalise();
+
+ // Do not allow any mag(val) < SMALL
+ if (mag(pfld[patchEdgei]) < SMALL)
+ {
+ pfld[patchEdgei] = vector::uniform(SMALL);
+ }
}
}
}
diff --git a/src/finiteArea/faMesh/faPatches/faPatch/faPatch.C b/src/finiteArea/faMesh/faPatches/faPatch/faPatch.C
index a6848fc4944..595ffb0702d 100644
--- a/src/finiteArea/faMesh/faPatches/faPatch/faPatch.C
+++ b/src/finiteArea/faMesh/faPatches/faPatch/faPatch.C
@@ -478,7 +478,19 @@ Foam::tmp<Foam::vectorField> Foam::faPatch::delta() const
{
// Use patch-normal delta for all non-coupled BCs
const vectorField nHat(edgeNormals());
- return nHat*(nHat & (edgeCentres() - edgeFaceCentres()));
+
+ vectorField edgePN(edgeCentres() - edgeFaceCentres());
+
+ // Do not allow any mag(val) < SMALL
+ for (vector& edgei : edgePN)
+ {
+ if (mag(edgei) < SMALL)
+ {
+ edgei = vector::uniform(SMALL);
+ }
+ }
+
+ return nHat*(nHat & edgePN);
}
diff --git a/src/finiteArea/finiteArea/gradSchemes/leastSquaresFaGrad/leastSquaresFaVectors.C b/src/finiteArea/finiteArea/gradSchemes/leastSquaresFaGrad/leastSquaresFaVectors.C
index 153d22fbedc..88e35b3d511 100644
--- a/src/finiteArea/finiteArea/gradSchemes/leastSquaresFaGrad/leastSquaresFaVectors.C
+++ b/src/finiteArea/finiteArea/gradSchemes/leastSquaresFaGrad/leastSquaresFaVectors.C
@@ -115,6 +115,13 @@ void Foam::leastSquaresFaVectors::makeLeastSquaresVectors() const
label nei = neighbour[facei];
vector d = C[nei] - C[own];
+
+ // Do not allow any mag(val) < SMALL
+ if (mag(d) < SMALL)
+ {
+ d = vector::uniform(SMALL);
+ }
+
symmTensor wdd = (1.0/magSqr(d))*sqr(d);
dd[own] += wdd;
@@ -169,6 +176,13 @@ void Foam::leastSquaresFaVectors::makeLeastSquaresVectors() const
label nei = neighbour[facei];
vector d = C[nei] - C[own];
+
+ // Do not allow any mag(val) < SMALL
+ if (mag(d) < SMALL)
+ {
+ d = vector::uniform(SMALL);
+ }
+
scalar magSfByMagSqrd = 1.0/magSqr(d);
lsP[facei] = magSfByMagSqrd*(invDd[own] & d);
diff --git a/src/finiteArea/interpolation/edgeInterpolation/edgeInterpolation.C b/src/finiteArea/interpolation/edgeInterpolation/edgeInterpolation.C
index ad65518493a..b57f0e75f23 100644
--- a/src/finiteArea/interpolation/edgeInterpolation/edgeInterpolation.C
+++ b/src/finiteArea/interpolation/edgeInterpolation/edgeInterpolation.C
@@ -248,6 +248,12 @@ void Foam::edgeInterpolation::makeLPN() const
);
lPNIn[edgeI] = (lPE + lEN);
+
+ // Do not allow any mag(val) < SMALL
+ if (mag(lPNIn[edgeI]) < SMALL)
+ {
+ lPNIn[edgeI] = SMALL;
+ }
}
@@ -287,7 +293,7 @@ void Foam::edgeInterpolation::makeWeights() const
false
),
mesh(),
- dimless
+ dimensionedScalar(dimless, 1)
);
edgeScalarField& weightingFactors = *weightingFactors_;
@@ -323,15 +329,12 @@ void Foam::edgeInterpolation::makeWeights() const
+ skewCorrEdge
);
- weightingFactorsIn[edgeI] =
- lEN
- /(
- lPE
-# ifdef BAD_MESH_STABILISATION
- + VSMALL
-# endif
- + lEN
- );
+ // weight = (0,1]
+ const scalar lPN = lPE + lEN;
+ if (mag(lPN) > SMALL)
+ {
+ weightingFactorsIn[edgeI] = lEN/lPN;
+ }
}
forAll(mesh().boundary(), patchI)
@@ -370,7 +373,7 @@ void Foam::edgeInterpolation::makeDeltaCoeffs() const
false
),
mesh(),
- dimless/dimLength
+ dimensionedScalar(dimless/dimLength, SMALL)
);
edgeScalarField& DeltaCoeffs = *differenceFactors_;
scalarField& dc = DeltaCoeffs.primitiveFieldRef();
@@ -427,11 +430,12 @@ void Foam::edgeInterpolation::makeDeltaCoeffs() const
// Edge normal - area tangent
edgeNormal = normalised(lengths[edgeI]);
- // Delta coefficient as per definition
-// dc[edgeI] = 1.0/(lPN*(unitDelta & edgeNormal));
-
- // Stabilised form for bad meshes. HJ, 23/Jul/2009
- dc[edgeI] = 1.0/max((lPN*(unitDelta & edgeNormal)), 0.05*lPN);
+ // Do not allow any mag(val) < SMALL
+ const scalar alpha = lPN*(unitDelta & edgeNormal);
+ if (mag(alpha) > SMALL)
+ {
+ dc[edgeI] = scalar(1)/max(alpha, 0.05*lPN);
+ }
}
@@ -478,7 +482,7 @@ void Foam::edgeInterpolation::makeCorrectionVectors() const
const edgeList& edges = mesh().edges();
const pointField& points = mesh().points();
- scalarField deltaCoeffs(owner.size());
+ scalarField deltaCoeffs(owner.size(), SMALL);
vectorField& CorrVecsIn = CorrVecs.primitiveFieldRef();
@@ -499,8 +503,12 @@ void Foam::edgeInterpolation::makeCorrectionVectors() const
// Edge normal - area tangent
edgeNormal = normalised(lengths[edgeI]);
- // Delta coeffs
- deltaCoeffs[edgeI] = 1.0/(unitDelta & edgeNormal);
+ // Do not allow any mag(val) < SMALL
+ const scalar alpha = unitDelta & edgeNormal;
+ if (mag(alpha) > SMALL)
+ {
+ deltaCoeffs[edgeI] = scalar(1)/alpha;
+ }
// Edge correction vector
CorrVecsIn[edgeI] =
@@ -570,7 +578,7 @@ void Foam::edgeInterpolation::makeSkewCorrectionVectors() const
false
),
mesh(),
- dimless
+ dimensionedVector(dimless, Zero)
);
edgeVectorField& SkewCorrVecs = *skewCorrectionVectors_;
@@ -590,19 +598,22 @@ void Foam::edgeInterpolation::makeSkewCorrectionVectors() const
const vector& P = C[owner[edgeI]];
const vector& N = C[neighbour[edgeI]];
const vector& S = points[edges[edgeI].start()];
- vector e = edges[edgeI].vec(points);
- const scalar beta = magSqr((N - P)^e);
+ // (T:Eq. 5.4)
+ const vector d(N - P);
+ const vector e(edges[edgeI].vec(points));
+ const vector de(d^e);
+ const scalar alpha = magSqr(de);
- if (beta < ROOTVSMALL)
+ if (alpha < SMALL)
{
// Too small - skew correction remains zero
continue;
}
+ const scalar beta = -((d^(S - P)) & de)/alpha;
- const scalar alpha = -(((N - P)^(S - P))&((N - P)^e))/beta;
-
- vector E(S + alpha*e);
+ // (T:Eq. 5.3)
+ const vector E(S + beta*e);
SkewCorrVecs[edgeI] = Ce[edgeI] - E;
}
@@ -630,28 +641,26 @@ void Foam::edgeInterpolation::makeSkewCorrectionVectors() const
const vector& P = C[edgeFaces[edgeI]];
const vector& N = ngbC[edgeI];
const vector& S = points[patchEdges[edgeI].start()];
- vector e = patchEdges[edgeI].vec(points);
- const scalar beta = magSqr((N - P)^e);
+ // (T:Eq. 5.4)
+ const vector d(N - P);
+ const vector e(patchEdges[edgeI].vec(points));
+ const vector de(d^e);
+ const scalar alpha = magSqr(de);
- if (beta < ROOTVSMALL)
+ if (alpha < SMALL)
{
// Too small - skew correction remains zero
continue;
}
+ const scalar beta = -((d^(S - P)) & de)/alpha;
- const scalar alpha = -(((N - P)^(S - P))&((N - P)^e))/beta;
-
- vector E(S + alpha*e);
+ const vector E(S + beta*e);
patchSkewCorrVecs[edgeI] =
Ce.boundaryField()[patchI][edgeI] - E;
}
}
- else
- {
- patchSkewCorrVecs = vector::zero;
- }
}
#ifdef FA_SKEW_CORRECTION
diff --git a/src/finiteArea/interpolation/edgeInterpolation/schemes/NVDscheme/faNVDscheme.C b/src/finiteArea/interpolation/edgeInterpolation/schemes/NVDscheme/faNVDscheme.C
index c659ad3f4f0..8b7c6850109 100644
--- a/src/finiteArea/interpolation/edgeInterpolation/schemes/NVDscheme/faNVDscheme.C
+++ b/src/finiteArea/interpolation/edgeInterpolation/schemes/NVDscheme/faNVDscheme.C
@@ -108,6 +108,12 @@ Foam::tmp<Foam::edgeScalarField> Foam::faNVDscheme<Type,NVDweight>::weights
d.normalise();
d *= mesh.edgeInterpolation::lPN().internalField()[edge];
+ // Do not allow any mag(val) < SMALL
+ if (mag(d) < SMALL)
+ {
+ d = vector::uniform(SMALL);
+ }
+
weights[edge] =
this->weight
(
@@ -189,6 +195,12 @@ Foam::tmp<Foam::edgeScalarField> Foam::faNVDscheme<Type,NVDweight>::weights
d.normalise();
d *= pLPN[edgeI];
+ // Do not allow any mag(val) < SMALL
+ if (mag(d) < SMALL)
+ {
+ d = vector::uniform(SMALL);
+ }
+
pWeights[edgeI] =
this->weight
(
--
GitLab