diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/TEqn.H b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/TEqn.H
index dbfc61739f..9a835792a4 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 35387f4ff4..20a05e5cd4 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 ebf68d409c..54519906a4 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 23d20cfa96..342af079d8 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 21be033f9b..60828e18dc 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 76cc4da8ba..477a322833 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 5664bb2154..50149ed360 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 d4878d063d..8c6a3f7671 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 0075c1e3ed..cb4c4d34f6 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 b8ac5595e4..e39d0bb803 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 3125cc3ffa..94537508b3 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 a1c3dbfed5..e625f122a3 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 06fa5d12a0..64cc110cec 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 cf23150332..d6fa9acee9 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();
-    p_rgh = p - rho*gh;
-
-    label p_rghRefCell = 0;
-    scalar p_rghRefValue = 0.0;
-    setRefCell
+    Info<< "Reading field p_rgh\n" << endl;
+    volScalarField p_rgh
     (
-        p_rgh,
-        mesh.solutionDict().subDict("SIMPLE"),
-        p_rghRefCell,
-        p_rghRefValue
+        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 pRefCell = 0;
     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)
-        );
-    }
+    setRefCell
+    (
+        p,
+        p_rgh,
+        mesh.solutionDict().subDict("SIMPLE"),
+        pRefCell,
+        pRefValue
+    );
 
     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 3088b42c69..f1a62d4770 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 0dbe80c67c..6c35536333 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 a81b0faaf3..36b3f1c3b0 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 d5ae07ff9b..68b29fd0c2 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 2d2ccf9a5e..167484a7c7 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,74 @@
             ).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 cef0122928..7869046fa0 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 9462dc2c39..317d377f1b 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 8491056ea3..cac7750e97 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 5198941366..2b6de83ca3 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 65467f8086..119eeb4fd7 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 1826a77217..b0a7f95912 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 e070537db2..16ba36f7dc 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 a264b68fe5..9bf3c4acc5 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 1f8cdb82a3..0000000000
--- 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 50a3011484..89aaec4737 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 86dd4344c1..b36cf89d34 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 f63e85458e..94fd01a2ba 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 be13e670c3..837305659e 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 ce8b1d0f40..05b26d1626 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 b02a7c61e8..5851cdcf95 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 5f373cc224..f32be6214b 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 ced6b64004..40eb26709e 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 b31ad45ed4..d3efde0b34 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/doc/changes/splitCyclic.txt b/doc/changes/splitCyclic.txt
index 814dfe0850..3d0ffb71cf 100644
--- a/doc/changes/splitCyclic.txt
+++ b/doc/changes/splitCyclic.txt
@@ -20,7 +20,8 @@ The disadvantages:
 - a patch-wise loop now might need to store data to go to the neighbour half
 since it is no longer handled in a single patch.
 - decomposed cyclics now require overlapping communications so will
-only work in non-blocking mode. Hence the underlying message passing library
+only work in 'nonBlocking' mode or 'blocking' (=buffered) mode but not
+in 'scheduled' mode. The underlying message passing library
 will require overlapping communications with message tags.
 - it is quite a code-change and there might be some oversights.
 - once converted (see foamUpgradeCyclics below) cases are not backwards
@@ -103,19 +104,14 @@ type 'processorCyclic'.
 
 
 - processor patches use overlapping communication using a different message
-tag. This maps straight through into the MPI message tag.
-See processorCyclicPolyPatch::tag(). This needs to be calculated the
-same on both sides so is calculated as
-        Pstream::nProcs()*max(myProcNo, neighbProcNo)
-      + min(myProcNo, neighbProcNo)
-which is
-- unique
-- commutative
-- does not interfere with the default tag (= 1)
+tag. This maps straight through into the MPI message tag. Each processor
+'interface' (processorPolyPatch, processorFvPatch, etc.) has a 'tag()'
+to use for communication.
 
 
 - when constructing a GeometricField from a dictionary it will explicitly
 check for non-existing entries for cyclic patches and exit with an error message
-warning to run foamUpgradeCyclics.
+warning to run foamUpgradeCyclics. (1.7.x will check if you are trying
+to run a case which has split cyclics)
 
 
diff --git a/src/OpenFOAM/Make/files b/src/OpenFOAM/Make/files
index d9548b67b1..94bc86f2ef 100644
--- a/src/OpenFOAM/Make/files
+++ b/src/OpenFOAM/Make/files
@@ -296,6 +296,7 @@ primitiveShapes = meshes/primitiveShapes
 $(primitiveShapes)/line/line.C
 $(primitiveShapes)/plane/plane.C
 $(primitiveShapes)/triangle/intersection.C
+$(primitiveShapes)/objectHit/pointIndexHitIOList.C
 
 meshShapes = meshes/meshShapes
 $(meshShapes)/edge/edge.C
diff --git a/src/OpenFOAM/meshes/meshShapes/face/face.H b/src/OpenFOAM/meshes/meshShapes/face/face.H
index 81ce90d8a9..dd5fc4a832 100644
--- a/src/OpenFOAM/meshes/meshShapes/face/face.H
+++ b/src/OpenFOAM/meshes/meshShapes/face/face.H
@@ -126,6 +126,14 @@ class face
 
 public:
 
+    //- Return types for classify
+    enum proxType
+    {
+        NONE,
+        POINT,  // Close to point
+        EDGE    // Close to edge
+    };
+
     // Static data members
 
         static const char* const typeName;
@@ -249,6 +257,20 @@ public:
             const pointField& meshPoints
         ) const;
 
+        //- Return nearest point to face and classify it:
+        //  + near point (nearType=POINT, nearLabel=0, 1, 2)
+        //  + near edge (nearType=EDGE, nearLabel=0, 1, 2)
+        //    Note: edges are counted from starting vertex so
+        //    e.g. edge n is from f[n] to f[0], where the face has n + 1
+        //    points
+        pointHit nearestPointClassify
+        (
+            const point& p,
+            const pointField& meshPoints,
+            label& nearType,
+            label& nearLabel
+        ) const;
+
         //- Return contact sphere diameter
         scalar contactSphereDiameter
         (
diff --git a/src/OpenFOAM/meshes/meshShapes/face/faceIntersection.C b/src/OpenFOAM/meshes/meshShapes/face/faceIntersection.C
index d4184e766a..f9f02eb612 100644
--- a/src/OpenFOAM/meshes/meshShapes/face/faceIntersection.C
+++ b/src/OpenFOAM/meshes/meshShapes/face/faceIntersection.C
@@ -181,6 +181,22 @@ Foam::pointHit Foam::face::nearestPoint
     const point& p,
     const pointField& meshPoints
 ) const
