diff --git a/applications/solvers/compressible/steadyCompressibleFoam/createFields.H b/applications/solvers/compressible/steadyCompressibleFoam/createFields.H
index e1d363a75..9bf40b3a6 100644
--- a/applications/solvers/compressible/steadyCompressibleFoam/createFields.H
+++ b/applications/solvers/compressible/steadyCompressibleFoam/createFields.H
@@ -9,7 +9,8 @@
     volScalarField& p = thermo.p();
     volScalarField& h = thermo.h();
     const volScalarField& T = thermo.T();
-    const volScalarField& psi = thermo.psi();
+    volScalarField psis("psi", thermo.psi()/thermo.Cp()*thermo.Cv());
+    psis.oldTime();
 
     volScalarField rho
     (
diff --git a/applications/solvers/compressible/steadyCompressibleFoam/hEqn.H b/applications/solvers/compressible/steadyCompressibleFoam/hEqn.H
index 780c9401d..91a8e1dc3 100644
--- a/applications/solvers/compressible/steadyCompressibleFoam/hEqn.H
+++ b/applications/solvers/compressible/steadyCompressibleFoam/hEqn.H
@@ -34,4 +34,5 @@
     h.correctBoundaryConditions();
 
     thermo.correct();
+    psis = thermo.psi()/thermo.Cp()*thermo.Cv();
 }
