diff --git a/applications/solvers/compressible/rhoPimpleFoam/createFields.H b/applications/solvers/compressible/rhoPimpleFoam/createFields.H
index a7ee3eca45..84f70f2cd0 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/createFields.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/createFields.H
@@ -2,11 +2,11 @@
 
 Info<< "Reading thermophysical properties\n" << endl;
 
-autoPtr<psiThermo> pThermo
+autoPtr<fluidThermo> pThermo
 (
-    psiThermo::New(mesh)
+    fluidThermo::New(mesh)
 );
-psiThermo& thermo = pThermo();
+fluidThermo& thermo = pThermo();
 thermo.validate(args.executable(), "h", "e");
 
 volScalarField& p = thermo.p();
@@ -40,27 +40,7 @@ volVectorField U
 
 #include "compressibleCreatePhi.H"
 
-dimensionedScalar rhoMax
-(
-    dimensionedScalar::lookupOrDefault
-    (
-        "rhoMax",
-        pimple.dict(),
-        dimDensity,
-        GREAT
-    )
-);
-
-dimensionedScalar rhoMin
-(
-    dimensionedScalar::lookupOrDefault
-    (
-        "rhoMin",
-        pimple.dict(),
-        dimDensity,
-        0
-    )
-);
+pressureControl pressureControl(p, rho, pimple.dict(), false);
 
 Info<< "Creating turbulence model\n" << endl;
 autoPtr<compressible::turbulenceModel> turbulence
diff --git a/applications/solvers/compressible/rhoPimpleFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
index ac7107acf0..73ecafd89f 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
@@ -1,8 +1,3 @@
-rho = thermo.rho();
-rho = max(rho, rhoMin);
-rho = min(rho, rhoMax);
-rho.relax();
-
 volScalarField rAU(1.0/UEqn.A());
 surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
 volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
@@ -12,55 +7,54 @@ if (pimple.nCorrPISO() <= 1)
     tUEqn.clear();
 }
 
+surfaceScalarField phiHbyA
+(
+    "phiHbyA",
+    (
+        fvc::flux(rho*HbyA)
+      + rhorAUf*fvc::ddtCorr(rho, U, phi)
+    )
+);
+
+MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
+
+// Update the pressure BCs to ensure flux consistency
+constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF);
+
 if (pimple.transonic())
 {
     surfaceScalarField phid
     (
         "phid",
-        fvc::interpolate(psi)
-       *(
-            fvc::flux(HbyA)
-          + rhorAUf*fvc::ddtCorr(rho, U, phi)/fvc::interpolate(rho)
-        )
+        (fvc::interpolate(psi)/fvc::interpolate(rho))*phiHbyA
     );
-
-    MRF.makeRelative(fvc::interpolate(psi), phid);
+    phiHbyA -= fvc::interpolate(p)*phid;
 
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
             fvm::ddt(psi, p)
+          + fvc::div(phiHbyA)
           + fvm::div(phid, p)
           - fvm::laplacian(rhorAUf, p)
          ==
             fvOptions(psi, p, rho.name())
         );
 
+        // Relax the pressure equation to ensure diagonal-dominance
+        pEqn.relax();
+
         pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
         if (pimple.finalNonOrthogonalIter())
         {
-            phi == pEqn.flux();
+            phi = phiHbyA + pEqn.flux();
         }
     }
 }
 else
 {
-    surfaceScalarField phiHbyA
-    (
-        "phiHbyA",
-        (
-            fvc::flux(rho*HbyA)
-          + rhorAUf*fvc::ddtCorr(rho, U, phi)
-        )
-    );
-
-    MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
-
-    // Update the pressure BCs to ensure flux consistency
-    constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF);
-
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
@@ -87,19 +81,20 @@ else
 // Explicitly relax pressure for momentum corrector
 p.relax();
 
-// Recalculate density from the relaxed pressure
-rho = thermo.rho();
-rho = max(rho, rhoMin);
-rho = min(rho, rhoMax);
-rho.relax();
-Info<< "rho max/min : " << max(rho).value()
-    << " " << min(rho).value() << endl;
-
 U = HbyA - rAU*fvc::grad(p);
 U.correctBoundaryConditions();
 fvOptions.correct(U);
 K = 0.5*magSqr(U);
 
+pressureControl.limit(p);
+p.correctBoundaryConditions();
+rho = thermo.rho();
+
+if (!pimple.transonic())
+{
+    rho.relax();
+}
+
 if (thermo.dpdt())
 {
     dpdt = fvc::ddt(p);
diff --git a/applications/solvers/compressible/rhoPimpleFoam/pcEqn.H b/applications/solvers/compressible/rhoPimpleFoam/pcEqn.H
index 713f443fc5..afbc2851e4 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/pcEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/pcEqn.H
@@ -1,8 +1,3 @@
-rho = thermo.rho();
-rho = max(rho, rhoMin);
-rho = min(rho, rhoMax);
-rho.relax();
-
 volScalarField rAU(1.0/UEqn.A());
 volScalarField rAtU(1.0/(1.0/rAU - UEqn.H1()));
 volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
@@ -12,72 +7,64 @@ if (pimple.nCorrPISO() <= 1)
     tUEqn.clear();
 }
 
