{ volScalarField rUA = 1.0/UEqn.A(); surfaceScalarField psisf = fvc::interpolate(psis); surfaceScalarField rhof = fvc::interpolate(rho); // Needs to be outside of loop since p is changing, but psi and rho are not. surfaceScalarField rhoReff = rhof - psisf*fvc::interpolate(p); for (int corr = 0; corr < nCorr; corr++) { U = rUA*UEqn.H(); phi = rhoReff*(fvc::interpolate(U) & mesh.Sf()); surfaceScalarField phid("phid", psisf/rhoReff*phi); // Make fluxes relative within the MRF zone mrfZones.relativeFlux(psisf, phid); mrfZones.relativeFlux(rhoReff, phi); p.storePrevIter(); for (int nonOrth = 0; nonOrth <= nNonOrthCorr; nonOrth++) { fvScalarMatrix pEqn ( fvm::ddt(psis, p) + fvm::div(phid, p) + fvc::div(phi) - fvm::laplacian(rho*rUA, p) ); // Retain the residual from the first pressure solution eqnResidual = pEqn.solve().initialResidual(); if (corr == 0 && nonOrth == 0) { maxResidual = max(eqnResidual, maxResidual); } // Calculate the flux if (nonOrth == nNonOrthCorr) { phi += pEqn.flux(); } } # include "compressibleContinuityErrs.H" // Explicitly relax the pressure for momentum corrector p.relax(); U -= rUA*fvc::grad(p); U.correctBoundaryConditions(); } // Bound the pressure if (min(p) < pMin || max(p) > pMax) { p.max(pMin); p.min(pMax); p.correctBoundaryConditions(); } // Bound the velocity volScalarField magU = mag(U); if (max(magU) > UMax) { volScalarField Ulimiter = pos(magU - UMax)*UMax/(magU + smallU) + neg(magU - UMax); Ulimiter.max(scalar(0)); Ulimiter.min(scalar(1)); U *= Ulimiter; U.correctBoundaryConditions(); } }