Steady-state PISO solver for compressible flow

I developed steady-state PISO solver for compressible flow.
On developing the solver, I used the algorithm published in [1].

In PISO method, there are two corrector steps.

First corrector step is identical to the corrector step of SIMPLE.
And PISO has an additional corrector step.

I named the new solver 'myRhoPisoSFoam' by applying the following changes to myRhoSimpleFoam:

1. Adds new pEqn2.H.

pEqn2.H

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

    volVectorField diffHbyA("diffHbyA", rAU*(UEqn().H() - H1));
    H1 = UEqn().H();

    UEqn.clear();

    //bool closedVolume = false;

    if (simple.transonic())
    {
        surfaceScalarField phid
        (
            "phid",
            fvc::interpolate(psi)
           *(fvc::interpolate(U) & mesh.Sf())
        );
        surfaceScalarField phidDiffHbyA
        (
            "phidDiffHbyA",
            fvc::interpolate(psi)
           *(fvc::interpolate(diffHbyA) & mesh.Sf())
        );

        fvOptions.makeRelative(fvc::interpolate(psi), phid);

        while (simple.correctNonOrthogonal())
        {
            fvScalarMatrix pPrimeEqn2
            (
                fvm::div(phid, pPrime)
              + fvm::div(phidDiffHbyA, pPrime)
              - fvm::laplacian(rho*rAU, pPrime)
              ==
                fvOptions(psi, pPrime, rho.name())
            );

            // Relax the pressure equation to ensure diagonal-dominance
            pPrimeEqn2.relax();

            fvOptions.constrain(pPrimeEqn2);

            //pPrimeEqn2.setReference(pRefCell, pRefValue);

            pPrimeEqn2.solve();

            if (simple.finalNonOrthogonalIter())
            {
                phi == pPrimeEqn2.flux();
            }
        }
    }
    else
    {
        surfaceScalarField phiUstar
        (
            "phiUstar",
            fvc::interpolate(rho*U) & mesh.Sf()
        );
        surfaceScalarField phiDiffHbyA
        (
            "phiDiffHbyA",
            fvc::interpolate(rho*diffHbyA) & mesh.Sf()
        );

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

        //closedVolume = adjustPhi(phiUstar, U, p);

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

            //pPrimeEqn2.setReference(pRefCell, pRefValue);

            fvOptions.constrain(pPrimeEqn2);

            pPrimeEqn2.solve();

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


    #include "incompressible/continuityErrs.H"

    // Explicitly relax pressure for momentum corrector
    pPrime.relax();
    
    p += pPrime;
    U += diffHbyA - rAU*fvc::grad(pPrime);
    p.correctBoundaryConditions();
    U.correctBoundaryConditions();
    fvOptions.correct(U);

    // For closed-volume cases adjust the pressure and density levels
    // to obey overall mass continuity
/*
    if (closedVolume)
    {
        Info<< "closedVolume" << endl;
        p += (initialMass - fvc::domainIntegrate(psi*p))
            /fvc::domainIntegrate(psi);
    }
*/
    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;
    Info<< "pPrime max/min : "
        << max(pPrime).value() << " "
        << min(pPrime).value() << endl;

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

myRhoPisoSFoam.C:

        // Pressure-velocity SIMPLE corrector
        {
            #include "UEqn.H"
            #include "EEqn.H"

            #include "pEqn.H"
            #include "pEqn2.H"
        }

References

1. N. W. Bressloff, A parallel pressure implicit plitting of operators algorithm applied to flows at all speeds, International Journal for Numerical Methods in Fluids 2001; 36:497–518(DOI: 10.1002/fld.140).