ENH: support exprField specification for SemiImplicitSource

- this allows more flexibility when defining the location or intensity
  of sources.

  For example,

  {
      type            scalarSemiImplicitSource;
      volumeMode      specific;
      selectionMode   all;

      sources
      {
          tracer0
          {
              explicit
              {
                  type       exprField;

                  functions<scalar>
                  {
                      square
                      {
                          type square;
                          scale 0.0025;
                          level 0.0025;
                          frequency 10;
                      }
                  }

                  expression
                  #{
                      (hypot(pos().x() + 0.025, pos().y()) < 0.01)
                    ? fn:square(time())
                    : 0
                  #};
              }
          }
      }
  }

ENH: SemiImplicitSource: handle "sources" with explicit/implicit entries

- essentially the same as injectionRateSuSp with Su/Sp,
  but potentially clearer in purpose.

ENH: add Function1 good() method to define if function can be evaluated

- for example, provides a programmatic means of avoiding the 'none'
  function
This commit is contained in:
Mark Olesen 2022-05-30 13:27:22 +02:00
parent ef743147ea
commit d2e10bca40
14 changed files with 575 additions and 107 deletions

View File

@ -243,6 +243,9 @@ public:
//- Is value constant (i.e. independent of x)
virtual bool constant() const { return false; }
//- Can function be evaluated?
virtual bool good() const { return true; }
// Evaluation

View File

@ -5,7 +5,7 @@
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2021 OpenCFD Ltd.
Copyright (C) 2021-2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -34,8 +34,8 @@ Description
\*---------------------------------------------------------------------------*/
#ifndef Function1Types_None_H
#define Function1Types_None_H
#ifndef Foam_Function1Types_None_H
#define Foam_Function1Types_None_H
#include "Function1.H"
@ -102,6 +102,9 @@ public:
//- Value is independent of x
virtual inline bool constant() const { return true; }
//- Function cannot be evaluated
virtual inline bool good() const { return false; }
//- Placeholder: generates an error if called
virtual Type value(const scalar x) const;

View File

@ -44,7 +44,7 @@ namespace Foam
void Foam::fa::option::resetApplied()
{
applied_.resize(fieldNames_.size());
applied_.resize_nocopy(fieldNames_.size());
applied_ = false;
}
@ -91,7 +91,7 @@ Foam::autoPtr<Foam::fa::option> Foam::fa::option::New
Info<< indent
<< "Selecting finite area options type " << modelType << endl;
const_cast<Time&>(patch.boundaryMesh().mesh().time()).libs().open
patch.boundaryMesh().mesh().time().libs().open
(
coeffs,
"libs",
@ -102,15 +102,16 @@ Foam::autoPtr<Foam::fa::option> Foam::fa::option::New
if (!ctorPtr)
{
FatalErrorInFunction
<< "Unknown faOption model type "
<< modelType << nl << nl
<< "Valid faOption types are:" << nl
<< dictionaryConstructorTablePtr_->sortedToc()
<< exit(FatalError);
FatalIOErrorInLookup
(
coeffs,
"faOption",
modelType,
*dictionaryConstructorTablePtr_
) << exit(FatalIOError);
}
return autoPtr<option>(ctorPtr(name, modelType, coeffs, patch));
return autoPtr<fa::option>(ctorPtr(name, modelType, coeffs, patch));
}

View File

@ -72,8 +72,8 @@ SourceFiles
\*---------------------------------------------------------------------------*/
#ifndef faOption_H
#define faOption_H
#ifndef Foam_faOption_H
#define Foam_faOption_H
#include "faMatricesFwd.H"
#include "areaFieldsFwd.H"
@ -263,7 +263,7 @@ public:
//- Return dictionary
inline const dictionary& coeffs() const noexcept;
//- Return const access to the source active flag
//- True if source is active
inline bool active() const noexcept;
//- Set the applied flag to true for field index fieldi

View File

