diff --git a/applications/solvers/compressible/rhoPimpleFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
index 6ea8cc76e3..81676f5fa7 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
@@ -1,9 +1,5 @@
 rho = thermo.rho();
-
-if (!pimple.transonic())
-{
-    rho.relax();
-}
+rho.relax();
 
 // Thermodynamic density needs to be updated by psi*d(p) after the
 // pressure solution
@@ -53,13 +49,14 @@ constrainPressure(p, rho, U, phiHbyA, rhorAAtUf, MRF);
 
 if (pimple.transonic())
 {
-    surfaceScalarField phid
+    const surfaceScalarField phid
     (
         "phid",
         (fvc::interpolate(psi)/fvc::interpolate(rho))*phiHbyA
     );
 
-    phiHbyA -= fvc::interpolate(psi*p)*phiHbyA/fvc::interpolate(rho);
+    const fvScalarMatrix divPhidp(fvm::div(phid, p));
+    phiHbyA -= divPhidp.flux();
 
     if (pimple.consistent())
     {
@@ -70,7 +67,7 @@ if (pimple.transonic())
     fvScalarMatrix pDDtEqn
     (
         fvc::ddt(rho) + psi*correction(fvm::ddt(p))
-      + fvc::div(phiHbyA) + fvm::div(phid, p)
+      + fvc::div(phiHbyA) + divPhidp
      ==
         fvModels.source(psi, p, rho.name())
     );
@@ -186,7 +183,7 @@ if (mesh.steady() || pimple.simpleRho() || adjustMass)
 // Correct rhoUf if the mesh is moving
 fvc::correctRhoUf(rhoUf, rho, U, phi);
 
-if ((mesh.steady() || pimple.simpleRho()) && !pimple.transonic())
+if (mesh.steady() || pimple.simpleRho())
 {
     rho.relax();
 }
diff --git a/applications/solvers/compressible/rhoSimpleFoam/pEqn.H b/applications/solvers/compressible/rhoSimpleFoam/pEqn.H
index 6ea12ab40e..8812504130 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/pEqn.H
@@ -1,10 +1,5 @@
 rho = thermo.rho();
-
-if (!simple.transonic())
-{
-    rho.relax();
-}
-
+rho.relax();
 
 const volScalarField rAU("rAU", 1.0/UEqn.A());
 const surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
@@ -45,13 +40,14 @@ constrainPressure(p, rho, U, phiHbyA, rhorAAtUf, MRF);
 
 if (simple.transonic())
 {
-    surfaceScalarField phid
+    const surfaceScalarField phid
     (
         "phid",
-        (fvc::interpolate(psi)/fvc::interpolate(rho))*phiHbyA
+        fvc::interpolate(psi)*phiHbyA/fvc::interpolate(rho)
     );
 
-    phiHbyA -= fvc::interpolate(psi*p)*phiHbyA/fvc::interpolate(rho);
+    const fvScalarMatrix divPhidp(fvm::div(phid, p));
+    phiHbyA -= divPhidp.flux();
 
     if (simple.consistent())
     {
@@ -63,8 +59,7 @@ if (simple.transonic())
     {
         fvScalarMatrix pEqn
         (
-            fvc::div(phiHbyA)
-          + fvm::div(phid, p)
+            fvc::div(phiHbyA) + divPhidp
           - fvm::laplacian(rhorAAtUf, p)
           ==
             fvModels.source(psi, p, rho.name())
@@ -143,8 +138,4 @@ if (closedVolume && !thermo.incompressible())
 }
 
 rho = thermo.rho();
-
-if (!simple.transonic())
-{
-    rho.relax();
-}
+rho.relax();
diff --git a/applications/solvers/heatTransfer/buoyantFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantFoam/pEqn.H
index f960c7831e..1e3497e755 100644
--- a/applications/solvers/heatTransfer/buoyantFoam/pEqn.H
+++ b/applications/solvers/heatTransfer/buoyantFoam/pEqn.H
@@ -1,9 +1,5 @@
 rho = thermo.rho();
-
-if (!pimple.transonic())
-{
-    rho.relax();
-}
+rho.relax();
 
 // Thermodynamic density needs to be updated by psi*d(p) after the
 // pressure solution
@@ -29,7 +25,7 @@ MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
 
 bool adjustMass = mesh.steady() && adjustPhi(phiHbyA, U, p_rgh);
 
-surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
+const surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
 
 phiHbyA += phig;
 
@@ -42,18 +38,19 @@ fvScalarMatrix p_rghEqn(p_rgh, dimMass/dimTime);
 
 if (pimple.transonic())
 {
-    surfaceScalarField phid
+    const surfaceScalarField phid
     (
         "phid",
         (fvc::interpolate(psi)/fvc::interpolate(rho))*phiHbyA
     );
 
-    phiHbyA -= fvc::interpolate(psi*p_rgh)*phiHbyA/fvc::interpolate(rho);
+    const fvScalarMatrix divPhidp(fvm::div(phid, p));
+    phiHbyA -= divPhidp.flux();
 
     fvScalarMatrix p_rghDDtEqn
     (
         fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
-      + fvc::div(phiHbyA) + fvm::div(phid, p_rgh)
+      + fvc::div(phiHbyA) + divPhidp
      ==
         fvModels.source(psi, p_rgh, rho.name())
     );
@@ -158,7 +155,7 @@ if (mesh.steady() || pimple.simpleRho() || adjustMass)
     rho = thermo.rho();
 }
 
-if ((mesh.steady() || pimple.simpleRho()) && !pimple.transonic())
+if (mesh.steady() || pimple.simpleRho())
 {
     rho.relax();
 }
diff --git a/tutorials/compressible/rhoSimpleFoam/aerofoilNACA0012/constant/physicalProperties b/tutorials/compressible/rhoSimpleFoam/aerofoilNACA0012/constant/physicalProperties
index 80135919d0..6a52d18984 100644
--- a/tutorials/compressible/rhoSimpleFoam/aerofoilNACA0012/constant/physicalProperties
+++ b/tutorials/compressible/rhoSimpleFoam/aerofoilNACA0012/constant/physicalProperties
@@ -22,7 +22,7 @@ thermoType
     thermo          hConst;
     equationOfState perfectGas;
     specie          specie;
-    energy          sensibleInternalEnergy;
+    energy          sensibleEnthalpy;
 }
 
 mixture // air at room temperature (293 K)
@@ -31,11 +31,13 @@ mixture // air at room temperature (293 K)
     {
         molWeight   28.9;
     }
+
     thermodynamics
     {
         Cp          1005;
         Hf          0;
     }
+
     transport
     {
         mu          1.82e-05;
diff --git a/tutorials/compressible/rhoSimpleFoam/aerofoilNACA0012/system/controlDict b/tutorials/compressible/rhoSimpleFoam/aerofoilNACA0012/system/controlDict
index 28b62e0f95..3fd5104525 100644
--- a/tutorials/compressible/rhoSimpleFoam/aerofoilNACA0012/system/controlDict
+++ b/tutorials/compressible/rhoSimpleFoam/aerofoilNACA0012/system/controlDict
@@ -15,13 +15,13 @@ FoamFile
 
 application     rhoSimpleFoam;
 
-startFrom       latestTime;
+startFrom       startTime;
 
 startTime       0;
 
 stopAt          endTime;
 
-endTime         2000;
+endTime         5000;
 
 deltaT          1;
 
diff --git a/tutorials/compressible/rhoSimpleFoam/aerofoilNACA0012/system/fvSchemes b/tutorials/compressible/rhoSimpleFoam/aerofoilNACA0012/system/fvSchemes
index cccaf6441f..392a31de7a 100644
--- a/tutorials/compressible/rhoSimpleFoam/aerofoilNACA0012/system/fvSchemes
+++ b/tutorials/compressible/rhoSimpleFoam/aerofoilNACA0012/system/fvSchemes
@@ -34,19 +34,14 @@ divSchemes
 
     div(phi,U)      bounded Gauss linearUpwind limited;
 
-    turbulence      bounded Gauss upwind;
     energy          bounded Gauss linearUpwind limited;
+    div(phi,h)      $energy;
+    div(phi,K)      $energy;
 
+    turbulence      bounded Gauss upwind;
     div(phi,k)      $turbulence;
     div(phi,omega)  $turbulence;
 
-    div(phi,e)      $energy;
-    div(phi,K)      $energy;
-    div(phi,Ekp)    $energy;
-
-    div(phid,p)     Gauss upwind;
-    div((phi|interpolate(rho)),p)  bounded Gauss upwind;
-
     div(((rho*nuEff)*dev2(T(grad(U)))))    Gauss linear;
 }
 
diff --git a/tutorials/compressible/rhoSimpleFoam/aerofoilNACA0012/system/fvSolution b/tutorials/compressible/rhoSimpleFoam/aerofoilNACA0012/system/fvSolution
index 12aeed15f4..eb5a5295ef 100644
--- a/tutorials/compressible/rhoSimpleFoam/aerofoilNACA0012/system/fvSolution
+++ b/tutorials/compressible/rhoSimpleFoam/aerofoilNACA0012/system/fvSolution
@@ -18,12 +18,13 @@ solvers
     p
     {
         solver          GAMG;
-        smoother        GaussSeidel;
+        smoother        DIC;
+
         tolerance       1e-8;
         relTol          0.01;
     }
 
-    "(U|k|omega|e)"
+    "(U|k|omega|h)"
     {
         solver          PBiCGStab;
         preconditioner  DILU;
@@ -38,7 +39,7 @@ SIMPLE
     {
         p               1e-6;
         U               1e-5;
-        "(k|omega|e)"   1e-5;
+        "(k|omega|h)"   1e-5;
     }
 
     nNonOrthogonalCorrectors 0;
@@ -51,10 +52,11 @@ relaxationFactors
         p               0.3;
         rho             0.01;
     }
+
     equations
     {
         U               0.7;
-        e               0.7;
+        h               0.7;
         "(k|omega)"     0.7;
     }
 }
diff --git a/tutorials/compressible/rhoSimpleFoam/angledDuctExplicitFixedCoeff/constant/physicalProperties b/tutorials/compressible/rhoSimpleFoam/angledDuctExplicitFixedCoeff/constant/physicalProperties
index bc04ce5d9d..5d778e9433 100644
--- a/tutorials/compressible/rhoSimpleFoam/angledDuctExplicitFixedCoeff/constant/physicalProperties
+++ b/tutorials/compressible/rhoSimpleFoam/angledDuctExplicitFixedCoeff/constant/physicalProperties
@@ -22,7 +22,7 @@ thermoType
     thermo          hConst;
     equationOfState perfectGas;
     specie          specie;
-    energy          sensibleInternalEnergy;
+    energy          sensibleEnthalpy;
 }
 
 mixture
diff --git a/tutorials/compressible/rhoSimpleFoam/angledDuctExplicitFixedCoeff/system/fvSchemes b/tutorials/compressible/rhoSimpleFoam/angledDuctExplicitFixedCoeff/system/fvSchemes
index c44866691c..a4b3853b38 100644
--- a/tutorials/compressible/rhoSimpleFoam/angledDuctExplicitFixedCoeff/system/fvSchemes
+++ b/tutorials/compressible/rhoSimpleFoam/angledDuctExplicitFixedCoeff/system/fvSchemes
@@ -29,11 +29,14 @@ divSchemes
     default         none;
 
     div(phi,U)      bounded Gauss upwind;
-    div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
-    div(phi,e)      bounded Gauss upwind;
+
+    div(phi,h)      bounded Gauss upwind;
+    div(phi,K)      bounded Gauss upwind;
+
     div(phi,epsilon) bounded Gauss upwind;
     div(phi,k)      bounded Gauss upwind;
-    div(phi,Ekp)    bounded Gauss upwind;
+
+    div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes
diff --git a/tutorials/compressible/rhoSimpleFoam/angledDuctExplicitFixedCoeff/system/fvSolution b/tutorials/compressible/rhoSimpleFoam/angledDuctExplicitFixedCoeff/system/fvSolution
index 3a52cc71cc..2b75e75d44 100644
--- a/tutorials/compressible/rhoSimpleFoam/angledDuctExplicitFixedCoeff/system/fvSolution
+++ b/tutorials/compressible/rhoSimpleFoam/angledDuctExplicitFixedCoeff/system/fvSolution
@@ -29,22 +29,24 @@ solvers
     {
         solver          smoothSolver;
         smoother        GaussSeidel;
-        nSweeps         2;
-        tolerance       1e-06;
+
+        tolerance       1e-6;
         relTol          0.1;
     }
 
-    e
+    h
     {
         solver          smoothSolver;
         smoother        symGaussSeidel;
-        tolerance       1e-06;
+
+        tolerance       1e-6;
         relTol          0.1;
     }
 
     "(k|epsilon)"
     {
         $U;
+
         tolerance       1e-07;
         relTol          0.1;
     }
@@ -59,9 +61,7 @@ SIMPLE
         p               1e-2;
         U               1e-4;
         e               1e-3;
-
-        // Possibly check turbulence fields
-        "(k|epsilon|omega)" 1e-3;
+        "(k|epsilon)"   1e-3;
     }
 }
 
