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-2022 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  sumTriMagSf_ = Zero;
165  sumTriMagSf_[Pstream::myProcNo() + 1] = sum(triMagSf);
166 
167  Pstream::listCombineReduce(sumTriMagSf_, maxEqOp<scalar>());
168 
169  for (label i = 1; i < triMagSf.size(); ++i)
170  {
171  triMagSf[i] += triMagSf[i-1];
172  }
173 
174  // Transfer to persistent storage
175  triFace_.transfer(triFace);
176  triToFace_.transfer(triToFace);
177  triCumulativeMagSf_.transfer(triMagSf);
178 
179  // Convert sumTriMagSf_ into cumulative sum of areas per proc
180  for (label i = 1; i < sumTriMagSf_.size(); ++i)
181  {
182  sumTriMagSf_[i] += sumTriMagSf_[i-1];
183  }
184 
185  // Global patch area (over all processors)
186  patchArea_ = sumTriMagSf_.last();
187 
188  // Local patch bounds (this processor)
189  patchBounds_ = boundBox(patch.localPoints(), false);
190  patchBounds_.inflate(0.1);
191 
192  // Determine if all eddies spawned from a single processor
193  singleProc_ = returnReduceOr
194  (
195  patch.size() == returnReduce(patch.size(), sumOp<label>())
196  );
197 }
198 
199 
200 void Foam::turbulentDFSEMInletFvPatchVectorField::initialiseEddyBox()
201 {
202  const scalarField& magSf = patch().magSf();
203 
204  const scalarField L(L_->value(db().time().timeOutputValue())/Lref_);
205 
206  // (PCF:Eq. 14)
207  const scalarField cellDx(max(Foam::sqrt(magSf), 2/patch().deltaCoeffs()));
208 
209  // Inialise eddy box extents
210  forAll(*this, faceI)
211  {
212  scalar& s = sigmax_[faceI];
213 
214  // Average length scale (SST:Eq. 24)
215  // Personal communication regarding (PCR:Eq. 14)
216  // - the min operator in Eq. 14 is a typo, and should be a max operator
217  s = min(mag(L[faceI]), kappa_*delta_);
218  s = max(s, nCellPerEddy_*cellDx[faceI]);
219  }
220 
221  // Maximum extent across all processors
222  maxSigmaX_ = gMax(sigmax_);
223 
224  // Eddy box volume
225  v0_ = 2*gSum(magSf)*maxSigmaX_;
226 
227  if (debug)
228  {
229  Info<< "Patch: " << patch().patch().name() << " eddy box:" << nl
230  << " volume : " << v0_ << nl
231  << " maxSigmaX : " << maxSigmaX_ << nl
232  << endl;
233  }
234 }
235 
236 
237 Foam::pointIndexHit Foam::turbulentDFSEMInletFvPatchVectorField::setNewPosition
238 (
239  const bool global
240 )
241 {
242  // Initialise to miss
243  pointIndexHit pos(false, vector::max, -1);
244 
245  const polyPatch& patch = this->patch().patch();
246  const pointField& points = patch.points();
247 
248  if (global)
249  {
250  const scalar areaFraction =
251  rndGen_.globalPosition<scalar>(0, patchArea_);
252 
253  // Determine which processor to use
254  label procI = 0;
255  forAllReverse(sumTriMagSf_, i)
256  {
257  if (areaFraction >= sumTriMagSf_[i])
258  {
259  procI = i;
260  break;
261  }
262  }
263 
264  if (Pstream::myProcNo() == procI)
265  {
266  // Find corresponding decomposed face triangle
267  label triI = 0;
268  const scalar offset = sumTriMagSf_[procI];
269  forAllReverse(triCumulativeMagSf_, i)
270  {
271  if (areaFraction > triCumulativeMagSf_[i] + offset)
272  {
273  triI = i;
274  break;
275  }
276  }
277 
278  // Find random point in triangle
279  const face& tf = triFace_[triI];
280  const triPointRef tri(points[tf[0]], points[tf[1]], points[tf[2]]);
281 
282  pos.hitPoint(tri.randomPoint(rndGen_));
283  pos.setIndex(triToFace_[triI]);
284  }
285  }
286  else
287  {
288  // Find corresponding decomposed face triangle on local processor
289  label triI = 0;
290  const scalar maxAreaLimit = triCumulativeMagSf_.last();
291  const scalar areaFraction = rndGen_.position<scalar>(0, maxAreaLimit);
292 
293  forAllReverse(triCumulativeMagSf_, i)
294  {
295  if (areaFraction > triCumulativeMagSf_[i])
296  {
297  triI = i;
298  break;
299  }
300  }
301 
302  // Find random point in triangle
303  const face& tf = triFace_[triI];
304  const triPointRef tri(points[tf[0]], points[tf[1]], points[tf[2]]);
305 
306  pos.hitPoint(tri.randomPoint(rndGen_));
307  pos.setIndex(triToFace_[triI]);
308  }
309 
310  return pos;
311 }
312 
313 
314 void Foam::turbulentDFSEMInletFvPatchVectorField::initialiseEddies()
315 {
316  const scalar t = db().time().timeOutputValue();
317  const symmTensorField R(R_->value(t)/sqr(Uref_));
318 
319  DynamicList<eddy> eddies(size());
320 
321  // Initialise eddy properties
322  scalar sumVolEddy = 0;
323  scalar sumVolEddyAllProc = 0;
324 
325  while (sumVolEddyAllProc/v0_ < d_)
326  {
327  bool search = true;
328  label iter = 0;
329 
330  while (search && iter++ < seedIterMax_)
331  {
332  // Get new parallel consistent position
333  pointIndexHit pos(setNewPosition(true));
334  const label patchFaceI = pos.index();
335 
336  // Note: only 1 processor will pick up this face
337  if (patchFaceI != -1)
338  {
339  eddy e
340  (
341  patchFaceI,
342  pos.hitPoint(),
343  rndGen_.position<scalar>(-maxSigmaX_, maxSigmaX_),
344  sigmax_[patchFaceI],
345  R[patchFaceI],
346  rndGen_
347  );
348 
349  // If eddy valid, patchFaceI is non-zero
350  if (e.patchFaceI() != -1)
351  {
352  eddies.append(e);
353  sumVolEddy += e.volume();
354  search = false;
355  }
356  }
357  // else eddy on remote processor
358 
359  reduce(search, andOp<bool>());
360  }
361 
362 
363  sumVolEddyAllProc = returnReduce(sumVolEddy, sumOp<scalar>());
364  }
365  eddies_.transfer(eddies);
366 
367  nEddy_ = eddies_.size();
368 
369  if (debug)
370  {
371  Pout<< "Patch:" << patch().patch().name();
372 
373  if (Pstream::parRun())
374  {
375  Pout<< " processor:" << Pstream::myProcNo();
376  }
377 
378  Pout<< " seeded:" << nEddy_ << " eddies" << endl;
379  }
380 
381  reduce(nEddy_, sumOp<label>());
382 
383  if (nEddy_ > 0)
384  {
385  Info<< "Turbulent DFSEM patch: " << patch().name()
386  << " seeded " << nEddy_ << " eddies with total volume "
387  << sumVolEddyAllProc
388  << endl;
389  }
390  else
391  {
393  << "Patch: " << patch().patch().name()
394  << " on field " << internalField().name()
395  << ": No eddies seeded - please check your set-up"
396  << endl;
397  }
398 }
399 
400 
401 void Foam::turbulentDFSEMInletFvPatchVectorField::convectEddies
402 (
403  const vector& UBulk,
404  const scalar deltaT
405 )
406 {
407  const scalar t = db().time().timeOutputValue();
408  const symmTensorField R(R_->value(t)/sqr(Uref_));
409 
410  // Note: all operations applied to local processor only
411 
412  label nRecycled = 0;
413 
414  forAll(eddies_, eddyI)
415  {
416  eddy& e = eddies_[eddyI];
417  e.move(deltaT*(UBulk & patchNormal_));
418 
419  const scalar position0 = e.x();
420 
421  // Check to see if eddy has exited downstream box plane
422  if (position0 > maxSigmaX_)
423  {
424  bool search = true;
425  label iter = 0;
426 
427  while (search && iter++ < seedIterMax_)
428  {
429  // Spawn new eddy with new random properties (intensity etc)
430  pointIndexHit pos(setNewPosition(false));
431  const label patchFaceI = pos.index();
432 
433  e = eddy
434  (
435  patchFaceI,
436  pos.hitPoint(),
437  position0 - floor(position0/maxSigmaX_)*maxSigmaX_,
438  sigmax_[patchFaceI],
439  R[patchFaceI],
440  rndGen_
441  );
442 
443  if (e.patchFaceI() != -1)
444  {
445  search = false;
446  }
447  }
448 
449  nRecycled++;
450  }
451  }
452 
453  if (debug)
454  {
455  reduce(nRecycled, sumOp<label>());
456 
457  if (nRecycled)
458  {
459  Info<< "Patch: " << patch().patch().name()
460  << " recycled " << nRecycled << " eddies"
461  << endl;
462  }
463  }
464 }
465 
466 
467 Foam::vector Foam::turbulentDFSEMInletFvPatchVectorField::uPrimeEddy
468 (
469  const List<eddy>& eddies,
470  const point& patchFaceCf
471 ) const
472 {
473  vector uPrime(Zero);
474 
475  forAll(eddies, k)
476  {
477  const eddy& e = eddies[k];
478  uPrime += e.uPrime(patchFaceCf, patchNormal_);
479  }
480 
481  return uPrime;
482 }
483 
484 
485 void Foam::turbulentDFSEMInletFvPatchVectorField::calcOverlappingProcEddies
486 (
487  List<List<eddy>>& overlappingEddies
488 ) const
489 {
490  int oldTag = UPstream::msgType();
491  UPstream::msgType() = oldTag + 1;
492 
493  List<boundBox> patchBBs(Pstream::nProcs());
494  patchBBs[Pstream::myProcNo()] = patchBounds_;
495  Pstream::allGatherList(patchBBs);
496 
497  // Per processor indices into all segments to send
498  List<DynamicList<label>> sendMap(UPstream::nProcs());
499 
500  // Collect overlapping eddies
501  forAll(eddies_, i)
502  {
503  const eddy& e = eddies_[i];
504 
505  // Eddy bounds
506  const point x(e.position(patchNormal_));
507  boundBox ebb(e.bounds());
508  ebb.min() += x;
509  ebb.max() += x;
510 
511  forAll(patchBBs, proci)
512  {
513  // Not including intersection with local patch
514  if (proci != Pstream::myProcNo())
515  {
516  if (ebb.overlaps(patchBBs[proci]))
517  {
518  sendMap[proci].push_back(i);
519  }
520  }
521  }
522  }
523 
524 
525  PstreamBuffers pBufs(UPstream::commsTypes::nonBlocking);
526 
527  forAll(sendMap, proci)
528  {
529  if (proci != Pstream::myProcNo() && !sendMap[proci].empty())
530  {
531  UOPstream os(proci, pBufs);
532 
533  os << UIndirectList<eddy>(eddies_, sendMap[proci]);
534  }
535  }
536 
537  pBufs.finishedSends();
538 
539  for (const int proci : pBufs.allProcs())
540  {
541  if (proci != Pstream::myProcNo() && pBufs.recvDataCount(proci))
542  {
543  UIPstream is(proci, pBufs);
544 
545  is >> overlappingEddies[proci];
546  }
547  }
548 
549  // Restore tag
550  UPstream::msgType() = oldTag;
551 }
552 
553 
554 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
555 
558 (
559  const fvPatch& p,
561 )
562 :
564  U_(nullptr),
565  R_(nullptr),
566  L_(nullptr),
567  delta_(1.0),
568  d_(1.0),
569  kappa_(0.41),
570  Uref_(1.0),
571  Lref_(1.0),
572  scale_(1.0),
573  m_(0.5),
574  nCellPerEddy_(5),
575 
576  patchArea_(-1),
577  triFace_(),
578  triToFace_(),
579  triCumulativeMagSf_(),
580  sumTriMagSf_(Pstream::nProcs() + 1, Zero),
581  patchNormal_(Zero),
582  patchBounds_(),
583 
584  eddies_(Zero),
585  v0_(Zero),
586  rndGen_(Pstream::myProcNo()),
587  sigmax_(size(), Zero),
588  maxSigmaX_(Zero),
589  nEddy_(0),
590  curTimeIndex_(-1),
591  singleProc_(false),
592  writeEddies_(false)
593 {}
594 
595 
598 (
600  const fvPatch& p,
602  const fvPatchFieldMapper& mapper
603 )
604 :
605  fixedValueFvPatchField<vector>(ptf, p, iF, mapper),
606  U_(ptf.U_.clone(patch().patch())),
607  R_(ptf.R_.clone(patch().patch())),
608  L_(ptf.L_.clone(patch().patch())),
609  delta_(ptf.delta_),
610  d_(ptf.d_),
611  kappa_(ptf.kappa_),
612  Uref_(ptf.Uref_),
613  Lref_(ptf.Lref_),
614  scale_(ptf.scale_),
615  m_(ptf.m_),
616  nCellPerEddy_(ptf.nCellPerEddy_),
617 
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_),
625 
626  eddies_(ptf.eddies_),
627  v0_(ptf.v0_),
628  rndGen_(ptf.rndGen_),
629  sigmax_(ptf.sigmax_, mapper),
630  maxSigmaX_(ptf.maxSigmaX_),
631  nEddy_(ptf.nEddy_),
632  curTimeIndex_(-1),
633  singleProc_(ptf.singleProc_),
634  writeEddies_(ptf.writeEddies_)
635 {}
636 
637 
640 (
641  const fvPatch& p,
643  const dictionary& dict
644 )
645 :
647  U_(PatchFunction1<vector>::New(patch().patch(), "U", dict)),
648  R_(PatchFunction1<symmTensor>::New(patch().patch(), "R", dict)),
649  L_(PatchFunction1<scalar>::New(patch().patch(), "L", dict)),
650  delta_(dict.getCheck<scalar>("delta", scalarMinMax::ge(0))),
651  d_(dict.getCheckOrDefault<scalar>("d", 1, scalarMinMax::ge(SMALL))),
652  kappa_(dict.getCheckOrDefault<scalar>("kappa", 0.41, scalarMinMax::ge(0))),
653  Uref_(dict.getCheckOrDefault<scalar>("Uref", 1, scalarMinMax::ge(SMALL))),
654  Lref_(dict.getCheckOrDefault<scalar>("Lref", 1, scalarMinMax::ge(SMALL))),
655  scale_(dict.getCheckOrDefault<scalar>("scale", 1, scalarMinMax::ge(0))),
656  m_(dict.getCheckOrDefault<scalar>("m", 0.5, scalarMinMax::ge(0))),
657  nCellPerEddy_(dict.getOrDefault<label>("nCellPerEddy", 5)),
658 
659  patchArea_(-1),
660  triFace_(),
661  triToFace_(),
662  triCumulativeMagSf_(),
663  sumTriMagSf_(Pstream::nProcs() + 1, Zero),
664  patchNormal_(Zero),
665  patchBounds_(),
666 
667  eddies_(),
668  v0_(Zero),
669  rndGen_(),
670  sigmax_(size(), Zero),
671  maxSigmaX_(Zero),
672  nEddy_(0),
673  curTimeIndex_(-1),
674  singleProc_(false),
675  writeEddies_(dict.getOrDefault("writeEddies", false))
676 {
677  eddy::debug = debug;
678 
679  const scalar t = db().time().timeOutputValue();
680  const symmTensorField R(R_->value(t)/sqr(Uref_));
682  checkStresses(R);
683 }
684 
685 
688 (
690 )
691 :
693  U_(ptf.U_.clone(patch().patch())),
694  R_(ptf.R_.clone(patch().patch())),
695  L_(ptf.L_.clone(patch().patch())),
696  delta_(ptf.delta_),
697  d_(ptf.d_),
698  kappa_(ptf.kappa_),
699  Uref_(ptf.Uref_),
700  Lref_(ptf.Lref_),
701  scale_(ptf.scale_),
702  m_(ptf.m_),
703  nCellPerEddy_(ptf.nCellPerEddy_),
704 
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_),
712 
713  eddies_(ptf.eddies_),
714  v0_(ptf.v0_),
715  rndGen_(ptf.rndGen_),
716  sigmax_(ptf.sigmax_),
717  maxSigmaX_(ptf.maxSigmaX_),
718  nEddy_(ptf.nEddy_),
719  curTimeIndex_(-1),
720  singleProc_(ptf.singleProc_),
721  writeEddies_(ptf.writeEddies_)
722 {}
723 
724 
727 (
730 )
731 :
732  fixedValueFvPatchField<vector>(ptf, iF),
733  U_(ptf.U_.clone(patch().patch())),
734  R_(ptf.R_.clone(patch().patch())),
735  L_(ptf.L_.clone(patch().patch())),
736  delta_(ptf.delta_),
737  d_(ptf.d_),
738  kappa_(ptf.kappa_),
739  Uref_(ptf.Uref_),
740  Lref_(ptf.Lref_),
741  scale_(ptf.scale_),
742  m_(ptf.m_),
743  nCellPerEddy_(ptf.nCellPerEddy_),
744 
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_),
752 
753  eddies_(ptf.eddies_),
754  v0_(ptf.v0_),
755  rndGen_(ptf.rndGen_),
756  sigmax_(ptf.sigmax_),
757  maxSigmaX_(ptf.maxSigmaX_),
758  nEddy_(ptf.nEddy_),
759  curTimeIndex_(-1),
760  singleProc_(ptf.singleProc_),
761  writeEddies_(ptf.writeEddies_)
762 {}
763 
764 
765 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
766 
768 (
769  const symmTensorField& R
770 )
771 {
772  constexpr label maxDiffs = 5;
773  label nDiffs = 0;
774 
775  // (S:Eq. 4a-4c)
776  forAll(R, i)
777  {
778  bool diff = false;
779 
780  if (maxDiffs < nDiffs)
781  {
782  Info<< "More than " << maxDiffs << " times"
783  << " Reynolds-stress realizability checks failed."
784  << " Skipping further comparisons." << endl;
785  return;
786  }
787 
788  const symmTensor& r = R[i];
789 
790  if (r.xx() < 0)
791  {
793  << "Reynolds stress " << r << " at index " << i
794  << " does not obey the constraint: Rxx >= 0"
795  << endl;
796  diff = true;
797  }
798 
799  if ((r.xx()*r.yy() - sqr(r.xy())) < 0)
800  {
802  << "Reynolds stress " << r << " at index " << i
803  << " does not obey the constraint: Rxx*Ryy - sqr(Rxy) >= 0"
804  << endl;
805  diff = true;
806  }
807 
808  if (det(r) < 0)
809  {
811  << "Reynolds stress " << r << " at index " << i
812  << " does not obey the constraint: det(R) >= 0"
813  << endl;
814  diff = true;
815  }
816 
817  if (diff)
818  {
819  ++nDiffs;
820  }
821  }
822 }
823 
824 
826 (
827  const scalarField& R
828 )
829 {
830  if (min(R) <= 0)
831  {
833  << "Reynolds stresses contain at least one "
834  << "nonpositive element. min(R) = " << min(R)
835  << exit(FatalError);
836  }
837 }
838 
839 
841 (
842  const fvPatchFieldMapper& m
843 )
844 {
846 
847  if (U_)
848  {
849  U_->autoMap(m);
850  }
851  if (R_)
852  {
853  R_->autoMap(m);
854  }
855  if (L_)
856  {
857  L_->autoMap(m);
858  }
859 
860  sigmax_.autoMap(m);
861 }
862 
863 
865 (
866  const fvPatchVectorField& ptf,
867  const labelList& addr
868 )
869 {
871 
872  const auto& dfsemptf =
873  refCast<const turbulentDFSEMInletFvPatchVectorField>(ptf);
874 
875  if (U_)
876  {
877  U_->rmap(dfsemptf.U_(), addr);
878  }
879  if (R_)
880  {
881  R_->rmap(dfsemptf.R_(), addr);
882  }
883  if (L_)
884  {
885  L_->rmap(dfsemptf.L_(), addr);
886  }
887 
888  sigmax_.rmap(dfsemptf.sigmax_, addr);
889 }
890 
891 
893 {
894  if (updated())
895  {
896  return;
897  }
898 
899  if (curTimeIndex_ == -1)
900  {
901  initialisePatch();
902 
903  initialiseEddyBox();
904 
905  initialiseEddies();
906  }
907 
908 
909  if (curTimeIndex_ != db().time().timeIndex())
910  {
911  tmp<vectorField> UMean =
912  U_->value(db().time().timeOutputValue())/Uref_;
913 
914  // (PCR:p. 522)
915  const vector UBulk
916  (
917  gSum(UMean()*patch().magSf())
918  /(gSum(patch().magSf()) + ROOTVSMALL)
919  );
920 
921  // Move eddies using bulk velocity
922  const scalar deltaT = db().time().deltaTValue();
923  convectEddies(UBulk, deltaT);
924 
925  // Set mean velocity
926  vectorField& U = *this;
927  U = UMean;
928 
929  // Apply second part of normalisation coefficient
930  const scalar c =
931  scale_*Foam::pow(10*v0_, m_)/Foam::sqrt(scalar(nEddy_));
932 
933  // In parallel, need to collect all eddies that will interact with
934  // local faces
935 
936  const pointField& Cf = patch().Cf();
937 
938  if (singleProc_ || !Pstream::parRun())
939  {
940  forAll(U, faceI)
941  {
942  U[faceI] += c*uPrimeEddy(eddies_, Cf[faceI]);
943  }
944  }
945  else
946  {
947  // Process local eddy contributions
948  forAll(U, faceI)
949  {
950  U[faceI] += c*uPrimeEddy(eddies_, Cf[faceI]);
951  }
952 
953  // Add contributions from overlapping eddies
954  List<List<eddy>> overlappingEddies(Pstream::nProcs());
955  calcOverlappingProcEddies(overlappingEddies);
956 
957  for (const List<eddy>& eddies : overlappingEddies)
958  {
959  if (eddies.size())
960  {
961  forAll(U, faceI)
962  {
963  U[faceI] += c*uPrimeEddy(eddies, Cf[faceI]);
964  }
965  }
966  }
967  }
968 
969  // Re-scale to ensure correct flow rate
970  const scalar fCorr =
971  gSum((UBulk & patchNormal_)*patch().magSf())
972  /gSum(U & -patch().Sf());
973 
974  U *= fCorr;
975 
976  curTimeIndex_ = db().time().timeIndex();
977 
978  if (writeEddies_)
979  {
980  writeEddyOBJ();
981  }
982 
983  if (debug)
984  {
985  Info<< "Magnitude of bulk velocity: " << UBulk << endl;
986 
987  Info<< "Number of eddies: "
988  << returnReduce(eddies_.size(), sumOp<label>())
989  << endl;
990 
991  Info<< "Patch:" << patch().patch().name()
992  << " min/max(U):" << gMin(U) << ", " << gMax(U)
993  << endl;
994 
995  if (db().time().writeTime())
996  {
997  writeLumleyCoeffs();
998  }
999  }
1000  }
1001 
1002  fixedValueFvPatchVectorField::updateCoeffs();
1003 }
1004 
1005 
1007 {
1009  os.writeEntry("delta", delta_);
1010  os.writeEntryIfDifferent<scalar>("d", 1.0, d_);
1011  os.writeEntryIfDifferent<scalar>("kappa", 0.41, kappa_);
1012  os.writeEntryIfDifferent<scalar>("Uref", 1.0, Uref_);
1013  os.writeEntryIfDifferent<scalar>("Lref", 1.0, Lref_);
1014  os.writeEntryIfDifferent<scalar>("scale", 1.0, scale_);
1015  os.writeEntryIfDifferent<scalar>("m", 0.5, m_);
1016  os.writeEntryIfDifferent<label>("nCellPerEddy", 5, nCellPerEddy_);
1017  os.writeEntryIfDifferent("writeEddies", false, writeEddies_);
1018  if (U_)
1019  {
1020  U_->writeData(os);
1021  }
1022  if (R_)
1023  {
1024  R_->writeData(os);
1025  }
1026  if (L_)
1027  {
1028  L_->writeData(os);
1029  }
1031 }
1032 
1033 
1034 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
1035 
1036 namespace Foam
1037 {
1039  (
1041  turbulentDFSEMInletFvPatchVectorField
1042  );
1043 }
1044 
1045 
1046 // ************************************************************************* //
Cmpt invariantII(const SymmTensor< Cmpt > &st)
Return the 2nd invariant of a SymmTensor.
Definition: SymmTensorI.H:581
scalar diff(const triad &A, const triad &B)
Return a quantity of the difference between two triads.
Definition: triad.C:373
This boundary condition supplies a fixed value constraint, and is the base class for a number of othe...
dictionary dict
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.
Definition: fvPatchField.C:236
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:125
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
error FatalError
Error stream (stdout output on all processes), with additional &#39;FOAM FATAL ERROR&#39; header text and sta...
A list of keyword definitions, which are a keyword followed by a number of values (eg...
Definition: dictionary.H:120
const vector L(dict.get< vector >("L"))
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:578
Type gMin(const FieldField< Field, Type > &f)
Tab [isspace].
Definition: token.H:126
void append(const T &val)
Append an element at the end of the list.
Definition: List.H:500
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:40
dimensionedSymmTensor sqr(const dimensionedVector &dv)
constexpr char nl
The newline &#39;\n&#39; character (0x0a)
Definition: Ostream.H:49
dimensionedScalar sqrt(const dimensionedScalar &ds)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:487
A finiteVolume patch using a polyPatch and a fvBoundaryMesh.
Definition: fvPatch.H:70
PointIndexHit< point > pointIndexHit
A PointIndexHit with a 3D point.
Definition: pointIndexHit.H:58
static bool & parRun() noexcept
Test if this a parallel run.
Definition: UPstream.H:1004
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.
Definition: Ostream.H:312
volVectorField UMean(UMeanHeader, mesh)
static int & msgType() noexcept
Message tag of standard messages.
Definition: UPstream.H:1184
This class describes the interaction of an object (often a face) and a point. It carries the info of ...
Definition: pointIndexHit.H:44
void writeValueEntry(Ostream &os) const
Write *this field as a "value" entry.
Definition: fvPatchField.H:372
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...
Definition: UPstream.H:1029
virtual void write(Ostream &) const
Write.
Definition: fvPatchField.C:372
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.
Definition: stdFoam.H:414
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. ...
Definition: UPstream.H:1020
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:38
const dimensionedScalar e
Elementary charge.
Definition: createFields.H:11
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.
Definition: exprTraits.C:52
const pointField & points
SymmTensor< scalar > symmTensor
SymmTensor of scalars, i.e. SymmTensor<scalar>.
Definition: symmTensor.H:55
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Cmpt invariantIII(const SymmTensor< Cmpt > &st)
Return the 3rd invariant of a SymmTensor.
Definition: SymmTensorI.H:595
Inter-processor communications stream.
Definition: Pstream.H:57
face triFace(3)
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.
Definition: fvPatchField.C:299
dimensionedScalar tr(const dimensionedSphericalTensor &dt)
Vector< scalar > vector
Definition: vector.H:57
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:26
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:327
static int debug
Flag to activate debug statements.
Definition: eddy.H:172
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())
labelList f(nPoints)
turbulentDFSEMInletFvPatchVectorField(const fvPatch &, const DimensionedField< vector, volMesh > &)
Construct from patch and internal field.
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
U
Definition: pEqn.H:72
vector point
Point is a vector.
Definition: point.H:37
#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...
Definition: areaFieldsFwd.H:42
const std::string patch
OpenFOAM patch number as a std::string.
triangle< point, const point & > triPointRef
A triangle using referred points.
Definition: triangleFwd.H:39
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.
Definition: fileName.C:640
#define forAllReverse(list, i)
Reverse loop across all elements in list.
Definition: stdFoam.H:430
List< label > labelList
A List of labels.
Definition: List.H:62
volScalarField & p
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())
Namespace for OpenFOAM.
label timeIndex
Definition: getTimeIndex.H:24
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.
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:133