+{
+    // Dummy labels
+    label nearType = -1;
+    label nearLabel = -1;
+
+    return nearestPointClassify(p, meshPoints, nearType, nearLabel);
+}
+
+
+Foam::pointHit Foam::face::nearestPointClassify
+(
+    const point& p,
+    const pointField& meshPoints,
+    label& nearType,
+    label& nearLabel
+) const
 {
     const face& f = *this;
     point ctr = centre(meshPoints);
@@ -188,6 +204,9 @@ Foam::pointHit Foam::face::nearestPoint
     // Initialize to miss, distance=GREAT
     pointHit nearest(p);
 
+    nearType = -1;
+    nearLabel = -1;
+
     label nPoints = f.size();
 
     point nextPoint = ctr;
@@ -196,8 +215,10 @@ Foam::pointHit Foam::face::nearestPoint
     {
         nextPoint = meshPoints[f[fcIndex(pI)]];
 
+        label tmpNearType = -1;
+        label tmpNearLabel = -1;
+
         // Note: for best accuracy, centre point always comes last
-        //
         triPointRef tri
         (
             meshPoints[f[pI]],
@@ -205,12 +226,42 @@ Foam::pointHit Foam::face::nearestPoint
             ctr
         );
 
-        pointHit curHit = tri.nearestPoint(p);
+        pointHit curHit = tri.nearestPointClassify
+        (
+            p,
+            tmpNearType,
+            tmpNearLabel
+        );
 
         if (Foam::mag(curHit.distance()) < Foam::mag(nearest.distance()))
         {
             nearest.setDistance(curHit.distance());
 
+            // Assume at first that the near type is NONE on the
+            // triangle (i.e. on the face of the triangle) then it is
+            // therefore also for the face.
+
+            nearType = NONE;
+
+            if (tmpNearType == triPointRef::EDGE && tmpNearLabel == 0)
+            {
+                // If the triangle edge label is 0, then this is also
+                // an edge of the face, if not, it is on the face
+
+                nearType = EDGE;
+
+                nearLabel = pI;
+            }
+            else if (tmpNearType == triPointRef::POINT && tmpNearLabel < 2)
+            {
+                // If the triangle point label is 0 or 1, then this is
+                // also a point of the face, if not, it is on the face
+
+                nearType = POINT;
+
+                nearLabel = pI + tmpNearLabel;
+            }
+
             if (curHit.hit())
             {
                 nearest.setHit();
diff --git a/src/meshTools/octree/PointIndexHit.H b/src/OpenFOAM/meshes/primitiveShapes/objectHit/PointIndexHit.H
similarity index 100%
rename from src/meshTools/octree/PointIndexHit.H
rename to src/OpenFOAM/meshes/primitiveShapes/objectHit/PointIndexHit.H
diff --git a/src/meshTools/octree/pointHitSort.H b/src/OpenFOAM/meshes/primitiveShapes/objectHit/pointHitSort.H
similarity index 100%
rename from src/meshTools/octree/pointHitSort.H
rename to src/OpenFOAM/meshes/primitiveShapes/objectHit/pointHitSort.H
diff --git a/src/meshTools/octree/pointIndexHit.H b/src/OpenFOAM/meshes/primitiveShapes/objectHit/pointIndexHit.H
similarity index 100%
rename from src/meshTools/octree/pointIndexHit.H
rename to src/OpenFOAM/meshes/primitiveShapes/objectHit/pointIndexHit.H
diff --git a/src/meshTools/octree/pointIndexHitIOList.C b/src/OpenFOAM/meshes/primitiveShapes/objectHit/pointIndexHitIOList.C
similarity index 100%
rename from src/meshTools/octree/pointIndexHitIOList.C
rename to src/OpenFOAM/meshes/primitiveShapes/objectHit/pointIndexHitIOList.C
diff --git a/src/meshTools/octree/pointIndexHitIOList.H b/src/OpenFOAM/meshes/primitiveShapes/objectHit/pointIndexHitIOList.H
similarity index 100%
rename from src/meshTools/octree/pointIndexHitIOList.H
rename to src/OpenFOAM/meshes/primitiveShapes/objectHit/pointIndexHitIOList.H
diff --git a/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H
index 64d9318306..abb5d7d8e4 100644
--- a/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H
+++ b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H
@@ -79,21 +79,6 @@ class triangle
 
         PointRef a_, b_, c_;
 
-    // Private Member Functions
-
-        //- Fast distance to triangle calculation. From
-        //  "Distance Between Point and Trangle in 3D"
-        //  David Eberly, Magic Software Inc. Aug. 2002.
-        //  Works on function Q giving distance to point and tries to
-        //  minimize this.
-        static pointHit nearestPoint
-        (
-            const Point& baseVertex,
-            const vector& E0,
-            const vector& E1,
-            const point& P
-        );
-
 
 public:
 
@@ -202,24 +187,27 @@ public:
                 const scalar tol = 0.0
             ) const;
 
-            //- Return nearest point to p on triangle
-            inline pointHit nearestPoint
-            (
-                const point& p
-            ) const;
-
-            //- Classify point in triangle plane w.r.t. triangle edges.
-            //  - inside (true returned)/outside (false returned)
-            //  - near point (nearType=POINT, nearLabel=0, 1, 2)
-            //  - near edge (nearType=EDGE, nearLabel=0, 1, 2)
+            //- Find the nearest point to p on the triangle and classify it:
+            //  + near point (nearType=POINT, nearLabel=0, 1, 2)
+            //  + near edge (nearType=EDGE, nearLabel=0, 1, 2)
             //    Note: edges are counted from starting
             //    vertex so e.g. edge 2 is from f[2] to f[0]
-            //  tol is fraction to account for truncation error. Is only used
-            //  when comparing normalized (0..1) numbers.
+            pointHit nearestPointClassify
+            (
+                const point& p,
+                label& nearType,
+                label& nearLabel
+            ) const;
+
+            //- Return nearest point to p on triangle
+            inline pointHit nearestPoint(const point& p) const;
+
+            //- Classify nearest point to p in triangle plane
+            //  w.r.t. triangle edges and points.  Returns inside
+            //  (true)/outside (false).
             bool classify
             (
                 const point& p,
-                const scalar tol,
                 label& nearType,
                 label& nearLabel
             ) const;
diff --git a/src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H
index 10267a923f..7840d8ee7b 100644
--- a/src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H
+++ b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H
@@ -32,158 +32,6 @@ License
 namespace Foam
 {
 
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-template<class Point, class PointRef>
-pointHit triangle<Point, PointRef>::nearestPoint
-(
-    const Point& baseVertex,
-    const vector& E0,
-    const vector& E1,
-    const point& P
-)
-{
-    // Distance vector
-    const vector D(baseVertex - P);
-
-    // Some geometrical factors
-    const scalar a = E0 & E0;
-    const scalar b = E0 & E1;
-    const scalar c = E1 & E1;
-
-    // Precalculate distance factors
-    const scalar d = E0 & D;
-    const scalar e = E1 & D;
-    const scalar f = D & D;
-
-    // Do classification
-    const scalar det = a*c - b*b;
-    scalar s = b*e - c*d;
-    scalar t = b*d - a*e;
-
-    bool inside = false;
-
-    if (s+t < det)
-    {
-        if (s < 0)
-        {
-            if (t < 0)
-            {
-                // Region 4
-                if (e > 0)
-                {
-                    // min on edge t = 0
-                    t = 0;
-                    s = (d >= 0 ? 0 : (-d >= a ? 1 : -d/a));
-                }
-                else
-                {
-                    // min on edge s=0
-                    s = 0;
-                    t = (e >= 0 ? 0 : (-e >= c ? 1 : -e/c));
-                }
-            }
-            else
-            {
-                // Region 3. Min on edge s = 0
-                s = 0;
-                t = (e >= 0 ? 0 : (-e >= c ? 1 : -e/c));
-            }
-        }
-        else if (t < 0)
-        {
-            // Region 5
-            t = 0;
-            s = (d >= 0 ? 0 : (-d >= a ? 1 : -d/a));
-        }
-        else
-        {
-            // Region 0
-            const scalar invDet = 1/det;
-            s *= invDet;
-            t *= invDet;
-
-            inside = true;
-        }
-    }
-    else
-    {
-        if (s < 0)
-        {
-            // Region 2
-            const scalar tmp0 = b + d;
-            const scalar tmp1 = c + e;
-            if (tmp1 > tmp0)
-            {
-                // min on edge s+t=1
-                const scalar numer = tmp1 - tmp0;
-                const scalar denom = a-2*b+c;
-                s = (numer >= denom ? 1 : numer/denom);
-                t = 1 - s;
-            }
-            else
-            {
-                // min on edge s=0
-                s = 0;
-                t = (tmp1 <= 0 ? 1 : (e >= 0 ? 0 : - e/c));
-            }
-        }
-        else if (t < 0)
-        {
-            // Region 6
-            const scalar tmp0 = b + d;
-            const scalar tmp1 = c + e;
-            if (tmp1 > tmp0)
-            {
-                // min on edge s+t=1
-                const scalar numer = tmp1 - tmp0;
-                const scalar denom = a-2*b+c;
-                s = (numer >= denom ? 1 : numer/denom);
-                t = 1 - s;
-            }
-            else
-            {
-                // min on edge t=0
-                t = 0;
-                s = (tmp1 <= 0 ? 1 : (d >= 0 ? 0 : - d/a));
-            }
-        }
-        else
-        {
-            // Region 1
-            const scalar numer = c+e-(b+d);
-            if (numer <= 0)
-            {
-                s = 0;
-            }
-            else
-            {
-                const scalar denom = a-2*b+c;
-                s = (numer >= denom ? 1 : numer/denom);
-            }
-        }
-
-        t = 1 - s;
-    }
-
-    // Calculate distance.
-    // Note: Foam::mag used since truncation error causes negative distances
-    // with points very close to one of the triangle vertices.
-    // (Up to -2.77556e-14 seen). Could use +SMALL but that not large enough.
-
-    return pointHit
-    (
-        inside,
-        baseVertex + s*E0 + t*E1,
-        Foam::sqrt
-        (
-            Foam::mag(a*s*s + 2*b*s*t + c*t*t + 2*d*s + 2*e*t + f)
-        ),
-        !inside
-    );
-}
-
-
 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
 
 template<class Point, class PointRef>
@@ -247,7 +95,7 @@ inline Point triangle<Point, PointRef>::centre() const
 template<class Point, class PointRef>
 inline scalar triangle<Point, PointRef>::mag() const
 {
-    return ::Foam::mag(normal());
+    return Foam::mag(normal());
 }
 
 
@@ -536,7 +384,7 @@ inline pointHit triangle<Point, PointRef>::ray
         inter.setMiss(eligible);
 
         // The miss point is the nearest point on the triangle
-        inter.setPoint(nearestPoint(a_, E0, E1, p).rawPoint());
+        inter.setPoint(nearestPoint(p).rawPoint());
 
         // The distance to the miss is the distance between the
         // original point and plane of intersection
@@ -633,18 +481,130 @@ inline pointHit triangle<Point, PointRef>::intersection
 }
 
 
+template<class Point, class PointRef>
+pointHit triangle<Point, PointRef>::nearestPointClassify
+(
+    const point& p,
+    label& nearType,
+    label& nearLabel
+) const
+{
+    // Adapted from:
+    // Real-time collision detection, Christer Ericson, 2005, 136-142
+
+    // Check if P in vertex region outside A
+    vector ab = b_ - a_;
+    vector ac = c_ - a_;
+    vector ap = p - a_;
+
+    scalar d1 = ab & ap;
+    scalar d2 = ac & ap;
+
+    if (d1 <= 0.0 && d2 <= 0.0)
+    {
+        // barycentric coordinates (1, 0, 0)
+
+        nearType = POINT;
+        nearLabel = 0;
+        return pointHit(false, a_, Foam::mag(a_ - p), true);
+    }
+
+    // Check if P in vertex region outside B
+    vector bp = p - b_;
+    scalar d3 = ab & bp;
+    scalar d4 = ac & bp;
+
+    if (d3 >= 0.0 && d4 <= d3)
+    {
+        // barycentric coordinates (0, 1, 0)
+
+        nearType = POINT;
+        nearLabel = 1;
+        return pointHit(false, b_, Foam::mag(b_ - p), true);
+    }
+
+    // Check if P in edge region of AB, if so return projection of P onto AB
+    scalar vc = d1*d4 - d3*d2;
+
+    if (vc <= 0.0 && d1 >= 0.0 && d3 <= 0.0)
+    {
+        // barycentric coordinates (1-v, v, 0)
+        scalar v = d1/(d1 - d3);
+
+        point nearPt =  a_ + v*ab;
+        nearType = EDGE;
+        nearLabel = 0;
+        return pointHit(false, nearPt, Foam::mag(nearPt - p), true);
+    }
+
+    // Check if P in vertex region outside C
+    vector cp = p - c_;
+    scalar d5 = ab & cp;
+    scalar d6 = ac & cp;
+
+    if (d6 >= 0.0 && d5 <= d6)
+    {
+        // barycentric coordinates (0, 0, 1)
+
+        nearType = POINT;
+        nearLabel = 2;
+        return pointHit(false, c_, Foam::mag(c_ - p), true);
+    }
+
+    // Check if P in edge region of AC, if so return projection of P onto AC
+    scalar vb = d5*d2 - d1*d6;
+
+    if (vb <= 0.0 && d2 >= 0.0 && d6 <= 0.0)
+    {
+        // barycentric coordinates (1-w, 0, w)
+        scalar w = d2/(d2 - d6);
+
+        point nearPt = a_ + w*ac;
+        nearType = EDGE;
+        nearLabel = 2;
+        return pointHit(false, nearPt, Foam::mag(nearPt - p), true);
+    }
+
+    // Check if P in edge region of BC, if so return projection of P onto BC
+    scalar va = d3*d6 - d5*d4;
+
+    if (va <= 0.0 && (d4 - d3) >= 0.0 && (d5 - d6) >= 0.0)
+    {
+        // barycentric coordinates (0, 1-w, w)
+        scalar w = (d4 - d3)/((d4 - d3) + (d5 - d6));
+
+        point nearPt = b_ + w*(c_ - b_);
+        nearType = EDGE;
+        nearLabel = 1;
+        return pointHit(false, nearPt, Foam::mag(nearPt - p), true);
+    }
+
+    // P inside face region. Compute Q through its barycentric
+    // coordinates (u, v, w)
+    scalar denom = 1.0/(va + vb + vc);
+    scalar v = vb * denom;
+    scalar w = vc * denom;
+
+    // = u*a + v*b + w*c, u = va*denom = 1.0 - v - w
+
+    point nearPt = a_ + ab*v + ac*w;
+    nearType = NONE,
+    nearLabel = -1;
+    return pointHit(true, nearPt, Foam::mag(nearPt - p), false);
+}
+
+
 template<class Point, class PointRef>
 inline pointHit triangle<Point, PointRef>::nearestPoint
 (
     const point& p
 ) const
 {
-    // Express triangle in terms of baseVertex (point a_) and
-    // two edge vectors
-    vector E0 = b_ - a_;
-    vector E1 = c_ - a_;
+    // Dummy labels
+    label nearType = -1;
+    label nearLabel = -1;
 
-    return nearestPoint(a_, E0, E1, p);
+    return nearestPointClassify(p, nearType, nearLabel);
 }
 
 
@@ -652,160 +612,14 @@ template<class Point, class PointRef>
 inline bool triangle<Point, PointRef>::classify
 (
     const point& p,
-    const scalar tol,
     label& nearType,
     label& nearLabel
 ) const
 {
-    const vector E0 = b_ - a_;
-    const vector E1 = c_ - a_;
-    const vector n = 0.5*(E0 ^ E1);
-
-    // Get largest component of normal
-    scalar magX = Foam::mag(n.x());
-    scalar magY = Foam::mag(n.y());
-    scalar magZ = Foam::mag(n.z());
-
-    label i0 = -1;
-    if ((magX >= magY) && (magX >= magZ))
-    {
-        i0 = 0;
-    }
-    else if ((magY >= magX) && (magY >= magZ))
-    {
-        i0 = 1;
-    }
-    else
-    {
-        i0 = 2;
-    }
-
-    // Get other components
-    label i1 = (i0 + 1) % 3;
-    label i2 = (i1 + 1) % 3;
-
-
-    scalar u1 = E0[i1];
-    scalar v1 = E0[i2];
-
-    scalar u2 = E1[i1];
-    scalar v2 = E1[i2];
-
-    scalar det = v2*u1 - u2*v1;
-
-    scalar u0 = p[i1] - a_[i1];
-    scalar v0 = p[i2] - a_[i2];
-
-    scalar alpha = 0;
-    scalar beta = 0;
-
-    bool hit = false;
-
-    if (Foam::mag(u1) < ROOTVSMALL)
-    {
-        beta = u0/u2;
-
-        alpha = (v0 - beta*v2)/v1;
-
-        hit = ((alpha >= 0) && ((alpha + beta) <= 1));
-    }
-    else
-    {
-        beta = (v0*u1 - u0*v1)/det;
-
-        alpha = (u0 - beta*u2)/u1;
-
-        hit = ((alpha >= 0) && ((alpha + beta) <= 1));
-    }
-
-    //
-    // Now alpha, beta are the coordinates in the triangle local coordinate
-    // system E0, E1
-    //
-
-    //Pout<< "alpha:" << alpha << endl;
-    //Pout<< "beta:" << beta << endl;
-    //Pout<< "hit:" << hit << endl;
-    //Pout<< "tol:" << tol << endl;
-
-    if (hit)
-    {
-        // alpha,beta might get negative due to precision errors
-        alpha = max(0.0, min(1.0, alpha));
-        beta = max(0.0, min(1.0, beta));
-    }
-
-    nearType = NONE;
-    nearLabel = -1;
-
-    if (Foam::mag(alpha+beta-1) <= tol)
-    {
-        // On edge between vert 1-2 (=edge 1)
-
-        if (Foam::mag(alpha) <= tol)
-        {
-            nearType = POINT;
-            nearLabel = 2;
-        }
-        else if (Foam::mag(beta) <= tol)
-        {
-            nearType = POINT;
-            nearLabel = 1;
-        }
-        else if ((alpha >= 0) && (alpha <= 1) && (beta >= 0) && (beta <= 1))
-        {
-            nearType = EDGE;
-            nearLabel = 1;
-        }
-    }
-    else if (Foam::mag(alpha) <= tol)
-    {
-        // On edge between vert 2-0 (=edge 2)
-
-        if (Foam::mag(beta) <= tol)
-        {
-            nearType = POINT;
-            nearLabel = 0;
-        }
-        else if (Foam::mag(beta-1) <= tol)
-        {
-            nearType = POINT;
-            nearLabel = 2;
-        }
-        else if ((beta >= 0) && (beta <= 1))
-        {
-            nearType = EDGE;
-            nearLabel = 2;
-        }
-    }
-    else if (Foam::mag(beta) <= tol)
-    {
-        // On edge between vert 0-1 (= edge 0)
-
-        if (Foam::mag(alpha) <= tol)
-        {
-            nearType = POINT;
-            nearLabel = 0;
-        }
-        else if (Foam::mag(alpha-1) <= tol)
-        {
-            nearType = POINT;
-            nearLabel = 1;
-        }
-        else if ((alpha >= 0) && (alpha <= 1))
-        {
-            nearType = EDGE;
-            nearLabel = 0;
-        }
-    }
-
-    return hit;
+    return nearestPointClassify(p, nearType, nearLabel).hit();
 }
 
 
-
-
-
 // * * * * * * * * * * * * * * * Ostream Operator  * * * * * * * * * * * * * //
 
 template<class point, class pointRef>
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/directMappedVelocityFluxFixedValue/directMappedVelocityFluxFixedValueFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/derived/directMappedVelocityFluxFixedValue/directMappedVelocityFluxFixedValueFvPatchField.C
index 716c5b1235..17fa6f71fb 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/directMappedVelocityFluxFixedValue/directMappedVelocityFluxFixedValueFvPatchField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/directMappedVelocityFluxFixedValue/directMappedVelocityFluxFixedValueFvPatchField.C
@@ -112,13 +112,6 @@ directMappedVelocityFluxFixedValueFvPatchField
             << " in file " << dimensionedInternalField().objectPath()
             << exit(FatalError);
     }
-
-    // Force calculation of schedule (uses parallel comms)
-    const directMappedPolyPatch& mpp = refCast<const directMappedPolyPatch>
-    (
-        this->patch().patch()
-    );
-    (void)mpp.map().schedule();
 }
 
 
diff --git a/src/finiteVolume/fvMesh/singleCellFvMesh/singleCellFvMesh.C b/src/finiteVolume/fvMesh/singleCellFvMesh/singleCellFvMesh.C
index 4955f49007..8eaf3559f5 100644
--- a/src/finiteVolume/fvMesh/singleCellFvMesh/singleCellFvMesh.C
+++ b/src/finiteVolume/fvMesh/singleCellFvMesh/singleCellFvMesh.C
@@ -42,24 +42,26 @@ void Foam::singleCellFvMesh::agglomerateMesh
     const polyBoundaryMesh& oldPatches = mesh.boundaryMesh();
 
     // Check agglomeration within patch face range and continuous
-    labelList nAgglom(oldPatches.size());
+    labelList nAgglom(oldPatches.size(), 0);
 
     forAll(oldPatches, patchI)
     {
         const polyPatch& pp = oldPatches[patchI];
-
-        nAgglom[patchI] = max(agglom[patchI])+1;
-
-        forAll(pp, i)
+        if (pp.size() > 0)
         {
-            if (agglom[patchI][i] < 0  || agglom[patchI][i] >= pp.size())
+            nAgglom[patchI] = max(agglom[patchI])+1;
+
+            forAll(pp, i)
             {
-                FatalErrorIn
-                (
-                    "singleCellFvMesh::agglomerateMesh(..)"
-                )   << "agglomeration on patch " << patchI
-                    << " is out of range 0.." << pp.size()-1
-                    << exit(FatalError);
+                if (agglom[patchI][i] < 0  || agglom[patchI][i] >= pp.size())
+                {
+                    FatalErrorIn
+                    (
+                        "singleCellFvMesh::agglomerateMesh(..)"
+                    )   << "agglomeration on patch " << patchI
+                        << " is out of range 0.." << pp.size()-1
+                        << exit(FatalError);
+                }
             }
         }
     }
@@ -155,6 +157,8 @@ void Foam::singleCellFvMesh::agglomerateMesh
 
     forAll(oldPatches, patchI)
     {
+        patchStarts[patchI] = coarseI;
+
         const polyPatch& pp = oldPatches[patchI];
 
         if (pp.size() > 0)
@@ -170,8 +174,6 @@ void Foam::singleCellFvMesh::agglomerateMesh
             // From agglomeration to compact patch face
             labelList agglomToFace(nAgglom[patchI], -1);
 
-            patchStarts[patchI] = coarseI;
-
             forAll(pp, i)
             {
                 label myAgglom = agglom[patchI][i];
@@ -223,9 +225,9 @@ void Foam::singleCellFvMesh::agglomerateMesh
                     );
                 }
             }
