DAIrkPimpleFoam.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2 
3  DAFoam : Discrete Adjoint with OpenFOAM
4  Version : v4
5 
6  This class is modified from OpenFOAM's source code
7  applications/solvers/incompressible/pimpleFoam
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 "DAIrkPimpleFoam.H"
31 
32 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
33 
34 namespace Foam
35 {
36 
37 defineTypeNameAndDebug(DAIrkPimpleFoam, 0);
38 addToRunTimeSelectionTable(DASolver, DAIrkPimpleFoam, dictionary);
39 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
40 
41 DAIrkPimpleFoam::DAIrkPimpleFoam(
42  char* argsAll,
43  PyObject* pyOptions)
44  : DASolver(argsAll, pyOptions),
45  // Radau23 coefficients and weights
46  D10(-2),
47  D11(3.0 / 2),
48  D12(1.0 / 2),
49  D20(2),
50  D21(-9.0 / 2),
51  D22(5.0 / 2),
52  w1(3.0 / 4),
53  w2(1.0 / 4),
54 
55  // SA-fv3 model coefficients
56  sigmaNut(0.66666),
57  kappa(0.41),
58  Cb1(0.1355),
59  Cb2(0.622),
60  Cw1(Cb1 / sqr(kappa) + (1.0 + Cb2) / sigmaNut),
61  Cw2(0.3),
62  Cw3(2.0),
63  Cv1(7.1),
64  Cv2(5.0)
65 {
66 }
67 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
68 
69 // Functions for SA-fv3 model
70 #include "mySAModel.H"
71 
72 // Some utilities, move to DAUtility later
73 #include "myUtilities.H"
74 
75 void DAIrkPimpleFoam::calcPriResIrkOrig(
76  volVectorField& U0, // oldTime U
77  volVectorField& U1, // 1st stage
78  volScalarField& p1,
79  surfaceScalarField& phi1,
80  volScalarField& nuTilda1,
81  volScalarField& nut1,
82  volVectorField& U2, // 2nd stage
83  volScalarField& p2,
84  surfaceScalarField& phi2,
85  volScalarField& nuTilda2,
86  volScalarField& nut2,
87  const volScalarField& nu,
88  const scalar& deltaT, // current dt
89  volVectorField& U1Res, // Residual for 1st stage
90  volScalarField& p1Res,
91  surfaceScalarField& phi1Res,
92  volVectorField& U2Res, // Residual for end stage
93  volScalarField& p2Res,
94  surfaceScalarField& phi2Res,
95  const scalar& relaxUEqn)
96 {
97  // Numerical settings
98  word divUScheme = "div(phi,U)";
99  word divGradUScheme = "div((nuEff*dev2(T(grad(U)))))";
100 
101  // Update boundaries
102  U0.correctBoundaryConditions(); // oldTime U
103  U1.correctBoundaryConditions(); // 1st stage
104  p1.correctBoundaryConditions();
105  nuTilda1.correctBoundaryConditions();
106  U2.correctBoundaryConditions(); // 2nd stage
107  p2.correctBoundaryConditions();
108  nuTilda2.correctBoundaryConditions();
109 
110  // --- 1st stage
111  {
112  // Get nuEff1
113  this->correctNut(nut1, nuTilda1, nu);
114  volScalarField nuEff1("nuEff1", nu + nut1);
115 
116  // Initialize U1Eqn w/o ddt term
117  fvVectorMatrix U1Eqn(
118  fvm::div(phi1, U1, divUScheme)
119  //+ turbulence->divDevReff(U)
120  - fvm::laplacian(nuEff1, U1)
121  - fvc::div(nuEff1 * dev2(T(fvc::grad(U1))), divGradUScheme));
122 
123  // Update U1Eqn with pseudo-spectral terms
124  forAll(U1, cellI)
125  {
126  scalar meshV = U1.mesh().V()[cellI];
127 
128  // Add D11 / halfDeltaT[i] * V() to diagonal
129  U1Eqn.diag()[cellI] += D11 / deltaT * meshV; // Use one seg for now: halfDeltaTList[segI]
130 
131  // Minus D10 / halfDeltaT[i] * T0 * V() to source term
132  U1Eqn.source()[cellI] -= D10 / deltaT * U0[cellI] * meshV;
133 
134  // Minus D12 / halfDeltaT[i] * T2 * V() to source term
135  U1Eqn.source()[cellI] -= D12 / deltaT * U2[cellI] * meshV;
136  }
137 
138  U1Res = (U1Eqn & U1) + fvc::grad(p1);
139 
140  // We use relaxation factor = 1.0, cannot skip the step below
141  U1Eqn.relax(relaxUEqn);
142 
143  volScalarField rAU1(1.0 / U1Eqn.A());
144  volVectorField HbyA1(constrainHbyA(rAU1 * U1Eqn.H(), U1, p1));
145  surfaceScalarField phiHbyA1(
146  "phiHbyA1",
147  fvc::flux(HbyA1));
148  tmp<volScalarField> rAtU1(rAU1);
149 
150  fvScalarMatrix p1Eqn(
151  fvm::laplacian(rAtU1(), p1) == fvc::div(phiHbyA1));
152 
153  p1Res = p1Eqn & p1;
154 
155  // Then do phiRes
156  phi1Res = phiHbyA1 - p1Eqn.flux() - phi1;
157  }
158 
159  // --- 2nd stage
160  {
161  // Get nuEff2
162  this->correctNut(nut2, nuTilda2, nu);
163  volScalarField nuEff2("nuEff2", nu + nut2);
164 
165  // Initialize U2Eqn w/o ddt term
166  fvVectorMatrix U2Eqn(
167  fvm::div(phi2, U2, divUScheme)
168  //+ turbulence->divDevReff(U)
169  - fvm::laplacian(nuEff2, U2)
170  - fvc::div(nuEff2 * dev2(T(fvc::grad(U2))), divGradUScheme));
171 
172  // Update U2Eqn with pseudo-spectral terms
173  forAll(U2, cellI)
174  {
175  scalar meshV = U2.mesh().V()[cellI];
176 
177  // Add D22 / halfDeltaT[i] * V() to diagonal
178  U2Eqn.diag()[cellI] += D22 / deltaT * meshV; // Use one seg for now: halfDeltaTList[segI]
179 
180  // Minus D20 / halfDeltaT[i] * T0 * V() to source term
181  U2Eqn.source()[cellI] -= D20 / deltaT * U0[cellI] * meshV;
182 
183  // Minus D21 / halfDeltaT[i] * T2 * V() to source term
184  U2Eqn.source()[cellI] -= D21 / deltaT * U1[cellI] * meshV;
185  }
186 
187  U2Res = (U2Eqn & U2) + fvc::grad(p2);
188 
189  // We use relaxation factor = 1.0, cannot skip the step below
190  U2Eqn.relax(relaxUEqn);
191 
192  volScalarField rAU2(1.0 / U2Eqn.A());
193  volVectorField HbyA2(constrainHbyA(rAU2 * U2Eqn.H(), U2, p2));
194  surfaceScalarField phiHbyA2(
195  "phiHbyA2",
196  fvc::flux(HbyA2));
197  tmp<volScalarField> rAtU2(rAU2);
198 
199  // Non-orthogonal pressure corrector loop
200  fvScalarMatrix p2Eqn(
201  fvm::laplacian(rAtU2(), p2) == fvc::div(phiHbyA2));
202 
203  p2Res = p2Eqn & p2;
204 
205  // Then do phiRes
206  phi2Res = phiHbyA2 - p2Eqn.flux() - phi2;
207  }
208 }
209 
210 void DAIrkPimpleFoam::calcPriSAResIrkOrig(
211  volScalarField& nuTilda0, // oldTime nuTilda
212  volVectorField& U1, // 1st stage
213  surfaceScalarField& phi1,
214  volScalarField& nuTilda1,
215  volVectorField& U2, // 2nd stage
216  surfaceScalarField& phi2,
217  volScalarField& nuTilda2,
218  volScalarField& y,
219  const volScalarField& nu,
220  const scalar& deltaT, // current dt
221  volScalarField& nuTilda1Res, // Residual for 1st stage
222  volScalarField& nuTilda2Res) // Residual for 2nd stage
223 {
224  // Numerical settings
225  word divNuTildaScheme = "div(phi,nuTilda)";
226 
227  // Update boundaries
228  nuTilda0.correctBoundaryConditions();
229  U1.correctBoundaryConditions();
230  nuTilda1.correctBoundaryConditions();
231  U2.correctBoundaryConditions();
232  nuTilda2.correctBoundaryConditions();
233 
234  // --- 1st stage
235  {
236  // Get chi1 and fv11
237  volScalarField chi1("chi1", chi(nuTilda1, nu));
238  volScalarField fv11("fv11", fv1(chi1));
239 
240  // Get Stilda1
241  volScalarField Stilda1(
242  "Stilda1",
243  fv3(chi1, fv11) * ::sqrt(2.0) * mag(skew(fvc::grad(U1))) + fv2(chi1, fv11) * nuTilda1 / sqr(kappa * y));
244 
245  // Construct nuTilda1Eqn w/o ddt term
246  fvScalarMatrix nuTilda1Eqn(
247  fvm::div(phi1, nuTilda1, divNuTildaScheme)
248  - fvm::laplacian(DnuTildaEff(nuTilda1, nu), nuTilda1)
249  - Cb2 / sigmaNut * magSqr(fvc::grad(nuTilda1))
250  == Cb1 * Stilda1 * nuTilda1 // If field inversion, beta should be multiplied here
251  - fvm::Sp(Cw1 * fw(Stilda1, nuTilda1, y) * nuTilda1 / sqr(y), nuTilda1));
252 
253  // Update nuTilda1Eqn with pseudo-spectral terms
254  forAll(nuTilda1, cellI)
255  {
256  scalar meshV = nuTilda1.mesh().V()[cellI];
257 
258  // Add D11 / halfDeltaT[i] * V() to diagonal
259  nuTilda1Eqn.diag()[cellI] += D11 / deltaT * meshV;
260 
261  // Minus D10 / halfDeltaT[i] * T0 * V() to source term
262  nuTilda1Eqn.source()[cellI] -= D10 / deltaT * nuTilda0[cellI] * meshV;
263 
264  // Minus D12 / halfDeltaT[i] * T2 * V() to source term
265  nuTilda1Eqn.source()[cellI] -= D12 / deltaT * nuTilda2[cellI] * meshV;
266  }
267 
268  nuTilda1Res = nuTilda1Eqn & nuTilda1;
269  }
270 
271  // --- 2nd stage
272  {
273  // Get chi2 and fv12
274  volScalarField chi2("chi2", chi(nuTilda2, nu));
275  volScalarField fv12("fv12", fv1(chi2));
276 
277  // Get Stilda2
278  volScalarField Stilda2(
279  "Stilda2",
280  fv3(chi2, fv12) * ::sqrt(2.0) * mag(skew(fvc::grad(U2))) + fv2(chi2, fv12) * nuTilda2 / sqr(kappa * y));
281 
282  // Construct nuTilda2Eqn w/o ddt term
283  fvScalarMatrix nuTilda2Eqn(
284  fvm::div(phi2, nuTilda2, divNuTildaScheme)
285  - fvm::laplacian(DnuTildaEff(nuTilda2, nu), nuTilda2)
286  - Cb2 / sigmaNut * magSqr(fvc::grad(nuTilda2))
287  == Cb1 * Stilda2 * nuTilda2 // If field inversion, beta should be multiplied here
288  - fvm::Sp(Cw1 * fw(Stilda2, nuTilda2, y) * nuTilda2 / sqr(y), nuTilda2));
289 
290  // Update nuTilda2Eqn with pseudo-spectral terms
291  forAll(nuTilda2, cellI)
292  {
293  scalar meshV = nuTilda2.mesh().V()[cellI];
294 
295  // Add D22 / halfDeltaT[i] * V() to diagonal
296  nuTilda2Eqn.diag()[cellI] += D22 / deltaT * meshV;
297 
298  // Minus D20 / halfDeltaT[i] * T0 * V() to source term
299  nuTilda2Eqn.source()[cellI] -= D20 / deltaT * nuTilda0[cellI] * meshV;
300 
301  // Minus D21 / halfDeltaT[i] * T2 * V() to source term
302  nuTilda2Eqn.source()[cellI] -= D21 / deltaT * nuTilda1[cellI] * meshV;
303  }
304 
305  nuTilda2Res = nuTilda2Eqn & nuTilda2;
306  }
307 }
308 
309 void DAIrkPimpleFoam::initSolver()
310 {
311  /*
312  Description:
313  Initialize variables for DASolver
314  */
315  daOptionPtr_.reset(new DAOption(meshPtr_(), pyOptions_));
316 }
317 
318 label DAIrkPimpleFoam::solvePrimal()
319 {
320  /*
321  Description:
322  Call the primal solver to get converged state variables
323 
324  Output:
325  state variable vector
326  */
327 
328  Foam::argList& args = argsPtr_();
329 #include "createTime.H"
330 #include "createMesh.H"
331 #include "initContinuityErrs.H"
332 #include "createControl.H"
333 #include "createFieldsIrkPimple.H"
334 #include "CourantNo.H"
335 
336  // Turbulence disabled
337  //turbulence->validate();
338 
339  // Get nu, nut, nuTilda, y
340  volScalarField& nu = const_cast<volScalarField&>(mesh.thisDb().lookupObject<volScalarField>("nu"));
341  volScalarField& nut = const_cast<volScalarField&>(mesh.thisDb().lookupObject<volScalarField>("nut"));
342  volScalarField& nuTilda = const_cast<volScalarField&>(mesh.thisDb().lookupObject<volScalarField>("nuTilda"));
343  volScalarField& y = const_cast<volScalarField&>(mesh.thisDb().lookupObject<volScalarField>("yWall"));
344 
345  // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
346 
347  Info << "\nStarting time loop\n"
348  << endl;
349 
350  // get IRKDict settings, default to Radau23 for now
351  IOdictionary IRKDict(
352  IOobject(
353  "IRKDict",
354  mesh.time().system(),
355  mesh,
356  IOobject::READ_IF_PRESENT,
357  IOobject::NO_WRITE));
358 
359  scalar relaxU = 1.0;
360  if (IRKDict.found("relaxU"))
361  {
362  if (IRKDict.getScalar("relaxU") > 0)
363  {
364  relaxU = IRKDict.getScalar("relaxU");
365  }
366  }
367 
368  scalar relaxP = 1.0;
369  if (IRKDict.found("relaxP"))
370  {
371  if (IRKDict.getScalar("relaxP") > 0)
372  {
373  relaxP = IRKDict.getScalar("relaxP");
374  }
375  }
376 
377  scalar relaxPhi = 1.0;
378  if (IRKDict.found("relaxPhi"))
379  {
380  if (IRKDict.getScalar("relaxPhi") > 0)
381  {
382  relaxPhi = IRKDict.getScalar("relaxPhi");
383  }
384  }
385 
386  scalar relaxNuTilda = 1.0;
387  if (IRKDict.found("relaxNuTilda"))
388  {
389  if (IRKDict.getScalar("relaxNuTilda") > 0)
390  {
391  relaxNuTilda = IRKDict.getScalar("relaxNuTilda");
392  }
393  }
394 
395  scalar relaxStage1 = 0.8;
396  if (IRKDict.found("relaxStage1"))
397  {
398  if (IRKDict.getScalar("relaxStage1") > 0)
399  {
400  relaxStage1 = IRKDict.getScalar("relaxStage1");
401  }
402  }
403 
404  scalar relaxStage2 = 0.8;
405  if (IRKDict.found("relaxStage2"))
406  {
407  if (IRKDict.getScalar("relaxStage2") > 0)
408  {
409  relaxStage2 = IRKDict.getScalar("relaxStage2");
410  }
411  }
412 
413  scalar relaxUEqn = 1.0;
414  scalar relaxNuTildaEqn = 1.0;
415  if (IRKDict.found("relaxNuTildaEqn"))
416  {
417  if (IRKDict.getScalar("relaxNuTildaEqn") > 0)
418  {
419  relaxNuTildaEqn = IRKDict.getScalar("relaxNuTildaEqn");
420  }
421  }
422 
423  label maxSweep = 10;
424  if (IRKDict.found("maxSweep"))
425  {
426  if (IRKDict.getLabel("maxSweep") > 0)
427  {
428  maxSweep = IRKDict.getLabel("maxSweep");
429  }
430  }
431 
432  // Duplicate state variables for stages
433  volVectorField U1("U1", U);
434  volVectorField U2("U2", U);
435  volScalarField p1("p1", p);
436  volScalarField p2("p2", p);
437  surfaceScalarField phi1("phi1", phi);
438  surfaceScalarField phi2("phi2", phi);
439  // SA turbulence model
440  volScalarField nuTilda1("nuTilda1", nuTilda);
441  volScalarField nuTilda2("nuTilda2", nuTilda);
442  volScalarField nut1("nut1", nut);
443  volScalarField nut2("nut2", nut);
444 
445  // Settings for stage pressure
446  mesh.setFluxRequired(p1.name());
447  mesh.setFluxRequired(p2.name());
448 
449  // Initialize primal residuals
450  volVectorField U1Res(
451  IOobject(
452  "U1Res",
453  runTime.timeName(),
454  mesh,
455  IOobject::NO_READ,
456  IOobject::NO_WRITE),
457  mesh,
458  dimensionedVector("U1Res", dimensionSet(0, 1, -2, 0, 0, 0, 0), vector::zero),
459  zeroGradientFvPatchField<vector>::typeName);
460  volVectorField U2Res("U2Res", U1Res);
461 
462  volScalarField p1Res(
463  IOobject(
464  "p1Res",
465  runTime.timeName(),
466  mesh,
467  IOobject::NO_READ,
468  IOobject::NO_WRITE),
469  mesh,
470  dimensionedScalar("p1Res", dimensionSet(0, 0, -1, 0, 0, 0, 0), 0.0),
471  zeroGradientFvPatchField<scalar>::typeName);
472  volScalarField p2Res("p2Res", p1Res);
473 
474  surfaceScalarField phi1Res(
475  IOobject(
476  "phi1Res",
477  runTime.timeName(),
478  mesh,
479  IOobject::NO_READ,
480  IOobject::NO_WRITE),
481  phi * 0.0);
482  surfaceScalarField phi2Res("phi2Res", phi1Res);
483 
484  volScalarField nuTilda1Res(
485  IOobject(
486  "nuTilda1Res",
487  runTime.timeName(),
488  mesh,
489  IOobject::NO_READ,
490  IOobject::NO_WRITE),
491  mesh,
492  dimensionedScalar("nuTilda1Res", dimensionSet(0, 2, -2, 0, 0, 0, 0), 0.0),
493  zeroGradientFvPatchField<scalar>::typeName);
494  volScalarField nuTilda2Res("nuTilda2Res", nuTilda1Res);
495 
496  // Initialize oldTime() for under-relaxation
497  U1.oldTime() = U1;
498  U2.oldTime() = U2;
499  p1.oldTime() = p1;
500  p2.oldTime() = p2;
501  phi1.oldTime() = phi1;
502  phi2.oldTime() = phi2;
503  nuTilda1.oldTime() = nuTilda1;
504  nuTilda2.oldTime() = nuTilda2;
505 
506  // Numerical settings
507  word divUScheme = "div(phi,U)";
508  word divGradUScheme = "div((nuEff*dev2(T(grad(U)))))";
509  word divNuTildaScheme = "div(phi,nuTilda)";
510 
511  const fvSolution& myFvSolution = mesh.thisDb().lookupObject<fvSolution>("fvSolution");
512  dictionary solverDictU = myFvSolution.subDict("solvers").subDict("U");
513  dictionary solverDictP = myFvSolution.subDict("solvers").subDict("p");
514  dictionary solverDictNuTilda = myFvSolution.subDict("solvers").subDict("nuTilda");
515 
516  while (runTime.run())
517  {
518 
519 #include "CourantNo.H"
520 
521  ++runTime;
522 
523  Info << "Time = " << runTime.timeName() << nl << endl;
524 
525  scalar deltaT = runTime.deltaTValue();
526 
527  // --- GS sweeps for IRK-PIMPLE
528  label sweepIndex = 0;
529  while (sweepIndex < maxSweep)
530  {
531  Info << "Block GS sweep = " << sweepIndex + 1 << endl;
532 
533  // --- 1st stage
534  {
535 #include "U1EqnIrkPimple.H"
536 
537  while (pimple.correct())
538  {
539 #include "p1EqnIrkPimple.H"
540  }
541 
542  // --- Correct turbulence, using our own SAFv3
543 #include "nuTilda1EqnIrkPimple.H"
544  }
545 
546  // --- 2nd stage
547  {
548 #include "U2EqnIrkPimple.H"
549 
550  while (pimple.correct())
551  {
552 #include "p2EqnIrkPimple.H"
553  }
554 
555  // --- Correct turbulence, using our own SAFv3
556 #include "nuTilda2EqnIrkPimple.H"
557  }
558 
559  this->calcPriResIrkOrig(U, U1, p1, phi1, nuTilda1, nut1, U2, p2, phi2, nuTilda2, nut2, nu, deltaT, U1Res, p1Res, phi1Res, U2Res, p2Res, phi2Res, relaxUEqn);
560  this->calcPriSAResIrkOrig(nuTilda, U1, phi1, nuTilda1, U2, phi2, nuTilda2, y, nu, deltaT, nuTilda1Res, nuTilda2Res);
561 
562  Info << "L2 norm of U1Res: " << this->L2norm(U1Res.primitiveField()) << endl;
563  Info << "L2 norm of U2Res: " << this->L2norm(U2Res.primitiveField()) << endl;
564  Info << "L2 norm of p1Res: " << this->L2norm(p1Res.primitiveField()) << endl;
565  Info << "L2 norm of p2Res: " << this->L2norm(p2Res.primitiveField()) << endl;
566  Info << "L2 norm of phi1Res: " << this->L2norm(phi1Res, phi1.mesh().magSf()) << endl;
567  Info << "L2 norm of phi2Res: " << this->L2norm(phi2Res, phi2.mesh().magSf()) << endl;
568  Info << "L2 norm of nuTilda1Res: " << this->L2norm(nuTilda1Res.primitiveField()) << endl;
569  Info << "L2 norm of nuTilda2Res: " << this->L2norm(nuTilda2Res.primitiveField()) << endl;
570 
571  sweepIndex++;
572  }
573 
574  // Update new step values before write-to-disk
575  U = U2;
576  U.correctBoundaryConditions();
577  p = p2;
578  p.correctBoundaryConditions();
579  phi = phi2;
580  nuTilda = nuTilda2;
581  nuTilda.correctBoundaryConditions();
582  nut = nut2;
583  nut.correctBoundaryConditions();
584 
585  // Write to disk
586  runTime.write(); // This writes U, p, phi, nuTilda, nut
587  // Also write internal stages to disk (Radau23)
588  U1.write();
589  p1.write();
590  phi1.write();
591  nuTilda1.write();
592  nut1.write();
593 
594  // Use old step as initial guess for the next step
595  U1 = U;
596  U1.correctBoundaryConditions();
597  p1 = p;
598  p1.correctBoundaryConditions();
599  phi1 = phi;
600  nuTilda1 = nuTilda;
601  nuTilda1.correctBoundaryConditions();
602  nut1 = nut;
603  nut1.correctBoundaryConditions();
604 
605  runTime.printExecutionTime(Info);
606  }
607 
608  Info << "End\n"
609  << endl;
610 
611  return 0;
612 }
613 
614 } // End namespace Foam
615 
616 // ************************************************************************* //
HbyA1
volVectorField HbyA1(constrainHbyA(rAU1 *U1Eqn.H(), U1, p1))
Foam::DAIrkPimpleFoam::Cb1
scalar Cb1
Definition: DAIrkPimpleFoam.H:70
Foam::DAIrkPimpleFoam::calcPriResIrkOrig
void calcPriResIrkOrig(volVectorField &U0, volVectorField &U1, volScalarField &p1, surfaceScalarField &phi1, volScalarField &nuTilda1, volScalarField &nut1, volVectorField &U2, volScalarField &p2, surfaceScalarField &phi2, volScalarField &nuTilda2, volScalarField &nut2, const volScalarField &nu, const scalar &deltaT, volVectorField &U1Res, volScalarField &p1Res, surfaceScalarField &phi1Res, volVectorField &U2Res, volScalarField &p2Res, surfaceScalarField &phi2Res, const scalar &relaxUEqn)
Definition: DAIrkPimpleFoam.C:75
U
volVectorField & U
Definition: createRefsPimpleDyM.H:7
Stilda2
volScalarField Stilda2("Stilda2", fv3(chi2, fv12) *::sqrt(2.0) *mag(skew(fvc::grad(U2)))+fv2(chi2, fv12) *nuTilda2/sqr(kappa *y))
phi2
phi2
Definition: p2EqnIrkPimple.H:19
Foam::DAIrkPimpleFoam::fw
tmp< volScalarField > fw(volScalarField &Stilda, volScalarField &nuTilda, volScalarField &y)
Definition: mySAModel.H:41
Foam::DAIrkPimpleFoam::D12
scalar D12
Definition: DAIrkPimpleFoam.H:60
fv11
volScalarField fv11("fv11", fv1(chi1))
Foam::DAIrkPimpleFoam::fv2
tmp< volScalarField > fv2(volScalarField &chi, volScalarField &fv1)
Definition: mySAModel.H:21
pimple
pimpleControlDF & pimple
Definition: createRefsPimpleDyM.H:5
phiHbyA2
surfaceScalarField phiHbyA2("phiHbyA2", fvc::flux(HbyA2))
Foam::DAIrkPimpleFoam::solvePrimal
virtual label solvePrimal()
solve the primal equations
Definition: DAIrkPimpleFoam.C:318
U2
U2
Definition: p2EqnIrkPimple.H:26
nuTilda1EqnIrkPimple.H
Foam::DAIrkPimpleFoam::L2norm
scalar L2norm(const List< scalar > &v)
Definition: myUtilities.H:11
p1EqnIrkPimple.H
nuTilda2EqnIrkPimple.H
Foam::DASolver
Definition: DASolver.H:54
nuEff2
volScalarField nuEff2("nuEff2", nu+nut2)
Foam::DAOption
Definition: DAOption.H:29
nuTilda2Eqn
fvScalarMatrix nuTilda2Eqn(fvm::div(phi2, nuTilda2, divNuTildaScheme) - fvm::laplacian(DnuTildaEff(nuTilda2, nu), nuTilda2) - Cb2/sigmaNut *magSqr(fvc::grad(nuTilda2))==Cb1 *Stilda2 *nuTilda2 - fvm::Sp(Cw1 *fw(Stilda2, nuTilda2, y) *nuTilda2/sqr(y), nuTilda2))
Foam::DAIrkPimpleFoam::Cb2
scalar Cb2
Definition: DAIrkPimpleFoam.H:71
Foam::addToRunTimeSelectionTable
addToRunTimeSelectionTable(DAFunction, DAFunctionForce, dictionary)
createFieldsIrkPimple.H
nuTilda1Eqn
fvScalarMatrix nuTilda1Eqn(fvm::div(phi1, nuTilda1, divNuTildaScheme) - fvm::laplacian(DnuTildaEff(nuTilda1, nu), nuTilda1) - Cb2/sigmaNut *magSqr(fvc::grad(nuTilda1))==Cb1 *Stilda1 *nuTilda1 - fvm::Sp(Cw1 *fw(Stilda1, nuTilda1, y) *nuTilda1/sqr(y), nuTilda1))
Foam::DAIrkPimpleFoam::chi
tmp< volScalarField > chi(volScalarField &nuTilda, const volScalarField &nu)
Definition: mySAModel.H:7
U1Eqn
fvVectorMatrix U1Eqn(fvm::div(phi1, U1, divUScheme) - fvm::laplacian(nuEff1, U1) - fvc::div(nuEff1 *dev2(T(fvc::grad(U1))), divGradUScheme))
Foam::DAIrkPimpleFoam::D10
scalar D10
Definition: DAIrkPimpleFoam.H:58
Foam::DASolver::daOptionPtr_
autoPtr< DAOption > daOptionPtr_
DAOption pointer.
Definition: DASolver.H:81
Foam::DASolver::pyOptions_
PyObject * pyOptions_
all options in DAFoam
Definition: DASolver.H:69
Foam::DAIrkPimpleFoam::fv1
tmp< volScalarField > fv1(volScalarField &chi)
Definition: mySAModel.H:14
Foam::DAIrkPimpleFoam::correctNut
void correctNut(volScalarField &nut, volScalarField &nuTilda, const volScalarField &nu)
Definition: mySAModel.H:70
Foam::DAIrkPimpleFoam::DnuTildaEff
tmp< volScalarField > DnuTildaEff(volScalarField &nuTilda, const volScalarField &nu)
Definition: mySAModel.H:62
mySAModel.H
U2Eqn
fvVectorMatrix U2Eqn(fvm::div(phi2, U2, divUScheme) - fvm::laplacian(nuEff2, U2) - fvc::div(nuEff2 *dev2(T(fvc::grad(U2))), divGradUScheme))
rAtU2
tmp< volScalarField > rAtU2(rAU2)
phiHbyA1
surfaceScalarField phiHbyA1("phiHbyA1", fvc::flux(HbyA1))
U1EqnIrkPimple.H
chi1
volScalarField chi1("chi1", chi(nuTilda1, nu))
Foam::DAIrkPimpleFoam::D22
scalar D22
Definition: DAIrkPimpleFoam.H:63
mesh
fvMesh & mesh
Definition: createRefsHeatTransfer.H:4
nuTilda2
nuTilda2
Definition: nuTilda2EqnIrkPimple.H:38
phi1
phi1
Definition: p1EqnIrkPimple.H:19
myUtilities.H
Foam::DAIrkPimpleFoam::Cw1
scalar Cw1
Definition: DAIrkPimpleFoam.H:72
Stilda1
volScalarField Stilda1("Stilda1", fv3(chi1, fv11) *::sqrt(2.0) *mag(skew(fvc::grad(U1)))+fv2(chi1, fv11) *nuTilda1/sqr(kappa *y))
p2EqnIrkPimple.H
rAU1
volScalarField rAU1(1.0/U1Eqn.A())
p
volScalarField & p
Definition: createRefsPimpleDyM.H:6
fv12
volScalarField fv12("fv12", fv1(chi2))
nuEff1
volScalarField nuEff1("nuEff1", nu+nut1)
HbyA2
volVectorField HbyA2(constrainHbyA(rAU2 *U2Eqn.H(), U2, p2))
Foam::defineTypeNameAndDebug
defineTypeNameAndDebug(DAFunction, 0)
Foam::DAIrkPimpleFoam::D11
scalar D11
Definition: DAIrkPimpleFoam.H:59
phi
surfaceScalarField & phi
Definition: createRefsPimpleDyM.H:8
Foam
Definition: checkGeometry.C:32
Foam::DAIrkPimpleFoam::kappa
scalar kappa
Definition: DAIrkPimpleFoam.H:69
Foam::DAIrkPimpleFoam::initSolver
virtual void initSolver()
initialize fields and variables
Definition: DAIrkPimpleFoam.C:309
T
volScalarField & T
Definition: createRefsHeatTransfer.H:5
Foam::DAIrkPimpleFoam::sigmaNut
scalar sigmaNut
Definition: DAIrkPimpleFoam.H:68
Foam::DAIrkPimpleFoam::calcPriSAResIrkOrig
void calcPriSAResIrkOrig(volScalarField &nuTilda0, volVectorField &U1, surfaceScalarField &phi1, volScalarField &nuTilda1, volVectorField &U2, surfaceScalarField &phi2, volScalarField &nuTilda2, volScalarField &y, const volScalarField &nu, const scalar &deltaT, volScalarField &nuTilda1Res, volScalarField &nuTilda2Res)
Definition: DAIrkPimpleFoam.C:210
nuTilda1
nuTilda1
Definition: nuTilda1EqnIrkPimple.H:38
Foam::DAIrkPimpleFoam::D20
scalar D20
Definition: DAIrkPimpleFoam.H:61
p1
p1
Definition: p1EqnIrkPimple.H:24
forAll
forAll(nuTilda1, cellI)
Definition: nuTilda1EqnIrkPimple.H:19
Foam::DAIrkPimpleFoam::DAIrkPimpleFoam
DAIrkPimpleFoam(char *argsAll, PyObject *pyOptions)
Definition: DAIrkPimpleFoam.C:41
U2EqnIrkPimple.H
rAtU1
tmp< volScalarField > rAtU1(rAU1)
Foam::DASolver::meshPtr_
autoPtr< fvMesh > meshPtr_
fvMesh pointer
Definition: DASolver.H:78
Foam::DAIrkPimpleFoam::D21
scalar D21
Definition: DAIrkPimpleFoam.H:62
Foam::DASolver::argsPtr_
autoPtr< argList > argsPtr_
args pointer
Definition: DASolver.H:72
chi2
volScalarField chi2("chi2", chi(nuTilda2, nu))
rAU2
volScalarField rAU2(1.0/U2Eqn.A())
p1Eqn
fvScalarMatrix p1Eqn(fvm::laplacian(rAtU1(), p1)==fvc::div(phiHbyA1))
args
Foam::argList & args
Definition: setRootCasePython.H:42
U1
U1
Definition: p1EqnIrkPimple.H:26
runTime
Time & runTime
Definition: createRefsHeatTransfer.H:1
p2Eqn
fvScalarMatrix p2Eqn(fvm::laplacian(rAtU2(), p2)==fvc::div(phiHbyA2))
DAIrkPimpleFoam.H
Foam::DAIrkPimpleFoam::fv3
tmp< volScalarField > fv3(volScalarField &chi, volScalarField &fv1)
Definition: mySAModel.H:28
p2
p2
Definition: p2EqnIrkPimple.H:24

Generated by   doxygen 1.8.17