diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/createFields.H b/applications/solvers/heatTransfer/buoyantPimpleFoam/createFields.H
index 8c6402f8f48366795ce1288ea6cb2aa0f1ffc5d8..af2b8926aa5498a525c50ca47741253cb003ed02 100644
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/createFields.H
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/createFields.H
@@ -90,3 +90,27 @@ volScalarField dpdt
 
 Info<< "Creating field kinetic energy K\n" << endl;
 volScalarField K("K", 0.5*magSqr(U));
+
+label pRefCell = 0;
+scalar pRefValue = 0.0;
+
+if (p_rgh.needReference())
+{
+    setRefCell
+    (
+        p,
+        p_rgh,
+        pimple.dict(),
+        pRefCell,
+        pRefValue
+    );
+
+    p += dimensionedScalar
+    (
+        "p",
+        p.dimensions(),
+        pRefValue - getRefCellValue(p, pRefCell)
+    );
+}
+
+dimensionedScalar initialMass("initialMass", fvc::domainIntegrate(rho));
diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H
index e2f99671a356a93fed292c9b3e83ce1d3d211ee6..76b4925917a3b5d61f1fb98403aa6a29b1161f32 100644
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H
@@ -1,4 +1,9 @@
 {
+    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
@@ -36,19 +41,27 @@
         )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
     );
 
-    fvScalarMatrix p_rghDDtEqn
+    tmp<fvScalarMatrix> p_rghDDtEqn
     (
-        fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
-      + fvc::div(phiHbyA)
-      ==
-        fvOptions(psi, p_rgh, rho.name())
+        new fvScalarMatrix(p_rgh, dimMass/dimTime)
     );
 
+    if (compressible)
+    {
+        p_rghDDtEqn =
+        (
+            fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
+        ==
+            fvOptions(psi, p_rgh, rho.name())
+        );
+    }
+
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix p_rghEqn
         (
-            p_rghDDtEqn
+            p_rghDDtEqn()
+          + fvc::div(phiHbyA)
           - fvm::laplacian(rAUf, p_rgh)
         );
 
@@ -81,6 +94,32 @@
         dpdt = fvc::ddt(p);
     }
 
-    #include "rhoEqn.H"
+    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;
 }
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H
index 002faddb59433b7bcf0c52e159d65c64a3cbc040..cd42cc00fdfd49722907ed4d585f06dfbf130e3e 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H
@@ -36,12 +36,21 @@
         )/(mesh.magSf().boundaryField()*rhorAUf.boundaryField())
     );
 
+    tmp<fvScalarMatrix> p_rghDDtEqn
+    (
+        new fvScalarMatrix(p_rgh, dimMass/dimTime)
+    );
+
     {
-        fvScalarMatrix p_rghDDtEqn
-        (
-            fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
-          + fvc::div(phiHbyA)
-        );
+        if (compressible)
+        {
+            p_rghDDtEqn =
+            (
+                fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
+             ==
+                fvOptions(psi, p_rgh, rho.name())
+            );
+        }
 
         // Thermodynamic density needs to be updated by psi*d(p) after the
         // pressure solution - done in 2 parts. Part 1:
@@ -52,6 +61,7 @@
             fvScalarMatrix p_rghEqn
             (
                 p_rghDDtEqn
+              + fvc::div(phiHbyA)
               - fvm::laplacian(rhorAUf, p_rgh)
             );
 
@@ -93,8 +103,11 @@
         dpdt = fvc::ddt(p);
     }
 
-    // Solve continuity
-    #include "rhoEqn.H"
+    if (compressible)
+    {
+        // Solve continuity
+        #include "rhoEqn.H"
+    }
 
     // Update continuity errors
     #include "compressibleContinuityErrors.H"
diff --git a/tutorials/heatTransfer/buoyantPimpleFoam/hotRoom/system/fvSolution b/tutorials/heatTransfer/buoyantPimpleFoam/hotRoom/system/fvSolution
index ef63b73506cfaee6e247870251cef7070abe8ee5..78d0b84613401eb84234fe8e18e6c2acceab4ea5 100644
--- a/tutorials/heatTransfer/buoyantPimpleFoam/hotRoom/system/fvSolution
+++ b/tutorials/heatTransfer/buoyantPimpleFoam/hotRoom/system/fvSolution
@@ -60,6 +60,8 @@ PIMPLE
     nOuterCorrectors 1;
     nCorrectors     2;
     nNonOrthogonalCorrectors 0;
+    pRefCell        0;
+    pRefValue       1e5;
 }