Modified steady-state SIMPLER solver for compressible flow

I modified the SIMPLER solver by applying the idea of modified SIMPLE solver.

Step 5-7 in the SIMPLER solver developed in the previous article are changed as follows:

5'. Solve pressure equation for p'

volScalarField rAU(1.0/UEqn().A());
//volVectorField HbyA("HbyA", U);
//HbyA = rAU*UEqn().H();

UEqn.clear();

bool closedVolume = false;

surfaceScalarField phiUstar
(
    "phiUstar",
    fvc::interpolate(rho*U) & mesh.Sf()
);

fvOptions.makeRelative(fvc::interpolate(rho), phiUstar);

closedVolume = adjustPhi(phiUstar, U, p);

while (simple.correctNonOrthogonal())
{
    fvScalarMatrix pPrimeEqn
    (
        fvc::div(phiUstar)
      - fvm::laplacian(rho*rAU, pPrime)
      ==
        fvOptions(psi, pPrime, rho.name())
    );

    pPrimeEqn.setReference(pRefCell, pRefValue);

    fvOptions.constrain(pPrimeEqn);

    pPrimeEqn.solve();

    if (simple.finalNonOrthogonalIter())
    {
        phi = phiUstar + pPrimeEqn.flux();
    }
}

#include "incompressible/continuityErrs.H"

// Explicitly relax pressure for momentum corrector
pPrime.relax();

6'. Calculate U = (Ustar - grad(pPrime))/A and don't correct p.

Info<< "Calculating U" << endl;
//p += pPrime
U -= rAU*fvc::grad(pPrime);
U.correctBoundaryConditions();
fvOptions.correct(U);

// For closed-volume cases adjust the pressure and density levels
// to obey overall mass continuity
/*
if (closedVolume)
{
    p += (initialMass - fvc::domainIntegrate(psi*p))
        /fvc::domainIntegrate(psi);
}
*/

// Recalculate density from the relaxed pressure
/*
rho = thermo.rho();
rho = max(rho, rhoMin);
rho = min(rho, rhoMax);

if (!simple.transonic())
{
    rho.relax();
}

Info<< "rho max/min : " << max(rho).value()

        << " " << min(rho).value() << endl;
*/

7'. Set pPrime = 0

pPrime = min
         (
             max
             (
                 pPrime, dimensionedScalar("0", dimPressure, scalar(0))
             ),
             dimensionedScalar("0", dimPressure, scalar(0))
         );