@@ -72,10 +72,11 @@ relaxationFactors
         p               0.3;
         rho             0.01;
     }
+
     equations
     {
         U               0.7;
-        e               0.7;
+        h               0.7;
         "(k|epsilon)"   0.7;
     }
 }
diff --git a/tutorials/compressible/rhoSimpleFoam/squareBend/0/U b/tutorials/compressible/rhoSimpleFoam/squareBend/0/U
index 2f43261bde..46c339f698 100644
--- a/tutorials/compressible/rhoSimpleFoam/squareBend/0/U
+++ b/tutorials/compressible/rhoSimpleFoam/squareBend/0/U
@@ -23,6 +23,7 @@ boundaryField
     {
         type            noSlip;
     }
+
     inlet
     {
         type            flowRateInletVelocity;
@@ -30,11 +31,11 @@ boundaryField
         rhoInlet        0.5;    // Guess for rho
         profile         turbulentBL;
     }
+
     outlet
     {
-        type            inletOutlet;
+        type            pressureInletOutletVelocity;
         value           uniform (0 0 0);
-        inletValue      uniform (0 0 0);
     }
 }
 
diff --git a/tutorials/compressible/rhoSimpleFoam/squareBend/constant/physicalProperties b/tutorials/compressible/rhoSimpleFoam/squareBend/constant/physicalProperties
index bc04ce5d9d..5d778e9433 100644
--- a/tutorials/compressible/rhoSimpleFoam/squareBend/constant/physicalProperties
+++ b/tutorials/compressible/rhoSimpleFoam/squareBend/constant/physicalProperties
@@ -22,7 +22,7 @@ thermoType
     thermo          hConst;
     equationOfState perfectGas;
     specie          specie;
