Development of steady-state SIMPLER solver for compressible flow

I developed a steady-state SIMPLER solver for compressible flow and named it 'myRhoSimplerFoam'.

Steps of SIMPLER algorithm and differences from rhoSimpleFoam is as follows.
The unchanged part is indicated in green, commented-out part in blue, and added part in red.

1. Compute momentum coefficients

tmp<fvVectorMatrix> UEqn
(
    fvm::div(phi, U)
  + turbulence->divDevRhoReff(U)
 ==
    fvOptions(rho, U)
);

UEqn().relax();

2. Calculate Uhat = H(U0)/A

/* Do not solve here
fvOptions.constrain(UEqn());

solve(UEqn() == -fvc::grad(p));

fvOptions.correct(U);
*/

// Uhat = H(U0)/A
Info<< "Calculating Uhat" << endl;
volVectorField Uhat("Uhat", U);
Uhat = rAU*UEqn().H();

3. Solve pressure equation in order to obtain pStar

/*

    The same as the original pEqn.H except that HbyA is replaced by Uhat 

*/

4.  Solve momentum equation using momentum coefficients calculated in Step 1 and pStar for Ustar

Info<< "Solving implicitly for Ustar" << endl;

fvOptions.constrain(UEqn());

solve(UEqn() == -fvc::grad(p));

fvOptions.correct(U);

5. Solve pressure equation again for p

// Store pStar

volScalarField pTmp(p);

/*

    The same as the original pEqn.H

*/

6. Calculate U = (H(Ustar) - grad(p))/A

// Explicitly relax pressure for momentum corrector

p.relax();

U = HbyA - rAU*fvc::grad(p);

U.correctBoundaryConditions();

fvOptions.correct(U);

7. Set p = pStar

Info<< "Restoring pStar" << endl;

p = pTmp;

8. Return to Step 1

・ Validation with a tutorial case

(angledDuctExplicitFixedCoeff)

pressure

temperature

velocity

Complete source codes are shared in the storage page.