zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-11-28 03:28:01 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

ENH: robuster finite-area least squares gradient

Browse Source 
- avoid division by zero for small edges
This commit is contained in:
Mark Olesen
2023-02-16 14:04:41 +01:00
parent 58e4cfbc8a
commit cbb153372e
1 changed files with 49 additions and 72 deletions
Download Patch File Download Diff File Expand all files Collapse all files

113
src/finiteArea/finiteArea/gradSchemes/leastSquaresFaGrad/leastSquaresFaVectors.C
View File

@ -64,7 +64,8 @@ Foam::leastSquaresFaVectors::~leastSquaresFaVectors()
void Foam::leastSquaresFaVectors::makeLeastSquaresVectors() const void Foam::leastSquaresFaVectors::makeLeastSquaresVectors() const
{ {
DebugInFunction DebugInFunction
<< "Constructing finite area least square gradient vectors" << nl; << "Constructing finite area (invDist) least square gradient vectors"
<< nl;
pVectorsPtr_ = new edgeVectorField pVectorsPtr_ = new edgeVectorField
( (
@ -75,7 +76,7 @@ void Foam::leastSquaresFaVectors::makeLeastSquaresVectors() const
mesh().thisDb(), mesh().thisDb(),
IOobject::NO_READ, IOobject::NO_READ,
IOobject::NO_WRITE, IOobject::NO_WRITE,
false IOobject::NO_REGISTER
), ),
mesh(), mesh(),
dimensionedVector(dimless/dimLength, Zero) dimensionedVector(dimless/dimLength, Zero)
@ -91,7 +92,7 @@ void Foam::leastSquaresFaVectors::makeLeastSquaresVectors() const
mesh().thisDb(), mesh().thisDb(),
IOobject::NO_READ, IOobject::NO_READ,
IOobject::NO_WRITE, IOobject::NO_WRITE,
false IOobject::NO_REGISTER
), ),
mesh(), mesh(),
dimensionedVector(dimless/dimLength, Zero) dimensionedVector(dimless/dimLength, Zero)
@ -103,63 +104,45 @@ void Foam::leastSquaresFaVectors::makeLeastSquaresVectors() const
const labelUList& neighbour = mesh().neighbour(); const labelUList& neighbour = mesh().neighbour();
const areaVectorField& C = mesh().areaCentres(); const areaVectorField& C = mesh().areaCentres();
const edgeScalarField& w = mesh().weights();
// Set up temporary storage for the dd tensor (before inversion) // Set up temporary storage for the dd tensor (before inversion)
symmTensorField dd(mesh().nFaces(), Zero); symmTensorField dd(mesh().nFaces(), Zero);
forAll(owner, facei) forAll(owner, facei)
{ {
label own = owner[facei]; const label own = owner[facei];
label nei = neighbour[facei]; const label nei = neighbour[facei];
vector d = C[nei] - C[own]; const vector d(C[nei] - C[own]);
// Do not allow any mag(val) < SMALL // No contribution when mag(val) < SMALL
if (d.magSqr() < ROOTSMALL) const scalar magSqrDist = d.magSqr();
if (magSqrDist > ROOTSMALL)
{ {
d = vector::uniform(SMALL); const symmTensor wdd(sqr(d)/magSqrDist);
}
symmTensor wdd = (1.0/d.magSqr())*sqr(d);
dd[own] += wdd; dd[own] += wdd;
dd[nei] += wdd; dd[nei] += wdd;
} }
}
forAll(lsP.boundaryField(), patchi) for (const auto& patchLsP : lsP.boundaryField())
{ {
const faePatchScalarField& pw = w.boundaryField()[patchi]; const faPatch& p = patchLsP.patch();
const faPatch& p = pw.patch();
const labelUList& edgeFaces = p.edgeFaces(); const labelUList& edgeFaces = p.edgeFaces();
// Build the d-vectors // Build the d-vectors
// HJ, reconsider deltas at the boundary, consistent with FVM
// Current implementation is good for fixedValue boundaries, but may
// cause problems with fixedGradient. HJ, 4/Oct/2010
const vectorField pd(p.delta()); const vectorField pd(p.delta());
if (p.coupled())
{
forAll(pd, patchFacei) forAll(pd, patchFacei)
{ {
const vector& d = pd[patchFacei]; const vector& d = pd[patchFacei];
dd[edgeFaces[patchFacei]] += // No contribution when mag(val) < SMALL
(1.0/magSqr(d))*sqr(d); const scalar magSqrDist = d.magSqr();
} if (magSqrDist > ROOTSMALL)
}
else
{ {
forAll(pd, patchFacei) dd[edgeFaces[patchFacei]] += sqr(d)/magSqrDist;
{
const vector& d = pd[patchFacei];
dd[edgeFaces[patchFacei]] +=
(1.0/magSqr(d))*sqr(d);
} }
} }
} }
@ -172,60 +155,54 @@ void Foam::leastSquaresFaVectors::makeLeastSquaresVectors() const
// Revisit all faces and calculate the lsP and lsN vectors // Revisit all faces and calculate the lsP and lsN vectors
forAll(owner, facei) forAll(owner, facei)
{ {
label own = owner[facei]; const label own = owner[facei];
label nei = neighbour[facei]; const label nei = neighbour[facei];
vector d = C[nei] - C[own]; const vector d(C[nei] - C[own]);
// Do not allow any mag(val) < SMALL // No contribution when mag(val) < SMALL
if (d.magSqr() < ROOTSMALL) const scalar magSqrDist = d.magSqr();
if (magSqrDist > ROOTSMALL)
{ {
d = vector::uniform(SMALL); lsP[facei] = (invDd[own] & d)/magSqrDist;
lsN[facei] = -(invDd[nei] & d)/magSqrDist;
}
// ... already zero from initialisation
// else
// {
// lsP[facei] = Zero;
// lsN[facei] = Zero;
// }
} }
scalar magSfByMagSqrd = 1.0/d.magSqr(); for (auto& patchLsP : lsP.boundaryFieldRef())
lsP[facei] = magSfByMagSqrd*(invDd[own] & d);
lsN[facei] = -magSfByMagSqrd*(invDd[nei] & d);
}
forAll(lsP.boundaryField(), patchi)
{ {
faePatchVectorField& patchLsP = lsP.boundaryFieldRef()[patchi]; const faPatch& p = patchLsP.patch();
const faePatchScalarField& pw = w.boundaryField()[patchi];
const faPatch& p = pw.patch();
const labelUList& edgeFaces = p.edgeFaces(); const labelUList& edgeFaces = p.edgeFaces();
// Build the d-vectors // Build the d-vectors
const vectorField pd(p.delta()); const vectorField pd(p.delta());
if (p.coupled())
{
forAll(pd, patchFacei) forAll(pd, patchFacei)
{ {
const vector& d = pd[patchFacei]; const vector& d = pd[patchFacei];
patchLsP[patchFacei] = // No contribution when mag(val) < SMALL
(1.0/magSqr(d)) const scalar magSqrDist = d.magSqr();
*(invDd[edgeFaces[patchFacei]] & d); if (magSqrDist > ROOTSMALL)
}
}
else
{ {
forAll(pd, patchFacei)
{
const vector& d = pd[patchFacei];
patchLsP[patchFacei] = patchLsP[patchFacei] =
(1.0/magSqr(d)) (invDd[edgeFaces[patchFacei]] & d)/magSqrDist;
*(invDd[edgeFaces[patchFacei]] & d);
} }
// ... already zero from initialisation
// else
// {
// patchLsP[patchFacei] = Zero;
// }
} }
} }
DebugInFunction DebugInfo
<< "Done constructing finite area least square gradient vectors" << nl; << "Done constructing finite area least square gradient vectors" << nl;
} }