@ -47,7 +47,7 @@ namespace Foam
void Foam::fv::option::resetApplied()
{
applied_.resize(fieldNames_.size());
applied_.resize_nocopy(fieldNames_.size());
applied_ = false;
}

View File

@ -68,14 +68,13 @@ SourceFiles
\*---------------------------------------------------------------------------*/
#ifndef fvOption_H
#define fvOption_H
#ifndef Foam_fvOption_H
#define Foam_fvOption_H
#include "fvMatricesFwd.H"
#include "primitiveFieldsFwd.H"
#include "volFieldsFwd.H"
#include "dictionary.H"
#include "Switch.H"
#include "runTimeSelectionTables.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -83,6 +82,7 @@ SourceFiles
namespace Foam
{
// Forward Declarations
class fvMesh;
namespace fv

View File

@ -100,14 +100,12 @@ Usage
type scalarSemiImplicitSource;
enabled true;
scalarSemiImplicitSourceCoeffs
selectionMode all;
volumeMode specific;
sources
{
selectionMode all;
volumeMode specific;
injectionRateSuSp
{
s (1 0);
}
s (1 0);
}
}
}

View File

@ -50,56 +50,159 @@ Foam::fv::SemiImplicitSource<Type>::volumeModeTypeNames_
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
template<class Type>
void Foam::fv::SemiImplicitSource<Type>::setFieldInjectionRates
void Foam::fv::SemiImplicitSource<Type>::setFieldCoeffs
(
const dictionary& dict
const dictionary& dict,
const word& keyExplicit,
const word& keyImplicit
)
{
label count = dict.size();
fieldNames_.resize_nocopy(count);
Su_.resize(count);
Sp_.resize(count);
Su_.clear();
Sp_.clear();
Su_.resize(2*count);
Sp_.resize(2*count);
driverSu_.clear();
driverSp_.clear();
driverSu_.resize(2*count);
driverSp_.resize(2*count);
valueExprSu_.clear();
valueExprSp_.clear();
valueExprSu_.resize(2*count);
valueExprSp_.resize(2*count);
fv::option::resetApplied();
word modelType;
Tuple2<Type, scalar> sourceRates;
count = 0;
for (const entry& dEntry : dict)
{
const word& fieldName = dEntry.keyword();
bool ok = false;
if (dEntry.isDict())
{
const dictionary& subdict = dEntry.dict();
Su_.set(count, Function1<Type>::New("Su", subdict, &mesh_));
Sp_.set(count, Function1<scalar>::New("Sp", subdict, &mesh_));
const dictionary& subDict = dEntry.dict();
const entry* eptr;
if
(
(eptr = subDict.findEntry(keyExplicit, keyType::LITERAL))
!= nullptr
)
{
ok = true;
if
(
eptr->isDict()
&& eptr->dict().readEntry("type", modelType, keyType::LITERAL)
&& (modelType == "exprField")
)
{
const dictionary& exprDict = eptr->dict();
valueExprSu_.emplace_set(fieldName);
valueExprSu_[fieldName].readEntry("expression", exprDict);
driverSu_.set
(
fieldName,
new expressions::volumeExprDriver(mesh_, exprDict)
);
}
else
{
Su_.set
(
fieldName,
Function1<Type>::New(keyExplicit, subDict, &mesh_)
);
}
}
if
(
(eptr = subDict.findEntry(keyImplicit, keyType::LITERAL))
!= nullptr
)
{
ok = true;
if
(
eptr->isDict()
&& eptr->dict().readEntry("type", modelType, keyType::LITERAL)
&& (modelType == "exprField")
)
{
const dictionary& exprDict = eptr->dict();
valueExprSp_.emplace_set(fieldName);
valueExprSp_[fieldName].readEntry("expression", exprDict);
driverSp_.set
(
fieldName,
new expressions::volumeExprDriver(mesh_, exprDict)
);
}
else
{
Sp_.set
(
fieldName,
Function1<scalar>::New(keyImplicit, subDict, &mesh_)
);
}
}
}
else
{
Tuple2<Type, scalar> injectionRate;
dEntry.readEntry(injectionRate);
// Non-dictionary form
dEntry.readEntry(sourceRates);
ok = true;
Su_.set
(
count,
fieldName,
new Function1Types::Constant<Type>
(
"Su",
injectionRate.first()
keyExplicit,
sourceRates.first()
)
);
Sp_.set
(
count,
fieldName,
new Function1Types::Constant<scalar>
(
"Sp",
injectionRate.second()
keyImplicit,
sourceRates.second()
)
);
}
if (!ok)
{
FatalIOErrorInFunction(dict)
<< "Require at least one of "
<< keyExplicit << '/' << keyImplicit << " entries for "
<< "field: " << fieldName << endl
<< exit(FatalIOError);
}
fieldNames_[count] = fieldName;
++count;
}
@ -154,6 +257,7 @@ void Foam::fv::SemiImplicitSource<Type>::addSup
const GeometricField<Type, fvPatchField, volMesh>& psi = eqn.psi();
// Note: field name may deviate from psi name
const word& fieldName = fieldNames_[fieldi];
const scalar tmVal = mesh_.time().timeOutputValue();
@ -161,67 +265,254 @@ void Foam::fv::SemiImplicitSource<Type>::addSup
const dimensionSet SuDims(eqn.dimensions()/dimVolume);
const dimensionSet SpDims(SuDims/psi.dimensions());
// Explicit source
{
const dimensioned<Type> SuValue
(
"Su",
SuDims,
Su_[fieldi].value(tmVal)/VDash_
);
const auto iter1 = valueExprSu_.cfind(fieldName);
const auto iter2 = Su_.cfind(fieldName);
if (mag(SuValue.value()) <= ROOTVSMALL)
tmp<DimensionedField<Type, volMesh>> tsu;
if (iter1.found())
{
// No-op
}
else if (this->useSubMesh())
{
auto tsu = DimensionedField<Type, volMesh>::New
const auto& valueExpr = iter1.val();
typedef
GeometricField<Type, fvPatchField, volMesh>
ExprResultType;
if (debug)
{
Info<< "Explicit expression source:" << nl
<< ">>>>" << nl
<< valueExpr.c_str() << nl
<< "<<<<" << nl;
}
auto& driver = *(driverSu_[fieldName]);
driver.clearVariables();
if (notNull(rho))
{
driver.addContextObject("rho", &rho);
}
// Switch dimension checking off
const bool oldDimChecking = dimensionSet::checking(false);
driver.parse(valueExpr);
// Restore dimension checking
dimensionSet::checking(oldDimChecking);
const ExprResultType* ptr = driver.isResultType<ExprResultType>();
if (!ptr)
{
FatalErrorInFunction
<< "Expression for Su " << fieldName
<< " evaluated to <" << driver.resultType()
<< "> but expected <" << ExprResultType::typeName
<< ">" << endl
<< exit(FatalError);
}
else if (ptr->size() != mesh_.nCells())
{
FatalErrorInFunction
<< "Expression for Su " << fieldName
<< " evaluated to " << ptr->size()
<< " instead of " << mesh_.nCells() << " values" << endl
<< exit(FatalError);
}
if (notNull(rho))
{
driver.removeContextObject(&rho);
}
const Field<Type>& exprFld = ptr->primitiveField();
tsu = DimensionedField<Type, volMesh>::New
(
name_ + fieldName + "Su",
mesh_,
dimensioned<Type>(SuDims, Zero)
);
UIndirectList<Type>(tsu.ref(), cells_) = SuValue.value();
eqn += tsu;
if (this->useSubMesh())
{
for (const label celli : cells_)
{
tsu.ref()[celli] = exprFld[celli]/VDash_;
}
}
else
{
tsu.ref().field() = exprFld;
if (!equal(VDash_, 1))
{
tsu.ref().field() /= VDash_;
}
}
}
else
else if (iter2.found() && iter2.val()->good())
{
eqn += SuValue;
const dimensioned<Type> SuValue
(
"Su",
SuDims,
iter2.val()->value(tmVal)/VDash_
);
if (mag(SuValue.value()) <= ROOTVSMALL)
{
// No-op
}
else if (this->useSubMesh())
{
tsu = DimensionedField<Type, volMesh>::New
(
name_ + fieldName + "Su",
mesh_,
dimensioned<Type>(SuDims, Zero)
);
UIndirectList<Type>(tsu.ref(), cells_) = SuValue.value();
}
else
{
eqn += SuValue;
}
}
if (tsu.valid())
{
eqn += tsu;
}
}
// Implicit source
{
const dimensioned<scalar> SpValue
(
"Sp",
SpDims,
Sp_[fieldi].value(tmVal)/VDash_
);
const auto iter1 = valueExprSp_.cfind(fieldName);
const auto iter2 = Sp_.cfind(fieldName);
if (mag(SpValue.value()) <= ROOTVSMALL)
tmp<DimensionedField<scalar, volMesh>> tsp;
if (iter1.found())
{
// No-op
}
else if (this->useSubMesh())
{
auto tsp = DimensionedField<scalar, volMesh>::New
const auto& valueExpr = iter1.val();
typedef volScalarField ExprResultType;
if (debug)
{
Info<< "Implicit expression source:" << nl
<< ">>>>" << nl
<< valueExpr.c_str() << nl
<< "<<<<" << nl;
}
auto& driver = *(driverSp_[fieldName]);
driver.clearVariables();
if (notNull(rho))
{
driver.addContextObject("rho", &rho);
}
// Switch dimension checking off
const bool oldDimChecking = dimensionSet::checking(false);
driver.parse(valueExpr);
// Restore dimension checking
dimensionSet::checking(oldDimChecking);
const ExprResultType* ptr = driver.isResultType<ExprResultType>();
if (!ptr)
{
FatalErrorInFunction
<< "Expression for Sp " << fieldName
<< " evaluated to <" << driver.resultType()
<< "> but expected <" << ExprResultType::typeName
<< ">" << endl
<< exit(FatalError);
}
else if (ptr->size() != mesh_.nCells())
{
FatalErrorInFunction
<< "Expression for Sp " << fieldName
<< " evaluated to " << ptr->size()
<< " instead of " << mesh_.nCells() << " values" << endl
<< exit(FatalError);
}
if (notNull(rho))
{
driver.removeContextObject(&rho);
}
const Field<scalar>& exprFld = ptr->primitiveField();
tsp = DimensionedField<scalar, volMesh>::New
(
name_ + fieldName + "Sp",
mesh_,
dimensioned<scalar>(SpDims, Zero)
);
UIndirectList<scalar>(tsp.ref(), cells_) = SpValue.value();
eqn += fvm::SuSp(tsp, psi);
if (this->useSubMesh())
{
for (const label celli : cells_)
{
tsp.ref()[celli] = exprFld[celli]/VDash_;
}
}
else
{
tsp.ref().field() = exprFld;
if (!equal(VDash_, 1))
{
tsp.ref().field() /= VDash_;
}
}
}
else
else if (iter2.found() && iter2.val()->good())
{
eqn += fvm::SuSp(SpValue, psi);
const dimensioned<scalar> SpValue
(
"Sp",
SpDims,
iter2.val()->value(tmVal)/VDash_
);
if (mag(SpValue.value()) <= ROOTVSMALL)
{
// No-op
}
else if (this->useSubMesh())
{
tsp = DimensionedField<scalar, volMesh>::New
(
name_ + fieldName + "Sp",
mesh_,
dimensioned<scalar>(SpDims, Zero)
);
UIndirectList<scalar>(tsp.ref(), cells_) = SpValue.value();
}
else
{
eqn += fvm::SuSp(SpValue, psi);
}
}
if (tsp.valid())
{
eqn += fvm::SuSp(tsp, psi);
}
}
}
@ -242,10 +533,26 @@ bool Foam::fv::SemiImplicitSource<Type>::read(const dictionary& dict)
VDash_ = V_;
}
// Compatibility (2112 and earlier)
const dictionary* injectDict =
coeffs_.findDict("injectionRateSuSp", keyType::LITERAL);
if (injectDict)
{
setFieldInjectionRates
setFieldCoeffs
(
coeffs_.subDict("injectionRateSuSp", keyType::LITERAL)
*injectDict,
"Su", // Su = explicit
"Sp" // Sp = implicit
);
}
else
{
setFieldCoeffs
(
coeffs_.subDict("sources", keyType::LITERAL),
"explicit",
"implicit"
);
}

