applications/solvers/compressible/rhoPorousSimpleFoam/pEqn.H

   1 if (pressureImplicitPorosity)
   2 {
   3     U = trTU()&UEqn().H();
   4 }
   5 else
   6 {
   7     U = trAU()*UEqn().H();
   8 }
   9
  10 UEqn.clear();
  11
  12 phi = fvc::interpolate(rho*U) & mesh.Sf();
  13 bool closedVolume = adjustPhi(phi, U, p);
  14
  15 for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
  16 {
  17     tmp<fvScalarMatrix> tpEqn;
  18
  19     if (pressureImplicitPorosity)
  20     {
  21         tpEqn = (fvm::laplacian(rho*trTU(), p) == fvc::div(phi));
  22     }
  23     else
  24     {
  25         tpEqn = (fvm::laplacian(rho*trAU(), p) == fvc::div(phi));
  26     }
  27
  28     tpEqn().setReference(pRefCell, pRefValue);
  29     // retain the residual from the first iteration
  30     if (nonOrth == 0)
  31     {
  32         eqnResidual = tpEqn().solve().initialResidual();
  33         maxResidual = max(eqnResidual, maxResidual);
  34     }
  35     else
  36     {
  37         tpEqn().solve();
  38     }
  39
  40     if (nonOrth == nNonOrthCorr)
  41     {
  42         phi -= tpEqn().flux();
  43     }
  44 }
  45
  46 #include "incompressible/continuityErrs.H"
  47
  48 // Explicitly relax pressure for momentum corrector
  49 p.relax();
  50
  51 if (pressureImplicitPorosity)
  52 {
  53     U -= trTU()&fvc::grad(p);
  54 }
  55 else
  56 {
  57     U -= trAU()*fvc::grad(p);
  58 }
  59
  60 U.correctBoundaryConditions();
  61
  62 bound(p, pMin);
  63
  64 // For closed-volume cases adjust the pressure and density levels
  65 // to obey overall mass continuity
  66 if (closedVolume)
  67 {
  68     p += (initialMass - fvc::domainIntegrate(psi*p))
  69         /fvc::domainIntegrate(psi);
  70 }
  71
  72 rho = thermo.rho();
  73 rho.relax();
  74 Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl;