diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/TEqn.H b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/TEqn.H
index dbfc61739f7d21e2d095ae55181831345fc05a75..9a835792a4ec379bccfe2d933d34999c35b0a19a 100644
--- a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/TEqn.H
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/TEqn.H
@@ -12,7 +12,7 @@
     );
 
     TEqn.relax();
-    TEqn.solve();
+    TEqn.solve(mesh.solver(T.select(finalIter)));
 
     rhok = 1.0 - beta*(T - TRef);
 }
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/UEqn.H b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/UEqn.H
index 35387f4ff40468426b6ebc31a6fb517553759fe4..20a05e5cd448366b05f8f9de882545190c31e8a0 100644
--- a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/UEqn.H
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/UEqn.H
@@ -18,9 +18,10 @@
             fvc::reconstruct
             (
                 (
-                    fvc::interpolate(rhok)*(g & mesh.Sf())
-                  - fvc::snGrad(p)*mesh.magSf()
-                )
-            )
+                  - ghf*fvc::snGrad(rhok)
+                  - fvc::snGrad(p_rgh)
+                )*mesh.magSf()
+            ),
+            mesh.solver(U.select(finalIter))
         );
     }
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/buoyantBoussinesqPimpleFoam.C b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/buoyantBoussinesqPimpleFoam.C
index ebf68d409c9965a3333fe337090125190ba60b83..54519906a481cc48c873364a7036aab81c60dd96 100644
--- a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/buoyantBoussinesqPimpleFoam.C
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/buoyantBoussinesqPimpleFoam.C
@@ -87,7 +87,7 @@ int main(int argc, char *argv[])
 
             if (nOuterCorr != 1)
             {
-                p.storePrevIter();
+                p_rgh.storePrevIter();
             }
 
             #include "UEqn.H"
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/createFields.H b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/createFields.H
index 23d20cfa96d16d4bdecfbd2778d46f39b6a8cd67..342af079d8d3f2711e2ab1fb6c841ebc63e82559 100644
--- a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/createFields.H
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/createFields.H
@@ -14,12 +14,12 @@
         mesh
     );
 
-    Info<< "Reading field p\n" << endl;
-    volScalarField p
+    Info<< "Reading field p_rgh\n" << endl;
+    volScalarField p_rgh
     (
         IOobject
         (
-            "p",
+            "p_rgh",
             runTime.timeName(),
             mesh,
             IOobject::MUST_READ,
@@ -52,6 +52,18 @@
         incompressible::RASModel::New(U, phi, laminarTransport)
     );
 
+    // Kinematic density for buoyancy force
+    volScalarField rhok
+    (
+        IOobject
+        (
+            "rhok",
+            runTime.timeName(),
+            mesh
+        ),
+        1.0 - beta*(T - TRef)
+    );
+
     // kinematic turbulent thermal thermal conductivity m2/s
     Info<< "Reading field kappat\n" << endl;
     volScalarField kappat
@@ -67,25 +79,41 @@
         mesh
     );
 
+    Info<< "Calculating field g.h\n" << endl;
+    volScalarField gh("gh", g & mesh.C());
+    surfaceScalarField ghf("ghf", g & mesh.Cf());
+
+    volScalarField p
+    (
+        IOobject
+        (
+            "p",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        p_rgh + rhok*gh
+    );
+
     label pRefCell = 0;
     scalar pRefValue = 0.0;
     setRefCell
     (
         p,
+        p_rgh,
         mesh.solutionDict().subDict("PIMPLE"),
         pRefCell,
         pRefValue
     );
 
-
-    // Kinematic density for buoyancy force
-    volScalarField rhok
-    (
-        IOobject
+    if (p_rgh.needReference())
+    {
+        p += dimensionedScalar
         (
-            "rhok",
-            runTime.timeName(),
-            mesh
-        ),
-        1.0 - beta*(T - TRef)
-    );
+            "p",
+            p.dimensions(),
+            pRefValue - getRefCellValue(p, pRefCell)
+        );
+    }
+
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/pEqn.H
index 21be033f9bc5a209848fb1198fddfd4f8cbf9cfa..60828e18dc544a08dba7a5827577d2cf0df76ac5 100644
--- a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/pEqn.H
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/pEqn.H
@@ -7,22 +7,23 @@
     phi = (fvc::interpolate(U) & mesh.Sf())
         + fvc::ddtPhiCorr(rUA, U, phi);
 
-    surfaceScalarField buoyancyPhi =
-        rUAf*fvc::interpolate(rhok)*(g & mesh.Sf());
-    phi += buoyancyPhi;
+    surfaceScalarField buoyancyPhi = rUAf*ghf*fvc::snGrad(rhok)*mesh.magSf();
+    phi -= buoyancyPhi;
 
     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
-        fvScalarMatrix pEqn
+        fvScalarMatrix p_rghEqn
         (
-            fvm::laplacian(rUAf, p) == fvc::div(phi)
+            fvm::laplacian(rUAf, p_rgh) == fvc::div(phi)
         );
 
-        pEqn.solve
+        p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
+
+        p_rghEqn.solve
         (
             mesh.solver
             (
-                p.select
+                p_rgh.select
                 (
                     (
                         finalIter
@@ -36,17 +37,30 @@
         if (nonOrth == nNonOrthCorr)
         {
             // Calculate the conservative fluxes
-            phi -= pEqn.flux();
+            phi -= p_rghEqn.flux();
 
             // Explicitly relax pressure for momentum corrector
-            p.relax();
+            p_rgh.relax();
 
             // Correct the momentum source with the pressure gradient flux
             // calculated from the relaxed pressure
-            U += rUA*fvc::reconstruct((buoyancyPhi - pEqn.flux())/rUAf);
+            U -= rUA*fvc::reconstruct((buoyancyPhi + p_rghEqn.flux())/rUAf);
             U.correctBoundaryConditions();
         }
     }
 
     #include "continuityErrs.H"
+
+    p = p_rgh + rhok*gh;
+
+    if (p_rgh.needReference())
+    {
+        p += dimensionedScalar
+        (
+            "p",
+            p.dimensions(),
+            pRefValue - getRefCellValue(p, pRefCell)
+        );
+        p_rgh = p - rhok*gh;
+    }
 }
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/createFields.H b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/createFields.H
index 76cc4da8bace0256e6abe6285a12b71deacb86a9..477a3228331f2f98218d1f764ed8e84f1a72133f 100644
--- a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/createFields.H
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/createFields.H
@@ -96,29 +96,23 @@
         p_rgh + rhok*gh
     );
 
-    label p_rghRefCell = 0;
-    scalar p_rghRefValue = 0.0;
+    label pRefCell = 0;
+    scalar pRefValue = 0.0;
     setRefCell
     (
+        p,
         p_rgh,
         mesh.solutionDict().subDict("SIMPLE"),
-        p_rghRefCell,
-        p_rghRefValue
+        pRefCell,
+        pRefValue
     );
 
-    scalar pRefValue = 0.0;
-
     if (p_rgh.needReference())
     {
-        pRefValue = readScalar
-        (
-            mesh.solutionDict().subDict("SIMPLE").lookup("pRefValue")
-        );
-
         p += dimensionedScalar
         (
             "p",
             p.dimensions(),
-            pRefValue - getRefCellValue(p, p_rghRefCell)
+            pRefValue - getRefCellValue(p, pRefCell)
         );
     }
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pEqn.H
index 5664bb2154b13c134cd17ced9962e661f8ba8ff6..50149ed36013ce28d445149b454076b546c1ee0b 100644
--- a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pEqn.H
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pEqn.H
@@ -18,17 +18,9 @@
             fvm::laplacian(rUAf, p_rgh) == fvc::div(phi)
         );
 
-        p_rghEqn.setReference(p_rghRefCell, p_rghRefValue);
+        p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
 
-        // retain the residual from the first iteration
-        if (nonOrth == 0)
-        {
-            p_rghEqn.solve();
-        }
-        else
-        {
-            p_rghEqn.solve();
-        }
+        p_rghEqn.solve();
 
         if (nonOrth == nNonOrthCorr)
         {
@@ -55,7 +47,8 @@
         (
             "p",
             p.dimensions(),
-            pRefValue - getRefCellValue(p, p_rghRefCell)
+            pRefValue - getRefCellValue(p, pRefCell)
         );
+        p_rgh = p - rhok*gh;
     }
 }
diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/UEqn.H b/applications/solvers/heatTransfer/buoyantPimpleFoam/UEqn.H
index d4878d063dab6ee8ed4243c61d529832d3dd9c6a..8c6a3f7671aafa72e1bd3dc53cf8644d1345a642 100644
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/UEqn.H
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/UEqn.H
@@ -17,8 +17,11 @@
          ==
             fvc::reconstruct
             (
-                fvc::interpolate(rho)*(g & mesh.Sf())
-              - fvc::snGrad(p)*mesh.magSf()
-            )
+                (
+                  - ghf*fvc::snGrad(rho)
+                  - fvc::snGrad(p_rgh)
+                )*mesh.magSf()
+            ),
+            mesh.solver(U.select(finalIter))
         );
     }
diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.C b/applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.C
index 0075c1e3edb5736b8ba0252f94b6149b4b8192d3..cb4c4d34f68fc5e40759e5142d6ecffc9eb86e3e 100644
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.C
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.C
@@ -80,7 +80,7 @@ int main(int argc, char *argv[])
 
             if (nOuterCorr != 1)
             {
-                p.storePrevIter();
+                p_rgh.storePrevIter();
             }
 
             #include "UEqn.H"
diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/createFields.H b/applications/solvers/heatTransfer/buoyantPimpleFoam/createFields.H
index b8ac5595e47cec77aa257c54a1710f91b08f20f7..e39d0bb80333b1cb16731cb0d63356c100aaf093 100644
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/createFields.H
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/createFields.H
@@ -53,15 +53,30 @@
         )
     );
 
+    Info<< "Calculating field g.h\n" << endl;
+    volScalarField gh("gh", g & mesh.C());
+    surfaceScalarField ghf("ghf", g & mesh.Cf());
+
+    Info<< "Reading field p_rgh\n" << endl;
+    volScalarField p_rgh
+    (
+        IOobject
+        (
+            "p_rgh",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    // Force p_rgh to be consistent with p
+    p_rgh = p - rho*gh;
+
     Info<< "Creating field DpDt\n" << endl;
     volScalarField DpDt
     (
         "DpDt",
         fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p)
     );
-
-    thermo.correct();
-
-    dimensionedScalar initialMass = fvc::domainIntegrate(rho);
-
-    dimensionedScalar totalVolume = sum(mesh.V());
diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/hEqn.H b/applications/solvers/heatTransfer/buoyantPimpleFoam/hEqn.H
index 3125cc3ffa86ce120e7dbbf774c9b46941105418..94537508b3725cc562118f196e2cca0de6664651 100644
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/hEqn.H
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/hEqn.H
@@ -9,7 +9,7 @@
     );
 
     hEqn.relax();
-    hEqn.solve();
+    hEqn.solve(mesh.solver(h.select(finalIter)));
 
     thermo.correct();
 }
diff --git a/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H
index a1c3dbfed59fb52e459c460e2c439290499f5300..e625f122a3b86d8fd2b076bb3a066b98c531d0fa 100644
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H
@@ -1,11 +1,9 @@
 {
-    bool closedVolume = p.needReference();
-
     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;
+    thermo.rho() -= psi*p_rgh;
 
     volScalarField rUA = 1.0/UEqn.A();
     surfaceScalarField rhorUAf("(rho*(1|A(U)))", fvc::interpolate(rho*rUA));
@@ -18,24 +16,23 @@
       + fvc::ddtPhiCorr(rUA, rho, U, phi)
     );
 
-    surfaceScalarField buoyancyPhi =
-        rhorUAf*fvc::interpolate(rho)*(g & mesh.Sf());
+    surfaceScalarField buoyancyPhi = -rhorUAf*ghf*fvc::snGrad(rho)*mesh.magSf();
     phi += buoyancyPhi;
 
     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
-        fvScalarMatrix pEqn
+        fvScalarMatrix p_rghEqn
         (
-            fvc::ddt(rho) + psi*correction(fvm::ddt(p))
+            fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
           + fvc::div(phi)
-          - fvm::laplacian(rhorUAf, p)
+          - fvm::laplacian(rhorUAf, p_rgh)
         );
 