View File

@ -33,7 +33,7 @@ Group
Description
Applies semi-implicit source within a specified region for \c Type,
where \c \<Type\>=Scalar/Vector/SphericalTensor/SymmTensor/Tensor.
The injection rate coefficients are specified
The source rate coefficients are specified
as pairs of Su-Sp coefficients, i.e.:
\f[
@ -57,27 +57,49 @@ Usage
// Mandatory entries (runtime modifiable)
volumeMode <volumeModeType>;
injectionRateSuSp
// Specification of sources (OpenFOAM-2206 and newer)
sources
{
// Specified as explicit(Su)/implicit(Sp) tuple:
k (30.7 0);
epsilon (1.5 0);
// The injectionRate can also be specified as a Function1
// by having dictionaries for the field entries instead
// Specified as Function1 or exprField
k
{
// Time-ramp from 0 to 30.7 at time 5
Su table
(
(0 0.0)
(5 30.7)
);
Sp 0.0;
explicit table ((0 0) (5 30.7));
implicit none;
}
epsilon
{
Su 1.5;
Sp 0.0;
explicit
{
type exprField:
expression "(mag(pos()) < 0.1) ? 1.5 : 0";
}
}
}
// Traditional specification of sources (OpenFOAM-2112 and older)
injectionRateSuSp
{
// Specified as explicit(Su)/implicit(Sp) tuple:
k (30.7 0);
epsilon (1.5 0);
// Specified as Function1
k
{
// Time-ramp from 0 to 30.7 at time 5
Su table ((0 0) (5 30.7));
Sp 0;
}
epsilon
{
Su 1.5;
Sp 0;
}
}
@ -92,7 +114,8 @@ Usage
type | Type name: \<Type\>SemiImplicitSource <!--
--> | word | yes | -
volumeMode | Volume mode type | word | yes | -
injectionRateSuSp | Injection rate settings | dictionary | yes | -
sources | Explicit/implicit sources | dict | cndtnl | -
injectionRateSuSp | Explicit/implicit sources | dict | cndtnl | -
\endtable
The inherited entries are elaborated in:
@ -105,6 +128,14 @@ Usage
specific | Values are given as \<quantity\>/m3
\endverbatim
Note
Missing explicit/implicit, Su/Sp entries are equivalent to constant values
of zero. However, at one entry must be supplied for the source terms.
Note
To use the \c exprField (expression fields) handling, the \em sources
dictionary form must be used.
See also
- Foam::fvOption
@ -119,6 +150,8 @@ SourceFiles
#include "cellSetOption.H"
#include "Enum.H"
#include "Function1.H"
#include "HashPtrTable.H"
#include "volumeExprDriver.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -127,12 +160,6 @@ namespace Foam
namespace fv
{
// Forward Declarations
template<class Type> class SemiImplicitSource;
template<class Type>
Ostream& operator<<(Ostream&, const SemiImplicitSource<Type>&);
/*---------------------------------------------------------------------------*\
Class SemiImplicitSource Declaration
\*---------------------------------------------------------------------------*/
@ -168,16 +195,34 @@ private:
scalar VDash_;
//- Explicit source contributions
PtrList<Function1<Type>> Su_;
HashPtrTable<Function1<Type>> Su_;
//- Linearised implicit contributions
PtrList<Function1<scalar>> Sp_;
HashPtrTable<Function1<scalar>> Sp_;
//- Expression to evaluate for explicit source contribution
HashTable<expressions::exprString> valueExprSu_;
//- Expression to evaluate for linearised implicit contribution
HashTable<expressions::exprString> valueExprSp_;
//- Expression driver for explicit sources
HashPtrTable<expressions::volumeExprDriver> driverSu_;
//- Expression driver for implicit sources
HashPtrTable<expressions::volumeExprDriver> driverSp_;
// Private Member Functions
//- Set the source coefficients from "injectionRateSuSp"
void setFieldInjectionRates(const dictionary& dict);
//- Set the source coefficients from "sources" (explicit/implicit)
//- or from "injectionRateSuSp" (Su/Sp)
void setFieldCoeffs
(
const dictionary& dict,
const word& keyExplicit,
const word& keyImplicit
);
public:

View File

@ -5,7 +5,7 @@
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2018-2021 OpenCFD Ltd.
Copyright (C) 2018-2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -218,6 +218,9 @@ public:
//- Is value constant (i.e. independent of x)
virtual bool constant() const { return false; }
//- Can function be evaluated?
virtual bool good() const { return true; }
//- Is value uniform (i.e. independent of coordinate)
virtual bool uniform() const = 0;

View File

@ -130,7 +130,13 @@ public:
//- Is value constant (i.e. independent of x)
virtual inline bool constant() const
{
return uniformValuePtr_ && uniformValuePtr_->constant();
return (uniformValuePtr_ && uniformValuePtr_->constant());
}
//- Can function be evaluated?
virtual inline bool good() const
{
return (uniformValuePtr_ && uniformValuePtr_->good());
}
//- Is value uniform (i.e. independent of coordinate)

View File

@ -0,0 +1,47 @@
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v2206 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
object tracer0;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 0 0 0 0 0 0];
internalField uniform 0;
boundaryField
{
#includeEtc "caseDicts/setConstraintTypes"
rotor
{
type zeroGradient;
}
stator
{
type zeroGradient;
}
front
{
type empty;
}
back
{
type empty;
}
}
// ************************************************************************* //

View File

@ -33,9 +33,50 @@ tracer0
}
}
injectionRateSuSp
sources
{
tracer0 (0.01 0);
tracer0
{
explicit constant 0.01;
}
}
}
source2
{
type scalarSemiImplicitSource;
volumeMode specific;
selectionMode all;
sources
{
tracer0
{
explicit
{
type exprField;
functions<scalar>
{
square
{
type square;
scale 0.0025;
level 0.0025;
frequency 10;
}
}
expression
#{
(hypot(pos().x() + 0.025, pos().y()) < 0.01)
? fn:square(time())
: 0
#};
}
}
}
}
}

View File

@ -49,10 +49,24 @@ IOdictionary fvOptions
IOobject::NO_WRITE
)
);
const dictionary& gradPDict =
fvOptions.subDict("momentumSource").subDict("injectionRateSuSp");
const scalar K =
gradPDict.get<Tuple2<vector, scalar>>("U").first().x();
scalar K(0);
// Get x() component from U source
{
const dictionary& momSource = fvOptions.subDict("momentumSource");
if (momSource.findDict("sources"))
{
K = momSource.subDict("sources")
.get<Tuple2<vector, scalar>>("U").first().x();
}
else
{
K = momSource.subDict("injectionRateSuSp")
.get<Tuple2<vector, scalar>>("U").first().x();
}
}
dictionary probes(IFstream(runTime.system()/"probes")());
const point location = pointField(probes.lookup("probeLocations"))[0];