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