From 7b825d0817e14e1dab4c75a4b33d200b7cd77e41 Mon Sep 17 00:00:00 2001
From: Henry Weller <http://cfd.direct>
Date: Fri, 12 May 2017 11:04:38 +0100
Subject: [PATCH 1/6] rhoPimpleFoam: Updated transonic option to be consistent
 with sonicFoam

Improves stability on start-up and allows running at slightly larger time-steps.
---
 .../solvers/compressible/rhoPimpleFoam/pEqn.H |  15 +-
 .../compressible/rhoPimpleFoam/pcEqn.H        |  14 +-
 .../solvers/compressible/rhoSimpleFoam/pEqn.H | 197 +++++++++---------
 .../compressible/rhoSimpleFoam/pcEqn.H        |  14 +-
 .../general/pressureControl/pressureControl.C |  22 +-
 .../general/pressureControl/pressureControl.H |   4 +-
 6 files changed, 145 insertions(+), 121 deletions(-)

diff --git a/applications/solvers/compressible/rhoPimpleFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
index 5aab536239..9eecd01d99 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
@@ -1,3 +1,5 @@
+rho = thermo.rho();
+
 volScalarField rAU(1.0/UEqn.A());
 surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
 volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
@@ -10,7 +12,7 @@ if (pimple.nCorrPISO() <= 1)
 surfaceScalarField phiHbyA
 (
     "phiHbyA",
-    fvc::flux(rho*HbyA)
+    fvc::interpolate(rho)*fvc::flux(HbyA)
   + rhorAUf*fvc::ddtCorr(rho, U, phi)
 );
 
@@ -26,7 +28,8 @@ if (pimple.transonic())
         "phid",
         (fvc::interpolate(psi)/fvc::interpolate(rho))*phiHbyA
     );
-    phiHbyA -= fvc::interpolate(p)*phid;
+
+    phiHbyA -= fvc::interpolate(psi*p)*phiHbyA/fvc::interpolate(rho);
 
     while (pimple.correctNonOrthogonal())
     {
@@ -84,9 +87,11 @@ U.correctBoundaryConditions();
 fvOptions.correct(U);
 K = 0.5*magSqr(U);
 
-pressureControl.limit(p);
-p.correctBoundaryConditions();
-rho = thermo.rho();
+if (pressureControl.limit(p))
+{
+    p.correctBoundaryConditions();
+    rho = thermo.rho();
+}
 
 if (!pimple.transonic())
 {
diff --git a/applications/solvers/compressible/rhoPimpleFoam/pcEqn.H b/applications/solvers/compressible/rhoPimpleFoam/pcEqn.H
index afbc2851e4..34657a1892 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/pcEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/pcEqn.H
@@ -1,3 +1,5 @@
+rho = thermo.rho();
+
 volScalarField rAU(1.0/UEqn.A());
 volScalarField rAtU(1.0/(1.0/rAU - UEqn.H1()));
 volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
@@ -11,7 +13,7 @@ surfaceScalarField phiHbyA
 (
     "phiHbyA",
     (
-        fvc::flux(rho*HbyA)
+        fvc::interpolate(rho)*fvc::flux(HbyA)
       + fvc::interpolate(rho*rAU)*fvc::ddtCorr(rho, U, phi)
     )
 );
@@ -33,7 +35,7 @@ if (pimple.transonic())
 
     phiHbyA +=
         fvc::interpolate(rho*(rAtU - rAU))*fvc::snGrad(p)*mesh.magSf()
-      - fvc::interpolate(p)*phid;
+      - fvc::interpolate(psi*p)*phiHbyA/fvc::interpolate(rho);
 
     HbyA -= (rAU - rAtU)*fvc::grad(p);
 
@@ -96,9 +98,11 @@ U.correctBoundaryConditions();
 fvOptions.correct(U);
 K = 0.5*magSqr(U);
 
-pressureControl.limit(p);
-p.correctBoundaryConditions();
-rho = thermo.rho();
+if (pressureControl.limit(p))
+{
+    p.correctBoundaryConditions();
+    rho = thermo.rho();
+}
 
 if (!pimple.transonic())
 {
diff --git a/applications/solvers/compressible/rhoSimpleFoam/pEqn.H b/applications/solvers/compressible/rhoSimpleFoam/pEqn.H
index c7edbb9e87..0f290ed273 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/pEqn.H
@@ -1,109 +1,110 @@
+volScalarField rAU(1.0/UEqn.A());
+surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
+volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
+tUEqn.clear();
+
+bool closedVolume = false;
+
+surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(rho)*fvc::flux(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())
 {
-    volScalarField rAU(1.0/UEqn.A());
-    surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
-    volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
-    tUEqn.clear();
+    surfaceScalarField phid
+    (
+        "phid",
+        (fvc::interpolate(psi)/fvc::interpolate(rho))*phiHbyA
+    );
 
-    bool closedVolume = false;
+    phiHbyA -= fvc::interpolate(psi*p)*phiHbyA/fvc::interpolate(rho);
 
-    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())
+    while (simple.correctNonOrthogonal())
     {
-        surfaceScalarField phid
+        fvScalarMatrix pEqn
         (
-            "phid",
-            (fvc::interpolate(psi)/fvc::interpolate(rho))*phiHbyA
+            fvc::div(phiHbyA)
+          + fvm::div(phid, p)
+          - fvm::laplacian(rhorAUf, p)
+          ==
+            fvOptions(psi, p, rho.name())
         );
-        phiHbyA -= fvc::interpolate(p)*phid;
 
-        while (simple.correctNonOrthogonal())
+        // Relax the pressure equation to ensure diagonal-dominance
+        pEqn.relax();
+
+        pEqn.setReference
+        (
+            pressureControl.refCell(),
+            pressureControl.refValue()
+        );
+
+        pEqn.solve();
+
+        if (simple.finalNonOrthogonalIter())
         {
-            fvScalarMatrix pEqn
-            (
-                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.setReference
-            (
-                pressureControl.refCell(),
-                pressureControl.refValue()
-            );
-
-            pEqn.solve();
-
-            if (simple.finalNonOrthogonalIter())
-            {
-                phi = phiHbyA + pEqn.flux();
-            }
+            phi = phiHbyA + pEqn.flux();
         }
     }
-    else
-    {
-        closedVolume = adjustPhi(phiHbyA, U, p);
-
-        while (simple.correctNonOrthogonal())
-        {
-            fvScalarMatrix pEqn
-            (
-                fvc::div(phiHbyA)
-              - fvm::laplacian(rhorAUf, p)
-              ==
-                fvOptions(psi, p, rho.name())
-            );
-
-            pEqn.setReference
-            (
-                pressureControl.refCell(),
-                pressureControl.refValue()
-            );
-
-            pEqn.solve();
-
-            if (simple.finalNonOrthogonalIter())
-            {
-                phi = phiHbyA + pEqn.flux();
-            }
-        }
-    }
-
-
-    #include "incompressible/continuityErrs.H"
-
-    // Explicitly relax pressure for momentum corrector
-    p.relax();
-
-    U = HbyA - rAU*fvc::grad(p);
-    U.correctBoundaryConditions();
-    fvOptions.correct(U);
-
-    pressureControl.limit(p);
-
-    // For closed-volume cases adjust the pressure and density levels
-    // to obey overall mass continuity
-    if (closedVolume)
-    {
-        p += (initialMass - fvc::domainIntegrate(psi*p))
-            /fvc::domainIntegrate(psi);
-    }
-
-    p.correctBoundaryConditions();
-
-    rho = thermo.rho();
-
-    if (!simple.transonic())
-    {
-        rho.relax();
-    }
+}
+else
+{
+    closedVolume = adjustPhi(phiHbyA, U, p);
+
+    while (simple.correctNonOrthogonal())
+    {
+        fvScalarMatrix pEqn
+        (
+            fvc::div(phiHbyA)
+          - fvm::laplacian(rhorAUf, p)
+          ==
+            fvOptions(psi, p, rho.name())
+        );
+
+        pEqn.setReference
+        (
+            pressureControl.refCell(),
+            pressureControl.refValue()
+        );
+
+        pEqn.solve();
+
+        if (simple.finalNonOrthogonalIter())
+        {
+            phi = phiHbyA + pEqn.flux();
+        }
+    }
+}
+
+#include "incompressible/continuityErrs.H"
+
+// Explicitly relax pressure for momentum corrector
+p.relax();
+
+U = HbyA - rAU*fvc::grad(p);
+U.correctBoundaryConditions();
+fvOptions.correct(U);
+
+bool pLimited = pressureControl.limit(p);
+
+// For closed-volume cases adjust the pressure and density levels
+// to obey overall mass continuity
+if (closedVolume)
+{
+    p += (initialMass - fvc::domainIntegrate(psi*p))
+        /fvc::domainIntegrate(psi);
+}
+
+if (pLimited || closedVolume)
+{
+    p.correctBoundaryConditions();
+}
+
+rho = thermo.rho();
+
+if (!simple.transonic())
+{
+    rho.relax();
 }
diff --git a/applications/solvers/compressible/rhoSimpleFoam/pcEqn.H b/applications/solvers/compressible/rhoSimpleFoam/pcEqn.H
index c7dc0f864d..51560aab5f 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/pcEqn.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/pcEqn.H
@@ -1,3 +1,5 @@
+rho = thermo.rho();
+
 volScalarField rAU(1.0/UEqn.A());
 volScalarField rAtU(1.0/(1.0/rAU - UEqn.H1()));
 volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
@@ -5,7 +7,7 @@ tUEqn.clear();
 
 bool closedVolume = false;
 
-surfaceScalarField phiHbyA("phiHbyA", fvc::flux(rho*HbyA));
+surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(rho)*fvc::flux(HbyA));
 MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
 
 volScalarField rhorAtU("rhorAtU", rho*rAtU);
@@ -23,7 +25,7 @@ if (simple.transonic())
 
     phiHbyA +=
         fvc::interpolate(rho*(rAtU - rAU))*fvc::snGrad(p)*mesh.magSf()
-      - fvc::interpolate(p)*phid;
+      - fvc::interpolate(psi*p)*phiHbyA/fvc::interpolate(rho);
 
     HbyA -= (rAU - rAtU)*fvc::grad(p);
 
@@ -98,7 +100,7 @@ U = HbyA - rAtU*fvc::grad(p);
 U.correctBoundaryConditions();
 fvOptions.correct(U);
 
-pressureControl.limit(p);
+bool pLimited = pressureControl.limit(p);
 
 // For closed-volume cases adjust the pressure and density levels
 // to obey overall mass continuity
@@ -108,9 +110,11 @@ if (closedVolume)
         /fvc::domainIntegrate(psi);
 }
 
-p.correctBoundaryConditions();
+if (pLimited || closedVolume)
+{
+    p.correctBoundaryConditions();
+}
 
-// Recalculate density from the relaxed pressure
 rho = thermo.rho();
 
 if (!simple.transonic())
diff --git a/src/finiteVolume/cfdTools/general/pressureControl/pressureControl.C b/src/finiteVolume/cfdTools/general/pressureControl/pressureControl.C
index aca22e7699..abd71e60ac 100644
--- a/src/finiteVolume/cfdTools/general/pressureControl/pressureControl.C
+++ b/src/finiteVolume/cfdTools/general/pressureControl/pressureControl.C
@@ -207,14 +207,24 @@ Foam::pressureControl::pressureControl
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
-void Foam::pressureControl::limit(volScalarField& p) const
+bool Foam::pressureControl::limit(volScalarField& p) const
 {
-    Info<< "pressureControl: p max/min "
-        << max(p).value() << " "
-        << min(p).value() << endl;
+    scalar pMax = max(p).value();
+    scalar pMin = min(p).value();
 
-    p = max(p, pMin_);
-    p = min(p, pMax_);
+    if (pMax > pMax_.value() || pMin < pMin_.value())
+    {
+        Info<< "pressureControl: p max/min " << pMax << " " << pMin << endl;
+
+        p = max(p, pMin_);
+        p = min(p, pMax_);
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
 }
 
 
diff --git a/src/finiteVolume/cfdTools/general/pressureControl/pressureControl.H b/src/finiteVolume/cfdTools/general/pressureControl/pressureControl.H
index 61f510ac60..0bb49c99a5 100644
--- a/src/finiteVolume/cfdTools/general/pressureControl/pressureControl.H
+++ b/src/finiteVolume/cfdTools/general/pressureControl/pressureControl.H
@@ -89,8 +89,8 @@ public:
         //- Return the pressure reference level
         inline scalar refValue() const;
 
-        //- Limit the pressure
-        void limit(volScalarField& p) const;
+        //- Limit the pressure if necessary and return true if so
+        bool limit(volScalarField& p) const;
 };
 
 

From b83af3b085abde7ce80d8213b609666c1c699804 Mon Sep 17 00:00:00 2001
From: Henry Weller <http://cfd.direct>
Date: Fri, 12 May 2017 18:34:00 +0100
Subject: [PATCH 2/6] rhoPimpleDyMFoam: Updated transonic formulation for
 consistency with sonicFoam

---
 .../rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H           | 13 +++++++++----
 1 file changed, 9 insertions(+), 4 deletions(-)

diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H
index 0f85562618..1540f601c1 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H
@@ -1,3 +1,5 @@
+rho = thermo.rho();
+
 volScalarField rAU(1.0/UEqn.A());
 surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
 volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
@@ -27,7 +29,8 @@ if (pimple.transonic())
         "phid",
         (fvc::interpolate(psi)/fvc::interpolate(rho))*phiHbyA
     );
