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
502 const eddy&
e = eddies_[i];
505 const point x(
e.position(patchNormal_));
506 boundBox ebb(
e.bounds());
515 if (ebb.overlaps(patchBBs[proci]))
517 sendMap[proci].push_back(i);
530 UOPstream
os(proci, pBufs);
532 os << UIndirectList<eddy>(eddies_, sendMap[proci]);
536 pBufs.finishedSends();
538 for (
const int proci : pBufs.allProcs())
542 UIPstream is(proci, pBufs);
544 is >> overlappingEddies[proci];
577 triCumulativeMagSf_(),
585 sigmax_(size(),
Zero),
614 nCellPerEddy_(ptf.nCellPerEddy_),
616 patchArea_(ptf.patchArea_),
617 triFace_(ptf.triFace_),
618 triToFace_(ptf.triToFace_),
619 triCumulativeMagSf_(ptf.triCumulativeMagSf_),
620 sumTriMagSf_(ptf.sumTriMagSf_),
621 patchNormal_(ptf.patchNormal_),
622 patchBounds_(ptf.patchBounds_),
624 eddies_(ptf.eddies_),
626 rndGen_(ptf.rndGen_),
627 sigmax_(ptf.sigmax_, mapper),
628 maxSigmaX_(ptf.maxSigmaX_),
631 singleProc_(ptf.singleProc_),
632 writeEddies_(ptf.writeEddies_)
655 nCellPerEddy_(
dict.getOrDefault<label>(
"nCellPerEddy", 5)),
660 triCumulativeMagSf_(),
668 sigmax_(size(),
Zero),
673 writeEddies_(
dict.getOrDefault(
"writeEddies", false))
677 const scalar t = db().time().timeOutputValue();
701 nCellPerEddy_(ptf.nCellPerEddy_),
703 patchArea_(ptf.patchArea_),
704 triFace_(ptf.triFace_),
705 triToFace_(ptf.triToFace_),
706 triCumulativeMagSf_(ptf.triCumulativeMagSf_),
707 sumTriMagSf_(ptf.sumTriMagSf_),
708 patchNormal_(ptf.patchNormal_),
709 patchBounds_(ptf.patchBounds_),
711 eddies_(ptf.eddies_),
713 rndGen_(ptf.rndGen_),
714 sigmax_(ptf.sigmax_),
715 maxSigmaX_(ptf.maxSigmaX_),
718 singleProc_(ptf.singleProc_),
719 writeEddies_(ptf.writeEddies_)
741 nCellPerEddy_(ptf.nCellPerEddy_),
743 patchArea_(ptf.patchArea_),
744 triFace_(ptf.triFace_),
745 triToFace_(ptf.triToFace_),
746 triCumulativeMagSf_(ptf.triCumulativeMagSf_),
747 sumTriMagSf_(ptf.sumTriMagSf_),
748 patchNormal_(ptf.patchNormal_),
749 patchBounds_(ptf.patchBounds_),
751 eddies_(ptf.eddies_),
753 rndGen_(ptf.rndGen_),
754 sigmax_(ptf.sigmax_),
755 maxSigmaX_(ptf.maxSigmaX_),
758 singleProc_(ptf.singleProc_),
759 writeEddies_(ptf.writeEddies_)
770 constexpr label maxDiffs = 5;
778 if (maxDiffs < nDiffs)
780 Info<<
"More than " << maxDiffs <<
" times" 781 <<
" Reynolds-stress realizability checks failed." 782 <<
" Skipping further comparisons." <<
endl;
791 <<
"Reynolds stress " << r <<
" at index " << i
792 <<
" does not obey the constraint: Rxx >= 0" 797 if ((r.xx()*r.yy() -
sqr(r.xy())) < 0)
800 <<
"Reynolds stress " << r <<
" at index " << i
801 <<
" does not obey the constraint: Rxx*Ryy - sqr(Rxy) >= 0" 809 <<
"Reynolds stress " << r <<
" at index " << i
810 <<
" does not obey the constraint: det(R) >= 0" 831 <<
"Reynolds stresses contain at least one " 832 <<
"nonpositive element. min(R) = " <<
min(
R)
840 const fvPatchFieldMapper& m
870 const auto& dfsemptf =
871 refCast<const turbulentDFSEMInletFvPatchVectorField>(ptf);
875 U_->rmap(dfsemptf.U_(), addr);
879 R_->rmap(dfsemptf.R_(), addr);
883 L_->rmap(dfsemptf.L_(), addr);
886 sigmax_.rmap(dfsemptf.sigmax_, addr);
897 if (curTimeIndex_ == -1)
907 if (curTimeIndex_ != db().time().
timeIndex())
909 tmp<vectorField>
UMean =
910 U_->value(db().time().timeOutputValue())/Uref_;
920 const scalar deltaT = db().time().deltaTValue();
921 convectEddies(UBulk, deltaT);
940 U[faceI] +=
c*uPrimeEddy(eddies_, Cf[faceI]);
948 U[faceI] +=
c*uPrimeEddy(eddies_, Cf[faceI]);
953 calcOverlappingProcEddies(overlappingEddies);
955 for (
const List<eddy>& eddies : overlappingEddies)
961 U[faceI] +=
c*uPrimeEddy(eddies, Cf[faceI]);
969 gSum((UBulk & patchNormal_)*
patch().magSf())
974 curTimeIndex_ = db().time().timeIndex();
983 Info<<
"Magnitude of bulk velocity: " << UBulk <<
endl;
985 Info<<
"Number of eddies: " 989 Info<<
"Patch:" <<
patch().patch().name()
990 <<
" min/max(U):" <<
gMin(
U) <<
", " <<
gMax(
U)
993 if (db().time().writeTime())
1000 fixedValueFvPatchVectorField::updateCoeffs();
1039 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)
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)
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 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.
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)