32 #include "surfaceInterpolate.H" 38 namespace functionObjects
44 stabilityBlendingFactor,
54 Foam::functionObjects::stabilityBlendingFactor::indicator()
56 const word fldName =
"blendedIndicator" +
fieldName_;
78 mesh_.objectRegistry::store(fldPtr);
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." 338 auto& magGradCC = tmagGradCC.ref();
347 "cc" + word(vector::componentNames[i]),
348 mesh_.time().timeName(),
358 cci.correctBoundaryConditions();
364 Log <<
" Max magGradCc : " <<
max(magGradCC.ref()).value()
377 (magGradCC - maxGradCc_)
378 / (minGradCc_ - maxGradCc_)
393 <<
"Co2 should be larger than Co1." 412 auto& Co = CoPtr.ref();
414 Co.primitiveFieldRef() =
415 mesh_.time().deltaT()*
mag(*UNamePtr)/
cbrt(mesh_.V());
421 clamp((Co - Co1_)/(Co2_ - Co1_), zero_one{})
426 Log <<
" Max Co : " <<
max(Co).value()
438 label nCellsScheme1 = 0;
439 label nCellsScheme2 = 0;
440 label nCellsBlended = 0;
442 calcStats(nCellsScheme1, nCellsScheme2, nCellsBlended);
445 <<
" scheme 1 cells : " << nCellsScheme1 <<
nl 446 <<
" scheme 2 cells : " << nCellsScheme2 <<
nl 447 <<
" blended cells : " << nCellsBlended <<
nl 451 indicator.correctBoundaryConditions();
452 indicator.clamp_range(zero_one{});
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))
552 const auto* nonOrthPtr =
555 if (nonOrthogonality_ && !nonOrthPtr)
559 nonOrthogonalityName_,
570 mesh_.objectRegistry::store(vfPtr);
575 <<
"Field : " << nonOrthogonalityName_ <<
" not found." 576 <<
" The function object will not be used" 582 const auto* faceWeightsPtr =
585 if (faceWeight_ && !faceWeightsPtr)
600 mesh_.objectRegistry::store(vfPtr);
605 <<
"Field : " << faceWeightName_ <<
" not found." 606 <<
" The function object will not be used" 613 if (skewness_ && !skewnessPtr)
628 mesh_.objectRegistry::store(vfPtr);
633 <<
"Field : " << skewnessName_ <<
" not found." 634 <<
" The function object will not be used" 662 dict.readEntry(
"switchNonOrtho", nonOrthogonality_);
663 dict.readEntry(
"switchGradCc", gradCc_);
664 dict.readEntry(
"switchResiduals", residuals_);
665 dict.readEntry(
"switchFaceWeight", faceWeight_);
666 dict.readEntry(
"switchSkewness", skewness_);
667 dict.readEntry(
"switchCo", Co_);
669 dict.readIfPresent(
"maxNonOrthogonality", maxNonOrthogonality_);
670 dict.readIfPresent(
"maxGradCc", maxGradCc_);
671 dict.readIfPresent(
"maxResidual", maxResidual_);
672 dict.readIfPresent(
"maxSkewness", maxSkewness_);
673 dict.readIfPresent(
"maxFaceWeight", maxFaceWeight_);
674 dict.readIfPresent(
"Co2", Co2_);
676 dict.readIfPresent(
"minFaceWeight", minFaceWeight_);
677 dict.readIfPresent(
"minNonOrthogonality", minNonOrthogonality_);
678 dict.readIfPresent(
"minGradCc", minGradCc_);
679 dict.readIfPresent(
"minSkewness", minSkewness_);
680 dict.readIfPresent(
"Co1", Co1_);
683 dict.readIfPresent(
"P", P_);
684 dict.readIfPresent(
"I", I_);
685 dict.readIfPresent(
"D", D_);
690 dict.readIfPresent(
"tolerance", tolerance_)
691 && (tolerance_ < 0 || tolerance_ > 1)
695 <<
"tolerance must be in the range 0 to 1. Supplied value: " 700 if (nonOrthogonality_)
702 Info<<
" Including nonOrthogonality between: " 703 << minNonOrthogonality_ <<
" and " << maxNonOrthogonality_
708 Info<<
" Including gradient between: " 709 << minGradCc_ <<
" and " << maxGradCc_ <<
endl;
713 Info<<
" Including residuals" <<
endl;
717 Info<<
" Including faceWeight between: " 718 << minFaceWeight_ <<
" and " << maxFaceWeight_ <<
endl;
722 Info<<
" Including skewness between: " 723 << minSkewness_ <<
" and " << maxSkewness_ <<
endl;
727 Info<<
" Including Co between: " 728 << Co2_ <<
" and " << Co1_ <<
endl;
740 label nCellsScheme1 = 0;
741 label nCellsScheme2 = 0;
742 label nCellsBlended = 0;
744 calcStats(nCellsScheme1, nCellsScheme2, nCellsBlended);
748 writeCurrentTime(file());
751 <<
tab << nCellsScheme1
752 <<
tab << nCellsScheme2
753 <<
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.
tmp< GeometricField< cmptType, PatchField, GeoMesh > > component(const direction) const
Return a component of the field.
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.
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.
#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.
static tmp< T > New(Args &&... args)
Construct tmp with forwarding arguments.
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)
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()
void reduce(const List< UPstream::commsStruct > &comms, T &value, const BinaryOp &bop, const int tag, const label comm)
Reduce inplace (cf. MPI Allreduce) using specified communication schedule.
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.
dimensionSet clamp(const dimensionSet &a, const dimensionSet &range)
static constexpr const zero Zero
Global zero (0)
const Time & time_
Reference to the time database.