32#include "surfaceInterpolate.H"
38namespace functionObjects
54Foam::functionObjects::stabilityBlendingFactor::indicator()
56 const word fldName =
"blendedIndicator" +
fieldName_;
74 zeroGradientFvPatchScalarField::typeName
77 mesh_.objectRegistry::store(fldPtr);
84void Foam::functionObjects::stabilityBlendingFactor::calcStats
91 const auto* indicatorPtr =
92 mesh_.cfindObject<
volScalarField>(
"blendedIndicator" + fieldName_);
97 <<
"Indicator field not set"
101 const auto& indicator = *indicatorPtr;
103 for (
const scalar i : indicator)
109 else if (i > (1 - tolerance_))
119 reduce(nCellsScheme1, sumOp<label>());
120 reduce(nCellsScheme2, sumOp<label>());
121 reduce(nCellsBlended, sumOp<label>());
131 writeCommented(
os,
"Time");
132 writeTabbed(
os,
"Scheme1");
133 writeTabbed(
os,
"Scheme2");
134 writeTabbed(
os,
"Blended");
139bool Foam::functionObjects::stabilityBlendingFactor::calc()
148 const auto* residualPtr = mesh_.findObject<
IOField<scalar>>(residualName_);
150 auto& indicator = this->indicator();
157 <<
"Could not find residual field : " << residualName_
158 <<
". The residual mode will not be " <<
nl
159 <<
" considered for the blended field in the stability "
160 <<
"blending factor. " <<
nl
161 <<
" Please add the corresponding 'solverInfo'"
162 <<
" function object with 'writeResidualFields true'." <<
nl
163 <<
" If the solverInfo function object is already enabled "
164 <<
"you might need to wait " <<
nl
165 <<
" for the first iteration."
170 scalar meanRes =
gAverage(mag(*residualPtr)) + VSMALL;
173 <<
" Average(mag(residuals)) : " << meanRes <<
endl;
176 oldErrorIntegral_ = errorIntegral_;
177 error_ = mag(meanRes - mag(*residualPtr));
178 errorIntegral_ = oldErrorIntegral_ + 0.5*(error_ + oldError_);
179 const scalarField errorDifferential(error_ - oldError_);
185 + D_*errorDifferential
194 mag(factorList - meanRes)/(maxResidual_*meanRes),
203 indicator[i] = indicatorResidual[i];
211 if (nonOrthogonality_)
213 if (maxNonOrthogonality_ >= minNonOrthogonality_)
216 <<
"minNonOrthogonality should be larger than "
217 <<
"maxNonOrthogonality."
230 (*nonOrthPtr - maxNonOrthogonality_)
231 / (minNonOrthogonality_ - maxNonOrthogonality_)
239 Log <<
" Max non-orthogonality : " <<
max(*nonOrthPtr).value()
244 const auto* skewnessPtr = mesh_.cfindObject<
volScalarField>(skewnessName_);
248 if (maxSkewness_ >= minSkewness_)
251 <<
"minSkewness should be larger than maxSkewness."
264 (*skewnessPtr - maxSkewness_)
265 / (minSkewness_ - maxSkewness_)
273 Log <<
" Max skewness : " <<
max(*skewnessPtr).value()
278 const auto* faceWeightsPtr =
283 if (maxFaceWeight_ >= minFaceWeight_)
286 <<
"minFaceWeight should be larger than maxFaceWeight."
299 (minFaceWeight_ - *faceWeightsPtr)
300 / (minFaceWeight_ - maxFaceWeight_)
308 Log <<
" Min face weight: " <<
min(*faceWeightsPtr).value()
316 if (maxGradCc_ >= minGradCc_)
319 <<
"minGradCc should be larger than maxGradCc."
335 zeroGradientFvPatchScalarField::typeName
337 auto& magGradCC = tmagGradCC.ref();
347 mesh_.time().timeName(),
354 zeroGradientFvPatchScalarField::typeName
356 cci = mesh_.C().component(i);
357 cci.correctBoundaryConditions();
363 Log <<
" Max magGradCc : " <<
max(magGradCC.ref()).value()
376 (magGradCC - maxGradCc_)
377 / (minGradCc_ - maxGradCc_)
392 <<
"Co2 should be larger than Co1."
408 zeroGradientFvPatchScalarField::typeName
411 auto& Co = CoPtr.ref();
413 Co.primitiveFieldRef() =
414 mesh_.time().deltaT()*
mag(*UNamePtr)/
cbrt(mesh_.V());
420 min(
max(scalar(0), (Co - Co1_)/(Co2_ - Co1_)), scalar(1))
425 Log <<
" Max Co : " <<
max(Co).value()
437 label nCellsScheme1 = 0;
438 label nCellsScheme2 = 0;
439 label nCellsBlended = 0;
441 calcStats(nCellsScheme1, nCellsScheme2, nCellsBlended);
444 <<
" scheme 1 cells : " << nCellsScheme1 <<
nl
445 <<
" scheme 2 cells : " << nCellsScheme2 <<
nl
446 <<
" blended cells : " << nCellsBlended <<
nl
450 indicator.correctBoundaryConditions();
474 nonOrthogonality_(
dict.getOrDefault<
Switch>(
"switchNonOrtho", false)),
475 gradCc_(
dict.getOrDefault<
Switch>(
"switchGradCc", false)),
476 residuals_(
dict.getOrDefault<
Switch>(
"switchResiduals", false)),
477 faceWeight_(
dict.getOrDefault<
Switch>(
"switchFaceWeight", false)),
478 skewness_(
dict.getOrDefault<
Switch>(
"switchSkewness", false)),
479 Co_(
dict.getOrDefault<
Switch>(
"switchCo", false)),
483 dict.getOrDefault<scalar>(
"maxNonOrthogonality", 20.0)
487 dict.getOrDefault<scalar>(
"minNonOrthogonality", 60.0)
489 maxGradCc_(
dict.getOrDefault<scalar>(
"maxGradCc", 3.0)),
490 minGradCc_(
dict.getOrDefault<scalar>(
"minGradCc", 4.0)),
491 maxResidual_(
dict.getOrDefault<scalar>(
"maxResidual", 10.0)),
492 minFaceWeight_(
dict.getOrDefault<scalar>(
"minFaceWeight", 0.3)),
493 maxFaceWeight_(
dict.getOrDefault<scalar>(
"maxFaceWeight", 0.2)),
494 maxSkewness_(
dict.getOrDefault<scalar>(
"maxSkewness", 2.0)),
495 minSkewness_(
dict.getOrDefault<scalar>(
"minSkewness", 3.0)),
496 Co1_(
dict.getOrDefault<scalar>(
"Co1", 1.0)),
497 Co2_(
dict.getOrDefault<scalar>(
"Co2", 10.0)),
499 nonOrthogonalityName_
501 dict.getOrDefault<
word>(
"nonOrthogonality",
"nonOrthoAngle")
505 dict.getOrDefault<
word>(
"faceWeight",
"faceWeight")
509 dict.getOrDefault<
word>(
"skewness",
"skewness")
516 IOobject::scopedName(
"initialResidual",
"p")
525 error_(mesh_.nCells(),
Zero),
526 errorIntegral_(mesh_.nCells(),
Zero),
527 oldError_(mesh_.nCells(),
Zero),
528 oldErrorIntegral_(mesh_.nCells(),
Zero),
529 P_(
dict.getOrDefault<scalar>(
"P", 3)),
530 I_(
dict.getOrDefault<scalar>(
"I", 0.0)),
531 D_(
dict.getOrDefault<scalar>(
"D", 0.25))
551 const auto* nonOrthPtr =
554 if (nonOrthogonality_ && !nonOrthPtr)
558 nonOrthogonalityName_,
568 mesh_.objectRegistry::store(vfPtr);
573 <<
"Field : " << nonOrthogonalityName_ <<
" not found."
574 <<
" The function object will not be used"
580 const auto* faceWeightsPtr =
583 if (faceWeight_ && !faceWeightsPtr)
597 mesh_.objectRegistry::store(vfPtr);
602 <<
"Field : " << faceWeightName_ <<
" not found."
603 <<
" The function object will not be used"
610 if (skewness_ && !skewnessPtr)
624 mesh_.objectRegistry::store(vfPtr);
629 <<
"Field : " << skewnessName_ <<
" not found."
630 <<
" The function object will not be used"
658 dict.readEntry(
"switchNonOrtho", nonOrthogonality_);
659 dict.readEntry(
"switchGradCc", gradCc_);
660 dict.readEntry(
"switchResiduals", residuals_);
661 dict.readEntry(
"switchFaceWeight", faceWeight_);
662 dict.readEntry(
"switchSkewness", skewness_);
663 dict.readEntry(
"switchCo", Co_);
665 dict.readIfPresent(
"maxNonOrthogonality", maxNonOrthogonality_);
666 dict.readIfPresent(
"maxGradCc", maxGradCc_);
667 dict.readIfPresent(
"maxResidual", maxResidual_);
668 dict.readIfPresent(
"maxSkewness", maxSkewness_);
669 dict.readIfPresent(
"maxFaceWeight", maxFaceWeight_);
670 dict.readIfPresent(
"Co2", Co2_);
672 dict.readIfPresent(
"minFaceWeight", minFaceWeight_);
673 dict.readIfPresent(
"minNonOrthogonality", minNonOrthogonality_);
674 dict.readIfPresent(
"minGradCc", minGradCc_);
675 dict.readIfPresent(
"minSkewness", minSkewness_);
676 dict.readIfPresent(
"Co1", Co1_);
679 dict.readIfPresent(
"P", P_);
680 dict.readIfPresent(
"I", I_);
681 dict.readIfPresent(
"D", D_);
686 dict.readIfPresent(
"tolerance", tolerance_)
687 && (tolerance_ < 0 || tolerance_ > 1)
691 <<
"tolerance must be in the range 0 to 1. Supplied value: "
696 if (nonOrthogonality_)
698 Info<<
" Including nonOrthogonality between: "
699 << minNonOrthogonality_ <<
" and " << maxNonOrthogonality_
704 Info<<
" Including gradient between: "
705 << minGradCc_ <<
" and " << maxGradCc_ <<
endl;
709 Info<<
" Including residuals" <<
endl;
713 Info<<
" Including faceWeight between: "
714 << minFaceWeight_ <<
" and " << maxFaceWeight_ <<
endl;
718 Info<<
" Including skewness between: "
719 << minSkewness_ <<
" and " << maxSkewness_ <<
endl;
723 Info<<
" Including Co between: "
724 << Co2_ <<
" and " << Co1_ <<
endl;
736 label nCellsScheme1 = 0;
737 label nCellsScheme2 = 0;
738 label nCellsBlended = 0;
740 calcStats(nCellsScheme1, nCellsScheme2, nCellsBlended);
744 writeCurrentTime(file());
747 <<
tab << nCellsScheme1
748 <<
tab << nCellsScheme2
749 <<
tab << nCellsBlended
Macros for easy insertion into run-time selection tables.
#define addToRunTimeSelectionTable(baseType, thisType, argNames)
Add to construction table with typeName as the key.
A primitive field of type <T> with automated input and output.
Defines the attributes of an object for which implicit objectRegistry management is supported,...
bool typeHeaderOk(const bool checkType=true, const bool search=true, const bool verbose=true)
Read header (uses typeFilePath to find file) and check its info.
writeOption writeOpt() const noexcept
The write option.
An Ostream is an abstract base class for all output systems (streams, files, token lists,...
virtual bool read()
Re-read model coefficients if they have changed.
A simple wrapper around bool so that it can be read as a word: true/false, on/off,...
const word & constant() const
Return constant name.
Class to control time during OpenFOAM simulations that is also the top-level objectRegistry.
static autoPtr< Time > New()
Construct (dummy) Time - no functionObjects or libraries.
static word timeName(const scalar t, const int precision=precision_)
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Abstract base-class for Time/database function objects.
Intermediate class for handling field expression function objects (e.g. blendingFactor etc....
word resultName_
Name of result field.
word fieldName_
Name of field to process.
void setResultName(const word &typeName, const word &defaultArg)
Set the name of result field.
const fvMesh & mesh_
Reference to the fvMesh.
Computes the natural logarithm of an input volScalarField.
bool store(word &fieldName, const tmp< ObjectType > &tfield, bool cacheable=false)
Store the field in the (sub) objectRegistry under the given name.
Computes the stabilityBlendingFactor to be used by the local blended convection scheme....
virtual void writeFileHeader(Ostream &os) const
Write the file header.
virtual bool write()
Write the stabilityBlendingFactor.
virtual bool read(const dictionary &)
Read the stabilityBlendingFactor data.
const Time & time_
Reference to the time database.
Base class for writing single files from the function objects.
bool writeToFile_
Flag to enable/disable writing to file.
virtual OFstream & file()
Return access to the file (if only 1)
const Time & time() const
Return the top-level database.
const Type * findObject(const word &name, const bool recursive=false) const
Return const pointer to the object of the given Type.
Type * getObjectPtr(const word &name, const bool recursive=false) const
static const char *const componentNames[]
static constexpr direction nComponents
Number of components in bool is 1.
A class for handling words, derived from Foam::string.
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
OBJstream os(runTime.globalPath()/outputName)
Calculate the gradient of the given field.
#define WarningInFunction
Report a warning using Foam::Warning.
static tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > interpolate(const GeometricField< Type, fvPatchField, volMesh > &tvf, const surfaceScalarField &faceFlux, Istream &schemeData)
Interpolate field onto faces using scheme given by Istream.
tmp< GeometricField< typename outerProduct< vector, Type >::type, fvPatchField, volMesh > > grad(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
GeometricField< vector, fvPatchField, volMesh > volVectorField
const dimensionSet dimless
Dimensionless.
const dimensionSet dimLength(0, 1, 0, 0, 0, 0, 0)
GeometricField< scalar, fvPatchField, volMesh > volScalarField
static void writeHeader(Ostream &os, const word &fieldName)
messageStream Info
Information stream (stdout output on master, null elsewhere)
GeometricField< scalar, fvsPatchField, surfaceMesh > surfaceScalarField
fileName::Type type(const fileName &name, const bool followLink=true)
Return the file type: DIRECTORY or FILE, normally following symbolic links.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Ostream & endl(Ostream &os)
Add newline and flush stream.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
void reduce(const List< UPstream::commsStruct > &comms, T &value, const BinaryOp &bop, const int tag, const label comm)
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
errorManip< error > abort(error &err)
Type gAverage(const FieldField< Field, Type > &f)
static constexpr const zero Zero
Global zero (0)
dimensionedScalar cbrt(const dimensionedScalar &ds)
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
errorManipArg< error, int > exit(error &err, const int errNo=1)
constexpr char nl
The newline '\n' character (0x0a)
constexpr char tab
The tab '\t' character(0x09)
#define forAll(list, i)
Loop across all elements in list.