+surfaceScalarField phiHbyA
+(
+    "phiHbyA",
+    (
+        fvc::flux(rho*HbyA)
+      + fvc::interpolate(rho*rAU)*fvc::ddtCorr(rho, U, phi)
+    )
+);
+
+MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
+
+volScalarField rhorAtU("rhorAtU", rho*rAtU);
+
+// Update the pressure BCs to ensure flux consistency
+constrainPressure(p, rho, U, phiHbyA, rhorAtU, MRF);
+
 if (pimple.transonic())
 {
     surfaceScalarField phid
     (
         "phid",
-        fvc::interpolate(psi)
-       *(
-            fvc::flux(HbyA)
-          + fvc::interpolate(rho*rAU)*fvc::ddtCorr(rho, U, phi)
-           /fvc::interpolate(rho)
-        )
+        (fvc::interpolate(psi)/fvc::interpolate(rho))*phiHbyA
     );
 
-    MRF.makeRelative(fvc::interpolate(psi), phid);
-
-    surfaceScalarField phic
-    (
-        "phic",
+    phiHbyA +=
         fvc::interpolate(rho*(rAtU - rAU))*fvc::snGrad(p)*mesh.magSf()
-    );
+      - fvc::interpolate(p)*phid;
 
     HbyA -= (rAU - rAtU)*fvc::grad(p);
 
-    volScalarField rhorAtU("rhorAtU", rho*rAtU);
-
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
             fvm::ddt(psi, p)
+          + fvc::div(phiHbyA)
           + fvm::div(phid, p)
-          + fvc::div(phic)
           - fvm::laplacian(rhorAtU, p)
          ==
             fvOptions(psi, p, rho.name())
         );
 
+        // Relax the pressure equation to ensure diagonal-dominance
+        pEqn.relax();
+
         pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
         if (pimple.finalNonOrthogonalIter())
         {
-            phi == phic + pEqn.flux();
+            phi = phiHbyA + pEqn.flux();
         }
     }
 }
 else
 {
-    surfaceScalarField phiHbyA
-    (
-        "phiHbyA",
-        (
-            fvc::flux(rho*HbyA)
-          + fvc::interpolate(rho*rAU)*fvc::ddtCorr(rho, U, phi)
-        )
-    );
-
-    MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
-
     phiHbyA += fvc::interpolate(rho*(rAtU - rAU))*fvc::snGrad(p)*mesh.magSf();
     HbyA -= (rAU - rAtU)*fvc::grad(p);
 
-    volScalarField rhorAtU("rhorAtU", rho*rAtU);
-
-    // Update the pressure BCs to ensure flux consistency
-    constrainPressure(p, rho, U, phiHbyA, rhorAtU, MRF);
-
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
@@ -109,19 +96,16 @@ U.correctBoundaryConditions();
 fvOptions.correct(U);
 K = 0.5*magSqr(U);
 
-if (thermo.dpdt())
-{
-    dpdt = fvc::ddt(p);
-}
-
-// Recalculate density from the relaxed pressure
+pressureControl.limit(p);
+p.correctBoundaryConditions();
 rho = thermo.rho();
-rho = max(rho, rhoMin);
-rho = min(rho, rhoMax);
 
 if (!pimple.transonic())
 {
     rho.relax();
 }
 
-Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl;
+if (thermo.dpdt())
+{
+    dpdt = fvc::ddt(p);
+}
diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H
index 7bd540df40..0f85562618 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H
@@ -1,8 +1,3 @@
-rho = thermo.rho();
-rho = max(rho, rhoMin);
-rho = min(rho, rhoMax);
-rho.relax();
-
 volScalarField rAU(1.0/UEqn.A());
 surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
 volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
@@ -12,55 +7,53 @@ if (pimple.nCorrPISO() <= 1)
     tUEqn.clear();
 }
 
+surfaceScalarField phiHbyA
+(
+    "phiHbyA",
+    fvc::flux(rho*HbyA)
+  + rhorAUf*fvc::ddtCorr(rho, U, rhoUf)
+);
+
+fvc::makeRelative(phiHbyA, rho, U);
+MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
+
+// Update the pressure BCs to ensure flux consistency
+constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF);
+
 if (pimple.transonic())
 {
     surfaceScalarField phid
     (
         "phid",
-        fvc::interpolate(psi)
-       *(
-            fvc::flux(HbyA)
-          + rhorAUf*fvc::ddtCorr(rho, U, rhoUf)/fvc::interpolate(rho)
-        )
+        (fvc::interpolate(psi)/fvc::interpolate(rho))*phiHbyA
     );
-
-    fvc::makeRelative(phid, psi, U);
-    MRF.makeRelative(fvc::interpolate(psi), phid);
+    phiHbyA -= fvc::interpolate(p)*phid;
 
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
             fvm::ddt(psi, p)
+          + fvc::div(phiHbyA)
           + fvm::div(phid, p)
           - fvm::laplacian(rhorAUf, p)
          ==
             fvOptions(psi, p, rho.name())
         );
 
