applications/solvers/compressible/sonicDyMFoam/pEqn.H

   1 {
   2     U = UEqn.H()/UEqn.A();
   3
   4 #   include "limitU.H"
   5
   6     for (int nonOrth = 0; nonOrth <= nNonOrthCorr; nonOrth++)
   7     {
   8         // Calculate phi for boundary conditions
   9         phi = rhof*
  10         (
  11             (fvc::interpolate(U) & mesh.Sf())
  12           - fvc::meshPhi(rho, U)
  13         );
  14
  15         surfaceScalarField phid2 = rhoReff/rhof*phi;
  16
  17         surfaceScalarField phid("phid", psisf/rhof*phi);
  18
  19         // Store pressure for under-relaxation
  20         p.storePrevIter();
  21
  22         volScalarField divPhid
  23         (
  24             "divPhid",
  25             fvc::div(phid)
  26         );
  27
  28         fvScalarMatrix pEqn
  29         (
  30             fvm::ddt(psis, p)
  31           + fvm::div(phid, p)
  32             // Convective flux relaxation terms
  33           + fvm::SuSp(-divPhid, p)
  34           + divPhid*p
  35           + fvc::div(phid2)
  36           - fvm::laplacian(rho*rUA, p)
  37         );
  38
  39         if
  40         (
  41 //              oCorr == nOuterCorr - 1
  42             corr == nCorr - 1
  43          && nonOrth == nNonOrthCorr
  44         )
  45         {
  46             pEqn.solve
  47             (
  48                 mesh.solutionDict().solver(p.name() + "Final")
  49             );
  50         }
  51         else
  52         {
  53             pEqn.solve(mesh.solutionDict().solver(p.name()));
  54         }
  55
  56         // Calculate the flux
  57         if (nonOrth == nNonOrthCorr)
  58         {
  59             phi = phid2 + pEqn.flux();
  60         }
  61
  62
  63         // Bound the pressure
  64         if (min(p) < pMin || max(p) > pMax)
  65         {
  66             p.max(pMin);
  67             p.min(pMax);
  68             p.correctBoundaryConditions();
  69         }
  70
  71         // Relax the pressure
  72         p.relax();
  73     }
  74
  75 #   include "compressibleContinuityErrs.H"
  76
  77     U -= fvc::grad(p)/UEqn.A();
  78     U.correctBoundaryConditions();
  79
  80 #   include "limitU.H"
  81 }