-
-            patchSizes[patchI] = coarseI-patchStarts[patchI];
         }
+
+        patchSizes[patchI] = coarseI-patchStarts[patchI];
     }
 
     //Pout<< "patchStarts:" << patchStarts << endl;
diff --git a/src/meshTools/Make/files b/src/meshTools/Make/files
index 14e561fcba..2c9f81195d 100644
--- a/src/meshTools/Make/files
+++ b/src/meshTools/Make/files
@@ -44,7 +44,6 @@ octree/octreeDataFace.C
 octree/treeBoundBox.C
 octree/treeNodeName.C
 octree/treeLeafName.C
-octree/pointIndexHitIOList.C
 
 indexedOctree/indexedOctreeName.C
 indexedOctree/treeDataCell.C
diff --git a/src/meshTools/indexedOctree/indexedOctree.C b/src/meshTools/indexedOctree/indexedOctree.C
index 3f85828c09..9971bd97b4 100644
--- a/src/meshTools/indexedOctree/indexedOctree.C
+++ b/src/meshTools/indexedOctree/indexedOctree.C
@@ -2794,6 +2794,30 @@ Foam::indexedOctree<Type>::getVolumeType
 }
 
 
+template <class Type>
+template <class CompareOp>
+void Foam::indexedOctree<Type>::findNear
+(
+    const scalar nearDist,
+    const indexedOctree<Type>& tree2,
+    CompareOp& cop
+) const
+{
+    findNear
+    (
+        nearDist,
+        true,
+        *this,
+        nodePlusOctant(0, 0),
+        bb(),
+        tree2,
+        nodePlusOctant(0, 0),
+        tree2.bb(),
+        cop
+    );
+}
+
+
 // Print contents of nodeI
 template <class Type>
 void Foam::indexedOctree<Type>::print
