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;
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];
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);
200 void Foam::turbulentDFSEMInletFvPatchVectorField::initialiseEddyBox()
204 const scalarField L(L_->value(db().time().timeOutputValue())/Lref_);
212 scalar&
s = sigmax_[faceI];
217 s =
min(
mag(
L[faceI]), kappa_*delta_);
218 s =
max(
s, nCellPerEddy_*cellDx[faceI]);
222 maxSigmaX_ =
gMax(sigmax_);
225 v0_ = 2*
gSum(magSf)*maxSigmaX_;
229 Info<<
"Patch: " <<
patch().patch().name() <<
" eddy box:" <<
nl 230 <<
" volume : " << v0_ <<
nl 231 <<
" maxSigmaX : " << maxSigmaX_ <<
nl 245 const polyPatch&
patch = this->
patch().patch();
250 const scalar areaFraction =
251 rndGen_.globalPosition<scalar>(0, patchArea_);
257 if (areaFraction >= sumTriMagSf_[i])
268 const scalar offset = sumTriMagSf_[procI];
271 if (areaFraction > triCumulativeMagSf_[i] + offset)
279 const face& tf = triFace_[triI];
282 pos.hitPoint(tri.randomPoint(rndGen_));
283 pos.setIndex(triToFace_[triI]);
290 const scalar maxAreaLimit = triCumulativeMagSf_.last();
291 const scalar areaFraction = rndGen_.position<scalar>(0, maxAreaLimit);
295 if (areaFraction > triCumulativeMagSf_[i])
303 const face& tf = triFace_[triI];
306 pos.hitPoint(tri.randomPoint(rndGen_));
307 pos.setIndex(triToFace_[triI]);
314 void 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" 401 void 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)
455 reduce(nRecycled, sumOp<label>());
459 Info<<
"Patch: " <<
patch().patch().name()
460 <<
" recycled " << nRecycled <<
" eddies" 467 Foam::vector Foam::turbulentDFSEMInletFvPatchVectorField::uPrimeEddy
469 const List<eddy>& eddies,
470 const point& patchFaceCf
477 const eddy&
e = eddies[
k];
478 uPrime +=
e.uPrime(patchFaceCf, patchNormal_);
485 void Foam::turbulentDFSEMInletFvPatchVectorField::calcOverlappingProcEddies
487 List<List<eddy>>& overlappingEddies
503 const eddy&
e = eddies_[i];
506 const point x(
e.position(patchNormal_));
507 boundBox ebb(
e.bounds());
516 if (ebb.overlaps(patchBBs[proci]))
518 sendMap[proci].push_back(i);
531 UOPstream
os(proci, pBufs);
533 os << UIndirectList<eddy>(eddies_, sendMap[proci]);
537 pBufs.finishedSends();
539 for (
const int proci : pBufs.allProcs())
543 UIPstream is(proci, pBufs);
545 is >> overlappingEddies[proci];
579 triCumulativeMagSf_(),
587 sigmax_(size(),
Zero),
616 nCellPerEddy_(ptf.nCellPerEddy_),
618 patchArea_(ptf.patchArea_),
619 triFace_(ptf.triFace_),
620 triToFace_(ptf.triToFace_),
621 triCumulativeMagSf_(ptf.triCumulativeMagSf_),
622 sumTriMagSf_(ptf.sumTriMagSf_),
623 patchNormal_(ptf.patchNormal_),
624 patchBounds_(ptf.patchBounds_),
626 eddies_(ptf.eddies_),
628 rndGen_(ptf.rndGen_),
629 sigmax_(ptf.sigmax_, mapper),
630 maxSigmaX_(ptf.maxSigmaX_),
633 singleProc_(ptf.singleProc_),
634 writeEddies_(ptf.writeEddies_)
657 nCellPerEddy_(
dict.getOrDefault<label>(
"nCellPerEddy", 5)),
662 triCumulativeMagSf_(),
670 sigmax_(size(),
Zero),
675 writeEddies_(
dict.getOrDefault(
"writeEddies", false))
679 const scalar t = db().time().timeOutputValue();
703 nCellPerEddy_(ptf.nCellPerEddy_),
705 patchArea_(ptf.patchArea_),
706 triFace_(ptf.triFace_),
707 triToFace_(ptf.triToFace_),
708 triCumulativeMagSf_(ptf.triCumulativeMagSf_),
709 sumTriMagSf_(ptf.sumTriMagSf_),
710 patchNormal_(ptf.patchNormal_),
711 patchBounds_(ptf.patchBounds_),
713 eddies_(ptf.eddies_),
715 rndGen_(ptf.rndGen_),
716 sigmax_(ptf.sigmax_),
717 maxSigmaX_(ptf.maxSigmaX_),
720 singleProc_(ptf.singleProc_),
721 writeEddies_(ptf.writeEddies_)
743 nCellPerEddy_(ptf.nCellPerEddy_),
745 patchArea_(ptf.patchArea_),
746 triFace_(ptf.triFace_),
747 triToFace_(ptf.triToFace_),
748 triCumulativeMagSf_(ptf.triCumulativeMagSf_),
749 sumTriMagSf_(ptf.sumTriMagSf_),
750 patchNormal_(ptf.patchNormal_),
751 patchBounds_(ptf.patchBounds_),
753 eddies_(ptf.eddies_),
755 rndGen_(ptf.rndGen_),
756 sigmax_(ptf.sigmax_),
757 maxSigmaX_(ptf.maxSigmaX_),
760 singleProc_(ptf.singleProc_),
761 writeEddies_(ptf.writeEddies_)
772 constexpr label maxDiffs = 5;
780 if (maxDiffs < nDiffs)
782 Info<<
"More than " << maxDiffs <<
" times" 783 <<
" Reynolds-stress realizability checks failed." 784 <<
" Skipping further comparisons." <<
endl;
793 <<
"Reynolds stress " << r <<
" at index " << i
794 <<
" does not obey the constraint: Rxx >= 0" 799 if ((r.xx()*r.yy() -
sqr(r.xy())) < 0)
802 <<
"Reynolds stress " << r <<
" at index " << i
803 <<
" does not obey the constraint: Rxx*Ryy - sqr(Rxy) >= 0" 811 <<
"Reynolds stress " << r <<
" at index " << i
812 <<
" does not obey the constraint: det(R) >= 0" 833 <<
"Reynolds stresses contain at least one " 834 <<
"nonpositive element. min(R) = " <<
min(
R)
842 const fvPatchFieldMapper& m
872 const auto& dfsemptf =
873 refCast<const turbulentDFSEMInletFvPatchVectorField>(ptf);
877 U_->rmap(dfsemptf.U_(), addr);
881 R_->rmap(dfsemptf.R_(), addr);
885 L_->rmap(dfsemptf.L_(), addr);
888 sigmax_.rmap(dfsemptf.sigmax_, addr);
899 if (curTimeIndex_ == -1)
909 if (curTimeIndex_ != db().time().
timeIndex())
911 tmp<vectorField>
UMean =
912 U_->value(db().time().timeOutputValue())/Uref_;
922 const scalar deltaT = db().time().deltaTValue();
923 convectEddies(UBulk, deltaT);
942 U[faceI] +=
c*uPrimeEddy(eddies_, Cf[faceI]);
950 U[faceI] +=
c*uPrimeEddy(eddies_, Cf[faceI]);
955 calcOverlappingProcEddies(overlappingEddies);
957 for (
const List<eddy>& eddies : overlappingEddies)
963 U[faceI] +=
c*uPrimeEddy(eddies, Cf[faceI]);
971 gSum((UBulk & patchNormal_)*
patch().magSf())
976 curTimeIndex_ = db().time().timeIndex();
985 Info<<
"Magnitude of bulk velocity: " << UBulk <<
endl;
987 Info<<
"Number of eddies: " 991 Info<<
"Patch:" <<
patch().patch().name()
992 <<
" min/max(U):" <<
gMin(
U) <<
", " <<
gMax(
U)
995 if (db().time().writeTime())
1002 fixedValueFvPatchVectorField::updateCoeffs();
1041 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...
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &f1)
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)
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)
void append(const T &val)
Append an element at the end of the list.
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.
static void allGatherList(List< T > &values, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Gather data, but keep individual values separate. Uses linear/tree communication. ...
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 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.
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.
triangle< point, const point & > triPointRef
A triangle using referred points.
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.
"nonBlocking" : (MPI_Isend, MPI_Irecv)
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.
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
uindirectPrimitivePatch pp(UIndirectList< face >(mesh.faces(), faceLabels), mesh.points())
static void listCombineReduce(List< T > &values, const CombineOp &cop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
After completion all processors have the same data.
virtual void write(Ostream &) const
Write.
static constexpr const zero Zero
Global zero (0)