-    phiHbyA -= fvc::interpolate(p)*phid;
+
+    phiHbyA -= fvc::interpolate(psi*p)*phiHbyA/fvc::interpolate(rho);
 
     while (pimple.correctNonOrthogonal())
     {
@@ -86,9 +89,11 @@ U.correctBoundaryConditions();
 fvOptions.correct(U);
 K = 0.5*magSqr(U);
 
-pressureControl.limit(p);
-p.correctBoundaryConditions();
-rho = thermo.rho();
+if (pressureControl.limit(p))
+{
+    p.correctBoundaryConditions();
+    rho = thermo.rho();
+}
 
 if (!pimple.transonic())
 {

From 39476bde1c8feefcd4e37c55663bca8fa43e422d Mon Sep 17 00:00:00 2001
From: Henry Weller <http://cfd.direct>
Date: Wed, 17 May 2017 17:05:43 +0100
Subject: [PATCH 3/6] GIT: Resolve conflict associated with cherry-pick of
 Foundation commit 79ff91350

79ff91350 - rhoPimpleFoam: Improved support for compressible liquids
(2017-05-17 17:05:43 +0100) <Henry Weller>
---
 .../combustion/reactingFoam/createFields.H    |  22 +--
 .../solvers/combustion/reactingFoam/pEqn.H    |  17 +-
 .../solvers/combustion/reactingFoam/pcEqn.H   |  21 +--
 .../combustion/reactingFoam/reactingFoam.C    |   5 +-
 .../rhoReactingBuoyantFoam/pEqn.H             | 130 +++++++------
 .../rhoReactingFoam/createFields.H            |   2 +
 .../reactingFoam/rhoReactingFoam/pEqn.H       | 109 -----------
 .../rhoReactingFoam/rhoReactingFoam.C         |  12 +-
 .../compressible/rhoPimpleFoam/createFields.H |   4 +-
 .../solvers/compressible/rhoPimpleFoam/pEqn.H |  52 +++---
 .../compressible/rhoPimpleFoam/pcEqn.H        |  52 +++---
 .../rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H     |   9 +-
 .../rhoPimpleFoam/rhoPimpleFoam.C             |   2 +
 .../heatTransfer/buoyantPimpleFoam/pEqn.H     | 173 +++++++----------
 .../chtMultiRegionFoam/fluid/pEqn.H           | 176 ++++++++----------
 .../coalChemistryFoam/coalChemistryFoam.C     |   3 +-
 .../lagrangian/reactingParcelFoam/pEqn.H      | 128 +++++++------
 .../simpleReactingParcelFoam/pEqn.H           |  96 +++++-----
 .../compressibleInterDyMFoam/pEqn.H           |   4 +-
 .../compressibleInterFoam/createFields.H      |   4 +-
 .../multiphase/compressibleInterFoam/pEqn.H   |   4 +-
 .../pimpleControl/pimpleControl.C             |   4 +-
 .../pimpleControl/pimpleControl.H             |  10 +-
 .../pimpleControl/pimpleControlI.H            |   8 +-
 .../basic/fluidThermo/fluidThermo.H           |   6 +-
 .../basic/psiThermo/psiThermo.C               |   6 +-
 .../basic/psiThermo/psiThermo.H               |   7 +-
 .../basic/rhoThermo/rhoThermo.C               |   8 +-
 .../basic/rhoThermo/rhoThermo.H               |   6 +-
 .../RAS/angledDuct/system/fvSolution          |   1 +
 .../rhoPimpleFoam/RAS/squareBendLiq/0/T       |  42 +++++
 .../rhoPimpleFoam/RAS/squareBendLiq/0/U       |  43 +++++
 .../rhoPimpleFoam/RAS/squareBendLiq/0/alphat  |  45 +++++
 .../rhoPimpleFoam/RAS/squareBendLiq/0/epsilon |  49 +++++
 .../rhoPimpleFoam/RAS/squareBendLiq/0/k       |  46 +++++
 .../rhoPimpleFoam/RAS/squareBendLiq/0/nut     |  47 +++++
 .../rhoPimpleFoam/RAS/squareBendLiq/0/p       |  41 ++++
 .../constant/thermophysicalProperties         |  31 +++
 .../constant/turbulenceProperties             |  30 +++
 .../RAS/squareBendLiq/system/blockMeshDict    | 127 +++++++++++++
 .../RAS/squareBendLiq/system/controlDict      |  51 +++++
 .../RAS/squareBendLiq/system/decomposeParDict |  45 +++++
 .../RAS/squareBendLiq/system/fvSchemes        |  60 ++++++
 .../RAS/squareBendLiq/system/fvSolution       |  63 +++++++
 44 files changed, 1194 insertions(+), 607 deletions(-)
 delete mode 100644 applications/solvers/combustion/reactingFoam/rhoReactingFoam/pEqn.H
 create mode 100644 tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/0/T
 create mode 100644 tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/0/U
 create mode 100644 tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/0/alphat
 create mode 100644 tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/0/epsilon
 create mode 100644 tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/0/k
 create mode 100644 tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/0/nut
 create mode 100644 tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/0/p
 create mode 100644 tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/constant/thermophysicalProperties
 create mode 100644 tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/constant/turbulenceProperties
 create mode 100644 tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/system/blockMeshDict
 create mode 100644 tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/system/controlDict
 create mode 100644 tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/system/decomposeParDict
 create mode 100644 tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/system/fvSchemes
 create mode 100644 tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/system/fvSolution

diff --git a/applications/solvers/combustion/reactingFoam/createFields.H b/applications/solvers/combustion/reactingFoam/createFields.H
index 6fcf83624c..bc09358193 100644
--- a/applications/solvers/combustion/reactingFoam/createFields.H
+++ b/applications/solvers/combustion/reactingFoam/createFields.H
@@ -52,27 +52,7 @@ volScalarField& p = thermo.p();
 
 #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);
 
 mesh.setFluxRequired(p.name());
 
diff --git a/applications/solvers/combustion/reactingFoam/pEqn.H b/applications/solvers/combustion/reactingFoam/pEqn.H
index ac7107acf0..81ab6424b4 100644
--- a/applications/solvers/combustion/reactingFoam/pEqn.H
+++ b/applications/solvers/combustion/reactingFoam/pEqn.H
@@ -1,7 +1,4 @@
 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));
@@ -87,19 +84,17 @@ 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);
 
+if (pressureControl.limit(p))
+{
+    p.correctBoundaryConditions();
+    rho = thermo.rho();
+}
+
 if (thermo.dpdt())
 {
     dpdt = fvc::ddt(p);
diff --git a/applications/solvers/combustion/reactingFoam/pcEqn.H b/applications/solvers/combustion/reactingFoam/pcEqn.H
index 713f443fc5..a1564b0f20 100644
--- a/applications/solvers/combustion/reactingFoam/pcEqn.H
+++ b/applications/solvers/combustion/reactingFoam/pcEqn.H
@@ -1,7 +1,4 @@
 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()));
@@ -109,19 +106,13 @@ U.correctBoundaryConditions();
 fvOptions.correct(U);
 K = 0.5*magSqr(U);
 
+if (pressureControl.limit(p))
+{
+    p.correctBoundaryConditions();
+    rho = thermo.rho();
+}
+
 if (thermo.dpdt())
 {
     dpdt = fvc::ddt(p);
 }
-
-// Recalculate density from the relaxed pressure
-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;
diff --git a/applications/solvers/combustion/reactingFoam/reactingFoam.C b/applications/solvers/combustion/reactingFoam/reactingFoam.C
index a06eef279c..d3e1b19ad3 100644
--- a/applications/solvers/combustion/reactingFoam/reactingFoam.C
+++ b/applications/solvers/combustion/reactingFoam/reactingFoam.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
@@ -37,6 +37,7 @@ Description
 #include "psiCombustionModel.H"
 #include "multivariateScheme.H"
 #include "pimpleControl.H"
+#include "pressureControl.H"
 #include "fvOptions.H"
 #include "localEulerDdtScheme.H"
 #include "fvcSmooth.H"
@@ -114,6 +115,8 @@ int main(int argc, char *argv[])
             }
         }
 
+        rho = thermo.rho();
+
         runTime.write();
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
diff --git a/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/pEqn.H b/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/pEqn.H
index 926e69f7c1..04296ddfc2 100644
--- a/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/pEqn.H
+++ b/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/pEqn.H
@@ -1,76 +1,74 @@
+rho = thermo.rho();
+
+// Thermodynamic density needs to be updated by psi*d(p) after the
+// pressure solution
+const volScalarField psip0(psi*p);
+
+volScalarField rAU(1.0/UEqn.A());
+surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
+volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
+
+surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
+
+surfaceScalarField phiHbyA
+(
+    "phiHbyA",
+    (
+        fvc::flux(rho*HbyA)
+      + rhorAUf*fvc::ddtCorr(rho, U, phi)
+    )
+  + phig
+);
+
+MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
+
+// Update the pressure BCs to ensure flux consistency
+constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
+
+fvScalarMatrix p_rghDDtEqn
+(
+    fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
+  + fvc::div(phiHbyA)
+  ==
+    fvOptions(psi, p_rgh, rho.name())
+);
+
+while (pimple.correctNonOrthogonal())
 {
-    rho = thermo.rho();
-
-    // Thermodynamic density needs to be updated by psi*d(p) after the
-    // pressure solution - done in 2 parts. Part 1:
-    thermo.rho() -= psi*p;
-
-    volScalarField rAU(1.0/UEqn.A());
-    surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
-    volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
-
-    surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
-
-    surfaceScalarField phiHbyA
+    fvScalarMatrix p_rghEqn
     (
-        "phiHbyA",
-        (
-            fvc::flux(rho*HbyA)
-          + rhorAUf*fvc::ddtCorr(rho, U, phi)
-        )
-      + phig
+        p_rghDDtEqn
+      - fvm::laplacian(rhorAUf, p_rgh)
     );
 
-    MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
+    p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
 
-    // Update the pressure BCs to ensure flux consistency
-    constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
-
-    fvScalarMatrix p_rghDDtEqn
-    (
-        fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
-      + fvc::div(phiHbyA)
-      ==
-        fvOptions(psi, p_rgh, rho.name())
-    );
-
-    while (pimple.correctNonOrthogonal())
+    if (pimple.finalNonOrthogonalIter())
     {
-        fvScalarMatrix p_rghEqn
-        (
-            p_rghDDtEqn
-          - fvm::laplacian(rhorAUf, p_rgh)
-        );
+        // Calculate the conservative fluxes
+        phi = phiHbyA + p_rghEqn.flux();
 
-        p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
+        // Explicitly relax pressure for momentum corrector
+        p_rgh.relax();
 
-        if (pimple.finalNonOrthogonalIter())
-        {
-            // Calculate the conservative fluxes
-            phi = phiHbyA + p_rghEqn.flux();
-
-            // Explicitly relax pressure for momentum corrector
-            p_rgh.relax();
-
-            // Correct the momentum source with the pressure gradient flux
-            // calculated from the relaxed pressure
-            U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rhorAUf);
-            U.correctBoundaryConditions();
-            fvOptions.correct(U);
-            K = 0.5*magSqr(U);
-        }
+        // Correct the momentum source with the pressure gradient flux
+        // calculated from the relaxed pressure
+        U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rhorAUf);
+        U.correctBoundaryConditions();
+        fvOptions.correct(U);
+        K = 0.5*magSqr(U);
     }
-
-    p = p_rgh + rho*gh;
-
-    // Second part of thermodynamic density update
-    thermo.rho() += psi*p;
-
-    if (thermo.dpdt())
-    {
-        dpdt = fvc::ddt(p);
-    }
-
-    #include "rhoEqn.H"
-    #include "compressibleContinuityErrs.H"
 }
+
+p = p_rgh + rho*gh;
+
+// Thermodynamic density update
+thermo.correctRho(psi*p - psip0);
+
+if (thermo.dpdt())
+{
+    dpdt = fvc::ddt(p);
+}
+
+#include "rhoEqn.H"
+#include "compressibleContinuityErrs.H"
diff --git a/applications/solvers/combustion/reactingFoam/rhoReactingFoam/createFields.H b/applications/solvers/combustion/reactingFoam/rhoReactingFoam/createFields.H
index 687dffcc13..e56035f24a 100644
--- a/applications/solvers/combustion/reactingFoam/rhoReactingFoam/createFields.H
+++ b/applications/solvers/combustion/reactingFoam/rhoReactingFoam/createFields.H
@@ -52,6 +52,8 @@ volScalarField& p = thermo.p();
 
 #include "compressibleCreatePhi.H"
 