-        pEqn.solve
+        p_rghEqn.solve
         (
             mesh.solver
             (
-                p.select
+                p_rgh.select
                 (
                     (
                         finalIter
@@ -49,34 +46,25 @@
         if (nonOrth == nNonOrthCorr)
         {
             // Calculate the conservative fluxes
-            phi += pEqn.flux();
+            phi += p_rghEqn.flux();
 
             // Explicitly relax pressure for momentum corrector
-            p.relax();
+            p_rgh.relax();
 
             // Correct the momentum source with the pressure gradient flux
             // calculated from the relaxed pressure
-            U += rUA*fvc::reconstruct((buoyancyPhi + pEqn.flux())/rhorUAf);
+            U += rUA*fvc::reconstruct((buoyancyPhi + p_rghEqn.flux())/rhorUAf);
             U.correctBoundaryConditions();
         }
     }
 
+    p = p_rgh + rho*gh;
+
     // Second part of thermodynamic density update
-    thermo.rho() += psi*p;
+    thermo.rho() += psi*p_rgh;
 
     DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
 
     #include "rhoEqn.H"
     #include "compressibleContinuityErrs.H"
-
-    // 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);
-        thermo.rho() = psi*p;
-        rho += (initialMass - fvc::domainIntegrate(rho))/totalVolume;
-    }
 }
diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C b/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C
index 06fa5d12a02e9e932ac71b099a956a085c420897..64cc110cecee14e2cde99619620b568e08cd5405 100644
--- a/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C
+++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C
@@ -62,10 +62,7 @@ int main(int argc, char *argv[])
         {
             #include "UEqn.H"
             #include "hEqn.H"
-            for (int i=0; i<3; i++)
-            {
             #include "pEqn.H"
-            }
         }
 
         turbulence->correct();
diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H b/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H
index cf231503323ee8eb7661bfb42e36bb08c5b7d830..d6fa9acee96ce088fe9c4b08bba8ee8b2f7de786 100644
--- a/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H
+++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H
@@ -23,20 +23,6 @@
     volScalarField& h = thermo.h();
     const volScalarField& psi = thermo.psi();
 
-    Info<< "Reading field p_rgh\n" << endl;
-    volScalarField p_rgh
-    (
-        IOobject
-        (
-            "p_rgh",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
     Info<< "Reading field U\n" << endl;
     volVectorField U
     (
@@ -53,7 +39,6 @@
 
     #include "compressibleCreatePhi.H"
 
-
     Info<< "Creating turbulence model\n" << endl;
     autoPtr<compressible::RASModel> turbulence
     (
@@ -66,40 +51,39 @@
         )
     );
 
+
     Info<< "Calculating field g.h\n" << endl;
     volScalarField gh("gh", g & mesh.C());
     surfaceScalarField ghf("ghf", g & mesh.Cf());
 
-    p = p_rgh + rho*gh;
-    thermo.correct();
-    rho = thermo.rho();
+    Info<< "Reading field p_rgh\n" << endl;
+    volScalarField p_rgh
+    (
+        IOobject
+        (
+            "p_rgh",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    // Force p_rgh to be consistent with p
     p_rgh = p - rho*gh;
 
-    label p_rghRefCell = 0;
-    scalar p_rghRefValue = 0.0;
+
+    label pRefCell = 0;
+    scalar pRefValue = 0.0;
     setRefCell
     (
+        p,
         p_rgh,
         mesh.solutionDict().subDict("SIMPLE"),
-        p_rghRefCell,
-        p_rghRefValue
+        pRefCell,
+        pRefValue
     );
 
-    scalar pRefValue = 0.0;
-
-    if (p_rgh.needReference())
-    {
-        pRefValue = readScalar
-        (
-            mesh.solutionDict().subDict("SIMPLE").lookup("pRefValue")
-        );
-
-        p += dimensionedScalar
-        (
-            "p",
-            p.dimensions(),
-            pRefValue - getRefCellValue(p, p_rghRefCell)
-        );
-    }
-
     dimensionedScalar initialMass = fvc::domainIntegrate(rho);
+    dimensionedScalar totalVolume = sum(mesh.V());
diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H
index 3088b42c697ce467ec0c6448654fba2445b26eef..f1a62d4770a6581921e1b9a3479d3d2dd4a6c10d 100644
--- a/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H
+++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H
@@ -1,11 +1,12 @@
 {
     rho = thermo.rho();
+    rho.relax();
 
     volScalarField rUA = 1.0/UEqn().A();
     surfaceScalarField rhorUAf("(rho*(1|A(U)))", fvc::interpolate(rho*rUA));
 
     U = rUA*UEqn().H();
-    //UEqn.clear();
+    UEqn.clear();
 
     phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf());
     bool closedVolume = adjustPhi(phi, U, p_rgh);
@@ -20,7 +21,7 @@
             fvm::laplacian(rhorUAf, p_rgh) == fvc::div(phi)
         );
 
-        p_rghEqn.setReference(p_rghRefCell, p_rghRefValue);
+        p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
         p_rghEqn.solve();
 
         if (nonOrth == nNonOrthCorr)
@@ -42,13 +43,13 @@
 
     p = p_rgh + rho*gh;
 
-    // For closed-volume cases adjust the pressure and density levels
+    // For closed-volume cases adjust the pressure level
     // to obey overall mass continuity
     if (closedVolume)
     {
         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/buoyantSimpleRadiationFoam/buoyantSimpleRadiationFoam.C b/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/buoyantSimpleRadiationFoam.C
index 0dbe80c67c0d408339bead48b579580d79bb6523..6c35536333ace6a10582d022ef26fecb80b351d8 100644
--- a/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/buoyantSimpleRadiationFoam.C
+++ b/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/buoyantSimpleRadiationFoam.C
@@ -58,7 +58,7 @@ int main(int argc, char *argv[])
 
         #include "readSIMPLEControls.H"
 
-        p.storePrevIter();
+        p_rgh.storePrevIter();
         rho.storePrevIter();
 
         // Pressure-velocity SIMPLE corrector
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C
index a81b0faaf3eb60890151e8d3a66af70de4d4bff5..36b3f1c3b0f5357b4ae28cca263b09d063c1ea18 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C
@@ -93,6 +93,8 @@ int main(int argc, char *argv[])
         // --- PIMPLE loop
         for (int oCorr=0; oCorr<nOuterCorr; oCorr++)
         {
+            bool finalIter = oCorr == nOuterCorr-1;
+
             forAll(fluidRegions, i)
             {
                 Info<< "\nSolving for fluid region "
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/UEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/UEqn.H
index d5ae07ff9be5ed6c95ec926f5a07c7d7ecde78a0..68b29fd0c2babc531af4fe730abd8a08c3ecf5ff 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/UEqn.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/UEqn.H
@@ -13,7 +13,9 @@
      ==
         fvc::reconstruct
         (
-            fvc::interpolate(rho)*(g & mesh.Sf())
-          - fvc::snGrad(p)*mesh.magSf()
+            (
+              - ghf*fvc::snGrad(rho)
+              - fvc::snGrad(p_rgh)
+            )*mesh.magSf()
         )
     );
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/createFluidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/createFluidFields.H
index 2d2ccf9a5eb4d443059ba97985a77b229407073e..8d8cad34777bdb49b8d539a0c4de8be4283fe103 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/createFluidFields.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/createFluidFields.H
@@ -6,12 +6,16 @@
     PtrList<surfaceScalarField> phiFluid(fluidRegions.size());
     PtrList<uniformDimensionedVectorField> gFluid(fluidRegions.size());
     PtrList<compressible::turbulenceModel> turbulence(fluidRegions.size());
-    PtrList<volScalarField> DpDtf(fluidRegions.size());
+    PtrList<volScalarField> p_rghFluid(fluidRegions.size());
+    PtrList<volScalarField> ghFluid(fluidRegions.size());
+    PtrList<surfaceScalarField> ghfFluid(fluidRegions.size());
 
     List<scalar> initialMassFluid(fluidRegions.size());
     List<label> pRefCellFluid(fluidRegions.size(),0);
     List<scalar> pRefValueFluid(fluidRegions.size(),0.0);
 
+    PtrList<dimensionedScalar> rhoMax(fluidRegions.size());
+    PtrList<dimensionedScalar> rhoMin(fluidRegions.size());
 
     // Populate fluid field pointer lists
     forAll(fluidRegions, i)
@@ -130,15 +134,68 @@
             ).ptr()
         );
 
+        Info<< "    Adding to ghFluid\n" << endl;
+        ghFluid.set
+        (
+            i,
+            new volScalarField("gh", gFluid[i] & fluidRegions[i].C())
+        );
+
+        Info<< "    Adding to ghfFluid\n" << endl;
+        ghfFluid.set
+        (
+            i,
+            new surfaceScalarField("ghf", gFluid[i] & fluidRegions[i].Cf())
+        );
+
+        p_rghFluid.set
+        (
+            i,
+            new volScalarField
+            (
+                IOobject
+                (
+                    "p_rgh",
+                    runTime.timeName(),
+                    fluidRegions[i],
+                    IOobject::MUST_READ,
+                    IOobject::AUTO_WRITE
+                ),
+                fluidRegions[i]
+            )
+        );
+
+        // Force p_rgh to be consistent with p
+        p_rghFluid[i] = thermoFluid[i].p() - rhoFluid[i]*ghFluid[i];
+
         initialMassFluid[i] = fvc::domainIntegrate(rhoFluid[i]).value();
 
         setRefCell
         (
             thermoFluid[i].p(),
+            p_rghFluid[i],
             fluidRegions[i].solutionDict().subDict("SIMPLE"),
             pRefCellFluid[i],
             pRefValueFluid[i]
         );
+
+        rhoMax.set
+        (
+            i,
+            new dimensionedScalar
+            (
+                fluidRegions[i].solutionDict().subDict("SIMPLE").lookup("rhoMax")
+            )
+        );
+
+        rhoMin.set
+        (
+            i,
+            new dimensionedScalar
+            (
+                fluidRegions[i].solutionDict().subDict("SIMPLE").lookup("rhoMin")
+            )
+        );
     }
 
 
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/hEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/hEqn.H
index cef0122928dae7e51242c2a77a1afe422421151e..7869046fa0213e4b92d30119fa40f2149489095b 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/hEqn.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/hEqn.H
@@ -5,8 +5,8 @@
       - fvm::Sp(fvc::div(phi), h)
       - fvm::laplacian(turb.alphaEff(), h)
      ==
-        fvc::div(phi/fvc::interpolate(rho)*fvc::interpolate(p))
-      - p*fvc::div(phi/fvc::interpolate(rho))
+        fvc::div(phi/fvc::interpolate(rho), rho/psi, "div(U,p)")
+      - (rho/psi)*fvc::div(phi/fvc::interpolate(rho))
     );
 
     hEqn.relax();
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H
index 9462dc2c39debaeb2fef992a2f42a02a1ea2fc82..317d377f1bd376021c3450cfee07f8ec085d7bb4 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H
@@ -1,7 +1,9 @@
 {
     // From buoyantSimpleFoam
-
     rho = thermo.rho();
+    rho = max(rho, rhoMin[i]);
+    rho = min(rho, rhoMax[i]);
+    rho.relax();
 
     volScalarField rUA = 1.0/UEqn().A();
     surfaceScalarField rhorUAf("(rho*(1|A(U)))", fvc::interpolate(rho*rUA));
@@ -10,59 +12,54 @@
     UEqn.clear();
 
     phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf());
-    bool closedVolume = adjustPhi(phi, U, p);
+    bool closedVolume = adjustPhi(phi, U, p_rgh);
 
-    surfaceScalarField buoyancyPhi =
-        rhorUAf*fvc::interpolate(rho)*(g & mesh.Sf());
-    phi += buoyancyPhi;
+    surfaceScalarField buoyancyPhi = rhorUAf*ghf*fvc::snGrad(rho)*mesh.magSf();
+    phi -= buoyancyPhi;
 
     // Solve pressure
     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
-        fvScalarMatrix pEqn
+        fvScalarMatrix p_rghEqn
         (
-            fvm::laplacian(rhorUAf, p) == fvc::div(phi)
+            fvm::laplacian(rhorUAf, p_rgh) == fvc::div(phi)
         );
 
-        pEqn.setReference(pRefCell, pRefValue);
+        p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
 
-        // retain the residual from the first iteration
-        if (nonOrth == 0)
-        {
-            pEqn.solve();
-        }
-        else
-        {
-            pEqn.solve();
-        }
+        p_rghEqn.solve();
 
         if (nonOrth == nNonOrthCorr)
         {
-            // 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);
-            }
-
             // Calculate the conservative fluxes
-            phi -= pEqn.flux();
+            phi -= p_rghEqn.flux();
 
             // Explicitly relax pressure for momentum corrector
-            p.relax();
+            p_rgh.relax();
 
             // Correct the momentum source with the pressure gradient flux
             // calculated from the relaxed pressure
-            U += rUA*fvc::reconstruct((buoyancyPhi - pEqn.flux())/rhorUAf);
+            U -= rUA*fvc::reconstruct((buoyancyPhi + p_rghEqn.flux())/rhorUAf);
             U.correctBoundaryConditions();
         }
     }
 
+    p = p_rgh + rho*gh;
 
     #include "continuityErrs.H"
 
+    // For closed-volume cases adjust the pressure level
+    // to obey overall mass continuity
+    if (closedVolume)
+    {
+        p += (initialMass - fvc::domainIntegrate(psi*p))
+            /fvc::domainIntegrate(psi);
+        p_rgh = p - rho*gh;
+    }
+
     rho = thermo.rho();
+    rho = max(rho, rhoMin[i]);
+    rho = min(rho, rhoMax[i]);
     rho.relax();
 
     Info<< "Min/max rho:" << min(rho).value() << ' '
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/setRegionFluidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/setRegionFluidFields.H
index 8491056ea325e7f20258da63a16e60c9ab2130e2..cac7750e972927250bbfd7f9aec2c5131cacf3de 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/setRegionFluidFields.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/setRegionFluidFields.H
@@ -5,7 +5,6 @@
     volScalarField& K = KFluid[i];
     volVectorField& U = UFluid[i];
     surfaceScalarField& phi = phiFluid[i];
-    const dimensionedVector& g = gFluid[i];
 
     compressible::turbulenceModel& turb = turbulence[i];
 
@@ -22,3 +21,7 @@
 
     const label pRefCell = pRefCellFluid[i];
     const scalar pRefValue = pRefValueFluid[i];
+
+    volScalarField& p_rgh = p_rghFluid[i];
+    const volScalarField& gh = ghFluid[i];
+    const surfaceScalarField& ghf = ghfFluid[i];
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/solveFluid.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/solveFluid.H
index 5198941366285201708ca870b3763390a07c5db8..2b6de83ca3dadda1628be14ca5862743135d57e5 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/solveFluid.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/solveFluid.H
@@ -1,6 +1,6 @@
 //  Pressure-velocity SIMPLE corrector
 
-    p.storePrevIter();
+    p_rgh.storePrevIter();
     rho.storePrevIter();
     {
         #include "UEqn.H"
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/UEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/UEqn.H
index 65467f80864609f9225e9d058b9568890c500988..119eeb4fd7d314bbd353934519093edb327e39c4 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/UEqn.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/UEqn.H
@@ -16,8 +16,11 @@
          ==
             fvc::reconstruct
             (
-                fvc::interpolate(rho)*(g & mesh.Sf())
-              - fvc::snGrad(p)*mesh.magSf()
-            )
+                (
+                  - ghf*fvc::snGrad(rho)
+                  - fvc::snGrad(p_rgh)
+                )*mesh.magSf()
+            ),
+            mesh.solver(U.select(finalIter))
         );
     }
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H
index 1826a77217e86ac89941941c2d8e452d783e94fb..b0a7f95912fa92648880d15a94dc49aacb6386e4 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H
@@ -6,6 +6,9 @@
     PtrList<surfaceScalarField> phiFluid(fluidRegions.size());
     PtrList<uniformDimensionedVectorField> gFluid(fluidRegions.size());
     PtrList<compressible::turbulenceModel> turbulence(fluidRegions.size());
+    PtrList<volScalarField> p_rghFluid(fluidRegions.size());
+    PtrList<volScalarField> ghFluid(fluidRegions.size());
+    PtrList<surfaceScalarField> ghfFluid(fluidRegions.size());
     PtrList<volScalarField> DpDtFluid(fluidRegions.size());
 
     List<scalar> initialMassFluid(fluidRegions.size());
@@ -129,6 +132,42 @@
             ).ptr()
         );
 
+        Info<< "    Adding to ghFluid\n" << endl;
+        ghFluid.set
+        (
+            i,
+            new volScalarField("gh", gFluid[i] & fluidRegions[i].C())
+        );
+
+        Info<< "    Adding to ghfFluid\n" << endl;
+        ghfFluid.set
+        (
+            i,
+            new surfaceScalarField("ghf", gFluid[i] & fluidRegions[i].Cf())
+        );
+
+        p_rghFluid.set
+        (
+            i,
+            new volScalarField
+            (
+                IOobject
+                (
+                    "p_rgh",
+                    runTime.timeName(),
+                    fluidRegions[i],
+                    IOobject::MUST_READ,
+                    IOobject::AUTO_WRITE
+                ),
+                fluidRegions[i]
+            )
+        );
+
+        // Force p_rgh to be consistent with p
+        p_rghFluid[i] = thermoFluid[i].p() - rhoFluid[i]*ghFluid[i];
+
+        initialMassFluid[i] = fvc::domainIntegrate(rhoFluid[i]).value();
+
         Info<< "    Adding to DpDtFluid\n" << endl;
         DpDtFluid.set
         (
@@ -147,6 +186,4 @@
                 )
             )
         );
-
-        initialMassFluid[i] = fvc::domainIntegrate(rhoFluid[i]).value();
     }
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/hEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/hEqn.H
index e070537db2c701e27baafd6cf64ef015e3898100..16ba36f7dc71513a6204dde059fb1cfe8868062d 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/hEqn.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/hEqn.H
@@ -7,16 +7,9 @@
      ==
         DpDt
     );
-    if (oCorr == nOuterCorr-1)
-    {
-        hEqn.relax();
-        hEqn.solve(mesh.solver("hFinal"));
-    }
-    else
-    {
-        hEqn.relax();
-        hEqn.solve();
-    }
+
+    hEqn.relax();
+    hEqn.solve(mesh.solver(h.select(finalIter)));
 
     thermo.correct();
 
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H
index a264b68fe5eab2388d38ca1aaf75a3f416db2330..9bf3c4acc50b7f60fa54786de4d187c90d9f677a 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H
@@ -1,5 +1,5 @@
 {
-    bool closedVolume = p.needReference();
+    bool closedVolume = p_rgh.needReference();
 
     rho = thermo.rho();
 
@@ -17,34 +17,35 @@
         )
     );
 
-    phi = phiU + fvc::interpolate(rho)*(g & mesh.Sf())*rhorUAf;
+    phi = phiU - rhorUAf*ghf*fvc::snGrad(rho)*mesh.magSf();
 
     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
-        fvScalarMatrix pEqn
+        fvScalarMatrix p_rghEqn
         (
-            fvm::ddt(psi, p)
+            fvm::ddt(psi, p_rgh) + fvc::ddt(psi, rho)*gh
           + fvc::div(phi)
-          - fvm::laplacian(rhorUAf, p)
+          - fvm::laplacian(rhorUAf, p_rgh)
         );
 