-    energy          sensibleInternalEnergy;
+    energy          sensibleEnthalpy;
 }
 
 mixture
diff --git a/tutorials/compressible/rhoSimpleFoam/squareBend/system/fvSchemes b/tutorials/compressible/rhoSimpleFoam/squareBend/system/fvSchemes
index f8c15ed765..38a7bf1b68 100644
--- a/tutorials/compressible/rhoSimpleFoam/squareBend/system/fvSchemes
+++ b/tutorials/compressible/rhoSimpleFoam/squareBend/system/fvSchemes
@@ -29,14 +29,16 @@ divSchemes
     default             none;
 
     div(phi,U)          bounded Gauss upwind;
-    div(((rho*nuEff)*dev2(T(grad(U)))))      Gauss linear;
-    div(phi,e)          bounded Gauss upwind;
+    div(phi,h)          bounded Gauss upwind;
+
     div(phi,epsilon)    bounded Gauss upwind;
     div(phi,k)          bounded Gauss upwind;
 
+    div(phi,K)          bounded Gauss upwind;
+
     div(phid,p)         Gauss upwind;
-    div(phi,Ekp)        bounded Gauss upwind;
-    div((phi|interpolate(rho)),p)  Gauss upwind;
+
+    div(((rho*nuEff)*dev2(T(grad(U)))))      Gauss linear;
 }
 
 laplacianSchemes