+pressureControl pressureControl(p, rho, pimple.dict(), false);
+
 mesh.setFluxRequired(p.name());
 
 
diff --git a/applications/solvers/combustion/reactingFoam/rhoReactingFoam/pEqn.H b/applications/solvers/combustion/reactingFoam/rhoReactingFoam/pEqn.H
deleted file mode 100644
index eff014cf43..0000000000
--- a/applications/solvers/combustion/reactingFoam/rhoReactingFoam/pEqn.H
+++ /dev/null
@@ -1,109 +0,0 @@
-{
-    rho = thermo.rho();
-
-    // Thermodynamic density needs to be updated by psi*d(p) after the
-    // pressure solution - done in 2 parts. Part 1:
-    thermo.rho() -= psi*p;
-
-    volScalarField rAU(1.0/UEqn.A());
-    surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
-    volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
-
-    if (pimple.transonic())
-    {
-        surfaceScalarField phiHbyA
-        (
-            "phiHbyA",
-            (
-                fvc::flux(HbyA)
-              + rhorAUf*fvc::ddtCorr(rho, U, phi)/fvc::interpolate(rho)
-            )
-        );
-
-        MRF.makeRelative(phiHbyA);
-
-        surfaceScalarField phid("phid", fvc::interpolate(thermo.psi())*phiHbyA);
-
-        phiHbyA *= fvc::interpolate(rho);
-
-        fvScalarMatrix pDDtEqn
-        (
-            fvc::ddt(rho) + fvc::div(phiHbyA)
-          + correction(psi*fvm::ddt(p) + fvm::div(phid, p))
-        );
-
-        while (pimple.correctNonOrthogonal())
-        {
-            fvScalarMatrix pEqn
-            (
-                pDDtEqn
-              - fvm::laplacian(rhorAUf, p)
-             ==
-                fvOptions(psi, p, rho.name())
-            );
-
-            pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
-
-            if (pimple.finalNonOrthogonalIter())
-            {
-                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);
-
-        fvScalarMatrix pDDtEqn
-        (
-            fvc::ddt(rho) + psi*correction(fvm::ddt(p))
-          + fvc::div(phiHbyA)
-         ==
-            fvOptions(psi, p, rho.name())
-        );
-
-        while (pimple.correctNonOrthogonal())
-        {
-            fvScalarMatrix pEqn
-            (
-                pDDtEqn
-              - fvm::laplacian(rhorAUf, p)
-            );
-
-            pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
-
-            if (pimple.finalNonOrthogonalIter())
-            {
-                phi = phiHbyA + pEqn.flux();
-            }
-        }
-    }
-
-    // Second part of thermodynamic density update
-    thermo.rho() += psi*p;
-
-    #include "rhoEqn.H"
-    #include "compressibleContinuityErrs.H"
-
-    U = HbyA - rAU*fvc::grad(p);
-    U.correctBoundaryConditions();
-    fvOptions.correct(U);
-    K = 0.5*magSqr(U);
-
-    if (thermo.dpdt())
-    {
-        dpdt = fvc::ddt(p);
-    }
-}
diff --git a/applications/solvers/combustion/reactingFoam/rhoReactingFoam/rhoReactingFoam.C b/applications/solvers/combustion/reactingFoam/rhoReactingFoam/rhoReactingFoam.C
index 55de274d83..5be897c6c7 100644
--- a/applications/solvers/combustion/reactingFoam/rhoReactingFoam/rhoReactingFoam.C
+++ b/applications/solvers/combustion/reactingFoam/rhoReactingFoam/rhoReactingFoam.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
@@ -38,6 +38,7 @@ Description
 #include "turbulentFluidThermoModel.H"
 #include "multivariateScheme.H"
 #include "pimpleControl.H"
+#include "pressureControl.H"
 #include "fvOptions.H"
 #include "localEulerDdtScheme.H"
 #include "fvcSmooth.H"
@@ -100,7 +101,14 @@ int main(int argc, char *argv[])
             // --- Pressure corrector loop
             while (pimple.correct())
             {
-                #include "pEqn.H"
+                if (pimple.consistent())
+                {
+                    #include "../../../compressible/rhoPimpleFoam/pcEqn.H"
+                }
+                else
+                {
+                    #include "../../../compressible/rhoPimpleFoam/pEqn.H"
+                }
             }
 
             if (pimple.turbCorr())
diff --git a/applications/solvers/compressible/rhoPimpleFoam/createFields.H b/applications/solvers/compressible/rhoPimpleFoam/createFields.H
index 84f70f2cd0..3ec620dc22 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/createFields.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/createFields.H
@@ -42,6 +42,8 @@ volVectorField U
 
 pressureControl pressureControl(p, rho, pimple.dict(), false);
 
+mesh.setFluxRequired(p.name());
+
 Info<< "Creating turbulence model\n" << endl;
 autoPtr<compressible::turbulenceModel> turbulence
 (
@@ -54,8 +56,6 @@ autoPtr<compressible::turbulenceModel> turbulence
     )
 );
 
-mesh.setFluxRequired(p.name());
-
 Info<< "Creating field dpdt\n" << endl;
 volScalarField dpdt
 (
diff --git a/applications/solvers/compressible/rhoPimpleFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
index 9eecd01d99..83eeb8fc68 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
@@ -1,4 +1,11 @@
-rho = thermo.rho();
+if (!pimple.SIMPLErho())
+{
+    rho = thermo.rho();
+}
+
+// Thermodynamic density needs to be updated by psi*d(p) after the
+// pressure solution
+const volScalarField psip0(psi*p);
 
 volScalarField rAU(1.0/UEqn.A());
 surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
@@ -31,17 +38,17 @@ if (pimple.transonic())
 
     phiHbyA -= fvc::interpolate(psi*p)*phiHbyA/fvc::interpolate(rho);
 
+    fvScalarMatrix pDDtEqn
+    (
+        fvc::ddt(rho) + psi*correction(fvm::ddt(p))
+      + fvc::div(phiHbyA) + fvm::div(phid, p)
+     ==
+        fvOptions(psi, p, rho.name())
+    );
+
     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())
-        );
+        fvScalarMatrix pEqn(pDDtEqn - fvm::laplacian(rhorAUf, p));
 
         // Relax the pressure equation to ensure diagonal-dominance
         pEqn.relax();
@@ -56,16 +63,17 @@ if (pimple.transonic())
 }
 else
 {
+    fvScalarMatrix pDDtEqn
+    (
+        fvc::ddt(rho) + psi*correction(fvm::ddt(p))
+      + fvc::div(phiHbyA)
+     ==
+        fvOptions(psi, p, rho.name())
+    );
+
     while (pimple.correctNonOrthogonal())
     {
-        fvScalarMatrix pEqn
-        (
-            fvm::ddt(psi, p)
-          + fvc::div(phiHbyA)
-          - fvm::laplacian(rhorAUf, p)
-         ==
-            fvOptions(psi, p, rho.name())
-        );
+        fvScalarMatrix pEqn(pDDtEqn - fvm::laplacian(rhorAUf, p));
 
         pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
@@ -76,6 +84,9 @@ else
     }
 }
 
+// Thermodynamic density update
+thermo.correctRho(psi*p - psip0);
+
 #include "rhoEqn.H"
 #include "compressibleContinuityErrs.H"
 
@@ -92,10 +103,9 @@ if (pressureControl.limit(p))
     p.correctBoundaryConditions();
     rho = thermo.rho();
 }
-
-if (!pimple.transonic())
+else if (pimple.SIMPLErho())
 {
-    rho.relax();
+    rho = thermo.rho();
 }
 
 if (thermo.dpdt())
diff --git a/applications/solvers/compressible/rhoPimpleFoam/pcEqn.H b/applications/solvers/compressible/rhoPimpleFoam/pcEqn.H
index 34657a1892..e9b849bc05 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/pcEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/pcEqn.H
@@ -1,4 +1,11 @@
-rho = thermo.rho();
+if (!pimple.SIMPLErho())
+{
+    rho = thermo.rho();
+}
+
+// Thermodynamic density needs to be updated by psi*d(p) after the
+// pressure solution
+const volScalarField psip0(psi*p);
 
 volScalarField rAU(1.0/UEqn.A());
 volScalarField rAtU(1.0/(1.0/rAU - UEqn.H1()));
@@ -39,17 +46,17 @@ if (pimple.transonic())
 
     HbyA -= (rAU - rAtU)*fvc::grad(p);
 
+    fvScalarMatrix pDDtEqn
+    (
+        fvc::ddt(rho) + psi*correction(fvm::ddt(p))
+      + fvc::div(phiHbyA) + fvm::div(phid, p)
+     ==
+        fvOptions(psi, p, rho.name())
+    );
+
     while (pimple.correctNonOrthogonal())
     {
-        fvScalarMatrix pEqn
-        (
-            fvm::ddt(psi, p)
-          + fvc::div(phiHbyA)
-          + fvm::div(phid, p)
-          - fvm::laplacian(rhorAtU, p)
-         ==
-            fvOptions(psi, p, rho.name())
-        );
+        fvScalarMatrix pEqn(pDDtEqn - fvm::laplacian(rhorAtU, p));
 
         // Relax the pressure equation to ensure diagonal-dominance
         pEqn.relax();
@@ -67,16 +74,17 @@ else
     phiHbyA += fvc::interpolate(rho*(rAtU - rAU))*fvc::snGrad(p)*mesh.magSf();
     HbyA -= (rAU - rAtU)*fvc::grad(p);
 
+    fvScalarMatrix pDDtEqn
+    (
+        fvc::ddt(rho) + psi*correction(fvm::ddt(p))
+      + fvc::div(phiHbyA)
+     ==
+        fvOptions(psi, p, rho.name())
+    );
+
     while (pimple.correctNonOrthogonal())
     {
-        fvScalarMatrix pEqn
-        (
-            fvm::ddt(psi, p)
-          + fvc::div(phiHbyA)
-          - fvm::laplacian(rhorAtU, p)
-         ==
-            fvOptions(psi, p, rho.name())
-        );
+        fvScalarMatrix pEqn(pDDtEqn - fvm::laplacian(rhorAtU, p));
 
         pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
@@ -87,6 +95,9 @@ else
     }
 }
 
+// Thermodynamic density update
+thermo.correctRho(psi*p - psip0);
+
 #include "rhoEqn.H"
 #include "compressibleContinuityErrs.H"
 
@@ -103,10 +114,9 @@ if (pressureControl.limit(p))
     p.correctBoundaryConditions();
     rho = thermo.rho();
 }
-
-if (!pimple.transonic())
+else if (pimple.SIMPLErho())
 {
-    rho.relax();
+    rho = thermo.rho();
 }
 
 if (thermo.dpdt())
diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H
index 1540f601c1..ed802ba6d2 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H
@@ -36,7 +36,7 @@ if (pimple.transonic())
     {
         fvScalarMatrix pEqn
         (
-            fvm::ddt(psi, p)
+            fvc::ddt(rho) + psi*correction(fvm::ddt(p))
           + fvc::div(phiHbyA)
           + fvm::div(phid, p)
           - fvm::laplacian(rhorAUf, p)
@@ -62,7 +62,7 @@ else
         // Pressure corrector
         fvScalarMatrix pEqn
         (
-            fvm::ddt(psi, p)
+            fvc::ddt(rho) + psi*correction(fvm::ddt(p))
           + fvc::div(phiHbyA)
           - fvm::laplacian(rhorAUf, p)
          ==
@@ -95,11 +95,6 @@ if (pressureControl.limit(p))
     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/rhoPimpleFoam.C b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C
index 0fbe7717e3..5d5ed3f5f4 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C
+++ b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C
@@ -122,6 +122,8 @@ int main(int argc, char *argv[])
             }
         }
 
+        rho = thermo.rho();
+
         runTime.write();
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H
index 15c87f9c85..c097118753 100644
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H
@@ -1,115 +1,74 @@
+rho = thermo.rho();
+
+// Thermodynamic density needs to be updated by psi*d(p) after the
+// pressure solution
+const volScalarField psip0(psi*p);
+
+volScalarField rAU(1.0/UEqn.A());
+surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
+volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh));
+
+surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
+
+surfaceScalarField phiHbyA
+(
+    "phiHbyA",
+    (
+        fvc::flux(rho*HbyA)
+      + rhorAUf*fvc::ddtCorr(rho, U, phi)
+    )
+  + phig
+);
+
+MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
+
+// Update the pressure BCs to ensure flux consistency
+constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
+
+fvScalarMatrix p_rghDDtEqn
+(
+    fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
+  + fvc::div(phiHbyA)
+  ==
+    fvOptions(psi, p_rgh, rho.name())
+);
+
+while (pimple.correctNonOrthogonal())
 {
-    bool closedVolume = p_rgh.needReference();
-
-    dimensionedScalar compressibility = fvc::domainIntegrate(psi);
-    bool compressible = (compressibility.value() > SMALL);
-
-    rho = thermo.rho();
-
-    // Thermodynamic density needs to be updated by psi*d(p) after the
-    // pressure solution - done in 2 parts. Part 1:
-    thermo.rho() -= psi*p_rgh;
-
-    volScalarField rAU(1.0/UEqn.A());
-    surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
-    volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh));
-
-    surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
-
-    surfaceScalarField phiHbyA
+    fvScalarMatrix p_rghEqn
     (
-        "phiHbyA",
-        (
-            fvc::flux(rho*HbyA)
-          + rhorAUf*fvc::ddtCorr(rho, U, phi)
-        )
-      + phig
+        p_rghDDtEqn
+      - fvm::laplacian(rhorAUf, p_rgh)
     );
 