-        if
+        p_rghEqn.solve
         (
-            oCorr == nOuterCorr-1
-         && corr == nCorr-1
-         && nonOrth == nNonOrthCorr
-        )
-        {
-            pEqn.solve(mesh.solver(p.name() + "Final"));
-        }
-        else
-        {
-            pEqn.solve(mesh.solver(p.name()));
-        }
+            mesh.solver
+            (
+                p_rgh.select
+                (
+                    (
+                        oCorr == nOuterCorr-1
+                     && corr == nCorr-1
+                     && nonOrth == nNonOrthCorr
+                    )
+                )
+            )
+        );
 
         if (nonOrth == nNonOrthCorr)
         {
-            phi += pEqn.flux();
+            phi += p_rghEqn.flux();
         }
     }
 
@@ -52,6 +53,8 @@
     U += rUA*fvc::reconstruct((phi - phiU)/rhorUAf);
     U.correctBoundaryConditions();
 
+    p = p_rgh + rho*gh;
+
     // Update pressure substantive derivative
     DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
 
@@ -65,9 +68,10 @@
     // to obey overall mass continuity
     if (closedVolume)
     {
-        p += (massIni - fvc::domainIntegrate(psi*p))
+        p += (initialMass - fvc::domainIntegrate(psi*p))
             /fvc::domainIntegrate(psi);
         rho = thermo.rho();
+        p_rgh = p - rho*gh;
     }
 
     // Update thermal conductivity
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/readFluidMultiRegionPISOControls.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/readFluidMultiRegionPISOControls.H
deleted file mode 100644
index 1f8cdb82a334c408603fed2745d4fbc224a7bfa4..0000000000000000000000000000000000000000
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/readFluidMultiRegionPISOControls.H
+++ /dev/null
@@ -1,17 +0,0 @@
-    const dictionary& piso = fluidRegions[i].solutionDict().subDict("PISO");
-
-    const int nOuterCorr =
-        piso.lookupOrDefault<int>("nOuterCorrectors", 1);
-
-    const int nCorr =
-        piso.lookupOrDefault<int>("nCorrectors", 1);
-
-    const int nNonOrthCorr =
-        piso.lookupOrDefault<int>("nNonOrthogonalCorrectors", 0);
-
-    const bool momentumPredictor =
-        piso.lookupOrDefault("momentumPredictor", true);
-
-    const bool transonic =
-        piso.lookupOrDefault("transonic", false);
-
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H
index 50a3011484689c25f6a26cce1e454bcbd20239b2..89aaec4737e014d64eb68bca51c0240c75231ba1 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H
@@ -1,11 +1,10 @@
-    const fvMesh& mesh = fluidRegions[i];
+    fvMesh& mesh = fluidRegions[i];
 
     basicPsiThermo& thermo = thermoFluid[i];
     volScalarField& rho = rhoFluid[i];
     volScalarField& K = KFluid[i];
     volVectorField& U = UFluid[i];
     surfaceScalarField& phi = phiFluid[i];
-    const dimensionedVector& g = gFluid[i];
 
     compressible::turbulenceModel& turb = turbulence[i];
     volScalarField& DpDt = DpDtFluid[i];
@@ -14,4 +13,13 @@
     const volScalarField& psi = thermo.psi();
     volScalarField& h = thermo.h();
 
-    const dimensionedScalar massIni("massIni", dimMass, initialMassFluid[i]);
+    volScalarField& p_rgh = p_rghFluid[i];
+    const volScalarField& gh = ghFluid[i];
+    const surfaceScalarField& ghf = ghfFluid[i];
+
+    const dimensionedScalar initialMass
+    (
+        "initialMass",
+        dimMass,
+        initialMassFluid[i]
+    );
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveFluid.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveFluid.H
index 86dd4344c15310f219bfb9bda406b22f2ffe7585..b36cf89d34e0821d116be4614d00d60708040da2 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveFluid.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/solveFluid.H
@@ -1,3 +1,8 @@
+if (finalIter)
+{
+    mesh.data::add("finalIteration", true);
+}
+
 if (oCorr == 0)
 {
     #include "rhoEqn.H"
@@ -16,3 +21,8 @@ for (int corr=0; corr<nCorr; corr++)
 turb.correct();
 
 rho = thermo.rho();
+
+if (finalIter)
+{
+    mesh.data::remove("finalIteration");
+}
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/storeOldFluidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/storeOldFluidFields.H
index f63e85458e25253209094f92995e0c3c91858fc5..94fd01a2baa30ffc9b16466ba1672f7e8a70dfd2 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/storeOldFluidFields.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/storeOldFluidFields.H
@@ -1,2 +1,2 @@
-    p.storePrevIter();
+    p_rgh.storePrevIter();
     rho.storePrevIter();
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/createSolidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/createSolidFields.H
index be13e670c326347f530dd991a4656e1dad22030f..837305659e6c088fa7048bb0e174a6da57fc418c 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/createSolidFields.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/createSolidFields.H
@@ -8,9 +8,5 @@
             << solidRegions[i].name() << nl << endl;
 
         Info<< "    Adding to thermos\n" << endl;
-        thermos.set
-        (
-            i,
-            basicSolidThermo::New(solidRegions[i])
-        );
+        thermos.set(i, basicSolidThermo::New(solidRegions[i]));
     }
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solveSolid.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solveSolid.H
index ce8b1d0f408d4ec033f88c37846a8c7ab80e61b0..05b26d1626d159c635fa1fd4b11d4348c1cc931c 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solveSolid.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solveSolid.H
@@ -1,3 +1,8 @@
+if (finalIter)
+{
+    mesh.data::add("finalIteration", true);
+}
+
 {
     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
@@ -7,8 +12,13 @@
           - fvm::laplacian(K, T)
         );
         TEqn().relax();
-        TEqn().solve();
+        TEqn().solve(mesh.solver(T.select(finalIter)));
     }
 
     Info<< "Min/max T:" << min(T) << ' ' << max(T) << endl;
 }
+
+if (finalIter)
+{
+    mesh.data::remove("finalIteration");
+}
diff --git a/applications/solvers/lagrangian/coalChemistryFoam/UEqn.H b/applications/solvers/lagrangian/coalChemistryFoam/UEqn.H
index b02a7c61e86523217ddf733841a2084e46cc2673..5851cdcf9510b6412b701c89090bf80a20297cb7 100644
--- a/applications/solvers/lagrangian/coalChemistryFoam/UEqn.H
+++ b/applications/solvers/lagrangian/coalChemistryFoam/UEqn.H
@@ -13,5 +13,5 @@
 
     if (momentumPredictor)
     {
-        solve(UEqn == -fvc::grad(p));
+        solve(UEqn == -fvc::grad(p), mesh.solver(U.select(finalIter)));
     }
diff --git a/applications/solvers/lagrangian/coalChemistryFoam/coalChemistryFoam.C b/applications/solvers/lagrangian/coalChemistryFoam/coalChemistryFoam.C
index 5f373cc2243d4542b6a804f29778182f4feb8ea1..f32be6214b2a213a0370145040b4682e15cea103 100644
--- a/applications/solvers/lagrangian/coalChemistryFoam/coalChemistryFoam.C
+++ b/applications/solvers/lagrangian/coalChemistryFoam/coalChemistryFoam.C
@@ -90,17 +90,28 @@ int main(int argc, char *argv[])
         #include "rhoEqn.H"
 
         // --- PIMPLE loop
-        for (int ocorr=1; ocorr<=nOuterCorr; ocorr++)
+        for (int oCorr=0; oCorr<nOuterCorr; oCorr++)
         {
+            bool finalIter = oCorr == nOuterCorr - 1;
+            if (finalIter)
+            {
+                mesh.data::add("finalIteration", true);
+            }
+
             #include "UEqn.H"
             #include "YEqn.H"
             #include "hsEqn.H"
 
             // --- PISO loop
-            for (int corr=1; corr<=nCorr; corr++)
+            for (int corr=0; corr<nCorr; corr++)
             {
                 #include "pEqn.H"
             }
+
+            if (finalIter)
+            {
+                mesh.data::remove("finalIteration");
+            }
         }
 
         turbulence->correct();
diff --git a/applications/solvers/lagrangian/coalChemistryFoam/hsEqn.H b/applications/solvers/lagrangian/coalChemistryFoam/hsEqn.H
index ced6b6400448a84182e7fba55a56bb06d18f14ab..40eb26709e14404510769f6aba0c77eeb62ef948 100644
--- a/applications/solvers/lagrangian/coalChemistryFoam/hsEqn.H
+++ b/applications/solvers/lagrangian/coalChemistryFoam/hsEqn.H
@@ -15,7 +15,7 @@
 
     hsEqn.relax();
 
-    hsEqn.solve();
+    hsEqn.solve(mesh.solver(hs.select(finalIter)));
 
     thermo.correct();
 
diff --git a/applications/solvers/lagrangian/coalChemistryFoam/pEqn.H b/applications/solvers/lagrangian/coalChemistryFoam/pEqn.H
index b31ad45ed4529b30d6d0dbc6c32c1b1e7f5ca6f0..d3efde0b349680e486d04e2097b64e7be9220efd 100644
--- a/applications/solvers/lagrangian/coalChemistryFoam/pEqn.H
+++ b/applications/solvers/lagrangian/coalChemistryFoam/pEqn.H
@@ -26,7 +26,20 @@ if (transonic)
             coalParcels.Srho()
         );
 
-        pEqn.solve();
+        pEqn.solve
+        (
+            mesh.solver
+            (
+                p.select
+                (
+                    (
+                        finalIter
+                     && corr == nCorr-1
+                     && nonOrth == nNonOrthCorr
+                    )
+                )
+            )
+        );
 
         if (nonOrth == nNonOrthCorr)
         {
@@ -54,7 +67,20 @@ else
             coalParcels.Srho()
         );
 
