turbulentDFSEMInletFvPatchVectorField.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | www.openfoam.com
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8  Copyright (C) 2015 OpenFOAM Foundation
9  Copyright (C) 2016-2021 OpenCFD Ltd.
10 -------------------------------------------------------------------------------
11 License
12  This file is part of OpenFOAM.
13 
14  OpenFOAM is free software: you can redistribute it and/or modify it
15  under the terms of the GNU General Public License as published by
16  the Free Software Foundation, either version 3 of the License, or
17  (at your option) any later version.
18 
19  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22  for more details.
23 
24  You should have received a copy of the GNU General Public License
25  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26 
27 \*---------------------------------------------------------------------------*/
28 
31 #include "momentOfInertia.H"
32 #include "OFstream.H"
33 
34 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
35 
36 Foam::label Foam::turbulentDFSEMInletFvPatchVectorField::seedIterMax_ = 1000;
37 
38 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
39 
40 void Foam::turbulentDFSEMInletFvPatchVectorField::writeEddyOBJ() const
41 {
42  {
43  // Output the bounding box
44  OFstream os(db().time().path()/"eddyBox.obj");
45 
46  const polyPatch& pp = this->patch().patch();
47  const labelList& boundaryPoints = pp.boundaryPoints();
48  const pointField& localPoints = pp.localPoints();
49 
50  const vector offset(patchNormal_*maxSigmaX_);
51  forAll(boundaryPoints, i)
52  {
53  point p = localPoints[boundaryPoints[i]];
54  p += offset;
55  os << "v " << p.x() << " " << p.y() << " " << p.z() << nl;
56  }
57 
58  forAll(boundaryPoints, i)
59  {
60  point p = localPoints[boundaryPoints[i]];
61  p -= offset;
62  os << "v " << p.x() << " " << p.y() << " " << p.z() << nl;
63  }
64  }
65 
66  {
67  const Time& time = db().time();
68  OFstream os
69  (
70  time.path()/"eddies_" + Foam::name(time.timeIndex()) + ".obj"
71  );
72 
73  label pointOffset = 0;
74  forAll(eddies_, eddyI)
75  {
76  const eddy& e = eddies_[eddyI];
77  pointOffset += e.writeSurfaceOBJ(pointOffset, patchNormal_, os);
78  }
79  }
80 }
81 
82 
83 void Foam::turbulentDFSEMInletFvPatchVectorField::writeLumleyCoeffs() const
84 {
85  // Output list of xi vs eta
86 
87  OFstream os(db().time().path()/"lumley_interpolated.out");
88 
89  os << "# xi" << token::TAB << "eta" << endl;
90 
91  const scalar t = db().time().timeOutputValue();
92  const symmTensorField R(R_->value(t)/sqr(Uref_));
93 
94  forAll(R, faceI)
95  {
96  // Normalised anisotropy tensor
97  const symmTensor devR(dev(R[faceI]/(tr(R[faceI]))));
98 
99  // Second tensor invariant
100  const scalar ii = min(0, invariantII(devR));
101 
102  // Third tensor invariant
103  const scalar iii = invariantIII(devR);
104 
105  // xi, eta
106  // See Pope - characterization of Reynolds-stress anisotropy
107  const scalar xi = cbrt(0.5*iii);
108  const scalar eta = sqrt(-ii/3.0);
109  os << xi << token::TAB << eta << token::TAB
110  << ii << token::TAB << iii << endl;
111  }
112 }
113 
114 
115 void Foam::turbulentDFSEMInletFvPatchVectorField::initialisePatch()
116 {
117  const vectorField nf(patch().nf());
118 
119  // Patch normal points into domain
120  patchNormal_ = -gAverage(nf);
121 
122  // Check that patch is planar
123  const scalar error = max(magSqr(patchNormal_ + nf));
124 
125  if (error > SMALL)
126  {
128  << "Patch " << patch().name() << " is not planar"
129  << endl;
130  }
131 
132  patchNormal_ /= mag(patchNormal_) + ROOTVSMALL;
133 
134 
135  // Decompose the patch faces into triangles to enable point search
136 
137  const polyPatch& patch = this->patch().patch();
138  const pointField& points = patch.points();
139 
140  // Triangulate the patch faces and create addressing
141  DynamicList<label> triToFace(2*patch.size());
142  DynamicList<scalar> triMagSf(2*patch.size());
143  DynamicList<face> triFace(2*patch.size());
144  DynamicList<face> tris(5);
145 
146  // Set zero value at the start of the tri area list
147  triMagSf.append(0.0);
148 
149  forAll(patch, faceI)
150  {
151  const face& f = patch[faceI];
152 
153  tris.clear();
154  f.triangles(points, tris);
155 
156  forAll(tris, i)
157  {
158  triToFace.append(faceI);
159  triFace.append(tris[i]);
160  triMagSf.append(tris[i].mag(points));
161  }
162  }
163 
164  for (auto& s : sumTriMagSf_)
165  {
166  s = 0.0;
167  }
168 
169  sumTriMagSf_[Pstream::myProcNo() + 1] = sum(triMagSf);
170 
171  Pstream::listCombineGather(sumTriMagSf_, maxEqOp<scalar>());
172  Pstream::listCombineScatter(sumTriMagSf_);
173 
174  for (label i = 1; i < triMagSf.size(); ++i)
175  {
176  triMagSf[i] += triMagSf[i-1];
177  }
178 
179  // Transfer to persistent storage
180  triFace_.transfer(triFace);
181  triToFace_.transfer(triToFace);
182  triCumulativeMagSf_.transfer(triMagSf);
183 
184  // Convert sumTriMagSf_ into cumulative sum of areas per proc
185  for (label i = 1; i < sumTriMagSf_.size(); ++i)
186  {
187  sumTriMagSf_[i] += sumTriMagSf_[i-1];
188  }
189 
190  // Global patch area (over all processors)
191  patchArea_ = sumTriMagSf_.last();
192 
193  // Local patch bounds (this processor)
194  patchBounds_ = boundBox(patch.localPoints(), false);
195  patchBounds_.inflate(0.1);
196 
197  // Determine if all eddies spawned from a single processor
198  singleProc_ = patch.size() == returnReduce(patch.size(), sumOp<label>());
199  reduce(singleProc_, orOp<bool>());
200 }
201 
202 
203 void Foam::turbulentDFSEMInletFvPatchVectorField::initialiseEddyBox()
204 {
205  const scalarField& magSf = patch().magSf();
206 
207  const scalarField L(L_->value(db().time().timeOutputValue())/Lref_);
208 
209  // (PCF:Eq. 14)
210  const scalarField cellDx(max(Foam::sqrt(magSf), 2/patch().deltaCoeffs()));
211 
212  // Inialise eddy box extents
213  forAll(*this, faceI)
214  {
215  scalar& s = sigmax_[faceI];
216 
217  // Average length scale (SST:Eq. 24)
218  // Personal communication regarding (PCR:Eq. 14)
219  // - the min operator in Eq. 14 is a typo, and should be a max operator
220  s = min(mag(L[faceI]), kappa_*delta_);
221  s = max(s, nCellPerEddy_*cellDx[faceI]);
222  }
223 
224  // Maximum extent across all processors
225  maxSigmaX_ = gMax(sigmax_);
226 
227  // Eddy box volume
228  v0_ = 2*gSum(magSf)*maxSigmaX_;
229 
230  if (debug)
231  {
232  Info<< "Patch: " << patch().patch().name() << " eddy box:" << nl
233  << " volume : " << v0_ << nl
234  << " maxSigmaX : " << maxSigmaX_ << nl
235  << endl;
236  }
237 }
238 
239 
240 Foam::pointIndexHit Foam::turbulentDFSEMInletFvPatchVectorField::setNewPosition
241 (
242  const bool global
243 )
244 {
245  // Initialise to miss
246  pointIndexHit pos(false, vector::max, -1);
247 
248  const polyPatch& patch = this->patch().patch();
249  const pointField& points = patch.points();
250 
251  if (global)
252  {
253  const scalar areaFraction =
254  rndGen_.globalPosition<scalar>(0, patchArea_);
255 
256  // Determine which processor to use
257  label procI = 0;
258  forAllReverse(sumTriMagSf_, i)
259  {
260  if (areaFraction >= sumTriMagSf_[i])
261  {
262  procI = i;
263  break;
264  }
265  }
266 
267  if (Pstream::myProcNo() == procI)
268  {
269  // Find corresponding decomposed face triangle
270  label triI = 0;
271  const scalar offset = sumTriMagSf_[procI];
272  forAllReverse(triCumulativeMagSf_, i)
273  {
274  if (areaFraction > triCumulativeMagSf_[i] + offset)
275  {
276  triI = i;
277  break;
278  }
279  }
280 
281  // Find random point in triangle
282  const face& tf = triFace_[triI];
283  const triPointRef tri(points[tf[0]], points[tf[1]], points[tf[2]]);
284 
285  pos.setHit();
286  pos.setIndex(triToFace_[triI]);
287  pos.rawPoint() = tri.randomPoint(rndGen_);
288  }
289  }
290  else
291  {
292  // Find corresponding decomposed face triangle on local processor
293  label triI = 0;
294  const scalar maxAreaLimit = triCumulativeMagSf_.last();
295  const scalar areaFraction = rndGen_.position<scalar>(0, maxAreaLimit);
296 
297  forAllReverse(triCumulativeMagSf_, i)
298  {
299  if (areaFraction > triCumulativeMagSf_[i])
300  {
301  triI = i;
302  break;
303  }
304  }
305 
306  // Find random point in triangle
307  const face& tf = triFace_[triI];
308  const triPointRef tri(points[tf[0]], points[tf[1]], points[tf[2]]);
309 
310  pos.setHit();
311  pos.setIndex(triToFace_[triI]);
312  pos.rawPoint() = tri.randomPoint(rndGen_);
313  }
314 
315  return pos;
316 }
317 
318 
319 void Foam::turbulentDFSEMInletFvPatchVectorField::initialiseEddies()
320 {
321  const scalar t = db().time().timeOutputValue();
322  const symmTensorField R(R_->value(t)/sqr(Uref_));
323 
324  DynamicList<eddy> eddies(size());
325 
326  // Initialise eddy properties
327  scalar sumVolEddy = 0;
328  scalar sumVolEddyAllProc = 0;
329 
330  while (sumVolEddyAllProc/v0_ < d_)
331  {
332  bool search = true;
333  label iter = 0;
334 
335  while (search && iter++ < seedIterMax_)
336  {
337  // Get new parallel consistent position
338  pointIndexHit pos(setNewPosition(true));
339  const label patchFaceI = pos.index();
340 
341  // Note: only 1 processor will pick up this face
342  if (patchFaceI != -1)
343  {
344  eddy e
345  (
346  patchFaceI,
347  pos.hitPoint(),
348  rndGen_.position<scalar>(-maxSigmaX_, maxSigmaX_),
349  sigmax_[patchFaceI],
350  R[patchFaceI],
351  rndGen_
352  );
353 
354  // If eddy valid, patchFaceI is non-zero
355  if (e.patchFaceI() != -1)
356  {
357  eddies.append(e);
358  sumVolEddy += e.volume();
359  search = false;
360  }
361  }
362  // else eddy on remote processor
363 
364  reduce(search, andOp<bool>());
365  }
366 
367 
368  sumVolEddyAllProc = returnReduce(sumVolEddy, sumOp<scalar>());
369  }
370  eddies_.transfer(eddies);
371 
372  nEddy_ = eddies_.size();
373 
374  if (debug)
375  {
376  Pout<< "Patch:" << patch().patch().name();
377 
378  if (Pstream::parRun())
379  {
380  Pout<< " processor:" << Pstream::myProcNo();
381  }
382 
383  Pout<< " seeded:" << nEddy_ << " eddies" << endl;
384  }
385 
386  reduce(nEddy_, sumOp<label>());
387 
388  if (nEddy_ > 0)
389  {
390  Info<< "Turbulent DFSEM patch: " << patch().name()
391  << " seeded " << nEddy_ << " eddies with total volume "
392  << sumVolEddyAllProc
393  << endl;
394  }
395  else
396  {
398  << "Patch: " << patch().patch().name()
399  << " on field " << internalField().name()
400  << ": No eddies seeded - please check your set-up"
401  << endl;
402  }
403 }
404 
405 
406 void Foam::turbulentDFSEMInletFvPatchVectorField::convectEddies
407 (
408  const vector& UBulk,
409  const scalar deltaT
410 )
411 {
412  const scalar t = db().time().timeOutputValue();
413  const symmTensorField R(R_->value(t)/sqr(Uref_));
414 
415  // Note: all operations applied to local processor only
416 
417  label nRecycled = 0;
418 
419  forAll(eddies_, eddyI)
420  {
421  eddy& e = eddies_[eddyI];
422  e.move(deltaT*(UBulk & patchNormal_));
423 
424  const scalar position0 = e.x();
425 
426  // Check to see if eddy has exited downstream box plane
427  if (position0 > maxSigmaX_)
428  {
429  bool search = true;
430  label iter = 0;
431 
432  while (search && iter++ < seedIterMax_)
433  {
434  // Spawn new eddy with new random properties (intensity etc)
435  pointIndexHit pos(setNewPosition(false));
436  const label patchFaceI = pos.index();
437 
438  e = eddy
439  (
440  patchFaceI,
441  pos.hitPoint(),
442  position0 - floor(position0/maxSigmaX_)*maxSigmaX_,
443  sigmax_[patchFaceI],
444  R[patchFaceI],
445  rndGen_
446  );
447 
448  if (e.patchFaceI() != -1)
449  {
450  search = false;
451  }
452  }
453 
454  nRecycled++;
455  }
456  }
457 
458  reduce(nRecycled, sumOp<label>());
459 
460  if (debug && nRecycled > 0)
461  {
462  Info<< "Patch: " << patch().patch().name()
463  << " recycled " << nRecycled << " eddies"
464  << endl;
465  }
466 }
467 
468 
469 Foam::vector Foam::turbulentDFSEMInletFvPatchVectorField::uPrimeEddy
470 (
471  const List<eddy>& eddies,
472  const point& patchFaceCf
473 ) const
474 {
475  vector uPrime(Zero);
476 
477  forAll(eddies, k)
478  {
479  const eddy& e = eddies[k];
480  uPrime += e.uPrime(patchFaceCf, patchNormal_);
481  }
482 
483  return uPrime;
484 }
485 
486 
487 void Foam::turbulentDFSEMInletFvPatchVectorField::calcOverlappingProcEddies
488 (
489  List<List<eddy>>& overlappingEddies
490 ) const
491 {
492  int oldTag = UPstream::msgType();
493  UPstream::msgType() = oldTag + 1;
494 
495  List<boundBox> patchBBs(Pstream::nProcs());
496  patchBBs[Pstream::myProcNo()] = patchBounds_;
497  Pstream::gatherList(patchBBs);
498  Pstream::scatterList(patchBBs);
499 
500  // Per processor indices into all segments to send
501  List<DynamicList<label>> dynSendMap(Pstream::nProcs());
502 
503  forAll(eddies_, i)
504  {
505  // Collect overlapping eddies
506  const eddy& e = eddies_[i];
507 
508  // Eddy bounds
509  const point x(e.position(patchNormal_));
510  boundBox ebb(e.bounds());
511  ebb.min() += x;
512  ebb.max() += x;
513 
514  forAll(patchBBs, procI)
515  {
516  // Not including intersection with local patch
517  if (procI != Pstream::myProcNo())
518  {
519  if (ebb.overlaps(patchBBs[procI]))
520  {
521  dynSendMap[procI].append(i);
522  }
523  }
524  }
525  }
526 
527  labelListList sendMap(Pstream::nProcs());
528  forAll(sendMap, procI)
529  {
530  sendMap[procI].transfer(dynSendMap[procI]);
531  }
532 
533  // Send the number of eddies for local processors to receive
534  labelListList sendSizes(Pstream::nProcs());
535  sendSizes[Pstream::myProcNo()].setSize(Pstream::nProcs());
536  forAll(sendMap, procI)
537  {
538  sendSizes[Pstream::myProcNo()][procI] = sendMap[procI].size();
539  }
540  Pstream::gatherList(sendSizes);
541  Pstream::scatterList(sendSizes);
542 
543  // Determine order of receiving
544  labelListList constructMap(Pstream::nProcs());
545 
546  // Local segment first
547  constructMap[Pstream::myProcNo()] = identity
548  (
549  sendMap[Pstream::myProcNo()].size()
550  );
551 
552  label segmentI = constructMap[Pstream::myProcNo()].size();
553  forAll(constructMap, procI)
554  {
555  if (procI != Pstream::myProcNo())
556  {
557  // What I need to receive is what other processor is sending to me
558  const label nRecv = sendSizes[procI][Pstream::myProcNo()];
559  constructMap[procI].setSize(nRecv);
560 
561  for (label i = 0; i < nRecv; ++i)
562  {
563  constructMap[procI][i] = segmentI++;
564  }
565  }
566  }
567 
568  mapDistribute map(segmentI, std::move(sendMap), std::move(constructMap));
569 
570  PstreamBuffers pBufs(Pstream::commsTypes::nonBlocking);
571 
572  for (const int domain : Pstream::allProcs())
573  {
574  const labelList& sendElems = map.subMap()[domain];
575 
576  if (domain != Pstream::myProcNo() && sendElems.size())
577  {
578  List<eddy> subEddies(UIndirectList<eddy>(eddies_, sendElems));
579 
580  UOPstream toDomain(domain, pBufs);
581 
582  toDomain<< subEddies;
583  }
584  }
585 
586  // Start receiving
587  pBufs.finishedSends();
588 
589  // Consume
590  for (const int domain : Pstream::allProcs())
591  {
592  const labelList& recvElems = map.constructMap()[domain];
593 
594  if (domain != Pstream::myProcNo() && recvElems.size())
595  {
596  UIPstream str(domain, pBufs);
597  {
598  str >> overlappingEddies[domain];
599  }
600  }
601  }
602 
603  // Restore tag
604  UPstream::msgType() = oldTag;
605 }
606 
607 
608 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
609 
612 (
613  const fvPatch& p,
615 )
616 :
618  U_(nullptr),
619  R_(nullptr),
620  L_(nullptr),
621  delta_(1.0),
622  d_(1.0),
623  kappa_(0.41),
624  Uref_(1.0),
625  Lref_(1.0),
626  scale_(1.0),
627  m_(0.5),
628  nCellPerEddy_(5),
629 
630  patchArea_(-1),
631  triFace_(),
632  triToFace_(),
633  triCumulativeMagSf_(),
634  sumTriMagSf_(Pstream::nProcs() + 1, Zero),
635  patchNormal_(Zero),
636  patchBounds_(boundBox::invertedBox),
637 
638  eddies_(Zero),
639  v0_(Zero),
640  rndGen_(Pstream::myProcNo()),
641  sigmax_(size(), Zero),
642  maxSigmaX_(Zero),
643  nEddy_(0),
644  curTimeIndex_(-1),
645  singleProc_(false),
646  writeEddies_(false)
647 {}
648 
649 
652 (
654  const fvPatch& p,
656  const fvPatchFieldMapper& mapper
657 )
658 :
659  fixedValueFvPatchField<vector>(ptf, p, iF, mapper),
660  U_(ptf.U_.clone(patch().patch())),
661  R_(ptf.R_.clone(patch().patch())),
662  L_(ptf.L_.clone(patch().patch())),
663  delta_(ptf.delta_),
664  d_(ptf.d_),
665  kappa_(ptf.kappa_),
666  Uref_(ptf.Uref_),
667  Lref_(ptf.Lref_),
668  scale_(ptf.scale_),
669  m_(ptf.m_),
670  nCellPerEddy_(ptf.nCellPerEddy_),
671 
672  patchArea_(ptf.patchArea_),
673  triFace_(ptf.triFace_),
674  triToFace_(ptf.triToFace_),
675  triCumulativeMagSf_(ptf.triCumulativeMagSf_),
676  sumTriMagSf_(ptf.sumTriMagSf_),
677  patchNormal_(ptf.patchNormal_),
678  patchBounds_(ptf.patchBounds_),
679 
680  eddies_(ptf.eddies_),
681  v0_(ptf.v0_),
682  rndGen_(ptf.rndGen_),
683  sigmax_(ptf.sigmax_, mapper),
684  maxSigmaX_(ptf.maxSigmaX_),
685  nEddy_(ptf.nEddy_),
686  curTimeIndex_(-1),
687  singleProc_(ptf.singleProc_),
688  writeEddies_(ptf.writeEddies_)
689 {}
690 
691 
694 (
695  const fvPatch& p,
697  const dictionary& dict
698 )
699 :
704  delta_(dict.getCheck<scalar>("delta", scalarMinMax::ge(0))),
705  d_(dict.getCheckOrDefault<scalar>("d", 1, scalarMinMax::ge(SMALL))),
706  kappa_(dict.getCheckOrDefault<scalar>("kappa", 0.41, scalarMinMax::ge(0))),
707  Uref_(dict.getCheckOrDefault<scalar>("Uref", 1, scalarMinMax::ge(SMALL))),
708  Lref_(dict.getCheckOrDefault<scalar>("Lref", 1, scalarMinMax::ge(SMALL))),
709  scale_(dict.getCheckOrDefault<scalar>("scale", 1, scalarMinMax::ge(0))),
710  m_(dict.getCheckOrDefault<scalar>("m", 0.5, scalarMinMax::ge(0))),
711  nCellPerEddy_(dict.getOrDefault<label>("nCellPerEddy", 5)),
712 
713  patchArea_(-1),
714  triFace_(),
715  triToFace_(),
716  triCumulativeMagSf_(),
717  sumTriMagSf_(Pstream::nProcs() + 1, Zero),
718  patchNormal_(Zero),
719  patchBounds_(boundBox::invertedBox),
720 
721  eddies_(),
722  v0_(Zero),
723  rndGen_(),
724  sigmax_(size(), Zero),
725  maxSigmaX_(Zero),
726  nEddy_(0),
727  curTimeIndex_(-1),
728  singleProc_(false),
729  writeEddies_(dict.getOrDefault("writeEddies", false))
730 {
731  eddy::debug = debug;
732 
733  const scalar t = db().time().timeOutputValue();
734  const symmTensorField R(R_->value(t)/sqr(Uref_));
735 
736  checkStresses(R);
737 }
738 
739 
742 (
744 )
745 :
747  U_(ptf.U_.clone(patch().patch())),
748  R_(ptf.R_.clone(patch().patch())),
749  L_(ptf.L_.clone(patch().patch())),
750  delta_(ptf.delta_),
751  d_(ptf.d_),
752  kappa_(ptf.kappa_),
753  Uref_(ptf.Uref_),
754  Lref_(ptf.Lref_),
755  scale_(ptf.scale_),
756  m_(ptf.m_),
757  nCellPerEddy_(ptf.nCellPerEddy_),
758 
759  patchArea_(ptf.patchArea_),
760  triFace_(ptf.triFace_),
761  triToFace_(ptf.triToFace_),
762  triCumulativeMagSf_(ptf.triCumulativeMagSf_),
763  sumTriMagSf_(ptf.sumTriMagSf_),
764  patchNormal_(ptf.patchNormal_),
765  patchBounds_(ptf.patchBounds_),
766 
767  eddies_(ptf.eddies_),
768  v0_(ptf.v0_),
769  rndGen_(ptf.rndGen_),
770  sigmax_(ptf.sigmax_),
771  maxSigmaX_(ptf.maxSigmaX_),
772  nEddy_(ptf.nEddy_),
773  curTimeIndex_(-1),
774  singleProc_(ptf.singleProc_),
775  writeEddies_(ptf.writeEddies_)
776 {}
777 
778 
781 (
784 )
785 :
787  U_(ptf.U_.clone(patch().patch())),
788  R_(ptf.R_.clone(patch().patch())),
789  L_(ptf.L_.clone(patch().patch())),
790  delta_(ptf.delta_),
791  d_(ptf.d_),
792  kappa_(ptf.kappa_),
793  Uref_(ptf.Uref_),
794  Lref_(ptf.Lref_),
795  scale_(ptf.scale_),
796  m_(ptf.m_),
797  nCellPerEddy_(ptf.nCellPerEddy_),
798 
799  patchArea_(ptf.patchArea_),
800  triFace_(ptf.triFace_),
801  triToFace_(ptf.triToFace_),
802  triCumulativeMagSf_(ptf.triCumulativeMagSf_),
803  sumTriMagSf_(ptf.sumTriMagSf_),
804  patchNormal_(ptf.patchNormal_),
805  patchBounds_(ptf.patchBounds_),
806 
807  eddies_(ptf.eddies_),
808  v0_(ptf.v0_),
809  rndGen_(ptf.rndGen_),
810  sigmax_(ptf.sigmax_),
811  maxSigmaX_(ptf.maxSigmaX_),
812  nEddy_(ptf.nEddy_),
813  curTimeIndex_(-1),
814  singleProc_(ptf.singleProc_),
815  writeEddies_(ptf.writeEddies_)
816 {}
817 
818 
819 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
820 
822 (
823  const symmTensorField& Rf
824 )
825 {
826  // Perform checks of the stress tensor based on Cholesky decomposition
827  // constraints
828 
829  forAll(Rf, facei)
830  {
831  const symmTensor& R = Rf[facei];
832 
833  if (R.xx() <= 0)
834  {
836  << "Reynolds stress " << R << " at face " << facei
837  << " does not obey the constraint: R_xx > 0"
838  << exit(FatalError);
839  }
840 
841  const scalar a_xx = sqrt(R.xx());
842 
843  const scalar a_xy = R.xy()/a_xx;
844 
845  const scalar a_yy_2 = R.yy() - sqr(a_xy);
846 
847  if (a_yy_2 < 0)
848  {
850  << "Reynolds stress " << R << " at face " << facei
851  << " leads to an invalid Cholesky decomposition due to the "
852  << "constraint R_yy - sqr(a_xy) >= 0"
853  << exit(FatalError);
854  }
855 
856  const scalar a_yy = Foam::sqrt(a_yy_2);
857 
858  const scalar a_xz = R.xz()/a_xx;
859 
860  const scalar a_yz = (R.yz() - a_xy*a_xz)/a_yy;
861 
862  const scalar a_zz_2 = R.zz() - sqr(a_xz) - sqr(a_yz);
863 
864  if (a_zz_2 < 0)
865  {
867  << "Reynolds stress " << R << " at face " << facei
868  << " leads to an invalid Cholesky decomposition due to the "
869  << "constraint R_zz - sqr(a_xz) - sqr(a_yz) >= 0"
870  << exit(FatalError);
871  }
872 
873  const scalar a_zz = Foam::sqrt(a_zz_2);
874 
875  if (debug)
876  {
877  Pout<< "R: " << R
878  << " a_xx:" << a_xx << " a_xy:" << a_xy << " a_xz:" << a_xy
879  << " a_yy:" << a_yy << " a_yz:" << a_yz
880  << " a_zz:" << a_zz
881  << endl;
882  }
883  }
884 
885  return true;
886 }
887 
888 
890 (
891  const fvPatchFieldMapper& m
892 )
893 {
895 
896  if (U_)
897  {
898  U_->autoMap(m);
899  }
900  if (R_)
901  {
902  R_->autoMap(m);
903  }
904  if (L_)
905  {
906  L_->autoMap(m);
907  }
908 
909  sigmax_.autoMap(m);
910 }
911 
912 
914 (
915  const fvPatchVectorField& ptf,
916  const labelList& addr
917 )
918 {
920 
921  const auto& dfsemptf =
922  refCast<const turbulentDFSEMInletFvPatchVectorField>(ptf);
923 
924  if (U_)
925  {
926  U_->rmap(dfsemptf.U_(), addr);
927  }
928  if (R_)
929  {
930  R_->rmap(dfsemptf.R_(), addr);
931  }
932  if (L_)
933  {
934  L_->rmap(dfsemptf.L_(), addr);
935  }
936 
937  sigmax_.rmap(dfsemptf.sigmax_, addr);
938 }
939 
940 
942 {
943  if (updated())
944  {
945  return;
946  }
947 
948  if (curTimeIndex_ == -1)
949  {
950  initialisePatch();
951 
952  initialiseEddyBox();
953 
954  initialiseEddies();
955  }
956 
957 
958  if (curTimeIndex_ != db().time().timeIndex())
959  {
961  U_->value(db().time().timeOutputValue())/Uref_;
962 
963  // (PCR:p. 522)
964  const vector UBulk
965  (
966  gSum(UMean()*patch().magSf())
967  /(gSum(patch().magSf()) + ROOTVSMALL)
968  );
969 
970  // Move eddies using bulk velocity
971  const scalar deltaT = db().time().deltaTValue();
972  convectEddies(UBulk, deltaT);
973 
974  // Set mean velocity
975  vectorField& U = *this;
976  U = UMean;
977 
978  // Apply second part of normalisation coefficient
979  const scalar c =
980  scale_*Foam::pow(10*v0_, m_)/Foam::sqrt(scalar(nEddy_));
981 
982  // In parallel, need to collect all eddies that will interact with
983  // local faces
984 
985  const pointField& Cf = patch().Cf();
986 
987  if (singleProc_ || !Pstream::parRun())
988  {
989  forAll(U, faceI)
990  {
991  U[faceI] += c*uPrimeEddy(eddies_, Cf[faceI]);
992  }
993  }
994  else
995  {
996  // Process local eddy contributions
997  forAll(U, faceI)
998  {
999  U[faceI] += c*uPrimeEddy(eddies_, Cf[faceI]);
1000  }
1001 
1002  // Add contributions from overlapping eddies
1003  List<List<eddy>> overlappingEddies(Pstream::nProcs());
1004  calcOverlappingProcEddies(overlappingEddies);
1005 
1006  forAll(overlappingEddies, procI)
1007  {
1008  const List<eddy>& eddies = overlappingEddies[procI];
1009 
1010  if (eddies.size())
1011  {
1012  forAll(U, faceI)
1013  {
1014  U[faceI] += c*uPrimeEddy(eddies, Cf[faceI]);
1015  }
1016  }
1017  }
1018  }
1019 
1020  // Re-scale to ensure correct flow rate
1021  const scalar fCorr =
1022  gSum((UBulk & patchNormal_)*patch().magSf())
1023  /gSum(U & -patch().Sf());
1024 
1025  U *= fCorr;
1026 
1027  curTimeIndex_ = db().time().timeIndex();
1028 
1029  if (writeEddies_)
1030  {
1031  writeEddyOBJ();
1032  }
1033 
1034  if (debug)
1035  {
1036  Info<< "Magnitude of bulk velocity: " << UBulk << endl;
1037 
1038  label n = eddies_.size();
1039  Info<< "Number of eddies: " << returnReduce(n, sumOp<label>())
1040  << endl;
1041 
1042  Info<< "Patch:" << patch().patch().name()
1043  << " min/max(U):" << gMin(U) << ", " << gMax(U)
1044  << endl;
1045 
1046  if (db().time().writeTime())
1047  {
1048  writeLumleyCoeffs();
1049  }
1050  }
1051  }
1052 
1053  fixedValueFvPatchVectorField::updateCoeffs();
1054 }
1055 
1056 
1058 {
1060  os.writeEntry("delta", delta_);
1061  os.writeEntryIfDifferent<scalar>("d", 1.0, d_);
1062  os.writeEntryIfDifferent<scalar>("kappa", 0.41, kappa_);
1063  os.writeEntryIfDifferent<scalar>("Uref", 1.0, Uref_);
1064  os.writeEntryIfDifferent<scalar>("Lref", 1.0, Lref_);
1065  os.writeEntryIfDifferent<scalar>("scale", 1.0, scale_);
1066  os.writeEntryIfDifferent<scalar>("m", 0.5, m_);
1067  os.writeEntryIfDifferent<label>("nCellPerEddy", 5, nCellPerEddy_);
1068  os.writeEntryIfDifferent("writeEddies", false, writeEddies_);
1069  if (U_)
1070  {
1071  U_->writeData(os);
1072  }
1073  if (R_)
1074  {
1075  R_->writeData(os);
1076  }
1077  if (L_)
1078  {
1079  L_->writeData(os);
1080  }
1081  writeEntry("value", os);
1082 }
1083 
1084 
1085 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
1086 
1087 namespace Foam
1088 {
1090  (
1093  );
1094 }
1095 
1096 
1097 // ************************************************************************* //
Foam::expressions::patchExpr::debug
int debug
Static debugging option.
Foam::fvPatchField< vector >
Foam::labelList
List< label > labelList
A List of labels.
Definition: List.H:67
Foam::pointField
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:44
Foam::fvPatchField::write
virtual void write(Ostream &) const
Write.
Definition: fvPatchField.C:384
Foam::Ostream::writeEntryIfDifferent
Ostream & writeEntryIfDifferent(const word &key, const T &value1, const T &value2)
Write a keyword/value entry only when the two values differ.
Definition: Ostream.H:248
Foam::scalarField
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Definition: primitiveFieldsFwd.H:52
Foam::turbulentDFSEMInletFvPatchVectorField::checkStresses
static bool checkStresses(const symmTensorField &Rf)
Return true if input Reynold stresses are valid.
Definition: turbulentDFSEMInletFvPatchVectorField.C:822
L
const vector L(dict.get< vector >("L"))
Foam::SymmTensor< scalar >
p
volScalarField & p
Definition: createFieldRefs.H:8
Foam::eddy::debug
static int debug
Flag to activate debug statements.
Definition: eddy.H:147
Foam::returnReduce
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
Definition: PstreamReduceOps.H:94
s
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))
Definition: gmvOutputSpray.H:25
Foam::invariantIII
Cmpt invariantIII(const SymmTensor< Cmpt > &st)
Return the 3rd invariant of a SymmTensor.
Definition: SymmTensorI.H:362
Foam::tmp
A class for managing temporary objects.
Definition: PtrList.H:61
Foam::Zero
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
momentOfInertia.H
Foam::gAverage
Type gAverage(const FieldField< Field, Type > &f)
Definition: FieldFieldFunctions.C:604
Foam::Pstream::scatterList
static void scatterList(const List< commsStruct > &comms, List< T > &Values, const int tag, const label comm)
Scatter data. Reverse of gatherList.
Definition: gatherScatterList.C:215
Foam::boundBox::invertedBox
static const boundBox invertedBox
A large inverted boundBox: min/max == +/- ROOTVGREAT.
Definition: boundBox.H:86
Foam::turbulentDFSEMInletFvPatchVectorField::write
virtual void write(Ostream &) const
Write.
Definition: turbulentDFSEMInletFvPatchVectorField.C:1057
Foam::List::append
void append(const T &val)
Append an element at the end of the list.
Definition: ListI.H:175
Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:369
Foam::gSum
Type gSum(const FieldField< Field, Type > &f)
Definition: FieldFieldFunctions.C:594
Foam::Pout
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
Foam::turbulentDFSEMInletFvPatchVectorField::turbulentDFSEMInletFvPatchVectorField
turbulentDFSEMInletFvPatchVectorField(const fvPatch &, const DimensionedField< vector, volMesh > &)
Construct from patch and internal field.
Definition: turbulentDFSEMInletFvPatchVectorField.C:612
Foam::min
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33
Foam::sumOp
Definition: ops.H:213
Foam::fixedValueFvPatchField
This boundary condition supplies a fixed value constraint, and is the base class for a number of othe...
Definition: fixedValueFvPatchField.H:80
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:296
OFstream.H
Foam::symmTensorField
Field< symmTensor > symmTensorField
Specialisation of Field<T> for symmTensor.
Definition: primitiveFieldsFwd.H:56
Foam::vectorField
Field< vector > vectorField
Specialisation of Field<T> for vector.
Definition: primitiveFieldsFwd.H:54
Foam::magSqr
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
n
label n
Definition: TABSMDCalcMethod2.H:31
Foam::reduce
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
Definition: PstreamReduceOps.H:51
R
#define R(A, B, C, D, E, F, K, M)
Foam::PointIndexHit
This class describes the interaction of (usually) a face and a point. It carries the info of a succes...
Definition: PointIndexHit.H:52
turbulentDFSEMInletFvPatchVectorField.H
Foam::Field< symmTensor >
Foam::Info
messageStream Info
Information stream (stdout output on master, null elsewhere)
Foam::fvPatch
A finiteVolume patch using a polyPatch and a fvBoundaryMesh.
Definition: fvPatch.H:65
Foam::symmTensor
SymmTensor< scalar > symmTensor
SymmTensor of scalars, i.e. SymmTensor<scalar>.
Definition: symmTensor.H:59
Foam::dictionary::getCheckOrDefault
T getCheckOrDefault(const word &keyword, const T &deflt, const Predicate &pred, enum keyType::option matchOpt=keyType::REGEX) const
Definition: dictionaryTemplates.C:209
Foam::pow
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
Definition: dimensionedScalar.C:75
Foam::max
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:47
dict
dictionary dict
Definition: searchingEngine.H:14
Foam::fvPatchField::rmap
virtual void rmap(const fvPatchField< Type > &, const labelList &)
Reverse map the given fvPatchField onto this fvPatchField.
Definition: fvPatchField.C:311
Foam::FatalError
error FatalError
Foam::dictionary
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:123
os
OBJstream os(runTime.globalPath()/outputName)
addToRunTimeSelectionTable.H
Macros for easy insertion into run-time selection tables.
Foam::PatchFunction1
Top level data entry class for use in dictionaries. Provides a mechanism to specify a variable as a c...
Definition: PatchFunction1.H:60
Foam::invariantII
Cmpt invariantII(const SymmTensor< Cmpt > &st)
Return the 2nd invariant of a SymmTensor.
Definition: SymmTensorI.H:350
Foam
Namespace for OpenFOAM.
Definition: atmBoundaryLayer.C:33
Foam::vector
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:51
Foam::exit
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:130
Foam::turbulentDFSEMInletFvPatchVectorField::rmap
virtual void rmap(const fvPatchVectorField &ptf, const labelList &addr)
Reverse map the given fvPatchField onto this fvPatchField.
Definition: turbulentDFSEMInletFvPatchVectorField.C:914
Foam::Pstream::listCombineGather
static void listCombineGather(const List< commsStruct > &comms, List< T > &Value, const CombineOp &cop, const int tag, const label comm)
Definition: combineGatherScatter.C:290
Foam::labelListList
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:56
U
U
Definition: pEqn.H:72
FatalErrorInFunction
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
Foam::Pstream::gatherList
static void gatherList(const List< commsStruct > &comms, List< T > &Values, const int tag, const label comm)
Gather data but keep individual values separate.
Definition: gatherScatterList.C:52
Foam::UPstream::allProcs
static rangeType allProcs(const label communicator=worldComm)
Range of process indices for all processes.
Definition: UPstream.H:508
Foam::sqr
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Definition: dimensionedSymmTensor.C:51
Foam::pointIndexHit
PointIndexHit< point > pointIndexHit
A PointIndexHit for 3D points.
Definition: pointIndexHit.H:46
Foam::UPstream::commsTypes::nonBlocking
Foam::UPstream::msgType
static int & msgType() noexcept
Message tag of standard messages.
Definition: UPstream.H:540
Foam::UPstream::myProcNo
static int myProcNo(const label communicator=worldComm)
Number of this process (starting from masterNo() = 0)
Definition: UPstream.H:463
Foam::nl
constexpr char nl
Definition: Ostream.H:404
Foam::VectorSpace< Vector< scalar >, scalar, 3 >::max
static const Vector< scalar > max
Definition: VectorSpace.H:117
f
labelList f(nPoints)
Foam::foamVersion::patch
const std::string patch
OpenFOAM patch number as a std::string.
Foam::MinMax::ge
static MinMax< T > ge(const T &minVal)
A semi-infinite range from minVal to the type max.
Definition: MinMaxI.H:31
Foam::Vector< scalar >
Foam::UPstream::parRun
static bool & parRun() noexcept
Test if this a parallel run.
Definition: UPstream.H:433
Foam::List< label >
Foam::sqrt
dimensionedScalar sqrt(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:144
Foam::token::TAB
Tab [isspace].
Definition: token.H:123
Foam::mag
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
Foam::turbulentDFSEMInletFvPatchVectorField
The turbulentDFSEMInlet is a synthesised-eddy based velocity inlet boundary condition to generate syn...
Definition: turbulentDFSEMInletFvPatchVectorField.H:258
points
const pointField & points
Definition: gmvOutputHeader.H:1
k
label k
Boltzmann constant.
Definition: LISASMDCalcMethod2.H:41
Foam::constant::electromagnetic::e
const dimensionedScalar e
Elementary charge.
Definition: createFields.H:11
Foam::identity
labelList identity(const label len, label start=0)
Create identity map of the given length with (map[i] == i)
Definition: labelList.C:38
Foam::sum
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
Definition: DimensionedFieldFunctions.C:327
path
fileName path(UMean.rootPath()/UMean.caseName()/"graphs"/UMean.instance())
Foam::search
fileName search(const word &file, const fileName &directory)
Recursively search the given directory for the file.
Definition: fileName.C:571
x
x
Definition: LISASMDCalcMethod2.H:52
Foam::dictionary::getCheck
T getCheck(const word &keyword, const Predicate &pred, enum keyType::option matchOpt=keyType::REGEX) const
Definition: dictionaryTemplates.C:120
Foam::Ostream::writeEntry
Ostream & writeEntry(const keyType &key, const T &value)
Write a keyword/value entry.
Definition: Ostream.H:236
Foam::turbulentDFSEMInletFvPatchVectorField::updateCoeffs
virtual void updateCoeffs()
Update the coefficients associated with the patch field.
Definition: turbulentDFSEMInletFvPatchVectorField.C:941
Foam::constant::universal::c
const dimensionedScalar c
Speed of light in a vacuum.
Foam::name
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
Definition: exprTraits.C:59
forAllReverse
#define forAllReverse(list, i)
Reverse loop across all elements in list.
Definition: stdFoam.H:309
Foam::turbulentDFSEMInletFvPatchVectorField::autoMap
virtual void autoMap(const fvPatchFieldMapper &m)
Map (and resize as needed) from self given a mapping object.
Definition: turbulentDFSEMInletFvPatchVectorField.C:890
Foam::fvPatchFieldMapper
Foam::fvPatchFieldMapper.
Definition: fvPatchFieldMapper.H:47
timeIndex
label timeIndex
Definition: getTimeIndex.H:30
Foam::tr
dimensionedScalar tr(const dimensionedSphericalTensor &dt)
Definition: dimensionedSphericalTensor.C:51
Foam::Ostream
An Ostream is an abstract base class for all output systems (streams, files, token lists,...
Definition: Ostream.H:56
Foam::makePatchTypeField
makePatchTypeField(fvPatchScalarField, atmBoundaryLayerInletEpsilonFvPatchScalarField)
Foam::gMin
Type gMin(const FieldField< Field, Type > &f)
Definition: FieldFieldFunctions.C:593
Foam::Pstream::listCombineScatter
static void listCombineScatter(const List< commsStruct > &comms, List< T > &Value, const int tag, const label comm)
Scatter data. Reverse of combineGather.
Definition: combineGatherScatter.C:432
Foam::point
vector point
Point is a vector.
Definition: point.H:43
Foam::cbrt
dimensionedScalar cbrt(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:155
Foam::dictionary::getOrDefault
T getOrDefault(const word &keyword, const T &deflt, enum keyType::option matchOpt=keyType::REGEX) const
Definition: dictionaryTemplates.C:148
UMean
volVectorField UMean(UMeanHeader, mesh)
triFace
face triFace(3)
WarningInFunction
#define WarningInFunction
Report a warning using Foam::Warning.
Definition: messageStream.H:328
Foam::fvPatchField::autoMap
virtual void autoMap(const fvPatchFieldMapper &)
Map (and resize as needed) from self given a mapping object.
Definition: fvPatchField.C:248
Foam::gMax
Type gMax(const FieldField< Field, Type > &f)
Definition: FieldFieldFunctions.C:592
Foam::UPstream::nProcs
static label nProcs(const label communicator=worldComm)
Number of processes in parallel run, and 1 for serial run.
Definition: UPstream.H:445
Foam::DimensionedField
Field with dimensions and associated with geometry type GeoMesh which is used to size the field and a...
Definition: DimensionedField.H:54
Foam::triPointRef
triangle< point, const point & > triPointRef
A triangle using referred points.
Definition: triangle.H:71
Foam::dev
dimensionedSymmTensor dev(const dimensionedSymmTensor &dt)
Definition: dimensionedSymmTensor.C:106
Foam::pos
dimensionedScalar pos(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:177