-    MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
+    p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
 
-    // Update the pressure BCs to ensure flux consistency
-    constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
-
-    tmp<fvScalarMatrix> p_rghDDtEqn
-    (
-        new fvScalarMatrix(p_rgh, dimMass/dimTime)
-    );
-
-    if (compressible)
+    if (pimple.finalNonOrthogonalIter())
     {
-        p_rghDDtEqn =
-        (
-            fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
-        );
+        // Calculate the conservative fluxes
+        phi = phiHbyA + p_rghEqn.flux();
+
+        // Explicitly relax pressure for momentum corrector
+        p_rgh.relax();
+
+        // Correct the momentum source with the pressure gradient flux
+        // calculated from the relaxed pressure
+        U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rhorAUf);
+        U.correctBoundaryConditions();
+        fvOptions.correct(U);
+        K = 0.5*magSqr(U);
     }
-
-    while (pimple.correctNonOrthogonal())
-    {
-        fvScalarMatrix p_rghEqn
-        (
-            p_rghDDtEqn()
-          + fvc::div(phiHbyA)
-          - fvm::laplacian(rhorAUf, p_rgh)
-        ==
-            fvOptions(psi, p_rgh, rho.name())
-        );
-
-        p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
-
-        if (pimple.finalNonOrthogonalIter())
-        {
-            // Calculate the conservative fluxes
-            phi = phiHbyA + p_rghEqn.flux();
-
-            // Explicitly relax pressure for momentum corrector
-            p_rgh.relax();
-
-            // Correct the momentum source with the pressure gradient flux
-            // calculated from the relaxed pressure
-            U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rhorAUf);
-            U.correctBoundaryConditions();
-            fvOptions.correct(U);
-            K = 0.5*magSqr(U);
-        }
-    }
-
-    p = p_rgh + rho*gh;
-
-    // Second part of thermodynamic density update
-    thermo.rho() += psi*p_rgh;
-
-    if (thermo.dpdt())
-    {
-        dpdt = fvc::ddt(p);
-    }
-
-    if (compressible)
-    {
-        #include "rhoEqn.H"
-    }
-    #include "compressibleContinuityErrs.H"
-
-    if (closedVolume)
-    {
-        if (!compressible)
-        {
-            p += dimensionedScalar
-            (
-                "p",
-                p.dimensions(),
-                pRefValue - getRefCellValue(p, pRefCell)
-            );
-        }
-        else
-        {
-            p += (initialMass - fvc::domainIntegrate(thermo.rho()))
-                /compressibility;
-            rho = thermo.rho();
-        }
-        p_rgh = p - rho*gh;
-    }
-
-    Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value()
-        << endl;
 }
+
+p = p_rgh + rho*gh;
+
+// Thermodynamic density update
+thermo.correctRho(psi*p - psip0);
+
+if (thermo.dpdt())
+{
+    dpdt = fvc::ddt(p);
+}
+
+#include "rhoEqn.H"
+#include "compressibleContinuityErrs.H"
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H
index 09d90906c0..bd955882b1 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H
@@ -1,114 +1,100 @@
+bool closedVolume = p_rgh.needReference();
+dimensionedScalar compressibility = fvc::domainIntegrate(psi);
+bool compressible = (compressibility.value() > SMALL);
+
+rho = thermo.rho();
+
+volScalarField rAU("rAU", 1.0/UEqn.A());
+surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
+volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh));
+
+surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
+
+surfaceScalarField phiHbyA
+(
+    "phiHbyA",
+    (
+        fvc::flux(rho*HbyA)
+      + rhorAUf*fvc::ddtCorr(rho, U, phi)
+    )
+  + phig
+);
+
+MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
+
+// Update the pressure BCs to ensure flux consistency
+constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
+
 {
-    bool closedVolume = p_rgh.needReference();
-    dimensionedScalar compressibility = fvc::domainIntegrate(psi);
-    bool compressible = (compressibility.value() > SMALL);
-
-    rho = thermo.rho();
-
-    volScalarField rAU("rAU", 1.0/UEqn.A());
-    surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
-    volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh));
-
-    surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
-
-    surfaceScalarField phiHbyA
+    fvScalarMatrix p_rghDDtEqn
     (
-        "phiHbyA",
-        (
-            fvc::flux(rho*HbyA)
-          + rhorAUf*fvc::ddtCorr(rho, U, phi)
-        )
-      + phig
+        fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
+      + fvc::div(phiHbyA)
     );
 
-    MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
-
-    // Update the pressure BCs to ensure flux consistency
-    constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
-
-    tmp<fvScalarMatrix> p_rghDDtEqn
-    (
-        new fvScalarMatrix(p_rgh, dimMass/dimTime)
-    );
+    // Thermodynamic density needs to be updated by psi*d(p) after the
+    // pressure solution
+    const volScalarField psip0(psi*p);
 
+    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
-        if (compressible)
-        {
-            p_rghDDtEqn =
-            (
-                fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
-             ==
-                fvOptions(psi, p_rgh, rho.name())
-            );
-        }
+        fvScalarMatrix p_rghEqn
+        (
+            p_rghDDtEqn
+          - fvm::laplacian(rhorAUf, p_rgh)
+        );
 
-        // Thermodynamic density needs to be updated by psi*d(p) after the
-        // pressure solution - done in 2 parts. Part 1:
-        thermo.rho() -= psi*p_rgh;
-
-        for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
-        {
-            fvScalarMatrix p_rghEqn
+        solvPerfp_rgh = p_rghEqn.solve
+        (
+            mesh.solver
             (
-                p_rghDDtEqn()
-              + fvc::div(phiHbyA)
-              - fvm::laplacian(rhorAUf, p_rgh)
-            );
-
-            p_rghEqn.solve
-            (
-                mesh.solver
+                p_rgh.select
                 (
-                    p_rgh.select
                     (
-                        (
-                           oCorr == nOuterCorr-1
-                        && corr == nCorr-1
-                        && nonOrth == nNonOrthCorr
-                        )
+                        oCorr == nOuterCorr-1
+                     && corr == nCorr-1
+                     && nonOrth == nNonOrthCorr
                     )
                 )
-            );
+            )
+        );
 
-            if (nonOrth == nNonOrthCorr)
-            {
-                phi = phiHbyA + p_rghEqn.flux();
-                U = HbyA
-                  + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rhorAUf);
-                U.correctBoundaryConditions();
-                fvOptions.correct(U);
-                K = 0.5*magSqr(U);
-            }
+        if (nonOrth == nNonOrthCorr)
+        {
+            phi = phiHbyA + p_rghEqn.flux();
+            U = HbyA
+              + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rhorAUf);
+            U.correctBoundaryConditions();
+            fvOptions.correct(U);
+            K = 0.5*magSqr(U);
         }
-
-        // Second part of thermodynamic density update
-        thermo.rho() += psi*p_rgh;
     }
 
     p = p_rgh + rho*gh;
 
-    // Update pressure time derivative if needed
-    if (thermo.dpdt())
-    {
-        dpdt = fvc::ddt(p);
-    }
-
-    if (compressible)
-    {
-        // Solve continuity
-        #include "rhoEqn.H"
-    }
-
-    // Update continuity errors
-    #include "compressibleContinuityErrors.H"
-
-    // For closed-volume cases adjust the pressure and density levels
-    // to obey overall mass continuity
-    if (closedVolume && compressible)
-    {
-        p += (initialMass - fvc::domainIntegrate(thermo.rho()))
-            /compressibility;
-        rho = thermo.rho();
-        p_rgh = p - rho*gh;
-    }
+    // Thermodynamic density update
+    thermo.correctRho(psi*p - psip0);
+}
+
+
+// Update pressure time derivative if needed
+if (thermo.dpdt())
+{
+    dpdt = fvc::ddt(p);
+}
+
+// Solve continuity
+#include "rhoEqn.H"
+
+// Update continuity errors
+#include "compressibleContinuityErrors.H"
+
+// For closed-volume cases adjust the pressure and density levels
+// to obey overall mass continuity
+if (closedVolume && compressible)
+{
+    p += (initialMass - fvc::domainIntegrate(thermo.rho()))
+        /compressibility;
+    rho = thermo.rho();
+    p_rgh = p - rho*gh;
 }
diff --git a/applications/solvers/lagrangian/coalChemistryFoam/coalChemistryFoam.C b/applications/solvers/lagrangian/coalChemistryFoam/coalChemistryFoam.C
index 544cccdae2..c7ea95929c 100644
--- a/applications/solvers/lagrangian/coalChemistryFoam/coalChemistryFoam.C
+++ b/applications/solvers/lagrangian/coalChemistryFoam/coalChemistryFoam.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
@@ -42,6 +42,7 @@ Description
 #include "radiationModel.H"
 #include "SLGThermo.H"
 #include "pimpleControl.H"
+#include "pressureControl.H"
 #include "localEulerDdtScheme.H"
 #include "fvcSmooth.H"
 
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H b/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H
index f6c8ceec21..f2cb294718 100644
--- a/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H
+++ b/applications/solvers/lagrangian/reactingParcelFoam/pEqn.H
@@ -1,73 +1,71 @@
+rho = thermo.rho();
+
+// Thermodynamic density needs to be updated by psi*d(p) after the
+// pressure solution
+const volScalarField psip0(psi*p);
+
+volScalarField rAU(1.0/UEqn.A());
+surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
+volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
+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);
+
+fvScalarMatrix pDDtEqn
+(
+    fvc::ddt(rho) + psi*correction(fvm::ddt(p))
+  + fvc::div(phiHbyA)
+ ==
+    parcels.Srho()
+  + fvOptions(psi, p, rho.name())
+);
+
+while (pimple.correctNonOrthogonal())
 {
-    rho = thermo.rho();
-
-    // Thermodynamic density needs to be updated by psi*d(p) after the
-    // pressure solution - done in 2 parts. Part 1:
-    thermo.rho() -= psi*p;
-
-    volScalarField rAU(1.0/UEqn.A());
-    surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
-    volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
-    surfaceScalarField phiHbyA
+    fvScalarMatrix pEqn
     (
-        "phiHbyA",
-        (
-            fvc::flux(rho*HbyA)
-          + rhorAUf*fvc::ddtCorr(rho, U, phi)
-        )
+        pDDtEqn
+      - fvm::laplacian(rhorAUf, p)
     );
 
-    MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
+    pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
-    // Update the pressure BCs to ensure flux consistency
-    constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF);
-
-    fvScalarMatrix pDDtEqn
-    (
-        fvc::ddt(rho) + psi*correction(fvm::ddt(p))
-      + fvc::div(phiHbyA)
-     ==
-        parcels.Srho()
-      + fvOptions(psi, p, rho.name())
-    );
-
-    while (pimple.correctNonOrthogonal())
+    if (pimple.finalNonOrthogonalIter())
     {
-        fvScalarMatrix pEqn
-        (
-            pDDtEqn
-          - fvm::laplacian(rhorAUf, p)
-        );
-
-        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
-
-        if (pimple.finalNonOrthogonalIter())
-        {
-            phi = phiHbyA + pEqn.flux();
-        }
+        phi = phiHbyA + pEqn.flux();
     }
-
-    p.relax();
-
-    // Second part of thermodynamic density update
-    thermo.rho() += psi*p;
-
-    #include "rhoEqn.H" // NOTE: flux and time scales now inconsistent
-    #include "compressibleContinuityErrs.H"
-
-    U = HbyA - rAU*fvc::grad(p);
-    U.correctBoundaryConditions();
-    fvOptions.correct(U);
-    K = 0.5*magSqr(U);
-
-    if (thermo.dpdt())
-    {
-        dpdt = fvc::ddt(p);
-    }
-
-    rho = thermo.rho();
-    rho = max(rho, rhoMin);
-    rho = min(rho, rhoMax);
-
-    Info<< "p min/max = " << min(p).value() << ", " << max(p).value() << endl;
 }
