diff --git a/src/functionObjects/field/stabilityBlendingFactor/stabilityBlendingFactor.C b/src/functionObjects/field/stabilityBlendingFactor/stabilityBlendingFactor.C
index bda9620c394b5c78f24aaef911d0866e54a5cfa4..827d882b79723c894afb5726a29d4b8b828c9b87 100644
--- a/src/functionObjects/field/stabilityBlendingFactor/stabilityBlendingFactor.C
+++ b/src/functionObjects/field/stabilityBlendingFactor/stabilityBlendingFactor.C
@@ -50,6 +50,37 @@ namespace functionObjects
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
+Foam::volScalarField&
+Foam::functionObjects::stabilityBlendingFactor::indicator()
+{
+ const word fldName = "blendedIndicator" + fieldName_;
+
+ auto* fldPtr = mesh_.getObjectPtr<volScalarField>(fldName);
+
+ if (!fldPtr)
+ {
+ fldPtr = new volScalarField
+ (
+ IOobject
+ (
+ "blendedIndicator" + fieldName_,
+ time_.timeName(),
+ mesh_,
+ IOobject::NO_READ,
+ IOobject::NO_WRITE
+ ),
+ mesh_,
+ dimensionedScalar(dimless, Zero),
+ zeroGradientFvPatchScalarField::typeName
+ );
+
+ mesh_.objectRegistry::store(fldPtr);
+ }
+
+ return *fldPtr;
+}
+
+
void Foam::functionObjects::stabilityBlendingFactor::calcStats
(
label& nCellsScheme1,
@@ -57,21 +88,31 @@ void Foam::functionObjects::stabilityBlendingFactor::calcStats
label& nCellsBlended
) const
{
- forAll(indicator_, celli)
+ const auto* indicatorPtr =
+ mesh_.cfindObject<volScalarField>("blendedIndicator" + fieldName_);
+
+ if (!indicatorPtr)
{
- scalar i = indicator_[celli];
+ FatalErrorInFunction
+ << "Indicator field not set"
+ << abort(FatalError);
+ }
+ const auto& indicator = *indicatorPtr;
+
+ for (const scalar i : indicator)
+ {
if (i < tolerance_)
{
- nCellsScheme2++;
+ ++nCellsScheme2;
}
else if (i > (1 - tolerance_))
{
- nCellsScheme1++;
+ ++nCellsScheme1;
}
else
{
- nCellsBlended++;
+ ++nCellsBlended;
}
}
@@ -106,13 +147,15 @@ bool Foam::functionObjects::stabilityBlendingFactor::init(bool first)
{
const auto* residualPtr = mesh_.findObject<IOField<scalar>>(residualName_);
+ auto& indicator = this->indicator();
+
if (residuals_)
{
if (!residualPtr)
{
WarningInFunction
<< "Could not find residual field : " << residualName_
- << ". The residual mode won't be " << nl
+ << ". The residual mode will not be " << nl
<< " considered for the blended field in the stability "
<< "blending factor. " << nl
<< " Please add the corresponding 'solverInfo'"
@@ -155,9 +198,9 @@ bool Foam::functionObjects::stabilityBlendingFactor::init(bool first)
)
);
- forAll(indicator_, i)
+ forAll(indicator, i)
{
- indicator_[i] = indicatorResidual[i];
+ indicator[i] = indicatorResidual[i];
}
}
}
@@ -170,26 +213,26 @@ bool Foam::functionObjects::stabilityBlendingFactor::init(bool first)
if (maxNonOrthogonality_ >= minNonOrthogonality_)
{
FatalErrorInFunction
- << " minNonOrthogonality should be larger than "
+ << "minNonOrthogonality should be larger than "
<< "maxNonOrthogonality."
<< exit(FatalError);
}
- indicator_ =
- max
- (
- indicator_,
- min
+ indicator =
+ max
(
- max
+ indicator,
+ min
(
- scalar(0),
- (*nonOrthPtr - maxNonOrthogonality_)
- / (minNonOrthogonality_ - maxNonOrthogonality_)
- ),
- scalar(1)
- )
- );
+ max
+ (
+ scalar(0),
+ (*nonOrthPtr - maxNonOrthogonality_)
+ / (minNonOrthogonality_ - maxNonOrthogonality_)
+ ),
+ scalar(1)
+ )
+ );
if (first)
{
@@ -198,33 +241,32 @@ bool Foam::functionObjects::stabilityBlendingFactor::init(bool first)
}
}
- const volScalarField* skewnessPtr =
- mesh_.findObject<volScalarField>(skewnessName_);
+ const auto* skewnessPtr = mesh_.cfindObject<volScalarField>(skewnessName_);
if (skewness_)
{
if (maxSkewness_ >= minSkewness_)
{
FatalErrorInFunction
- << " minSkewness should be larger than maxSkewness."
+ << "minSkewness should be larger than maxSkewness."
<< exit(FatalError);
}
- indicator_ =
- max
- (
- indicator_,
- min
+ indicator =
+ max
(
- max
+ indicator,
+ min
(
- scalar(0),
- (*skewnessPtr - maxSkewness_)
- / (minSkewness_ - maxSkewness_)
- ),
- scalar(1)
- )
- );
+ max
+ (
+ scalar(0),
+ (*skewnessPtr - maxSkewness_)
+ / (minSkewness_ - maxSkewness_)
+ ),
+ scalar(1)
+ )
+ );
if (first)
{
@@ -233,22 +275,22 @@ bool Foam::functionObjects::stabilityBlendingFactor::init(bool first)
}
}
- const volScalarField* faceWeightsPtr =
- mesh_.findObject<volScalarField>(faceWeightName_);
+ const auto* faceWeightsPtr =
+ mesh_.cfindObject<volScalarField>(faceWeightName_);
if (faceWeight_)
{
if (maxFaceWeight_ >= minFaceWeight_)
{
FatalErrorInFunction
- << " minFaceWeight should be larger than maxFaceWeight."
+ << "minFaceWeight should be larger than maxFaceWeight."
<< exit(FatalError);
}
- indicator_ =
+ indicator =
max
(
- indicator_,
+ indicator,
min
(
max
@@ -274,7 +316,7 @@ bool Foam::functionObjects::stabilityBlendingFactor::init(bool first)
if (maxGradCc_ >= minGradCc_)
{
FatalErrorInFunction
- << " minGradCc should be larger than maxGradCc."
+ << "minGradCc should be larger than maxGradCc."
<< exit(FatalError);
}
@@ -322,10 +364,10 @@ bool Foam::functionObjects::stabilityBlendingFactor::init(bool first)
<< endl;
}
- indicator_ =
+ indicator =
max
(
- indicator_,
+ indicator,
min
(
max
@@ -340,15 +382,14 @@ bool Foam::functionObjects::stabilityBlendingFactor::init(bool first)
}
- const volVectorField* UNamePtr =
- mesh_.findObject<volVectorField>(UName_);
+ const auto* UNamePtr = mesh_.findObject<volVectorField>(UName_);
if (Co_)
{
if (Co1_ >= Co2_)
{
FatalErrorInFunction
- << " Co2 should be larger than Co1."
+ << "Co2 should be larger than Co1."
<< exit(FatalError);
}
@@ -372,10 +413,10 @@ bool Foam::functionObjects::stabilityBlendingFactor::init(bool first)
Co.primitiveFieldRef() =
mesh_.time().deltaT()*mag(*UNamePtr)/cbrt(mesh_.V());
- indicator_ =
+ indicator =
max
(
- indicator_,
+ indicator,
min(max(scalar(0), (Co - Co1_)/(Co2_ - Co1_)), scalar(1))
);
@@ -406,15 +447,14 @@ bool Foam::functionObjects::stabilityBlendingFactor::init(bool first)
<< endl;
}
- indicator_.correctBoundaryConditions();
- indicator_.min(1.0);
- indicator_.max(0.0);
+ indicator.correctBoundaryConditions();
+ indicator.min(1.0);
+ indicator.max(0.0);
// Update the blended surface field
- surfaceScalarField& surBlended =
- mesh_.lookupObjectRef<surfaceScalarField>(resultName_);
+ auto& surBlended = mesh_.lookupObjectRef<surfaceScalarField>(resultName_);
- surBlended = fvc::interpolate(indicator_);
+ surBlended = fvc::interpolate(indicator);
return true;
}
@@ -431,20 +471,6 @@ Foam::functionObjects::stabilityBlendingFactor::stabilityBlendingFactor
:
fieldExpression(name, runTime, dict),
writeFile(obr_, name, typeName, dict),
- indicator_
- (
- IOobject
- (
- "blendedIndicator" + fieldName_,
- time_.timeName(),
- mesh_,
- IOobject::NO_READ,
- IOobject::NO_WRITE
- ),
- mesh_,
- dimensionedScalar(dimless, Zero),
- zeroGradientFvPatchScalarField::typeName
- ),
nonOrthogonality_(dict.getOrDefault<Switch>("switchNonOrtho", false)),
gradCc_(dict.getOrDefault<Switch>("switchGradCc", false)),
residuals_(dict.getOrDefault<Switch>("switchResiduals", false)),
@@ -522,104 +548,93 @@ Foam::functionObjects::stabilityBlendingFactor::stabilityBlendingFactor
);
store(resultName_, faceBlendedPtr);
- const volScalarField* nonOrthPtr =
+ const auto* nonOrthPtr =
mesh_.findObject<volScalarField>(nonOrthogonalityName_);
- if (nonOrthogonality_)
+ if (nonOrthogonality_ && !nonOrthPtr)
{
- if (!nonOrthPtr)
- {
- IOobject fieldHeader
- (
- nonOrthogonalityName_,
- mesh_.time().constant(),
- mesh_,
- IOobject::MUST_READ,
- IOobject::NO_WRITE
- );
+ IOobject fieldHeader
+ (
+ nonOrthogonalityName_,
+ mesh_.time().constant(),
+ mesh_,
+ IOobject::MUST_READ,
+ IOobject::NO_WRITE
+ );
- if (fieldHeader.typeHeaderOk<volScalarField>(true, true, false))
- {
- volScalarField* vfPtr(new volScalarField(fieldHeader, mesh_));
- mesh_.objectRegistry::store(vfPtr);
- }
- else
- {
- FatalErrorInFunction
- << " Field : " << nonOrthogonalityName_ << " not found."
- << " The function object will not be used"
- << exit(FatalError);
- }
+ if (fieldHeader.typeHeaderOk<volScalarField>(true, true, false))
+ {
+ volScalarField* vfPtr = new volScalarField(fieldHeader, mesh_);
+ mesh_.objectRegistry::store(vfPtr);
+ }
+ else
+ {
+ FatalErrorInFunction
+ << "Field : " << nonOrthogonalityName_ << " not found."
+ << " The function object will not be used"
+ << exit(FatalError);
}
}
- const volScalarField* faceWeightsPtr =
+ const auto* faceWeightsPtr =
mesh_.findObject<volScalarField>(faceWeightName_);
- if (faceWeight_)
+ if (faceWeight_ && !faceWeightsPtr)
{
- if (!faceWeightsPtr)
- {
- IOobject fieldHeader
- (
- faceWeightName_,
- mesh_.time().constant(),
- mesh_,
- IOobject::MUST_READ,
- IOobject::NO_WRITE
- );
+ IOobject fieldHeader
+ (
+ faceWeightName_,
+ mesh_.time().constant(),
+ mesh_,
+ IOobject::MUST_READ,
+ IOobject::NO_WRITE
+ );
- if (fieldHeader.typeHeaderOk<volScalarField>(true, true, false))
- {
- volScalarField* vfPtr(new volScalarField(fieldHeader, mesh_));
- mesh_.objectRegistry::store(vfPtr);
- }
- else
- {
- FatalErrorInFunction
- << " Field : " << faceWeightName_ << " not found."
- << " The function object will not be used"
- << exit(FatalError);
- }
+ if (fieldHeader.typeHeaderOk<volScalarField>(true, true, false))
+ {
+ volScalarField* vfPtr = new volScalarField(fieldHeader, mesh_);
+ mesh_.objectRegistry::store(vfPtr);
+ }
+ else
+ {
+ FatalErrorInFunction
+ << "Field : " << faceWeightName_ << " not found."
+ << " The function object will not be used"
+ << exit(FatalError);
}
}
- const volScalarField* skewnessPtr =
- mesh_.findObject<volScalarField>(skewnessName_);
+ const auto* skewnessPtr = mesh_.findObject<volScalarField>(skewnessName_);
- if (skewness_)
+ if (skewness_ && !skewnessPtr)
{
- if (!skewnessPtr)
- {
- IOobject fieldHeader
- (
- skewnessName_,
- mesh_.time().constant(),
- mesh_,
- IOobject::MUST_READ,
- IOobject::NO_WRITE
- );
+ IOobject fieldHeader
+ (
+ skewnessName_,
+ mesh_.time().constant(),
+ mesh_,
+ IOobject::MUST_READ,
+ IOobject::NO_WRITE
+ );
- if (fieldHeader.typeHeaderOk<volScalarField>(true, true, false))
- {
- volScalarField* vfPtr(new volScalarField(fieldHeader, mesh_));
- mesh_.objectRegistry::store(vfPtr);
- }
- else
- {
- FatalErrorInFunction
- << " Field : " << skewnessName_ << " not found."
- << " The function object will not be used"
- << exit(FatalError);
- }
+ if (fieldHeader.typeHeaderOk<volScalarField>(true, true, false))
+ {
+ volScalarField* vfPtr(new volScalarField(fieldHeader, mesh_));
+ mesh_.objectRegistry::store(vfPtr);
+ }
+ else
+ {
+ FatalErrorInFunction
+ << "Field : " << skewnessName_ << " not found."
+ << " The function object will not be used"
+ << exit(FatalError);
}
}
if (log)
{
- indicator_.writeOpt(IOobject::AUTO_WRITE);
-
+ indicator().writeOpt(IOobject::AUTO_WRITE);
}
if (writeToFile_)
diff --git a/src/functionObjects/field/stabilityBlendingFactor/stabilityBlendingFactor.H b/src/functionObjects/field/stabilityBlendingFactor/stabilityBlendingFactor.H
index 579ad13c0ba388ab54356f9fb3219fc6f95f344c..5823e654951f0c5c0094b35882ff2d0332c9a3b8 100644
--- a/src/functionObjects/field/stabilityBlendingFactor/stabilityBlendingFactor.H
+++ b/src/functionObjects/field/stabilityBlendingFactor/stabilityBlendingFactor.H
@@ -340,9 +340,6 @@ class stabilityBlendingFactor
{
// Private Member Data
- //- Cell based blended indicator
- volScalarField indicator_;
-
// Switches
//- Switch for non-orthogonality
@@ -443,6 +440,9 @@ class stabilityBlendingFactor
//- Init fields
bool init(bool first);
+ //- Return access to the indicator field
+ volScalarField& indicator();
+
//- Calculate statistics
void calcStats(label&, label&, label&) const ;