diff --git a/tutorials/compressible/rhoSimpleFoam/squareBend/system/fvSolution b/tutorials/compressible/rhoSimpleFoam/squareBend/system/fvSolution
index 3a9721daf9..1491724c91 100644
--- a/tutorials/compressible/rhoSimpleFoam/squareBend/system/fvSolution
+++ b/tutorials/compressible/rhoSimpleFoam/squareBend/system/fvSolution
@@ -19,19 +19,19 @@ solvers
     p
     {
         solver          GAMG;
+        smoother        DILU;
+
         tolerance       1e-08;
         relTol          0.1;
-        smoother        GaussSeidel;
-        nCellsInCoarsestLevel 20;
     }
 
-    "(U|e|k|epsilon)"
+    "(U|h|k|epsilon)"
     {
-        solver          GAMG;
+        solver          PBiCGStab;
+        preconditioner  DILU;
+
         tolerance       1e-08;
         relTol          0.1;
-        smoother        GaussSeidel;
-        nCellsInCoarsestLevel 20;
     }
 }
 
@@ -46,10 +46,8 @@ SIMPLE
     {
         p               1e-3;
         U               1e-4;
-        e               1e-3;
-
-        // possibly check turbulence fields
-        "(k|epsilon|omega)" 1e-3;
+        h               1e-3;
+        "(k|epsilon)"   1e-3;
     }
 }
 