+
+p.relax();
+
+// Thermodynamic density update
+thermo.correctRho(psi*p - psip0);
+
+#include "rhoEqn.H" // NOTE: flux and time scales now inconsistent
+#include "compressibleContinuityErrs.H"
+
+U = HbyA - rAU*fvc::grad(p);
+U.correctBoundaryConditions();
+fvOptions.correct(U);
+K = 0.5*magSqr(U);
+
+if (thermo.dpdt())
+{
+    dpdt = fvc::ddt(p);
+}
+
+rho = thermo.rho();
+rho = max(rho, rhoMin);
+rho = min(rho, rhoMax);
+
+Info<< "p min/max = " << min(p).value() << ", " << max(p).value() << endl;
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/pEqn.H b/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/pEqn.H
index 055eff6f0b..1ce5db0ec8 100644
--- a/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/pEqn.H
+++ b/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/pEqn.H
@@ -1,57 +1,55 @@
-{
-    // Thermodynamic density needs to be updated by psi*d(p) after the
-    // pressure solution - done in 2 parts. Part 1:
-    thermo.rho() -= psi*p;
+// Thermodynamic density needs to be updated by psi*d(p) after the
+// pressure solution
+const volScalarField psip0(psi*p);
 
-    volScalarField rAU(1.0/UEqn.A());
-    surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
-    volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
-    tUEqn.clear();
-    surfaceScalarField phiHbyA
+volScalarField rAU(1.0/UEqn.A());
+surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
+volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
+tUEqn.clear();
+surfaceScalarField phiHbyA
+(
+    "phiHbyA",
+    fvc::interpolate(rho)*fvc::flux(HbyA)
+);
+
+MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
+
+// Update the pressure BCs to ensure flux consistency
+constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF);
+
+while (simple.correctNonOrthogonal())
+{
+    fvScalarMatrix pEqn
     (
-        "phiHbyA",
-        fvc::interpolate(rho)*fvc::flux(HbyA)
+        fvc::div(phiHbyA)
+      - fvm::laplacian(rhorAUf, p)
+     ==
+        parcels.Srho()
+      + fvOptions(psi, p, rho.name())
     );
 
-    MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
+    pEqn.solve();
 
-    // Update the pressure BCs to ensure flux consistency
-    constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF);
-
-    while (simple.correctNonOrthogonal())
+    if (simple.finalNonOrthogonalIter())
     {
-        fvScalarMatrix pEqn
-        (
-            fvc::div(phiHbyA)
-          - fvm::laplacian(rhorAUf, p)
-         ==
-            parcels.Srho()
-          + fvOptions(psi, p, rho.name())
-        );
-
-        pEqn.solve();
-
-        if (simple.finalNonOrthogonalIter())
-        {
-            phi = phiHbyA + pEqn.flux();
-        }
+        phi = phiHbyA + pEqn.flux();
     }
-
-    p.relax();
-
-    // Second part of thermodynamic density update
-    thermo.rho() += psi*p;
-
-    #include "compressibleContinuityErrs.H"
-
-    U = HbyA - rAU*fvc::grad(p);
-    U.correctBoundaryConditions();
-    fvOptions.correct(U);
-
-    rho = thermo.rho();
-    rho = max(rho, rhoMin);
-    rho = min(rho, rhoMax);
-    rho.relax();
-
-    Info<< "p min/max = " << min(p).value() << ", " << max(p).value() << endl;
 }
+
+p.relax();
+
+// Thermodynamic density update
+thermo.correctRho(psi*p - psip0);
+
+#include "compressibleContinuityErrs.H"
+
+U = HbyA - rAU*fvc::grad(p);
+U.correctBoundaryConditions();
+fvOptions.correct(U);
+
+rho = thermo.rho();
+rho = max(rho, rhoMin);
+rho = min(rho, rhoMax);
+rho.relax();
+
+Info<< "p min/max = " << min(p).value() << ", " << max(p).value() << endl;
diff --git a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H
index 7cd3875e85..e84c19453f 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/pEqn.H
@@ -104,8 +104,8 @@
     }
 
     // Update densities from change in p_rgh
-    rho1 += psi1*(p_rgh - p_rgh_0);
-    rho2 += psi2*(p_rgh - p_rgh_0);
+    mixture.thermo1().correctRho(psi1*(p_rgh - p_rgh_0));
+    mixture.thermo2().correctRho(psi2*(p_rgh - p_rgh_0));
 
     rho = alpha1*rho1 + alpha2*rho2;
 
diff --git a/applications/solvers/multiphase/compressibleInterFoam/createFields.H b/applications/solvers/multiphase/compressibleInterFoam/createFields.H
index 46168e6b87..0616b80ea8 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/createFields.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/createFields.H
@@ -38,8 +38,8 @@ volScalarField& alpha2(mixture.alpha2());
 
 Info<< "Reading thermophysical properties\n" << endl;
 
-volScalarField& rho1 = mixture.thermo1().rho();
-volScalarField& rho2 = mixture.thermo2().rho();
+const volScalarField& rho1 = mixture.thermo1().rho();
+const volScalarField& rho2 = mixture.thermo2().rho();
 
 volScalarField rho
 (
diff --git a/applications/solvers/multiphase/compressibleInterFoam/pEqn.H b/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
index 51c58ba536..93412c9dba 100644
--- a/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
@@ -104,8 +104,8 @@
     }
 
     // Update densities from change in p_rgh
-    rho1 += psi1*(p_rgh - p_rgh_0);
-    rho2 += psi2*(p_rgh - p_rgh_0);
+    mixture.thermo1().correctRho(psi1*(p_rgh - p_rgh_0));
+    mixture.thermo2().correctRho(psi2*(p_rgh - p_rgh_0));
 
     rho = alpha1*rho1 + alpha2*rho2;
 
diff --git a/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.C b/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.C
index e500ab7587..6d031f05b6 100644
--- a/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.C
+++ b/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.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
@@ -45,6 +45,7 @@ void Foam::pimpleControl::read()
     solveFlow_ = pimpleDict.lookupOrDefault<Switch>("solveFlow", true);
     nCorrPIMPLE_ = pimpleDict.lookupOrDefault<label>("nOuterCorrectors", 1);
     nCorrPISO_ = pimpleDict.lookupOrDefault<label>("nCorrectors", 1);
+    SIMPLErho_ = pimpleDict.lookupOrDefault<Switch>("SIMPLErho", false);
     turbOnFinalIterOnly_ =
         pimpleDict.lookupOrDefault<Switch>("turbOnFinalIterOnly", true);
 }
@@ -128,6 +129,7 @@ Foam::pimpleControl::pimpleControl(fvMesh& mesh, const word& dictName)
     nCorrPIMPLE_(0),
     nCorrPISO_(0),
     corrPISO_(0),
+    SIMPLErho_(false),
     turbOnFinalIterOnly_(true),
     converged_(false)
 {
diff --git a/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.H b/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.H
index 033ce84d2f..5e87b20ea3 100644
--- a/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.H
+++ b/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.H
@@ -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
@@ -81,6 +81,10 @@ protected:
             //- Current PISO corrector
             label corrPISO_;
 
+            //- Flag to indicate whether to update density in SIMPLE
+            //  rather than PISO mode
+            bool SIMPLErho_;
+
             //- Flag to indicate whether to only solve turbulence on final iter
             bool turbOnFinalIterOnly_;
 
@@ -128,6 +132,10 @@ public:
             //- Current PISO corrector index
             inline label corrPISO() const;
 
+            //- Flag to indicate whether to update density in SIMPLE
+            //  rather than PISO mode
+            inline bool SIMPLErho() const;
+
 
         // Solution control
 
diff --git a/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControlI.H b/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControlI.H
index f6c6e0ce35..deed2958c0 100644
--- a/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControlI.H
+++ b/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControlI.H
@@ -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
@@ -43,6 +43,12 @@ inline Foam::label Foam::pimpleControl::corrPISO() const
 }
 
 
+inline bool Foam::pimpleControl::SIMPLErho() const
+{
+    return SIMPLErho_;
+}
+
+
 inline bool Foam::pimpleControl::correct()
 {
     corrPISO_++;
diff --git a/src/thermophysicalModels/basic/fluidThermo/fluidThermo.H b/src/thermophysicalModels/basic/fluidThermo/fluidThermo.H
index e18cda7343..995b1e9ad5 100644
--- a/src/thermophysicalModels/basic/fluidThermo/fluidThermo.H
+++ b/src/thermophysicalModels/basic/fluidThermo/fluidThermo.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -103,6 +103,10 @@ public:
 
         // Access to thermodynamic state variables
 
+            //- Add the given density correction to the density field.
+            //  Used to update the density field following pressure solution
+            virtual void correctRho(const volScalarField& deltaRho) = 0;
+
             //- Compressibility [s^2/m^2]
             virtual const volScalarField& psi() const = 0;
 
diff --git a/src/thermophysicalModels/basic/psiThermo/psiThermo.C b/src/thermophysicalModels/basic/psiThermo/psiThermo.C
index 630402280a..fff1f8dfd6 100644
--- a/src/thermophysicalModels/basic/psiThermo/psiThermo.C
+++ b/src/thermophysicalModels/basic/psiThermo/psiThermo.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -102,6 +102,10 @@ Foam::tmp<Foam::scalarField> Foam::psiThermo::rho(const label patchi) const
 }
 
 
+void Foam::psiThermo::correctRho(const Foam::volScalarField& deltaRho)
+{}
+
+
 const Foam::volScalarField& Foam::psiThermo::psi() const
 {
     return psi_;
diff --git a/src/thermophysicalModels/basic/psiThermo/psiThermo.H b/src/thermophysicalModels/basic/psiThermo/psiThermo.H
index 76404a026b..454eb5674f 100644
--- a/src/thermophysicalModels/basic/psiThermo/psiThermo.H
+++ b/src/thermophysicalModels/basic/psiThermo/psiThermo.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -115,6 +115,11 @@ public:
 
         // Fields derived from thermodynamic state variables
 
+            //- Add the given density correction to the density field.
+            //  Used to update the density field following pressure solution.
+            //  For psiThermo does nothing.
+            virtual void correctRho(const volScalarField& deltaRho);
+
             //- Density [kg/m^3] - uses current value of pressure
             virtual tmp<volScalarField> rho() const;
 
diff --git a/src/thermophysicalModels/basic/rhoThermo/rhoThermo.C b/src/thermophysicalModels/basic/rhoThermo/rhoThermo.C
index ada7716c4a..2f2c5d893f 100644
--- a/src/thermophysicalModels/basic/rhoThermo/rhoThermo.C
+++ b/src/thermophysicalModels/basic/rhoThermo/rhoThermo.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -173,6 +173,12 @@ Foam::volScalarField& Foam::rhoThermo::rho()
 }
 
 
+void Foam::rhoThermo::correctRho(const Foam::volScalarField& deltaRho)
+{
+    rho_ += deltaRho;
+}
+
+
 const Foam::volScalarField& Foam::rhoThermo::psi() const
 {
     return psi_;
diff --git a/src/thermophysicalModels/basic/rhoThermo/rhoThermo.H b/src/thermophysicalModels/basic/rhoThermo/rhoThermo.H
index 680094b90f..a4abe62a53 100644
--- a/src/thermophysicalModels/basic/rhoThermo/rhoThermo.H
+++ b/src/thermophysicalModels/basic/rhoThermo/rhoThermo.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -136,6 +136,10 @@ public:
             //- Return non-const access to the local density field [kg/m^3]
             virtual volScalarField& rho();
 
+            //- Add the given density correction to the density field.
+            //  Used to update the density field following pressure solution
+            virtual void correctRho(const volScalarField& deltaRho);
+
             //- Compressibility [s^2/m^2]
             virtual const volScalarField& psi() const;
 
diff --git a/tutorials/compressible/rhoPimpleFoam/RAS/angledDuct/system/fvSolution b/tutorials/compressible/rhoPimpleFoam/RAS/angledDuct/system/fvSolution
index 5dc858a4d4..828b16db4f 100644
--- a/tutorials/compressible/rhoPimpleFoam/RAS/angledDuct/system/fvSolution
+++ b/tutorials/compressible/rhoPimpleFoam/RAS/angledDuct/system/fvSolution
@@ -54,6 +54,7 @@ PIMPLE
     nCorrectors         1;
     nNonOrthogonalCorrectors 0;
     consistent          yes;
+    SIMPLErho           yes;
 
     pMaxFactor          1.5;
     pMinFactor          0.9;
diff --git a/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/0/T b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/0/T
new file mode 100644
index 0000000000..2df93ba37e
--- /dev/null
+++ b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/0/T
@@ -0,0 +1,42 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      T;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 1 0 0 0];
+
+internalField   uniform 300;
+
+boundaryField
+{
+    Default_Boundary_Region
+    {
+        type            fixedValue;
+        value           uniform 350;
+    }
+
+    inlet
+    {
+        type            fixedValue;
+        value           $internalField;
+    }
+
+    outlet
+    {
+        type            inletOutlet;
+        inletValue      $internalField;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/0/U b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/0/U
new file mode 100644
index 0000000000..03d1d6a8d8
--- /dev/null
+++ b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/0/U
@@ -0,0 +1,43 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volVectorField;
+    object      U;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    Default_Boundary_Region
+    {
+        type            noSlip;
+    }
+
+    inlet
+    {
+        type            flowRateInletVelocity;
+        massFlowRate    constant 5;
+        rhoInlet        1000;    // Guess for rho
+    }
+
+    outlet
+    {
+        type            pressureInletOutletVelocity;
+        value           $internalField;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/0/alphat b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/0/alphat
new file mode 100644
index 0000000000..61043c15b1
--- /dev/null
+++ b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/0/alphat
@@ -0,0 +1,45 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      alphat;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -1 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    Default_Boundary_Region
+    {
+        type            compressible::alphatWallFunction;
+        Prt             0.85;
+        value           uniform 0;
+    }
+
+    inlet
+    {
+        type            calculated;
+        value           uniform 0;
+    }
+
+    outlet
+    {
+        type            calculated;
+        value           uniform 0;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/0/epsilon b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/0/epsilon
new file mode 100644
index 0000000000..20d9e01d3f
--- /dev/null
+++ b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/0/epsilon
@@ -0,0 +1,49 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      epsilon;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -3 0 0 0 0];
+
+internalField   uniform 200;
+
+boundaryField
+{
+    Default_Boundary_Region
+    {
+        type            epsilonWallFunction;
+        Cmu             0.09;
+        kappa           0.41;
+        E               9.8;
+        value           uniform 200;
+    }
+
+    inlet
+    {
+        type            turbulentMixingLengthDissipationRateInlet;
+        mixingLength    0.005;
+        value           uniform 200;
+    }
+
+    outlet
+    {
+        type            inletOutlet;
+        inletValue      uniform 200;
+        value           uniform 200;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/0/k b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/0/k
new file mode 100644
index 0000000000..4bddfd6596
--- /dev/null
+++ b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/0/k
@@ -0,0 +1,46 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      k;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -2 0 0 0 0];
+
+internalField   uniform 1;
+
+boundaryField
+{
+    Default_Boundary_Region
+    {
+        type            kqRWallFunction;
+        value           uniform 1;
+    }
+
+    inlet
+    {
+        type            turbulentIntensityKineticEnergyInlet;
+        intensity       0.05;
+        value           uniform 1;
+    }
+
+    outlet
+    {
+        type            inletOutlet;
+        inletValue      uniform 1;
+        value           uniform 1;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/0/nut b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/0/nut
new file mode 100644
index 0000000000..a7ffddb820
--- /dev/null
+++ b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/0/nut
@@ -0,0 +1,47 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      nut;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -1 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    Default_Boundary_Region
+    {
+        type            nutkWallFunction;
+        Cmu             0.09;
+        kappa           0.41;
+        E               9.8;
+        value           uniform 0;
+    }
+
+    inlet
+    {
+        type            calculated;
+        value           uniform 0;
+    }
+
+    outlet
+    {
+        type            calculated;
+        value           uniform 0;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/0/p b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/0/p
new file mode 100644
index 0000000000..b4c015f222
--- /dev/null
+++ b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/0/p
@@ -0,0 +1,41 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      p;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -2 0 0 0 0];
+
+internalField   uniform 1e5;
+
+boundaryField
+{
+    Default_Boundary_Region
+    {
+        type            zeroGradient;
+    }
+
+    inlet
+    {
+        type            zeroGradient;
+    }
+
+    outlet
+    {
+        type            fixedValue;
+        value           $internalField;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/constant/thermophysicalProperties b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/constant/thermophysicalProperties
new file mode 100644
index 0000000000..426ab72116
--- /dev/null
+++ b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/constant/thermophysicalProperties
@@ -0,0 +1,31 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      thermophysicalProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+thermoType
+{
+    type            heRhoThermo;
+    mixture         pureMixture;
+    properties      liquid;
+    energy          sensibleInternalEnergy;
+}
+
+mixture
+{
+    H2O;
+}
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/constant/turbulenceProperties b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/constant/turbulenceProperties
new file mode 100644
index 0000000000..cd2daf8229
--- /dev/null
+++ b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/constant/turbulenceProperties
@@ -0,0 +1,30 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      turbulenceProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+simulationType  RAS;
+
+RAS
+{
+    RASModel        kEpsilon;
+
+    turbulence      on;
+
+    printCoeffs     on;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/system/blockMeshDict b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/system/blockMeshDict
new file mode 100644
index 0000000000..9bec6a9547
--- /dev/null
+++ b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/system/blockMeshDict
@@ -0,0 +1,127 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+convertToMeters 0.001;
+
+vertices
+(
+    // front-plane: z = +25mm
+    // inlet region
+    (   -50  25   25)           // pt 0
+    (     0  25   25)           // pt 1
+    (   -50  75   25)           // pt 2
+    (     0  75   25)           // pt 3
+    // outlet region
+    (  -500 -75   25)           // pt 4
+    (     0 -75   25)           // pt 5
+    (  -500 -25   25)           // pt 6
+    (     0 -25   25)           // pt 7
+    // bend mid-points
+    (    25   0   25)           // pt 8
+    (    75   0   25)           // pt 9
+    // back-plane: z = -25mm
+    // inlet region
+    (   -50  25   -25)          // pt 0 + 10
+    (     0  25   -25)          // pt 1 + 10
+    (   -50  75   -25)          // pt 2 + 10
+    (     0  75   -25)          // pt 3 + 10
+    // outlet region
+    (  -500 -75   -25)          // pt 4 + 10
+    (     0 -75   -25)          // pt 5 + 10
+    (  -500 -25   -25)          // pt 7 + 10
+    (     0 -25   -25)          // pt 8 + 10
+    // bend mid-points
+    (    25   0   -25)          // pt 8 + 10
+    (    75   0   -25)          // pt 9 + 10
+);
+
+blocks
+(
+    hex (0 1 11 10  2 3 13 12) inlet  ( 20 20 20)  simpleGrading (1 1 1)
+    hex (4 5 15 14  6 7 17 16) outlet (200 20 20)  simpleGrading (1 1 1)
+
+    hex (1 8 18 11  3 9 19 13) bend1  ( 30 20 20)  simpleGrading (1 1 1)
+    hex (5 9 19 15  7 8 18 17) bend2  ( 30 20 20)  simpleGrading (1 1 1)
+);
+
+edges
+(
+   // block 2
+   arc  1  8  ( 17.678  17.678  25)
+   arc 11 18  ( 17.678  17.678 -25)
+   arc  3  9  ( 53.033  53.033  25)
+   arc 13 19  ( 53.033  53.033 -25)
+   // block 3
+   arc  7  8  ( 17.678  -17.678  25)
+   arc 17 18  ( 17.678  -17.678 -25)
+   arc  5  9  ( 53.033  -53.033  25)
+   arc 15 19  ( 53.033  -53.033 -25)
+);
+
+boundary
+(
+        // is there no way of defining all my 'defaultFaces' to be 'wall'?
+    Default_Boundary_Region
+    {
+        type wall;
+        faces
+        (
+            // block0
+            ( 0 1 3 2 )
+            ( 11 10 12 13 )
+            ( 0 10 11 1 )
+            ( 2 3 13 12 )
+            // block1
+            ( 4 5 7 6 )
+            ( 15 14 16 17 )
+            ( 4 14 15 5 )
+            ( 6 7 17 16 )
+            // block2
+            ( 1 8 9 3 )
+            ( 18 11 13 19 )
+            ( 3 9 19 13 )
+            ( 1 11 18 8 )
+            // block3
+            ( 5 9 8 7 )
+            ( 19 15 17 18 )
+            ( 5 15 19 9 )
+            ( 7 8 18 17 )
+        );
+    }
+    inlet
+    {
+        type patch;
+        faces
+        (
+            (0 2 12 10)
+        );
+    }
+
+    outlet
+    {
+        type patch;
+        faces
+        (
+            (4 6 16 14)
+        );
+    }
+);
+
+mergePatchPairs
+(
+);
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/system/controlDict b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/system/controlDict
new file mode 100644
index 0000000000..07384c85ac
--- /dev/null
+++ b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/system/controlDict
@@ -0,0 +1,51 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+application     rhoPimpleFoam;
+
+startFrom       startTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         0.5;
+
+deltaT          2e-3;
+
+writeControl    timeStep;
+
+writeInterval   10;
+
+purgeWrite      0;
+
+writeFormat     ascii;
+
+writePrecision  6;
+
+writeCompression off;
+
+timeFormat      general;
+
+timePrecision   6;
+
+graphFormat     raw;
+
+runTimeModifiable true;
+
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/system/decomposeParDict b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/system/decomposeParDict
new file mode 100644
index 0000000000..9e844a62ba
--- /dev/null
+++ b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/system/decomposeParDict
@@ -0,0 +1,45 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      decomposeParDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+numberOfSubdomains 8;
+
+method          hierarchical;
+
+simpleCoeffs
+{
+    n               (8 1 1);
+    delta           0.001;
+}
+
+hierarchicalCoeffs
+{
+    n               (4 2 1);
+    delta           0.001;
+    order           xyz;
+}
+
+manualCoeffs
+{
+    dataFile        "";
+}
+
+distributed     no;
+
+roots           ( );
+
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/system/fvSchemes b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/system/fvSchemes
new file mode 100644
index 0000000000..8c9c2b631b
--- /dev/null
+++ b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/system/fvSchemes
@@ -0,0 +1,60 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+    default             Euler;
+}
+
+gradSchemes
+{
+    default             Gauss linear;
+
+    limited             cellLimited Gauss linear 1;
+}
+
+divSchemes
+{
+    default             none;
+
+    div(phi,U)          Gauss linearUpwind limited;
+    div(phi,e)          Gauss linearUpwind limited;
+    div(phi,epsilon)    Gauss linearUpwind limited;
+    div(phi,k)          Gauss linearUpwind limited;
+    div(phi,K)          Gauss linearUpwind limited;
+    div(phiv,p)         Gauss linearUpwind limited;
+
+    div(((rho*nuEff)*dev2(T(grad(U)))))      Gauss linear;
+}
+
+laplacianSchemes
+{
+    default         Gauss linear corrected;
+}
+
+interpolationSchemes
+{
+    default         linear;
+}
+
+snGradSchemes
+{
+    default         corrected;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/system/fvSolution b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/system/fvSolution
new file mode 100644
index 0000000000..173adf6c59
--- /dev/null
+++ b/tutorials/compressible/rhoPimpleFoam/RAS/squareBendLiq/system/fvSolution
@@ -0,0 +1,63 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+    "rho.*"
+    {
+        solver nthn;
+    }
+
+    "p.*"
+    {
+        solver          GAMG;
+        smoother        GaussSeidel;
+
+        tolerance       1e-7;
+        relTol          0;
+    }
+
+    "(U|e|k|epsilon).*"
+    {
+        solver          PBiCGStab;
+        preconditioner  DILU;
+
+        tolerance       1e-7;
+        relTol          0;
+    }
+}
+
+PIMPLE
+{
+    nCorrectors     3;
+    nNonOrthogonalCorrectors 0;
+    pMinFactor      0.1;
+    pMaxFactor      1.5;
+
+    transonic       no;
+    consistent      no;
+}
+
+relaxationFactors
+{
+    equations
+    {
+        ".*"            1;
+    }
+}
+
+// ************************************************************************* //

From 111ec17cb5c0d44af76c3f840e2f44c9f36fce7c Mon Sep 17 00:00:00 2001
From: Andrew Heather <a.heather@opencfd.co.uk>
Date: Tue, 30 May 2017 12:53:59 +0100
Subject: [PATCH 4/6] GIT: Updated after cherry-picks from third party code

---
 .../solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H        | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H
index bd955882b1..d515c6a183 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H
@@ -44,7 +44,7 @@ constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
           - fvm::laplacian(rhorAUf, p_rgh)
         );
 
-        solvPerfp_rgh = p_rghEqn.solve
+        p_rghEqn.solve
         (
             mesh.solver
             (

From d1b651533f9295aa4ac1656d87a5a5e45b79f37c Mon Sep 17 00:00:00 2001
From: sergio <sergio>
Date: Thu, 1 Jun 2017 12:39:28 -0700
Subject: [PATCH 5/6] Modification on rhoPimpleFoam pEq's for handling rho
 thermo and incompressible EoS. Adding rho limiters if p is limited. This is
 important when LTS stepping or large Co number is used.

Updating rhoBuoyantPimpleFoam to handle closed domain for rho thermo and incompressible Eos.
Consolidating chtMultiRegionSimpleFoam and chtMultiRegionFoam pEqs to use the same formulation as rhoBuoyantPimpleFoam and
rhoBuoyantSimpleFoam
---
 .../compressible/rhoPimpleFoam/createFields.H | 23 ++++++++
 .../solvers/compressible/rhoPimpleFoam/pEqn.H | 11 +++-
 .../compressible/rhoPimpleFoam/pcEqn.H        | 12 +++-
 .../heatTransfer/buoyantPimpleFoam/pEqn.H     | 46 +++++++++++++--
 .../buoyantSimpleFoam/createFields.H          | 23 ++++++++
 .../heatTransfer/buoyantSimpleFoam/pEqn.H     | 17 +++---
 .../chtMultiRegionSimpleFoam/fluid/pEqn.H     | 20 ++++++-
 .../fluid/createFluidFields.H                 | 19 +++++++
 .../chtMultiRegionFoam/fluid/pEqn.H           | 57 ++++++++++++++-----
 .../fluid/setRegionFluidFields.H              |  3 +
 .../basic/fluidThermo/fluidThermo.H           |  7 ++-
 .../basic/psiThermo/psiThermo.C               |  8 ++-
 .../basic/psiThermo/psiThermo.H               |  7 ++-
 .../basic/rhoThermo/rhoThermo.C               | 10 +++-
 .../basic/rhoThermo/rhoThermo.H               |  8 ++-
 15 files changed, 229 insertions(+), 42 deletions(-)

diff --git a/applications/solvers/compressible/rhoPimpleFoam/createFields.H b/applications/solvers/compressible/rhoPimpleFoam/createFields.H
index 3ec620dc22..8079428225 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/createFields.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/createFields.H
@@ -73,3 +73,26 @@ Info<< "Creating field kinetic energy K\n" << endl;
 volScalarField K("K", 0.5*magSqr(U));
 
 #include "createMRF.H"
+
+
+dimensionedScalar rhoMax
+(
+    dimensionedScalar::lookupOrDefault
+    (
+        "rhoMax",
+        pimple.dict(),
+        dimDensity,
+        GREAT
+    )
+);
+
+dimensionedScalar rhoMin
+(
+    dimensionedScalar::lookupOrDefault
+    (
+        "rhoMin",
+        pimple.dict(),
+        dimDensity,
+        0
+    )
+);
diff --git a/applications/solvers/compressible/rhoPimpleFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
index 83eeb8fc68..3a0ca9d042 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
@@ -84,9 +84,6 @@ else
     }
 }
 
-// Thermodynamic density update
-thermo.correctRho(psi*p - psip0);
-
 #include "rhoEqn.H"
 #include "compressibleContinuityErrs.H"
 
@@ -101,12 +98,20 @@ K = 0.5*magSqr(U);
 if (pressureControl.limit(p))
 {
     p.correctBoundaryConditions();
+    thermo.correctRho(psi*p - psip0, rhoMin, rhoMax);
     rho = thermo.rho();
 }
 else if (pimple.SIMPLErho())
 {
+    thermo.correctRho(psi*p - psip0, rhoMin, rhoMax);
     rho = thermo.rho();
 }
+else
+{
+    thermo.correctRho(psi*p - psip0, rhoMin, rhoMax) ;
+}
+
+Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl;
 
 if (thermo.dpdt())
 {
diff --git a/applications/solvers/compressible/rhoPimpleFoam/pcEqn.H b/applications/solvers/compressible/rhoPimpleFoam/pcEqn.H
index e9b849bc05..546dfaf251 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/pcEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/pcEqn.H
@@ -95,9 +95,6 @@ else
     }
 }
 
-// Thermodynamic density update
-thermo.correctRho(psi*p - psip0);
-
 #include "rhoEqn.H"
 #include "compressibleContinuityErrs.H"
 
@@ -112,12 +109,21 @@ K = 0.5*magSqr(U);
 if (pressureControl.limit(p))
 {
     p.correctBoundaryConditions();
+    thermo.correctRho(psi*p - psip0, rhoMin, rhoMax);
     rho = thermo.rho();
 }
+
 else if (pimple.SIMPLErho())
 {
+    thermo.correctRho(psi*p - psip0, rhoMin, rhoMax);
     rho = thermo.rho();
 }
+else
+{
+    thermo.correctRho(psi*p - psip0, rhoMin, rhoMax);
+}
+
+Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl;
 
 if (thermo.dpdt())
 {
diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H
index c097118753..0f53379eef 100644
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H
@@ -1,3 +1,7 @@
+
+dimensionedScalar compressibility = fvc::domainIntegrate(psi);
+bool compressible = (compressibility.value() > SMALL);
+
 rho = thermo.rho();
 
 // Thermodynamic density needs to be updated by psi*d(p) after the
@@ -41,6 +45,12 @@ while (pimple.correctNonOrthogonal())
       - fvm::laplacian(rhorAUf, p_rgh)
     );
 
+    p_rghEqn.setReference
+    (
+        pRefCell,
+        compressible ? getRefCellValue(p_rgh, pRefCell) : pRefValue
+    );
+
     p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
 
     if (pimple.finalNonOrthogonalIter())
@@ -63,12 +73,40 @@ while (pimple.correctNonOrthogonal())
 p = p_rgh + rho*gh;
 
 // Thermodynamic density update
-thermo.correctRho(psi*p - psip0);
+//thermo.correctRho(psi*p - psip0);
+
+#include "rhoEqn.H"
+#include "compressibleContinuityErrs.H"
+
+if (p_rgh.needReference())
+{
+    if (!compressible)
+    {
+        p += dimensionedScalar
+        (
+            "p",
+            p.dimensions(),
+            pRefValue - getRefCellValue(p, pRefCell)
+        );
+    }
+    else
+    {
+        p += (initialMass - fvc::domainIntegrate(psi*p))
+            /compressibility;
+        thermo.correctRho(psi*p - psip0);
+        rho = thermo.rho();
+
+    }
+    p_rgh = p - rho*gh;
+}
+else
+{
+    thermo.correctRho(psi*p - psip0);
+}
+
+rho = thermo.rho();
 
 if (thermo.dpdt())
 {
     dpdt = fvc::ddt(p);
 }
-
-#include "rhoEqn.H"
-#include "compressibleContinuityErrs.H"
diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H b/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H
index b9f275407f..64bc8558c8 100644
--- a/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H
+++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H
@@ -88,3 +88,26 @@ dimensionedScalar totalVolume = sum(mesh.V());
 
 #include "createMRF.H"
 #include "createRadiationModel.H"
+
+dimensionedScalar rhoMax
+(
+    dimensionedScalar::lookupOrDefault
+    (
+        "rhoMax",
+        simple.dict(),
+        dimDensity,
+        GREAT
+    )
+);
+
+dimensionedScalar rhoMin
+(
+    dimensionedScalar::lookupOrDefault
+    (
+        "rhoMin",
+        simple.dict(),
+        dimDensity,
+        0
+    )
+);
+
diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H
index 647cf954bb..aa3ad0f163 100644
--- a/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H
+++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H
@@ -1,7 +1,4 @@
 {
-    rho = thermo.rho();
-    rho.relax();
-
     volScalarField rAU("rAU", 1.0/UEqn.A());
     surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
     volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh));
@@ -24,6 +21,8 @@
     // Update the pressure BCs to ensure flux consistency
     constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
 
+    dimensionedScalar compressibility = fvc::domainIntegrate(psi);
+    bool compressible = (compressibility.value() > SMALL);
 
     while (simple.correctNonOrthogonal())
     {
@@ -32,7 +31,12 @@
             fvm::laplacian(rhorAUf, p_rgh) == fvc::div(phiHbyA)
         );
 
-        p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
+        p_rghEqn.setReference
+        (
+            pRefCell,
+            compressible ? getRefCellValue(p_rgh, pRefCell) : pRefValue
+        );
+
         p_rghEqn.solve();
 
         if (simple.finalNonOrthogonalIter())
@@ -51,12 +55,9 @@
         }
     }
 
-    #include "continuityErrs.H"
-
     p = p_rgh + rho*gh;
 
-    dimensionedScalar compressibility = fvc::domainIntegrate(psi);
-    bool compressible = (compressibility.value() > SMALL);
+    #include "continuityErrs.H"
 
     // For closed-volume cases adjust the pressure level
     // to obey overall mass continuity
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H
index 7fe544d6a5..5b8286ca4b 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H
@@ -1,8 +1,10 @@
 {
+    /*
     rho = thermo.rho();
     rho = max(rho, rhoMin[i]);
     rho = min(rho, rhoMax[i]);
     rho.relax();
+    */
 
     volScalarField rAU("rAU", 1.0/UEqn.A());
     surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
@@ -67,10 +69,22 @@
 
     // For closed-volume cases adjust the pressure level
     // to obey overall mass continuity
-    if (closedVolume && compressible)
+    if (closedVolume)
     {
-        p += (initialMass - fvc::domainIntegrate(thermo.rho()))
-            /compressibility;
+        if (!compressible)
+        {
+            p += dimensionedScalar
+            (
+                "p",
+                p.dimensions(),
+                pRefValue - getRefCellValue(p, pRefCell)
+            );
+        }
+        else
+        {
+            p += (initialMass - fvc::domainIntegrate(psi*p))
+                /compressibility;
+        }
         p_rgh = p - rho*gh;
     }
 
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H
index be3ec79a2c..87e372fcae 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H
@@ -19,6 +19,9 @@ List<bool> frozenFlowFluid(fluidRegions.size(), false);
 PtrList<IOMRFZoneList> MRFfluid(fluidRegions.size());
 PtrList<fv::options> fluidFvOptions(fluidRegions.size());
 
+List<label> refCellFluid(fluidRegions.size());
+List<scalar> refValueFluid(fluidRegions.size());
+
 // Populate fluid field pointer lists
 forAll(fluidRegions, i)
 {
@@ -248,4 +251,20 @@ forAll(fluidRegions, i)
     );
 
     turbulence[i].validate();
+
+    refCellFluid[i] = 0;
+    refValueFluid[i] = 0.0;
+
+    if (p_rghFluid[i].needReference())
+    {
+        setRefCell
+        (
+            thermoFluid[i].p(),
+            p_rghFluid[i],
+            pimpleDict,
+            refCellFluid[i],
+            refValueFluid[i]
+        );
+    }
+
 }
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H
index d515c6a183..0c9195ee9d 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H
@@ -4,6 +4,10 @@ bool compressible = (compressibility.value() > SMALL);
 
 rho = thermo.rho();
 
+// Thermodynamic density needs to be updated by psi*d(p) after the
+// pressure solution
+const volScalarField psip0(psi*p);
+
 volScalarField rAU("rAU", 1.0/UEqn.A());
 surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
 volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh));
@@ -32,10 +36,6 @@ constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
       + fvc::div(phiHbyA)
     );
 