+        // Relax the pressure equation to ensure diagonal-dominance
+        pEqn.relax();
+
         pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
         if (pimple.finalNonOrthogonalIter())
         {
-            phi == pEqn.flux();
+            phi = phiHbyA + pEqn.flux();
         }
     }
 }
 else
 {
-    surfaceScalarField phiHbyA
-    (
-        "phiHbyA",
-        fvc::flux(rho*HbyA)
-      + rhorAUf*fvc::ddtCorr(rho, U, rhoUf)
-    );
-
-    fvc::makeRelative(phiHbyA, rho, U);
-    MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
-
-    // Update the pressure BCs to ensure flux consistency
-    constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF);
-
     while (pimple.correctNonOrthogonal())
     {
         // Pressure corrector
@@ -88,19 +81,20 @@ else
 // Explicitly relax pressure for momentum corrector
 p.relax();
 
-// Recalculate density from the relaxed pressure
-rho = thermo.rho();
-rho = max(rho, rhoMin);
-rho = min(rho, rhoMax);
-rho.relax();
-Info<< "rho max/min : " << max(rho).value()
-    << " " << min(rho).value() << endl;
-
 U = HbyA - rAU*fvc::grad(p);
 U.correctBoundaryConditions();
 fvOptions.correct(U);
 K = 0.5*magSqr(U);
 
+pressureControl.limit(p);
+p.correctBoundaryConditions();
+rho = thermo.rho();
+
+if (!pimple.transonic())
+{
+    rho.relax();
+}
+
 {
     rhoUf = fvc::interpolate(rho*U);
     surfaceVectorField n(mesh.Sf()/mesh.magSf());
diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/rhoPimpleDyMFoam.C b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/rhoPimpleDyMFoam.C
index c7046bb9b4..1d20c3f10f 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/rhoPimpleDyMFoam.C
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/rhoPimpleDyMFoam.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -22,10 +22,7 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Application
-    rhoPimpleFoam
-
-Group
-    grpCompressibleSolvers grpMovingMeshSolvers
+    rhoPimpleDyMFoam
 
 Description
     Transient solver for turbulent flow of compressible fluids for HVAC and
@@ -38,10 +35,11 @@ Description
 
 #include "fvCFD.H"
 #include "dynamicFvMesh.H"
-#include "psiThermo.H"
+#include "fluidThermo.H"
 #include "turbulentFluidThermoModel.H"
 #include "bound.H"
 #include "pimpleControl.H"
+#include "pressureControl.H"
 #include "CorrectPhi.H"
 #include "fvOptions.H"
 #include "localEulerDdtScheme.H"
diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C
index 108ad9aa0e..f30dae76a3 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -34,10 +34,11 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "psiThermo.H"
+#include "fluidThermo.H"
 #include "turbulentFluidThermoModel.H"
 #include "bound.H"
 #include "pimpleControl.H"
+#include "pressureControl.H"
 #include "fvOptions.H"
 #include "localEulerDdtScheme.H"
 #include "fvcSmooth.H"
diff --git a/applications/solvers/compressible/rhoSimpleFoam/createFields.H b/applications/solvers/compressible/rhoSimpleFoam/createFields.H
index 4671347b66..c7914f89d3 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/createFields.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/createFields.H
@@ -7,6 +7,8 @@ autoPtr<fluidThermo> pThermo
 fluidThermo& thermo = pThermo();
 thermo.validate(args.executable(), "h", "e");
 
+volScalarField& p = thermo.p();
+
 volScalarField rho
 (
     IOobject
@@ -20,8 +22,6 @@ volScalarField rho
     thermo.rho()
 );
 
-volScalarField& p = thermo.p();
-
 Info<< "Reading field U\n" << endl;
 volVectorField U
 (
diff --git a/applications/solvers/compressible/rhoSimpleFoam/pEqn.H b/applications/solvers/compressible/rhoSimpleFoam/pEqn.H
index 5f092ada2f..c7edbb9e87 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/pEqn.H
@@ -6,16 +6,18 @@
 
     bool closedVolume = false;
 
+    surfaceScalarField phiHbyA("phiHbyA", fvc::flux(rho*HbyA));
+    MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
+
+    // Update the pressure BCs to ensure flux consistency
+    constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF);
+
     if (simple.transonic())
     {
-        surfaceScalarField phiHbyA("phiHbyA", fvc::flux(rho*HbyA));
-        surfaceScalarField rhof(fvc::interpolate(rho));
-        MRF.makeRelative(rhof, phiHbyA);
-
         surfaceScalarField phid
         (
             "phid",
-            (fvc::interpolate(psi)/rhof)*phiHbyA
+            (fvc::interpolate(psi)/fvc::interpolate(rho))*phiHbyA
         );
         phiHbyA -= fvc::interpolate(p)*phid;
 
@@ -49,14 +51,8 @@
     }
     else
     {
-        surfaceScalarField phiHbyA("phiHbyA", fvc::flux(rho*HbyA));
-        MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
-
         closedVolume = adjustPhi(phiHbyA, U, p);
 
-        // Update the pressure BCs to ensure flux consistency
-        constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF);
-
         while (simple.correctNonOrthogonal())
         {
             fvScalarMatrix pEqn
diff --git a/applications/solvers/compressible/rhoSimpleFoam/pcEqn.H b/applications/solvers/compressible/rhoSimpleFoam/pcEqn.H
index af2005febd..c7dc0f864d 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/pcEqn.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/pcEqn.H
@@ -5,16 +5,20 @@ tUEqn.clear();
 
 bool closedVolume = false;
 
+surfaceScalarField phiHbyA("phiHbyA", fvc::flux(rho*HbyA));
+MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
+
+volScalarField rhorAtU("rhorAtU", rho*rAtU);
+
+// Update the pressure BCs to ensure flux consistency
+constrainPressure(p, rho, U, phiHbyA, rhorAtU, MRF);
+
 if (simple.transonic())
 {
-    surfaceScalarField phiHbyA("phiHbyA", fvc::flux(rho*HbyA));
-    surfaceScalarField rhof(fvc::interpolate(rho));
-    MRF.makeRelative(rhof, phiHbyA);
-
     surfaceScalarField phid
     (
         "phid",
-        (fvc::interpolate(psi)/rhof)*phiHbyA
+        (fvc::interpolate(psi)/fvc::interpolate(rho))*phiHbyA
     );
 
     phiHbyA +=
@@ -23,14 +27,12 @@ if (simple.transonic())
 
     HbyA -= (rAU - rAtU)*fvc::grad(p);
 
-    volScalarField rhorAtU("rhorAtU", rho*rAtU);
-
     while (simple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
-            fvm::div(phid, p)
-          + fvc::div(phiHbyA)
+            fvc::div(phiHbyA)
+          + fvm::div(phid, p)
           - fvm::laplacian(rhorAtU, p)
          ==
             fvOptions(psi, p, rho.name())
@@ -55,19 +57,11 @@ if (simple.transonic())
 }
 else
 {
-    surfaceScalarField phiHbyA("phiHbyA", fvc::flux(rho*HbyA));
-    MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
-
     closedVolume = adjustPhi(phiHbyA, U, p);
 
     phiHbyA += fvc::interpolate(rho*(rAtU - rAU))*fvc::snGrad(p)*mesh.magSf();
     HbyA -= (rAU - rAtU)*fvc::grad(p);
 
-    volScalarField rhorAtU("rhorAtU", rho*rAtU);
-
-    // Update the pressure BCs to ensure flux consistency
-    constrainPressure(p, rho, U, phiHbyA, rhorAtU, MRF);
-
     while (simple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousSimpleFoam/createFields.H b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousSimpleFoam/createFields.H
deleted file mode 100644
index 4671347b66..0000000000
--- a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousSimpleFoam/createFields.H
+++ /dev/null
@@ -1,59 +0,0 @@
-Info<< "Reading thermophysical properties\n" << endl;
-
-autoPtr<fluidThermo> pThermo
-(
-    fluidThermo::New(mesh)
-);
-fluidThermo& thermo = pThermo();
-thermo.validate(args.executable(), "h", "e");
-
-volScalarField rho
-(
-    IOobject
-    (
-        "rho",
-        runTime.timeName(),
-        mesh,
-        IOobject::READ_IF_PRESENT,
-        IOobject::AUTO_WRITE
-    ),
-    thermo.rho()
-);
-
-volScalarField& p = thermo.p();
-
-Info<< "Reading field U\n" << endl;
-volVectorField U
-(
-    IOobject
-    (
-        "U",
-        runTime.timeName(),
-        mesh,
-        IOobject::MUST_READ,
-        IOobject::AUTO_WRITE
-    ),
-    mesh
-);
-
-#include "compressibleCreatePhi.H"
-
-pressureControl pressureControl(p, rho, simple.dict());
-
-mesh.setFluxRequired(p.name());
-
-Info<< "Creating turbulence model\n" << endl;
-autoPtr<compressible::turbulenceModel> turbulence
-(
-    compressible::turbulenceModel::New
-    (
-        rho,
-        U,
-        phi,
-        thermo
-    )
-);
-
-dimensionedScalar initialMass = fvc::domainIntegrate(rho);
-
-#include "createMRF.H"
diff --git a/src/finiteVolume/cfdTools/general/pressureControl/pressureControl.C b/src/finiteVolume/cfdTools/general/pressureControl/pressureControl.C
index 101126a5e8..aca22e7699 100644
--- a/src/finiteVolume/cfdTools/general/pressureControl/pressureControl.C
+++ b/src/finiteVolume/cfdTools/general/pressureControl/pressureControl.C
@@ -32,105 +32,171 @@ Foam::pressureControl::pressureControl
 (
     const volScalarField& p,
     const volScalarField& rho,
-    const dictionary& dict
+    const dictionary& dict,
+    const bool pRefRequired
 )
 :
-    refCell_(0),
+    refCell_(-1),
     refValue_(0),
     pMax_("pMax", dimPressure, 0),
     pMin_("pMin", dimPressure, GREAT)
 {
-    if (setRefCell(p, dict, refCell_, refValue_))
+    bool pLimits = false;
+
+    // Set the reference cell and value for closed domain simulations
+    if (pRefRequired && setRefCell(p, dict, refCell_, refValue_))
     {
+        pLimits = true;
+
         pMax_.value() = refValue_;
         pMin_.value() = refValue_;
     }
 
-    const volScalarField::Boundary& pbf = p.boundaryField();
-    const volScalarField::Boundary& rhobf = rho.boundaryField();
-
-    scalar rhoRefMax = -GREAT;
-    scalar rhoRefMin = GREAT;
-    bool rhoLimits = false;
-
-    forAll(pbf, patchi)
-    {
-        if (pbf[patchi].fixesValue())
-        {
-            rhoLimits = true;
-
-            pMax_.value() = max(pMax_.value(), max(pbf[patchi]));
-            pMin_.value() = min(pMin_.value(), min(pbf[patchi]));
-
-            rhoRefMax = max(rhoRefMax, max(rhobf[patchi]));
-            rhoRefMin = min(rhoRefMin, min(rhobf[patchi]));
-        }
-    }
-
-    if (dict.found("pMax"))
+    if (dict.found("pMax") && dict.found("pMin"))
     {
         pMax_.value() = readScalar(dict.lookup("pMax"));
-    }
-    else if (dict.found("pMaxFactor"))
-    {
-        const scalar pMaxFactor(readScalar(dict.lookup("pMaxFactor")));
-        pMax_ *= pMaxFactor;
-    }
-    else if (dict.found("rhoMax"))
-    {
-        // For backward-compatibility infer the pMax from rhoMax
-
-        IOWarningInFunction(dict)
-            << "'rhoMax' specified rather than 'pMax' or 'pMaxFactor'" << nl
-            << "    This is supported for backward-compatibility but "
-               "'pMax' or 'pMaxFactor' are more reliable." << endl;
-
-        if (!rhoLimits)
-        {
-            FatalIOErrorInFunction(dict)
-                << "'rhoMax' specified rather than 'pMaxFactor'" << nl
-                << "    but the corresponding reference density cannot"
-                   " be evaluated from the boundary conditions." << nl
-                << "Please specify 'pMaxFactor' rather than 'rhoMax'"
-                << exit(FatalError);
-        }
-
-        dimensionedScalar rhoMax("rhoMax", dimDensity, dict);
-
-        pMax_ *= max(rhoMax.value()/rhoRefMax, 1);
-    }
-
-    if (dict.found("pMin"))
-    {
         pMin_.value() = readScalar(dict.lookup("pMin"));
     }
-    else if (dict.found("pMinFactor"))
+    else
     {
-        const scalar pMinFactor(readScalar(dict.lookup("pMinFactor")));
-        pMin_ *= pMinFactor;
-    }
-    else if (dict.found("rhoMin"))
-    {
-        // For backward-compatibility infer the pMin from rhoMin
+        const volScalarField::Boundary& pbf = p.boundaryField();
+        const volScalarField::Boundary& rhobf = rho.boundaryField();
 
-        IOWarningInFunction(dict)
-            << "'rhoMin' specified rather than 'pMin' or 'pMinFactor'" << nl
-            << "    This is supported for backward-compatibility but"
-               "'pMin' or 'pMinFactor' are more reliable." << endl;
+        scalar rhoRefMax = -GREAT;
+        scalar rhoRefMin = GREAT;
+        bool rhoLimits = false;
 
-        if (!rhoLimits)
+        forAll(pbf, patchi)
         {
-            FatalIOErrorInFunction(dict)
-                << "'rhoMin' specified rather than 'pMinFactor'" << nl
-                << "    but the corresponding reference density cannot"
-                   " be evaluated from the boundary conditions." << nl
-                << "Please specify 'pMinFactor' rather than 'rhoMin'"
-                << exit(FatalError);
+            if (pbf[patchi].fixesValue())
+            {
+                pLimits = true;
+                rhoLimits = true;
+
+                pMax_.value() = max(pMax_.value(), max(pbf[patchi]));
+                pMin_.value() = min(pMin_.value(), min(pbf[patchi]));
+
+                rhoRefMax = max(rhoRefMax, max(rhobf[patchi]));
+                rhoRefMin = min(rhoRefMin, min(rhobf[patchi]));
+            }
         }
 
-        dimensionedScalar rhoMin("rhoMin", dimDensity, dict);
+        reduce(rhoLimits, andOp<bool>());
+        if (rhoLimits)
+        {
+            reduce(pMax_.value(), maxOp<scalar>());
+            reduce(pMin_.value(), minOp<scalar>());
 
-        pMin_ *= min(rhoMin.value()/rhoRefMin, 1);
+            reduce(rhoRefMax, maxOp<scalar>());
+            reduce(rhoRefMin, minOp<scalar>());
+        }
+
+        if (dict.found("pMax"))
+        {
+            pMax_.value() = readScalar(dict.lookup("pMax"));
+        }
+        else if (dict.found("pMaxFactor"))
+        {
+            if (!pLimits)
+            {
+                FatalIOErrorInFunction(dict)
+                    << "'pMaxFactor' specified rather than 'pMax'" << nl
+                    << "    but the corresponding reference pressure cannot"
+                       " be evaluated from the boundary conditions." << nl
+                    << "    Please specify 'pMax' rather than 'pMaxFactor'"
+                    << exit(FatalIOError);
+            }
+
+            const scalar pMaxFactor(readScalar(dict.lookup("pMaxFactor")));
+            pMax_ *= pMaxFactor;
+        }
+        else if (dict.found("rhoMax"))
+        {
+            // For backward-compatibility infer the pMax from rhoMax
+
+            IOWarningInFunction(dict)
+                 << "'rhoMax' specified rather than 'pMax' or 'pMaxFactor'"
+                 << nl
+                 << "    This is supported for backward-compatibility but "
+                    "'pMax' or 'pMaxFactor' are more reliable." << endl;
+
+            if (!pLimits)
+            {
+                FatalIOErrorInFunction(dict)
+                    << "'rhoMax' specified rather than 'pMax'" << nl
+                    << "    but the corresponding reference pressure cannot"
+                       " be evaluated from the boundary conditions." << nl
+                    << "    Please specify 'pMax' rather than 'rhoMax'"
+                    << exit(FatalIOError);
+            }
+
+            if (!rhoLimits)
+            {
+                FatalIOErrorInFunction(dict)
+                    << "'rhoMax' specified rather than 'pMaxFactor'" << nl
+                    << "    but the corresponding reference density cannot"
+                       " be evaluated from the boundary conditions." << nl
+                    << "    Please specify 'pMaxFactor' rather than 'rhoMax'"
+                    << exit(FatalIOError);
+            }
+
+            dimensionedScalar rhoMax("rhoMax", dimDensity, dict);
+
+            pMax_ *= max(rhoMax.value()/rhoRefMax, 1);
+        }
+
+        if (dict.found("pMin"))
+        {
+            pMin_.value() = readScalar(dict.lookup("pMin"));
+        }
+        else if (dict.found("pMinFactor"))
+        {
+            if (!pLimits)
+            {
+                FatalIOErrorInFunction(dict)
+                    << "'pMinFactor' specified rather than 'pMin'" << nl
+                    << "    but the corresponding reference pressure cannot"
+                       " be evaluated from the boundary conditions." << nl
+                    << "    Please specify 'pMin' rather than 'pMinFactor'"
+                    << exit(FatalIOError);
+            }
+
+            const scalar pMinFactor(readScalar(dict.lookup("pMinFactor")));
+            pMin_ *= pMinFactor;
+        }
+        else if (dict.found("rhoMin"))
+        {
+            // For backward-compatibility infer the pMin from rhoMin
+
+            IOWarningInFunction(dict)
+                << "'rhoMin' specified rather than 'pMin' or 'pMinFactor'" << nl
+                << "    This is supported for backward-compatibility but"
+                   "'pMin' or 'pMinFactor' are more reliable." << endl;
+
+            if (!pLimits)
+            {
+                FatalIOErrorInFunction(dict)
+                    << "'rhoMin' specified rather than 'pMin'" << nl
+                    << "    but the corresponding reference pressure cannot"
+                       " be evaluated from the boundary conditions." << nl
+                    << "    Please specify 'pMin' rather than 'rhoMin'"
+                    << exit(FatalIOError);
+            }
+
+            if (!rhoLimits)
+            {
+                FatalIOErrorInFunction(dict)
+                    << "'rhoMin' specified rather than 'pMinFactor'" << nl
+                    << "    but the corresponding reference density cannot"
+                       " be evaluated from the boundary conditions." << nl
+                    << "    Please specify 'pMinFactor' rather than 'rhoMin'"
+                    << exit(FatalIOError);
+            }
+
+            dimensionedScalar rhoMin("rhoMin", dimDensity, dict);
+
+            pMin_ *= min(rhoMin.value()/rhoRefMin, 1);
+        }
     }
 
     Info<< "pressureControl" << nl
diff --git a/src/finiteVolume/cfdTools/general/pressureControl/pressureControl.H b/src/finiteVolume/cfdTools/general/pressureControl/pressureControl.H
index 19bd053d43..61f510ac60 100644
--- a/src/finiteVolume/cfdTools/general/pressureControl/pressureControl.H
+++ b/src/finiteVolume/cfdTools/general/pressureControl/pressureControl.H
@@ -76,7 +76,8 @@ public:
         (
             const volScalarField& p,
             const volScalarField& rho,
-            const dictionary& dict
+            const dictionary& dict,
+            const bool pRefRequired = true
         );
 
 
diff --git a/tutorials/compressible/rhoPimpleDyMFoam/annularThermalMixer/system/fvSolution b/tutorials/compressible/rhoPimpleDyMFoam/annularThermalMixer/system/fvSolution
index 66786f3772..c771e4cfe3 100644
--- a/tutorials/compressible/rhoPimpleDyMFoam/annularThermalMixer/system/fvSolution
+++ b/tutorials/compressible/rhoPimpleDyMFoam/annularThermalMixer/system/fvSolution
@@ -63,8 +63,9 @@ PIMPLE
     nOuterCorrectors    3;
     nCorrectors         1;
     nNonOrthogonalCorrectors 0;
-    rhoMin          0.5;
-    rhoMax          2.0;
+
+    pMaxFactor          1.2;
+    pMinFactor          0.8;
 }
 
 relaxationFactors
diff --git a/tutorials/compressible/rhoPimpleFoam/LES/pitzDaily/system/fvSolution b/tutorials/compressible/rhoPimpleFoam/LES/pitzDaily/system/fvSolution
index 4eddb38721..b2dc2e505e 100644
--- a/tutorials/compressible/rhoPimpleFoam/LES/pitzDaily/system/fvSolution
+++ b/tutorials/compressible/rhoPimpleFoam/LES/pitzDaily/system/fvSolution
@@ -52,8 +52,9 @@ PIMPLE
     nOuterCorrectors 3;
     nCorrectors     1;
     nNonOrthogonalCorrectors 0;
-    rhoMin          0.5;
-    rhoMax          2.0;
+
+    pMinFactor      0.5;
+    pMaxFactor      2.0;
 }
 
 relaxationFactors
diff --git a/tutorials/compressible/rhoPimpleFoam/RAS/angledDuct/constant/thermophysicalProperties b/tutorials/compressible/rhoPimpleFoam/RAS/angledDuct/constant/thermophysicalProperties
index 32542b344e..7fc31fc96b 100644
--- a/tutorials/compressible/rhoPimpleFoam/RAS/angledDuct/constant/thermophysicalProperties
+++ b/tutorials/compressible/rhoPimpleFoam/RAS/angledDuct/constant/thermophysicalProperties
@@ -23,7 +23,7 @@ thermoType
     thermo          hConst;
     equationOfState perfectGas;
     specie          specie;
-    energy          sensibleEnthalpy;
+    energy          sensibleInternalEnergy;
 }
 
 mixture
diff --git a/tutorials/compressible/rhoPimpleFoam/RAS/angledDuct/system/fvSchemes b/tutorials/compressible/rhoPimpleFoam/RAS/angledDuct/system/fvSchemes
index 359d6ad3e1..c0522b83d6 100644
--- a/tutorials/compressible/rhoPimpleFoam/RAS/angledDuct/system/fvSchemes
+++ b/tutorials/compressible/rhoPimpleFoam/RAS/angledDuct/system/fvSchemes
@@ -30,8 +30,9 @@ divSchemes
     default         none;
     div(phi,U)      Gauss upwind;
     div(phid,p)     Gauss upwind;
+    div(phiv,p)     Gauss linear;
     div(phi,K)      Gauss linear;
-    div(phi,h)      Gauss upwind;
+    div(phi,e)      Gauss upwind;
     div(phi,k)      Gauss upwind;
     div(phi,epsilon) Gauss upwind;
     div(phi,R)      Gauss upwind;
diff --git a/tutorials/compressible/rhoPimpleFoam/RAS/angledDuct/system/fvSolution b/tutorials/compressible/rhoPimpleFoam/RAS/angledDuct/system/fvSolution
index ff24406781..3c7d9238f7 100644
--- a/tutorials/compressible/rhoPimpleFoam/RAS/angledDuct/system/fvSolution
+++ b/tutorials/compressible/rhoPimpleFoam/RAS/angledDuct/system/fvSolution
@@ -28,11 +28,10 @@ solvers
     pFinal
     {
         $p;
-        tolerance       1e-07;
         relTol          0;
     }
 
-    "(rho|U|h|k|epsilon|omega)"
+    "(rho|U|e|k|epsilon|omega)"
     {
         solver          smoothSolver;
         smoother        symGaussSeidel;
@@ -40,10 +39,9 @@ solvers
         relTol          0.1;
     }
 
-    "(rho|U|h|k|epsilon|omega)Final"
+    "(rho|U|e|k|epsilon|omega)Final"
     {
         $U;
-        tolerance       1e-06;
         relTol          0;
     }
 }
@@ -57,8 +55,8 @@ PIMPLE
     nNonOrthogonalCorrectors 0;
     consistent          yes;
 
-    rhoMin          0.4;
-    rhoMax          2.0;
+    pMaxFactor          1.5;
+    pMinFactor          0.9;
 
     residualControl
     {
@@ -76,13 +74,13 @@ relaxationFactors
 {
     fields
     {
-        "p.*"           0.9;
+        "p.*"           1;
         "rho.*"         1;
     }
     equations
     {
         "U.*"           0.9;
-        "h.*"           0.7;
+        "e.*"           0.7;
         "(k|epsilon|omega).*" 0.8;
     }
 }
diff --git a/tutorials/compressible/rhoPimpleFoam/RAS/angledDuctLTS/system/fvSolution b/tutorials/compressible/rhoPimpleFoam/RAS/angledDuctLTS/system/fvSolution
index 83e599258f..8b4068e0d6 100644
--- a/tutorials/compressible/rhoPimpleFoam/RAS/angledDuctLTS/system/fvSolution
+++ b/tutorials/compressible/rhoPimpleFoam/RAS/angledDuctLTS/system/fvSolution
@@ -56,8 +56,8 @@ PIMPLE
     nCorrectors       1;
     nNonOrthogonalCorrectors 0;
 
-    rhoMin            0.5;
-    rhoMax            2.0;
+    pMaxFactor          1.5;
+    pMinFactor          0.9;
 
     maxCo             0.2;
     rDeltaTSmoothingCoeff 0.1;
diff --git a/tutorials/compressible/rhoPimpleFoam/RAS/cavity/system/fvSolution b/tutorials/compressible/rhoPimpleFoam/RAS/cavity/system/fvSolution
index e7c64b955d..17ed812b9d 100644
--- a/tutorials/compressible/rhoPimpleFoam/RAS/cavity/system/fvSolution
+++ b/tutorials/compressible/rhoPimpleFoam/RAS/cavity/system/fvSolution
@@ -59,8 +59,9 @@ PIMPLE
     nOuterCorrectors 1;
     nCorrectors     2;
     nNonOrthogonalCorrectors 0;
-    rhoMin          0.5;
-    rhoMax          2.0;
+
+    pMax                1.2e5;
+    pMin                0.8e5;
 }
 
 
diff --git a/tutorials/compressible/rhoPimpleFoam/RAS/mixerVessel2D/system/fvSolution b/tutorials/compressible/rhoPimpleFoam/RAS/mixerVessel2D/system/fvSolution
index d4a55d13c4..02ea295883 100644
--- a/tutorials/compressible/rhoPimpleFoam/RAS/mixerVessel2D/system/fvSolution
+++ b/tutorials/compressible/rhoPimpleFoam/RAS/mixerVessel2D/system/fvSolution
@@ -84,14 +84,13 @@ solvers
 
 PIMPLE
 {
-    nOuterCorrectors 1;
-    nCorrectors     2;
-    nNonOrthogonalCorrectors 0;
-    momentumPredictor yes;
-    rhoMin          0.5;
-    rhoMax          2.0;
-    pRefCell        0;
-    pRefValue       1e5;
+    nOuterCorrectors    1;
+    nCorrectors         2;
+    nNonOrthogonalCorrectors  0;
+    momentumPredictor   yes;
+
+    pMax                1.2e5;
+    pMin                0.8e5;
 }
 
 relaxationFactors
diff --git a/tutorials/compressible/rhoPimpleFoam/laminar/helmholtzResonance/system/fvSolution b/tutorials/compressible/rhoPimpleFoam/laminar/helmholtzResonance/system/fvSolution
index 4eddb38721..b2dc2e505e 100644
--- a/tutorials/compressible/rhoPimpleFoam/laminar/helmholtzResonance/system/fvSolution
+++ b/tutorials/compressible/rhoPimpleFoam/laminar/helmholtzResonance/system/fvSolution
@@ -52,8 +52,9 @@ PIMPLE
     nOuterCorrectors 3;
     nCorrectors     1;
     nNonOrthogonalCorrectors 0;
-    rhoMin          0.5;
-    rhoMax          2.0;
+
+    pMinFactor      0.5;
+    pMaxFactor      2.0;
 }
 
 relaxationFactors