@@ -59,11 +57,12 @@ relaxationFactors
     {
         p               1;
     }
+
     equations
     {
         p               1;
         U               0.9;
-        e               0.8;
+        h               0.9;
         k               0.9;
         epsilon         0.9;
     }
diff --git a/tutorials/compressible/rhoSimpleFoam/squareBendLiq/constant/physicalProperties b/tutorials/compressible/rhoSimpleFoam/squareBendLiq/constant/physicalProperties
index b3238afa33..c0b11b5c3a 100644
--- a/tutorials/compressible/rhoSimpleFoam/squareBendLiq/constant/physicalProperties
+++ b/tutorials/compressible/rhoSimpleFoam/squareBendLiq/constant/physicalProperties
@@ -19,7 +19,7 @@ thermoType
     type            heRhoThermo;
     mixture         pureMixture;
     properties      liquid;
-    energy          sensibleInternalEnergy;
+    energy          sensibleEnthalpy;
 }
 
 mixture
diff --git a/tutorials/compressible/rhoSimpleFoam/squareBendLiq/system/fvSchemes b/tutorials/compressible/rhoSimpleFoam/squareBendLiq/system/fvSchemes
index ad0ae211ec..6b871313a1 100644
--- a/tutorials/compressible/rhoSimpleFoam/squareBendLiq/system/fvSchemes
+++ b/tutorials/compressible/rhoSimpleFoam/squareBendLiq/system/fvSchemes
@@ -31,12 +31,12 @@ divSchemes
     default             none;
 
     div(phi,U)          bounded Gauss linearUpwindV limited;
+
     div(phi,h)          bounded Gauss linearUpwind limited;
-    div(phi,e)          bounded Gauss linearUpwind limited;
+    div(phi,K)          bounded Gauss linearUpwind limited;
+
     div(phi,epsilon)    bounded Gauss linearUpwind limited;
     div(phi,k)          bounded Gauss linearUpwind limited;
-    div(phi,Ekp)        bounded Gauss linearUpwind limited;
-    div(phi,K)          bounded Gauss linearUpwind limited;
 
     div(((rho*nuEff)*dev2(T(grad(U)))))      Gauss linear;
 }
diff --git a/tutorials/compressible/rhoSimpleFoam/squareBendLiq/system/fvSolution b/tutorials/compressible/rhoSimpleFoam/squareBendLiq/system/fvSolution
index 911b47bac9..53026f793e 100644
--- a/tutorials/compressible/rhoSimpleFoam/squareBendLiq/system/fvSolution
+++ b/tutorials/compressible/rhoSimpleFoam/squareBendLiq/system/fvSolution
@@ -25,7 +25,7 @@ solvers
         relTol          0.1;
     }
 
-    "(U|h|e|k|epsilon)"
+    "(U|h|k|epsilon)"
     {
         solver          smoothSolver;
         smoother        symGaussSeidel;
@@ -46,7 +46,7 @@ SIMPLE
     {
         p               1e-4;
         U               1e-4;
-        "(k|omega|epsilon|h|e)" 1e-3;
+        "(k|omega|epsilon|h)" 1e-3;
     }
 }
 
@@ -56,11 +56,11 @@ relaxationFactors
     {
         p               0.3;
     }
+
     equations
     {
         U               0.7;
         h               0.7;
-        e               0.7;
         k               0.7;
         epsilon         0.7;
     }