-        pEqn.solve();
+        pEqn.solve
+        (
+            mesh.solver
+            (
+                p.select
+                (
+                    (
+                        finalIter
+                     && corr == nCorr-1
+                     && nonOrth == nNonOrthCorr
+                    )
+                )
+            )
+        );
 
         if (nonOrth == nNonOrthCorr)
         {
diff --git a/tutorials/heatTransfer/buoyantBoussinesqPimpleFoam/hotRoom/0/T.org b/tutorials/heatTransfer/buoyantBoussinesqPimpleFoam/hotRoom/0/T.org
index 98cc191af7c9e6234287f45ebeab0cdea6ab8b62..77b58a1a00e6421c4d9c3f72a36b3d190e41ef7f 100644
--- a/tutorials/heatTransfer/buoyantBoussinesqPimpleFoam/hotRoom/0/T.org
+++ b/tutorials/heatTransfer/buoyantBoussinesqPimpleFoam/hotRoom/0/T.org
@@ -23,411 +23,7 @@ boundaryField
     floor
     {
         type            fixedValue;
-        value           nonuniform List<scalar> 
-400
-(
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-600
-600
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-600
-600
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-)
-;
+        value           uniform 300; 
     }
     ceiling
     {
diff --git a/tutorials/heatTransfer/buoyantBoussinesqPimpleFoam/hotRoom/0/p b/tutorials/heatTransfer/buoyantBoussinesqPimpleFoam/hotRoom/0/p
index b247cdcbaf97dc148fad58ea63f67e45e15d6d72..ed297deb8eed1122515449b3984ffe2cc178e15c 100644
--- a/tutorials/heatTransfer/buoyantBoussinesqPimpleFoam/hotRoom/0/p
+++ b/tutorials/heatTransfer/buoyantBoussinesqPimpleFoam/hotRoom/0/p
@@ -22,23 +22,20 @@ boundaryField
 {
     floor
     {
-        type            buoyantPressure;
-        rho             rhok;
-        value           uniform 0;
+        type            calculated;
+        value           $internalField;
     }
 
     ceiling
     {
-        type            buoyantPressure;
-        rho             rhok;
-        value           uniform 0;
+        type            calculated;
+        value           $internalField;
     }
 
     fixedWalls
     {
-        type            buoyantPressure;
-        rho             rhok;
-        value           uniform 0;
+        type            calculated;
+        value           $internalField;
     }
 }
 
diff --git a/tutorials/heatTransfer/buoyantBoussinesqPimpleFoam/hotRoom/0/p_rgh b/tutorials/heatTransfer/buoyantBoussinesqPimpleFoam/hotRoom/0/p_rgh
new file mode 100644
index 0000000000000000000000000000000000000000..270f01455806dd8d9937dc89b463b5e082fe41c9
--- /dev/null
+++ b/tutorials/heatTransfer/buoyantBoussinesqPimpleFoam/hotRoom/0/p_rgh
@@ -0,0 +1,45 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      p_rgh;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -2 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    floor
+    {
+        type            buoyantPressure;
+        rho             rhok;
+        value           uniform 0;
+    }
+
+    ceiling
+    {
+        type            buoyantPressure;
+        rho             rhok;
+        value           uniform 0;
+    }
+
+    fixedWalls
+    {
+        type            buoyantPressure;
+        rho             rhok;
+        value           uniform 0;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/heatTransfer/buoyantBoussinesqPimpleFoam/hotRoom/Allclean b/tutorials/heatTransfer/buoyantBoussinesqPimpleFoam/hotRoom/Allclean
index 25243dd42f8aff1c422994b6303af87d894424c3..fcc6e9c9cd8501447ff8d319f904432911616575 100755
--- a/tutorials/heatTransfer/buoyantBoussinesqPimpleFoam/hotRoom/Allclean
+++ b/tutorials/heatTransfer/buoyantBoussinesqPimpleFoam/hotRoom/Allclean
@@ -5,6 +5,6 @@ cd ${0%/*} || exit 1    # run from this directory
 . $WM_PROJECT_DIR/bin/tools/CleanFunctions
 
 cleanCase
-cp 0/T.org 0/T
+rm -f 0/T
 
 # ----------------------------------------------------------------- end-of-file
diff --git a/tutorials/heatTransfer/buoyantBoussinesqPimpleFoam/hotRoom/Allrun b/tutorials/heatTransfer/buoyantBoussinesqPimpleFoam/hotRoom/Allrun
index cb0be1e6e06e5244e2bddc492b8cb56b1ca04acd..800bb5731994d0f3a70b9f5ff7865f9659c1e4f7 100755
--- a/tutorials/heatTransfer/buoyantBoussinesqPimpleFoam/hotRoom/Allrun
+++ b/tutorials/heatTransfer/buoyantBoussinesqPimpleFoam/hotRoom/Allrun
@@ -8,6 +8,7 @@ application=`getApplication`
 
 compileApplication ../../buoyantPimpleFoam/hotRoom/setHotRoom
 runApplication blockMesh
+cp 0/T.org 0/T
 runApplication setHotRoom
 runApplication $application
 
diff --git a/tutorials/heatTransfer/buoyantBoussinesqPimpleFoam/hotRoom/system/fvSchemes b/tutorials/heatTransfer/buoyantBoussinesqPimpleFoam/hotRoom/system/fvSchemes
index fef36c6b1496b2039f68515cb96bd1a7667d41ad..764772f38ec2e85d1b5ed6c80bb02bad97ec37e6 100644
--- a/tutorials/heatTransfer/buoyantBoussinesqPimpleFoam/hotRoom/system/fvSchemes
+++ b/tutorials/heatTransfer/buoyantBoussinesqPimpleFoam/hotRoom/system/fvSchemes
@@ -40,12 +40,12 @@ divSchemes
 laplacianSchemes
 {
     default         none;
-    laplacian(nuEff,U) Gauss linear corrected;
-    laplacian((1|A(U)),p) Gauss linear corrected;
-    laplacian(kappaEff,T) Gauss linear corrected;
-    laplacian(DkEff,k) Gauss linear corrected;
-    laplacian(DepsilonEff,epsilon) Gauss linear corrected;
-    laplacian(DREff,R) Gauss linear corrected;
+    laplacian(nuEff,U) Gauss linear uncorrected;
+    laplacian((1|A(U)),p_rgh) Gauss linear uncorrected;
+    laplacian(kappaEff,T) Gauss linear uncorrected;
+    laplacian(DkEff,k) Gauss linear uncorrected;
+    laplacian(DepsilonEff,epsilon) Gauss linear uncorrected;
+    laplacian(DREff,R) Gauss linear uncorrected;
 }
 
 interpolationSchemes
@@ -55,13 +55,13 @@ interpolationSchemes
 
 snGradSchemes
 {
-    default         corrected;
+    default         uncorrected;
 }
 
 fluxRequired
 {
     default         no;
-    p               ;
+    p_rgh;
 }
 
 
diff --git a/tutorials/heatTransfer/buoyantBoussinesqPimpleFoam/hotRoom/system/fvSolution b/tutorials/heatTransfer/buoyantBoussinesqPimpleFoam/hotRoom/system/fvSolution
index 4c0b29ac70ad13ef23ad10f3daf5c0a46cb10d88..dccd8e082d3693f9f1ffeb6421104fc5fcb6996f 100644
--- a/tutorials/heatTransfer/buoyantBoussinesqPimpleFoam/hotRoom/system/fvSolution
+++ b/tutorials/heatTransfer/buoyantBoussinesqPimpleFoam/hotRoom/system/fvSolution
@@ -17,17 +17,17 @@ FoamFile
 
 solvers
 {
-    p
+    p_rgh
     {
         solver          PCG;
         preconditioner  DIC;
         tolerance       1e-8;
-        relTol          0.1;
+        relTol          0.01;
     }
 
-    pFinal
+    p_rghFinal
     {
-        $p;
+        $p_rgh;
         relTol          0;
     }
 
diff --git a/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/hotRoom/0/p b/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/hotRoom/0/p
index b247cdcbaf97dc148fad58ea63f67e45e15d6d72..ed297deb8eed1122515449b3984ffe2cc178e15c 100644
--- a/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/hotRoom/0/p
+++ b/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/hotRoom/0/p
@@ -22,23 +22,20 @@ boundaryField
 {
     floor
     {
-        type            buoyantPressure;
-        rho             rhok;
-        value           uniform 0;
+        type            calculated;
+        value           $internalField;
     }
 
     ceiling
     {
-        type            buoyantPressure;
-        rho             rhok;
-        value           uniform 0;
+        type            calculated;
+        value           $internalField;
     }
 
     fixedWalls
     {
-        type            buoyantPressure;
-        rho             rhok;
-        value           uniform 0;
+        type            calculated;
+        value           $internalField;
     }
 }
 
diff --git a/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/hotRoom/Allclean b/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/hotRoom/Allclean
index 25243dd42f8aff1c422994b6303af87d894424c3..fcc6e9c9cd8501447ff8d319f904432911616575 100755
--- a/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/hotRoom/Allclean
+++ b/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/hotRoom/Allclean
@@ -5,6 +5,6 @@ cd ${0%/*} || exit 1    # run from this directory
 . $WM_PROJECT_DIR/bin/tools/CleanFunctions
 
 cleanCase
-cp 0/T.org 0/T
+rm -f 0/T
 
 # ----------------------------------------------------------------- end-of-file
diff --git a/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/hotRoom/Allrun b/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/hotRoom/Allrun
index cb0be1e6e06e5244e2bddc492b8cb56b1ca04acd..800bb5731994d0f3a70b9f5ff7865f9659c1e4f7 100755
--- a/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/hotRoom/Allrun
+++ b/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/hotRoom/Allrun
@@ -8,6 +8,7 @@ application=`getApplication`
 
 compileApplication ../../buoyantPimpleFoam/hotRoom/setHotRoom
 runApplication blockMesh
+cp 0/T.org 0/T
 runApplication setHotRoom
 runApplication $application
 
diff --git a/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/hotRoom/system/controlDict b/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/hotRoom/system/controlDict
index dc6fb2083038ddd652546158951f808baa343ac6..820670013786a49ad59768fc77430307f6b2c088 100644
--- a/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/hotRoom/system/controlDict
+++ b/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/hotRoom/system/controlDict
@@ -45,5 +45,25 @@ timePrecision   6;
 
 runTimeModifiable true;
 
+functions
+{
+    residualControl1
+    {
+        type            residualControl;
+        functionObjectLibs ( "libjobControl.so" );
+        outputControl   timeStep;
+        outputInterval  1;
+
+        maxResiduals
+        {
+            p_rgh   1e-2;
+            U       1e-4;
+            T       1e-3;
+
+            // possibly check turbulence fields
+            "(k|epsilon|omega)" 1e-3;
+        }
+    }
+}
 
 // ************************************************************************* //
diff --git a/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/hotRoom/system/fvSolution b/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/hotRoom/system/fvSolution
index 68da3948e5ca1d13e30c97424ef5ef1a1eb8c1cf..c51949b80579958618776e4c7ce21eae0171eb65 100644
--- a/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/hotRoom/system/fvSolution
+++ b/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/hotRoom/system/fvSolution
@@ -37,17 +37,15 @@ solvers
 SIMPLE
 {
     nNonOrthogonalCorrectors 0;
-    p_rghRefCell    0;
-    p_rghRefValue   0;
+    pRefCell        0;
     pRefValue       0;
 }
 
 relaxationFactors
 {
-    rho             1;
     p_rgh           0.7;
     U               0.2;
-    T               0.7;
+    T               0.5;
     "(k|epsilon|R)" 0.7;
 }
 
diff --git a/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/iglooWithFridges/0/p b/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/iglooWithFridges/0/p
index 5486bbf9fcedf55174c7dc3c9f25ed67ccdc7e4b..3f7a67aa5e3d7442ddeeb510a4298b05deeea011 100644
--- a/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/iglooWithFridges/0/p
+++ b/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/iglooWithFridges/0/p
@@ -22,30 +22,26 @@ boundaryField
 {
     ground
     {
-        type            buoyantPressure;
-        rho             rhok;
-        value           uniform 0;
+        type            calculated;
+        value           $internalField;
     }
 
     igloo_region0
     {
-        type            buoyantPressure;
-        rho             rhok;
-        value           uniform 0;
+        type            calculated;
+        value           $internalField;
     }
 
     twoFridgeFreezers_seal_0
     {
-        type            buoyantPressure;
-        rho             rhok;
-        value           uniform 0;
+        type            calculated;
+        value           $internalField;
     }
 
     twoFridgeFreezers_herring_1
     {
-        type            buoyantPressure;
-        rho             rhok;
-        value           uniform 0;
+        type            calculated;
+        value           $internalField;
     }
 }
 
diff --git a/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/iglooWithFridges/constant/polyMesh/boundary b/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/iglooWithFridges/constant/polyMesh/boundary
index e3e3c88c206fc91de134a934be1f72721a804725..2876f3adad56b802155c54691fe41731d7fe90d3 100644
--- a/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/iglooWithFridges/constant/polyMesh/boundary
+++ b/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/iglooWithFridges/constant/polyMesh/boundary
@@ -21,55 +21,55 @@ FoamFile
     {
         type            empty;
         nFaces          0;
-        startFace       60456;
+        startFace       60336;
     }
     minX
     {
         type            empty;
         nFaces          0;
-        startFace       60456;
+        startFace       60336;
     }
     maxX
     {
         type            empty;
         nFaces          0;
-        startFace       60456;
+        startFace       60336;
     }
     minY
     {
         type            empty;
         nFaces          0;
-        startFace       60456;
+        startFace       60336;
     }
     ground
     {
         type            wall;
         nFaces          590;
-        startFace       60456;
+        startFace       60336;
     }
     maxZ
     {
         type            empty;
         nFaces          0;
-        startFace       61046;
+        startFace       60926;
     }
     igloo_region0
     {
         type            wall;
         nFaces          2260;
-        startFace       61046;
+        startFace       60926;
     }
     twoFridgeFreezers_seal_0
     {
         type            wall;
         nFaces          1344;
-        startFace       63306;
+        startFace       63186;
     }
     twoFridgeFreezers_herring_1
     {
         type            wall;
         nFaces          1116;
-        startFace       64650;
+        startFace       64530;
     }
 )
 
diff --git a/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/iglooWithFridges/system/controlDict b/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/iglooWithFridges/system/controlDict
index 73cfe8ada43473f1bb82beb55a0584a57aae0025..bf11b4f154badfe3e89c5fa7ff6f674294b40800 100644
--- a/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/iglooWithFridges/system/controlDict
+++ b/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/iglooWithFridges/system/controlDict
@@ -45,5 +45,25 @@ timePrecision   6;
 
 runTimeModifiable true;
 
+functions
+{
+    residualControl1
+    {
+        type            residualControl;
+        functionObjectLibs ( "libjobControl.so" );
+        outputControl   timeStep;
+        outputInterval  1;
+
+        maxResiduals
+        {
+            p_rgh   1e-2;
+            U       1e-4;
+            T       1e-3;
+
+            // possibly check turbulence fields
+            "(k|epsilon|omega)" 1e-3;
+        }
+    }
+}
 
 // ************************************************************************* //
diff --git a/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/iglooWithFridges/system/fvSolution b/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/iglooWithFridges/system/fvSolution
index 8df126a48366dc10aea3e57600f533461b45a65d..edea946f9ac02addf005716ef56d84c750cf2b06 100644
--- a/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/iglooWithFridges/system/fvSolution
+++ b/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/iglooWithFridges/system/fvSolution
@@ -37,16 +37,15 @@ solvers
 SIMPLE
 {
     nNonOrthogonalCorrectors 0;
-    p_rghRefCell    0;
-    p_rghRefValue   0;
+    pRefCell        0;
     pRefValue       0;
 }
 
 relaxationFactors
 {
-    p_rgh           0.8;
+    p_rgh           0.7;
     U               0.2;
-    T               0.7;
+    T               0.5;
     "(k|epsilon)"   0.7;
 }
 
diff --git a/tutorials/heatTransfer/buoyantPimpleFoam/hotRoom/0/p b/tutorials/heatTransfer/buoyantPimpleFoam/hotRoom/0/p
index 9854db7add58f4a75054fe9cff734f99d9304e8c..d6ab24211f697f7dfee98b8bfc5f8359c5bed1ef 100644
--- a/tutorials/heatTransfer/buoyantPimpleFoam/hotRoom/0/p
+++ b/tutorials/heatTransfer/buoyantPimpleFoam/hotRoom/0/p
@@ -22,20 +22,20 @@ boundaryField
 {
     floor
     {
-        type            buoyantPressure;
-        value           uniform 1e5;
+        type            calculated;
+        value           $internalField;
     }
 
     ceiling
     {
-        type            buoyantPressure;
-        value           uniform 1e5;
+        type            calculated;
+        value           $internalField;
     }
 
     fixedWalls
     {
-        type            buoyantPressure;
-        value           uniform 1e5;
+        type            calculated;
+        value           $internalField;
     }
 }
 
diff --git a/tutorials/heatTransfer/buoyantPimpleFoam/hotRoom/0/p_rgh b/tutorials/heatTransfer/buoyantPimpleFoam/hotRoom/0/p_rgh
new file mode 100644
index 0000000000000000000000000000000000000000..371d93c38aace14c5187b795c88f08d8f9f36bfa
--- /dev/null
+++ b/tutorials/heatTransfer/buoyantPimpleFoam/hotRoom/0/p_rgh
@@ -0,0 +1,42 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      p_rgh;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -2 0 0 0 0];
+
+internalField   uniform 1e5;
+
+boundaryField
+{
+    floor
+    {
+        type            buoyantPressure;
+        value           uniform 1e5;
+    }
+
+    ceiling
+    {
+        type            buoyantPressure;
+        value           uniform 1e5;
+    }
+
+    fixedWalls
+    {
+        type            buoyantPressure;
+        value           uniform 1e5;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/heatTransfer/buoyantPimpleFoam/hotRoom/system/fvSchemes b/tutorials/heatTransfer/buoyantPimpleFoam/hotRoom/system/fvSchemes
index 8caf680281b70920a18d0744a606c78262c44e3d..b238c08e08add804d21bb294586ad8a118c9cd8d 100644
--- a/tutorials/heatTransfer/buoyantPimpleFoam/hotRoom/system/fvSchemes
+++ b/tutorials/heatTransfer/buoyantPimpleFoam/hotRoom/system/fvSchemes
@@ -42,7 +42,7 @@ laplacianSchemes
 {
     default         none;
     laplacian(muEff,U) Gauss linear corrected;
-    laplacian((rho*(1|A(U))),p) Gauss linear corrected;
+    laplacian((rho*(1|A(U))),p_rgh) Gauss linear corrected;
     laplacian(alphaEff,h) Gauss linear corrected;
     laplacian(DkEff,k) Gauss linear corrected;
     laplacian(DepsilonEff,epsilon) Gauss linear corrected;
@@ -62,7 +62,7 @@ snGradSchemes
 fluxRequired
 {
     default         no;
-    p               ;
+    p_rgh;
 }
 
 
diff --git a/tutorials/heatTransfer/buoyantPimpleFoam/hotRoom/system/fvSolution b/tutorials/heatTransfer/buoyantPimpleFoam/hotRoom/system/fvSolution
index 9ad9aaa4ceaa48d329158cb6eeeb993a71e22214..0473bd2bf3b60c51713f03f3f70c26156b6933d3 100644
--- a/tutorials/heatTransfer/buoyantPimpleFoam/hotRoom/system/fvSolution
+++ b/tutorials/heatTransfer/buoyantPimpleFoam/hotRoom/system/fvSolution
@@ -25,17 +25,17 @@ solvers
         relTol          0;
     }
 
-    p
+    p_rgh
     {
         solver          PCG;
         preconditioner  DIC;
         tolerance       1e-8;
-        relTol          0.1;
+        relTol          0.01;
     }
 
-    pFinal
+    p_rghFinal
     {
-        $p;
+        $p_rgh;
         relTol          0;
     }
 
@@ -44,7 +44,7 @@ solvers
         solver          PBiCG;
         preconditioner  DILU;
         tolerance       1e-6;
-        relTol          0;
+        relTol          0.1;
     }
 
     "(U|h|k|epsilon|R)Final"
@@ -56,7 +56,7 @@ solvers
 
 PIMPLE
 {
-    momentumPredictor no;
+    momentumPredictor yes;
     nOuterCorrectors 1;
     nCorrectors     2;
     nNonOrthogonalCorrectors 0;
diff --git a/tutorials/heatTransfer/buoyantSimpleFoam/buoyantCavity/0/p b/tutorials/heatTransfer/buoyantSimpleFoam/buoyantCavity/0/p
index 0961828a5b013e1a4ef3778d971a1b10656fcc3f..9d710bdf49ee77f16275c99e13aa7727ecb93a51 100644
--- a/tutorials/heatTransfer/buoyantSimpleFoam/buoyantCavity/0/p
+++ b/tutorials/heatTransfer/buoyantSimpleFoam/buoyantCavity/0/p
@@ -23,26 +23,26 @@ boundaryField
 {
     frontAndBack 
     {
-        type            buoyantPressure;
-        value           uniform 1e5;
+        type            calculated;
+        value           $internalField;
     }
 
     topAndBottom
     {
-        type            buoyantPressure;
-        value           uniform 1e5;
+        type            calculated;
+        value           $internalField;
     }
 
     hot
     {
-        type            buoyantPressure;
-        value           uniform 1e5;
+        type            calculated;
+        value           $internalField;
     }
 
     cold
     {
-        type            buoyantPressure;
-        value           uniform 1e5;
+        type            calculated;
+        value           $internalField;
     }
 }
 
diff --git a/tutorials/heatTransfer/buoyantSimpleFoam/buoyantCavity/system/controlDict b/tutorials/heatTransfer/buoyantSimpleFoam/buoyantCavity/system/controlDict
index 8ca46687fba9ebfce54d100a24d5b415d7dd4e91..599204b449e686d742d5b57fd479705514e8930d 100644
--- a/tutorials/heatTransfer/buoyantSimpleFoam/buoyantCavity/system/controlDict
+++ b/tutorials/heatTransfer/buoyantSimpleFoam/buoyantCavity/system/controlDict
@@ -44,5 +44,25 @@ timePrecision   6;
 
 runTimeModifiable true;
 
+functions
+{
+    residualControl1
+    {
+        type            residualControl;
+        functionObjectLibs ( "libjobControl.so" );
+        outputControl   timeStep;
+        outputInterval  1;
+
+        maxResiduals
+        {
+            p_rgh   1e-2;
+            U       1e-4;
+            T       1e-3;
+
+            // possibly check turbulence fields
+            "(k|epsilon|omega)" 1e-3;
+        }
+    }
+}
 
 // ************************************************************************* //
diff --git a/tutorials/heatTransfer/buoyantSimpleFoam/buoyantCavity/system/fvSolution b/tutorials/heatTransfer/buoyantSimpleFoam/buoyantCavity/system/fvSolution
index 992508c411d48ac26bc83378a684c3205bb2cfb4..c453b27d81419a3accdc50084dcf14be08be2cbb 100644
--- a/tutorials/heatTransfer/buoyantSimpleFoam/buoyantCavity/system/fvSolution
+++ b/tutorials/heatTransfer/buoyantSimpleFoam/buoyantCavity/system/fvSolution
@@ -44,17 +44,16 @@ SIMPLE
 {
     momentumPredictor yes;
     nNonOrthogonalCorrectors 0;
-    p_rghRefCell    0;
-    p_rghRefValue   100000;
-    pRefValue       100000;
-    convergence     1e-04;
+    pRefCell        0;
+    pRefValue       0;
 }
 
 relaxationFactors
 {
-    p_rgh           0.9;
+    rho             1.0;
+    p_rgh           0.7;
     U               0.3;
-    h               0.7;
+    h               0.3;
     "(k|epsilon|omega)" 0.7;
 }
 
diff --git a/tutorials/heatTransfer/buoyantSimpleFoam/hotRoom/0/T b/tutorials/heatTransfer/buoyantSimpleFoam/hotRoom/0/T
index 98cc191af7c9e6234287f45ebeab0cdea6ab8b62..77b58a1a00e6421c4d9c3f72a36b3d190e41ef7f 100644
--- a/tutorials/heatTransfer/buoyantSimpleFoam/hotRoom/0/T
+++ b/tutorials/heatTransfer/buoyantSimpleFoam/hotRoom/0/T
@@ -23,411 +23,7 @@ boundaryField
     floor
     {
         type            fixedValue;
-        value           nonuniform List<scalar> 
-400
-(
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-600
-600
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-600
-600
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-)
-;
+        value           uniform 300; 
     }
     ceiling
     {
diff --git a/tutorials/heatTransfer/buoyantSimpleFoam/hotRoom/0/T.org b/tutorials/heatTransfer/buoyantSimpleFoam/hotRoom/0/T.org
index 98cc191af7c9e6234287f45ebeab0cdea6ab8b62..77b58a1a00e6421c4d9c3f72a36b3d190e41ef7f 100644
--- a/tutorials/heatTransfer/buoyantSimpleFoam/hotRoom/0/T.org
+++ b/tutorials/heatTransfer/buoyantSimpleFoam/hotRoom/0/T.org
@@ -23,411 +23,7 @@ boundaryField
     floor
     {
         type            fixedValue;
-        value           nonuniform List<scalar> 
-400
-(
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-600
-600
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-600
-600
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-300
-)
-;
+        value           uniform 300; 
     }
     ceiling
     {
diff --git a/tutorials/heatTransfer/buoyantSimpleFoam/hotRoom/0/p b/tutorials/heatTransfer/buoyantSimpleFoam/hotRoom/0/p
index 9854db7add58f4a75054fe9cff734f99d9304e8c..d6ab24211f697f7dfee98b8bfc5f8359c5bed1ef 100644
--- a/tutorials/heatTransfer/buoyantSimpleFoam/hotRoom/0/p
+++ b/tutorials/heatTransfer/buoyantSimpleFoam/hotRoom/0/p
@@ -22,20 +22,20 @@ boundaryField
 {
     floor
     {
-        type            buoyantPressure;
-        value           uniform 1e5;
+        type            calculated;
+        value           $internalField;
     }
 
     ceiling
     {
-        type            buoyantPressure;
-        value           uniform 1e5;
+        type            calculated;
+        value           $internalField;
     }
 
     fixedWalls
     {
-        type            buoyantPressure;
-        value           uniform 1e5;
+        type            calculated;
+        value           $internalField;
     }
 }
 
diff --git a/tutorials/heatTransfer/buoyantSimpleFoam/hotRoom/0/p_rgh b/tutorials/heatTransfer/buoyantSimpleFoam/hotRoom/0/p_rgh
new file mode 100644
index 0000000000000000000000000000000000000000..371d93c38aace14c5187b795c88f08d8f9f36bfa
--- /dev/null
+++ b/tutorials/heatTransfer/buoyantSimpleFoam/hotRoom/0/p_rgh
@@ -0,0 +1,42 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      p_rgh;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -2 0 0 0 0];
+
+internalField   uniform 1e5;
+
+boundaryField
+{
+    floor
+    {
+        type            buoyantPressure;
+        value           uniform 1e5;
+    }
+
+    ceiling
+    {
+        type            buoyantPressure;
+        value           uniform 1e5;
+    }
+
+    fixedWalls
+    {
+        type            buoyantPressure;
+        value           uniform 1e5;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/heatTransfer/buoyantSimpleFoam/hotRoom/system/controlDict b/tutorials/heatTransfer/buoyantSimpleFoam/hotRoom/system/controlDict
index f9d959e2dc8b1c2493a4b0e826157c7d4d34e8df..242675a288a20edab069535d3aba9bc598a47f83 100644
--- a/tutorials/heatTransfer/buoyantSimpleFoam/hotRoom/system/controlDict
+++ b/tutorials/heatTransfer/buoyantSimpleFoam/hotRoom/system/controlDict
@@ -45,5 +45,25 @@ timePrecision   6;
 
 runTimeModifiable true;
 
+functions
+{
+    residualControl1
+    {
+        type            residualControl;
+        functionObjectLibs ( "libjobControl.so" );
+        outputControl   timeStep;
+        outputInterval  1;
+
+        maxResiduals
+        {
+            p_rgh   1e-2;
+            U       1e-4;
+            T       1e-3;
+
+            // possibly check turbulence fields
+            "(k|epsilon|omega)" 1e-3;
+        }
+    }
+}
 
 // ************************************************************************* //
diff --git a/tutorials/heatTransfer/buoyantSimpleFoam/hotRoom/system/fvSchemes b/tutorials/heatTransfer/buoyantSimpleFoam/hotRoom/system/fvSchemes
index 328bd3ed59418bf59f3a1fad1e6aebee33ef49ce..ab84aca2c4ab2ac7d880eefbb91b2ebf1b99d40a 100644
--- a/tutorials/heatTransfer/buoyantSimpleFoam/hotRoom/system/fvSchemes
+++ b/tutorials/heatTransfer/buoyantSimpleFoam/hotRoom/system/fvSchemes
@@ -40,12 +40,12 @@ divSchemes
 laplacianSchemes
 {
     default         none;
-    laplacian(muEff,U) Gauss linear corrected;
-    laplacian((rho*(1|A(U))),p) Gauss linear corrected;
-    laplacian(alphaEff,h) Gauss linear corrected;
-    laplacian(DkEff,k) Gauss linear corrected;
-    laplacian(DepsilonEff,epsilon) Gauss linear corrected;
-    laplacian(DREff,R) Gauss linear corrected;
+    laplacian(muEff,U) Gauss linear uncorrected;
+    laplacian((rho*(1|A(U))),p_rgh) Gauss linear uncorrected;
+    laplacian(alphaEff,h) Gauss linear uncorrected;
+    laplacian(DkEff,k) Gauss linear uncorrected;
+    laplacian(DepsilonEff,epsilon) Gauss linear uncorrected;
+    laplacian(DREff,R) Gauss linear uncorrected;
 }
 
 interpolationSchemes
@@ -61,7 +61,7 @@ snGradSchemes
 fluxRequired
 {
     default         no;
-    p               ;
+    p_rgh;
 }
 
 
diff --git a/tutorials/heatTransfer/buoyantSimpleFoam/hotRoom/system/fvSolution b/tutorials/heatTransfer/buoyantSimpleFoam/hotRoom/system/fvSolution
index c9a423a0631042e51b594d56419b28bb1d66bd85..aa535a28eb93c623863c7191fa7c33ccbb26fe2e 100644
--- a/tutorials/heatTransfer/buoyantSimpleFoam/hotRoom/system/fvSolution
+++ b/tutorials/heatTransfer/buoyantSimpleFoam/hotRoom/system/fvSolution
@@ -17,12 +17,12 @@ FoamFile
 
 solvers
 {
-    p
+    p_rgh
     {
         solver          PCG;
         preconditioner  DIC;
         tolerance       1e-08;
-        relTol          0;
+        relTol          0.01;
     }
 
     "(U|h|k|epsilon|R)"
@@ -30,7 +30,7 @@ solvers
         solver          PBiCG;
         preconditioner  DILU;
         tolerance       1e-05;
-        relTol          0;
+        relTol          0.1;
     }
 }
 
@@ -38,16 +38,16 @@ SIMPLE
 {
     nNonOrthogonalCorrectors 0;
     pRefCell        0;
-    pRefValue       100000;
+    pRefValue       0;
 }
 
 relaxationFactors
 {
-    rho             1;
-    p               0.7;
+    rho             1.0;
+    p_rgh           0.7;
     U               0.2;
-    h               0.7;
-    "(k|epsilon|R)" 0.7;
+    h               0.2;
+    "(k|epsilon|R)" 0.5;
 }
 
 
diff --git a/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoom/0/G b/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoom/0/G
index d49a1b712b778f241fd43b9bf23c315b3d772fc3..928afaeb56586af05a00b38f3c6c33ffbec6cf1b 100644
--- a/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoom/0/G
+++ b/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoom/0/G
@@ -25,7 +25,7 @@ boundaryField
         type            MarshakRadiation;
         T               T;
         emissivity      1;
-        value           uniform 0;
+//        value           uniform 0;
     }
 
     fixedWalls
@@ -33,7 +33,7 @@ boundaryField
         type            MarshakRadiation;
         T               T;
         emissivity      1;
-        value           uniform 0;
+//        value           uniform 0;
     }
 
     ceiling
@@ -41,7 +41,7 @@ boundaryField
         type            MarshakRadiation;
         T               T;
         emissivity      1;
-        value           uniform 0;
+ //       value           uniform 0;
     }
 
     box
@@ -49,7 +49,7 @@ boundaryField
         type            MarshakRadiation;
         T               T;
         emissivity      1;
-        value           uniform 0;
+//        value           uniform 0;
     }
 }
 
diff --git a/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoom/0/p b/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoom/0/p
index a1da943200cf2bd68718b4eb8f1252e550968663..fb9be497bfd16e3bb979bd490191158bbc635d74 100644
--- a/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoom/0/p
+++ b/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoom/0/p
@@ -22,26 +22,26 @@ boundaryField
 {
     floor
     {
-        type            buoyantPressure;
-        value           uniform 100000;
+        type            calculated;
+        value           $internalField;
     }
 
     ceiling
     {
-        type            buoyantPressure;
-        value           uniform 100000;
+        type            calculated;
+        value           $internalField;
     }
 
     fixedWalls
     {
-        type            buoyantPressure;
-        value           uniform 100000;
+        type            calculated;
+        value           $internalField;
     }
 
     box
     {
-        type            buoyantPressure;
-        value           uniform 100000;
+        type            calculated;
+        value           $internalField;
     }
 }
 
diff --git a/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoom/0/p_rgh b/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoom/0/p_rgh
new file mode 100644
index 0000000000000000000000000000000000000000..87dd59413431afaeed37d495039e21fcbddaf000
--- /dev/null
+++ b/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoom/0/p_rgh
@@ -0,0 +1,48 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      p_rgh;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -2 0 0 0 0];
+
+internalField   uniform 100000;
+
+boundaryField
+{
+    floor
+    {
+        type            buoyantPressure;
+        value           uniform 100000;
+    }
+
+    ceiling
+    {
+        type            buoyantPressure;
+        value           uniform 100000;
+    }
+
+    fixedWalls
+    {
+        type            buoyantPressure;
+        value           uniform 100000;
+    }
+
+    box
+    {
+        type            buoyantPressure;
+        value           uniform 100000;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoom/system/controlDict b/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoom/system/controlDict
index 25a67762b362324fe7617c9716cc0fb200611ba4..80041f81dc9b7cc6e5ca7e51bc58c76f3534103e 100644
--- a/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoom/system/controlDict
+++ b/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoom/system/controlDict
@@ -45,5 +45,25 @@ timePrecision   6;
 
 runTimeModifiable true;
 
+functions
+{
+    residualControl1
+    {
+        type            residualControl;
+        functionObjectLibs ( "libjobControl.so" );
+        outputControl   timeStep;
+        outputInterval  1;
+
+        maxResiduals
+        {
+            p_rgh   1e-2;
+            U       1e-4;
+            T       1e-3;
+
+            // possibly check turbulence fields
+            "(k|epsilon|omega)" 1e-3;
+        }
+    }
+}
 
 // ************************************************************************* //
diff --git a/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoom/system/fvSchemes b/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoom/system/fvSchemes
index 8d30fa0a3790c1431883202074d0c266a8b74b36..2dfb6798b87c78561c73b8f93271895718f7403b 100644
--- a/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoom/system/fvSchemes
+++ b/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoom/system/fvSchemes
@@ -41,7 +41,7 @@ laplacianSchemes
 {
     default         none;
     laplacian(muEff,U) Gauss linear corrected;
-    laplacian((rho*(1|A(U))),p) Gauss linear corrected;
+    laplacian((rho*(1|A(U))),p_rgh) Gauss linear corrected;
     laplacian(alphaEff,h) Gauss linear corrected;
     laplacian(DkEff,k) Gauss linear corrected;
     laplacian(DepsilonEff,epsilon) Gauss linear corrected;
@@ -62,7 +62,7 @@ snGradSchemes
 fluxRequired
 {
     default         no;
-    p;
+    p_rgh;
 }
 
 
diff --git a/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoom/system/fvSolution b/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoom/system/fvSolution
index c490be6e8f4f101d157ad2dd9dd71219ef9137b0..fa8d3d533976867e6584d6a9e8df8f0b12d8c776 100644
--- a/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoom/system/fvSolution
+++ b/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoom/system/fvSolution
@@ -17,7 +17,7 @@ FoamFile
 
 solvers
 {
-    p
+    p_rgh
     {
         solver          PCG;
         preconditioner  DIC;
@@ -35,7 +35,7 @@ solvers
 
     G
     {
-        $p;
+        $p_rgh;
         tolerance       1e-05;
         relTol          0.1;
     }
@@ -50,11 +50,11 @@ SIMPLE
 
 relaxationFactors
 {
-    rho             1;
-    p               0.3;
-    U               0.7;
-    h               0.7;
-    "(k|epsilon)"   0.7;
+    rho             1.0;
+    p_rgh           0.7;
+    U               0.2;
+    h               0.2;
+    "(k|epsilon|R)" 0.5;
     G               0.7;
 }
 
diff --git a/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoomFvDOM/0/IDefault b/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoomFvDOM/0/IDefault
index 3bee6a9b107a580ac654fc8574bd487d17655a9a..94f365296e9262ca775182d6bba125ae7bd39c16 100644
--- a/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoomFvDOM/0/IDefault
+++ b/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoomFvDOM/0/IDefault
@@ -24,7 +24,7 @@ boundaryField
     {
         type            greyDiffusiveRadiation;
         T               T;
-        emissivity      0.5;
+        emissivity      1.0;
         value           uniform 0;
     }
 }
diff --git a/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoomFvDOM/0/p b/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoomFvDOM/0/p
index a1da943200cf2bd68718b4eb8f1252e550968663..48fe8a009fd10c7b86b76f7de6160ccb31d5045c 100644
--- a/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoomFvDOM/0/p
+++ b/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoomFvDOM/0/p
@@ -22,25 +22,25 @@ boundaryField
 {
     floor
     {
-        type            buoyantPressure;
+        type            calculated;
         value           uniform 100000;
     }
 
     ceiling
     {
-        type            buoyantPressure;
+        type            calculated;
         value           uniform 100000;
     }
 
     fixedWalls
     {
-        type            buoyantPressure;
+        type            calculated;
         value           uniform 100000;
     }
 
     box
     {
-        type            buoyantPressure;
+        type            calculated;
         value           uniform 100000;
     }
 }
diff --git a/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoomFvDOM/0/p_rgh b/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoomFvDOM/0/p_rgh
new file mode 100644
index 0000000000000000000000000000000000000000..87dd59413431afaeed37d495039e21fcbddaf000
--- /dev/null
+++ b/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoomFvDOM/0/p_rgh
@@ -0,0 +1,48 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      p_rgh;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -2 0 0 0 0];
+
+internalField   uniform 100000;
+
+boundaryField
+{
+    floor
+    {
+        type            buoyantPressure;
+        value           uniform 100000;
+    }
+
+    ceiling
+    {
+        type            buoyantPressure;
+        value           uniform 100000;
+    }
+
+    fixedWalls
+    {
+        type            buoyantPressure;
+        value           uniform 100000;
+    }
+
+    box
+    {
+        type            buoyantPressure;
+        value           uniform 100000;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoomFvDOM/system/controlDict b/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoomFvDOM/system/controlDict
index 1cca62ded7f545918bfa54e88eb08c241fbeaba5..64496570c2a733a336e0b21789310937b7b0e857 100644
--- a/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoomFvDOM/system/controlDict
+++ b/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoomFvDOM/system/controlDict
@@ -45,5 +45,25 @@ timePrecision   6;
 
 runTimeModifiable true;
 
+functions
+{
+    residualControl1
+    {
+        type            residualControl;
+        functionObjectLibs ( "libjobControl.so" );
+        outputControl   timeStep;
+        outputInterval  1;
+
+        maxResiduals
+        {
+            p_rgh   1e-2;
+            U       1e-4;
+            T       1e-3;
+
+            // possibly check turbulence fields
+            "(k|epsilon|omega)" 1e-3;
+        }
+    }
+}
 
 // ************************************************************************* //
diff --git a/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoomFvDOM/system/fvSchemes b/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoomFvDOM/system/fvSchemes
index 58fe62417e93f0e2cc0fa7c802c471f3aa200e66..c9900e695fcb668a7fd7f64cc89c1d7c3080aedf 100644
--- a/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoomFvDOM/system/fvSchemes
+++ b/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoomFvDOM/system/fvSchemes
@@ -42,7 +42,7 @@ laplacianSchemes
 {
     default         none;
     laplacian(muEff,U) Gauss linear corrected;
-    laplacian((rho*(1|A(U))),p) Gauss linear corrected;
+    laplacian((rho*(1|A(U))),p_rgh) Gauss linear corrected;
     laplacian(alphaEff,h) Gauss linear corrected;
     laplacian(DkEff,k) Gauss linear corrected;
     laplacian(DepsilonEff,epsilon) Gauss linear corrected;
@@ -63,7 +63,7 @@ snGradSchemes
 fluxRequired
 {
     default         no;
-    p;
+    p_rgh;
 }
 
 
diff --git a/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoomFvDOM/system/fvSolution b/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoomFvDOM/system/fvSolution
index c219a76e03fea5e0e6ed24f300a4da969c3c71fc..daa2528426bab6d94d3859b6323754407f9fb261 100644
--- a/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoomFvDOM/system/fvSolution
+++ b/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoomFvDOM/system/fvSolution
@@ -17,7 +17,7 @@ FoamFile
 
 solvers
 {
-    p
+    p_rgh
     {
         solver          GAMG;
         tolerance       1e-06;
@@ -47,11 +47,12 @@ SIMPLE
 
 relaxationFactors
 {
-    rho             1;
-    p               0.3;
+    rho             1.0;
+    p_rgh           0.3;
     U               0.7;
     h               0.7;
-    "(k|epsilon)"   0.7;
+    k               0.7;
+    epsilon         0.7;
     "ILambda.*"     0.7;
 }
 
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/0/epsilon b/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/0/epsilon
index 6034a1d7bfc3f328fece0218588466e1deb72932..f4e27ad95fb84fff1121e59c121f6d2c844d1b43 100644
--- a/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/0/epsilon
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/0/epsilon
@@ -10,7 +10,6 @@ FoamFile
     version     2.0;
     format      ascii;
     class       volScalarField;
-    location    "0.001";
     object      epsilon;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/0/k b/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/0/k
index eb1abac09b21d07f72329ef894699f330741de42..b2509f967076c366a4ab4e410f070916983536e4 100644
--- a/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/0/k
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/0/k
@@ -10,7 +10,6 @@ FoamFile
     version     2.0;
     format      ascii;
     class       volScalarField;
-    location    "0.001";
     object      k;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/0/p_rgh b/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/0/p_rgh
new file mode 100644
index 0000000000000000000000000000000000000000..638c8c7d7cd472cf343b4fc5091d9a8f1bf038ae
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/0/p_rgh
@@ -0,0 +1,29 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      p_rgh;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -2 0 0 0 0];
+
+internalField   uniform 1e5;
+
+boundaryField
+{
+    ".*"
+    {
+        type            calculated;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/bottomAir/changeDictionaryDict b/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/bottomAir/changeDictionaryDict
index c36ea6b8039bbf390c437d1875aa6eebe84cf63f..6faab27ee56705411dc2abcf18838952a6f25fa1 100644
--- a/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/bottomAir/changeDictionaryDict
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/bottomAir/changeDictionaryDict
@@ -92,16 +92,30 @@ dictionaryReplacement
         }
     }
 
-    p
+    p_rgh
     {
-        internalField   uniform 100000;
+        internalField   uniform 1e5;
 
         boundaryField
         {
             ".*"
             {
                 type            buoyantPressure;
-                value           1e5;
+                value           uniform 1e5;
+            }
+        }
+    }
+
+    p
+    {
+        internalField   uniform 1e5;
+
+        boundaryField
+        {
+            ".*"
+            {
+                type            calculated;
+                value           uniform 1e5;
             }
         }
     }
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/bottomAir/fvSchemes b/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/bottomAir/fvSchemes
index 90364fcbf2309c83fb6dbf1002215f43b97d23d5..7d3e57864ac0660fb7bd646e2850a67eafb1d49d 100644
--- a/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/bottomAir/fvSchemes
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/bottomAir/fvSchemes
@@ -17,13 +17,11 @@ FoamFile
 ddtSchemes
 {
     default Euler;
-    //default CoEuler phi rho 0.1;
 }
 
 gradSchemes
 {
     default         Gauss linear;
-//    grad(U)         cellLimited Gauss linear 1;
 }
 
 divSchemes
@@ -43,7 +41,7 @@ laplacianSchemes
 {
     default         none;
     laplacian(muEff,U) Gauss linear limited 0.333;
-    laplacian((rho*(1|A(U))),p) Gauss linear limited 0.333;
+    laplacian((rho*(1|A(U))),p_rgh) Gauss linear limited 0.333;
     laplacian(alphaEff,h) Gauss linear limited 0.333;
     laplacian(DkEff,k) Gauss linear limited 0.333;
     laplacian(DepsilonEff,epsilon) Gauss linear limited 0.333;
@@ -63,7 +61,7 @@ snGradSchemes
 fluxRequired
 {
     default         no;
-    p;
+    p_rgh;
 }
 
 // ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/bottomAir/fvSolution b/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/bottomAir/fvSolution
index 16665b05a12a2e01e1211ecd4534473e4d42960b..1e968f0aea32e538282ed8a628684928be7acdfc 100644
--- a/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/bottomAir/fvSolution
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/bottomAir/fvSolution
@@ -24,7 +24,7 @@ solvers
         relTol          0;
     }
 
-    p
+    p_rgh
     {
         solver           GAMG;
         tolerance        1e-8;
@@ -38,31 +38,26 @@ solvers
         mergeLevels      1;
     }
 
-    pFinal
+    p_rghFinal
     {
-        $p;
+        $p_rgh;
         tolerance        1e-8;
         relTol           0;
     }
 
-    U
+    "(U|h|k|epsilon|R)"
     {
         solver           PBiCG;
         preconditioner   DILU;
         tolerance        1e-08;
-        relTol           0;
-    }
-
-    h
-    {
-        $U;
-        tolerance        1e-08;
         relTol           0.1;
     }
 
-    "(hFinal|k|epsilon|R)"
+    "(U|h|k|epsilon|R)Final"
     {
         $U;
+        tolerance        1e-08;
+        relTol           0;
     }
 }
 
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/controlDict b/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/controlDict
index ec9f3455bd051faa1da52f6088d977381ec873b1..590a47c55a6124c258398c7e6ce82f659aed0aea 100644
--- a/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/controlDict
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/controlDict
@@ -35,9 +35,9 @@ writeInterval   50;
 
 purgeWrite      0;
 
-writeFormat     ascii;
+writeFormat     binary;
 
-writePrecision  7;
+writePrecision  8;
 
 writeCompression off;
 
@@ -45,11 +45,12 @@ timeFormat      general;
 
 timePrecision   6;
 
-runTimeModifiable true;
+runTimeModifiable yes;
 
 maxCo           0.3;
 
-maxDi          10.0;
+// Maximum diffusion number
+maxDi           10.0;
 
 adjustTimeStep  yes;
 
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/heater/decomposeParDict b/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/heater/decomposeParDict
index e53a523c9d58e8314f6eeea8b310d284b54f6cd0..8f6ced4ac38fc5f058b8c4370ed4002c036d61f6 100644
--- a/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/heater/decomposeParDict
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/heater/decomposeParDict
@@ -76,7 +76,7 @@ manualCoeffs
 //// Is the case distributed
 //distributed     yes;
 //// Per slave (so nProcs-1 entries) the directory above the case.
-//roots           
+//roots
 //(
 //    "/tmp"
 //    "/tmp"
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/heater/fvSolution b/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/heater/fvSolution
index b773399f729bcb6ee8fc6b7c458c3d82e9ea2a8f..7c00e466f4b508dc216b21874718108befe65e16 100644
--- a/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/heater/fvSolution
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/heater/fvSolution
@@ -21,6 +21,12 @@ solvers
         solver           PCG;
         preconditioner   DIC;
         tolerance        1E-06;
+        relTol           0.1;
+    }
+    TFinal
+    {
+        $T;
+        tolerance        1E-06;
         relTol           0;
     }
 }
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/leftSolid/fvSolution b/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/leftSolid/fvSolution
index b773399f729bcb6ee8fc6b7c458c3d82e9ea2a8f..7c00e466f4b508dc216b21874718108befe65e16 100644
--- a/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/leftSolid/fvSolution
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/leftSolid/fvSolution
@@ -21,6 +21,12 @@ solvers
         solver           PCG;
         preconditioner   DIC;
         tolerance        1E-06;
+        relTol           0.1;
+    }
+    TFinal
+    {
+        $T;
+        tolerance        1E-06;
         relTol           0;
     }
 }
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/rightSolid/fvSolution b/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/rightSolid/fvSolution
index b773399f729bcb6ee8fc6b7c458c3d82e9ea2a8f..7c00e466f4b508dc216b21874718108befe65e16 100644
--- a/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/rightSolid/fvSolution
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/rightSolid/fvSolution
@@ -21,6 +21,12 @@ solvers
         solver           PCG;
         preconditioner   DIC;
         tolerance        1E-06;
+        relTol           0.1;
+    }
+    TFinal
+    {
+        $T;
+        tolerance        1E-06;
         relTol           0;
     }
 }
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/topAir/changeDictionaryDict b/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/topAir/changeDictionaryDict
index 71114003ec106a4a9af29c5abfa86f94cf185c56..ff4ba100119b83b40959e27389c224519da39a2e 100644
--- a/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/topAir/changeDictionaryDict
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/topAir/changeDictionaryDict
@@ -18,7 +18,7 @@ dictionaryReplacement
 {
     U
     {
-        internalField   uniform (0.01 0 0);
+        internalField   uniform (0.1 0 0);
 
         boundaryField
         {
@@ -30,13 +30,13 @@ dictionaryReplacement
             minX
             {
                 type            fixedValue;
-                value           uniform ( 0.01 0 0 );
+                value           uniform ( 0.1 0 0 );
             }
             maxX
             {
                 type            inletOutlet;
                 inletValue      uniform ( 0 0 0 );
-                value           uniform ( 0 0 0 );
+                value           uniform ( 0.1 0 0 );
             }
         }
     }
@@ -127,25 +127,42 @@ dictionaryReplacement
         }
     }
 
-    p
+    p_rgh
     {
-        internalField   uniform 100000;
+        internalField   uniform 1e5;
 
         boundaryField
         {
             ".*"
             {
                 type            buoyantPressure;
-                value           1e5;
+                value           uniform 1e5;
+            }
+
+            maxX
+            {
+                type            fixedValue;
+                value           uniform 1e5;
+            }
+        }
+    }
+
+    p
+    {
+        internalField   uniform 1e5;
+
+        boundaryField
+        {
+            ".*"
+            {
+                type            calculated;
+                value           uniform 1e5;
             }
 
             maxX
             {
-                type            waveTransmissive;
-                gamma           1.4;
-                fieldInf        100000;
-                lInf            0.4;
-                value           uniform 100000;
+                type            calculated;
+                value           uniform 1e5;
             }
         }
     }
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/topAir/fvSchemes b/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/topAir/fvSchemes
index 90364fcbf2309c83fb6dbf1002215f43b97d23d5..7d3e57864ac0660fb7bd646e2850a67eafb1d49d 100644
--- a/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/topAir/fvSchemes
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/topAir/fvSchemes
@@ -17,13 +17,11 @@ FoamFile
 ddtSchemes
 {
     default Euler;
-    //default CoEuler phi rho 0.1;
 }
 
 gradSchemes
 {
     default         Gauss linear;
-//    grad(U)         cellLimited Gauss linear 1;
 }
 
 divSchemes
@@ -43,7 +41,7 @@ laplacianSchemes
 {
     default         none;
     laplacian(muEff,U) Gauss linear limited 0.333;
-    laplacian((rho*(1|A(U))),p) Gauss linear limited 0.333;
+    laplacian((rho*(1|A(U))),p_rgh) Gauss linear limited 0.333;
     laplacian(alphaEff,h) Gauss linear limited 0.333;
     laplacian(DkEff,k) Gauss linear limited 0.333;
     laplacian(DepsilonEff,epsilon) Gauss linear limited 0.333;
@@ -63,7 +61,7 @@ snGradSchemes
 fluxRequired
 {
     default         no;
-    p;
+    p_rgh;
 }
 
 // ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/topAir/fvSolution b/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/topAir/fvSolution
index 28c1a63b7aa3986ac043fd9eea7b5975e68d7c5d..d2d44df4fc5f85c3a91d97b0f5402c4649e7102f 100644
--- a/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/topAir/fvSolution
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/topAir/fvSolution
@@ -20,14 +20,14 @@ solvers
     {
         solver          PCG
         preconditioner  DIC;
-        tolerance       1e-8;
+        tolerance       1e-7;
         relTol          0;
     }
 
-    p
+    p_rgh
     {
         solver           GAMG;
-        tolerance        1e-8;
+        tolerance        1e-7;
         relTol           0.01;
 
         smoother         GaussSeidel;
@@ -36,47 +36,43 @@ solvers
         nCellsInCoarsestLevel 10;
         agglomerator     faceAreaPair;
         mergeLevels      1;
+
+        maxIter          100;
     }
 
-    pFinal
+    p_rghFinal
     {
-        $p;
-        tolerance        1e-8;
+        $p_rgh;
+        tolerance        1e-7;
         relTol           0;
     }
 
-    U
+    "(U|h|k|epsilon|R)"
     {
         solver           PBiCG;
         preconditioner   DILU;
-        tolerance        1e-08;
-        relTol           0;
-    }
-
-    h
-    {
-        $U;
-        tolerance        1e-08;
+        tolerance        1e-7;
         relTol           0.1;
     }
 
-    "(hFinal|k|epsilon|R)"
+    "(U|h|k|epsilon|R)Final"
     {
         $U;
+        tolerance        1e-07;
+        relTol           0;
     }
 }
 
 PISO
 {
-    momentumPredictor   off;
+    momentumPredictor    off;
     nOuterCorrectors     1;
-    nCorrectors     2;
+    nCorrectors          2;
     nNonOrthogonalCorrectors 1;
-    pRefPoint       (-0.081 -0.0257 8.01);
-    pRefValue       1e5;
+    pRefPoint            (-0.081 -0.0257 8.01);
+    pRefValue            1e5;
 }
 
-
 PIMPLE
 {
     momentumPredictor   on;
@@ -84,11 +80,10 @@ PIMPLE
     nNonOrthogonalCorrectors 0;
 }
 
-
 relaxationFactors
 {
-//    h               0.9;
-//    U               0.9;
+    h               1;
+    U               1;
 }
 
 // ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/0/alphat b/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/0/alphat
index 5afdd3b3cc00a9e5b503e967ee5c0685f489464b..5dba622286ec20e5ed7cf843857f7587026841fc 100644
--- a/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/0/alphat
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/0/alphat
@@ -10,7 +10,6 @@ FoamFile
     version     2.0;
     format      ascii;
     class       volScalarField;
-    location    "0.001";
     object      alphat;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/0/epsilon b/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/0/epsilon
index 6034a1d7bfc3f328fece0218588466e1deb72932..f4e27ad95fb84fff1121e59c121f6d2c844d1b43 100644
--- a/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/0/epsilon
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/0/epsilon
@@ -10,7 +10,6 @@ FoamFile
     version     2.0;
     format      ascii;
     class       volScalarField;
-    location    "0.001";
     object      epsilon;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/0/k b/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/0/k
index eb1abac09b21d07f72329ef894699f330741de42..b2509f967076c366a4ab4e410f070916983536e4 100644
--- a/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/0/k
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/0/k
@@ -10,7 +10,6 @@ FoamFile
     version     2.0;
     format      ascii;
     class       volScalarField;
-    location    "0.001";
     object      k;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/0/p_rgh b/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/0/p_rgh
new file mode 100644
index 0000000000000000000000000000000000000000..638c8c7d7cd472cf343b4fc5091d9a8f1bf038ae
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/0/p_rgh
@@ -0,0 +1,29 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      p_rgh;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -2 0 0 0 0];
+
+internalField   uniform 1e5;
+
+boundaryField
+{
+    ".*"
+    {
+        type            calculated;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/bottomAir/changeDictionaryDict b/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/bottomAir/changeDictionaryDict
index 1e2a7dd206395a389f4621e82d7cd22813dbbaea..fc51ab1ad718189e9d837a7d71fd329a12f39a47 100644
--- a/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/bottomAir/changeDictionaryDict
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/bottomAir/changeDictionaryDict
@@ -94,16 +94,30 @@ dictionaryReplacement
         }
     }
 
-    p
+    p_rgh
     {
-        internalField   uniform 100000;
+        internalField   uniform 1e5;
 
         boundaryField
         {
             ".*"
             {
                 type            buoyantPressure;
-                value           1e5;
+                value           uniform 1e5;
+            }
+        }
+    }
+
+    p
+    {
+        internalField   uniform 1e5;
+
+        boundaryField
+        {
+            ".*"
+            {
+                type            calculated;
+                value           uniform 1e5;
             }
         }
     }
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/bottomAir/decomposeParDict b/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/bottomAir/decomposeParDict
index e53a523c9d58e8314f6eeea8b310d284b54f6cd0..8f6ced4ac38fc5f058b8c4370ed4002c036d61f6 100644
--- a/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/bottomAir/decomposeParDict
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/bottomAir/decomposeParDict
@@ -76,7 +76,7 @@ manualCoeffs
 //// Is the case distributed
 //distributed     yes;
 //// Per slave (so nProcs-1 entries) the directory above the case.
-//roots           
+//roots
 //(
 //    "/tmp"
 //    "/tmp"
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/bottomAir/fvSchemes b/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/bottomAir/fvSchemes
index 90364fcbf2309c83fb6dbf1002215f43b97d23d5..7d3e57864ac0660fb7bd646e2850a67eafb1d49d 100644
--- a/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/bottomAir/fvSchemes
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/bottomAir/fvSchemes
@@ -17,13 +17,11 @@ FoamFile
 ddtSchemes
 {
     default Euler;
-    //default CoEuler phi rho 0.1;
 }
 
 gradSchemes
 {
     default         Gauss linear;
-//    grad(U)         cellLimited Gauss linear 1;
 }
 
 divSchemes
@@ -43,7 +41,7 @@ laplacianSchemes
 {
     default         none;
     laplacian(muEff,U) Gauss linear limited 0.333;
-    laplacian((rho*(1|A(U))),p) Gauss linear limited 0.333;
+    laplacian((rho*(1|A(U))),p_rgh) Gauss linear limited 0.333;
     laplacian(alphaEff,h) Gauss linear limited 0.333;
     laplacian(DkEff,k) Gauss linear limited 0.333;
     laplacian(DepsilonEff,epsilon) Gauss linear limited 0.333;
@@ -63,7 +61,7 @@ snGradSchemes
 fluxRequired
 {
     default         no;
-    p;
+    p_rgh;
 }
 
 // ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/bottomAir/fvSolution b/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/bottomAir/fvSolution
index ad1f5492fb46f3751c4eaef9ee2665b31fef1e0f..ed5885da9a3aeebf717fead8cd15050bb82fdd23 100644
--- a/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/bottomAir/fvSolution
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/bottomAir/fvSolution
@@ -20,14 +20,14 @@ solvers
     {
         solver          PCG
         preconditioner  DIC;
-        tolerance       1e-8;
+        tolerance       1e-7;
         relTol          0;
     }
 
-    p
+    p_rgh
     {
         solver           GAMG;
-        tolerance        1e-8;
+        tolerance        1e-7;
         relTol           0.01;
 
         smoother         GaussSeidel;
@@ -38,46 +38,39 @@ solvers
         mergeLevels      1;
     }
 
-    pFinal
+    p_rghFinal
     {
-        $p;
-        tolerance        1e-8;
+        $p_rgh;
+        tolerance        1e-7;
         relTol           0;
     }
 
-    U
+    "(U|h|k|epsilon|R)"
     {
         solver           PBiCG;
         preconditioner   DILU;
-        tolerance        1e-08;
-        relTol           0;
-    }
-
-    h
-    {
-        $U;
-        tolerance        1e-08;
+        tolerance        1e-7;
         relTol           0.1;
     }
 
-
-    "(hFinal|k|epsilon|R)"
+    "(U|h|k|epsilon|R)Final"
     {
         $U;
+        tolerance        1e-07;
+        relTol           0;
     }
 }
 
 PISO
 {
-    momentumPredictor   off;
+    momentumPredictor    off;
     nOuterCorrectors     1;
-    nCorrectors     2;
+    nCorrectors          2;
     nNonOrthogonalCorrectors 1;
-    pRefPoint       (-0.081 -0.0257 8.01);
-    pRefValue       1e5;
+    pRefPoint            (-0.081 -0.0257 8.01);
+    pRefValue            1e5;
 }
 
-
 PIMPLE
 {
     momentumPredictor   on;
@@ -85,11 +78,10 @@ PIMPLE
     nNonOrthogonalCorrectors 0;
 }
 
-
 relaxationFactors
 {
-//    h               0.9;
-//    U               0.9;
+    h               1;
+    U               1;
 }
 
 // ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/heater/decomposeParDict b/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/heater/decomposeParDict
index e53a523c9d58e8314f6eeea8b310d284b54f6cd0..8f6ced4ac38fc5f058b8c4370ed4002c036d61f6 100644
--- a/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/heater/decomposeParDict
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/heater/decomposeParDict
@@ -76,7 +76,7 @@ manualCoeffs
 //// Is the case distributed
 //distributed     yes;
 //// Per slave (so nProcs-1 entries) the directory above the case.
-//roots           
+//roots
 //(
 //    "/tmp"
 //    "/tmp"
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/heater/fvSolution b/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/heater/fvSolution
index b773399f729bcb6ee8fc6b7c458c3d82e9ea2a8f..7c00e466f4b508dc216b21874718108befe65e16 100644
--- a/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/heater/fvSolution
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/heater/fvSolution
@@ -21,6 +21,12 @@ solvers
         solver           PCG;
         preconditioner   DIC;
         tolerance        1E-06;
+        relTol           0.1;
+    }
+    TFinal
+    {
+        $T;
+        tolerance        1E-06;
         relTol           0;
     }
 }
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/leftSolid/decomposeParDict b/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/leftSolid/decomposeParDict
index e53a523c9d58e8314f6eeea8b310d284b54f6cd0..8f6ced4ac38fc5f058b8c4370ed4002c036d61f6 100644
--- a/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/leftSolid/decomposeParDict
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/leftSolid/decomposeParDict
@@ -76,7 +76,7 @@ manualCoeffs
 //// Is the case distributed
 //distributed     yes;
 //// Per slave (so nProcs-1 entries) the directory above the case.
-//roots           
+//roots
 //(
 //    "/tmp"
 //    "/tmp"
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/leftSolid/fvSolution b/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/leftSolid/fvSolution
index b773399f729bcb6ee8fc6b7c458c3d82e9ea2a8f..18b2ebb7ebf2761aad98a3b88bb550724bbd8a1e 100644
--- a/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/leftSolid/fvSolution
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/leftSolid/fvSolution
@@ -21,6 +21,13 @@ solvers
         solver           PCG;
         preconditioner   DIC;
         tolerance        1E-06;
+        relTol           0.1;
+    }
+
+    TFinal
+    {
+        $T;
+        tolerance        1E-06;
         relTol           0;
     }
 }
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/rightSolid/decomposeParDict b/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/rightSolid/decomposeParDict
index e53a523c9d58e8314f6eeea8b310d284b54f6cd0..8f6ced4ac38fc5f058b8c4370ed4002c036d61f6 100644
--- a/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/rightSolid/decomposeParDict
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/rightSolid/decomposeParDict
@@ -76,7 +76,7 @@ manualCoeffs
 //// Is the case distributed
 //distributed     yes;
 //// Per slave (so nProcs-1 entries) the directory above the case.
-//roots           
+//roots
 //(
 //    "/tmp"
 //    "/tmp"
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/rightSolid/fvSolution b/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/rightSolid/fvSolution
index b773399f729bcb6ee8fc6b7c458c3d82e9ea2a8f..18b2ebb7ebf2761aad98a3b88bb550724bbd8a1e 100644
--- a/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/rightSolid/fvSolution
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/rightSolid/fvSolution
@@ -21,6 +21,13 @@ solvers
         solver           PCG;
         preconditioner   DIC;
         tolerance        1E-06;
+        relTol           0.1;
+    }
+
+    TFinal
+    {
+        $T;
+        tolerance        1E-06;
         relTol           0;
     }
 }
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/topAir/changeDictionaryDict b/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/topAir/changeDictionaryDict
index 71114003ec106a4a9af29c5abfa86f94cf185c56..ff4ba100119b83b40959e27389c224519da39a2e 100644
--- a/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/topAir/changeDictionaryDict
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/topAir/changeDictionaryDict
@@ -18,7 +18,7 @@ dictionaryReplacement
 {
     U
     {
-        internalField   uniform (0.01 0 0);
+        internalField   uniform (0.1 0 0);
 
         boundaryField
         {
@@ -30,13 +30,13 @@ dictionaryReplacement
             minX
             {
                 type            fixedValue;
-                value           uniform ( 0.01 0 0 );
+                value           uniform ( 0.1 0 0 );
             }
             maxX
             {
                 type            inletOutlet;
                 inletValue      uniform ( 0 0 0 );
-                value           uniform ( 0 0 0 );
+                value           uniform ( 0.1 0 0 );
             }
         }
     }
@@ -127,25 +127,42 @@ dictionaryReplacement
         }
     }
 
-    p
+    p_rgh
     {
-        internalField   uniform 100000;
+        internalField   uniform 1e5;
 
         boundaryField
         {
             ".*"
             {
                 type            buoyantPressure;
-                value           1e5;
+                value           uniform 1e5;
+            }
+
+            maxX
+            {
+                type            fixedValue;
+                value           uniform 1e5;
+            }
+        }
+    }
+
+    p
+    {
+        internalField   uniform 1e5;
+
+        boundaryField
+        {
+            ".*"
+            {
+                type            calculated;
+                value           uniform 1e5;
             }
 
             maxX
             {
-                type            waveTransmissive;
-                gamma           1.4;
-                fieldInf        100000;
-                lInf            0.4;
-                value           uniform 100000;
+                type            calculated;
+                value           uniform 1e5;
             }
         }
     }
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/topAir/decomposeParDict b/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/topAir/decomposeParDict
index e53a523c9d58e8314f6eeea8b310d284b54f6cd0..8f6ced4ac38fc5f058b8c4370ed4002c036d61f6 100644
--- a/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/topAir/decomposeParDict
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/topAir/decomposeParDict
@@ -76,7 +76,7 @@ manualCoeffs
 //// Is the case distributed
 //distributed     yes;
 //// Per slave (so nProcs-1 entries) the directory above the case.
-//roots           
+//roots
 //(
 //    "/tmp"
 //    "/tmp"
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/topAir/fvSchemes b/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/topAir/fvSchemes
index 90364fcbf2309c83fb6dbf1002215f43b97d23d5..7d3e57864ac0660fb7bd646e2850a67eafb1d49d 100644
--- a/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/topAir/fvSchemes
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/topAir/fvSchemes
@@ -17,13 +17,11 @@ FoamFile
 ddtSchemes
 {
     default Euler;
-    //default CoEuler phi rho 0.1;
 }
 
 gradSchemes
 {
     default         Gauss linear;
-//    grad(U)         cellLimited Gauss linear 1;
 }
 
 divSchemes
@@ -43,7 +41,7 @@ laplacianSchemes
 {
     default         none;
     laplacian(muEff,U) Gauss linear limited 0.333;
-    laplacian((rho*(1|A(U))),p) Gauss linear limited 0.333;
+    laplacian((rho*(1|A(U))),p_rgh) Gauss linear limited 0.333;
     laplacian(alphaEff,h) Gauss linear limited 0.333;
     laplacian(DkEff,k) Gauss linear limited 0.333;
     laplacian(DepsilonEff,epsilon) Gauss linear limited 0.333;
@@ -63,7 +61,7 @@ snGradSchemes
 fluxRequired
 {
     default         no;
-    p;
+    p_rgh;
 }
 
 // ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/topAir/fvSolution b/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/topAir/fvSolution
index 28c1a63b7aa3986ac043fd9eea7b5975e68d7c5d..d2d44df4fc5f85c3a91d97b0f5402c4649e7102f 100644
--- a/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/topAir/fvSolution
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/topAir/fvSolution
@@ -20,14 +20,14 @@ solvers
     {
         solver          PCG
         preconditioner  DIC;
-        tolerance       1e-8;
+        tolerance       1e-7;
         relTol          0;
     }
 
-    p
+    p_rgh
     {
         solver           GAMG;
-        tolerance        1e-8;
+        tolerance        1e-7;
         relTol           0.01;
 
         smoother         GaussSeidel;
@@ -36,47 +36,43 @@ solvers
         nCellsInCoarsestLevel 10;
         agglomerator     faceAreaPair;
         mergeLevels      1;
+
+        maxIter          100;
     }
 
-    pFinal
+    p_rghFinal
     {
-        $p;
-        tolerance        1e-8;
+        $p_rgh;
+        tolerance        1e-7;
         relTol           0;
     }
 
-    U
+    "(U|h|k|epsilon|R)"
     {
         solver           PBiCG;
         preconditioner   DILU;
-        tolerance        1e-08;
-        relTol           0;
-    }
-
-    h
-    {
-        $U;
-        tolerance        1e-08;
+        tolerance        1e-7;
         relTol           0.1;
     }
 
-    "(hFinal|k|epsilon|R)"
+    "(U|h|k|epsilon|R)Final"
     {
         $U;
+        tolerance        1e-07;
+        relTol           0;
     }
 }
 
 PISO
 {
-    momentumPredictor   off;
+    momentumPredictor    off;
     nOuterCorrectors     1;
-    nCorrectors     2;
+    nCorrectors          2;
     nNonOrthogonalCorrectors 1;
-    pRefPoint       (-0.081 -0.0257 8.01);
-    pRefValue       1e5;
+    pRefPoint            (-0.081 -0.0257 8.01);
+    pRefValue            1e5;
 }
 
-
 PIMPLE
 {
     momentumPredictor   on;
@@ -84,11 +80,10 @@ PIMPLE
     nNonOrthogonalCorrectors 0;
 }
 
-
 relaxationFactors
 {
-//    h               0.9;
-//    U               0.9;
+    h               1;
+    U               1;
 }
 
 // ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/0/U b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/0/U
index 62058ec7acf331fb2b61fa165da191f4776f7eee..8e47fbcddf8a6d77315c7deb728c6eb1215eb68e 100644
--- a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/0/U
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/0/U
@@ -16,7 +16,7 @@ FoamFile
 
 dimensions      [0 1 -1 0 0 0 0];
 
-internalField   uniform (0.01 0 0);
+internalField   uniform (0.1 0 0);
 
 boundaryField
 {
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/0/epsilon b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/0/epsilon
index 6034a1d7bfc3f328fece0218588466e1deb72932..f4e27ad95fb84fff1121e59c121f6d2c844d1b43 100644
--- a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/0/epsilon
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/0/epsilon
@@ -10,7 +10,6 @@ FoamFile
     version     2.0;
     format      ascii;
     class       volScalarField;
-    location    "0.001";
     object      epsilon;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/0/k b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/0/k
index eb1abac09b21d07f72329ef894699f330741de42..b2509f967076c366a4ab4e410f070916983536e4 100644
--- a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/0/k
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/0/k
@@ -10,7 +10,6 @@ FoamFile
     version     2.0;
     format      ascii;
     class       volScalarField;
-    location    "0.001";
     object      k;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/0/p_rgh b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/0/p_rgh
new file mode 100644
index 0000000000000000000000000000000000000000..638c8c7d7cd472cf343b4fc5091d9a8f1bf038ae
--- /dev/null
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/0/p_rgh
@@ -0,0 +1,29 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      p_rgh;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -2 0 0 0 0];
+
+internalField   uniform 1e5;
+
+boundaryField
+{
+    ".*"
+    {
+        type            calculated;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/bottomAir/changeDictionaryDict b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/bottomAir/changeDictionaryDict
index c36ea6b8039bbf390c437d1875aa6eebe84cf63f..6faab27ee56705411dc2abcf18838952a6f25fa1 100644
--- a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/bottomAir/changeDictionaryDict
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/bottomAir/changeDictionaryDict
@@ -92,16 +92,30 @@ dictionaryReplacement
         }
     }
 
-    p
+    p_rgh
     {
-        internalField   uniform 100000;
+        internalField   uniform 1e5;
 
         boundaryField
         {
             ".*"
             {
                 type            buoyantPressure;
-                value           1e5;
+                value           uniform 1e5;
+            }
+        }
+    }
+
+    p
+    {
+        internalField   uniform 1e5;
+
+        boundaryField
+        {
+            ".*"
+            {
+                type            calculated;
+                value           uniform 1e5;
             }
         }
     }
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/bottomAir/decomposeParDict b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/bottomAir/decomposeParDict
index e53a523c9d58e8314f6eeea8b310d284b54f6cd0..8f6ced4ac38fc5f058b8c4370ed4002c036d61f6 100644
--- a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/bottomAir/decomposeParDict
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/bottomAir/decomposeParDict
@@ -76,7 +76,7 @@ manualCoeffs
 //// Is the case distributed
 //distributed     yes;
 //// Per slave (so nProcs-1 entries) the directory above the case.
-//roots           
+//roots
 //(
 //    "/tmp"
 //    "/tmp"
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/bottomAir/fvSchemes b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/bottomAir/fvSchemes
index a06836c2b3bcf34f7c2d7df8a6bb39f7f0b3f5a0..c136ffea3116ad29d3f67eb00de18abd99b3ef5b 100644
--- a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/bottomAir/fvSchemes
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/bottomAir/fvSchemes
@@ -28,7 +28,7 @@ divSchemes
 {
     default         none;
     div(phi,U)      Gauss upwind;
-    div(phiU,p)     Gauss linear;
+    div(U,p)        Gauss linear;
     div(phi,h)      Gauss upwind;
     div(phi,k)      Gauss upwind;
     div(phi,epsilon) Gauss upwind;
@@ -41,7 +41,7 @@ laplacianSchemes
 {
     default         none;
     laplacian(muEff,U) Gauss linear limited 0.333;
-    laplacian((rho*(1|A(U))),p) Gauss linear limited 0.333;
+    laplacian((rho*(1|A(U))),p_rgh) Gauss linear limited 0.333;
     laplacian(alphaEff,h) Gauss linear limited 0.333;
     laplacian(DkEff,k) Gauss linear limited 0.333;
     laplacian(DepsilonEff,epsilon) Gauss linear limited 0.333;
@@ -61,7 +61,7 @@ snGradSchemes
 fluxRequired
 {
     default         no;
-    p;
+    p_rgh;
 }
 
 // ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/bottomAir/fvSolution b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/bottomAir/fvSolution
index 842419ac72ba1d61e74173ee300cafd9fb171ee5..eff94d930fcb9875b4853c30e7592634920cd851 100644
--- a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/bottomAir/fvSolution
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/bottomAir/fvSolution
@@ -20,14 +20,14 @@ solvers
     {
         solver          PCG
         preconditioner  DIC;
-        tolerance       1e-8;
+        tolerance       1e-7;
         relTol          0;
     }
 
-    p
+    p_rgh
     {
         solver           GAMG;
-        tolerance        1e-8;
+        tolerance        1e-7;
         relTol           0.01;
 
         smoother         GaussSeidel;
@@ -42,22 +42,25 @@ solvers
     {
         solver           PBiCG;
         preconditioner   DILU;
-        tolerance        1e-08;
+        tolerance        1e-7;
         relTol           0.1;
     }
 }
 
 SIMPLE
 {
+    momentumPredictor on;
     nNonOrthogonalCorrectors 0;
     pRefCell                0;
     pRefValue               100000;
+    rhoMin      rhoMin [1 -3 0 0 0] 0.2;
+    rhoMax      rhoMax [1 -3 0 0 0] 2;
 }
 
 relaxationFactors
 {
     rho         1;
-    p           0.7;
+    p_rgh       0.7;
     U           0.3;
     h           0.7;
     nuTilda     0.7;
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/controlDict b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/controlDict
index 0c78e811c64e188c09e6d22ea16aa332c7024503..32e60a1dd1db695e64f63501af2c5e5de42616cc 100644
--- a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/controlDict
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/controlDict
@@ -17,7 +17,7 @@ FoamFile
 
 libs            ("libcompressibleTurbulenceModel.so" "libcompressibleRASModels.so");
 
-application     chtMultiRegionFoam;
+application     chtMultiRegionSimpleFoam;
 
 startFrom       startTime;
 
@@ -25,14 +25,15 @@ startTime       0;
 
 stopAt          endTime;
 
-endTime         1000;
+endTime         2000;
 
 deltaT          1;
 
 writeControl    timeStep;
+
 writeInterval   100;
 
-purgeWrite      0;
+purgeWrite      5;
 
 writeFormat     ascii;
 
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/decomposeParDict b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/decomposeParDict
index e53a523c9d58e8314f6eeea8b310d284b54f6cd0..8f6ced4ac38fc5f058b8c4370ed4002c036d61f6 100644
--- a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/decomposeParDict
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/decomposeParDict
@@ -76,7 +76,7 @@ manualCoeffs
 //// Is the case distributed
 //distributed     yes;
 //// Per slave (so nProcs-1 entries) the directory above the case.
-//roots           
+//roots
 //(
 //    "/tmp"
 //    "/tmp"
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/heater/decomposeParDict b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/heater/decomposeParDict
index e53a523c9d58e8314f6eeea8b310d284b54f6cd0..8f6ced4ac38fc5f058b8c4370ed4002c036d61f6 100644
--- a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/heater/decomposeParDict
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/heater/decomposeParDict
@@ -76,7 +76,7 @@ manualCoeffs
 //// Is the case distributed
 //distributed     yes;
 //// Per slave (so nProcs-1 entries) the directory above the case.
-//roots           
+//roots
 //(
 //    "/tmp"
 //    "/tmp"
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/leftSolid/decomposeParDict b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/leftSolid/decomposeParDict
index e53a523c9d58e8314f6eeea8b310d284b54f6cd0..8f6ced4ac38fc5f058b8c4370ed4002c036d61f6 100644
--- a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/leftSolid/decomposeParDict
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/leftSolid/decomposeParDict
@@ -76,7 +76,7 @@ manualCoeffs
 //// Is the case distributed
 //distributed     yes;
 //// Per slave (so nProcs-1 entries) the directory above the case.
-//roots           
+//roots
 //(
 //    "/tmp"
 //    "/tmp"
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/rightSolid/decomposeParDict b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/rightSolid/decomposeParDict
index e53a523c9d58e8314f6eeea8b310d284b54f6cd0..8f6ced4ac38fc5f058b8c4370ed4002c036d61f6 100644
--- a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/rightSolid/decomposeParDict
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/rightSolid/decomposeParDict
@@ -76,7 +76,7 @@ manualCoeffs
 //// Is the case distributed
 //distributed     yes;
 //// Per slave (so nProcs-1 entries) the directory above the case.
-//roots           
+//roots
 //(
 //    "/tmp"
 //    "/tmp"
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/topAir/changeDictionaryDict b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/topAir/changeDictionaryDict
index cf0e2f131c205d67073d5777a4cd6f6f500677a4..ff4ba100119b83b40959e27389c224519da39a2e 100644
--- a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/topAir/changeDictionaryDict
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/topAir/changeDictionaryDict
@@ -18,7 +18,7 @@ dictionaryReplacement
 {
     U
     {
-        internalField   uniform (0.01 0 0);
+        internalField   uniform (0.1 0 0);
 
         boundaryField
         {
@@ -30,13 +30,13 @@ dictionaryReplacement
             minX
             {
                 type            fixedValue;
-                value           uniform ( 0.01 0 0 );
+                value           uniform ( 0.1 0 0 );
             }
             maxX
             {
                 type            inletOutlet;
                 inletValue      uniform ( 0 0 0 );
-                value           uniform ( 0 0 0 );
+                value           uniform ( 0.1 0 0 );
             }
         }
     }
@@ -127,22 +127,42 @@ dictionaryReplacement
         }
     }
 
-    p
+    p_rgh
     {
-        internalField   uniform 100000;
+        internalField   uniform 1e5;
 
         boundaryField
         {
             ".*"
             {
                 type            buoyantPressure;
-                value           1e5;
+                value           uniform 1e5;
             }
 
             maxX
             {
                 type            fixedValue;
-                value           uniform 100000;
+                value           uniform 1e5;
+            }
+        }
+    }
+
+    p
+    {
+        internalField   uniform 1e5;
+
+        boundaryField
+        {
+            ".*"
+            {
+                type            calculated;
+                value           uniform 1e5;
+            }
+
+            maxX
+            {
+                type            calculated;
+                value           uniform 1e5;
             }
         }
     }
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/topAir/decomposeParDict b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/topAir/decomposeParDict
index e53a523c9d58e8314f6eeea8b310d284b54f6cd0..8f6ced4ac38fc5f058b8c4370ed4002c036d61f6 100644
--- a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/topAir/decomposeParDict
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/topAir/decomposeParDict
@@ -76,7 +76,7 @@ manualCoeffs
 //// Is the case distributed
 //distributed     yes;
 //// Per slave (so nProcs-1 entries) the directory above the case.
-//roots           
+//roots
 //(
 //    "/tmp"
 //    "/tmp"
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/topAir/fvSchemes b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/topAir/fvSchemes
index a06836c2b3bcf34f7c2d7df8a6bb39f7f0b3f5a0..c136ffea3116ad29d3f67eb00de18abd99b3ef5b 100644
--- a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/topAir/fvSchemes
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/topAir/fvSchemes
@@ -28,7 +28,7 @@ divSchemes
 {
     default         none;
     div(phi,U)      Gauss upwind;
-    div(phiU,p)     Gauss linear;
+    div(U,p)        Gauss linear;
     div(phi,h)      Gauss upwind;
     div(phi,k)      Gauss upwind;
     div(phi,epsilon) Gauss upwind;
@@ -41,7 +41,7 @@ laplacianSchemes
 {
     default         none;
     laplacian(muEff,U) Gauss linear limited 0.333;
-    laplacian((rho*(1|A(U))),p) Gauss linear limited 0.333;
+    laplacian((rho*(1|A(U))),p_rgh) Gauss linear limited 0.333;
     laplacian(alphaEff,h) Gauss linear limited 0.333;
     laplacian(DkEff,k) Gauss linear limited 0.333;
     laplacian(DepsilonEff,epsilon) Gauss linear limited 0.333;
@@ -61,7 +61,7 @@ snGradSchemes
 fluxRequired
 {
     default         no;
-    p;
+    p_rgh;
 }
 
 // ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/topAir/fvSolution b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/topAir/fvSolution
index 842419ac72ba1d61e74173ee300cafd9fb171ee5..6e31656bb367a276b1ca28c2b355e60e10725968 100644
--- a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/topAir/fvSolution
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/topAir/fvSolution
@@ -20,14 +20,14 @@ solvers
     {
         solver          PCG
         preconditioner  DIC;
-        tolerance       1e-8;
+        tolerance       1e-7;
         relTol          0;
     }
 
-    p
+    p_rgh
     {
         solver           GAMG;
-        tolerance        1e-8;
+        tolerance        1e-7;
         relTol           0.01;
 
         smoother         GaussSeidel;
@@ -36,28 +36,33 @@ solvers
         nCellsInCoarsestLevel 10;
         agglomerator     faceAreaPair;
         mergeLevels      1;
+
+        maxIter          100;
     }
 
     "(U|h|k|epsilon)"
     {
         solver           PBiCG;
         preconditioner   DILU;
-        tolerance        1e-08;
+        tolerance        1e-7;
         relTol           0.1;
     }
 }
 
 SIMPLE
 {
+    momentumPredictor on;
     nNonOrthogonalCorrectors 0;
     pRefCell                0;
     pRefValue               100000;
+    rhoMin      rhoMin [1 -3 0 0 0] 0.2;
+    rhoMax      rhoMax [1 -3 0 0 0] 2;
 }
 
 relaxationFactors
 {
     rho         1;
-    p           0.7;
+    p_rgh       0.7;
     U           0.3;
     h           0.7;
     nuTilda     0.7;