39 Foam::label Foam::turbulentDFSEMInletFvPatchVectorField::seedIterMax_ = 1000;
43 void Foam::turbulentDFSEMInletFvPatchVectorField::writeEddyOBJ()
const
47 OFstream os(db().time().
path()/
"eddyBox.obj");
49 const polyPatch& pp = this->
patch().patch();
50 const labelList& boundaryPoints = pp.boundaryPoints();
51 const pointField& localPoints = pp.localPoints();
53 vector offset = patchNormal_*maxSigmaX_;
56 point p = localPoints[boundaryPoints[i]];
58 os <<
"v " <<
p.x() <<
" " <<
p.y() <<
" " <<
p.z() <<
nl;
63 point p = localPoints[boundaryPoints[i]];
65 os <<
"v " <<
p.x() <<
" " <<
p.y() <<
" " <<
p.z() <<
nl;
88 const Time& time = db().time();
91 time.path()/
"eddies_" +
Foam::name(time.timeIndex()) +
".obj"
94 label pointOffset = 0;
97 const eddy&
e = eddies_[eddyI];
98 pointOffset +=
e.writeSurfaceOBJ(pointOffset, patchNormal_, os);
104 void Foam::turbulentDFSEMInletFvPatchVectorField::writeLumleyCoeffs()
const
118 /db().time().caseConstant()
127 autoPtr<ISstream> isPtr
135 Field<symmTensor> Rexp(isPtr());
137 OFstream os(db().time().
path()/
"lumley_input.out");
154 scalar xi =
cbrt(0.5*iii);
155 scalar eta =
sqrt(-ii/3.0);
163 OFstream os(db().time().
path()/
"lumley_interpolated.out");
180 scalar xi =
cbrt(0.5*iii);
181 scalar eta =
sqrt(-ii/3.0);
190 Foam::turbulentDFSEMInletFvPatchVectorField::patchMapper()
const
193 if (mapperPtr_.empty())
196 fileName samplePointsFile
198 this->db().time().
path()
199 /this->db().time().caseConstant()
205 pointField samplePoints((IFstream(samplePointsFile)()));
210 <<
" Read " << samplePoints.size() <<
" sample points from "
211 << samplePointsFile <<
endl;
219 && mapMethod_ !=
"planarInterpolation"
225 new pointToPointPlanarInterpolation
239 void Foam::turbulentDFSEMInletFvPatchVectorField::initialisePatch()
247 scalar error =
max(
magSqr(patchNormal_ + nf));
252 <<
"Patch " <<
patch().name() <<
" is not planar"
256 patchNormal_ /=
mag(patchNormal_) + ROOTVSMALL;
261 const polyPatch&
patch = this->
patch().patch();
265 DynamicList<label> triToFace(2*
patch.size());
266 DynamicList<scalar> triMagSf(2*
patch.size());
268 DynamicList<face> tris(5);
271 triMagSf.append(0.0);
275 const face&
f =
patch[faceI];
282 triToFace.append(faceI);
290 sumTriMagSf_[i] = 0.0;
298 for (
label i = 1; i < triMagSf.size(); i++)
300 triMagSf[i] += triMagSf[i-1];
305 triToFace_.transfer(triToFace);
306 triCumulativeMagSf_.transfer(triMagSf);
309 for (
label i = 1; i < sumTriMagSf_.size(); i++)
311 sumTriMagSf_[i] += sumTriMagSf_[i-1];
315 patchArea_ = sumTriMagSf_.last();
318 patchBounds_ = boundBox(
patch.localPoints(),
false);
319 patchBounds_.inflate(0.1);
323 reduce(singleProc_, orOp<bool>());
327 void Foam::turbulentDFSEMInletFvPatchVectorField::initialiseEddyBox()
337 scalar&
s = sigmax_[faceI];
341 s =
min(
s, kappa_*delta_);
346 s =
max(
s, nCellPerEddy_*cellDx[faceI]);
350 maxSigmaX_ =
gMax(sigmax_);
353 v0_ = 2*
gSum(magSf)*maxSigmaX_;
357 Info<<
"Patch: " <<
patch().patch().name() <<
" eddy box:" <<
nl
358 <<
" volume : " << v0_ <<
nl
359 <<
" maxSigmaX : " << maxSigmaX_ <<
nl
373 const polyPatch&
patch = this->
patch().patch();
378 scalar areaFraction = rndGen_.globalPosition<scalar>(0, patchArea_);
384 if (areaFraction >= sumTriMagSf_[i])
395 scalar offset = sumTriMagSf_[procI];
398 if (areaFraction > triCumulativeMagSf_[i] + offset)
406 const face& tf = triFace_[triI];
410 pos.setIndex(triToFace_[triI]);
411 pos.rawPoint() = tri.randomPoint(rndGen_);
418 scalar maxAreaLimit = triCumulativeMagSf_.last();
419 scalar areaFraction = rndGen_.position<scalar>(0, maxAreaLimit);
423 if (areaFraction > triCumulativeMagSf_[i])
431 const face& tf = triFace_[triI];
435 pos.setIndex(triToFace_[triI]);
436 pos.rawPoint() = tri.randomPoint(rndGen_);
443 void Foam::turbulentDFSEMInletFvPatchVectorField::initialiseEddies()
445 DynamicList<eddy> eddies(size());
448 scalar sumVolEddy = 0;
449 scalar sumVolEddyAllProc = 0;
451 while (sumVolEddyAllProc/v0_ < d_)
456 while (
search && iter++ < seedIterMax_)
469 rndGen_.position<scalar>(-maxSigmaX_, maxSigmaX_),
476 if (
e.patchFaceI() != -1)
479 sumVolEddy +=
e.volume();
489 sumVolEddyAllProc =
returnReduce(sumVolEddy, sumOp<scalar>());
491 eddies_.transfer(eddies);
493 nEddy_ = eddies_.size();
497 Pout<<
"Patch:" <<
patch().patch().name();
504 Pout<<
" seeded:" << nEddy_ <<
" eddies" <<
endl;
507 reduce(nEddy_, sumOp<label>());
511 Info<<
"Turbulent DFSEM patch: " <<
patch().name()
512 <<
" seeded " << nEddy_ <<
" eddies with total volume "
519 <<
"Patch: " <<
patch().patch().name()
520 <<
" on field " << internalField().name()
521 <<
": No eddies seeded - please check your set-up" <<
endl;
526 void Foam::turbulentDFSEMInletFvPatchVectorField::convectEddies
537 eddy&
e = eddies_[eddyI];
538 e.move(deltaT*(UMean_ & patchNormal_));
540 const scalar position0 =
e.x();
543 if (position0 > maxSigmaX_)
548 while (
search && iter++ < seedIterMax_)
558 position0 - floor(position0/maxSigmaX_)*maxSigmaX_,
564 if (
e.patchFaceI() != -1)
574 reduce(nRecycled, sumOp<label>());
576 if (
debug && nRecycled > 0)
578 Info<<
"Patch: " <<
patch().patch().name() <<
" recycled "
579 << nRecycled <<
" eddies" <<
endl;
584 Foam::vector Foam::turbulentDFSEMInletFvPatchVectorField::uDashEddy
586 const List<eddy>& eddies,
587 const point& patchFaceCf
594 const eddy&
e = eddies[
k];
595 uDash +=
e.uDash(patchFaceCf, patchNormal_);
602 void Foam::turbulentDFSEMInletFvPatchVectorField::calcOverlappingProcEddies
604 List<List<eddy>>& overlappingEddies
621 const eddy&
e = eddies_[i];
624 point x =
e.position(patchNormal_);
625 boundBox ebb =
e.bounds();
634 if (ebb.overlaps(patchBBs[procI]))
636 dynSendMap[procI].append(i);
645 sendMap[procI].transfer(dynSendMap[procI]);
668 forAll(constructMap, procI)
674 constructMap[procI].setSize(nRecv);
676 for (
label i = 0; i < nRecv; i++)
678 constructMap[procI][i] = segmentI++;
683 mapDistribute map(segmentI, std::move(sendMap), std::move(constructMap));
689 const labelList& sendElems = map.subMap()[domain];
693 List<eddy> subEddies(UIndirectList<eddy>(eddies_, sendElems));
695 UOPstream toDomain(domain, pBufs);
697 toDomain<< subEddies;
702 pBufs.finishedSends();
707 const labelList& recvElems = map.constructMap()[domain];
711 UIPstream str(domain, pBufs);
713 str >> overlappingEddies[domain];
738 mapMethod_(
"nearestCell"),
740 interpolateR_(
false),
741 interpolateL_(
false),
742 interpolateU_(
false),
751 triCumulativeMagSf_(),
759 sigmax_(size(),
Zero),
783 perturb_(ptf.perturb_),
784 mapMethod_(ptf.mapMethod_),
786 interpolateR_(ptf.interpolateR_),
787 interpolateL_(ptf.interpolateL_),
788 interpolateU_(ptf.interpolateU_),
794 patchArea_(ptf.patchArea_),
795 triFace_(ptf.triFace_),
796 triToFace_(ptf.triToFace_),
797 triCumulativeMagSf_(ptf.triCumulativeMagSf_),
798 sumTriMagSf_(ptf.sumTriMagSf_),
800 eddies_(ptf.eddies_),
801 nCellPerEddy_(ptf.nCellPerEddy_),
802 patchNormal_(ptf.patchNormal_),
804 rndGen_(ptf.rndGen_),
805 sigmax_(ptf.sigmax_, mapper),
806 maxSigmaX_(ptf.maxSigmaX_),
809 patchBounds_(ptf.patchBounds_),
810 singleProc_(ptf.singleProc_),
811 writeEddies_(ptf.writeEddies_)
824 delta_(
dict.
get<scalar>(
"delta")),
834 R_(interpolateOrRead<symmTensor>(
"R",
dict, interpolateR_)),
835 L_(interpolateOrRead<scalar>(
"L",
dict, interpolateL_)),
836 U_(interpolateOrRead<vector>(
"U",
dict, interpolateU_)),
842 triCumulativeMagSf_(),
850 sigmax_(size(),
Zero),
878 perturb_(ptf.perturb_),
879 mapMethod_(ptf.mapMethod_),
881 interpolateR_(ptf.interpolateR_),
882 interpolateL_(ptf.interpolateL_),
883 interpolateU_(ptf.interpolateU_),
889 patchArea_(ptf.patchArea_),
890 triFace_(ptf.triFace_),
891 triToFace_(ptf.triToFace_),
892 triCumulativeMagSf_(ptf.triCumulativeMagSf_),
893 sumTriMagSf_(ptf.sumTriMagSf_),
895 eddies_(ptf.eddies_),
896 nCellPerEddy_(ptf.nCellPerEddy_),
897 patchNormal_(ptf.patchNormal_),
899 rndGen_(ptf.rndGen_),
900 sigmax_(ptf.sigmax_),
901 maxSigmaX_(ptf.maxSigmaX_),
904 patchBounds_(ptf.patchBounds_),
905 singleProc_(ptf.singleProc_),
906 writeEddies_(ptf.writeEddies_)
922 perturb_(ptf.perturb_),
923 mapMethod_(ptf.mapMethod_),
925 interpolateR_(ptf.interpolateR_),
926 interpolateL_(ptf.interpolateL_),
927 interpolateU_(ptf.interpolateU_),
933 patchArea_(ptf.patchArea_),
934 triFace_(ptf.triFace_),
935 triToFace_(ptf.triToFace_),
936 triCumulativeMagSf_(ptf.triCumulativeMagSf_),
937 sumTriMagSf_(ptf.sumTriMagSf_),
939 eddies_(ptf.eddies_),
940 nCellPerEddy_(ptf.nCellPerEddy_),
941 patchNormal_(ptf.patchNormal_),
943 rndGen_(ptf.rndGen_),
944 sigmax_(ptf.sigmax_),
945 maxSigmaX_(ptf.maxSigmaX_),
948 patchBounds_(ptf.patchBounds_),
949 singleProc_(ptf.singleProc_),
950 writeEddies_(ptf.writeEddies_)
971 <<
"Reynolds stress " <<
R <<
" at face " << facei
972 <<
" does not obey the constraint: R_xx > 0"
976 scalar a_xx =
sqrt(
R.xx());
978 scalar a_xy =
R.xy()/a_xx;
980 scalar a_yy_2 =
R.yy() -
sqr(a_xy);
985 <<
"Reynolds stress " <<
R <<
" at face " << facei
986 <<
" leads to an invalid Cholesky decomposition due to the "
987 <<
"constraint R_yy - sqr(a_xy) >= 0"
993 scalar a_xz =
R.xz()/a_xx;
995 scalar a_yz = (
R.yz() - a_xy*a_xz)*a_yy;
997 scalar a_zz_2 =
R.zz() -
sqr(a_xz) -
sqr(a_yz);
1002 <<
"Reynolds stress " <<
R <<
" at face " << facei
1003 <<
" leads to an invalid Cholesky decomposition due to the "
1004 <<
"constraint R_zz - sqr(a_xz) - sqr(a_yz) >= 0"
1013 <<
" a_xx:" << a_xx <<
" a_xy:" << a_xy <<
" a_xz:" << a_xy
1014 <<
" a_yy:" << a_yy <<
" a_yz:" << a_yz
1050 refCast<const turbulentDFSEMInletFvPatchVectorField>(ptf);
1052 R_.rmap(dfsemptf.R_, addr);
1053 L_.rmap(dfsemptf.L_, addr);
1054 U_.rmap(dfsemptf.U_, addr);
1059 sigmax_.rmap(dfsemptf.sigmax_, addr);
1070 if (curTimeIndex_ == -1)
1074 initialiseEddyBox();
1080 if (curTimeIndex_ != db().time().
timeIndex())
1084 label n = eddies_.size();
1089 const scalar deltaT = db().time().deltaTValue();
1092 convectEddies(deltaT);
1103 const scalar FACTOR = 2;
1113 U[faceI] +=
c*uDashEddy(eddies_, Cf[faceI]);
1121 U[faceI] +=
c*uDashEddy(eddies_, Cf[faceI]);
1126 calcOverlappingProcEddies(overlappingEddies);
1128 forAll(overlappingEddies, procI)
1130 const List<eddy>& eddies = overlappingEddies[procI];
1139 U[faceI] +=
c*uDashEddy(eddies, Cf[faceI]);
1147 gSum((UMean_ & patchNormal_)*
patch().magSf())
1154 Info<<
"Patch:" <<
patch().patch().name()
1158 curTimeIndex_ = db().time().timeIndex();
1165 if (
debug && db().time().writeTime())
1167 writeLumleyCoeffs();
1171 fixedValueFvPatchVectorField::updateCoeffs();
1178 writeEntry(
"value", os);
1201 if (!mapMethod_.empty())