34#include "surfaceInterpolate.H"
65Foam::fv::directionalPressureGradientExplicitSource::pressureDropModelNames_
67 { pressureDropModel::pVolumetricFlowRateTable,
"volumetricFlowRateTable" },
68 { pressureDropModel::pConstant,
"constant" },
69 { pressureDropModel::pDarcyForchheimer,
"DarcyForchheimer" },
75void Foam::fv::directionalPressureGradientExplicitSource::initialise()
80 facePatchId_.
setSize(fZone.size());
85 const label faceI = fZone[i];
88 label facePatchId = -1;
98 const auto* cpp = isA<coupledPolyPatch>(pp);
102 faceId = (cpp->owner() ? pp.whichFace(faceI) : -1);
104 else if (!isA<emptyPolyPatch>(pp))
106 faceId = pp.whichFace(faceI);
117 facePatchId_[
count] = facePatchId;
127void Foam::fv::directionalPressureGradientExplicitSource::writeProps
133 if (mesh_.time().writeTime())
139 name_ +
"Properties",
140 mesh_.time().timeName(),
155Foam::fv::directionalPressureGradientExplicitSource::
156directionalPressureGradientExplicitSource
158 const word& sourceName,
159 const word& modelType,
165 model_(pressureDropModelNames_.get(
"model", coeffs_)),
166 gradP0_(cells_.size(),
Zero),
167 dGradP_(cells_.size(),
Zero),
168 gradPporous_(cells_.size(),
Zero),
169 flowDir_(coeffs_.get<
vector>(
"flowDir")),
176 faceZoneName_(coeffs_.get<
word>(
"faceZone")),
177 zoneID_(mesh_.faceZones().findZoneID(faceZoneName_)),
180 relaxationFactor_(coeffs_.getOrDefault<scalar>(
"relaxationFactor", 0.3)),
181 cellFaceMap_(cells_.size(), -1)
185 flowDir_.normalise();
190 <<
"Source can only be applied to a single field. Current "
197 <<
type() <<
" " << this->
name() <<
": "
198 <<
" Unknown face zone name: " << faceZoneName_
220 <<
"Did not find mode " << model_ <<
nl
221 <<
"Please set 'model' to one of "
222 << pressureDropModelNames_
234 if (propsFile.
good())
236 Info<<
" Reading pressure gradient from file" <<
endl;
241 Info<<
" Initial pressure gradient = " << gradP0_ <<
nl <<
endl;
262 case pDarcyForchheimer:
266 const auto& turbModel =
274 gradPporous_ = -flowDir_*(D_*
nu + I_*0.5*magUn)*magUn*length_;
278 const auto& turbModel =
289 - flowDir_*(D_*
mu + I_*0.5*
rho*magUn)*magUn*length_;
295 gradPporous_ = -flowDir_*pressureDrop_;
299 case pVolumetricFlowRateTable:
301 scalar volFlowRate = 0;
306 label faceI = faceId_[i];
307 if (facePatchId_[i] != -1)
309 label patchI = facePatchId_[i];
310 totalphi +=
phi.boundaryField()[patchI][faceI];
314 totalphi +=
phi[faceI];
321 volFlowRate =
mag(totalphi);
325 const auto& turbModel =
333 volFlowRate =
mag(totalphi)/rhoAve;
336 gradPporous_ = -flowDir_*flowRate_(volFlowRate);
341 const faceZone& fZone = mesh_.faceZones()[zoneID_];
343 labelList meshToLocal(mesh_.nCells(), -1);
346 meshToLocal[cells_[i]] = i;
355 label masterCellI = mc[i];
357 if (meshToLocal[masterCellI] != -1 && masterCellI != -1)
359 faceToCellIndex[i] = meshToLocal[masterCellI];
361 else if (meshToLocal[masterCellI] == -1)
364 <<
"Did not find cell " << masterCellI
365 <<
"in cellZone :" << zoneName()
378 forAll(
U.boundaryField(), patchI)
386 else if (!isA<emptyFvPatchScalarField>(pf))
394 label faceI = fZone[i];
395 label
cellId = faceToCellIndex[i];
399 label sourceCellId = sc[i];
400 if (mesh_.isInternalFace(faceI))
402 scalar w = mesh_.magSf()[faceI];
403 UfCells[
cellId] +=
U[sourceCellId]*w;
404 UfCellWeights[
cellId] += w;
408 label patchI = pbm.
patchID()[faceI-mesh_.nInternalFaces()];
409 label localFaceI = pbm[patchI].whichFace(faceI);
411 scalar w = mesh_.magSf().boundaryField()[patchI][localFaceI];
413 if (upwindValues.
set(patchI))
415 UfCells[
cellId] += upwindValues[patchI][localFaceI]*w;
416 UfCellWeights[
cellId] += w;
422 UfCells /= UfCellWeights;
426 label cellI = cells_[i];
428 const vector Ufnorm(UfCells[i]/(
mag(UfCells[i]) + SMALL));
435 (
D & UfCells[i]) -
U[cellI]
441 Info<<
"Difference mag(U) = "
442 <<
mag(UfCells[i]) -
mag(
U[cellI])
444 Info<<
"Pressure drop in flowDir direction : "
445 << gradPporous_[i] <<
endl;
446 Info<<
"UfCell:= " << UfCells[i] <<
"U : " <<
U[cellI] <<
endl;
449 writeProps(gradP0_ + dGradP_);
463 name_ + fieldNames_[fieldI] +
"Sup",
464 mesh_.time().timeName(),
486 this->addSup(eqn, fieldI);
505 mesh_.time().timeName(),
516 invAPtr_() = 1.0/eqn.
A();
544 flowDir_.normalise();
546 if (model_ == pConstant)
548 coeffs.
readEntry(
"pressureDrop", pressureDrop_);
550 else if (model_ == pDarcyForchheimer)
Macros for easy insertion into run-time selection tables.
#define addToRunTimeSelectionTable(baseType, thisType, argNames)
Add to construction table with typeName as the key.
Field with dimensions and associated with geometry type GeoMesh which is used to size the field and a...
Enum is a wrapper around a list of names/values that represent particular enumeration (or int) values...
A field of fields is a PtrList of fields with reference counting.
Input from file stream, using an ISstream.
Defines the attributes of an object for which implicit objectRegistry management is supported,...
bool good() const noexcept
True if next operation might succeed.
Templated abstract base class for single-phase incompressible turbulence models.
void setSize(const label n)
Alias for resize()
An Ostream is an abstract base class for all output systems (streams, files, token lists,...
const T * set(const label i) const
virtual void correct()
Solve the turbulence equations and correct the turbulence viscosity.
virtual bool read()
Re-read model coefficients if they have changed.
Templated wrapper class to provide compressible turbulence models thermal diffusivity based thermal t...
fileName timePath() const
Return current time path.
A List with indirect addressing. Like IndirectList but does not store addressing.
void size(const label n)
Older name for setAddressableSize.
label size() const noexcept
The number of elements in the list.
wordList names() const
A list of the zone names.
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
T getOrDefault(const word &keyword, const T &deflt, enum keyType::option matchOpt=keyType::REGEX) const
bool readEntry(const word &keyword, T &val, enum keyType::option matchOpt=keyType::REGEX, bool mandatory=true) const
A subset of mesh faces organised as a primitive patch.
const labelList & masterCells() const
const boolList & flipMap() const noexcept
Return face flip map.
const labelList & slaveCells() const
Return labels of slave cells.
A special matrix type and solver, designed for finite volume solutions of scalar equations....
tmp< volScalarField > A() const
Return the central coefficient.
const dimensionSet & dimensions() const noexcept
Mesh data needed to do the Finite Volume discretisation.
const Time & time() const
Return the top-level database.
virtual tmp< Field< Type > > patchNeighbourField() const
Return patchField on the opposite patch of a coupled patch.
virtual bool coupled() const
Return true if this patch field is coupled.
Intermediate abstract class for handling cell-set options for the derived fvOptions.
Applies an explicit pressure gradient source in such a way to deflect the flow towards an specific di...
virtual void writeData(Ostream &os) const
Write the source properties.
virtual void addSup(fvMatrix< vector > &eqn, const label fieldI)
Add explicit contribution to momentum equation.
pressureDropModel
Modes of pressure drop.
@ pVolumetricFlowRateTable
virtual void constrain(fvMatrix< vector > &eqn, const label fieldI)
Set 1/A coefficient.
Base abstract class for handling finite volume options (i.e. fvOption).
const word & name() const noexcept
Return const access to the source name.
const fvMesh & mesh_
Reference to the mesh database.
wordList fieldNames_
Field names to apply source to - populated by derived models.
dictionary coeffs_
Dictionary containing source coefficients.
void resetApplied()
Resize/reset applied flag list for all fieldNames_ entries.
const word name_
Source name.
An interpolation/look-up table of scalar vs <Type> values. The reference scalar values must be monoto...
A polyBoundaryMesh is a polyPatch list with additional search methods and registered IO.
label whichPatch(const label faceIndex) const
Return patch index for a given face label.
const labelList & patchID() const
Per boundary face label the patch index.
const faceZoneMesh & faceZones() const noexcept
Return face zone mesh.
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
bool isInternalFace(const label faceIndex) const noexcept
Return true if given face label is internal to the mesh.
static const word propertiesName
Default name of the turbulence properties dictionary.
A class for handling words, derived from Foam::string.
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
const volScalarField & mu
volVectorField gradP(fvc::grad(p))
#define NotImplemented
Issue a FatalErrorIn for a function not currently implemented.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
OBJstream os(runTime.globalPath()/outputName)
A special matrix type and solver, designed for finite volume solutions of scalar equations.
tmp< volScalarField > rAU
unsigned int count(const UList< bool > &bools, const bool val=true)
Count number of 'true' entries.
scalarProduct< Type, Type >::type gSumProd(const UList< Type > &f1, const UList< Type > &f2, const label comm)
Type gSum(const FieldField< Field, Type > &f)
tensor rotationTensor(const vector &n1, const vector &n2)
Rotational transformation tensor from vector n1 to n2.
const dimensionSet dimArea(sqr(dimLength))
const dimensionSet dimVelocity
messageStream Info
Information stream (stdout output on master, null elsewhere)
fileName::Type type(const fileName &name, const bool followLink=true)
Return the file type: DIRECTORY or FILE, normally following symbolic links.
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)
Field< vector > vectorField
Specialisation of Field<T> for vector.
static constexpr const zero Zero
Global zero (0)
const dimensionSet dimVolume(pow3(dimLength))
errorManipArg< error, int > exit(error &err, const int errNo=1)
constexpr char nl
The newline '\n' character (0x0a)
const dimensionedScalar & D
#define forAll(list, i)
Loop across all elements in list.
IOdictionary propsDict(dictIO)