32 #include "surfaceInterpolate.H" 38 namespace functionObjects
44 stabilityBlendingFactor,
54 Foam::functionObjects::stabilityBlendingFactor::indicator()
56 const word fldName =
"blendedIndicator" +
fieldName_;
85 void Foam::functionObjects::stabilityBlendingFactor::calcStats
92 const auto* indicatorPtr =
93 mesh_.cfindObject<
volScalarField>(
"blendedIndicator" + fieldName_);
98 <<
"Indicator field not set" 102 const auto& indicator = *indicatorPtr;
104 for (
const scalar i : indicator)
110 else if (i > (1 - tolerance_))
132 writeCommented(
os,
"Time");
133 writeTabbed(
os,
"Scheme1");
134 writeTabbed(
os,
"Scheme2");
135 writeTabbed(
os,
"Blended");
140 bool Foam::functionObjects::stabilityBlendingFactor::calc()
147 bool Foam::functionObjects::stabilityBlendingFactor::init(
bool first)
149 const auto* residualPtr = mesh_.findObject<IOField<scalar>>(residualName_);
151 auto& indicator = this->indicator();
158 <<
"Could not find residual field : " << residualName_
159 <<
". The residual mode will not be " <<
nl 160 <<
" considered for the blended field in the stability " 161 <<
"blending factor. " <<
nl 162 <<
" Please add the corresponding 'solverInfo'" 163 <<
" function object with 'writeResidualFields true'." <<
nl 164 <<
" If the solverInfo function object is already enabled " 165 <<
"you might need to wait " <<
nl 166 <<
" for the first iteration." 171 scalar meanRes =
gAverage(
mag(*residualPtr)) + VSMALL;
174 <<
" Average(mag(residuals)) : " << meanRes <<
endl;
177 oldErrorIntegral_ = errorIntegral_;
178 error_ =
mag(meanRes -
mag(*residualPtr));
179 errorIntegral_ = oldErrorIntegral_ + 0.5*(error_ + oldError_);
180 const scalarField errorDifferential(error_ - oldError_);
186 + D_*errorDifferential
195 mag(factorList - meanRes)/(maxResidual_*meanRes),
204 indicator[i] = indicatorResidual[i];
212 if (nonOrthogonality_)
214 if (maxNonOrthogonality_ >= minNonOrthogonality_)
217 <<
"minNonOrthogonality should be larger than " 218 <<
"maxNonOrthogonality." 231 (*nonOrthPtr - maxNonOrthogonality_)
232 / (minNonOrthogonality_ - maxNonOrthogonality_)
240 Log <<
" Max non-orthogonality : " <<
max(*nonOrthPtr).value()
245 const auto* skewnessPtr = mesh_.cfindObject<
volScalarField>(skewnessName_);
249 if (maxSkewness_ >= minSkewness_)
252 <<
"minSkewness should be larger than maxSkewness." 265 (*skewnessPtr - maxSkewness_)
266 / (minSkewness_ - maxSkewness_)
274 Log <<
" Max skewness : " <<
max(*skewnessPtr).value()
279 const auto* faceWeightsPtr =
284 if (maxFaceWeight_ >= minFaceWeight_)
287 <<
"minFaceWeight should be larger than maxFaceWeight." 300 (minFaceWeight_ - *faceWeightsPtr)
301 / (minFaceWeight_ - maxFaceWeight_)
309 Log <<
" Min face weight: " <<
min(*faceWeightsPtr).value()
317 if (maxGradCc_ >= minGradCc_)
320 <<
"minGradCc should be larger than maxGradCc." 332 auto& magGradCC = tmagGradCC.ref();
339 "cc" + word(vector::componentNames[i]),
345 auto& cci = tcci.ref();
347 cci = mesh_.C().component(i);
348 cci.correctBoundaryConditions();
354 Log <<
" Max magGradCc : " <<
max(magGradCC.ref()).value()
367 (magGradCC - maxGradCc_)
368 / (minGradCc_ - maxGradCc_)
383 <<
"Co2 should be larger than Co1." 395 auto& Co = CoPtr.ref();
397 Co.primitiveFieldRef() =
398 mesh_.time().deltaT()*
mag(*UNamePtr)/
cbrt(mesh_.V());
404 clamp((Co - Co1_)/(Co2_ - Co1_), zero_one{})
409 Log <<
" Max Co : " <<
max(Co).value()
421 label nCellsScheme1 = 0;
422 label nCellsScheme2 = 0;
423 label nCellsBlended = 0;
425 calcStats(nCellsScheme1, nCellsScheme2, nCellsBlended);
428 <<
" scheme 1 cells : " << nCellsScheme1 <<
nl 429 <<
" scheme 2 cells : " << nCellsScheme2 <<
nl 430 <<
" blended cells : " << nCellsBlended <<
nl 434 indicator.correctBoundaryConditions();
435 indicator.clamp_range(zero_one{});
457 nonOrthogonality_(
dict.getOrDefault<
Switch>(
"switchNonOrtho", false)),
458 gradCc_(
dict.getOrDefault<
Switch>(
"switchGradCc", false)),
459 residuals_(
dict.getOrDefault<
Switch>(
"switchResiduals", false)),
460 faceWeight_(
dict.getOrDefault<
Switch>(
"switchFaceWeight", false)),
461 skewness_(
dict.getOrDefault<
Switch>(
"switchSkewness", false)),
462 Co_(
dict.getOrDefault<
Switch>(
"switchCo", false)),
466 dict.getOrDefault<scalar>(
"maxNonOrthogonality", 20.0)
470 dict.getOrDefault<scalar>(
"minNonOrthogonality", 60.0)
472 maxGradCc_(
dict.getOrDefault<scalar>(
"maxGradCc", 3.0)),
473 minGradCc_(
dict.getOrDefault<scalar>(
"minGradCc", 4.0)),
474 maxResidual_(
dict.getOrDefault<scalar>(
"maxResidual", 10.0)),
475 minFaceWeight_(
dict.getOrDefault<scalar>(
"minFaceWeight", 0.3)),
476 maxFaceWeight_(
dict.getOrDefault<scalar>(
"maxFaceWeight", 0.2)),
477 maxSkewness_(
dict.getOrDefault<scalar>(
"maxSkewness", 2.0)),
478 minSkewness_(
dict.getOrDefault<scalar>(
"minSkewness", 3.0)),
479 Co1_(
dict.getOrDefault<scalar>(
"Co1", 1.0)),
480 Co2_(
dict.getOrDefault<scalar>(
"Co2", 10.0)),
482 nonOrthogonalityName_
484 dict.getOrDefault<
word>(
"nonOrthogonality",
"nonOrthoAngle")
488 dict.getOrDefault<
word>(
"faceWeight",
"faceWeight")
492 dict.getOrDefault<
word>(
"skewness",
"skewness")
499 IOobject::scopedName(
"initialResidual",
"p")
508 error_(mesh_.nCells(),
Zero),
509 errorIntegral_(mesh_.nCells(),
Zero),
510 oldError_(mesh_.nCells(),
Zero),
511 oldErrorIntegral_(mesh_.nCells(),
Zero),
512 P_(
dict.getOrDefault<scalar>(
"P", 3)),
513 I_(
dict.getOrDefault<scalar>(
"I", 0.0)),
514 D_(
dict.getOrDefault<scalar>(
"D", 0.25))
528 const auto* nonOrthPtr =
531 if (nonOrthogonality_ && !nonOrthPtr)
535 nonOrthogonalityName_,
551 <<
"Field : " << nonOrthogonalityName_ <<
" not found." 552 <<
" The function object will not be used" 558 const auto* faceWeightsPtr =
561 if (faceWeight_ && !faceWeightsPtr)
581 <<
"Field : " << faceWeightName_ <<
" not found." 582 <<
" The function object will not be used" 589 if (skewness_ && !skewnessPtr)
609 <<
"Field : " << skewnessName_ <<
" not found." 610 <<
" The function object will not be used" 638 dict.readEntry(
"switchNonOrtho", nonOrthogonality_);
639 dict.readEntry(
"switchGradCc", gradCc_);
640 dict.readEntry(
"switchResiduals", residuals_);
641 dict.readEntry(
"switchFaceWeight", faceWeight_);
642 dict.readEntry(
"switchSkewness", skewness_);
643 dict.readEntry(
"switchCo", Co_);
645 dict.readIfPresent(
"maxNonOrthogonality", maxNonOrthogonality_);
646 dict.readIfPresent(
"maxGradCc", maxGradCc_);
647 dict.readIfPresent(
"maxResidual", maxResidual_);
648 dict.readIfPresent(
"maxSkewness", maxSkewness_);
649 dict.readIfPresent(
"maxFaceWeight", maxFaceWeight_);
650 dict.readIfPresent(
"Co2", Co2_);
652 dict.readIfPresent(
"minFaceWeight", minFaceWeight_);
653 dict.readIfPresent(
"minNonOrthogonality", minNonOrthogonality_);
654 dict.readIfPresent(
"minGradCc", minGradCc_);
655 dict.readIfPresent(
"minSkewness", minSkewness_);
656 dict.readIfPresent(
"Co1", Co1_);
659 dict.readIfPresent(
"P", P_);
660 dict.readIfPresent(
"I", I_);
661 dict.readIfPresent(
"D", D_);
666 dict.readIfPresent(
"tolerance", tolerance_)
667 && (tolerance_ < 0 || tolerance_ > 1)
671 <<
"tolerance must be in the range 0 to 1. Supplied value: " 676 if (nonOrthogonality_)
678 Info<<
" Including nonOrthogonality between: " 679 << minNonOrthogonality_ <<
" and " << maxNonOrthogonality_
684 Info<<
" Including gradient between: " 685 << minGradCc_ <<
" and " << maxGradCc_ <<
endl;
689 Info<<
" Including residuals" <<
endl;
693 Info<<
" Including faceWeight between: " 694 << minFaceWeight_ <<
" and " << maxFaceWeight_ <<
endl;
698 Info<<
" Including skewness between: " 699 << minSkewness_ <<
" and " << maxSkewness_ <<
endl;
703 Info<<
" Including Co between: " 704 << Co2_ <<
" and " << Co1_ <<
endl;
716 label nCellsScheme1 = 0;
717 label nCellsScheme2 = 0;
718 label nCellsBlended = 0;
720 calcStats(nCellsScheme1, nCellsScheme2, nCellsBlended);
724 writeCurrentTime(file());
727 <<
tab << nCellsScheme1
728 <<
tab << nCellsScheme2
729 <<
tab << nCellsBlended
static const word & zeroGradientType() noexcept
The type name for zeroGradient patch fields.
virtual OFstream & file()
Return access to the file (if only 1)
tmp< GeometricField< typename outerProduct< vector, Type >::type, fvPatchField, volMesh >> grad(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
defineTypeNameAndDebug(ObukhovLength, 0)
word resultName_
Name of result field.
static void writeHeader(Ostream &os, const word &fieldName)
dimensionedScalar log(const dimensionedScalar &ds)
errorManipArg< error, int > exit(error &err, const int errNo=1)
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
error FatalError
Error stream (stdout output on all processes), with additional 'FOAM FATAL ERROR' header text and sta...
A list of keyword definitions, which are a keyword followed by a number of values (eg...
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
word fieldName_
Name of field to process.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
constexpr char nl
The newline '\n' character (0x0a)
Ostream & endl(Ostream &os)
Add newline and flush stream.
virtual bool read(const dictionary &dict)
Read the fieldExpression data.
A simple wrapper around bool so that it can be read as a word: true/false, on/off, yes/no, any/none. Also accepts 0/1 as a string and shortcuts t/f, y/n.
virtual void writeFileHeader(Ostream &os) const
Write the file header.
bool store()
Register object with its registry and transfer ownership to the registry.
constexpr char tab
The tab '\t' character(0x09)
const Type * findObject(const word &name, const bool recursive=false) const
Return const pointer to the object of the given Type.
::Foam::direction nComponents(const expressions::valueTypeCode) noexcept
The number of components associated with given valueTypeCode.
Ignore writing from objectRegistry::writeObject()
const dimensionSet dimless
Dimensionless.
const Time & time() const
Return the top-level database.
GeometricField< vector, fvPatchField, volMesh > volVectorField
Class to control time during OpenFOAM simulations that is also the top-level objectRegistry.
Macros for easy insertion into run-time selection tables.
virtual const objectRegistry & thisDb() const
Return the object registry - resolve conflict polyMesh/lduMesh.
#define forAll(list, i)
Loop across all elements in list.
GeometricField< scalar, fvPatchField, volMesh > volScalarField
void setResultName(const word &typeName, const word &defaultArg)
Set the name of result field.
virtual bool read(const dictionary &)
Read the stabilityBlendingFactor data.
fileName::Type type(const fileName &name, const bool followLink=true)
Return the file type: DIRECTORY or FILE, normally following symbolic links.
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for expressions::valueTypeCode::INVALID.
Calculate the gradient of the given field.
A class for handling words, derived from Foam::string.
Type * getObjectPtr(const word &name, const bool recursive=false) const
Return non-const pointer to the object of the given Type, using a const-cast to have it behave like a...
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
dimensionedScalar cbrt(const dimensionedScalar &ds)
bool store(word &fieldName, const tmp< ObjectType > &tfield, bool cacheable=false)
Store the field in the (sub) objectRegistry under the given name.
stabilityBlendingFactor(const word &name, const Time &runTime, const dictionary &dict)
Construct from Time and dictionary.
bool log
Flag to write log into Info.
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
errorManip< error > abort(error &err)
static tmp< GeometricField< scalar, fvPatchField, volMesh > > New(const word &name, IOobjectOption::registerOption regOpt, const Mesh &mesh, const dimensionSet &dims, const word &patchFieldType=fvPatchField< scalar >::calculatedType())
Return tmp field (NO_READ, NO_WRITE) from name, mesh, dimensions and patch type. [Takes current timeN...
An Ostream is an abstract base class for all output systems (streams, files, token lists...
static word timeName(const scalar t, const int precision=precision_)
Return a time name for the given scalar time value formatted with the given precision.
const word & constant() const noexcept
Return constant name.
OBJstream os(runTime.globalPath()/outputName)
addToRunTimeSelectionTable(functionObject, ObukhovLength, dictionary)
virtual bool read(const dictionary &dict)
Read.
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.
bool writeToFile_
Flag to enable/disable writing to file.
#define WarningInFunction
Report a warning using Foam::Warning.
const dimensionSet dimLength(0, 1, 0, 0, 0, 0, 0)
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
Type gAverage(const FieldField< Field, Type > &f)
Intermediate class for handling field expression function objects (e.g. blendingFactor etc...
Automatically write from objectRegistry::writeObject()
messageStream Info
Information stream (stdout output on master, null elsewhere)
virtual bool write()
Write the stabilityBlendingFactor.
GeometricField< scalar, fvsPatchField, surfaceMesh > surfaceScalarField
Base class for writing single files from the function objects.
Defines the attributes of an object for which implicit objectRegistry management is supported...
Request registration (bool: true)
const fvMesh & mesh_
Reference to the fvMesh.
Do not request registration (bool: false)
void reduce(T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Reduce inplace (cf. MPI Allreduce) using linear/tree communication schedule.
dimensionSet clamp(const dimensionSet &a, const dimensionSet &range)
static constexpr const zero Zero
Global zero (0)
const Time & time_
Reference to the time database.