diff --git a/src/meshTools/indexedOctree/treeDataTriSurface.C b/src/meshTools/indexedOctree/treeDataTriSurface.C
index e0f6d1a9df..f380167174 100644
--- a/src/meshTools/indexedOctree/treeDataTriSurface.C
+++ b/src/meshTools/indexedOctree/treeDataTriSurface.C
@@ -35,145 +35,145 @@ defineTypeNameAndDebug(Foam::treeDataTriSurface, 0);
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-// Fast distance to triangle calculation. From
-// "Distance Between Point and Trangle in 3D"
-// David Eberly, Magic Software Inc. Aug. 2003.
-// Works on function Q giving distance to point and tries to minimize this.
-Foam::scalar Foam::treeDataTriSurface::nearestCoords
-(
-    const point& base,
-    const point& E0,
-    const point& E1,
-    const scalar a,
-    const scalar b,
-    const scalar c,
-    const point& P,
-    scalar& s,
-    scalar& t
-)
-{
-    // distance vector
-    const vector D(base - P);
+// // Fast distance to triangle calculation. From
+// // "Distance Between Point and Triangle in 3D"
+// // David Eberly, Magic Software Inc. Aug. 2003.
+// // Works on function Q giving distance to point and tries to minimize this.
+// Foam::scalar Foam::treeDataTriSurface::nearestCoords
+// (
+//     const point& base,
+//     const point& E0,
+//     const point& E1,
+//     const scalar a,
+//     const scalar b,
+//     const scalar c,
+//     const point& P,
+//     scalar& s,
+//     scalar& t
+// )
+// {
+//     // distance vector
+//     const vector D(base - P);
 
-    // Precalculate distance factors.
-    const scalar d = E0 & D;
-    const scalar e = E1 & D;
+//     // Precalculate distance factors.
+//     const scalar d = E0 & D;
+//     const scalar e = E1 & D;
 
-    // Do classification
-    const scalar det = a*c - b*b;
+//     // Do classification
+//     const scalar det = a*c - b*b;
 
-    s = b*e - c*d;
-    t = b*d - a*e;
+//     s = b*e - c*d;
+//     t = b*d - a*e;
 
-    if (s+t < det)
-    {
-        if (s < 0)
-        {
-            if (t < 0)
-            {
-                //region 4
-                if (e > 0)
-                {
-                    //min on edge t = 0
-                    t = 0;
-                    s = (d >= 0 ? 0 : (-d >= a ? 1 : -d/a));
-                }
-                else
-                {
-                    //min on edge s=0
-                    s = 0;
-                    t = (e >= 0 ? 0 : (-e >= c ? 1 : -e/c));
-                }
-            }
-            else
-            {
-                //region 3. Min on edge s = 0
-                s = 0;
-                t = (e >= 0 ? 0 : (-e >= c ? 1 : -e/c));
-            }
-        }
-        else if (t < 0)
-        {
-            //region 5
-            t = 0;
-            s = (d >= 0 ? 0 : (-d >= a ? 1 : -d/a));
-        }
-        else
-        {
-            //region 0
-            const scalar invDet = 1/det;
-            s *= invDet;
-            t *= invDet;
-        }
-    }
-    else
-    {
-        if (s < 0)
-        {
-            //region 2
-            const scalar tmp0 = b + d;
-            const scalar tmp1 = c + e;
-            if (tmp1 > tmp0)
-            {
-                //min on edge s+t=1
-                const scalar numer = tmp1 - tmp0;
-                const scalar denom = a-2*b+c;
-                s = (numer >= denom ? 1 : numer/denom);
-                t = 1 - s;
-            }
-            else
-            {
-                //min on edge s=0
-                s = 0;
-                t = (tmp1 <= 0 ? 1 : (e >= 0 ? 0 : - e/c));
-            }
-        }
-        else if (t < 0)
-        {
-            //region 6
-            const scalar tmp0 = b + d;
-            const scalar tmp1 = c + e;
-            if (tmp1 > tmp0)
-            {
-                //min on edge s+t=1
-                const scalar numer = tmp1 - tmp0;
-                const scalar denom = a-2*b+c;
-                s = (numer >= denom ? 1 : numer/denom);
-                t = 1 - s;
-            }
-            else
-            {
-                //min on edge t=0
-                t = 0;
-                s = (tmp1 <= 0 ? 1 : (d >= 0 ? 0 : - d/a));
-            }
-        }
-        else
-        {
-            //region 1
-            const scalar numer = c+e-(b+d);
-            if (numer <= 0)
-            {
-                s = 0;
-            }
-            else
-            {
-                const scalar denom = a-2*b+c;
-                s = (numer >= denom ? 1 : numer/denom);
-            }
-        }
-        t = 1 - s;
-    }
+//     if (s+t < det)
+//     {
+//         if (s < 0)
+//         {
+//             if (t < 0)
+//             {
+//                 //region 4
+//                 if (e > 0)
+//                 {
+//                     //min on edge t = 0
+//                     t = 0;
+//                     s = (d >= 0 ? 0 : (-d >= a ? 1 : -d/a));
+//                 }
+//                 else
+//                 {
+//                     //min on edge s=0
+//                     s = 0;
+//                     t = (e >= 0 ? 0 : (-e >= c ? 1 : -e/c));
+//                 }
+//             }
+//             else
+//             {
+//                 //region 3. Min on edge s = 0
+//                 s = 0;
+//                 t = (e >= 0 ? 0 : (-e >= c ? 1 : -e/c));
+//             }
+//         }
+//         else if (t < 0)
+//         {
+//             //region 5
+//             t = 0;
+//             s = (d >= 0 ? 0 : (-d >= a ? 1 : -d/a));
+//         }
+//         else
+//         {
+//             //region 0
+//             const scalar invDet = 1/det;
+//             s *= invDet;
+//             t *= invDet;
+//         }
+//     }
+//     else
+//     {
+//         if (s < 0)
+//         {
+//             //region 2
+//             const scalar tmp0 = b + d;
+//             const scalar tmp1 = c + e;
+//             if (tmp1 > tmp0)
+//             {
+//                 //min on edge s+t=1
+//                 const scalar numer = tmp1 - tmp0;
+//                 const scalar denom = a-2*b+c;
+//                 s = (numer >= denom ? 1 : numer/denom);
+//                 t = 1 - s;
+//             }
+//             else
+//             {
+//                 //min on edge s=0
+//                 s = 0;
+//                 t = (tmp1 <= 0 ? 1 : (e >= 0 ? 0 : - e/c));
+//             }
+//         }
+//         else if (t < 0)
+//         {
+//             //region 6
+//             const scalar tmp0 = b + d;
+//             const scalar tmp1 = c + e;
+//             if (tmp1 > tmp0)
+//             {
+//                 //min on edge s+t=1
+//                 const scalar numer = tmp1 - tmp0;
+//                 const scalar denom = a-2*b+c;
+//                 s = (numer >= denom ? 1 : numer/denom);
+//                 t = 1 - s;
+//             }
+//             else
+//             {
+//                 //min on edge t=0
+//                 t = 0;
+//                 s = (tmp1 <= 0 ? 1 : (d >= 0 ? 0 : - d/a));
+//             }
+//         }
+//         else
+//         {
+//             //region 1
+//             const scalar numer = c+e-(b+d);
+//             if (numer <= 0)
+//             {
+//                 s = 0;
+//             }
+//             else
+//             {
+//                 const scalar denom = a-2*b+c;
+//                 s = (numer >= denom ? 1 : numer/denom);
+//             }
+//         }
+//         t = 1 - s;
+//     }
 
 
-    // Calculate distance.
-    // Note: abs should not be needed but truncation error causes problems
-    // with points very close to one of the triangle vertices.
-    // (seen up to -9e-15). Alternatively add some small value.
+//     // Calculate distance.
+//     // Note: abs should not be needed but truncation error causes problems
+//     // with points very close to one of the triangle vertices.
+//     // (seen up to -9e-15). Alternatively add some small value.
 
-    const scalar f = D & D;
-    return Foam::mag(a*s*s + 2*b*s*t + c*t*t + 2*d*s + 2*e*t + f);
-}
+//     const scalar f = D & D;
+//     return Foam::mag(a*s*s + 2*b*s*t + c*t*t + 2*d*s + 2*e*t + f);
+// }
 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
