DAkEpsilon.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2 
3  DAFoam : Discrete Adjoint with OpenFOAM
4  Version : v4
5 
6  This file is modified from OpenFOAM's source code
7  src/TurbulenceModels/turbulenceModels/RAS/kEpsilon/kEpsilon.H
8 
9  OpenFOAM: The Open Source CFD Toolbox
10 
11  Copyright (C): 2011-2016 OpenFOAM Foundation
12 
13  OpenFOAM License:
14 
15  OpenFOAM is free software: you can redistribute it and/or modify it
16  under the terms of the GNU General Public License as published by
17  the Free Software Foundation, either version 3 of the License, or
18  (at your option) any later version.
19 
20  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
21  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
22  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
23  for more details.
24 
25  You should have received a copy of the GNU General Public License
26  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
27 
28 \*---------------------------------------------------------------------------*/
29 
30 #include "DAkEpsilon.H"
31 
32 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
33 
34 namespace Foam
35 {
36 
37 defineTypeNameAndDebug(DAkEpsilon, 0);
38 addToRunTimeSelectionTable(DATurbulenceModel, DAkEpsilon, dictionary);
39 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
40 
42  const word modelType,
43  const fvMesh& mesh,
44  const DAOption& daOption)
45  : DATurbulenceModel(modelType, mesh, daOption),
46  // KE parameters
47  Cmu_(dimensioned<scalar>::lookupOrAddToDict(
48  "Cmu",
49  this->coeffDict_,
50  0.09)),
51  C1_(dimensioned<scalar>::lookupOrAddToDict(
52  "C1",
53  this->coeffDict_,
54  1.44)),
55  C2_(dimensioned<scalar>::lookupOrAddToDict(
56  "C2",
57  this->coeffDict_,
58  1.92)),
59  C3_(dimensioned<scalar>::lookupOrAddToDict(
60  "C3",
61  this->coeffDict_,
62  0)),
63  sigmak_(dimensioned<scalar>::lookupOrAddToDict(
64  "sigmak",
65  this->coeffDict_,
66  1.0)),
67  sigmaEps_(dimensioned<scalar>::lookupOrAddToDict(
68  "sigmaEps",
69  this->coeffDict_,
70  1.3)),
71  // Augmented variables
72  epsilon_(const_cast<volScalarField&>(
73  mesh_.thisDb().lookupObject<volScalarField>("epsilon"))),
74  epsilonRes_(
75  IOobject(
76  "epsilonRes",
77  mesh.time().timeName(),
78  mesh,
79  IOobject::NO_READ,
80  IOobject::NO_WRITE),
81  mesh,
82  dimensionedScalar("epsilonRes", dimensionSet(0, 0, 0, 0, 0, 0, 0), 0.0),
83  zeroGradientFvPatchField<scalar>::typeName),
84  k_(const_cast<volScalarField&>(
85  mesh_.thisDb().lookupObject<volScalarField>("k"))),
86  kRes_(
87  IOobject(
88  "kRes",
89  mesh.time().timeName(),
90  mesh,
91  IOobject::NO_READ,
92  IOobject::NO_WRITE),
93  mesh,
94  dimensionedScalar("kRes", dimensionSet(0, 0, 0, 0, 0, 0, 0), 0.0),
95  zeroGradientFvPatchField<scalar>::typeName),
96  GPtr_(nullptr),
97  betaFIK_(
98  IOobject(
99  "betaFIK",
100  mesh.time().timeName(),
101  mesh,
102  IOobject::READ_IF_PRESENT,
103  IOobject::AUTO_WRITE),
104  mesh,
105  dimensionedScalar("betaFIK", dimensionSet(0, 0, 0, 0, 0, 0, 0), 1.0),
106  "zeroGradient"),
107  betaFIEpsilon_(
108  IOobject(
109  "betaFIEpsilon",
110  mesh.time().timeName(),
111  mesh,
112  IOobject::READ_IF_PRESENT,
113  IOobject::AUTO_WRITE),
114  mesh,
115  dimensionedScalar("betaFIEpsilon", dimensionSet(0, 0, 0, 0, 0, 0, 0), 1.0),
116  "zeroGradient")
117 {
118 
119  // we need to reset the nuTildaRes's dimension based on the model type
120  if (turbModelType_ == "incompressible")
121  {
122  epsilonRes_.dimensions().reset(dimensionSet(0, 2, -4, 0, 0, 0, 0));
123  kRes_.dimensions().reset(dimensionSet(0, 2, -3, 0, 0, 0, 0));
124  }
125 
126  if (turbModelType_ == "compressible")
127  {
128  epsilonRes_.dimensions().reset(dimensionSet(1, -1, -4, 0, 0, 0, 0));
129  kRes_.dimensions().reset(dimensionSet(1, -1, -3, 0, 0, 0, 0));
130  }
131 
132  // calculate the size of epsilonWallFunction faces
133  label nWallFaces = 0;
134  forAll(epsilon_.boundaryField(), patchI)
135  {
136  if (epsilon_.boundaryField()[patchI].type() == "epsilonWallFunction"
137  and epsilon_.boundaryField()[patchI].size() > 0)
138  {
139  forAll(epsilon_.boundaryField()[patchI], faceI)
140  {
141  nWallFaces++;
142  //Info<<"patchI: "<<patchI<<" faceI: "<<faceI<<endl;
143  }
144  }
145  }
146 
147  // initialize epsilonNearWall
148  epsilonNearWall_.setSize(nWallFaces);
149 
150  // initialize the G field
151  tmp<volTensorField> tgradU = fvc::grad(U_);
152  GPtr_.reset(new volScalarField::Internal(
153  "kEpsilon:G",
154  nut_.v() * (dev(twoSymm(tgradU().v())) && tgradU().v())));
155 }
156 
157 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
158 
159 // KE member functions
160 tmp<fvScalarMatrix> DAkEpsilon::kSource() const
161 {
162  return tmp<fvScalarMatrix>(
163  new fvScalarMatrix(
164  k_,
165  dimVolume * this->rhoDimensions() * k_.dimensions()
166  / dimTime));
167 }
168 
169 tmp<fvScalarMatrix> DAkEpsilon::epsilonSource() const
170 {
171  return tmp<fvScalarMatrix>(
172  new fvScalarMatrix(
173  epsilon_,
174  dimVolume * this->rhoDimensions() * epsilon_.dimensions()
175  / dimTime));
176 }
177 
178 tmp<volScalarField> DAkEpsilon::DkEff() const
179 {
180  return tmp<volScalarField>(
181  new volScalarField(
182  "DkEff",
183  (nut_ / sigmak_ + this->nu())));
184 }
185 
186 tmp<volScalarField> DAkEpsilon::DepsilonEff() const
187 {
188  return tmp<volScalarField>(
189  new volScalarField(
190  "DepsilonEff",
191  (nut_ / sigmaEps_ + this->nu())));
192 }
193 
194 // Augmented functions
195 void DAkEpsilon::correctModelStates(wordList& modelStates) const
196 {
197  /*
198  Description:
199  Update the name in modelStates based on the selected physical model at runtime
200 
201  Example:
202  In DAStateInfo, if the modelStates reads:
203 
204  modelStates = {"nut"}
205 
206  then for the SA model, calling correctModelStates(modelStates) will give:
207 
208  modelStates={"nuTilda"}
209 
210  while calling correctModelStates(modelStates) for the SST model will give
211 
212  modelStates={"k","epsilon"}
213 
214  We don't udpate the names for the radiation model because users are
215  supposed to set modelStates={"G"}
216  */
217 
218  // For SST model, we need to replace nut with k, epsilon
219 
220  forAll(modelStates, idxI)
221  {
222  word stateName = modelStates[idxI];
223  if (stateName == "nut")
224  {
225  modelStates[idxI] = "epsilon";
226  modelStates.append("k");
227  }
228  }
229 }
230 
232 {
233  /*
234  Description:
235  Update nut based on other turbulence variables and update the BCs
236  Also update alphat if is present
237  */
238 
239  nut_ = Cmu_ * sqr(k_) / epsilon_;
240 
241  nut_.correctBoundaryConditions(); // nutkWallFunction: update wall face nut based on k
242 
243  // this is basically BasicTurbulenceModel::correctNut();
244  this->correctAlphat();
245 
246  return;
247 }
248 
250 {
251  /*
252  Description:
253  Update turbulence variable boundary values
254  */
255 
256  // correct the BCs for the perturbed fields
257  // kqWallFunction is a zero-gradient BC
258  k_.correctBoundaryConditions();
259 }
260 
262 {
263  /*
264  Description:
265  this is a special treatment for epsilon BC because we cant directly call epsilon.
266  correctBoundaryConditions() because it will modify the internal epsilon and G that
267  are right next to walls. This will mess up adjoint Jacobians
268  To solve this issue,
269  1. we store the near wall epsilon before calling epsilon.correctBoundaryConditions()
270  2. call epsilon.correctBoundaryConditions()
271  3. Assign the stored near wall epsilon back
272  4. Apply a zeroGradient BC for epsilon at the wall patches
273  *********** NOTE *************
274  this treatment will obviously downgrade the accuracy of adjoint derivative since it is
275  not 100% consistent with what is used for the flow solver; however, our observation
276  shows that the impact is not very large.
277  */
278 
279  // save the perturbed epsilon at the wall
280  this->saveEpsilonNearWall();
281  // correct epsilon boundary conditions, this includes updating wall face and near wall epsilon values,
282  // updating the inter-proc BCs
283  epsilon_.correctBoundaryConditions();
284  // reset the corrected epsilon near wall cell to its perturbed value
285  this->setEpsilonNearWall();
286 }
287 
289 {
290  /*
291  Description:
292  Save the current near wall epsilon values to epsilonNearWall_
293  */
294  label counterI = 0;
295  forAll(epsilon_.boundaryField(), patchI)
296  {
297  if (epsilon_.boundaryField()[patchI].type() == "epsilonWallFunction"
298  and epsilon_.boundaryField()[patchI].size() > 0)
299  {
300  const UList<label>& faceCells = mesh_.boundaryMesh()[patchI].faceCells();
301  forAll(faceCells, faceI)
302  {
303  //Info<<"Near Wall cellI: "<<faceCells[faceI]<<endl;
304  epsilonNearWall_[counterI] = epsilon_[faceCells[faceI]];
305  counterI++;
306  }
307  }
308  }
309  return;
310 }
311 
313 {
314  /*
315  Description:
316  Set the current near wall epsilon values from epsilonNearWall_
317  Here we also apply a zeroGradient BC to the wall faces
318  */
319  label counterI = 0;
320  forAll(epsilon_.boundaryField(), patchI)
321  {
322  if (epsilon_.boundaryField()[patchI].type() == "epsilonWallFunction"
323  && epsilon_.boundaryField()[patchI].size() > 0)
324  {
325  const UList<label>& faceCells = mesh_.boundaryMesh()[patchI].faceCells();
326  forAll(faceCells, faceI)
327  {
328  epsilon_[faceCells[faceI]] = epsilonNearWall_[counterI];
329  // zeroGradient BC
330  epsilon_.boundaryFieldRef()[patchI][faceI] = epsilon_[faceCells[faceI]];
331  counterI++;
332  }
333  }
334  }
335  return;
336 }
337 
339 {
340  /*
341  Description:
342  Update nut based on nuTilda. Note: we need to update nut and its BC since we
343  may have perturbed other turbulence vars that affect the nut values
344  */
345 
346  this->correctNut();
347 }
348 
349 void DAkEpsilon::correctStateResidualModelCon(List<List<word>>& stateCon) const
350 {
351  /*
352  Description:
353  Update the original variable connectivity for the adjoint state
354  residuals in stateCon. Basically, we modify/add state variables based on the
355  original model variables defined in stateCon.
356 
357  Input:
358 
359  stateResCon: the connectivity levels for a state residual, defined in Foam::DAJacCon
360 
361  Example:
362  If stateCon reads:
363  stateCon=
364  {
365  {"U", "p", "nut"},
366  {"p"}
367  }
368 
369  For the SA turbulence model, calling this function for will get a new stateCon
370  stateCon=
371  {
372  {"U", "p", "nuTilda"},
373  {"p"}
374  }
375 
376  For the SST turbulence model, calling this function will give
377  stateCon=
378  {
379  {"U", "p", "k", "omega"},
380  {"p", "U"}
381  }
382  ***NOTE***: we add a extra level of U connectivity because nut is
383  related to grad(U), k, and omega in SST!
384  */
385 
386  forAll(stateCon, idxI)
387  {
388  forAll(stateCon[idxI], idxJ)
389  {
390  word conStateName = stateCon[idxI][idxJ];
391  if (conStateName == "nut")
392  {
393  stateCon[idxI][idxJ] = "epsilon"; // replace nut with epsilon
394  stateCon[idxI].append("k"); // also add k for that level
395  }
396  }
397  }
398 }
399 
400 void DAkEpsilon::addModelResidualCon(HashTable<List<List<word>>>& allCon) const
401 {
402  /*
403  Description:
404  Add the connectivity levels for all physical model residuals to allCon
405 
406  Input:
407  allCon: the connectivity levels for all state residual, defined in DAJacCon
408 
409  Example:
410  If stateCon reads:
411  allCon=
412  {
413  "URes":
414  {
415  {"U", "p", "nut"},
416  {"p"}
417  }
418  }
419 
420  For the SA turbulence model, calling this function for will get a new stateCon,
421  something like this:
422  allCon=
423  {
424  "URes":
425  {
426  {"U", "p", "nuTilda"},
427  {"p"}
428  },
429  "nuTildaRes":
430  {
431  {"U", "phi", "nuTilda"},
432  {"U"}
433  }
434  }
435 
436  */
437 
438  word pName;
439 
440  if (mesh_.thisDb().foundObject<volScalarField>("p"))
441  {
442  pName = "p";
443  }
444  else if (mesh_.thisDb().foundObject<volScalarField>("p_rgh"))
445  {
446  pName = "p_rgh";
447  }
448  else
449  {
450  FatalErrorIn(
451  "Neither p nor p_rgh was found in mesh.thisDb()!"
452  "addModelResidualCon failed to setup turbulence residuals!")
453  << exit(FatalError);
454  }
455 
456  // NOTE: for compressible flow, it depends on rho so we need to add T and p
457  if (turbModelType_ == "incompressible")
458  {
459  allCon.set(
460  "epsilonRes",
461  {
462  {"U", "epsilon", "k", "phi"}, // lv0
463  {"U", "epsilon", "k"}, // lv1
464  {"U", "epsilon", "k"} // lv2
465  });
466  allCon.set(
467  "kRes",
468  {
469  {"U", "epsilon", "k", "phi"}, // lv0
470  {"U", "epsilon", "k"}, // lv1
471  {"U", "epsilon", "k"} // lv2
472  });
473  }
474  else if (turbModelType_ == "compressible")
475  {
476  allCon.set(
477  "epsilonRes",
478  {
479  {"U", "T", pName, "epsilon", "k", "phi"}, // lv0
480  {"U", "T", pName, "epsilon", "k"}, // lv1
481  {"U", "T", pName, "epsilon", "k"} // lv2
482  });
483  allCon.set(
484  "kRes",
485  {
486  {"U", "T", pName, "epsilon", "k", "phi"}, // lv0
487  {"U", "T", pName, "epsilon", "k"}, // lv1
488  {"U", "T", pName, "epsilon", "k"} // lv2
489  });
490  }
491 }
492 
493 void DAkEpsilon::correct(label printToScreen)
494 {
495  /*
496  Descroption:
497  Solve the residual equations and update the state. This function will be called
498  by the DASolver. It is needed because we want to control the output frequency
499  of the residual convergence every 100 steps. If using the correct from turbulence
500  it will output residual every step which will be too much of information.
501  */
502 
503  // We set the flag solveTurbState_ to 1 such that in the calcResiduals function
504  // we will solve and update nuTilda
505  solveTurbState_ = 1;
506  dictionary dummyOptions;
507  dummyOptions.set("printToScreen", printToScreen);
508  this->calcResiduals(dummyOptions);
509  // after it, we reset solveTurbState_ = 0 such that calcResiduals will not
510  // update nuTilda when calling from the adjoint class, i.e., solveAdjoint from DASolver.
511  solveTurbState_ = 0;
512 }
513 
514 void DAkEpsilon::calcResiduals(const dictionary& options)
515 {
516  /*
517  Descroption:
518  If solveTurbState_ == 1, this function solve and update k and epsilon, and
519  is the same as calling turbulence.correct(). If solveTurbState_ == 0,
520  this function compute residuals for turbulence variables, e.g., nuTildaRes_
521 
522  Input:
523  options.isPC: 1 means computing residuals for preconditioner matrix.
524  This essentially use the first order scheme for div(phi,nuTilda)
525 
526  p_, U_, phi_, etc: State variables in OpenFOAM
527 
528  Output:
529  kRes_/epsilonRes_: If solveTurbState_ == 0, update the residual field variable
530 
531  k_/epsilon_: If solveTurbState_ == 1, update them
532  */
533 
534  // Copy and modify based on the "correct" function
535 
536  word divKScheme = "div(phi,k)";
537  word divEpsilonScheme = "div(phi,epsilon)";
538 
539  label isPC = 0;
540 
541  label printToScreen = options.lookupOrDefault("printToScreen", 0);
542 
543  if (!solveTurbState_)
544  {
545  isPC = options.getLabel("isPC");
546 
547  if (isPC)
548  {
549  divKScheme = "div(pc)";
550  divEpsilonScheme = "div(pc)";
551  }
552  }
553 
554  volScalarField rho = this->rho();
555 
556  // Note: for compressible flow, the "this->phi()" function divides phi by fvc:interpolate(rho),
557  // while for the incompresssible "this->phi()" returns phi only
558  // see src/TurbulenceModels/compressible/compressibleTurbulenceModel.C line 62 to 73
559  volScalarField::Internal divU(
560  fvc::div(fvc::absolute(phi_ / fvc::interpolate(rho), U_))().v());
561 
562  tmp<volTensorField> tgradU = fvc::grad(U_);
563  volScalarField::Internal& G = const_cast<volScalarField::Internal&>(GPtr_());
564  G = nut_.v() * (dev(twoSymm(tgradU().v())) && tgradU().v());
565  tgradU.clear();
566 
567  if (solveTurbState_)
568  {
569  // Update epsilon and G at the wall
570  epsilon_.boundaryFieldRef().updateCoeffs();
571  }
572  else
573  {
574  // special treatment for epsilon BC
576  }
577 
578  // Dissipation equation
579  tmp<fvScalarMatrix> epsEqn(
580  fvm::ddt(phase_, rho, epsilon_)
581  + fvm::div(phaseRhoPhi_, epsilon_, divEpsilonScheme)
582  - fvm::laplacian(phase_ * rho * DepsilonEff(), epsilon_)
583  == C1_ * phase_() * rho() * G * epsilon_() / k_() * betaFIEpsilon_()
584  - fvm::SuSp((scalar(2.0 / 3.0) * C1_ - C3_) * phase_() * rho() * divU, epsilon_)
585  - fvm::Sp(C2_ * phase_() * rho() * epsilon_() / k_(), epsilon_)
586  + epsilonSource());
587 
588  epsEqn.ref().relax();
589 
590  epsEqn.ref().boundaryManipulate(epsilon_.boundaryFieldRef());
591 
592  if (solveTurbState_)
593  {
594  // get the solver performance info such as initial
595  // and final residuals
596  SolverPerformance<scalar> solverEpsilon = solve(epsEqn);
597  DAUtility::primalResidualControl(solverEpsilon, printToScreen, "epsilon", daGlobalVar_.primalMaxRes);
598 
599  DAUtility::boundVar(allOptions_, epsilon_, printToScreen);
600  }
601  else
602  {
603  // reset the corrected omega near wall cell to its perturbed value
604  this->setEpsilonNearWall();
605 
606  // calculate residuals
607  epsilonRes_ = epsEqn() & epsilon_;
608  // need to normalize residuals
609  normalizeResiduals(epsilonRes);
610  }
611 
612  // Turbulent kinetic energy equation
613  tmp<fvScalarMatrix> kEqn(
614  fvm::ddt(phase_, rho, k_)
615  + fvm::div(phaseRhoPhi_, k_, divKScheme)
616  - fvm::laplacian(phase_ * rho * DkEff(), k_)
617  == phase_() * rho() * G * betaFIK_()
618  - fvm::SuSp((2.0 / 3.0) * phase_() * rho() * divU, k_)
619  - fvm::Sp(phase_() * rho() * epsilon_() / k_(), k_)
620  + kSource());
621 
622  kEqn.ref().relax();
623 
624  if (solveTurbState_)
625  {
626 
627  // get the solver performance info such as initial
628  // and final residuals
629  SolverPerformance<scalar> solverK = solve(kEqn);
630  DAUtility::primalResidualControl(solverK, printToScreen, "k", daGlobalVar_.primalMaxRes);
631 
632  DAUtility::boundVar(allOptions_, k_, printToScreen);
633 
634  this->correctNut();
635  }
636  else
637  {
638  // calculate residuals
639  kRes_ = kEqn() & k_;
640  // need to normalize residuals
641  normalizeResiduals(kRes);
642  }
643 
644  return;
645 }
646 
648  const word varName,
649  scalarField& diag,
650  scalarField& upper,
651  scalarField& lower)
652 {
653  /*
654  Description:
655  return the diag(), upper(), and lower() scalarFields from the turbulence model's fvMatrix
656  this will be use to compute the preconditioner matrix
657  */
658 
659  if (varName != "k" && varName != "epsilon")
660  {
661  FatalErrorIn(
662  "varName not valid. It has to be k or epsilon")
663  << exit(FatalError);
664  }
665 
666  volScalarField rho = this->rho();
667 
668  // Note: for compressible flow, the "this->phi()" function divides phi by fvc:interpolate(rho),
669  // while for the incompresssible "this->phi()" returns phi only
670  // see src/TurbulenceModels/compressible/compressibleTurbulenceModel.C line 62 to 73
671  volScalarField::Internal divU(
672  fvc::div(fvc::absolute(phi_ / fvc::interpolate(rho), U_))().v());
673 
674  tmp<volTensorField> tgradU = fvc::grad(U_);
675  volScalarField::Internal& G = const_cast<volScalarField::Internal&>(GPtr_());
676  G = nut_.v() * (dev(twoSymm(tgradU().v())) && tgradU().v());
677  tgradU.clear();
678 
679  // special treatment for epsilon BC
681 
682  if (varName == "epsilon")
683  {
684  // Dissipation equation
685  fvScalarMatrix epsEqn(
686  fvm::ddt(phase_, rho, epsilon_)
687  + fvm::div(phaseRhoPhi_, epsilon_, "div(pc)")
688  - fvm::laplacian(phase_ * rho * DepsilonEff(), epsilon_)
689  == C1_ * phase_() * rho() * G * epsilon_() / k_() * betaFIEpsilon_()
690  - fvm::SuSp((scalar(2.0 / 3.0) * C1_ - C3_) * phase_() * rho() * divU, epsilon_)
691  - fvm::Sp(C2_ * phase_() * rho() * epsilon_() / k_(), epsilon_)
692  + epsilonSource());
693 
694  epsEqn.relax();
695 
696  // reset the corrected omega near wall cell to its perturbed value
697  this->setEpsilonNearWall();
698 
699  diag = epsEqn.D();
700  upper = epsEqn.upper();
701  lower = epsEqn.lower();
702  }
703  else if (varName == "k")
704  {
705  fvScalarMatrix kEqn(
706  fvm::ddt(phase_, rho, k_)
707  + fvm::div(phaseRhoPhi_, k_, "div(pc)")
708  - fvm::laplacian(phase_ * rho * DkEff(), k_)
709  == phase_() * rho() * G * betaFIK_()
710  - fvm::SuSp((2.0 / 3.0) * phase_() * rho() * divU, k_)
711  - fvm::Sp(phase_() * rho() * epsilon_() / k_(), k_)
712  + kSource());
713 
714  kEqn.relax();
715 
716  diag = kEqn.D();
717  upper = kEqn.upper();
718  lower = kEqn.lower();
719  }
720 }
721 
722 void DAkEpsilon::getTurbProdOverDestruct(volScalarField& PoD) const
723 {
724  /*
725  Description:
726  Return the value of the production over destruction term from the turbulence model
727  */
728  tmp<volTensorField> tgradU = fvc::grad(U_);
729  volScalarField::Internal& G = const_cast<volScalarField::Internal&>(GPtr_());
730  G = nut_.v() * (dev(twoSymm(tgradU().v())) && tgradU().v());
731 
732  volScalarField rho = this->rho();
733 
734  volScalarField::Internal P = phase_() * rho() * G;
735  volScalarField::Internal D = phase_() * rho() * epsilon_();
736 
737  forAll(P, cellI)
738  {
739  PoD[cellI] = P[cellI] / (D[cellI] + P[cellI] + 1e-16);
740  }
741 }
742 
743 void DAkEpsilon::getTurbConvOverProd(volScalarField& CoP) const
744 {
745  /*
746  Description:
747  Return the value of the convective over production term from the turbulence model
748  */
749 
750  tmp<volTensorField> tgradU = fvc::grad(U_);
751  volScalarField::Internal& G = const_cast<volScalarField::Internal&>(GPtr_());
752  G = nut_.v() * (dev(twoSymm(tgradU().v())) && tgradU().v());
753 
754  volScalarField rho = this->rho();
755 
756  volScalarField::Internal P = phase_() * rho() * G;
757  volScalarField C = fvc::div(phaseRhoPhi_, k_);
758 
759  forAll(P, cellI)
760  {
761  CoP[cellI] = C[cellI] / (P[cellI] + C[cellI] + 1e-16);
762  }
763 }
764 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
765 
766 } // End namespace Foam
767 
768 // ************************************************************************* //
Foam::DAkEpsilon::DkEff
tmp< volScalarField > DkEff() const
Definition: DAkEpsilon.C:178
Foam::DAkEpsilon::epsilonNearWall_
scalarList epsilonNearWall_
Definition: DAkEpsilon.H:104
Foam::DAkEpsilon::getTurbProdOverDestruct
virtual void getTurbProdOverDestruct(volScalarField &PoD) const
return the value of the destruction term from the turbulence model
Definition: DAkEpsilon.C:722
Foam::DAkEpsilon::Cmu_
dimensionedScalar Cmu_
Definition: DAkEpsilon.H:55
Foam::DATurbulenceModel::phi_
surfaceScalarField & phi_
face flux
Definition: DATurbulenceModel.H:83
Foam::DAkEpsilon::addModelResidualCon
virtual void addModelResidualCon(HashTable< List< List< word >>> &allCon) const
add the model residual connectivity to stateCon
Definition: DAkEpsilon.C:400
Foam::DAkEpsilon::epsilon_
volScalarField & epsilon_
Definition: DAkEpsilon.H:86
solve
nuTilda1Eqn solve(solverDictNuTilda)
D
volVectorField & D
Definition: createRefsSolidDisplacement.H:8
Foam::DAkEpsilon::getTurbConvOverProd
virtual void getTurbConvOverProd(volScalarField &CoP) const
return the value of the convective over production term from the turbulence model
Definition: DAkEpsilon.C:743
Foam::DAOption
Definition: DAOption.H:29
Foam::addToRunTimeSelectionTable
addToRunTimeSelectionTable(DAFunction, DAFunctionForce, dictionary)
Foam::DAkEpsilon::correctBoundaryConditions
virtual void correctBoundaryConditions()
update turbulence variable boundary values
Definition: DAkEpsilon.C:249
Foam::DAkEpsilon::betaFIK_
volScalarField betaFIK_
beta field for field inversion
Definition: DAkEpsilon.H:97
Foam::DAkEpsilon::calcResiduals
virtual void calcResiduals(const dictionary &options)
compute the turbulence residuals
Definition: DAkEpsilon.C:514
normalizeResiduals
#define normalizeResiduals(resName)
Definition: DAMacroFunctions.H:28
Foam::DATurbulenceModel::phaseRhoPhi_
surfaceScalarField & phaseRhoPhi_
phase*phi*density field
Definition: DATurbulenceModel.H:89
Foam::DATurbulenceModel::allOptions_
const dictionary & allOptions_
all DAFoam options
Definition: DATurbulenceModel.H:71
Foam::DATurbulenceModel::rho
tmp< volScalarField > rho() const
get the density field
Definition: DATurbulenceModel.C:294
mesh
fvMesh & mesh
Definition: createRefsHeatTransfer.H:4
Foam::DAUtility::boundVar
static void boundVar(const dictionary &allOptions, volScalarField &var, const label printToScreen)
bound a volScalar variable based on parametes defined in DAOption::allOptions_
Definition: DAUtility.C:475
Foam::DAkEpsilon::getFvMatrixFields
virtual void getFvMatrixFields(const word varName, scalarField &diag, scalarField &upper, scalarField &lower)
return the diag(), upper(), and lower() scalarFields from the turbulence model's fvMatrix
Definition: DAkEpsilon.C:647
Foam::DAUtility::primalResidualControl
static void primalResidualControl(const SolverPerformance< scalar > &solverP, const label printToScreen, const word varName, scalar &primalMaxRes)
control when to print the residual and also compute the maxInitRes
Definition: DAUtility.C:734
Foam::DATurbulenceModel::daGlobalVar_
DAGlobalVar & daGlobalVar_
global variable
Definition: DATurbulenceModel.H:74
Foam::DAkEpsilon::correctModelStates
virtual void correctModelStates(wordList &modelStates) const
update the turbulence state for DAStateInfo::regStates_
Definition: DAkEpsilon.C:195
Foam::DATurbulenceModel::U_
volVectorField & U_
velocity
Definition: DATurbulenceModel.H:80
Foam::DAkEpsilon::correct
virtual void correct(label printToScreen)
solve the residual equations and update the state
Definition: DAkEpsilon.C:493
Foam::DATurbulenceModel::nut_
volScalarField & nut_
turbulence viscosity
Definition: DATurbulenceModel.H:77
Foam::DAkEpsilon::sigmaEps_
dimensionedScalar sigmaEps_
Definition: DAkEpsilon.H:60
Foam::DAGlobalVar::primalMaxRes
scalar primalMaxRes
the maximal residual for primal this var will be used in multiple classes
Definition: DAGlobalVar.H:62
Foam::defineTypeNameAndDebug
defineTypeNameAndDebug(DAFunction, 0)
Foam::DAkEpsilon::betaFIEpsilon_
volScalarField betaFIEpsilon_
Definition: DAkEpsilon.H:98
Foam::DAkEpsilon::epsilonRes_
volScalarField epsilonRes_
Definition: DAkEpsilon.H:87
Foam
Definition: checkGeometry.C:32
Foam::DAkEpsilon::C1_
dimensionedScalar C1_
Definition: DAkEpsilon.H:56
Foam::DAkEpsilon::GPtr_
autoPtr< volScalarField::Internal > GPtr_
Definition: DAkEpsilon.H:94
Foam::DATurbulenceModel::rhoDimensions
dimensionSet rhoDimensions() const
return the dimension of rho
Definition: DATurbulenceModel.C:321
Foam::DATurbulenceModel::nu
tmp< volScalarField > nu() const
get the nu field
Definition: DATurbulenceModel.C:267
Foam::DAkEpsilon::kSource
tmp< fvScalarMatrix > kSource() const
Definition: DAkEpsilon.C:160
Foam::DAkEpsilon::correctEpsilonBoundaryConditions
void correctEpsilonBoundaryConditions()
specially treatment to correct epsilon BC
Definition: DAkEpsilon.C:261
Foam::DAkEpsilon::setEpsilonNearWall
void setEpsilonNearWall()
set epsilonNearWall_ to near wall epsilon values
Definition: DAkEpsilon.C:312
forAll
forAll(nuTilda1, cellI)
Definition: nuTilda1EqnIrkPimple.H:19
Foam::DAkEpsilon::C3_
dimensionedScalar C3_
Definition: DAkEpsilon.H:58
Foam::DATurbulenceModel::correctAlphat
void correctAlphat()
update alphat
Definition: DATurbulenceModel.C:203
Foam::DAkEpsilon::correctNut
virtual void correctNut()
update nut based on other turbulence variables and update the BCs
Definition: DAkEpsilon.C:231
Foam::DAkEpsilon::k_
volScalarField & k_
Definition: DAkEpsilon.H:88
Foam::DAkEpsilon::epsilonSource
tmp< fvScalarMatrix > epsilonSource() const
Definition: DAkEpsilon.C:169
Foam::DAkEpsilon::kRes_
volScalarField kRes_
Definition: DAkEpsilon.H:89
Foam::DATurbulenceModel
Definition: DATurbulenceModel.H:52
Foam::DAkEpsilon::correctStateResidualModelCon
virtual void correctStateResidualModelCon(List< List< word >> &stateCon) const
update the original variable connectivity for the adjoint state residuals in stateCon
Definition: DAkEpsilon.C:349
Foam::DATurbulenceModel::mesh_
const fvMesh & mesh_
fvMesh
Definition: DATurbulenceModel.H:65
Foam::DAkEpsilon::DAkEpsilon
DAkEpsilon(const word modelType, const fvMesh &mesh, const DAOption &daOption)
Definition: DAkEpsilon.C:41
Foam::DAkEpsilon::saveEpsilonNearWall
void saveEpsilonNearWall()
save near wall epsilon values to epsilonNearWall_
Definition: DAkEpsilon.C:288
Foam::DAkEpsilon::updateIntermediateVariables
virtual void updateIntermediateVariables()
update any intermediate variables that are dependent on state variables and are used in calcResiduals
Definition: DAkEpsilon.C:338
Foam::DAkEpsilon::solveTurbState_
label solveTurbState_
whether to solve for turb states
Definition: DAkEpsilon.H:107
DAkEpsilon.H
Foam::DAkEpsilon::sigmak_
dimensionedScalar sigmak_
Definition: DAkEpsilon.H:59
Foam::DAkEpsilon::DepsilonEff
tmp< volScalarField > DepsilonEff() const
Definition: DAkEpsilon.C:186
Foam::DAkEpsilon::C2_
dimensionedScalar C2_
Definition: DAkEpsilon.H:57
Foam::DATurbulenceModel::turbModelType_
word turbModelType_
whether the turbulence model is incompressible or compressible
Definition: DATurbulenceModel.H:95
Foam::DATurbulenceModel::phase_
volScalarField phase_
phase field
Definition: DATurbulenceModel.H:86