Go to the documentation of this file.
32 #include "surfaceInterpolate.H"
38 namespace functionObjects
44 stabilityBlendingFactor,
53 void Foam::functionObjects::stabilityBlendingFactor::calcStats
62 scalar i = indicator_[celli];
68 else if (i > (1 - tolerance_))
78 reduce(nCellsScheme1, sumOp<label>());
79 reduce(nCellsScheme2, sumOp<label>());
80 reduce(nCellsBlended, sumOp<label>());
90 writeCommented(os,
"Time");
91 writeTabbed(os,
"Scheme1");
92 writeTabbed(os,
"Scheme2");
93 writeTabbed(os,
"Blended");
98 bool Foam::functionObjects::stabilityBlendingFactor::calc()
105 bool Foam::functionObjects::stabilityBlendingFactor::init(
bool first)
107 const auto* residualPtr = mesh_.findObject<
IOField<scalar>>(residualName_);
114 <<
"Could not find residual field : " << residualName_
115 <<
". The residual mode won't be " <<
nl
116 <<
" considered for the blended field in the stability "
117 <<
"blending factor. " <<
nl
118 <<
" Please add the corresponding 'solverInfo'"
119 <<
" function object with 'writeResidualFields true'." <<
nl
120 <<
" If the solverInfo function object is already enabled "
121 <<
"you might need to wait " <<
nl
122 <<
" for the first iteration."
127 scalar meanRes =
gAverage(
mag(*residualPtr)) + VSMALL;
130 <<
" Average(mag(residuals)) : " << meanRes <<
endl;
133 oldErrorIntegral_ = errorIntegral_;
134 error_ =
mag(meanRes -
mag(*residualPtr));
135 errorIntegral_ = oldErrorIntegral_ + 0.5*(error_ + oldError_);
136 const scalarField errorDifferential(error_ - oldError_);
142 + D_*errorDifferential
151 mag(factorList - meanRes)/(maxResidual_*meanRes),
160 indicator_[i] = indicatorResidual[i];
168 if (nonOrthogonality_)
170 if (maxNonOrthogonality_ >= minNonOrthogonality_)
173 <<
" minNonOrthogonality should be larger than "
174 <<
"maxNonOrthogonality."
187 (*nonOrthPtr - maxNonOrthogonality_)
188 / (minNonOrthogonality_ - maxNonOrthogonality_)
196 Log <<
" Max non-orthogonality : " <<
max(*nonOrthPtr).value()
206 if (maxSkewness_ >= minSkewness_)
209 <<
" minSkewness should be larger than maxSkewness."
222 (*skewnessPtr - maxSkewness_)
223 / (minSkewness_ - maxSkewness_)
231 Log <<
" Max skewness : " <<
max(*skewnessPtr).value()
241 if (maxFaceWeight_ >= minFaceWeight_)
244 <<
" minFaceWeight should be larger than maxFaceWeight."
257 (minFaceWeight_ - *faceWeightsPtr)
258 / (minFaceWeight_ - maxFaceWeight_)
266 Log <<
" Min face weight: " <<
min(*faceWeightsPtr).value()
274 if (maxGradCc_ >= minGradCc_)
277 <<
" minGradCc should be larger than maxGradCc."
293 zeroGradientFvPatchScalarField::typeName
295 auto& magGradCC = tmagGradCC.ref();
305 mesh_.time().timeName(),
312 zeroGradientFvPatchScalarField::typeName
315 cci.correctBoundaryConditions();
321 Log <<
" Max magGradCc : " <<
max(magGradCC.ref()).value()
334 (magGradCC - maxGradCc_)
335 / (minGradCc_ - maxGradCc_)
351 <<
" Co2 should be larger than Co1."
367 zeroGradientFvPatchScalarField::typeName
370 auto& Co = CoPtr.ref();
372 Co.primitiveFieldRef() =
373 mesh_.time().deltaT()*
mag(*UNamePtr)/
cbrt(mesh_.V());
379 min(
max(scalar(0), (Co - Co1_)/(Co2_ - Co1_)), scalar(1))
384 Log <<
" Max Co : " <<
max(Co).value()
396 label nCellsScheme1 = 0;
397 label nCellsScheme2 = 0;
398 label nCellsBlended = 0;
400 calcStats(nCellsScheme1, nCellsScheme2, nCellsBlended);
403 <<
" scheme 1 cells : " << nCellsScheme1 <<
nl
404 <<
" scheme 2 cells : " << nCellsScheme2 <<
nl
405 <<
" blended cells : " << nCellsBlended <<
nl
409 indicator_.correctBoundaryConditions();
438 "blendedIndicator" + fieldName_,
446 zeroGradientFvPatchScalarField::typeName
448 nonOrthogonality_(
dict.getOrDefault<
Switch>(
"switchNonOrtho",
false)),
449 gradCc_(
dict.getOrDefault<
Switch>(
"switchGradCc",
false)),
450 residuals_(
dict.getOrDefault<
Switch>(
"switchResiduals",
false)),
451 faceWeight_(
dict.getOrDefault<
Switch>(
"switchFaceWeight",
false)),
452 skewness_(
dict.getOrDefault<
Switch>(
"switchSkewness",
false)),
453 Co_(
dict.getOrDefault<
Switch>(
"switchCo",
false)),
457 dict.getOrDefault<scalar>(
"maxNonOrthogonality", 20.0)
461 dict.getOrDefault<scalar>(
"minNonOrthogonality", 60.0)
463 maxGradCc_(
dict.getOrDefault<scalar>(
"maxGradCc", 3.0)),
464 minGradCc_(
dict.getOrDefault<scalar>(
"minGradCc", 4.0)),
465 maxResidual_(
dict.getOrDefault<scalar>(
"maxResidual", 10.0)),
466 minFaceWeight_(
dict.getOrDefault<scalar>(
"minFaceWeight", 0.3)),
467 maxFaceWeight_(
dict.getOrDefault<scalar>(
"maxFaceWeight", 0.2)),
468 maxSkewness_(
dict.getOrDefault<scalar>(
"maxSkewness", 2.0)),
469 minSkewness_(
dict.getOrDefault<scalar>(
"minSkewness", 3.0)),
470 Co1_(
dict.getOrDefault<scalar>(
"Co1", 1.0)),
471 Co2_(
dict.getOrDefault<scalar>(
"Co2", 10.0)),
473 nonOrthogonalityName_
475 dict.getOrDefault<
word>(
"nonOrthogonality",
"nonOrthoAngle")
479 dict.getOrDefault<
word>(
"faceWeight",
"faceWeight")
483 dict.getOrDefault<
word>(
"skewness",
"skewness")
487 dict.getOrDefault<
word>(
"residual",
"initialResidual:p")
495 error_(mesh_.nCells(),
Zero),
496 errorIntegral_(mesh_.nCells(),
Zero),
497 oldError_(mesh_.nCells(),
Zero),
498 oldErrorIntegral_(mesh_.nCells(),
Zero),
499 P_(
dict.getOrDefault<scalar>(
"P", 3)),
500 I_(
dict.getOrDefault<scalar>(
"I", 0.0)),
501 D_(
dict.getOrDefault<scalar>(
"D", 0.25))
504 setResultName(typeName,
"");
519 store(resultName_, faceBlendedPtr);
524 if (nonOrthogonality_)
530 nonOrthogonalityName_,
531 mesh_.time().constant(),
540 mesh_.objectRegistry::store(vfPtr);
545 <<
" Field : " << nonOrthogonalityName_ <<
" not found."
546 <<
" The function object will not be used"
563 mesh_.time().constant(),
572 mesh_.objectRegistry::store(vfPtr);
577 <<
" Field : " << faceWeightName_ <<
" not found."
578 <<
" The function object will not be used"
594 mesh_.time().constant(),
603 mesh_.objectRegistry::store(vfPtr);
608 <<
" Field : " << skewnessName_ <<
" not found."
609 <<
" The function object will not be used"
623 writeFileHeader(file());
639 dict.readEntry(
"switchNonOrtho", nonOrthogonality_);
640 dict.readEntry(
"switchGradCc", gradCc_);
641 dict.readEntry(
"switchResiduals", residuals_);
642 dict.readEntry(
"switchFaceWeight", faceWeight_);
643 dict.readEntry(
"switchSkewness", skewness_);
644 dict.readEntry(
"switchCo", Co_);
646 dict.readIfPresent(
"maxNonOrthogonality", maxNonOrthogonality_);
647 dict.readIfPresent(
"maxGradCc", maxGradCc_);
648 dict.readIfPresent(
"maxResidual", maxResidual_);
649 dict.readIfPresent(
"maxSkewness", maxSkewness_);
650 dict.readIfPresent(
"maxFaceWeight", maxFaceWeight_);
651 dict.readIfPresent(
"Co2", Co2_);
653 dict.readIfPresent(
"minFaceWeight", minFaceWeight_);
654 dict.readIfPresent(
"minNonOrthogonality", minNonOrthogonality_);
655 dict.readIfPresent(
"minGradCc", minGradCc_);
656 dict.readIfPresent(
"minSkewness", minSkewness_);
657 dict.readIfPresent(
"Co1", Co1_);
660 dict.readIfPresent(
"P", P_);
661 dict.readIfPresent(
"I", I_);
662 dict.readIfPresent(
"D", D_);
667 dict.readIfPresent(
"tolerance", tolerance_)
668 && (tolerance_ < 0 || tolerance_ > 1)
672 <<
"tolerance must be in the range 0 to 1. Supplied value: "
677 if (nonOrthogonality_)
679 Info<<
" Including nonOrthogonality between: "
680 << minNonOrthogonality_ <<
" and " << maxNonOrthogonality_
685 Info<<
" Including gradient between: "
686 << minGradCc_ <<
" and " << maxGradCc_ <<
endl;
690 Info<<
" Including residuals" <<
endl;
694 Info<<
" Including faceWeight between: "
695 << minFaceWeight_ <<
" and " << maxFaceWeight_ <<
endl;
699 Info<<
" Including skewness between: "
700 << minSkewness_ <<
" and " << maxSkewness_ <<
endl;
704 Info<<
" Including Co between: "
705 << Co2_ <<
" and " << Co1_ <<
endl;
717 label nCellsScheme1 = 0;
718 label nCellsScheme2 = 0;
719 label nCellsBlended = 0;
721 calcStats(nCellsScheme1, nCellsScheme2, nCellsBlended);
725 writeCurrentTime(file());
728 <<
tab << nCellsScheme1
729 <<
tab << nCellsScheme2
730 <<
tab << nCellsBlended
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Defines the attributes of an object for which implicit objectRegistry management is supported,...
tmp< GeometricField< cmptType, PatchField, GeoMesh > > component(const direction) const
Return a component of the field.
A simple wrapper around bool so that it can be read as a word: true/false, on/off,...
const dimensionSet dimless(0, 0, 0, 0, 0, 0, 0)
Dimensionless.
tmp< GeometricField< typename outerProduct< vector, Type >::type, fvPatchField, volMesh >> grad(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Class to control time during OpenFOAM simulations that is also the top-level objectRegistry.
A class for handling words, derived from Foam::string.
A primitive field of type <T> with automated input and output.
static const char *const componentNames[]
const dimensionSet dimLength(0, 1, 0, 0, 0, 0, 0)
static constexpr const zero Zero
Global zero (0)
Type gAverage(const FieldField< Field, Type > &f)
virtual bool write()
Write the stabilityBlendingFactor.
bool read(const char *buf, int32_t &val)
Same as readInt32.
virtual bool read(const dictionary &dict)
Read the fieldExpression data.
Ostream & endl(Ostream &os)
Add newline and flush stream.
static void writeHeader(Ostream &os, const word &fieldName)
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
#define forAll(list, i)
Loop across all elements in list.
virtual bool read(const dictionary &dict)
Read.
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
messageStream Info
Information stream (uses stdout - output is on the master only)
word name(const complex &c)
Return string representation of complex.
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
GeometricField< scalar, fvPatchField, volMesh > volScalarField
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
dimensionedScalar log(const dimensionedScalar &ds)
Macros for easy insertion into run-time selection tables.
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.
errorManipArg< error, int > exit(error &err, const int errNo=1)
GeometricField< vector, fvPatchField, volMesh > volVectorField
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Intermediate class for handling field expression function objects (e.g. blendingFactor etc....
GeometricField< scalar, fvsPatchField, surfaceMesh > surfaceScalarField
addToRunTimeSelectionTable(functionObject, ObukhovLength, dictionary)
Calculate the gradient of the given field.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
fileName::Type type(const fileName &name, const bool followLink=true)
Return the file type: DIRECTORY or FILE, normally following symbolic links.
static tmp< T > New(Args &&... args)
Construct tmp of T with forwarding arguments.
defineTypeNameAndDebug(ObukhovLength, 0)
Computes the natural logarithm of an input volScalarField.
Base class for writing single files from the function objects.
stabilityBlendingFactor(const word &name, const Time &runTime, const dictionary &dict)
Construct from Time and dictionary.
An Ostream is an abstract base class for all output systems (streams, files, token lists,...
virtual bool read(const dictionary &)
Read the stabilityBlendingFactor data.
dimensionedScalar cbrt(const dimensionedScalar &ds)
#define WarningInFunction
Report a warning using Foam::Warning.
static constexpr direction nComponents
Number of components in this vector space.
virtual void writeFileHeader(Ostream &os) const
Write the file header.