songtianlun/openfoam
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

ENH: finiteArea: reduce various code footprint

Browse Source
This commit is contained in:
Kutalmis Bercin 2024-01-23 15:49:41 +00:00 committed by Andrew Heather
parent abfe30454a
commit d03a225061
6 changed files with 181 additions and 248 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

27
src/finiteArea/finiteArea/d2dt2Schemes/EulerFaD2dt2Scheme/EulerFaD2dt2Scheme.C
View File

@ -82,10 +82,10 @@ EulerFaD2dt2Scheme<Type>::facD2dt2
if (mesh().moving())
{
scalar halfRdeltaT2 = rDeltaT2.value()/2.0;
scalar halfRdeltaT2 = 0.5*rDeltaT2.value();
scalarField SS0 = mesh().S() + mesh().S0();
scalarField S0S00 = mesh().S0() + mesh().S00();
scalarField SS0(mesh().S() + mesh().S0());
scalarField S0S00(mesh().S0() + mesh().S00());
tmp<GeometricField<Type, faPatchField, areaMesh>> tdt2dt2
(
@ -148,10 +148,10 @@ EulerFaD2dt2Scheme<Type>::facD2dt2
if (mesh().moving())
{
scalar halfRdeltaT2 = rDeltaT2.value()/2.0;
scalar halfRdeltaT2 = 0.5*rDeltaT2.value();
scalarField SS0 = mesh().S() + mesh().S0();
scalarField S0S00 = mesh().S0() + mesh().S00();
scalarField SS0(mesh().S() + mesh().S0());
scalarField S0S00(mesh().S0() + mesh().S00());
return tmp<GeometricField<Type, faPatchField, areaMesh>>
(
@ -228,13 +228,15 @@ EulerFaD2dt2Scheme<Type>::facD2dt2
{
scalar halfRdeltaT2 = 0.5*rDeltaT2.value();
scalarField SS0rhoRho0 =
(mesh().S() + mesh().S0())
*rho.value();
scalarField SS0rhoRho0
(
(mesh().S() + mesh().S0())*rho.value()
);
scalarField S0S00rho0Rho00 =
(mesh().S0() + mesh().S00())
*rho.value();
scalarField S0S00rho0Rho00
(
(mesh().S0() + mesh().S00())*rho.value()
);
return tmp<GeometricField<Type, faPatchField, areaMesh>>
(
@ -489,7 +491,6 @@ EulerFaD2dt2Scheme<Type>::famD2dt2
scalar halfRdeltaT2 = 0.5*rDeltaT2;
scalarField SS0(mesh().S() + mesh().S0());
scalarField S0S00(mesh().S0() + mesh().S00());
fam.diag() = rho.value()*(coefft*halfRdeltaT2)*SS0;

6
src/finiteArea/finiteArea/ddtSchemes/boundedBackwardFaDdtScheme/boundedBackwardFaDdtScheme.H
View File

@ -76,10 +76,8 @@ class boundedBackwardFaDdtScheme
{
return GREAT;
}
else
{
return deltaT0_();
}
return deltaT0_();
}

42
src/finiteArea/finiteArea/gradSchemes/leastSquaresFaGrad/leastSquaresFaGrad.C
View File

@ -91,10 +91,10 @@ leastSquaresFaGrad<Type>::calcGrad
forAll(own, edgei)
{
label ownEdgeI = own[edgei];
label neiEdgeI = nei[edgei];
const label ownEdgeI = own[edgei];
const label neiEdgeI = nei[edgei];
Type deltaVsf = vsf[neiEdgeI] - vsf[ownEdgeI];
const Type deltaVsf(vsf[neiEdgeI] - vsf[ownEdgeI]);
lsGrad[ownEdgeI] += ownLs[edgei]*deltaVsf;
lsGrad[neiEdgeI] -= neiLs[edgei]*deltaVsf;
@ -103,35 +103,23 @@ leastSquaresFaGrad<Type>::calcGrad
// Boundary edges
forAll(vsf.boundaryField(), patchi)
{
const faePatchVectorField& patchOwnLs = ownLs.boundaryField()[patchi];
const faPatchField<Type>& bf = vsf.boundaryField()[patchi];
const Field<Type>& vsfp =
(
bf.coupled()
? bf.patchNeighbourField().cref()
: const_cast<faPatchField<Type>&>(bf)
);
const faePatchVectorField& ownLsp = ownLs.boundaryField()[patchi];
const labelUList& edgeFaces =
lsGrad.boundaryField()[patchi].patch().edgeFaces();
if (vsf.boundaryField()[patchi].coupled())
forAll(vsfp, pEdgei)
{
Field<Type> neiVsf
(
vsf.boundaryField()[patchi].patchNeighbourField()
);
forAll(neiVsf, patchEdgeI)
{
lsGrad[edgeFaces[patchEdgeI]] +=
patchOwnLs[patchEdgeI]
*(neiVsf[patchEdgeI] - vsf[edgeFaces[patchEdgeI]]);
}
}
else
{
const faPatchField<Type>& patchVsf = vsf.boundaryField()[patchi];
forAll(patchVsf, patchEdgeI)
{
lsGrad[edgeFaces[patchEdgeI]] +=
patchOwnLs[patchEdgeI]
*(patchVsf[patchEdgeI] - vsf[edgeFaces[patchEdgeI]]);
}
lsGrad[edgeFaces[pEdgei]] +=
ownLsp[pEdgei]*(vsfp[pEdgei] - vsf[edgeFaces[pEdgei]]);
}
}

216
src/finiteArea/finiteArea/gradSchemes/limitedGradSchemes/edgeLimitedFaGrad/edgeLimitedFaGrads.C
View File

@ -91,22 +91,22 @@ tmp<areaVectorField> edgeLimitedGrad<scalar>::calcGrad
const areaVectorField& C = mesh.areaCentres();
const edgeVectorField& Cf = mesh.edgeCentres();
// create limiter
// Create limiter field
scalarField limiter(vsf.internalField().size(), 1.0);
scalar rk = (1.0/k_ - 1.0);
const scalar rk = (1.0/k_ - 1.0);
forAll(owner, edgei)
{
label own = owner[edgei];
label nei = neighbour[edgei];
const label own = owner[edgei];
const label nei = neighbour[edgei];
scalar vsfOwn = vsf[own];
scalar vsfNei = vsf[nei];
const scalar vsfOwn = vsf[own];
const scalar vsfNei = vsf[nei];
scalar maxEdge = max(vsfOwn, vsfNei);
scalar minEdge = min(vsfOwn, vsfNei);
scalar maxMinEdge = rk*(maxEdge - minEdge);
const scalar maxMinEdge = rk*(maxEdge - minEdge);
maxEdge += maxMinEdge;
minEdge -= maxMinEdge;
@ -114,7 +114,8 @@ tmp<areaVectorField> edgeLimitedGrad<scalar>::calcGrad
limitEdge
(
limiter[own],
maxEdge - vsfOwn, minEdge - vsfOwn,
maxEdge - vsfOwn,
minEdge - vsfOwn,
(Cf[edgei] - C[own]) & g[own]
);
@ -122,67 +123,53 @@ tmp<areaVectorField> edgeLimitedGrad<scalar>::calcGrad
limitEdge
(
limiter[nei],
maxEdge - vsfNei, minEdge - vsfNei,
maxEdge - vsfNei,
minEdge - vsfNei,
(Cf[edgei] - C[nei]) & g[nei]
);
}
const areaScalarField::Boundary& bsf = vsf.boundaryField();
// Lambda expression to update limiter for boundary edges
auto updateLimiter = [&](const label patchi, const scalarField& fld) -> void
{
const labelUList& pOwner = mesh.boundary()[patchi].edgeFaces();
const vectorField& pCf = Cf.boundaryField()[patchi];
forAll(pOwner, edgei)
{
const label own = pOwner[edgei];
const scalar vsfOwn = vsf[own];
const scalar vsfNei = fld[edgei];
scalar maxEdge = max(vsfOwn, vsfNei);
scalar minEdge = min(vsfOwn, vsfNei);
const scalar maxMinEdge = rk*(maxEdge - minEdge);
maxEdge += maxMinEdge;
minEdge -= maxMinEdge;
limitEdge
(
limiter[own],
maxEdge - vsfOwn,
minEdge - vsfOwn,
(pCf[edgei] - C[own]) & g[own]
);
}
};
const areaScalarField::Boundary& bsf = vsf.boundaryField();
forAll(bsf, patchi)
{
const faPatchScalarField& psf = bsf[patchi];
const labelUList& pOwner = mesh.boundary()[patchi].edgeFaces();
const vectorField& pCf = Cf.boundaryField()[patchi];
if (psf.coupled())
{
const scalarField psfNei(psf.patchNeighbourField());
forAll(pOwner, pEdgei)
{
label own = pOwner[pEdgei];
scalar vsfOwn = vsf[own];
scalar vsfNei = psfNei[pEdgei];
scalar maxEdge = max(vsfOwn, vsfNei);
scalar minEdge = min(vsfOwn, vsfNei);
scalar maxMinEdge = rk*(maxEdge - minEdge);
maxEdge += maxMinEdge;
minEdge -= maxMinEdge;
limitEdge
(
limiter[own],
maxEdge - vsfOwn, minEdge - vsfOwn,
(pCf[pEdgei] - C[own]) & g[own]
);
}
updateLimiter(patchi, psf.patchNeighbourField());
}
else if (psf.fixesValue())
{
forAll(pOwner, pEdgei)
{
label own = pOwner[pEdgei];
scalar vsfOwn = vsf[own];
scalar vsfNei = psf[pEdgei];
scalar maxEdge = max(vsfOwn, vsfNei);
scalar minEdge = min(vsfOwn, vsfNei);
scalar maxMinEdge = rk*(maxEdge - minEdge);
maxEdge += maxMinEdge;
minEdge -= maxMinEdge;
limitEdge
(
limiter[own],
maxEdge - vsfOwn, minEdge - vsfOwn,
(pCf[pEdgei] - C[own]) & g[own]
);
}
updateLimiter(patchi, psf);
}
}
@ -226,35 +213,33 @@ tmp<areaTensorField> edgeLimitedGrad<vector>::calcGrad
const areaVectorField& C = mesh.areaCentres();
const edgeVectorField& Cf = mesh.edgeCentres();
// create limiter
// Create limiter
scalarField limiter(vvf.internalField().size(), 1.0);
scalar rk = (1.0/k_ - 1.0);
const scalar rk = (1.0/k_ - 1.0);
forAll(owner, edgei)
{
label own = owner[edgei];
label nei = neighbour[edgei];
vector vvfOwn = vvf[own];
vector vvfNei = vvf[nei];
const label own = owner[edgei];
const label nei = neighbour[edgei];
// owner side
vector gradf = (Cf[edgei] - C[own]) & g[own];
vector gradf((Cf[edgei] - C[own]) & g[own]);
scalar vsfOwn = gradf & vvfOwn;
scalar vsfNei = gradf & vvfNei;
scalar vsfOwn = gradf & vvf[own];
scalar vsfNei = gradf & vvf[nei];
scalar maxEdge = max(vsfOwn, vsfNei);
scalar minEdge = min(vsfOwn, vsfNei);
scalar maxMinEdge = rk*(maxEdge - minEdge);
const scalar maxMinEdge = rk*(maxEdge - minEdge);
maxEdge += maxMinEdge;
minEdge -= maxMinEdge;
limitEdge
(
limiter[own],
maxEdge - vsfOwn, minEdge - vsfOwn,
maxEdge - vsfOwn,
minEdge - vsfOwn,
magSqr(gradf)
);
@ -262,8 +247,8 @@ tmp<areaTensorField> edgeLimitedGrad<vector>::calcGrad
// neighbour side
gradf = (Cf[edgei] - C[nei]) & g[nei];
vsfOwn = gradf & vvfOwn;
vsfNei = gradf & vvfNei;
vsfOwn = gradf & vvf[own];
vsfNei = gradf & vvf[nei];
maxEdge = max(vsfOwn, vsfNei);
minEdge = min(vsfOwn, vsfNei);
@ -271,78 +256,57 @@ tmp<areaTensorField> edgeLimitedGrad<vector>::calcGrad
limitEdge
(
limiter[nei],
maxEdge - vsfNei, minEdge - vsfNei,
maxEdge - vsfNei,
minEdge - vsfNei,
magSqr(gradf)
);
}
const areaVectorField::Boundary& bvf = vvf.boundaryField();
// Lambda expression to update limiter for boundary edges
auto updateLimiter = [&](const label patchi, const vectorField& fld) -> void
{
const labelUList& pOwner = mesh.boundary()[patchi].edgeFaces();
const vectorField& pCf = Cf.boundaryField()[patchi];
forAll(pOwner, edgei)
{
const label own = pOwner[edgei];
const vector gradf((pCf[edgei] - C[own]) & g[own]);
const scalar vsfOwn = gradf & vvf[own];
const scalar vsfNei = gradf & fld[edgei];
scalar maxEdge = max(vsfOwn, vsfNei);
scalar minEdge = min(vsfOwn, vsfNei);
const scalar maxMinEdge = rk*(maxEdge - minEdge);
maxEdge += maxMinEdge;
minEdge -= maxMinEdge;
limitEdge
(
limiter[own],
maxEdge - vsfOwn,
minEdge - vsfOwn,
magSqr(gradf)
);
}
};
const areaVectorField::Boundary& bvf = vvf.boundaryField();
forAll(bvf, patchi)
{
const faPatchVectorField& psf = bvf[patchi];
const labelUList& pOwner = mesh.boundary()[patchi].edgeFaces();
const vectorField& pCf = Cf.boundaryField()[patchi];
if (psf.coupled())
{
const vectorField psfNei(psf.patchNeighbourField());
forAll(pOwner, pEdgei)
{
label own = pOwner[pEdgei];
vector vvfOwn = vvf[own];
vector vvfNei = psfNei[pEdgei];
vector gradf = (pCf[pEdgei] - C[own]) & g[own];
scalar vsfOwn = gradf & vvfOwn;
scalar vsfNei = gradf & vvfNei;
scalar maxEdge = max(vsfOwn, vsfNei);
scalar minEdge = min(vsfOwn, vsfNei);
scalar maxMinEdge = rk*(maxEdge - minEdge);
maxEdge += maxMinEdge;
minEdge -= maxMinEdge;
limitEdge
(
limiter[own],
maxEdge - vsfOwn, minEdge - vsfOwn,
magSqr(gradf)
);
}
updateLimiter(patchi, psf.patchNeighbourField());
}
else if (psf.fixesValue())
{
forAll(pOwner, pEdgei)
{
label own = pOwner[pEdgei];
vector vvfOwn = vvf[own];
vector vvfNei = psf[pEdgei];
vector gradf = (pCf[pEdgei] - C[own]) & g[own];
scalar vsfOwn = gradf & vvfOwn;
scalar vsfNei = gradf & vvfNei;
scalar maxEdge = max(vsfOwn, vsfNei);
scalar minEdge = min(vsfOwn, vsfNei);
scalar maxMinEdge = rk*(maxEdge - minEdge);
maxEdge += maxMinEdge;
minEdge -= maxMinEdge;
limitEdge
(
limiter[own],
maxEdge - vsfOwn, minEdge - vsfOwn,
magSqr(gradf)
);
}
updateLimiter(patchi, psf);
}
}

107
src/finiteArea/finiteArea/gradSchemes/limitedGradSchemes/faceLimitedFaGrad/faceLimitedFaGrads.C
View File

@ -118,11 +118,11 @@ tmp<areaVectorField> faceLimitedGrad<scalar>::calcGrad
forAll(owner, facei)
{
label own = owner[facei];
label nei = neighbour[facei];
const label own = owner[facei];
const label nei = neighbour[facei];
scalar vsfOwn = vsf[own];
scalar vsfNei = vsf[nei];
const scalar vsfOwn = vsf[own];
const scalar vsfNei = vsf[nei];
maxVsf[own] = max(maxVsf[own], vsfNei);
minVsf[own] = min(minVsf[own], vsfNei);
@ -132,37 +132,33 @@ tmp<areaVectorField> faceLimitedGrad<scalar>::calcGrad
}
const areaScalarField::Boundary& bsf = vsf.boundaryField();
// Lambda expression to update limiter for boundary edges
auto updateLimiter = [&](const label patchi, const scalarField& fld) -> void
{
const labelUList& pOwner = mesh.boundary()[patchi].edgeFaces();
forAll(pOwner, facei)
{
const label own = pOwner[facei];
const scalar vsf = fld[facei];
maxVsf[own] = max(maxVsf[own], vsf);
minVsf[own] = min(minVsf[own], vsf);
}
};
const areaScalarField::Boundary& bsf = vsf.boundaryField();
forAll(bsf, patchi)
{
const faPatchScalarField& psf = bsf[patchi];
const labelUList& pOwner = mesh.boundary()[patchi].edgeFaces();
if (psf.coupled())
{
scalarField psfNei(psf.patchNeighbourField());
forAll(pOwner, pFacei)
{
label own = pOwner[pFacei];
scalar vsfNei = psfNei[pFacei];
maxVsf[own] = max(maxVsf[own], vsfNei);
minVsf[own] = min(minVsf[own], vsfNei);
}
updateLimiter(patchi, psf.patchNeighbourField());
}
else
{
forAll(pOwner, pFacei)
{
label own = pOwner[pFacei];
scalar vsfNei = psf[pFacei];
maxVsf[own] = max(maxVsf[own], vsfNei);
minVsf[own] = min(minVsf[own], vsfNei);
}
updateLimiter(patchi, psf);
}
}
@ -171,19 +167,19 @@ tmp<areaVectorField> faceLimitedGrad<scalar>::calcGrad
if (k_ < 1.0)
{
scalarField maxMinVsf((1.0/k_ - 1.0)*(maxVsf - minVsf));
const scalarField maxMinVsf((1.0/k_ - 1.0)*(maxVsf - minVsf));
maxVsf += maxMinVsf;
minVsf -= maxMinVsf;
}
// create limiter
// Create limiter
scalarField limiter(vsf.internalField().size(), 1.0);
forAll(owner, facei)
{
label own = owner[facei];
label nei = neighbour[facei];
const label own = owner[facei];
const label nei = neighbour[facei];
// owner side
limitEdge
@ -211,7 +207,7 @@ tmp<areaVectorField> faceLimitedGrad<scalar>::calcGrad
forAll(pOwner, pFacei)
{
label own = pOwner[pFacei];
const label own = pOwner[pFacei];
limitEdge
(
@ -268,8 +264,8 @@ tmp<areaTensorField> faceLimitedGrad<vector>::calcGrad
forAll(owner, facei)
{
label own = owner[facei];
label nei = neighbour[facei];
const label own = owner[facei];
const label nei = neighbour[facei];
const vector& vsfOwn = vsf[own];
const vector& vsfNei = vsf[nei];
@ -282,36 +278,33 @@ tmp<areaTensorField> faceLimitedGrad<vector>::calcGrad
}
const areaVectorField::Boundary& bsf = vsf.boundaryField();
// Lambda expression to update limiter for boundary edges
auto updateLimiter = [&](const label patchi, const vectorField& fld) -> void
{
const labelUList& pOwner = mesh.boundary()[patchi].edgeFaces();
forAll(pOwner, facei)
{
const label own = pOwner[facei];
const vector& vsf = fld[facei];
maxVsf[own] = max(maxVsf[own], vsf);
minVsf[own] = min(minVsf[own], vsf);
}
};
const areaVectorField::Boundary& bsf = vsf.boundaryField();
forAll(bsf, patchi)
{
const faPatchVectorField& psf = bsf[patchi];
const labelUList& pOwner = mesh.boundary()[patchi].edgeFaces();
if (psf.coupled())
{
vectorField psfNei(psf.patchNeighbourField());
forAll(pOwner, pFacei)
{
label own = pOwner[pFacei];
const vector& vsfNei = psfNei[pFacei];
maxVsf[own] = max(maxVsf[own], vsfNei);
minVsf[own] = min(minVsf[own], vsfNei);
}
updateLimiter(patchi, psf.patchNeighbourField());
}
else
{
forAll(pOwner, pFacei)
{
label own = pOwner[pFacei];
const vector& vsfNei = psf[pFacei];
maxVsf[own] = max(maxVsf[own], vsfNei);
minVsf[own] = min(minVsf[own], vsfNei);
}
updateLimiter(patchi, psf);
}
}
@ -320,7 +313,7 @@ tmp<areaTensorField> faceLimitedGrad<vector>::calcGrad
if (k_ < 1.0)
{
vectorField maxMinVsf((1.0/k_ - 1.0)*(maxVsf - minVsf));
const vectorField maxMinVsf((1.0/k_ - 1.0)*(maxVsf - minVsf));
maxVsf += maxMinVsf;
minVsf -= maxMinVsf;
@ -329,13 +322,13 @@ tmp<areaTensorField> faceLimitedGrad<vector>::calcGrad
}
// create limiter
// Create limiter
vectorField limiter(vsf.internalField().size(), vector::one);
forAll(owner, facei)
{
label own = owner[facei];
label nei = neighbour[facei];
const label own = owner[facei];
const label nei = neighbour[facei];
// owner side
limitEdge
@ -363,7 +356,7 @@ tmp<areaTensorField> faceLimitedGrad<vector>::calcGrad
forAll(pOwner, pFacei)
{
label own = pOwner[pFacei];
const label own = pOwner[pFacei];
limitEdge
(

31
src/finiteArea/interpolation/edgeInterpolation/schemes/Gamma/Gamma.H
View File

@ -76,7 +76,7 @@ public:
// Rescale k_ to be >= 0 and <= 0.5 (TVD conformant)
// and avoid the /0 when k_ = 0
k_ = max(k_/2.0, SMALL);
k_ = max(0.5*k_, SMALL);
}
@ -91,30 +91,19 @@ public:
const vector& d
) const
{
scalar magd = mag(d);
vector dHat = d/mag(d);
const vector dHat(normalised(d));
scalar gradf = (phiN - phiP)/magd;
scalar gradcf;
scalar udWeight;
if (faceFlux > 0)
{
gradcf = dHat & gradcP;
udWeight = 1;
}
else
{
gradcf = dHat & gradcN;
udWeight = 0;
}
// Choose gradc based on faceFlux
const vector& gradcPN = (faceFlux > 0) ? gradcP : gradcN;
const scalar udWeight = (faceFlux > 0) ? 1 : 0;
// Stabilise for division
gradcf = stabilise(gradcf, SMALL);
const scalar gradcf = stabilise(dHat & gradcPN, SMALL);
scalar phict = 1 - 0.5*gradf/gradcf;
scalar limiter = clamp(phict/k_, zero_one{});
const scalar gradf = (phiN - phiP)/mag(d);
const scalar phict = 1 - 0.5*gradf/gradcf;
const scalar limiter = clamp(phict/k_, zero_one{});
return lerp(udWeight, cdWeight, limiter);
}