diff --git a/applications/solvers/compressible/steadyCompressibleFoam/pEqn.H b/applications/solvers/compressible/steadyCompressibleFoam/pEqn.H
index 5a214ca34..a85bad7be 100644
--- a/applications/solvers/compressible/steadyCompressibleFoam/pEqn.H
+++ b/applications/solvers/compressible/steadyCompressibleFoam/pEqn.H
@@ -1,24 +1,22 @@
 {
     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();
 
-        // Execute ddtPhiCorr before recalculating flux
-        // HJ, 27/Apr/2010
-        surfaceScalarField phid
-        (
-            "phid",
-            fvc::interpolate(psi)*
-            (
-                (fvc::interpolate(U) & mesh.Sf())
-              + fvc::ddtPhiCorr(rUA, rho, U, phi)
-            )
-        );
-
         // Calculate phi for boundary conditions
-        phi = fvc::interpolate(rho*U) & mesh.Sf();
+        phi = rhof*fvc::interpolate(U) & mesh.Sf();
+
+        surfaceScalarField phid2 = rhoReff/rhof*phi;
+
+        surfaceScalarField phid("phid", psisf/rhof*phi);
 
         p.storePrevIter();
 
@@ -26,8 +24,9 @@
         {
             fvScalarMatrix pEqn
             (
-                fvm::ddt(psi, p)
+                fvm::ddt(psis, p)
               + fvm::div(phid, p)
+              + fvc::div(phid2)
               - fvm::laplacian(rho*rUA, p)
             );
 
@@ -42,7 +41,7 @@
             // Calculate the flux
             if (nonOrth == nNonOrthCorr)
             {
-                phi = pEqn.flux();
+                phi = phid2 + pEqn.flux();
             }
         }
 
diff --git a/applications/solvers/compressible/steadyCompressibleMRFFoam/pEqn.H b/applications/solvers/compressible/steadyCompressibleMRFFoam/pEqn.H
index 1cd8b9efb..b45e1573d 100644
--- a/applications/solvers/compressible/steadyCompressibleMRFFoam/pEqn.H
+++ b/applications/solvers/compressible/steadyCompressibleMRFFoam/pEqn.H
@@ -11,13 +11,17 @@
     {
         U = rUA*UEqn.H();
 
-        phi = rhoReff*(fvc::interpolate(U) & mesh.Sf());
+        // Calculate phi for boundary conditions
+        phi = rhof*fvc::interpolate(U) & mesh.Sf();
 
-        surfaceScalarField phid("phid", psisf/rhoReff*phi);
+        surfaceScalarField phid2 = rhoReff/rhof*phi;
+
+        surfaceScalarField phid("phid", psisf/rhof*phi);
 
         // Make fluxes relative within the MRF zone
-        mrfZones.relativeFlux(psisf, phid);
         mrfZones.relativeFlux(rhoReff, phi);
+        mrfZones.relativeFlux(psisf, phid);
+        mrfZones.relativeFlux(rhoReff, phid2);
 
         p.storePrevIter();
 
@@ -27,7 +31,7 @@
             (
                 fvm::ddt(psis, p)
               + fvm::div(phid, p)
-              + fvc::div(phi)
+              + fvc::div(phid2)
               - fvm::laplacian(rho*rUA, p)
             );
 
@@ -42,13 +46,13 @@
             // Calculate the flux
             if (nonOrth == nNonOrthCorr)
             {
-                phi += pEqn.flux();
+                phi = phid2 + pEqn.flux();
             }
         }
 
 #       include "compressibleContinuityErrs.H"
 
-        // Explicitly relax the pressure for momentum corrector
+        // Relax the pressure
         p.relax();
 
         U -= rUA*fvc::grad(p);
diff --git a/applications/solvers/compressible/steadyCompressibleSRFFoam/createFields.H b/applications/solvers/compressible/steadyCompressibleSRFFoam/createFields.H
index 49b4b0b97..6b84b46ba 100644
--- a/applications/solvers/compressible/steadyCompressibleSRFFoam/createFields.H
+++ b/applications/solvers/compressible/steadyCompressibleSRFFoam/createFields.H
@@ -9,7 +9,8 @@
     volScalarField& p = thermo.p();
     volScalarField& h = thermo.h();
     const volScalarField& T = thermo.T();
-    const volScalarField& psi = thermo.psi();
+    volScalarField psis("psi", thermo.psi()/thermo.Cp()*thermo.Cv());
+    psis.oldTime();
 
     volScalarField rho
     (
diff --git a/applications/solvers/compressible/steadyCompressibleSRFFoam/iEqn.H b/applications/solvers/compressible/steadyCompressibleSRFFoam/iEqn.H
index c80d715ca..d91a1ffb3 100644
--- a/applications/solvers/compressible/steadyCompressibleSRFFoam/iEqn.H
+++ b/applications/solvers/compressible/steadyCompressibleSRFFoam/iEqn.H
@@ -1,5 +1,5 @@
 {
-    // Solve the enthalpy equation
+    // Solve the rothalpy equation
     T.storePrevIter();
 
     // Calculate face velocity from flux
@@ -44,4 +44,5 @@
     i = h - 0.5*magSqr(Urot);
 
     thermo.correct();
+    psis = thermo.psi()/thermo.Cp()*thermo.Cv();
 }
diff --git a/applications/solvers/compressible/steadyCompressibleSRFFoam/pEqn.H b/applications/solvers/compressible/steadyCompressibleSRFFoam/pEqn.H
index ed734fc7d..f161a3c4d 100644
--- a/applications/solvers/compressible/steadyCompressibleSRFFoam/pEqn.H
+++ b/applications/solvers/compressible/steadyCompressibleSRFFoam/pEqn.H
@@ -1,24 +1,22 @@
 {
     volScalarField rUrelA = 1.0/UrelEqn.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++)
     {
         Urel = rUrelA*UrelEqn.H();
 
-        // Execute ddtPhiCorr before recalculating flux
-        // HJ, 27/Apr/2010
-        surfaceScalarField phid
-        (
-            "phid",
-            fvc::interpolate(thermo.psi())*
-            (
-                (fvc::interpolate(Urel) & mesh.Sf())
-              + fvc::ddtPhiCorr(rUrelA, rho, Urel, phi)
-            )
-        );
-
         // Calculate phi for boundary conditions
-        phi = fvc::interpolate(rho*Urel) & mesh.Sf();
+        phi = rhof*fvc::interpolate(Urel) & mesh.Sf();
+
+        surfaceScalarField phid2 = rhoReff/rhof*phi;
+
+        surfaceScalarField phid("phid", psisf/rhof*phi);
 
         p.storePrevIter();
 
@@ -26,8 +24,9 @@
         {
             fvScalarMatrix pEqn
             (
-                fvm::ddt(psi, p)
+                fvm::ddt(psis, p)
               + fvm::div(phid, p)
+              + fvc::div(phid2)
               - fvm::laplacian(rho*rUrelA, p)
             );
 
@@ -42,14 +41,13 @@
             // Calculate the flux
             if (nonOrth == nNonOrthCorr)
             {
-                phi = pEqn.flux();
+                phi = phid2 + pEqn.flux();
             }
         }
 
 #       include "compressibleContinuityErrs.H"
 
-
-        // Explicitly relax the pressure for momentum corrector
+        // Relax the pressure
         p.relax();
 
         Urel -= rUrelA*fvc::grad(p);
diff --git a/tutorials/compressible/steadyCompressibleFoam/2bump/system/fvSolution b/tutorials/compressible/steadyCompressibleFoam/2bump/system/fvSolution
index c567ebbed..cf8cf1e2c 100644
--- a/tutorials/compressible/steadyCompressibleFoam/2bump/system/fvSolution
+++ b/tutorials/compressible/steadyCompressibleFoam/2bump/system/fvSolution
@@ -58,11 +58,9 @@ PIMPLE
 relaxationFactors
 {
     // Note: under-relaxation factors used in wave-transmissive schemes
-    U               0.5;
-    p               0.2;
-    h               0.5;
-    rho             0.5;
-    T               1;
+    U               0.7;
+    p               0.7;
+    h               0.7;
 }
 
 fieldBounds
diff --git a/tutorials/compressible/steadyCompressibleFoam/transonicBump/system/fvSchemes b/tutorials/compressible/steadyCompressibleFoam/transonicBump/system/fvSchemes
index 1bb526f98..8c0c582e2 100644
--- a/tutorials/compressible/steadyCompressibleFoam/transonicBump/system/fvSchemes
+++ b/tutorials/compressible/steadyCompressibleFoam/transonicBump/system/fvSchemes
@@ -20,7 +20,7 @@ ddtSchemes
 
     ddt(rho,U)      steadyState;
     ddt(rho,h)      steadyState;
-    ddt(psi,p)      steadyInertial phi rho 0.5;
+    ddt(psi,p)      steadyInertial phi rho 0.95;
 
     U               steadyState;
     T               steadyState;
diff --git a/tutorials/compressible/steadyCompressibleFoam/transonicBump/system/fvSolution b/tutorials/compressible/steadyCompressibleFoam/transonicBump/system/fvSolution
index 5a67d8053..e499e4b91 100644
--- a/tutorials/compressible/steadyCompressibleFoam/transonicBump/system/fvSolution
+++ b/tutorials/compressible/steadyCompressibleFoam/transonicBump/system/fvSolution
@@ -76,8 +76,6 @@ relaxationFactors
     U               0.1;
     p               0.25;
     h               0.1;
-    rho             0.5;
-    T               0.5;
 }
 
 fieldBounds