34 #include "surfaceInterpolate.H" 52 directionalPressureGradientExplicitSource,
65 Foam::fv::directionalPressureGradientExplicitSource::pressureDropModelNames_
67 { pressureDropModel::pVolumetricFlowRateTable,
"volumetricFlowRateTable" },
68 { pressureDropModel::pConstant,
"constant" },
69 { pressureDropModel::pDarcyForchheimer,
"DarcyForchheimer" },
75 void Foam::fv::directionalPressureGradientExplicitSource::initialise()
80 label numFaces = fZone.
size();
91 const label meshFacei = fZone[i];
95 label facePatchId = -1;
103 if (isA<emptyPolyPatch>(
pp))
108 const auto* cpp = isA<coupledPolyPatch>(
pp);
110 if (cpp && !cpp->owner())
120 faceId_[numFaces] =
faceId;
121 facePatchId_[numFaces] = facePatchId;
130 facePatchId_.
resize(numFaces);
134 void Foam::fv::directionalPressureGradientExplicitSource::writeProps
140 if (mesh_.time().writeTime())
146 name_ +
"Properties",
147 mesh_.time().timeName(),
165 const word& sourceName,
166 const word& modelType,
172 model_(pressureDropModelNames_.
get(
"model", coeffs_)),
173 gradP0_(cells_.size(),
Zero),
174 dGradP_(cells_.size(),
Zero),
175 gradPporous_(cells_.size(),
Zero),
176 flowDir_(coeffs_.
get<
vector>(
"flowDir")),
183 faceZoneName_(coeffs_.
get<
word>(
"faceZone")),
184 zoneID_(mesh_.faceZones().findZoneID(faceZoneName_)),
187 relaxationFactor_(coeffs_.getOrDefault<scalar>(
"relaxationFactor", 0.3)),
188 cellFaceMap_(cells_.size(), -1)
192 flowDir_.normalise();
197 <<
"Source can only be applied to a single field. Current " 204 <<
type() <<
" " << this->
name() <<
": " 205 <<
" Unknown face zone name: " << faceZoneName_
212 flowRate_ = interpolationTable<scalar>(
coeffs_);
227 <<
"Did not find mode " << model_ <<
nl 228 <<
"Please set 'model' to one of " 229 << pressureDropModelNames_
241 if (propsFile.good())
243 Info<<
" Reading pressure gradient from file" <<
endl;
245 propsDict.readEntry(
"gradient", gradP0_);
248 Info<<
" Initial pressure gradient = " << gradP0_ <<
nl <<
endl;
269 case pDarcyForchheimer:
273 const auto& turbModel =
281 gradPporous_ = -flowDir_*(D_*
nu + I_*0.5*magUn)*magUn*length_;
285 const auto& turbModel =
296 - flowDir_*(D_*
mu + I_*0.5*
rho*magUn)*magUn*length_;
302 gradPporous_ = -flowDir_*pressureDrop_;
306 case pVolumetricFlowRateTable:
308 scalar volFlowRate = 0;
313 label faceI = faceId_[i];
314 if (facePatchId_[i] != -1)
316 label patchI = facePatchId_[i];
317 totalphi +=
phi.boundaryField()[patchI][faceI];
321 totalphi +=
phi[faceI];
324 reduce(totalphi, sumOp<scalar>());
328 volFlowRate =
mag(totalphi);
332 const auto& turbModel =
340 volFlowRate =
mag(totalphi)/rhoAve;
343 gradPporous_ = -flowDir_*flowRate_(volFlowRate);
348 const faceZone& fZone = mesh_.faceZones()[zoneID_];
350 labelList meshToLocal(mesh_.nCells(), -1);
353 meshToLocal[cells_[i]] = i;
356 labelList faceToCellIndex(fZone.size(), -1);
357 const labelList& mc = fZone.masterCells();
358 const labelList& sc = fZone.slaveCells();
362 label masterCellI = mc[i];
364 if (meshToLocal[masterCellI] != -1 && masterCellI != -1)
366 faceToCellIndex[i] = meshToLocal[masterCellI];
368 else if (meshToLocal[masterCellI] == -1)
371 <<
"Did not find cell " << masterCellI
372 <<
"in cellZone :" << zoneName()
381 const polyBoundaryMesh&
pbm = mesh_.boundaryMesh();
383 FieldField<Field, vector> upwindValues(
pbm.size());
385 forAll(
U.boundaryField(), patchI)
391 upwindValues.set(patchI, pf.patchNeighbourField());
393 else if (!isA<emptyFvPatchScalarField>(pf))
401 label faceI = fZone[i];
402 label
cellId = faceToCellIndex[i];
406 label sourceCellId = sc[i];
407 if (mesh_.isInternalFace(faceI))
409 scalar w = mesh_.magSf()[faceI];
410 UfCells[
cellId] +=
U[sourceCellId]*w;
411 UfCellWeights[
cellId] += w;
413 else if (fZone.flipMap()[i])
415 const label patchI =
pbm.patchID(faceI);
416 label localFaceI =
pbm[patchI].whichFace(faceI);
418 scalar w = mesh_.magSf().boundaryField()[patchI][localFaceI];
420 if (upwindValues.set(patchI))
422 UfCells[
cellId] += upwindValues[patchI][localFaceI]*w;
423 UfCellWeights[
cellId] += w;
429 UfCells /= UfCellWeights;
433 label cellI = cells_[i];
435 const vector Ufnorm(UfCells[i]/(
mag(UfCells[i]) + SMALL));
442 (
D & UfCells[i]) -
U[cellI]
448 Info<<
"Difference mag(U) = " 449 <<
mag(UfCells[i]) -
mag(
U[cellI])
451 Info<<
"Pressure drop in flowDir direction : " 452 << gradPporous_[i] <<
endl;
453 Info<<
"UfCell:= " << UfCells[i] <<
"U : " <<
U[cellI] <<
endl;
456 writeProps(gradP0_ + dGradP_);
462 fvMatrix<vector>& eqn,
466 DimensionedField<vector, volMesh>
Su 470 name_ + fieldNames_[fieldI] +
"Sup",
471 mesh_.time().timeName(),
493 this->addSup(eqn, fieldI);
512 mesh_.time().timeName(),
523 invAPtr_() = 1.0/eqn.
A();
550 coeffs.readEntry(
"flowDir", flowDir_);
551 flowDir_.normalise();
553 if (model_ == pConstant)
555 coeffs.readEntry(
"pressureDrop", pressureDrop_);
557 else if (model_ == pDarcyForchheimer)
559 coeffs.readEntry(
"D", D_);
560 coeffs.readEntry(
"I", I_);
561 coeffs.readEntry(
"length", length_);
List< ReturnType > get(const UPtrList< T > &list, const AccessOp &aop)
List of values generated by applying the access operation to each list item.
const polyBoundaryMesh & pbm
void size(const label n)
Older name for setAddressableSize.
fvPatchField< vector > fvPatchVectorField
errorManipArg< error, int > exit(error &err, const int errNo=1)
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
wordList fieldNames_
Field names to apply source to - populated by derived models.
void resize(const label len)
Adjust allocated size of list.
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...
defineTypeNameAndDebug(atmAmbientTurbSource, 0)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
fileName timePath() const
Return current time path.
const fvMesh & mesh_
Reference to the mesh database.
pressureDropModel
Modes of pressure drop.
constexpr char nl
The newline '\n' character (0x0a)
scalarProduct< Type, Type >::type gSumProd(const UList< Type > &f1, const UList< Type > &f2, const label comm)
dimensioned< vector > dimensionedVector
Dimensioned vector obtained from generic dimensioned type.
ThermalDiffusivity< CompressibleTurbulenceModel< fluidThermo > > turbulenceModel
Ostream & endl(Ostream &os)
Add newline and flush stream.
void resize_nocopy(const label len)
Adjust allocated size of list without necessarily.
Ignore writing from objectRegistry::writeObject()
Templated wrapper class to provide compressible turbulence models thermal diffusivity based thermal t...
const Time & time() const
Return the top-level database.
tensor rotationTensor(const vector &n1, const vector &n2)
Rotational transformation tensor from vector n1 to n2.
tmp< volScalarField > rAU
Macros for easy insertion into run-time selection tables.
bool isInternalFace(const label faceIndex) const noexcept
Return true if given face label is internal to the mesh.
IOdictionary propsDict(dictIO)
virtual void correct(volVectorField &U)
Correct the pressure gradient.
#define forAll(list, i)
Loop across all elements in list.
volVectorField gradP(fvc::grad(p))
bool readEntry(const word &keyword, T &val, enum keyType::option matchOpt=keyType::REGEX, IOobjectOption::readOption readOpt=IOobjectOption::MUST_READ) const
Find entry and assign to T val. FatalIOError if it is found and the number of tokens is incorrect...
Templated abstract base class for single-phase incompressible turbulence models.
fileName::Type type(const fileName &name, const bool followLink=true)
Return the file type: DIRECTORY or FILE, normally following symbolic links.
const dimensionSet dimVolume(pow3(dimLength))
Type gSum(const FieldField< Field, Type > &f)
const word name_
Source name.
static const word propertiesName
Default name of the turbulence properties dictionary.
const polyBoundaryMesh & boundaryMesh() const noexcept
Return boundary mesh.
A class for handling words, derived from Foam::string.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
label size() const noexcept
The number of entries in the list.
A special matrix type and solver, designed for finite volume solutions of scalar equations. Face addressing is used to make all matrix assembly and solution loops vectorise.
const word & name() const noexcept
Return const access to the source name.
label whichPatch(const label meshFacei) const
Return patch index for a given mesh face index. Uses binary search.
addToRunTimeSelectionTable(option, atmAmbientTurbSource, dictionary)
virtual void writeData(Ostream &os) const
Write the source properties.
An Ostream is an abstract base class for all output systems (streams, files, token lists...
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
int debug
Static debugging option.
OBJstream os(runTime.globalPath()/outputName)
const faceZoneMesh & faceZones() const noexcept
Return face zone mesh.
virtual bool read(const dictionary &dict)
Read source dictionary.
const dimensionedScalar mu
Atomic mass unit.
virtual void constrain(fvMatrix< vector > &eqn, const label fieldI)
Set 1/A coefficient.
tmp< volScalarField > A() const
Return the central coefficient.
Enum is a wrapper around a list of names/values that represent particular enumeration (or int) values...
directionalPressureGradientExplicitSource(const word &sourceName, const word &modelType, const dictionary &dict, const fvMesh &mesh)
Construct from explicit source name and mesh.
Mesh data needed to do the Finite Volume discretisation.
A List with indirect addressing. Like IndirectList but does not store addressing. ...
void resetApplied()
Resize/reset applied flag list for all fieldNames_ entries.
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0)
A special matrix type and solver, designed for finite volume solutions of scalar equations.
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.
wordList names() const
A list of the zone names.
messageStream Info
Information stream (stdout output on master, null elsewhere)
Field< vector > vectorField
Specialisation of Field<T> for vector.
Intermediate abstract class for handling cell-set options for the derived fvOptions.
const dimensionedScalar & D
T getOrDefault(const word &keyword, const T &deflt, enum keyType::option matchOpt=keyType::REGEX) const
Find and return a T, or return the given default value. FatalIOError if it is found and the number of...
List< label > labelList
A List of labels.
#define NotImplemented
Issue a FatalErrorIn for a function not currently implemented.
Defines the attributes of an object for which implicit objectRegistry management is supported...
dictionary coeffs_
Dictionary containing source coefficients.
uindirectPrimitivePatch pp(UIndirectList< face >(mesh.faces(), faceLabels), mesh.points())
static constexpr const zero Zero
Global zero (0)
virtual void addSup(fvMatrix< vector > &eqn, const label fieldI)
Add explicit contribution to momentum equation.