36Foam::label Foam::turbulentDFSEMInletFvPatchVectorField::seedIterMax_ = 1000;
40void 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();
48 const pointField& localPoints = pp.localPoints();
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);
83void 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);
115void Foam::turbulentDFSEMInletFvPatchVectorField::initialisePatch()
123 const scalar error =
max(
magSqr(patchNormal_ + nf));
128 <<
"Patch " <<
patch().name() <<
" is not planar"
132 patchNormal_ /=
mag(patchNormal_) + ROOTVSMALL;
137 const polyPatch&
patch = this->patch().patch();
141 DynamicList<label> triToFace(2*
patch.size());
142 DynamicList<scalar> triMagSf(2*
patch.size());
144 DynamicList<face> tris(5);
147 triMagSf.append(0.0);
151 const face&
f =
patch[faceI];
158 triToFace.append(faceI);
169 for (label i = 1; i < triMagSf.size(); ++i)
171 triMagSf[i] += triMagSf[i-1];
175 triFace_.transfer(triFace);
176 triToFace_.transfer(triToFace);
177 triCumulativeMagSf_.transfer(triMagSf);
180 for (label i = 1; i < sumTriMagSf_.size(); ++i)
182 sumTriMagSf_[i] += sumTriMagSf_[i-1];
186 patchArea_ = sumTriMagSf_.last();
189 patchBounds_ = boundBox(
patch.localPoints(),
false);
190 patchBounds_.inflate(0.1);
194 reduce(singleProc_, orOp<bool>());
198void Foam::turbulentDFSEMInletFvPatchVectorField::initialiseEddyBox()
202 const scalarField L(L_->value(db().time().timeOutputValue())/Lref_);
210 scalar&
s = sigmax_[faceI];
215 s =
min(
mag(
L[faceI]), kappa_*delta_);
216 s =
max(
s, nCellPerEddy_*cellDx[faceI]);
220 maxSigmaX_ =
gMax(sigmax_);
223 v0_ = 2*
gSum(magSf)*maxSigmaX_;
227 Info<<
"Patch: " <<
patch().patch().name() <<
" eddy box:" <<
nl
228 <<
" volume : " << v0_ <<
nl
229 <<
" maxSigmaX : " << maxSigmaX_ <<
nl
243 const polyPatch&
patch = this->patch().patch();
248 const scalar areaFraction =
249 rndGen_.globalPosition<scalar>(0, patchArea_);
255 if (areaFraction >= sumTriMagSf_[i])
266 const scalar offset = sumTriMagSf_[procI];
269 if (areaFraction > triCumulativeMagSf_[i] + offset)
277 const face& tf = triFace_[triI];
281 pos.setIndex(triToFace_[triI]);
282 pos.rawPoint() = tri.randomPoint(rndGen_);
289 const scalar maxAreaLimit = triCumulativeMagSf_.last();
290 const scalar areaFraction = rndGen_.position<scalar>(0, maxAreaLimit);
294 if (areaFraction > triCumulativeMagSf_[i])
302 const face& tf = triFace_[triI];
306 pos.setIndex(triToFace_[triI]);
307 pos.rawPoint() = tri.randomPoint(rndGen_);
314void Foam::turbulentDFSEMInletFvPatchVectorField::initialiseEddies()
316 const scalar t = db().time().timeOutputValue();
319 DynamicList<eddy> eddies(size());
322 scalar sumVolEddy = 0;
323 scalar sumVolEddyAllProc = 0;
325 while (sumVolEddyAllProc/v0_ < d_)
330 while (
search && iter++ < seedIterMax_)
334 const label patchFaceI =
pos.index();
337 if (patchFaceI != -1)
343 rndGen_.position<scalar>(-maxSigmaX_, maxSigmaX_),
350 if (
e.patchFaceI() != -1)
353 sumVolEddy +=
e.volume();
363 sumVolEddyAllProc =
returnReduce(sumVolEddy, sumOp<scalar>());
365 eddies_.transfer(eddies);
367 nEddy_ = eddies_.size();
371 Pout<<
"Patch:" <<
patch().patch().name();
378 Pout<<
" seeded:" << nEddy_ <<
" eddies" <<
endl;
381 reduce(nEddy_, sumOp<label>());
385 Info<<
"Turbulent DFSEM patch: " <<
patch().name()
386 <<
" seeded " << nEddy_ <<
" eddies with total volume "
393 <<
"Patch: " <<
patch().patch().name()
394 <<
" on field " << internalField().name()
395 <<
": No eddies seeded - please check your set-up"
401void Foam::turbulentDFSEMInletFvPatchVectorField::convectEddies
407 const scalar t = db().time().timeOutputValue();
416 eddy&
e = eddies_[eddyI];
417 e.move(deltaT*(UBulk & patchNormal_));
419 const scalar position0 =
e.x();
422 if (position0 > maxSigmaX_)
427 while (
search && iter++ < seedIterMax_)
431 const label patchFaceI =
pos.index();
437 position0 - floor(position0/maxSigmaX_)*maxSigmaX_,
443 if (
e.patchFaceI() != -1)
453 reduce(nRecycled, sumOp<label>());
455 if (debug && nRecycled > 0)
457 Info<<
"Patch: " <<
patch().patch().name()
458 <<
" recycled " << nRecycled <<
" eddies"
464Foam::vector Foam::turbulentDFSEMInletFvPatchVectorField::uPrimeEddy
466 const List<eddy>& eddies,
467 const point& patchFaceCf
474 const eddy&
e = eddies[
k];
475 uPrime +=
e.uPrime(patchFaceCf, patchNormal_);
482void Foam::turbulentDFSEMInletFvPatchVectorField::calcOverlappingProcEddies
484 List<List<eddy>>& overlappingEddies
500 const eddy&
e = eddies_[i];
503 const point x(
e.position(patchNormal_));
504 boundBox ebb(
e.bounds());
513 if (ebb.overlaps(patchBBs[procI]))
515 dynSendMap[procI].append(i);
524 sendMap[procI].transfer(dynSendMap[procI]);
546 forAll(constructMap, procI)
552 constructMap[procI].setSize(nRecv);
554 for (label i = 0; i < nRecv; ++i)
556 constructMap[procI][i] = segmentI++;
561 mapDistribute map(segmentI, std::move(sendMap), std::move(constructMap));
567 const labelList& sendElems = map.subMap()[domain];
571 List<eddy> subEddies(UIndirectList<eddy>(eddies_, sendElems));
573 UOPstream toDomain(domain, pBufs);
575 toDomain<< subEddies;
580 pBufs.finishedSends();
585 const labelList& recvElems = map.constructMap()[domain];
589 UIPstream str(domain, pBufs);
591 str >> overlappingEddies[domain];
603Foam::turbulentDFSEMInletFvPatchVectorField::
604turbulentDFSEMInletFvPatchVectorField
626 triCumulativeMagSf_(),
629 patchBounds_(
boundBox::invertedBox),
634 sigmax_(size(),
Zero),
643Foam::turbulentDFSEMInletFvPatchVectorField::
644turbulentDFSEMInletFvPatchVectorField
653 U_(ptf.U_.clone(patch().patch())),
654 R_(ptf.R_.clone(patch().patch())),
655 L_(ptf.L_.clone(patch().patch())),
663 nCellPerEddy_(ptf.nCellPerEddy_),
665 patchArea_(ptf.patchArea_),
666 triFace_(ptf.triFace_),
667 triToFace_(ptf.triToFace_),
668 triCumulativeMagSf_(ptf.triCumulativeMagSf_),
669 sumTriMagSf_(ptf.sumTriMagSf_),
670 patchNormal_(ptf.patchNormal_),
671 patchBounds_(ptf.patchBounds_),
673 eddies_(ptf.eddies_),
675 rndGen_(ptf.rndGen_),
676 sigmax_(ptf.sigmax_, mapper),
677 maxSigmaX_(ptf.maxSigmaX_),
680 singleProc_(ptf.singleProc_),
681 writeEddies_(ptf.writeEddies_)
685Foam::turbulentDFSEMInletFvPatchVectorField::
686turbulentDFSEMInletFvPatchVectorField
704 nCellPerEddy_(
dict.getOrDefault<label>(
"nCellPerEddy", 5)),
709 triCumulativeMagSf_(),
712 patchBounds_(
boundBox::invertedBox),
717 sigmax_(size(),
Zero),
722 writeEddies_(
dict.getOrDefault(
"writeEddies", false))
726 const scalar t = db().time().timeOutputValue();
733Foam::turbulentDFSEMInletFvPatchVectorField::
734turbulentDFSEMInletFvPatchVectorField
740 U_(ptf.U_.clone(patch().patch())),
741 R_(ptf.R_.clone(patch().patch())),
742 L_(ptf.L_.clone(patch().patch())),
750 nCellPerEddy_(ptf.nCellPerEddy_),
752 patchArea_(ptf.patchArea_),
753 triFace_(ptf.triFace_),
754 triToFace_(ptf.triToFace_),
755 triCumulativeMagSf_(ptf.triCumulativeMagSf_),
756 sumTriMagSf_(ptf.sumTriMagSf_),
757 patchNormal_(ptf.patchNormal_),
758 patchBounds_(ptf.patchBounds_),
760 eddies_(ptf.eddies_),
762 rndGen_(ptf.rndGen_),
763 sigmax_(ptf.sigmax_),
764 maxSigmaX_(ptf.maxSigmaX_),
767 singleProc_(ptf.singleProc_),
768 writeEddies_(ptf.writeEddies_)
772Foam::turbulentDFSEMInletFvPatchVectorField::
773turbulentDFSEMInletFvPatchVectorField
780 U_(ptf.U_.clone(patch().patch())),
781 R_(ptf.R_.clone(patch().patch())),
782 L_(ptf.L_.clone(patch().patch())),
790 nCellPerEddy_(ptf.nCellPerEddy_),
792 patchArea_(ptf.patchArea_),
793 triFace_(ptf.triFace_),
794 triToFace_(ptf.triToFace_),
795 triCumulativeMagSf_(ptf.triCumulativeMagSf_),
796 sumTriMagSf_(ptf.sumTriMagSf_),
797 patchNormal_(ptf.patchNormal_),
798 patchBounds_(ptf.patchBounds_),
800 eddies_(ptf.eddies_),
802 rndGen_(ptf.rndGen_),
803 sigmax_(ptf.sigmax_),
804 maxSigmaX_(ptf.maxSigmaX_),
807 singleProc_(ptf.singleProc_),
808 writeEddies_(ptf.writeEddies_)
819 constexpr label maxDiffs = 5;
827 if (maxDiffs < nDiffs)
829 Info<<
"More than " << maxDiffs <<
" times"
830 <<
" Reynolds-stress realizability checks failed."
831 <<
" Skipping further comparisons." <<
endl;
840 <<
"Reynolds stress " << r <<
" at index " << i
841 <<
" does not obey the constraint: Rxx >= 0"
849 <<
"Reynolds stress " << r <<
" at index " << i
850 <<
" does not obey the constraint: Rxx*Ryy - sqr(Rxy) >= 0"
858 <<
"Reynolds stress " << r <<
" at index " << i
859 <<
" does not obey the constraint: det(R) >= 0"
880 <<
"Reynolds stresses contain at least one "
881 <<
"nonpositive element. min(R) = " <<
min(
R)
919 const auto& dfsemptf =
920 refCast<const turbulentDFSEMInletFvPatchVectorField>(ptf);
924 U_->rmap(dfsemptf.U_(), addr);
928 R_->rmap(dfsemptf.R_(), addr);
932 L_->rmap(dfsemptf.L_(), addr);
935 sigmax_.rmap(dfsemptf.sigmax_, addr);
946 if (curTimeIndex_ == -1)
956 if (curTimeIndex_ != db().time().
timeIndex())
959 U_->value(db().time().timeOutputValue())/Uref_;
965 /(
gSum(patch().magSf()) + ROOTVSMALL)
969 const scalar deltaT = db().time().deltaTValue();
970 convectEddies(UBulk, deltaT);
989 U[faceI] += c*uPrimeEddy(eddies_, Cf[faceI]);
997 U[faceI] += c*uPrimeEddy(eddies_, Cf[faceI]);
1002 calcOverlappingProcEddies(overlappingEddies);
1004 forAll(overlappingEddies, procI)
1006 const List<eddy>& eddies = overlappingEddies[procI];
1012 U[faceI] += c*uPrimeEddy(eddies, Cf[faceI]);
1019 const scalar fCorr =
1020 gSum((UBulk & patchNormal_)*patch().magSf())
1021 /
gSum(
U & -patch().Sf());
1025 curTimeIndex_ = db().time().timeIndex();
1034 Info<<
"Magnitude of bulk velocity: " << UBulk <<
endl;
1036 label
n = eddies_.size();
1040 Info<<
"Patch:" << patch().patch().name()
1041 <<
" min/max(U):" <<
gMin(
U) <<
", " <<
gMax(
U)
1044 if (db().time().writeTime())
1046 writeLumleyCoeffs();
1051 fixedValueFvPatchVectorField::updateCoeffs();
1079 writeEntry(
"value",
os);
#define R(A, B, C, D, E, F, K, M)
Macros for easy insertion into run-time selection tables.
volVectorField UMean(UMeanHeader, mesh)
Field with dimensions and associated with geometry type GeoMesh which is used to size the field and a...
void append(const T &val)
Append an element at the end of the list.
An Ostream is an abstract base class for all output systems (streams, files, token lists,...
Ostream & writeEntry(const keyType &key, const T &value)
Write a keyword/value entry.
Ostream & writeEntryIfDifferent(const word &key, const T &value1, const T &value2)
Write a keyword/value entry only when the two values differ.
Top level data entry class for use in dictionaries. Provides a mechanism to specify a variable as a c...
This class describes the interaction of (usually) a face and a point. It carries the info of a succes...
UPstream::rangeType allProcs() const noexcept
Range of ranks indices associated with PstreamBuffers.
label nProcs() const noexcept
Number of ranks associated with PstreamBuffers.
Inter-processor communications stream.
static void allGatherList(const List< commsStruct > &comms, List< T > &values, const int tag, const label comm)
static void listCombineAllGather(const List< commsStruct > &comms, List< T > &values, const CombineOp &cop, const int tag, const label comm)
After completion all processors have the same data.
void size(const label n)
Older name for setAddressableSize.
@ nonBlocking
"nonBlocking"
static int & msgType() noexcept
Message tag of standard messages.
static bool & parRun() noexcept
Test if this a parallel run.
void autoMap(const fvPatchFieldMapper &)
Map (and resize as needed) from self given a mapping object.
void rmap(const atmBoundaryLayer &, const labelList &)
Reverse map the given fvPatchField onto this fvPatchField.
A bounding box defined in terms of min/max extrema points.
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
static int debug
Flag to activate debug statements.
This boundary condition supplies a fixed value constraint, and is the base class for a number of othe...
virtual bool write()
Write the output fields.
A FieldMapper for finite-volume patch fields.
A finiteVolume patch using a polyPatch and a fvBoundaryMesh.
static const complex max
complex (VGREAT,VGREAT)
int myProcNo() const noexcept
Return processor number.
A class for managing temporary objects.
The turbulentDFSEMInlet is a synthesised-eddy based velocity inlet boundary condition to generate syn...
static void checkStresses(const symmTensorField &R)
Check if input Reynolds stresses are valid.
virtual void rmap(const fvPatchVectorField &ptf, const labelList &addr)
Reverse map the given fvPatchField onto this fvPatchField.
virtual void autoMap(const fvPatchFieldMapper &m)
Map (and resize as needed) from self given a mapping object.
virtual void updateCoeffs()
Update the coefficients associated with the patch field.
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
fileName path(UMean.rootPath()/UMean.caseName()/"graphs"/UMean.instance())
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
OBJstream os(runTime.globalPath()/outputName)
#define makePatchTypeField(PatchTypeField, typePatchTypeField)
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))
#define WarningInFunction
Report a warning using Foam::Warning.
const std::string patch
OpenFOAM patch number as a std::string.
dimensionedSymmTensor dev(const dimensionedSymmTensor &dt)
dimensionedScalar det(const dimensionedSphericalTensor &dt)
dimensionedScalar pos(const dimensionedScalar &ds)
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Type gSum(const FieldField< Field, Type > &f)
labelList identity(const label len, label start=0)
Return an identity map of the given length with (map[i] == i)
dimensionedScalar tr(const dimensionedSphericalTensor &dt)
List< label > labelList
A List of labels.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
PointIndexHit< point > pointIndexHit
A PointIndexHit for 3D points.
Cmpt invariantII(const SymmTensor< Cmpt > &st)
Return the 2nd invariant of a SymmTensor.
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
SymmTensor< scalar > symmTensor
SymmTensor of scalars, i.e. SymmTensor<scalar>.
messageStream Info
Information stream (stdout output on master, null elsewhere)
vectorField pointField
pointField is a vectorField.
vector point
Point is a vector.
Cmpt invariantIII(const SymmTensor< Cmpt > &st)
Return the 3rd invariant of a SymmTensor.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
Ostream & endl(Ostream &os)
Add newline and flush stream.
dimensionedScalar sqrt(const dimensionedScalar &ds)
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
List< labelList > labelListList
A List of labelList.
void reduce(const List< UPstream::commsStruct > &comms, T &value, const BinaryOp &bop, const int tag, const label comm)
triangle< point, const point & > triPointRef
A triangle using referred points.
scalar diff(const triad &A, const triad &B)
Return a quantity of the difference between two triads.
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Field< vector > vectorField
Specialisation of Field<T> for vector.
Field< symmTensor > symmTensorField
Specialisation of Field<T> for symmTensor.
Type gAverage(const FieldField< Field, Type > &f)
static constexpr const zero Zero
Global zero (0)
tmp< DimensionedField< TypeR, GeoMesh > > New(const tmp< DimensionedField< TypeR, GeoMesh > > &tdf1, const word &name, const dimensionSet &dimensions)
Global function forwards to reuseTmpDimensionedField::New.
dimensionedScalar cbrt(const dimensionedScalar &ds)
Type gMin(const FieldField< Field, Type > &f)
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
T returnReduce(const T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Reduce (copy) and return value.
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Type gMax(const FieldField< Field, Type > &f)
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
fileName search(const word &file, const fileName &directory)
Recursively search the given directory for the file.
constexpr char nl
The newline '\n' character (0x0a)
#define forAll(list, i)
Loop across all elements in list.
#define forAllReverse(list, i)
Reverse loop across all elements in list.
const vector L(dict.get< vector >("L"))