36 Foam::label Foam::turbulentDFSEMInletFvPatchVectorField::seedIterMax_ = 1000;
40 void Foam::turbulentDFSEMInletFvPatchVectorField::writeEddyOBJ()
const 44 OFstream
os(db().time().
path()/
"eddyBox.obj");
46 const polyPatch&
pp = this->
patch().patch();
47 const labelList& boundaryPoints =
pp.boundaryPoints();
50 const vector offset(patchNormal_*maxSigmaX_);
53 point p = localPoints[boundaryPoints[i]];
55 os <<
"v " <<
p.x() <<
" " <<
p.y() <<
" " <<
p.z() <<
nl;
60 point p = localPoints[boundaryPoints[i]];
62 os <<
"v " <<
p.x() <<
" " <<
p.y() <<
" " <<
p.z() <<
nl;
67 const Time& time = db().time();
70 time.path()/
"eddies_" +
Foam::name(time.timeIndex()) +
".obj" 73 label pointOffset = 0;
76 const eddy&
e = eddies_[eddyI];
77 pointOffset +=
e.writeSurfaceOBJ(pointOffset, patchNormal_,
os);
83 void Foam::turbulentDFSEMInletFvPatchVectorField::writeLumleyCoeffs()
const 87 OFstream
os(db().time().
path()/
"lumley_interpolated.out");
91 const scalar t = db().time().timeOutputValue();
107 const scalar xi =
cbrt(0.5*iii);
108 const scalar eta =
sqrt(-ii/3.0);
115 void Foam::turbulentDFSEMInletFvPatchVectorField::initialisePatch()
123 const scalar error =
max(
magSqr(patchNormal_ + nf));
128 <<
"Patch " <<
patch().name() <<
" is not planar" 132 patchNormal_ /=
mag(patchNormal_) + ROOTVSMALL;
135 const polyPatch&
patch = this->
patch().patch();
141 for (
const face&
f :
patch)
143 nTris +=
f.nTriangles();
146 DynamicList<labelledTri> dynTriFace(nTris);
147 DynamicList<face> tris(8);
151 const face&
f =
patch[facei];
156 for (
const auto& t : tris)
163 triFace_.transfer(dynTriFace);
172 triCumulativeMagSf_.resize_nocopy(triFace_.size()+1);
174 auto iter = triCumulativeMagSf_.begin();
177 scalar patchArea = 0;
181 for (
const auto& t : triFace_)
183 patchArea += t.mag(
points);
187 sumTriMagSf_.resize_nocopy(numProc+1);
192 slice[myProci] = patchArea;
197 for (label i = 1; i < sumTriMagSf_.size(); ++i)
199 sumTriMagSf_[i] += sumTriMagSf_[i-1];
204 patchArea_ = sumTriMagSf_.back();
207 patchBounds_ = boundBox(
patch.localPoints(),
false);
208 patchBounds_.inflate(0.1);
218 void Foam::turbulentDFSEMInletFvPatchVectorField::initialiseEddyBox()
222 const scalarField L(L_->value(db().time().timeOutputValue())/Lref_);
230 scalar&
s = sigmax_[faceI];
235 s =
min(
mag(
L[faceI]), kappa_*delta_);
236 s =
max(
s, nCellPerEddy_*cellDx[faceI]);
240 maxSigmaX_ =
gMax(sigmax_);
243 v0_ = 2*
gSum(magSf)*maxSigmaX_;
247 Info<<
"Patch: " <<
patch().patch().name() <<
" eddy box:" <<
nl 248 <<
" volume : " << v0_ <<
nl 249 <<
" maxSigmaX : " << maxSigmaX_ <<
nl 260 const polyPatch&
patch = this->
patch().patch();
267 const scalar areaFraction =
268 rndGen_.globalPosition<scalar>(0, patchArea_);
274 if (areaFraction >= sumTriMagSf_[i])
285 const scalar offset = sumTriMagSf_[proci];
288 if (areaFraction > triCumulativeMagSf_[i] + offset)
300 const scalar maxAreaLimit = triCumulativeMagSf_.back();
301 const scalar areaFraction = rndGen_.position<scalar>(0, maxAreaLimit);
305 if (areaFraction > triCumulativeMagSf_[i])
320 triFace_[triI].tri(
points).randomPoint(rndGen_),
321 triFace_[triI].index()
330 void Foam::turbulentDFSEMInletFvPatchVectorField::initialiseEddies()
332 const scalar t = db().time().timeOutputValue();
335 DynamicList<eddy> eddies(size());
338 scalar sumVolEddy = 0;
339 scalar sumVolEddyAllProc = 0;
341 while (sumVolEddyAllProc/v0_ < d_)
346 while (
search && iter++ < seedIterMax_)
350 const label patchFaceI =
pos.index();
353 if (patchFaceI != -1)
359 rndGen_.position<scalar>(-maxSigmaX_, maxSigmaX_),
366 if (
e.patchFaceI() != -1)
369 sumVolEddy +=
e.volume();
379 sumVolEddyAllProc =
returnReduce(sumVolEddy, sumOp<scalar>());
381 eddies_.transfer(eddies);
383 nEddy_ = eddies_.size();
387 Pout<<
"Patch:" <<
patch().patch().name();
394 Pout<<
" seeded:" << nEddy_ <<
" eddies" <<
endl;
397 reduce(nEddy_, sumOp<label>());
401 Info<<
"Turbulent DFSEM patch: " <<
patch().name()
402 <<
" seeded " << nEddy_ <<
" eddies with total volume " 409 <<
"Patch: " <<
patch().patch().name()
410 <<
" on field " << internalField().name()
411 <<
": No eddies seeded - please check your set-up" 417 void Foam::turbulentDFSEMInletFvPatchVectorField::convectEddies
423 const scalar t = db().time().timeOutputValue();
432 eddy&
e = eddies_[eddyI];
433 e.move(deltaT*(UBulk & patchNormal_));
435 const scalar position0 =
e.x();
438 if (position0 > maxSigmaX_)
443 while (
search && iter++ < seedIterMax_)
447 const label patchFaceI =
pos.index();
453 position0 - floor(position0/maxSigmaX_)*maxSigmaX_,
459 if (
e.patchFaceI() != -1)
471 reduce(nRecycled, sumOp<label>());
475 Info<<
"Patch: " <<
patch().patch().name()
476 <<
" recycled " << nRecycled <<
" eddies" 483 Foam::vector Foam::turbulentDFSEMInletFvPatchVectorField::uPrimeEddy
485 const List<eddy>& eddies,
486 const point& patchFaceCf
493 const eddy&
e = eddies[
k];
494 uPrime +=
e.uPrime(patchFaceCf, patchNormal_);
501 void Foam::turbulentDFSEMInletFvPatchVectorField::calcOverlappingProcEddies
503 List<List<eddy>>& overlappingEddies
518 const eddy&
e = eddies_[i];
521 const point x(
e.position(patchNormal_));
522 boundBox ebb(
e.bounds());
531 if (ebb.overlaps(patchBBs[proci]))
533 sendMap[proci].push_back(i);
540 PstreamBuffers pBufs;
546 UOPstream
os(proci, pBufs);
548 os << UIndirectList<eddy>(eddies_, sendMap[proci]);
552 pBufs.finishedSends();
554 for (
const int proci : pBufs.allProcs())
558 UIPstream is(proci, pBufs);
560 is >> overlappingEddies[proci];
592 triCumulativeMagSf_(),
600 sigmax_(size(),
Zero),
629 nCellPerEddy_(ptf.nCellPerEddy_),
631 patchArea_(ptf.patchArea_),
632 triFace_(ptf.triFace_),
633 triCumulativeMagSf_(ptf.triCumulativeMagSf_),
634 sumTriMagSf_(ptf.sumTriMagSf_),
635 patchNormal_(ptf.patchNormal_),
636 patchBounds_(ptf.patchBounds_),
638 eddies_(ptf.eddies_),
640 rndGen_(ptf.rndGen_),
641 sigmax_(ptf.sigmax_, mapper),
642 maxSigmaX_(ptf.maxSigmaX_),
645 singleProc_(ptf.singleProc_),
646 writeEddies_(ptf.writeEddies_)
669 nCellPerEddy_(
dict.getOrDefault<label>(
"nCellPerEddy", 5)),
673 triCumulativeMagSf_(),
681 sigmax_(size(),
Zero),
686 writeEddies_(
dict.getOrDefault(
"writeEddies", false))
690 const scalar t = db().time().timeOutputValue();
715 nCellPerEddy_(ptf.nCellPerEddy_),
717 patchArea_(ptf.patchArea_),
718 triFace_(ptf.triFace_),
719 triCumulativeMagSf_(ptf.triCumulativeMagSf_),
720 sumTriMagSf_(ptf.sumTriMagSf_),
721 patchNormal_(ptf.patchNormal_),
722 patchBounds_(ptf.patchBounds_),
724 eddies_(ptf.eddies_),
726 rndGen_(ptf.rndGen_),
727 sigmax_(ptf.sigmax_),
728 maxSigmaX_(ptf.maxSigmaX_),
731 singleProc_(ptf.singleProc_),
732 writeEddies_(ptf.writeEddies_)
753 constexpr label maxDiffs = 5;
761 if (maxDiffs < nDiffs)
763 Info<<
"More than " << maxDiffs <<
" times" 764 <<
" Reynolds-stress realizability checks failed." 765 <<
" Skipping further comparisons." <<
endl;
774 <<
"Reynolds stress " << r <<
" at index " << i
775 <<
" does not obey the constraint: Rxx >= 0" 780 if ((r.xx()*r.yy() -
sqr(r.xy())) < 0)
783 <<
"Reynolds stress " << r <<
" at index " << i
784 <<
" does not obey the constraint: Rxx*Ryy - sqr(Rxy) >= 0" 792 <<
"Reynolds stress " << r <<
" at index " << i
793 <<
" does not obey the constraint: det(R) >= 0" 814 <<
"Reynolds stresses contain at least one " 815 <<
"nonpositive element. min(R) = " <<
min(
R)
823 const fvPatchFieldMapper& m
853 const auto& dfsemptf =
854 refCast<const turbulentDFSEMInletFvPatchVectorField>(ptf);
858 U_->rmap(dfsemptf.U_(), addr);
862 R_->rmap(dfsemptf.R_(), addr);
866 L_->rmap(dfsemptf.L_(), addr);
869 sigmax_.rmap(dfsemptf.sigmax_, addr);
880 if (curTimeIndex_ == -1)
890 if (curTimeIndex_ != db().time().
timeIndex())
892 tmp<vectorField>
UMean =
893 U_->value(db().time().timeOutputValue())/Uref_;
903 const scalar deltaT = db().time().deltaTValue();
904 convectEddies(UBulk, deltaT);
923 U[faceI] +=
c*uPrimeEddy(eddies_, Cf[faceI]);
931 U[faceI] +=
c*uPrimeEddy(eddies_, Cf[faceI]);
936 calcOverlappingProcEddies(overlappingEddies);
938 for (
const List<eddy>& eddies : overlappingEddies)
944 U[faceI] +=
c*uPrimeEddy(eddies, Cf[faceI]);
952 gSum((UBulk & patchNormal_)*
patch().magSf())
957 curTimeIndex_ = db().time().timeIndex();
966 Info<<
"Magnitude of bulk velocity: " << UBulk <<
endl;
968 Info<<
"Number of eddies: " 972 Info<<
"Patch:" <<
patch().patch().name()
973 <<
" min/max(U):" <<
gMin(
U) <<
", " <<
gMax(
U)
976 if (db().time().writeTime())
983 fixedValueFvPatchVectorField::updateCoeffs();
1022 turbulentDFSEMInletFvPatchVectorField
Cmpt invariantII(const SymmTensor< Cmpt > &st)
Return the 2nd invariant of a SymmTensor.
scalar diff(const triad &A, const triad &B)
Return a quantity of the difference between two triads.
This boundary condition supplies a fixed value constraint, and is the base class for a number of othe...
Field< symmTensor > symmTensorField
Specialisation of Field<T> for symmTensor.
fvPatchField< vector > fvPatchVectorField
virtual void autoMap(const fvPatchFieldMapper &)
Map (and resize as needed) from self given a mapping object.
errorManipArg< error, int > exit(error &err, const int errNo=1)
static int incrMsgType(int val=1) noexcept
Increment the message tag for standard messages.
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...
const vector L(dict.get< vector >("L"))
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Type gMin(const FieldField< Field, Type > &f)
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
constexpr char nl
The newline '\n' character (0x0a)
dimensionedScalar sqrt(const dimensionedScalar &ds)
Ostream & endl(Ostream &os)
Add newline and flush stream.
A finiteVolume patch using a polyPatch and a fvBoundaryMesh.
PointIndexHit< point > pointIndexHit
A PointIndexHit with a 3D point.
static bool & parRun() noexcept
Test if this a parallel run.
tmp< DimensionedField< TypeR, GeoMesh > > New(const tmp< DimensionedField< TypeR, GeoMesh >> &tf1, const word &name, const dimensionSet &dimensions, const bool initCopy=false)
Global function forwards to reuseTmpDimensionedField::New.
Ostream & writeEntry(const keyType &key, const T &value)
Write a keyword/value entry.
volVectorField UMean(UMeanHeader, mesh)
static int & msgType() noexcept
Message tag of standard messages.
This class describes the interaction of an object (often a face) and a point. It carries the info of ...
void writeValueEntry(Ostream &os) const
Write *this field as a "value" entry.
label k
Boltzmann constant.
dimensionedScalar det(const dimensionedSphericalTensor &dt)
static int myProcNo(const label communicator=worldComm)
Rank of this process in the communicator (starting from masterNo()). Can be negative if the process i...
virtual void write(Ostream &) const
Write.
virtual void rmap(const fvPatchVectorField &ptf, const labelList &addr)
Reverse map the given fvPatchField onto this fvPatchField.
T returnReduce(const T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Perform reduction on a copy, using specified binary operation.
Macros for easy insertion into run-time selection tables.
virtual void updateCoeffs()
Update the coefficients associated with the patch field.
#define forAll(list, i)
Loop across all elements in list.
dimensionedScalar pos(const dimensionedScalar &ds)
static label nProcs(const label communicator=worldComm)
Number of ranks in parallel run (for given communicator). It is 1 for serial run. ...
vectorField pointField
pointField is a vectorField.
const dimensionedScalar e
Elementary charge.
dimensionedSymmTensor dev(const dimensionedSymmTensor &dt)
Type gSum(const FieldField< Field, Type > &f)
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for expressions::valueTypeCode::INVALID.
SymmTensor< scalar > symmTensor
SymmTensor of scalars, i.e. SymmTensor<scalar>.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Cmpt invariantIII(const SymmTensor< Cmpt > &st)
Return the 3rd invariant of a SymmTensor.
Inter-processor communications stream.
A FieldMapper for finite-volume patch fields.
dimensionedScalar cbrt(const dimensionedScalar &ds)
virtual void rmap(const fvPatchField< Type > &, const labelList &)
Reverse map the given fvPatchField onto this fvPatchField.
dimensionedScalar tr(const dimensionedSphericalTensor &dt)
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Ostream & writeEntryIfDifferent(const word &key, const T &value1, const T &value2)
Write a keyword/value entry only when the two values differ.
static int debug
Flag to activate debug statements.
Top level data entry class for use in dictionaries. Provides a mechanism to specify a variable as a c...
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
virtual void autoMap(const fvPatchFieldMapper &m)
Map (and resize as needed) from self given a mapping object.
int debug
Static debugging option.
Type gMax(const FieldField< Field, Type > &f)
The turbulentDFSEMInlet is a synthesised-eddy based velocity inlet boundary condition to generate syn...
OBJstream os(runTime.globalPath()/outputName)
fileName path(UMean.rootPath()/UMean.caseName()/"graphs"/UMean.instance())
turbulentDFSEMInletFvPatchVectorField(const fvPatch &, const DimensionedField< vector, volMesh > &)
Construct from patch and internal field.
T & emplace_back(Args &&... args)
Construct an element at the end of the list, return reference to the new list element.
static void allGatherList(UList< T > &values, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Gather data, but keep individual values separate. Uses MPI_Allgather or manual linear/tree communicat...
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
vector point
Point is a vector.
#define R(A, B, C, D, E, F, K, M)
#define WarningInFunction
Report a warning using Foam::Warning.
Type gAverage(const FieldField< Field, Type > &f)
const dimensionedScalar c
Speed of light in a vacuum.
Field with dimensions and associated with geometry type GeoMesh which is used to size the field and a...
const std::string patch
OpenFOAM patch number as a std::string.
messageStream Info
Information stream (stdout output on master, null elsewhere)
Field< vector > vectorField
Specialisation of Field<T> for vector.
fileName search(const word &file, const fileName &directory)
Recursively search the given directory for the file.
#define forAllReverse(list, i)
Reverse loop across all elements in list.
List< label > labelList
A List of labels.
static void checkStresses(const symmTensorField &R)
Check if input Reynolds stresses are valid.
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
bool returnReduceOr(const bool value, const label comm=UPstream::worldComm)
Perform logical (or) MPI Allreduce on a copy. Uses UPstream::reduceOr.
makePatchTypeField(fvPatchScalarField, atmBoundaryLayerInletEpsilonFvPatchScalarField)
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
void reduce(T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Reduce inplace (cf. MPI Allreduce) using linear/tree communication schedule.
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
uindirectPrimitivePatch pp(UIndirectList< face >(mesh.faces(), faceLabels), mesh.points())
virtual void write(Ostream &) const
Write.
static constexpr const zero Zero
Global zero (0)