-    // Thermodynamic density needs to be updated by psi*d(p) after the
-    // pressure solution
-    const volScalarField psip0(psi*p);
-
     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
         fvScalarMatrix p_rghEqn
@@ -44,6 +44,12 @@ constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
           - fvm::laplacian(rhorAUf, p_rgh)
         );
 
+        p_rghEqn.setReference
+        (
+            pRefCell,
+            compressible ? getRefCellValue(p_rgh, pRefCell) : pRefValue
+        );
+
         p_rghEqn.solve
         (
             mesh.solver
@@ -62,6 +68,9 @@ constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
         if (nonOrth == nNonOrthCorr)
         {
             phi = phiHbyA + p_rghEqn.flux();
+
+            p_rgh.relax();
+
             U = HbyA
               + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rhorAUf);
             U.correctBoundaryConditions();
@@ -73,15 +82,10 @@ constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
     p = p_rgh + rho*gh;
 
     // Thermodynamic density update
-    thermo.correctRho(psi*p - psip0);
+    //thermo.correctRho(psi*p - psip0);
 }
 
 
-// Update pressure time derivative if needed
-if (thermo.dpdt())
-{
-    dpdt = fvc::ddt(p);
-}
 
 // Solve continuity
 #include "rhoEqn.H"
@@ -91,10 +95,35 @@ if (thermo.dpdt())
 
 // For closed-volume cases adjust the pressure and density levels
 // to obey overall mass continuity
-if (closedVolume && compressible)
+if (closedVolume)
 {
-    p += (initialMass - fvc::domainIntegrate(thermo.rho()))
-        /compressibility;
-    rho = thermo.rho();
+    if (!compressible)
+    {
+        p += dimensionedScalar
+        (
+            "p",
+            p.dimensions(),
+            pRefValue - getRefCellValue(p, pRefCell)
+        );
+    }
+    else
+    {
+        p += (initialMass - fvc::domainIntegrate(psi*p))
+            /compressibility;
+        thermo.correctRho(psi*p - psip0);
+        rho = thermo.rho();
+    }
     p_rgh = p - rho*gh;
 }
+else
+{
+    thermo.correctRho(psi*p - psip0);
+}
+
+rho = thermo.rho();
+
+// Update pressure time derivative if needed
+if (thermo.dpdt())
+{
+    dpdt = fvc::ddt(p);
+}
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H
index 8d3a291bef..5edd27ac5e 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H
@@ -33,3 +33,6 @@
     );
 
     const bool frozenFlow = frozenFlowFluid[i];