@@ -234,9 +234,7 @@ Foam::label Foam::treeDataTriSurface::getVolumeType
     (
         surface_,
         sample,
-        pHit.index(),
-        pHit.hitPoint(),
-        indexedOctree<treeDataTriSurface>::perturbTol()
+        pHit.index()
     );
 
     if (t == triSurfaceTools::UNKNOWN)
@@ -353,39 +351,43 @@ void Foam::treeDataTriSurface::findNearest
         //    )
         //)
         {
-            // Get spanning vectors of triangle
-            vector base(p1);
-            vector E0(p0 - p1);
-            vector E1(p2 - p1);
+            // // Get spanning vectors of triangle
+            // vector base(p1);
+            // vector E0(p0 - p1);
+            // vector E1(p2 - p1);
 
-            scalar a(E0& E0);
-            scalar b(E0& E1);
-            scalar c(E1& E1);
+            // scalar a(E0& E0);
+            // scalar b(E0& E1);
+            // scalar c(E1& E1);
 
-            // Get nearest point in s,t coordinates (s is along E0, t is along
-            // E1)
-            scalar s;
-            scalar t;
+            // // Get nearest point in s,t coordinates (s is along E0, t
+            // // is along E1)
+            // scalar s;
+            // scalar t;
 
-            scalar distSqr = nearestCoords
-            (
-                base,
-                E0,
-                E1,
-                a,
-                b,
-                c,
-                sample,
+            // scalar distSqr = nearestCoords
+            // (
+            //     base,
+            //     E0,
+            //     E1,
+            //     a,
+            //     b,
+            //     c,
+            //     sample,
 
-                s,
-                t
-            );
+            //     s,
+            //     t
+            // );
+
+            pointHit pHit = triPointRef(p0, p1, p2).nearestPoint(sample);
+
+            scalar distSqr = sqr(pHit.distance());
 
             if (distSqr < nearestDistSqr)
             {
                 nearestDistSqr = distSqr;
                 minIndex = index;
-                nearestPoint = base + s*E0 + t*E1;
+                nearestPoint = pHit.rawPoint();
             }
         }
     }
