turbulentDigitalFilterInletFvPatchVectorField.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) 2019-2020 OpenCFD Ltd.
9 -------------------------------------------------------------------------------
10 License
11  This file is part of OpenFOAM.
12 
13  OpenFOAM is free software: you can redistribute it and/or modify it
14  under the terms of the GNU General Public License as published by
15  the Free Software Foundation, either version 3 of the License, or
16  (at your option) any later version.
17 
18  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
19  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
20  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
21  for more details.
22 
23  You should have received a copy of the GNU General Public License
24  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
25 
26 \*---------------------------------------------------------------------------*/
27 
29 #include "volFields.H"
30 #include "mathematicalConstants.H"
32 
33 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
34 
36 Foam::turbulentDigitalFilterInletFvPatchVectorField::patchMapper() const
37 {
38  // Initialise interpolation (2D planar interpolation by triangulation)
39  if (!mapperPtr_)
40  {
41  // Reread values and interpolate
42  const fileName samplePointsFile
43  (
44  this->db().time().globalPath()
45  /this->db().time().constant()
46  /"boundaryData"
47  /this->patch().name()
48  /"points"
49  );
50 
51  IOobject io
52  (
53  samplePointsFile, // absolute path
54  this->db().time(),
57  false, // no need to register
58  true // is global object (currently not used)
59  );
60 
61  // Read data
62  const rawIOField<point> samplePoints(io, false);
63 
64  // tbd: run-time selection
65  bool nearestOnly =
66  (
67  !mapMethod_.empty()
68  && mapMethod_ != "planarInterpolation"
69  );
70 
71  // Allocate the interpolator
72  mapperPtr_.reset
73  (
74  new pointToPointPlanarInterpolation
75  (
76  samplePoints,
77  this->patch().patch().faceCentres(),
78  perturb_,
79  nearestOnly
80  )
81  );
82  }
83 
84  return *mapperPtr_;
85 }
86 
87 
89 Foam::turbulentDigitalFilterInletFvPatchVectorField::indexPairs()
90 {
91  // Get patch normal direction into the domain
92  const vector nf((-gAverage(patch().nf())).normalise());
93 
94  // Find the second local coordinate direction
95  direction minCmpt = 0;
96  scalar minMag = mag(nf[minCmpt]);
97  for (direction cmpt = 1; cmpt < pTraits<vector>::nComponents; ++cmpt)
98  {
99  const scalar s = mag(nf[cmpt]);
100  if (s < minMag)
101  {
102  minMag = s;
103  minCmpt = cmpt;
104  }
105  }
106 
107  // Create the second local coordinate direction
108  vector e2(Zero);
109  e2[minCmpt] = 1;
110 
111  // Remove normal component
112  e2 -= (nf&e2)*nf;
113 
114  // Create the local coordinate system
115  coordSystem::cartesian cs
116  (
117  Zero, // origin
118  nf, // normal
119  e2 // 0-axis
120  );
121 
122  // Convert patch points into local coordinate system
123  const pointField localPos
124  (
125  cs.localPosition
126  (
127  pointField
128  (
129  patch().patch().points(),
130  patch().patch().meshPoints()
131  )
132  )
133  );
134 
135  const boundBox bb(localPos);
136 
137  // Normalise to (i, j) coordinates
138  const Vector2D<label> n(n_.first(), n_.second());
139  invDelta_[0] = n[0]/bb.span()[0];
140  invDelta_[1] = n[1]/bb.span()[1];
141  const point localMinPt(bb.min());
142 
143  // Compute virtual-actual patch index pairs
144  List<Pair<label>> indexPairs(this->size(), Pair<label>(Zero, Zero));
145 
146  forAll(*this, facei)
147  {
148  const scalar& centre0 = localPos[facei][0];
149  const scalar& centre1 = localPos[facei][1];
150 
151  // Virtual turbulence plane indices
152  const label j = label((centre0 - localMinPt[0])*invDelta_[0]);
153  const label k = label((centre1 - localMinPt[1])*invDelta_[1]);
154 
155  indexPairs[facei] = Pair<label>(facei, k*n[0] + j);
156  }
157 
158  return indexPairs;
159 }
160 
161 
162 void Foam::turbulentDigitalFilterInletFvPatchVectorField::checkR() const
163 {
164  const vectorField& faceCentres = this->patch().patch().faceCentres();
165 
166  forAll(R_, facei)
167  {
168  if (R_[facei].xx() <= 0)
169  {
171  << "Reynolds stress tensor component Rxx cannot be negative"
172  << " or zero, where Rxx = " << R_[facei].xx()
173  << " at the face centre = " << faceCentres[facei]
174  << exit(FatalError);
175  }
176 
177  if (R_[facei].yy() < 0 || R_[facei].zz() < 0)
178  {
180  << "Reynolds stress tensor components Ryy or Rzz cannot be"
181  << " negative where Ryy = " << R_[facei].yy()
182  << ", and Rzz = " << R_[facei].zz()
183  << " at the face centre = " << faceCentres[facei]
184  << exit(FatalError);
185  }
186 
187  const scalar x0 = R_[facei].xx()*R_[facei].yy() - sqr(R_[facei].xy());
188 
189  if (x0 <= 0)
190  {
192  << "Reynolds stress tensor component group, Rxx*Ryy - Rxy^2"
193  << " cannot be negative or zero"
194  << " at the face centre = " << faceCentres[facei]
195  << exit(FatalError);
196  }
197 
198  const scalar x1 = R_[facei].zz() - sqr(R_[facei].xz())/R_[facei].xx();
199  const scalar x2 = sqr(R_[facei].yz() - R_[facei].xy()*R_[facei].xz()
200  /(R_[facei].xx()*x0));
201  const scalar x3 = x1 - x2;
202 
203  if (x3 < 0)
204  {
206  << "Reynolds stress tensor component group, "
207  << "Rzz - Rxz^2/Rxx - (Ryz - Rxy*Rxz/(Rxx*(Rxx*Ryy - Rxy^2)))^2"
208  << " cannot be negative at the face centre = "
209  << faceCentres[facei]
210  << exit(FatalError);
211  }
212  }
213 
214  Info<< " # Reynolds stress tensor on patch is consistent #" << endl;
215 }
216 
217 
218 Foam::symmTensorField Foam::turbulentDigitalFilterInletFvPatchVectorField::
219 calcLund() const
220 {
221  checkR();
222 
223  symmTensorField Lund(symmTensorField(R_.size()));
224 
225  forAll(R_, facei)
226  {
227  const symmTensor& R = R_[facei];
228  symmTensor& lws = Lund[facei];
229 
230  // (KSJ:Eq. 5)
231  lws.xx() = Foam::sqrt(R.xx());
232  lws.xy() = R.xy()/lws.xx();
233  lws.xz() = R.xz()/lws.xx();
234  lws.yy() = Foam::sqrt(R.yy() - sqr(lws.xy()));
235  lws.yz() = (R.yz() - lws.xy()*lws.xz())/lws.yy();
236  lws.zz() = Foam::sqrt(R.zz() - sqr(lws.xz()) - sqr(lws.yz()));
237  }
238 
239  return Lund;
240 }
241 
242 
243 Foam::scalar Foam::turbulentDigitalFilterInletFvPatchVectorField::
244 calcFlowRate() const
245 {
246  const vector patchNormal((-gAverage(patch().nf())).normalise());
247  return gSum((UMean_ & patchNormal)*patch().magSf());
248 }
249 
250 
251 Foam::tensor Foam::turbulentDigitalFilterInletFvPatchVectorField::
252 meterToCell
253 (
254  const tensor& L
255 ) const
256 {
257  tensor Ls(L);
258 
259  // Convert streamwise integral length scales to integral
260  // time scales by using Taylor's frozen turbulence hypothesis
261  forAll(Ls.x(), i)
262  {
263  Ls[i] /= Ubulk_;
264  }
265 
266  const scalar invDeltaT = 1.0/db().time().deltaTValue();
267 
268  // (KSJ:Eq. 13)
269  Ls.row(0, Ls.x()*invDeltaT);
270  Ls.row(1, Ls.y()*invDelta_[0]);
271  Ls.row(2, Ls.z()*invDelta_[1]);
272 
273  return Ls;
274 }
275 
276 
277 Foam::List<Foam::label> Foam::turbulentDigitalFilterInletFvPatchVectorField::
278 initBox() const
279 {
280  label initValue = 0;
281  label rangeModifier = 0;
282 
283  if (fsm_)
284  {
285  // Initialise with 1 since streamwise-dir possesses 1 cell in FSM
286  initValue = 1;
287  rangeModifier = pTraits<vector>::nComponents;
288  }
289 
290  List<label> szBox(pTraits<tensor>::nComponents, initValue);
291  Vector<label> szPlane
292  (
293  1,
294  n_.first(),
295  n_.second()
296  );
297 
298  for
299  (
300  const label i
301  : labelRange(rangeModifier, pTraits<tensor>::nComponents - rangeModifier)
302  )
303  {
304  // Slicing index
305  const label slicei = label(i/pTraits<vector>::nComponents);
306 
307  // Refer to "calcFilterCoeffs()"
308  const label n = ceil(L_[i]);
309  const label twiceN = 4*n;
310 
311  // Size of random-number sets
312  szBox[i] = szPlane[slicei] + twiceN;
313  }
314 
315  return szBox;
316 }
317 
318 
319 Foam::List<Foam::label> Foam::turbulentDigitalFilterInletFvPatchVectorField::
320 initFactors2D() const
321 {
322  List<label> boxFactors2D(pTraits<vector>::nComponents);
323 
324  for (direction dir = 0; dir < pTraits<vector>::nComponents; ++dir)
325  {
326  boxFactors2D[dir] =
327  szBox_[pTraits<vector>::nComponents + dir]
328  *szBox_[pTraits<symmTensor>::nComponents + dir];
329  }
330 
331  return boxFactors2D;
332 }
333 
334 
335 Foam::List<Foam::label> Foam::turbulentDigitalFilterInletFvPatchVectorField::
336 initFactors3D() const
337 {
338  List<label> boxFactors3D(pTraits<vector>::nComponents);
339 
340  for (direction dir = 0; dir < pTraits<vector>::nComponents; ++dir)
341  {
342  boxFactors3D[dir] = boxFactors2D_[dir]*szBox_[dir];
343  }
344 
345  return boxFactors3D;
346 }
347 
348 
349 Foam::List<Foam::label> Foam::turbulentDigitalFilterInletFvPatchVectorField::
350 initPlaneFactors() const
351 {
352  List<label> planeFactors(pTraits<vector>::nComponents);
353 
354  for (direction dir = 0; dir < pTraits<vector>::nComponents; ++dir)
355  {
356  planeFactors[dir] = boxFactors2D_[dir]*(szBox_[dir] - 1);
357  }
358 
359  return planeFactors;
360 }
361 
362 
364 Foam::turbulentDigitalFilterInletFvPatchVectorField::fillBox()
365 {
366  List<List<scalar>> box(pTraits<vector>::nComponents, List<scalar>());
367 
368  // Initialise: Remaining convenience factors for (e1 e2 e3)
369  for (direction dir = 0; dir < pTraits<vector>::nComponents; ++dir)
370  {
371  // Initialise: Random-box content with random-number sets
372  box[dir] = randomSet<List<scalar>, scalar>(boxFactors3D_[dir]);
373  }
374 
375  return box;
376 }
377 
378 
380 Foam::turbulentDigitalFilterInletFvPatchVectorField::calcFilterCoeffs() const
381 {
382  List<List<scalar>> filterCoeffs
383  (
384  pTraits<tensor>::nComponents,
385  List<scalar>(1, 1.0)
386  );
387 
388  label initValue = 0;
389 
390  if (fsm_)
391  {
392  initValue = pTraits<vector>::nComponents;
393  }
394 
395  for (direction dir = initValue; dir < pTraits<tensor>::nComponents; ++dir)
396  {
397  // The smallest filter width larger than length scale
398  // (KSJ:'n' in Eq. 15)
399  const label n = ceil(L_[dir]);
400 
401  // Full filter-width (except mid-zero) according to the condition
402  // (KSJ:Eq. 15 whereat N is minimum =2n)
403  const label twiceN = 4*n;
404 
405  // Initialise filter-coeffs containers with full filter-width size
406  // Extra elem is mid-zero within [-N, N]
407  filterCoeffs[dir] = List<scalar>(twiceN + 1, Zero);
408 
409  // First element: -N within [-N, N]
410  const scalar initElem = -2*scalar(n);
411 
412  // Container initialised with [-N, N] (KSJ:p. 658, item-b)
413  std::iota
414  (
415  filterCoeffs[dir].begin(),
416  filterCoeffs[dir].end(),
417  initElem
418  );
419 
420  // Compute filter-coeffs (KSJ:Eq. 14 (Gaussian))
421  List<scalar> fTemp(filterCoeffs[dir]);
422  scalar fSum = 0;
423  const scalar nSqr = n*n;
424 
425  if (Gaussian_)
426  {
427  fTemp = sqr(exp(C1_*sqr(fTemp)/nSqr));
428  fSum = Foam::sqrt(sum(fTemp));
429 
430  filterCoeffs[dir] =
431  exp(C1_*sqr(filterCoeffs[dir])/nSqr)/fSum;
432  }
433  else
434  {
435  fTemp = sqr(exp(C1_*mag(fTemp)/n));
436  fSum = Foam::sqrt(sum(fTemp));
437 
438  filterCoeffs[dir] = exp(C1_*mag(filterCoeffs[dir])/n)/fSum;
439  }
440  }
441 
442  return filterCoeffs;
443 }
444 
445 
446 void Foam::turbulentDigitalFilterInletFvPatchVectorField::shiftRefill()
447 {
448  forAll(box_, dir)
449  {
450  List<scalar>& r = box_[dir];
451 
452  // Shift forward from the back to the front / First Out
453  inplaceRotateList(r, boxFactors2D_[dir]);
454 
455  // Refill the back with a new random-number set / First In
456  for (label i = 0; i < boxFactors2D_[dir]; ++i)
457  {
458  r[i] = rndGen_.GaussNormal<scalar>();
459  }
460  }
461 }
462 
463 
464 void Foam::turbulentDigitalFilterInletFvPatchVectorField::mapFilteredBox
465 (
466  vectorField& U
467 )
468 {
469  for (const auto& x : indexPairs_)
470  {
471  const label xf = x.first();
472  const label xs = x.second();
473  U[xf][0] = filteredBox_[0][xs];
474  U[xf][1] = filteredBox_[1][xs];
475  U[xf][2] = filteredBox_[2][xs];
476  }
477 }
478 
479 
480 void Foam::turbulentDigitalFilterInletFvPatchVectorField::onePointCorrs
481 (
482  vectorField& U
483 ) const
484 {
485  forAll(Lund_, facei)
486  {
487  vector& Us = U[facei];
488  const symmTensor& lws = Lund_[facei];
489 
490  // (KSJ:p. 658, item-e)
491  Us.z() = Us.x()*lws.xz() + Us.y()*lws.yz() + Us.z()*lws.zz();
492  Us.y() = Us.x()*lws.xy() + Us.y()*lws.yy();
493  Us.x() = Us.x()*lws.xx();
494  }
495 }
496 
497 
498 void Foam::turbulentDigitalFilterInletFvPatchVectorField::twoPointCorrs()
499 {
500  for (direction dir = 0; dir < pTraits<vector>::nComponents; ++dir)
501  {
502  List<scalar>& in = box_[dir];
503  List<scalar>& out = filteredBox_[dir];
504  const List<scalar>& filter1 = filterCoeffs_[dir];
505  const List<scalar>& filter2 = filterCoeffs_[3 + dir];
506  const List<scalar>& filter3 = filterCoeffs_[6 + dir];
507 
508  const label sz1 = szBox_[dir];
509  const label sz2 = szBox_[3 + dir];
510  const label sz3 = szBox_[6 + dir];
511  const label szfilter1 = filterCoeffs_[dir].size();
512  const label szfilter2 = filterCoeffs_[3 + dir].size();
513  const label szfilter3 = filterCoeffs_[6 + dir].size();
514  const label sz23 = boxFactors2D_[dir];
515  const label sz123 = boxFactors3D_[dir];
516  const label validSlice2 = sz2 - (szfilter2 - label(1));
517  const label validSlice3 = sz3 - (szfilter3 - label(1));
518 
519  // Convolution summation - Along 1st direction
520  scalarField tmp(sz123);
521  label filterCentre = label(szfilter2/label(2));
522  label endIndex = sz2 - filterCentre;
523  label i0 = 0;
524  label i1 = 0;
525  label i2 = 0;
526 
527  for (label i = 0; i < sz1; ++i)
528  {
529  for (label j = 0; j < sz3; ++j)
530  {
531  i1 += filterCentre;
532 
533  for (label k = filterCentre; k < endIndex; ++k)
534  {
535  tmp[i1] = 0.0;
536  label q = 0;
537 
538  for (label p = szfilter2 - 1; p >= 0; --p, ++q)
539  {
540  tmp[i1] += in[i0 + q]*filter2[p];
541  }
542  ++i0;
543  ++i1;
544 
545  }
546  i0 += (filterCentre + filterCentre);
547  i1 += filterCentre;
548  }
549  }
550 
551  // Convolution summation - Along 2nd direction
552  scalarField tmp2(tmp);
553  filterCentre = label(szfilter3/label(2));
554  endIndex = sz3 - filterCentre;
555  i1 = 0;
556  i2 = 0;
557 
558  for (label i = 0; i < sz1; ++i)
559  {
560  label sl = i*sz23;
561 
562  for (label j = 0; j < sz2; ++j)
563  {
564  i1 = j + sl;
565  i2 = i1;
566 
567  for (label k = 0; k < endIndex - filterCentre; ++k)
568  {
569  tmp[i1] = 0.0;
570  label q = 0;
571 
572  for (label p = szfilter3 - 1; p >= 0; --p, ++q)
573  {
574  tmp[i1] += tmp2[i2 + q*sz2]*filter3[p];
575  }
576  i1 += sz2;
577  i2 += sz2;
578 
579  }
580  i1 += (sz2 + filterCentre);
581  i2 += (sz2 + filterCentre);
582  }
583  }
584 
585  // Convolution summation - Along 3rd direction
586  filterCentre = label(szfilter1/label(2));
587  endIndex = sz1 - filterCentre;
588  i1 = (szfilter2 - label(1))/label(2);
589  i2 = (szfilter2 - label(1))/label(2);
590  label i3 = 0;
591 
592  for (label i = 0; i < validSlice3; ++i)
593  {
594  for (label j = 0; j < validSlice2; ++j)
595  {
596  scalar sum = 0.0;
597  i1 = i2 + j;
598 
599  for (label k = szfilter1 - 1; k >= 0; --k)
600  {
601  sum += tmp[i1]*filter1[k];
602  i1 += sz23;
603  }
604  out[i3] = sum;
605  ++i3;
606 
607  }
608  i2 += sz2;
609  }
610  }
611 }
612 
613 
614 Foam::List<Foam::scalar> Foam::turbulentDigitalFilterInletFvPatchVectorField::
615 calcCoeffs1FSM() const
616 {
617  List<scalar> coeffs1FSM(pTraits<vector>::nComponents);
618 
619  forAll(L_.x(), i)
620  {
621  coeffs1FSM[i] = exp(C1FSM_/L_[i]);
622  }
623 
624  return coeffs1FSM;
625 }
626 
627 
628 Foam::List<Foam::scalar> Foam::turbulentDigitalFilterInletFvPatchVectorField::
629 calcCoeffs2FSM() const
630 {
631  List<scalar> coeffs2FSM(pTraits<vector>::nComponents);
632 
633  forAll(L_.x(), i)
634  {
635  coeffs2FSM[i] = Foam::sqrt(1.0 - exp(C2FSM_/L_[i]));
636  }
637 
638  return coeffs2FSM;
639 }
640 
641 
642 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
643 
646 (
647  const fvPatch& p,
649 )
650 :
652  mapperPtr_(nullptr),
653  fsm_(false),
654  Gaussian_(true),
655  fixSeed_(true),
656  continuous_(false),
657  correctFlowRate_(true),
658  interpR_(false),
659  interpUMean_(false),
660  mapMethod_("nearestCell"),
661  curTimeIndex_(-1),
662  Ubulk_(0.0),
663  C1_(-0.5*constant::mathematical::pi),
664  perturb_(1e-5),
665  flowRate_(1.0),
666  rndGen_(1234),
667  n_(0, 0),
668  invDelta_(Zero),
669  indexPairs_(Zero),
670  R_(Zero),
671  Lund_(Zero),
672  UMean_(Zero),
673  Lbak_(Zero),
674  L_(Zero),
675  C1FSM_(Zero),
676  C2FSM_(Zero),
677  coeffs1FSM_(Zero),
678  coeffs2FSM_(Zero),
679  szBox_(Zero),
680  boxFactors2D_(Zero),
681  boxFactors3D_(Zero),
682  iNextToLastPlane_(Zero),
683  box_(Zero),
684  filterCoeffs_(Zero),
685  filteredBox_(Zero),
686  U0_(Zero)
687 {}
688 
689 
692 (
694  const fvPatch& p,
696  const fvPatchFieldMapper& mapper
697 )
698 :
699  fixedValueFvPatchField<vector>(ptf, p, iF, mapper),
700  mapperPtr_(nullptr),
701  fsm_(ptf.fsm_),
702  Gaussian_(ptf.Gaussian_),
703  fixSeed_(ptf.fixSeed_),
704  continuous_(ptf.continuous_),
705  correctFlowRate_(ptf.correctFlowRate_),
706  interpR_(ptf.interpR_),
707  interpUMean_(ptf.interpUMean_),
708  mapMethod_(ptf.mapMethod_),
709  curTimeIndex_(-1),
710  Ubulk_(ptf.Ubulk_),
711  C1_(ptf.C1_),
712  perturb_(ptf.perturb_),
713  flowRate_(ptf.flowRate_),
714  rndGen_(ptf.rndGen_),
715  n_(ptf.n_),
716  invDelta_(ptf.invDelta_),
717  indexPairs_(ptf.indexPairs_),
718  R_(ptf.R_),
719  Lund_(ptf.Lund_),
720  UMean_(ptf.UMean_),
721  Lbak_(ptf.Lbak_),
722  L_(ptf.L_),
723  C1FSM_(ptf.C1FSM_),
724  C2FSM_(ptf.C2FSM_),
725  coeffs1FSM_(ptf.coeffs1FSM_),
726  coeffs2FSM_(ptf.coeffs2FSM_),
727  szBox_(ptf.szBox_),
728  boxFactors2D_(ptf.boxFactors2D_),
729  boxFactors3D_(ptf.boxFactors3D_),
730  iNextToLastPlane_(ptf.iNextToLastPlane_),
731  box_(ptf.box_),
732  filterCoeffs_(ptf.filterCoeffs_),
733  filteredBox_(ptf.filteredBox_),
734  U0_(ptf.U0_)
735 {}
736 
737 
740 (
741  const fvPatch& p,
743  const dictionary& dict
744 )
745 :
747  mapperPtr_(nullptr),
748  fsm_(dict.getOrDefault("fsm", false)),
749  Gaussian_(fsm_ ? false : dict.getOrDefault("Gaussian", true)),
750  fixSeed_(dict.getOrDefault("fixSeed", true)),
751  continuous_(dict.getOrDefault("continuous", false)),
752  correctFlowRate_(dict.getOrDefault("correctFlowRate", true)),
753  interpR_(false),
754  interpUMean_(false),
755  mapMethod_(dict.getOrDefault<word>("mapMethod", "nearestCell")),
756  curTimeIndex_(-1),
757  Ubulk_(dict.getCheck<scalar>("Ubulk", scalarMinMax::ge(ROOTSMALL))),
758  C1_
759  (
760  dict.getOrDefault<scalar>("C1", -0.5*constant::mathematical::pi)
761  ),
762  perturb_
763  (
764  dict.getCheckOrDefault<scalar>("perturb", 1e-5, scalarMinMax::ge(SMALL))
765  ),
766  flowRate_(1.0),
767  rndGen_(fixSeed_ ? 1234 : time(0)),
768  n_
769  (
771  (
772  "n",
773  [&](const Tuple2<label, label>& x)
774  {
775  return min(x.first(), x.second()) > 0 ? true : false;
776  }
777  )
778  ),
779  invDelta_(Zero),
780  indexPairs_(indexPairs()),
781  R_(interpolateOrRead<symmTensor>("R", dict, interpR_)),
782  Lund_(calcLund()),
783  UMean_(interpolateOrRead<vector>("UMean", dict, interpUMean_)),
784  Lbak_
785  (
787  (
788  "L",
789  [&](const tensor& x){return cmptMin(x) > ROOTSMALL ? true : false;}
790  )
791  ),
792  L_(meterToCell(Lbak_)),
793  C1FSM_
794  (
795  fsm_
796  ? dict.getOrDefault<scalar>
797  (
798  "C1FSM",
800  )
801  : 0.0
802  ),
803  C2FSM_
804  (
805  fsm_
806  ? dict.getOrDefault<scalar>
807  (
808  "C2FSM",
810  )
811  : 0.0
812  ),
813  coeffs1FSM_(fsm_ ? calcCoeffs1FSM() : List<scalar>()),
814  coeffs2FSM_(fsm_ ? calcCoeffs2FSM() : List<scalar>()),
815  szBox_(initBox()),
816  boxFactors2D_(initFactors2D()),
817  boxFactors3D_(initFactors3D()),
818  iNextToLastPlane_(initPlaneFactors()),
819  box_
820  (
821  (continuous_ && Pstream::master())
823  (
824  "box",
825  fillBox() // First time-step
826  )
827  :
828  fillBox()
829  ),
830  filterCoeffs_(calcFilterCoeffs()),
831  filteredBox_
832  (
834  List<scalar>(n_.first()*n_.second(), Zero)
835  ),
836  U0_
837  (
838  fsm_
839  ? randomSet<vectorField, vector>(patch().size())
840  : vectorField()
841  )
842 {
843  if (correctFlowRate_)
844  {
845  flowRate_ = calcFlowRate();
846  }
847 
848  // Check if varying or fixed time-step computation
849  if (db().time().isAdjustTimeStep())
850  {
852  << " # Varying time-step computations are not fully supported #"
853  << endl;
854  }
855 
856  Info<< " # turbulentDigitalFilterInlet initialisation completed #" << endl;
857 }
858 
859 
862 (
864 )
865 :
867  mapperPtr_(nullptr),
868  fsm_(ptf.fsm_),
869  Gaussian_(ptf.Gaussian_),
870  fixSeed_(ptf.fixSeed_),
871  continuous_(ptf.continuous_),
872  correctFlowRate_(ptf.correctFlowRate_),
873  interpR_(ptf.interpR_),
874  interpUMean_(ptf.interpUMean_),
875  mapMethod_(ptf.mapMethod_),
876  curTimeIndex_(ptf.curTimeIndex_),
877  Ubulk_(ptf.Ubulk_),
878  C1_(ptf.C1_),
879  perturb_(ptf.perturb_),
880  flowRate_(ptf.flowRate_),
881  rndGen_(ptf.rndGen_),
882  n_(ptf.n_),
883  invDelta_(ptf.invDelta_),
884  indexPairs_(ptf.indexPairs_),
885  R_(ptf.R_),
886  Lund_(ptf.Lund_),
887  UMean_(ptf.UMean_),
888  Lbak_(ptf.Lbak_),
889  L_(ptf.L_),
890  C1FSM_(ptf.C1FSM_),
891  C2FSM_(ptf.C2FSM_),
892  coeffs1FSM_(ptf.coeffs1FSM_),
893  coeffs2FSM_(ptf.coeffs2FSM_),
894  szBox_(ptf.szBox_),
895  boxFactors2D_(ptf.boxFactors2D_),
896  boxFactors3D_(ptf.boxFactors3D_),
897  iNextToLastPlane_(ptf.iNextToLastPlane_),
898  box_(ptf.box_),
899  filterCoeffs_(ptf.filterCoeffs_),
900  filteredBox_(ptf.filteredBox_),
901  U0_(ptf.U0_)
902 {}
903 
904 
907 (
910 )
911 :
913  mapperPtr_(nullptr),
914  fsm_(ptf.fsm_),
915  Gaussian_(ptf.Gaussian_),
916  fixSeed_(ptf.fixSeed_),
917  continuous_(ptf.continuous_),
918  correctFlowRate_(ptf.correctFlowRate_),
919  interpR_(ptf.interpR_),
920  interpUMean_(ptf.interpUMean_),
921  mapMethod_(ptf.mapMethod_),
922  curTimeIndex_(-1),
923  Ubulk_(ptf.Ubulk_),
924  C1_(ptf.C1_),
925  perturb_(ptf.perturb_),
926  flowRate_(ptf.flowRate_),
927  rndGen_(ptf.rndGen_),
928  n_(ptf.n_),
929  invDelta_(ptf.invDelta_),
930  indexPairs_(ptf.indexPairs_),
931  R_(ptf.R_),
932  Lund_(ptf.Lund_),
933  UMean_(ptf.UMean_),
934  Lbak_(ptf.Lbak_),
935  L_(ptf.L_),
936  C1FSM_(ptf.C1FSM_),
937  C2FSM_(ptf.C2FSM_),
938  coeffs1FSM_(ptf.coeffs1FSM_),
939  coeffs2FSM_(ptf.coeffs2FSM_),
940  szBox_(ptf.szBox_),
941  boxFactors2D_(ptf.boxFactors2D_),
942  boxFactors3D_(ptf.boxFactors3D_),
943  iNextToLastPlane_(ptf.iNextToLastPlane_),
944  box_(ptf.box_),
945  filterCoeffs_(ptf.filterCoeffs_),
946  filteredBox_(ptf.filteredBox_),
947  U0_(ptf.U0_)
948 {}
949 
950 
951 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
952 
954 {
955  if (updated())
956  {
957  return;
958  }
959 
960  if (curTimeIndex_ != db().time().timeIndex())
961  {
962  vectorField& U = *this;
963 
964  if (Pstream::master())
965  {
966  twoPointCorrs();
967  shiftRefill();
968  }
969 
970  Pstream::scatter(filteredBox_);
971 
972  mapFilteredBox(U);
973 
974  if (fsm_)
975  {
976  //(XC:Eq. 14)
977  for (direction dir = 0; dir < pTraits<vector>::nComponents; ++dir)
978  {
979  U.component(dir) =
980  U0_.component(dir)*coeffs1FSM_[dir]
981  + U.component(dir)*coeffs2FSM_[dir];
982  }
983 
984  U0_ = U;
985  }
986 
987  onePointCorrs(U);
988 
989  U += UMean_;
990 
991  if (correctFlowRate_)
992  {
993  // (KCX:Eq. 8)
994  U *= (flowRate_/gSum(U & -patch().Sf()));
995  }
996 
997  curTimeIndex_ = db().time().timeIndex();
998  }
999 
1000  fixedValueFvPatchVectorField::updateCoeffs();
1001 }
1002 
1003 
1006  Ostream& os
1007 ) const
1008 {
1010  os.writeEntry("n", n_);
1011  os.writeEntry("L", Lbak_);
1012 
1013  if (!interpR_)
1014  {
1015  R_.writeEntry("R", os);
1016  }
1017 
1018  if (!interpUMean_)
1019  {
1020  UMean_.writeEntry("UMean", os);
1021  }
1022 
1023  os.writeEntryIfDifferent<bool>("fsm", false, fsm_);
1024  os.writeEntryIfDifferent<bool>("Gaussian", true, Gaussian_);
1025  os.writeEntryIfDifferent<bool>("fixSeed", true, fixSeed_);
1026  os.writeEntryIfDifferent<bool>("continuous", false, continuous_);
1027  os.writeEntryIfDifferent<bool>("correctFlowRate", true, correctFlowRate_);
1028 
1029  if (!mapMethod_.empty())
1030  {
1032  (
1033  "mapMethod",
1034  "nearestCell",
1035  mapMethod_
1036  );
1037  }
1038 
1039  os.writeEntry("Ubulk", Ubulk_);
1040  os.writeEntry("C1", C1_);
1041  os.writeEntry("perturb", perturb_);
1042 
1043  if (fsm_)
1044  {
1045  os.writeEntry("C1FSM", C1FSM_);
1046  os.writeEntry("C2FSM", C2FSM_);
1047  }
1048 
1049  if (continuous_ && Pstream::master())
1050  {
1051  os.writeEntry("box", box_);
1052  }
1053 }
1054 
1055 
1058  const fvPatchFieldMapper& m
1059 )
1060 {
1062 }
1063 
1064 
1067  const fvPatchVectorField& ptf,
1068  const labelList& addr
1069 )
1070 {
1072 }
1073 
1074 
1075 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
1076 
1077 namespace Foam
1078 {
1080  (
1083  );
1084 }
1085 
1086 
1087 // ************************************************************************* //
Foam::fvPatchField< vector >
Foam::IOobject::NO_WRITE
Definition: IOobject.H:195
volFields.H
Foam::pointField
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:44
Foam::Tensor< scalar >
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
L
const vector L(dict.get< vector >("L"))
mathematicalConstants.H
p
volScalarField & p
Definition: createFieldRefs.H:8
Foam::word
A class for handling words, derived from Foam::string.
Definition: word.H:65
stdFoam::begin
constexpr auto begin(C &c) -> decltype(c.begin())
Return iterator to the beginning of the container c.
Definition: stdFoam.H:97
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::Zero
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
Foam::gAverage
Type gAverage(const FieldField< Field, Type > &f)
Definition: FieldFieldFunctions.C:604
Foam::inplaceRotateList
void inplaceRotateList(ListType< DataType > &list, label n)
Inplace reversal of a list using the Reversal Block Swapping algorithm.
Definition: ListOpsTemplates.C:1053
Foam::UPstream::master
static bool master(const label communicator=worldComm)
Am I the master process.
Definition: UPstream.H:457
Foam::Pstream::scatter
static void scatter(const List< commsStruct > &comms, T &Value, const int tag, const label comm)
Scatter data. Distribute without modification. Reverse of gather.
Definition: gatherScatter.C:150
Foam::pointToPointPlanarInterpolation
Interpolates between two sets of unstructured points using 2D Delaunay triangulation....
Definition: pointToPointPlanarInterpolation.H:53
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::exp
dimensionedScalar exp(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:261
Foam::turbulentDigitalFilterInletFvPatchVectorField::write
virtual void write(Ostream &) const
Write.
Definition: turbulentDigitalFilterInletFvPatchVectorField.C:1005
Foam::fixedValueFvPatchField
This boundary condition supplies a fixed value constraint, and is the base class for a number of othe...
Definition: fixedValueFvPatchField.H:80
turbulentDigitalFilterInletFvPatchVectorField.H
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:296
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::fixedValueFvPatchField::write
virtual void write(Ostream &) const
Write.
Definition: fixedValueFvPatchField.C:156
Foam::SymmTensor::xx
const Cmpt & xx() const
Definition: SymmTensorI.H:84
n
label n
Definition: TABSMDCalcMethod2.H:31
Foam::turbulentDigitalFilterInletFvPatchVectorField::autoMap
virtual void autoMap(const fvPatchFieldMapper &m)
Map (and resize as needed) from self given a mapping object.
Definition: turbulentDigitalFilterInletFvPatchVectorField.C:1057
R
#define R(A, B, C, D, E, F, K, M)
Foam::Field< symmTensor >
Foam::Info
messageStream Info
Information stream (stdout output on master, null elsewhere)
Foam::turbulentDigitalFilterInletFvPatchVectorField::updateCoeffs
virtual void updateCoeffs()
Update the coefficients associated with the patch field.
Definition: turbulentDigitalFilterInletFvPatchVectorField.C:953
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
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.
stdFoam::end
constexpr auto end(C &c) -> decltype(c.end())
Return iterator to the end of the container c.
Definition: stdFoam.H:121
Foam
Namespace for OpenFOAM.
Definition: atmBoundaryLayer.C:33
Foam::turbulentDigitalFilterInletFvPatchVectorField::turbulentDigitalFilterInletFvPatchVectorField
turbulentDigitalFilterInletFvPatchVectorField(const fvPatch &, const DimensionedField< vector, volMesh > &)
Construct from patch and internal field.
Definition: turbulentDigitalFilterInletFvPatchVectorField.C:646
Foam::turbulentDigitalFilterInletFvPatchVectorField
Definition: turbulentDigitalFilterInletFvPatchVectorField.H:349
Foam::vector
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:51
Foam::turbulentDigitalFilterInletFvPatchVectorField::rmap
virtual void rmap(const fvPatchVectorField &ptf, const labelList &addr)
Reverse map the given fvPatchField onto this fvPatchField.
Definition: turbulentDigitalFilterInletFvPatchVectorField.C:1066
Foam::exit
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:130
U
U
Definition: pEqn.H:72
FatalErrorInFunction
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
Foam::sqr
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Definition: dimensionedSymmTensor.C:51
Foam::constant::mathematical::pi
constexpr scalar pi(M_PI)
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::List
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
Definition: BitOps.H:63
Foam::sqrt
dimensionedScalar sqrt(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:144
Foam::pTraits
A traits class, which is primarily used for primitives.
Definition: pTraits.H:56
Us
Us
Definition: createFaFields.H:51
Foam::mag
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
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::sum
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
Definition: DimensionedFieldFunctions.C:327
Foam::direction
uint8_t direction
Definition: direction.H:52
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::name
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
Definition: exprTraits.C:59
Foam::fvPatchFieldMapper
Foam::fvPatchFieldMapper.
Definition: fvPatchFieldMapper.H:47
timeIndex
label timeIndex
Definition: getTimeIndex.H:30
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::Tuple2< label, label >
Foam::point
vector point
Point is a vector.
Definition: point.H:43
Foam::dictionary::getOrDefault
T getOrDefault(const word &keyword, const T &deflt, enum keyType::option matchOpt=keyType::REGEX) const
Definition: dictionaryTemplates.C:148
constant
constant condensation/saturation model.
WarningInFunction
#define WarningInFunction
Report a warning using Foam::Warning.
Definition: messageStream.H:328
Foam::tensor
Tensor< scalar > tensor
Tensor of scalars, i.e. Tensor<scalar>.
Definition: symmTensor.H:61
Foam::fvPatchField::autoMap
virtual void autoMap(const fvPatchFieldMapper &)
Map (and resize as needed) from self given a mapping object.
Definition: fvPatchField.C:248
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::IOobject::MUST_READ
Definition: IOobject.H:185