+
+    const label pRefCell = refCellFluid[i];
+    const scalar pRefValue = refValueFluid[i];
diff --git a/src/thermophysicalModels/basic/fluidThermo/fluidThermo.H b/src/thermophysicalModels/basic/fluidThermo/fluidThermo.H
index 995b1e9ad5..bb0364dc79 100644
--- a/src/thermophysicalModels/basic/fluidThermo/fluidThermo.H
+++ b/src/thermophysicalModels/basic/fluidThermo/fluidThermo.H
@@ -105,7 +105,12 @@ public:
 
             //- Add the given density correction to the density field.
             //  Used to update the density field following pressure solution
-            virtual void correctRho(const volScalarField& deltaRho) = 0;
+            virtual void correctRho
+            (
+                const volScalarField& deltaRho,
+                const dimensionedScalar& rhoMin,
+                const dimensionedScalar& rhoMax
+            ) = 0;
 
             //- Compressibility [s^2/m^2]
             virtual const volScalarField& psi() const = 0;
diff --git a/src/thermophysicalModels/basic/psiThermo/psiThermo.C b/src/thermophysicalModels/basic/psiThermo/psiThermo.C
index fff1f8dfd6..b3e65677d9 100644
--- a/src/thermophysicalModels/basic/psiThermo/psiThermo.C
+++ b/src/thermophysicalModels/basic/psiThermo/psiThermo.C
@@ -102,10 +102,14 @@ Foam::tmp<Foam::scalarField> Foam::psiThermo::rho(const label patchi) const
 }
 
 
-void Foam::psiThermo::correctRho(const Foam::volScalarField& deltaRho)
+void Foam::psiThermo::correctRho
+(
+    const Foam::volScalarField& deltaRho,
+    const dimensionedScalar& rhoMin,
+    const dimensionedScalar& rhoMax
+)
 {}
 
-
 const Foam::volScalarField& Foam::psiThermo::psi() const
 {
     return psi_;
diff --git a/src/thermophysicalModels/basic/psiThermo/psiThermo.H b/src/thermophysicalModels/basic/psiThermo/psiThermo.H
index 454eb5674f..b2e05a9e9e 100644
--- a/src/thermophysicalModels/basic/psiThermo/psiThermo.H
+++ b/src/thermophysicalModels/basic/psiThermo/psiThermo.H
@@ -118,7 +118,12 @@ public:
             //- Add the given density correction to the density field.
             //  Used to update the density field following pressure solution.
             //  For psiThermo does nothing.
-            virtual void correctRho(const volScalarField& deltaRho);
+            virtual void correctRho
+            (
+                const volScalarField& deltaRho,
+                const dimensionedScalar& rhoMin,
+                const dimensionedScalar& rhoMax
+            );
 
             //- Density [kg/m^3] - uses current value of pressure
             virtual tmp<volScalarField> rho() const;
diff --git a/src/thermophysicalModels/basic/rhoThermo/rhoThermo.C b/src/thermophysicalModels/basic/rhoThermo/rhoThermo.C
index 2f2c5d893f..bbc1a64ba5 100644
--- a/src/thermophysicalModels/basic/rhoThermo/rhoThermo.C
+++ b/src/thermophysicalModels/basic/rhoThermo/rhoThermo.C
@@ -173,12 +173,18 @@ Foam::volScalarField& Foam::rhoThermo::rho()
 }
 
 
-void Foam::rhoThermo::correctRho(const Foam::volScalarField& deltaRho)
+void Foam::rhoThermo::correctRho
+(
+    const Foam::volScalarField& deltaRho,
+    const dimensionedScalar& rhoMin,
+    const dimensionedScalar& rhoMax
+)
 {
     rho_ += deltaRho;
+    rho_ = max(rho_, rhoMin);
+    rho_ = min(rho_, rhoMax);
 }
 
-
 const Foam::volScalarField& Foam::rhoThermo::psi() const
 {
     return psi_;
diff --git a/src/thermophysicalModels/basic/rhoThermo/rhoThermo.H b/src/thermophysicalModels/basic/rhoThermo/rhoThermo.H
index a4abe62a53..fd8bb6b81c 100644
--- a/src/thermophysicalModels/basic/rhoThermo/rhoThermo.H
+++ b/src/thermophysicalModels/basic/rhoThermo/rhoThermo.H
@@ -138,7 +138,13 @@ public:
 
             //- Add the given density correction to the density field.
             //  Used to update the density field following pressure solution
-            virtual void correctRho(const volScalarField& deltaRho);
+            //  Limit thermo rho between rhoMin and rhoMax
+            virtual void correctRho
+            (
+                const volScalarField& deltaRho,
+                const dimensionedScalar& rhoMin,
+                const dimensionedScalar& rhoMax
+            );
 
             //- Compressibility [s^2/m^2]
             virtual const volScalarField& psi() const;

From c44322ebb8f6c1119bc64dbeac3c2b7fb88b52cf Mon Sep 17 00:00:00 2001
From: sergio <sergio>
Date: Wed, 7 Jun 2017 12:11:01 +0100
Subject: [PATCH 6/6] Adding correctRho plus rhoMin and rhoMax limiters to
 thermos. Update pEq for close domains.

---
 applications/solvers/compressible/rhoPimpleFoam/pEqn.H   | 2 --
 applications/solvers/compressible/rhoPimpleFoam/pcEqn.H  | 3 ---
 .../solvers/heatTransfer/buoyantPimpleFoam/pEqn.H        | 6 +-----
 .../solvers/heatTransfer/buoyantSimpleFoam/pEqn.H        | 2 +-
 .../chtMultiRegionSimpleFoam/fluid/pEqn.H                | 2 +-
 .../solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H | 2 +-
 src/thermophysicalModels/basic/fluidThermo/fluidThermo.H | 9 +++++++++
 src/thermophysicalModels/basic/psiThermo/psiThermo.C     | 6 ++++++
 src/thermophysicalModels/basic/psiThermo/psiThermo.H     | 7 +++++++
 src/thermophysicalModels/basic/rhoThermo/rhoThermo.C     | 6 ++++++
 src/thermophysicalModels/basic/rhoThermo/rhoThermo.H     | 4 ++++
 11 files changed, 36 insertions(+), 13 deletions(-)

diff --git a/applications/solvers/compressible/rhoPimpleFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
index 3a0ca9d042..ab89fcd811 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
@@ -111,8 +111,6 @@ else
     thermo.correctRho(psi*p - psip0, rhoMin, rhoMax) ;
 }
 
-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/pcEqn.H b/applications/solvers/compressible/rhoPimpleFoam/pcEqn.H
index 546dfaf251..55be311749 100644
--- a/applications/solvers/compressible/rhoPimpleFoam/pcEqn.H
+++ b/applications/solvers/compressible/rhoPimpleFoam/pcEqn.H
@@ -112,7 +112,6 @@ if (pressureControl.limit(p))
     thermo.correctRho(psi*p - psip0, rhoMin, rhoMax);
     rho = thermo.rho();
 }
-
 else if (pimple.SIMPLErho())
 {
     thermo.correctRho(psi*p - psip0, rhoMin, rhoMax);
@@ -123,8 +122,6 @@ else
     thermo.correctRho(psi*p - psip0, rhoMin, rhoMax);
 }
 
-Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl;
-
 if (thermo.dpdt())
 {
     dpdt = fvc::ddt(p);
diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H
index 0f53379eef..677fabffd5 100644
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H
@@ -72,9 +72,6 @@ while (pimple.correctNonOrthogonal())
 
 p = p_rgh + rho*gh;
 
-// Thermodynamic density update
-//thermo.correctRho(psi*p - psip0);
-
 #include "rhoEqn.H"
 #include "compressibleContinuityErrs.H"
 
@@ -95,9 +92,8 @@ if (p_rgh.needReference())
             /compressibility;
         thermo.correctRho(psi*p - psip0);
         rho = thermo.rho();
-
+        p_rgh = p - rho*gh;
     }
-    p_rgh = p - rho*gh;
 }
 else
 {
diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H
index aa3ad0f163..13ffff11b6 100644
--- a/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H
+++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H
@@ -76,8 +76,8 @@
         {
             p += (initialMass - fvc::domainIntegrate(psi*p))
                 /fvc::domainIntegrate(psi);
+            p_rgh = p - rho*gh;
         }
-        p_rgh = p - rho*gh;
     }
 
     rho = thermo.rho();
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H
index 5b8286ca4b..410fc7166d 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H
@@ -84,8 +84,8 @@
         {
             p += (initialMass - fvc::domainIntegrate(psi*p))
                 /compressibility;
+            p_rgh = p - rho*gh;
         }
-        p_rgh = p - rho*gh;
     }
 
     rho = thermo.rho();
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H
index 0c9195ee9d..dc67ab09e2 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H
@@ -112,8 +112,8 @@ if (closedVolume)
             /compressibility;
         thermo.correctRho(psi*p - psip0);
         rho = thermo.rho();
+        p_rgh = p - rho*gh;
     }
-    p_rgh = p - rho*gh;
 }
 else
 {
diff --git a/src/thermophysicalModels/basic/fluidThermo/fluidThermo.H b/src/thermophysicalModels/basic/fluidThermo/fluidThermo.H
index bb0364dc79..79a313b208 100644
--- a/src/thermophysicalModels/basic/fluidThermo/fluidThermo.H
+++ b/src/thermophysicalModels/basic/fluidThermo/fluidThermo.H
@@ -105,6 +105,7 @@ public:
 
             //- Add the given density correction to the density field.
             //  Used to update the density field following pressure solution
+            //  Limit thermo rho between rhoMin and rhoMax
             virtual void correctRho
             (
                 const volScalarField& deltaRho,
@@ -112,6 +113,14 @@ public:
                 const dimensionedScalar& rhoMax
             ) = 0;
 
+
+            //- Add the given density correction to the density field.
+            //  Used to update the density field following pressure solution
+            virtual void correctRho
+            (
+                const volScalarField& deltaRho
+            ) = 0;
+
             //- Compressibility [s^2/m^2]
             virtual const volScalarField& psi() const = 0;
 
diff --git a/src/thermophysicalModels/basic/psiThermo/psiThermo.C b/src/thermophysicalModels/basic/psiThermo/psiThermo.C
index b3e65677d9..3dafd068a1 100644
--- a/src/thermophysicalModels/basic/psiThermo/psiThermo.C
+++ b/src/thermophysicalModels/basic/psiThermo/psiThermo.C
@@ -110,6 +110,12 @@ void Foam::psiThermo::correctRho
 )
 {}
 
+void Foam::psiThermo::correctRho
+(
+    const Foam::volScalarField& deltaRho
+)
+{}
+
 const Foam::volScalarField& Foam::psiThermo::psi() const
 {
     return psi_;
diff --git a/src/thermophysicalModels/basic/psiThermo/psiThermo.H b/src/thermophysicalModels/basic/psiThermo/psiThermo.H
index b2e05a9e9e..5b24276e0e 100644
--- a/src/thermophysicalModels/basic/psiThermo/psiThermo.H
+++ b/src/thermophysicalModels/basic/psiThermo/psiThermo.H
@@ -125,6 +125,13 @@ public:
                 const dimensionedScalar& rhoMax
             );
 
+            //- Add the given density correction to the density field.
+            //  Used to update the density field following pressure solution
+            virtual void correctRho
+            (
+                const volScalarField& deltaRho
+            );
+
             //- Density [kg/m^3] - uses current value of pressure
             virtual tmp<volScalarField> rho() const;
 
diff --git a/src/thermophysicalModels/basic/rhoThermo/rhoThermo.C b/src/thermophysicalModels/basic/rhoThermo/rhoThermo.C
index bbc1a64ba5..e7b70a4fc9 100644
--- a/src/thermophysicalModels/basic/rhoThermo/rhoThermo.C
+++ b/src/thermophysicalModels/basic/rhoThermo/rhoThermo.C
@@ -185,6 +185,12 @@ void Foam::rhoThermo::correctRho
     rho_ = min(rho_, rhoMax);
 }
 
+void Foam::rhoThermo::correctRho(const Foam::volScalarField& deltaRho)
+{
+    rho_ += deltaRho;
+}
+
+
 const Foam::volScalarField& Foam::rhoThermo::psi() const
 {
     return psi_;
diff --git a/src/thermophysicalModels/basic/rhoThermo/rhoThermo.H b/src/thermophysicalModels/basic/rhoThermo/rhoThermo.H
index fd8bb6b81c..2728cf02e7 100644
--- a/src/thermophysicalModels/basic/rhoThermo/rhoThermo.H
+++ b/src/thermophysicalModels/basic/rhoThermo/rhoThermo.H
@@ -146,6 +146,10 @@ public:
                 const dimensionedScalar& rhoMax
             );
 
+            //- Add the given density correction to the density field.
+            //  Used to update the density field following pressure solution
+            virtual void correctRho(const volScalarField& deltaRho);
+
             //- Compressibility [s^2/m^2]
             virtual const volScalarField& psi() const;