diff --git a/src/meshTools/indexedOctree/treeDataTriSurface.H b/src/meshTools/indexedOctree/treeDataTriSurface.H
index f602783afe..dc959f78ff 100644
--- a/src/meshTools/indexedOctree/treeDataTriSurface.H
+++ b/src/meshTools/indexedOctree/treeDataTriSurface.H
@@ -60,20 +60,20 @@ class treeDataTriSurface
 
     // Private Member Functions
 
-        //- fast triangle nearest point calculation. Returns point in E0, E1
-        //  coordinate system:  base + s*E0 + t*E1
-        static scalar nearestCoords
-        (
-            const point& base,
-            const point& E0,
-            const point& E1,
-            const scalar a,
-            const scalar b,
-            const scalar c,
-            const point& P,
-            scalar& s,
-            scalar& t
-        );
+        // //- fast triangle nearest point calculation. Returns point in E0, E1
+        // //  coordinate system:  base + s*E0 + t*E1
+        // static scalar nearestCoords
+        // (
+        //     const point& base,
+        //     const point& E0,
+        //     const point& E1,
+        //     const scalar a,
+        //     const scalar b,
+        //     const scalar c,
+        //     const point& P,
+        //     scalar& s,
+        //     scalar& t
+        // );
 
 public:
 
diff --git a/src/meshTools/triSurface/booleanOps/surfaceIntersection/edgeIntersections.C b/src/meshTools/triSurface/booleanOps/surfaceIntersection/edgeIntersections.C
index 110f9dcd66..c8955f7470 100644
--- a/src/meshTools/triSurface/booleanOps/surfaceIntersection/edgeIntersections.C
+++ b/src/meshTools/triSurface/booleanOps/surfaceIntersection/edgeIntersections.C
@@ -430,10 +430,7 @@ bool Foam::edgeIntersections::offsetPerturb
 
         point ctr = tri.centre();
 
-        // Get measure for tolerance.
-        scalar tolDim = 0.001*mag(tri.a() - ctr);
-
-        tri.classify(pHit.hitPoint(), tolDim, nearType, nearLabel);
+        tri.classify(pHit.hitPoint(), nearType, nearLabel);
 
         if (nearType == triPointRef::POINT || nearType == triPointRef::EDGE)
         {
diff --git a/src/meshTools/triSurface/booleanOps/surfaceIntersection/surfaceIntersection.C b/src/meshTools/triSurface/booleanOps/surfaceIntersection/surfaceIntersection.C
index 6602fbcf43..58c206aa75 100644
--- a/src/meshTools/triSurface/booleanOps/surfaceIntersection/surfaceIntersection.C
+++ b/src/meshTools/triSurface/booleanOps/surfaceIntersection/surfaceIntersection.C
@@ -315,7 +315,7 @@ void Foam::surfaceIntersection::classifyHit
         surf2Pts[f2[0]],
         surf2Pts[f2[1]],
         surf2Pts[f2[2]]
-    ).classify(pHit.hitPoint(), tolDim, nearType, nearLabel);
+    ).classify(pHit.hitPoint(), nearType, nearLabel);
 
     // Classify points on edge of surface1
     label edgeEnd =
diff --git a/src/meshTools/triSurface/octreeData/octreeDataTriSurface.C b/src/meshTools/triSurface/octreeData/octreeDataTriSurface.C
index c80991819e..55e60a3111 100644
--- a/src/meshTools/triSurface/octreeData/octreeDataTriSurface.C
+++ b/src/meshTools/triSurface/octreeData/octreeDataTriSurface.C
@@ -43,7 +43,7 @@ Foam::scalar Foam::octreeDataTriSurface::tol(1E-6);
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 // Fast distance to triangle calculation. From
-// "Distance Between Point and Trangle in 3D"
+// "Distance Between Point and Triangle in 3D"
 // David Eberly, Magic Software Inc. Aug. 2003.
 // Works on function Q giving distance to point and tries to minimize this.
 void Foam::octreeDataTriSurface::nearestCoords
diff --git a/src/meshTools/triSurface/orientedSurface/orientedSurface.C b/src/meshTools/triSurface/orientedSurface/orientedSurface.C
index 6462df4cb3..38a7afd4b8 100644
--- a/src/meshTools/triSurface/orientedSurface/orientedSurface.C
+++ b/src/meshTools/triSurface/orientedSurface/orientedSurface.C
@@ -234,9 +234,7 @@ void Foam::orientedSurface::propagateOrientation
     (
         s,
         samplePoint,
-        nearestFaceI,
-        nearestPt,
-        10*SMALL
+        nearestFaceI
     );
 
     if (side == triSurfaceTools::UNKNOWN)
diff --git a/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C b/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C
index 5a0cc03623..af30b3bcc8 100644
--- a/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C
+++ b/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C
@@ -2121,12 +2121,13 @@ Foam::vector Foam::triSurfaceTools::surfaceNormal
 
     label nearType;
     label nearLabel;
+
     triPointRef
     (
         points[f[0]],
         points[f[1]],
         points[f[2]]
-    ).classify(nearestPt, 1E-6, nearType, nearLabel);
+    ).classify(nearestPt, nearType, nearLabel);
 
     if (nearType == triPointRef::NONE)
     {
@@ -2199,28 +2200,61 @@ Foam::triSurfaceTools::sideType Foam::triSurfaceTools::surfaceSide
 (
     const triSurface& surf,
     const point& sample,
-    const label nearestFaceI,   // nearest face
-    const point& nearestPoint,  // nearest point on nearest face
-    const scalar tol
+    const label nearestFaceI
 )
 {
     const labelledTri& f = surf[nearestFaceI];
     const pointField& points = surf.points();
 
-    // Find where point is on triangle. Note tolerance needed. Is relative
-    // one (used in comparing normalized [0..1] triangle coordinates).
+    // Find where point is on triangle.
     label nearType, nearLabel;
-    triPointRef
+
+    pointHit pHit = triPointRef
     (
         points[f[0]],
         points[f[1]],
         points[f[2]]
-    ).classify(nearestPoint, tol, nearType, nearLabel);
+    ).nearestPointClassify(sample, nearType, nearLabel);
+
+    const point& nearestPoint(pHit.rawPoint());
 
     if (nearType == triPointRef::NONE)
     {
+        vector sampleNearestVec = (sample - nearestPoint);
+
         // Nearest to face interior. Use faceNormal to determine side
-        scalar c = (sample - nearestPoint) & surf.faceNormals()[nearestFaceI];
+        scalar c = sampleNearestVec & surf.faceNormals()[nearestFaceI];
+
+        // // If the sample is essentially on the face, do not check for
+        // // it being perpendicular.
+
+        // scalar magSampleNearestVec = mag(sampleNearestVec);
+
+        // if (magSampleNearestVec > SMALL)
+        // {
+        //     c /= magSampleNearestVec*mag(surf.faceNormals()[nearestFaceI]);
+
+        //     if (mag(c) < 0.99)
+        //     {
+        //         FatalErrorIn("triSurfaceTools::surfaceSide")
+        //             << "nearestPoint identified as being on triangle face "
+        //             << "but vector from nearestPoint to sample is not "
+        //             << "perpendicular to the normal." << nl
+        //             << "sample: " << sample << nl
+        //             << "nearestPoint: " << nearestPoint << nl
+        //             << "sample - nearestPoint: "
+        //             << sample - nearestPoint << nl
+        //             << "normal: " << surf.faceNormals()[nearestFaceI] << nl
+        //             << "mag(sample - nearestPoint): "
+        //             << mag(sample - nearestPoint) << nl
+        //             << "normalised dot product: " << c << nl
+        //             << "triangle vertices: " << nl
+        //             << "    " << points[f[0]] << nl
+        //             << "    " << points[f[1]] << nl
+        //             << "    " << points[f[2]] << nl
+        //             << abort(FatalError);
+        //     }
+        // }
 
         if (c > 0)
         {
@@ -2239,13 +2273,13 @@ Foam::triSurfaceTools::sideType Foam::triSurfaceTools::surfaceSide
         // Get the edge. Assume order of faceEdges same as face vertices.
         label edgeI = surf.faceEdges()[nearestFaceI][nearLabel];
 
-        //if (debug)
-        //{
+        // if (debug)
+        // {
         //    // Check order of faceEdges same as face vertices.
         //    const edge& e = surf.edges()[edgeI];
         //    const labelList& meshPoints = surf.meshPoints();
         //    const edge meshEdge(meshPoints[e[0]], meshPoints[e[1]]);
-        //
+
         //    if
         //    (
         //        meshEdge
@@ -2259,7 +2293,7 @@ Foam::triSurfaceTools::sideType Foam::triSurfaceTools::surfaceSide
         //            << " in face " << f
         //            << abort(FatalError);
         //    }
-        //}
+        // }
 
         return edgeSide(surf, sample, nearestPoint, edgeI);
     }
@@ -2717,7 +2751,14 @@ void Foam::triSurfaceTools::calcInterpolationWeights
 
             triPointRef tri(f.tri(points));
 
-            pointHit nearest = tri.nearestPoint(samplePt);
+            label nearType, nearLabel;
+
+            pointHit nearest = tri.nearestPointClassify
+            (
+                samplePt,
+                nearType,
+                nearLabel
+            );
 
             if (nearest.hit())
             {
@@ -2741,14 +2782,6 @@ void Foam::triSurfaceTools::calcInterpolationWeights
                 minDistance = nearest.distance();
 
                 // Outside triangle. Store nearest.
-                label nearType, nearLabel;
-                tri.classify
-                (
-                    nearest.rawPoint(),
-                    1E-6,               // relative tolerance
-                    nearType,
-                    nearLabel
-                );
 
                 if (nearType == triPointRef::POINT)
                 {
@@ -2779,12 +2812,12 @@ void Foam::triSurfaceTools::calcInterpolationWeights
                         max
                         (
                             0,
-                            mag(nearest.rawPoint()-p0)/mag(p1-p0)
+                            mag(nearest.rawPoint() - p0)/mag(p1 - p0)
                         )
                     );
 
                     // Interpolate
-                    weights[0] = 1-s;
+                    weights[0] = 1 - s;
                     weights[1] = s;
                     weights[2] = -GREAT;
 
@@ -2830,7 +2863,6 @@ Foam::surfaceLocation Foam::triSurfaceTools::classify
     triPointRef(s[triI].tri(s.points())).classify
     (
         trianglePoint,
-        1E-6,
         elemType,
         index
     );
diff --git a/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.H b/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.H
index 62d66d031f..744a12217d 100644
--- a/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.H
+++ b/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.H
@@ -458,9 +458,7 @@ public:
         (
             const triSurface& surf,
             const point& sample,
-            const label nearestFaceI,   // nearest face
-            const point& nearestPt,     // nearest point on nearest face
-            const scalar tol            // tolerance for nearness test.
+            const label nearestFaceI
         );
 
     // Triangulation of faces
diff --git a/src/postProcessing/functionObjects/field/streamLine/streamLineParticle.C b/src/postProcessing/functionObjects/field/streamLine/streamLineParticle.C
index ce90b7a55d..56b05c88ca 100644
--- a/src/postProcessing/functionObjects/field/streamLine/streamLineParticle.C
+++ b/src/postProcessing/functionObjects/field/streamLine/streamLineParticle.C
@@ -168,7 +168,6 @@ bool Foam::streamLineParticle::move(streamLineParticle::trackData& td)
         td.keepParticle
     && !td.switchProcessor
     && lifeTime_ > 0
-    && tEnd > ROOTVSMALL
     )
     {
         // TBD: implement subcycling so step through cells in more than
@@ -191,6 +190,12 @@ bool Foam::streamLineParticle::move(streamLineParticle::trackData& td)
 
         tEnd -= dt;
         stepFraction() = 1.0 - tEnd/deltaT;
+
+        if (tEnd <= ROOTVSMALL)
+        {
+            // Force removal
+            lifeTime_ = 0;
+        }
     }
 
     if (!td.keepParticle || lifeTime_ == 0)
diff --git a/src/sampling/sampledSurface/distanceSurface/distanceSurface.C b/src/sampling/sampledSurface/distanceSurface/distanceSurface.C
index b777edbdb3..4a4c9219d6 100644
--- a/src/sampling/sampledSurface/distanceSurface/distanceSurface.C
+++ b/src/sampling/sampledSurface/distanceSurface/distanceSurface.C
@@ -89,20 +89,29 @@ void Foam::distanceSurface::createGeometry()
 
         if (signed_)
         {
-            vectorField normal;
-            surfPtr_().getNormal(nearest, normal);
+            List<searchableSurface::volumeType> volType;
 
-            forAll(nearest, i)
+            surfPtr_().getVolumeType(cc, volType);
+
+            forAll(volType, i)
             {
-                vector d(cc[i]-nearest[i].hitPoint());
+                searchableSurface::volumeType vT = volType[i];
 
-                if ((d&normal[i]) > 0)
+                if (vT == searchableSurface::OUTSIDE)
                 {
-                    fld[i] = Foam::mag(d);
+                    fld[i] = Foam::mag(cc[i] - nearest[i].hitPoint());
+                }
+                else if (vT == searchableSurface::INSIDE)
+                {
+                    fld[i] = -Foam::mag(cc[i] - nearest[i].hitPoint());
                 }
                 else
                 {
-                    fld[i] = -Foam::mag(d);
+                    FatalErrorIn
+                    (
+                        "void Foam::distanceSurface::createGeometry()"
+                    )   << "getVolumeType failure, neither INSIDE or OUTSIDE"
+                        << exit(FatalError);
                 }
             }
         }
@@ -132,20 +141,30 @@ void Foam::distanceSurface::createGeometry()
 
             if (signed_)
             {
-                vectorField normal;
-                surfPtr_().getNormal(nearest, normal);
+                List<searchableSurface::volumeType> volType;
 
-                forAll(nearest, i)
+                surfPtr_().getVolumeType(cc, volType);
+
+                forAll(volType, i)
                 {
-                    vector d(cc[i]-nearest[i].hitPoint());
+                    searchableSurface::volumeType vT = volType[i];
 
-                    if ((d&normal[i]) > 0)
+                    if (vT == searchableSurface::OUTSIDE)
                     {
-                        fld[i] = Foam::mag(d);
+                        fld[i] = Foam::mag(cc[i] - nearest[i].hitPoint());
+                    }
+                    else if (vT == searchableSurface::INSIDE)
+                    {
+                        fld[i] = -Foam::mag(cc[i] - nearest[i].hitPoint());
                     }
                     else
                     {
-                        fld[i] = -Foam::mag(d);
+                        FatalErrorIn
+                        (
+                            "void Foam::distanceSurface::createGeometry()"
+                        )   << "getVolumeType failure, "
+                            << "neither INSIDE or OUTSIDE"
+                            << exit(FatalError);
                     }
                 }
             }
@@ -179,20 +198,31 @@ void Foam::distanceSurface::createGeometry()
 
         if (signed_)
         {
-            vectorField normal;
-            surfPtr_().getNormal(nearest, normal);
+            List<searchableSurface::volumeType> volType;
 
-            forAll(nearest, i)
+            surfPtr_().getVolumeType(pts, volType);
+
+            forAll(volType, i)
             {
-                vector d(pts[i]-nearest[i].hitPoint());
+                searchableSurface::volumeType vT = volType[i];
 
-                if ((d&normal[i]) > 0)
+                if (vT == searchableSurface::OUTSIDE)
                 {
-                    pointDistance_[i] = Foam::mag(d);
+                    pointDistance_[i] =
+                        Foam::mag(pts[i] - nearest[i].hitPoint());
+                }
+                else if (vT == searchableSurface::INSIDE)
+                {
+                    pointDistance_[i] =
+                        -Foam::mag(pts[i] - nearest[i].hitPoint());
                 }
                 else
                 {
-                    pointDistance_[i] = -Foam::mag(d);
+                    FatalErrorIn
+                    (
+                        "void Foam::distanceSurface::createGeometry()"
+                    )   << "getVolumeType failure, neither INSIDE or OUTSIDE"
+                        << exit(FatalError);
                 }
             }
         }
diff --git a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/kappatWallFunctions/kappatJayatillekeWallFunction/kappatJayatillekeWallFunctionFvPatchScalarField.C b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/kappatWallFunctions/kappatJayatillekeWallFunction/kappatJayatillekeWallFunctionFvPatchScalarField.C
index 7aa9fa2fe9..6efef7364b 100644
--- a/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/kappatWallFunctions/kappatJayatillekeWallFunction/kappatJayatillekeWallFunctionFvPatchScalarField.C
+++ b/src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/kappatWallFunctions/kappatJayatillekeWallFunction/kappatJayatillekeWallFunctionFvPatchScalarField.C
@@ -231,19 +231,17 @@ void kappatJayatillekeWallFunctionFvPatchScalarField::updateCoeffs()
         scalar P = Psmooth(Prat);
         scalar yPlusTherm = this->yPlusTherm(P, Prat);
 
-        // Evaluate new effective thermal diffusivity
-        scalar kappaEff = 0.0;
-        if (yPlus < yPlusTherm)
+        // Update turbulent thermal conductivity
+        if (yPlus > yPlusTherm)
         {
-            kappaEff = Pr*yPlus;
+            scalar nu = nuw[faceI];
+            scalar kt = nu*(yPlus/(Prt_/kappa_*log(E_*yPlusTherm) + P) - 1/Pr);
+            kappatw[faceI] = max(0.0, kt);
         }
         else
         {
-            kappaEff = nuw[faceI]*yPlus/(Prt_/kappa_*log(E_*yPlusTherm) + P);
+            kappatw[faceI] = 0.0;
         }
-
-        // Update turbulent thermal diffusivity
-        kappatw[faceI] = max(0.0, kappaEff - nuw[faceI]/Pr);
     }
 
     fixedValueFvPatchField<scalar>::updateCoeffs();
diff --git a/tutorials/heatTransfer/buoyantBoussinesqPimpleFoam/hotRoom/0/T.org b/tutorials/heatTransfer/buoyantBoussinesqPimpleFoam/hotRoom/0/T.org
index 98cc191af7..77b58a1a00 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 b247cdcbaf..ed297deb8e 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 0000000000..270f014558
--- /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 25243dd42f..fcc6e9c9cd 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 cb0be1e6e0..800bb57319 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 fef36c6b14..764772f38e 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 4c0b29ac70..dccd8e082d 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 b247cdcbaf..ed297deb8e 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 25243dd42f..fcc6e9c9cd 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 cb0be1e6e0..800bb57319 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 dc6fb20830..8206700137 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 68da3948e5..c51949b805 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 5486bbf9fc..3f7a67aa5e 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 e3e3c88c20..2876f3adad 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 73cfe8ada4..bf11b4f154 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 8df126a483..edea946f9a 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 9854db7add..d6ab24211f 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 0000000000..371d93c38a
--- /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 8caf680281..b238c08e08 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 9ad9aaa4ce..0473bd2bf3 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 0961828a5b..9d710bdf49 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 8ca46687fb..599204b449 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 992508c411..c453b27d81 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 98cc191af7..77b58a1a00 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 98cc191af7..77b58a1a00 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 9854db7add..d6ab24211f 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 0000000000..371d93c38a
--- /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 f9d959e2dc..242675a288 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 328bd3ed59..ab84aca2c4 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 c9a423a063..aa535a28eb 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 d49a1b712b..928afaeb56 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 a1da943200..fb9be497bf 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 0000000000..87dd594134
--- /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 25a67762b3..80041f81dc 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 8d30fa0a37..2dfb6798b8 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 c490be6e8f..fa8d3d5339 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 3bee6a9b10..94f365296e 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 a1da943200..48fe8a009f 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 0000000000..87dd594134
--- /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 1cca62ded7..64496570c2 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 58fe62417e..c9900e695f 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 c219a76e03..daa2528426 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 6034a1d7bf..f4e27ad95f 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 eb1abac09b..b2509f9670 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 0000000000..638c8c7d7c
--- /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 c36ea6b803..6faab27ee5 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 90364fcbf2..7d3e57864a 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 16665b05a1..1e968f0aea 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 ec9f3455bd..590a47c55a 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 e53a523c9d..8f6ced4ac3 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 b773399f72..7c00e466f4 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 b773399f72..7c00e466f4 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 b773399f72..7c00e466f4 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 71114003ec..ff4ba10011 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            waveTransmissive;
-                gamma           1.4;
-                fieldInf        100000;
-                lInf            0.4;
-                value           uniform 100000;
+                type            fixedValue;
+                value           uniform 1e5;
+            }
+        }
+    }
+
+    p
+    {
+        internalField   uniform 1e5;
+
+        boundaryField
+        {
+            ".*"
+            {
+                type            calculated;
+                value           uniform 1e5;
+            }
+
+            maxX
+            {
+                type            calculated;
+                value           uniform 1e5;
             }
         }
     }
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/topAir/fvSchemes b/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/topAir/fvSchemes
index 90364fcbf2..7d3e57864a 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 28c1a63b7a..d2d44df4fc 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 5afdd3b3cc..5dba622286 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 6034a1d7bf..f4e27ad95f 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 eb1abac09b..b2509f9670 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 0000000000..638c8c7d7c
--- /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 1e2a7dd206..fc51ab1ad7 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 e53a523c9d..8f6ced4ac3 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 90364fcbf2..7d3e57864a 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 ad1f5492fb..ed5885da9a 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 e53a523c9d..8f6ced4ac3 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 b773399f72..7c00e466f4 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 e53a523c9d..8f6ced4ac3 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 b773399f72..18b2ebb7eb 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 e53a523c9d..8f6ced4ac3 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 b773399f72..18b2ebb7eb 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 71114003ec..ff4ba10011 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            waveTransmissive;
-                gamma           1.4;
-                fieldInf        100000;
-                lInf            0.4;
-                value           uniform 100000;
+                type            fixedValue;
+                value           uniform 1e5;
+            }
+        }
+    }
+
+    p
+    {
+        internalField   uniform 1e5;
+
+        boundaryField
+        {
+            ".*"
+            {
+                type            calculated;
+                value           uniform 1e5;
+            }
+
+            maxX
+            {
+                type            calculated;
+                value           uniform 1e5;
             }
         }
     }
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/topAir/decomposeParDict b/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/topAir/decomposeParDict
index e53a523c9d..8f6ced4ac3 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 90364fcbf2..7d3e57864a 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 28c1a63b7a..d2d44df4fc 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 62058ec7ac..8e47fbcddf 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 6034a1d7bf..f4e27ad95f 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 eb1abac09b..b2509f9670 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 0000000000..638c8c7d7c
--- /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 c36ea6b803..6faab27ee5 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 e53a523c9d..8f6ced4ac3 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 a06836c2b3..c136ffea31 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 842419ac72..eff94d930f 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 0c78e811c6..32e60a1dd1 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 e53a523c9d..8f6ced4ac3 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 e53a523c9d..8f6ced4ac3 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 e53a523c9d..8f6ced4ac3 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 e53a523c9d..8f6ced4ac3 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 cf0e2f131c..ff4ba10011 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 e53a523c9d..8f6ced4ac3 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 a06836c2b3..c136ffea31 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 842419ac72..6e31656bb3 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;