diff --git a/ReleaseNotes-dev b/ReleaseNotes-dev
index 0a6041c55b..d35d1b518e 100644
--- a/ReleaseNotes-dev
+++ b/ReleaseNotes-dev
@@ -39,13 +39,43 @@
 *** DSMC
 *** Dynamic Mesh
 *** Numerics
+*** *Updated* command line help, e.g. `snappyHexMesh -help' now gives:
+
+    Usage: snappyHexMesh [OPTIONS]
+    options:
+      -case <dir>       specify alternate case directory, default is the cwd
+      -overwrite        overwrite existing mesh/results files
+      -parallel         run in parallel
+      -srcDoc           display source code in browser
+      -doc              display application documentation in browser
+      -help             print the usage
+
+*** *New* Surface film library
+    + Creation of films by particle addition, or initial film distribution
+    + Coupled with the lagrangian/intermediate cloud hierarchy library
+    + Hierarchical design, consisting of
+      + kinematic film: mass, momentum
+        + constant thermodynamic properties
+      + thermodynamic film: mass, momentum and enthalpy
+        + constant, or temperature dependant thermodynamic properties
+    + Sub-models:
+      + detachment/dripping whereby particles (re)enter the originating cloud
+        + particle sizes set according to PDF
+        + other properties set to ensure mass, momentum and energy conservation
+      + heat transfer to/from walls and film surface
+      + film evaporation and boiling
+    + Additional wall functions for primary region momentum and temperature
+      taking film into account
+    + Parallel aware
+*** *Updated* particle tracking algorithm
 
 * Solvers
   A number of new solvers have been developed for a range of engineering
   applications.  There has been a set of improvements to certain classes of
   solver that are introduced in this release.
 *** *New* Solvers
-    + ...
+    + =reactingParcelFilmFoam=: Lagrangian cloud and film transport in a
+      reacting gas phase system
 *** Modifications to multiphase and buoyant solvers
     + ...
 *** Modifications to solvers for sensible enthalpy
@@ -58,23 +88,43 @@
 * Boundary conditions
   New boundary conditions have been introduced to support new applications in
   OpenFOAM.
-  + ...
+  + *New* wall functions:
+    + kappatJayatillekeWallFunction: incompressible RAS thermal wall function
 
 * Utilities
   There have been some utilities added and updated in this release.
 *** *New* utilities
-    + ...
+    + =extrudeToRegionMesh=: Extrude faceZones into separate mesh (as a
+      different region)
+      + used to e.g. extrude baffles (extrude internal faces) or create
+        liquid film regions
+      + if extruding internal faces:
+         + create baffles in original mesh with directMappedWall patches
+      + if extruding boundary faces:
+        + convert boundary faces to directMappedWall patches
+      + extrude edges of faceZone as a <zone>_sidePatch
+      + extrude edges inbetween different faceZones as a
+        (nonuniformTransform)cyclic <zoneA>_<zoneB>
+      + extrudes into master direction (i.e. away from the owner cell
+        if flipMap is false)
 *** Updated utilities
     + ...
 
 * Post-processing
   + =foamToEnsight=: new =-nodeValues= option to generate and output nodal
-    field data. This is optional and saves memory (compared to letting EnSight
-    do the interpolation) in case of meshes with large
-    amounts of polyhedral cells. On typical snappyHexMesh generated meshes
-    we've seen differences of 0.4 v.s. 2.5Gb memory usage.
+    field data.
+  + Function objects:
+    + =residualControl=: new function object to allow users to terminate steady
+    state calculations when the defined residual levels are achieved
+    + =abortCalculation=: watches for presence of the named file in the
+    $FOAM_CASE directory and aborts the calculation if it is present
+    + =timeActivatedFileUpdate=: performs a file copy/replacement once a
+    specified time has been reached, e.g. to automagically change fvSchemes and
+    fvSolution during a calculation
+    + =streamLine=: generate streamlines; ouputs both trajectory and field data
 
 * New tutorials
   There is a large number of new tutorials to support the new solvers in the
   release.
-  + ...
+  + =reactingParcelFilmFoam= tutorials:
+    + multipleBoxes, hotBoxes, panel, evaporationTest
diff --git a/applications/solvers/compressible/rhoPorousSimpleFoam/createFields.H b/applications/solvers/compressible/rhoPorousSimpleFoam/createFields.H
index d7a6b97ac1..360687a60d 100644
--- a/applications/solvers/compressible/rhoPorousSimpleFoam/createFields.H
+++ b/applications/solvers/compressible/rhoPorousSimpleFoam/createFields.H
@@ -45,14 +45,9 @@
     scalar pRefValue = 0.0;
     setRefCell(p, mesh.solutionDict().subDict("SIMPLE"), pRefCell, pRefValue);
 
-    dimensionedScalar rhoMax
+    dimensionedScalar pMin
     (
-        mesh.solutionDict().subDict("SIMPLE").lookup("rhoMax")
-    );
-
-    dimensionedScalar rhoMin
-    (
-        mesh.solutionDict().subDict("SIMPLE").lookup("rhoMin")
+        mesh.solutionDict().subDict("SIMPLE").lookup("pMin")
     );
 
     Info<< "Creating turbulence model\n" << endl;
@@ -70,16 +65,29 @@
     dimensionedScalar initialMass = fvc::domainIntegrate(rho);
 
     thermalPorousZones pZones(mesh);
+    Switch pressureImplicitPorosity(false);
 
     // nUCorrectors used for pressureImplicitPorosity
     int nUCorr = 0;
-    const bool pressureImplicitPorosity =
-    (
-        pZones.size()
-     && mesh.solutionDict().subDict("SIMPLE").readIfPresent
-        (
-            "nUCorrectors",
-            nUCorr
-        )
-     && (nUCorr > 0)
-    );
+    if (pZones.size())
+    {
+        // nUCorrectors for pressureImplicitPorosity
+        if (mesh.solutionDict().subDict("SIMPLE").found("nUCorrectors"))
+        {
+            nUCorr = readInt
+            (
+                mesh.solutionDict().subDict("SIMPLE").lookup("nUCorrectors")
+            );
+        }
+
+        if (nUCorr > 0)
+        {
+            pressureImplicitPorosity = true;
+            Info<< "Using pressure implicit porosity" << endl;
+        }
+        else
+        {
+            Info<< "Using pressure explicit porosity" << endl;
+        }
+    }
+
diff --git a/applications/solvers/compressible/rhoPorousSimpleFoam/hEqn.H b/applications/solvers/compressible/rhoPorousSimpleFoam/hEqn.H
index ff0b91bcb7..9e5138de8a 100644
--- a/applications/solvers/compressible/rhoPorousSimpleFoam/hEqn.H
+++ b/applications/solvers/compressible/rhoPorousSimpleFoam/hEqn.H
@@ -5,8 +5,8 @@
       - fvm::Sp(fvc::div(phi), h)
       - fvm::laplacian(turbulence->alphaEff(), h)
      ==
-        fvc::div(phi/fvc::interpolate(rho), rho/psi, "div(U,p)")
-      - (rho/psi)*fvc::div(phi/fvc::interpolate(rho))
+        fvc::div(phi/fvc::interpolate(rho)*fvc::interpolate(p, "div(U,p)"))
+      - p*fvc::div(phi/fvc::interpolate(rho))
     );
 
     pZones.addEnthalpySource(thermo, rho, hEqn);
diff --git a/applications/solvers/compressible/rhoPorousSimpleFoam/pEqn.H b/applications/solvers/compressible/rhoPorousSimpleFoam/pEqn.H
index 0386899612..40b7cb3969 100644
--- a/applications/solvers/compressible/rhoPorousSimpleFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoPorousSimpleFoam/pEqn.H
@@ -1,8 +1,3 @@
-rho = thermo.rho();
-rho = max(rho, rhoMin);
-rho = min(rho, rhoMax);
-rho.relax();
-
 if (pressureImplicitPorosity)
 {
     U = trTU()&UEqn().H();
@@ -63,6 +58,8 @@ else
 
 U.correctBoundaryConditions();
 
+bound(p, pMin);
+
 // For closed-volume cases adjust the pressure and density levels
 // to obey overall mass continuity
 if (closedVolume)
@@ -72,7 +69,5 @@ if (closedVolume)
 }
 
 rho = thermo.rho();
-rho = max(rho, rhoMin);
-rho = min(rho, rhoMax);
 rho.relax();
 Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl;
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/applications/utilities/postProcessing/sampling/sample/sample.C b/applications/utilities/postProcessing/sampling/sample/sample.C
index a68ab0b8fb..2b4f7ab008 100644
--- a/applications/utilities/postProcessing/sampling/sample/sample.C
+++ b/applications/utilities/postProcessing/sampling/sample/sample.C
@@ -116,7 +116,7 @@ int main(int argc, char *argv[])
         sampledSets::typeName,
         mesh,
         sampleDict,
-        IOobject::MUST_READ,
+        IOobject::MUST_READ_IF_MODIFIED,
         true
     );
 
@@ -125,7 +125,7 @@ int main(int argc, char *argv[])
         sampledSurfaces::typeName,
         mesh,
         sampleDict,
-        IOobject::MUST_READ,
+        IOobject::MUST_READ_IF_MODIFIED,
         true
     );
 
diff --git a/src/lagrangian/basic/Cloud/CloudIO.C b/src/lagrangian/basic/Cloud/CloudIO.C
index bc0e7defd2..87b23e4538 100644
--- a/src/lagrangian/basic/Cloud/CloudIO.C
+++ b/src/lagrangian/basic/Cloud/CloudIO.C
@@ -121,12 +121,16 @@ void Foam::Cloud<ParticleType>::initCloud(const bool checkClass)
     }
     else
     {
-        WarningIn("Cloud<ParticleType>::initCloud(const bool checkClass)")
-            << "Cannot read particle positions file " << nl
+        Pout<< "Cannot read particle positions file " << nl
             << "    " << ioP.path() << nl
             << "    assuming the initial cloud contains 0 particles." << endl;
     }
 
+    // Ask for the tetBasePtIs to trigger all processors to build
+    // them, otherwise, if some processors have no particles then
+    // there is a comms mismatch.
+    polyMesh_.tetBasePtIs();
+
     forAllIter(typename Cloud<ParticleType>, *this, pIter)
     {
         ParticleType& p = pIter();
diff --git a/src/thermophysicalModels/radiation/derivedFvPatchFields/MarshakRadiation/MarshakRadiationMixedFvPatchScalarField.C b/src/thermophysicalModels/radiation/derivedFvPatchFields/MarshakRadiation/MarshakRadiationMixedFvPatchScalarField.C
index 9ae7b435ca..9bc80d5c9f 100644
--- a/src/thermophysicalModels/radiation/derivedFvPatchFields/MarshakRadiation/MarshakRadiationMixedFvPatchScalarField.C
+++ b/src/thermophysicalModels/radiation/derivedFvPatchFields/MarshakRadiation/MarshakRadiationMixedFvPatchScalarField.C
@@ -77,24 +77,19 @@ Foam::MarshakRadiationFvPatchScalarField::MarshakRadiationFvPatchScalarField
 {
     if (dict.found("value"))
     {
-        fvPatchScalarField::operator=
-        (
-            scalarField("value", dict, p.size())
-        );
-        refValue() = scalarField("refValue", dict, p.size());
-        refGrad() = scalarField("refGradient", dict, p.size());
-        valueFraction() = scalarField("valueFraction", dict, p.size());
+        refValue() = scalarField("value", dict, p.size());
     }
     else
     {
-        const scalarField& Tp =
-            patch().lookupPatchField<volScalarField, scalar>(TName_);
-
-        refValue() = 4.0*constant::physicoChemical::sigma.value()*pow4(Tp);
-        refGrad() = 0.0;
-
-        fvPatchScalarField::operator=(refValue());
+        refValue() = 0.0;
     }
+
+    // zero gradient
+    refGrad() = 0.0;
+
+    valueFraction() = 1.0;
+
+    fvPatchScalarField::operator=(refValue());
 }
 
 
diff --git a/src/thermophysicalModels/radiation/derivedFvPatchFields/MarshakRadiationFixedT/MarshakRadiationFixedTMixedFvPatchScalarField.C b/src/thermophysicalModels/radiation/derivedFvPatchFields/MarshakRadiationFixedT/MarshakRadiationFixedTMixedFvPatchScalarField.C
index 99b5c68267..cb8cc7db97 100644
--- a/src/thermophysicalModels/radiation/derivedFvPatchFields/MarshakRadiationFixedT/MarshakRadiationFixedTMixedFvPatchScalarField.C
+++ b/src/thermophysicalModels/radiation/derivedFvPatchFields/MarshakRadiationFixedT/MarshakRadiationFixedTMixedFvPatchScalarField.C
@@ -78,24 +78,15 @@ MarshakRadiationFixedTMixedFvPatchScalarField
     Trad_("Trad", dict, p.size()),
     emissivity_(readScalar(dict.lookup("emissivity")))
 {
-    if (dict.found("value"))
-    {
-        fvPatchScalarField::operator=
-        (
-            scalarField("value", dict, p.size())
-        );
-        refValue() = scalarField("refValue", dict, p.size());
-        refGrad() = scalarField("refGradient", dict, p.size());
-        valueFraction() = scalarField("valueFraction", dict, p.size());
-    }
-    else
-    {
-        refValue() = 4.0*constant::physicoChemical::sigma.value()*pow4(Trad_);
-        refGrad() = 0.0;
-        valueFraction() = 1.0;
+    // refValue updated on each call to updateCoeffs()
+    refValue() = 4.0*constant::physicoChemical::sigma.value()*pow4(Trad_);
 
-        fvPatchScalarField::operator=(refValue());
-    }
+    // zero gradient
+    refGrad() = 0.0;
+
+    valueFraction() = 1.0;
+
+    fvPatchScalarField::operator=(refValue());
 }
 
 
diff --git a/src/thermophysicalModels/thermalPorousZone/Make/files b/src/thermophysicalModels/thermalPorousZone/Make/files
index 046355feaf..32e0e201f5 100644
--- a/src/thermophysicalModels/thermalPorousZone/Make/files
+++ b/src/thermophysicalModels/thermalPorousZone/Make/files
@@ -1,4 +1,9 @@
 thermalPorousZone/thermalPorousZone.C
 thermalPorousZone/thermalPorousZones.C
 
+thermalModel/thermalModel/thermalModel.C
+thermalModel/thermalModel/thermalModelNew.C
+thermalModel/fixedTemperature/fixedTemperature.C
+thermalModel/noThermalModel/noThermalModel.C
+
 LIB = $(FOAM_LIBBIN)/libthermalPorousZone
diff --git a/src/thermophysicalModels/thermalPorousZone/thermalModel/fixedTemperature/fixedTemperature.C b/src/thermophysicalModels/thermalPorousZone/thermalModel/fixedTemperature/fixedTemperature.C
new file mode 100644
index 0000000000..512746722c
--- /dev/null
+++ b/src/thermophysicalModels/thermalPorousZone/thermalModel/fixedTemperature/fixedTemperature.C
@@ -0,0 +1,98 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fixedTemperature.H"
+#include "addToRunTimeSelectionTable.H"
+#include "basicThermo.H"
+#include "volFields.H"
+#include "fvMatrices.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(fixedTemperature, 0);
+
+    addToRunTimeSelectionTable
+    (
+        thermalModel,
+        fixedTemperature,
+        pZone
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::fixedTemperature::fixedTemperature(const porousZone& pZone)
+:
+    thermalModel(pZone, typeName),
+    T_(readScalar(coeffDict_.lookup("T")))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor    * * * * * * * * * * * * * * //
+
+Foam::fixedTemperature::~fixedTemperature()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::fixedTemperature::addEnthalpySource
+(
+    const basicThermo& thermo,
+    const volScalarField& rho,
+    fvScalarMatrix& hEqn
+) const
+{
+    const labelList& zones = pZone_.zoneIds();
+    if (zones.empty() || T_ < 0.0)
+    {
+        return;
+    }
+
+    const fvMesh& mesh = pZone_.mesh();
+    const scalarField& V = mesh.V();
+    scalarField& hDiag = hEqn.diag();
+    scalarField& hSource = hEqn.source();
+
+    // TODO: generalize for non-fixedTemperature methods
+    const scalar rate = 1e6;
+
+    forAll(zones, zoneI)
+    {
+        const labelList& cells = mesh.cellZones()[zones[zoneI]];
+
+        forAll(cells, i)
+        {
+            hDiag[cells[i]] += rate*V[cells[i]]*rho[cells[i]];
+            hSource[cells[i]] += rate*V[cells[i]]*rho[cells[i]]*T_;
+        }
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/thermophysicalModels/thermalPorousZone/thermalModel/fixedTemperature/fixedTemperature.H b/src/thermophysicalModels/thermalPorousZone/thermalModel/fixedTemperature/fixedTemperature.H
new file mode 100644
index 0000000000..11c1c7c442
--- /dev/null
+++ b/src/thermophysicalModels/thermalPorousZone/thermalModel/fixedTemperature/fixedTemperature.H
@@ -0,0 +1,98 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::fixedTemperature
+
+Description
+    Fixed temperature model
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef fixedTemperature_H
+#define fixedTemperature_H
+
+#include "thermalModel.H"
+#include "autoPtr.H"
+#include "runTimeSelectionTables.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                      Class fixedTemperature Declaration
+\*---------------------------------------------------------------------------*/
+
+class fixedTemperature
+:
+    public thermalModel
+{
+
+protected:
+
+    // Protected data
+
+        //- Fixed temperature
+        const scalar T_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("fixedTemperature");
+
+
+    // Constructors
+
+        //- Construct from porous zone
+        fixedTemperature(const porousZone& pZone);
+
+
+    //- Destructor
+    virtual ~fixedTemperature();
+
+
+    // Member Functions
+
+        //- Add the thermal source to the enthalpy equation
+        virtual void addEnthalpySource
+        (
+            const basicThermo& thermo,
+            const volScalarField& rho,
+            fvScalarMatrix& hEqn
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
+
diff --git a/src/thermophysicalModels/thermalPorousZone/thermalModel/noThermalModel/noThermalModel.C b/src/thermophysicalModels/thermalPorousZone/thermalModel/noThermalModel/noThermalModel.C
new file mode 100644
index 0000000000..a7c3dfd017
--- /dev/null
+++ b/src/thermophysicalModels/thermalPorousZone/thermalModel/noThermalModel/noThermalModel.C
@@ -0,0 +1,74 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "noThermalModel.H"
+#include "addToRunTimeSelectionTable.H"
+#include "basicThermo.H"
+#include "volFields.H"
+#include "fvMatrices.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(noThermalModel, 0);
+
+    addToRunTimeSelectionTable
+    (
+        thermalModel,
+        noThermalModel,
+        pZone
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::noThermalModel::noThermalModel(const porousZone& pZone)
+:
+    thermalModel(pZone)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor    * * * * * * * * * * * * * * //
+
+Foam::noThermalModel::~noThermalModel()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::noThermalModel::addEnthalpySource
+(
+    const basicThermo&,
+    const volScalarField&,
+    fvScalarMatrix&
+) const
+{
+    // do nothing
+}
+
+
+// ************************************************************************* //
diff --git a/src/thermophysicalModels/thermalPorousZone/thermalModel/noThermalModel/noThermalModel.H b/src/thermophysicalModels/thermalPorousZone/thermalModel/noThermalModel/noThermalModel.H
new file mode 100644
index 0000000000..af23d09feb
--- /dev/null
+++ b/src/thermophysicalModels/thermalPorousZone/thermalModel/noThermalModel/noThermalModel.H
@@ -0,0 +1,90 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::noThermalModel
+
+Description
+    Dummy model for 'none' option
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef noThermalModel_H
+#define noThermalModel_H
+
+#include "thermalModel.H"
+#include "autoPtr.H"
+#include "runTimeSelectionTables.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                      Class noThermalModel Declaration
+\*---------------------------------------------------------------------------*/
+
+class noThermalModel
+:
+    public thermalModel
+{
+
+public:
+
+    //- Runtime type information
+    TypeName("none");
+
+
+    // Constructors
+
+        //- Construct from porous zone
+        noThermalModel(const porousZone& pZone);
+
+
+    //- Destructor
+    virtual ~noThermalModel();
+
+
+    // Member Functions
+
+        //- Add the thermal source to the enthalpy equation
+        virtual void addEnthalpySource
+        (
+            const basicThermo& thermo,
+            const volScalarField& rho,
+            fvScalarMatrix& hEqn
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
+
diff --git a/src/thermophysicalModels/thermalPorousZone/thermalModel/thermalModel/thermalModel.C b/src/thermophysicalModels/thermalPorousZone/thermalModel/thermalModel/thermalModel.C
new file mode 100644
index 0000000000..e8ba858b8c
--- /dev/null
+++ b/src/thermophysicalModels/thermalPorousZone/thermalModel/thermalModel/thermalModel.C
@@ -0,0 +1,63 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "thermalModel.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(thermalModel, 0);
+    defineRunTimeSelectionTable(thermalModel, pZone);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::thermalModel::thermalModel(const porousZone& pZone)
+:
+    pZone_(pZone),
+    coeffDict_(dictionary::null)
+{}
+
+
+Foam::thermalModel::thermalModel
+(
+    const porousZone& pZone,
+    const word& modelType
+)
+:
+    pZone_(pZone),
+    coeffDict_(pZone_.dict().subDict(modelType + "Coeffs"))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor    * * * * * * * * * * * * * * //
+
+Foam::thermalModel::~thermalModel()
+{}
+
+
+// ************************************************************************* //
diff --git a/src/thermophysicalModels/thermalPorousZone/thermalModel/thermalModel/thermalModel.H b/src/thermophysicalModels/thermalPorousZone/thermalModel/thermalModel/thermalModel.H
new file mode 100644
index 0000000000..c6bd861919
--- /dev/null
+++ b/src/thermophysicalModels/thermalPorousZone/thermalModel/thermalModel/thermalModel.H
@@ -0,0 +1,122 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::thermalModel
+
+Description
+    Base class for selecting the temperature specification models
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef thermalModel_H
+#define thermalModel_H
+
+#include "porousZone.H"
+#include "autoPtr.H"
+#include "runTimeSelectionTables.H"
+#include "volFieldsFwd.H"
+#include "fvMatricesFwd.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class basicThermo;
+
+/*---------------------------------------------------------------------------*\
+                      Class thermalModel Declaration
+\*---------------------------------------------------------------------------*/
+
+class thermalModel
+{
+
+protected:
+
+    // Protected data
+
+        //- Reference to the porous zone
+        const porousZone& pZone_;
+
+        //- Sub-model coefficients dictionary
+        const dictionary coeffDict_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("thermalModel");
+
+    //- Declare runtime constructor selection table
+    declareRunTimeSelectionTable
+    (
+        autoPtr,
+        thermalModel,
+        pZone,
+        (
+            const porousZone& pZone
+        ),
+        (pZone)
+    );
+
+
+    // Constructors
+
+        //- Construct null from porous zone
+        thermalModel(const porousZone& pZone);
+
+        //- Construct from porous zone and model type name
+        thermalModel(const porousZone& pZone, const word& modelType);
+
+
+    //- Destructor
+    virtual ~thermalModel();
+
+
+    //- Selector
+    static autoPtr<thermalModel> New(const porousZone& pZone);
+
+
+    // Member Functions
+
+        //- Add the thermal source to the enthalpy equation
+        virtual void addEnthalpySource
+        (
+            const basicThermo& thermo,
+            const volScalarField& rho,
+            fvScalarMatrix& hEqn
+        ) const = 0;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
+
diff --git a/src/thermophysicalModels/thermalPorousZone/thermalModel/thermalModel/thermalModelNew.C b/src/thermophysicalModels/thermalPorousZone/thermalModel/thermalModel/thermalModelNew.C
new file mode 100644
index 0000000000..e97685a3fc
--- /dev/null
+++ b/src/thermophysicalModels/thermalPorousZone/thermalModel/thermalModel/thermalModelNew.C
@@ -0,0 +1,58 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "thermalModel.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::thermalModel> Foam::thermalModel::New
+(
+    const porousZone& pZone
+)
+{
+    const word modelType(pZone.dict().lookup("thermalModel"));
+
+    Info<< "Selecting thermalModel " << modelType << endl;
+
+    pZoneConstructorTable::iterator cstrIter =
+        pZoneConstructorTablePtr_->find(modelType);
+
+    if (cstrIter == pZoneConstructorTablePtr_->end())
+    {
+        FatalErrorIn
+        (
+            "thermalModel::New(const porousZone&)"
+        )   << "Unknown thermalModel type "
+            << modelType << nl << nl
+            << "Valid thermalModel types are :" << endl
+            << pZoneConstructorTablePtr_->sortedToc()
+            << abort(FatalError);
+    }
+
+    return autoPtr<thermalModel>(cstrIter()(pZone));
+}
+
+
+// ************************************************************************* //
diff --git a/src/thermophysicalModels/thermalPorousZone/thermalPorousZone/thermalPorousZone.C b/src/thermophysicalModels/thermalPorousZone/thermalPorousZone/thermalPorousZone.C
index df6a32d84b..6e19d7a177 100644
--- a/src/thermophysicalModels/thermalPorousZone/thermalPorousZone/thermalPorousZone.C
+++ b/src/thermophysicalModels/thermalPorousZone/thermalPorousZone/thermalPorousZone.C
@@ -24,9 +24,9 @@ License
 \*----------------------------------------------------------------------------*/
 
 #include "thermalPorousZone.H"
-#include "fvMesh.H"
-#include "fvMatrices.H"
 #include "basicThermo.H"
+#include "volFields.H"
+#include "fvMatrices.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -38,29 +38,8 @@ Foam::thermalPorousZone::thermalPorousZone
 )
 :
     porousZone(key, mesh, dict),
-    T_("T", dimTemperature, -GREAT)
-{
-    if (const dictionary* dictPtr = dict.subDictPtr("thermalModel"))
-    {
-        const word thermalModel(dictPtr->lookup("type"));
-
-        if (thermalModel == "fixedTemperature")
-        {
-            dictPtr->lookup("T") >> T_;
-        }
-        else
-        {
-            FatalIOErrorIn
-            (
-                "thermalPorousZone::thermalPorousZone"
-                "(const keyType&, const fvMesh&, const dictionary&)",
-                *dictPtr
-            )   << "thermalModel " << thermalModel << " is not supported" << nl
-                << "    Supported thermalModels are: fixedTemperature"
-                << exit(FatalIOError);
-        }
-    }
-}
+    model_(thermalModel::New(*this))
+{}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
@@ -72,29 +51,7 @@ void Foam::thermalPorousZone::addEnthalpySource
     fvScalarMatrix& hEqn
 ) const
 {
-    const labelList& zones = this->zoneIds();
-    if (zones.empty() || T_.value() < 0.0)
-    {
-        return;
-    }
-
-    const scalarField& V = mesh().V();
-    scalarField& hDiag = hEqn.diag();
-    scalarField& hSource = hEqn.source();
-
-    // TODO: generalize for non-fixedTemperature methods
-    const scalar rate = 1e6;
-
-    forAll(zones, zoneI)
-    {
-        const labelList& cells = mesh().cellZones()[zones[zoneI]];
-
-        forAll(cells, i)
-        {
-            hDiag[cells[i]] += rate*V[cells[i]]*rho[cells[i]];
-            hSource[cells[i]] += rate*V[cells[i]]*rho[cells[i]]*T_.value();
-        }
-    }
+    model_->addEnthalpySource(thermo, rho, hEqn);
 }
 
 
diff --git a/src/thermophysicalModels/thermalPorousZone/thermalPorousZone/thermalPorousZone.H b/src/thermophysicalModels/thermalPorousZone/thermalPorousZone/thermalPorousZone.H
index 7278f4ec98..3264194939 100644
--- a/src/thermophysicalModels/thermalPorousZone/thermalPorousZone/thermalPorousZone.H
+++ b/src/thermophysicalModels/thermalPorousZone/thermalPorousZone/thermalPorousZone.H
@@ -29,7 +29,8 @@ Description
     equations.
 
 See Also
-    porousZone, thermalPorousZones and coordinateSystems
+    porousZone, thermalPorousZones and coordinateSystems with run-time
+    selectable thermal model
 
 SourceFiles
     thermalPorousZone.C
@@ -40,6 +41,7 @@ SourceFiles
 #define thermalPorousZone_H
 
 #include "porousZone.H"
+#include "thermalModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -50,7 +52,7 @@ class fvMesh;
 class basicThermo;
 
 /*---------------------------------------------------------------------------*\
-                        Class thermalPorousZone Declaration
+                       Class thermalPorousZone Declaration
 \*---------------------------------------------------------------------------*/
 
 class thermalPorousZone
@@ -59,23 +61,27 @@ class thermalPorousZone
 {
     // Private data
 
-        //- Temperature in the porous-zone
-        dimensionedScalar T_;
-
-
         //- Disallow default bitwise copy construct
         thermalPorousZone(const thermalPorousZone&);
 
         //- Disallow default bitwise assignment
         void operator=(const thermalPorousZone&);
 
+        //- Thermal model
+        autoPtr<thermalModel> model_;
+
 
 public:
 
     // Constructors
 
         //- Construct from components
-        thermalPorousZone(const keyType& key, const fvMesh&, const dictionary&);
+        thermalPorousZone
+        (
+            const keyType& key,
+            const fvMesh& mesh,
+            const dictionary& dict
+        );
 
         //- Return clone
         autoPtr<thermalPorousZone> clone() const
@@ -118,20 +124,6 @@ public:
 
     // Member Functions
 
-        // Access
-
-            //- Return the temperature in the porous-zone
-            const dimensionedScalar& T() const
-            {
-                return T_;
-            }
-
-            //- Edit access to the temperature in the porous-zone
-            dimensionedScalar& T()
-            {
-                return T_;
-            }
-
         //- Add the thermal source to the enthalpy equation
         void addEnthalpySource
         (
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/combustion/dieselFoam/aachenBomb/constant/sprayProperties b/tutorials/combustion/dieselFoam/aachenBomb/constant/sprayProperties
index 878439295f..c8cfa62f0b 100644
--- a/tutorials/combustion/dieselFoam/aachenBomb/constant/sprayProperties
+++ b/tutorials/combustion/dieselFoam/aachenBomb/constant/sprayProperties
@@ -127,8 +127,8 @@ hollowConeInjectorCoeffs
         {
             minValue        1e-06;
             maxValue        0.00015;
-            d               ( 0.00015 );
-            n               ( 3 );
+            d               0.00015;
+            n               3;
         }
 
         exponentialPDF
diff --git a/tutorials/combustion/fireFoam/les/smallPoolFire2D/system/fvSolution b/tutorials/combustion/fireFoam/les/smallPoolFire2D/system/fvSolution
index e71c1e85db..e97f1de618 100644
--- a/tutorials/combustion/fireFoam/les/smallPoolFire2D/system/fvSolution
+++ b/tutorials/combustion/fireFoam/les/smallPoolFire2D/system/fvSolution
@@ -20,10 +20,17 @@ solvers
 {
     rho
     {
-        solver          PCG;
-        preconditioner  DIC;
-        tolerance       0;
-        relTol          0;
+        solver              PCG;
+        preconditioner      DIC;
+        tolerance           1e-7;
+        relTol              0.1;
+    };
+
+    rhoFinal
+    {
+        $rho;
+        tolerance           1e-7;
+        relTol              0;
     };
 
     p
@@ -40,37 +47,13 @@ solvers
 
     pFinal
     {
-        solver              GAMG;
+        $p;
         tolerance           1e-7;
         relTol              0;
-        smoother            GaussSeidel;
-        cacheAgglomeration  true;
-        nCellsInCoarsestLevel 10;
-        agglomerator     faceAreaPair;
-        mergeLevels      1;
     };
 
 
-    ft
-    {
-        solver          smoothSolver;
-        smoother        GaussSeidel;
-        tolerance       1e-7;
-        relTol          0;
-        nSweeps         1;
-    };
-
-    fu
-    {
-        solver          smoothSolver;
-        smoother        GaussSeidel;
-        tolerance       1e-7;
-        relTol          0;
-        nSweeps         1;
-    };
-
-
-    U
+    "(U|ft|fu|k|hs)"
     {
         solver          smoothSolver;
         smoother        GaussSeidel;
@@ -79,31 +62,13 @@ solvers
         nSweeps         1;
     };
 
-    UFinal
+    "(U|ft|fu|k|hs)Final"
     {
-        solver          smoothSolver;
-        smoother        GaussSeidel;
+        $U;
         tolerance       1e-7;
         relTol          0;
-        nSweeps         1;
     };
 
-    k
-    {
-        solver          smoothSolver;
-        smoother        GaussSeidel;
-        tolerance       1e-7;
-        relTol          0;
-        nSweeps         1;
-    };
-
-    hs
-    {
-        solver          PBiCG;
-        preconditioner  DILU;
-        tolerance       1e-7;
-        relTol          0;
-    };
 
     Ii
     {
diff --git a/tutorials/combustion/reactingFoam/ras/counterFlowFlame2D/system/fvSolution b/tutorials/combustion/reactingFoam/ras/counterFlowFlame2D/system/fvSolution
index ad29357a94..38b8f5dbcc 100755
--- a/tutorials/combustion/reactingFoam/ras/counterFlowFlame2D/system/fvSolution
+++ b/tutorials/combustion/reactingFoam/ras/counterFlowFlame2D/system/fvSolution
@@ -32,7 +32,7 @@ solvers
         relTol           0.0;
     }
 
-    "(U|Yi|h|k|epsilon)"
+    "(U|Yi|hs|k|epsilon)"
     {
         solver          PBiCG;
         preconditioner  DILU;
diff --git a/tutorials/compressible/rhoCentralFoam/biconic25-55Run35/constant/thermophysicalProperties b/tutorials/compressible/rhoCentralFoam/biconic25-55Run35/constant/thermophysicalProperties
index 04ef086ac6..e30c05b865 100644
--- a/tutorials/compressible/rhoCentralFoam/biconic25-55Run35/constant/thermophysicalProperties
+++ b/tutorials/compressible/rhoCentralFoam/biconic25-55Run35/constant/thermophysicalProperties
@@ -15,11 +15,16 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-Pr              Pr [ 0 0 0 0 0 0 0 ] 0.72;
+Pr              0.72;
 
 thermoType      ePsiThermo<pureMixture<sutherlandTransport<specieThermo<janafThermo<perfectGas>>>>>;
 
-mixture         N2 1 28.01348 100 10000 1000 2.9525407 0.0013968838 -4.9262577e-07 7.8600091e-11 -4.6074978e-15 -923.93753 5.8718221 3.5309628 -0.0001236595 -5.0299339e-07 2.4352768e-09 -1.4087954e-12 -1046.9637 2.9674391 1.458e-06 110;
+mixture         N2
+     1 28.01348
+     100 10000 1000
+     2.9525407 0.0013968838 -4.9262577e-07 7.8600091e-11 -4.6074978e-15 -923.93753 5.8718221
+     3.5309628 -0.0001236595 -5.0299339e-07 2.4352768e-09 -1.4087954e-12 -1046.9637 2.9674391
+     1.458e-06 110;
 
 
 // ************************************************************************* //
diff --git a/tutorials/compressible/rhoPimpleFoam/ras/angledDuct/Allrun b/tutorials/compressible/rhoPimpleFoam/ras/angledDuct/Allrun
new file mode 100755
index 0000000000..13c8195c59
--- /dev/null
+++ b/tutorials/compressible/rhoPimpleFoam/ras/angledDuct/Allrun
@@ -0,0 +1,10 @@
+cd ${0%/*} || exit 1    # run from this directory                              
+
+m4 constant/polyMesh/blockMeshDict.m4 > constant/polyMesh/blockMeshDict
+
+# Source tutorial run functions
+. $WM_PROJECT_DIR/bin/tools/RunFunctions
+
+runApplication blockMesh
+runApplication rhoPimpleFoam
+
diff --git a/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/system/fvSolution b/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/system/fvSolution
index 4ebcefc07d..2e74c6af29 100644
--- a/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/system/fvSolution
+++ b/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/system/fvSolution
@@ -22,6 +22,13 @@ solvers
         solver          PCG;
         preconditioner  DIC;
         tolerance       1e-5;
+        relTol          0.1;
+    }
+
+    rhoFinal
+    {
+        $rho;
+        tolerance       1e-5;
         relTol          0;
     }
 
@@ -65,7 +72,7 @@ solvers
         mergeLevels     1;
     }
 
-    h
+    "(U|h|k|epsilon)"
     {
         solver          PBiCG;
         preconditioner  DILU;
@@ -73,26 +80,9 @@ solvers
         relTol          0.1;
     }
 
-    hFinal
+    "(U|h|k|epsilon)Final"
     {
-        solver          PBiCG;
-        preconditioner  DILU;
-        tolerance       1e-5;
-        relTol          0;
-    }
-
-    k
-    {
-        solver          PBiCG;
-        preconditioner  DILU;
-        tolerance       1e-5;
-        relTol          0;
-    }
-
-    epsilon
-    {
-        solver          PBiCG;
-        preconditioner  DILU;
+        $U;
         tolerance       1e-5;
         relTol          0;
     }
diff --git a/tutorials/compressible/rhoPorousSimpleFoam/angledDuctExplicit/Allrun b/tutorials/compressible/rhoPorousSimpleFoam/angledDuctExplicit/Allrun
new file mode 120000
index 0000000000..09efeaf930
--- /dev/null
+++ b/tutorials/compressible/rhoPorousSimpleFoam/angledDuctExplicit/Allrun
@@ -0,0 +1 @@
+../angledDuctImplicit/Allrun
\ No newline at end of file
diff --git a/tutorials/compressible/rhoPorousSimpleFoam/angledDuctExplicit/system/controlDict b/tutorials/compressible/rhoPorousSimpleFoam/angledDuctExplicit/system/controlDict
index 83de517420..4b8c2697ad 100644
--- a/tutorials/compressible/rhoPorousSimpleFoam/angledDuctExplicit/system/controlDict
+++ b/tutorials/compressible/rhoPorousSimpleFoam/angledDuctExplicit/system/controlDict
@@ -47,5 +47,25 @@ graphFormat     raw;
 
 runTimeModifiable true;
 
+functions
+{
+    residualControl1
+    {
+        type            residualControl;
+        functionObjectLibs ( "libjobControl.so" );
+        outputControl   timeStep;
+        outputInterval  1;
+
+        maxResiduals
+        {
+            p       1e-2;
+            U       1e-4;
+            T       1e-3;
+
+            // possibly check turbulence fields
+            "(k|epsilon|omega)" 1e-3;
+        }
+    }
+}
 
 // ************************************************************************* //
diff --git a/tutorials/compressible/rhoPorousSimpleFoam/angledDuctExplicit/system/fvSchemes b/tutorials/compressible/rhoPorousSimpleFoam/angledDuctExplicit/system/fvSchemes
index a1376603cf..4d843cfccb 100644
--- a/tutorials/compressible/rhoPorousSimpleFoam/angledDuctExplicit/system/fvSchemes
+++ b/tutorials/compressible/rhoPorousSimpleFoam/angledDuctExplicit/system/fvSchemes
@@ -34,7 +34,6 @@ divSchemes
     div(phi,h)      Gauss upwind;
     div(phi,epsilon) Gauss upwind;
     div(phi,k)      Gauss upwind;
-    div(phid,p)     Gauss linear;
 }
 
 laplacianSchemes
diff --git a/tutorials/compressible/rhoPorousSimpleFoam/angledDuctImplicit/Allrun b/tutorials/compressible/rhoPorousSimpleFoam/angledDuctImplicit/Allrun
new file mode 100755
index 0000000000..ad77d33e3d
--- /dev/null
+++ b/tutorials/compressible/rhoPorousSimpleFoam/angledDuctImplicit/Allrun
@@ -0,0 +1,10 @@
+cd ${0%/*} || exit 1    # run from this directory                              
+
+m4 constant/polyMesh/blockMeshDict.m4 > constant/polyMesh/blockMeshDict
+
+# Source tutorial run functions
+. $WM_PROJECT_DIR/bin/tools/RunFunctions
+
+runApplication blockMesh
+runApplication rhoPorousSimpleFoam
+
diff --git a/tutorials/compressible/rhoPorousSimpleFoam/angledDuctImplicit/constant/porousZones b/tutorials/compressible/rhoPorousSimpleFoam/angledDuctImplicit/constant/porousZones
index 4a14bc6eee..e4d251bf9e 100644
--- a/tutorials/compressible/rhoPorousSimpleFoam/angledDuctImplicit/constant/porousZones
+++ b/tutorials/compressible/rhoPorousSimpleFoam/angledDuctImplicit/constant/porousZones
@@ -31,11 +31,11 @@ FoamFile
             f   f [0 -1 0 0 0 0 0] (0 0 0);
         }
 
-        thermalModel
-        {
-            type fixedTemperature;
+        thermalModel none; // fixedTemperature;
 
-            T  T [0 0 0 1 0] 350;
+        fixedTemperatureCoeffs
+        {
+            T              350;
         }
     }
 )
diff --git a/tutorials/compressible/rhoPorousSimpleFoam/angledDuctImplicit/system/controlDict b/tutorials/compressible/rhoPorousSimpleFoam/angledDuctImplicit/system/controlDict
index 3b3db3d164..fd58f436b9 100644
--- a/tutorials/compressible/rhoPorousSimpleFoam/angledDuctImplicit/system/controlDict
+++ b/tutorials/compressible/rhoPorousSimpleFoam/angledDuctImplicit/system/controlDict
@@ -47,5 +47,25 @@ graphFormat     raw;
 
 runTimeModifiable true;
 
+functions
+{
+    residualControl1
+    {
+        type            residualControl;
+        functionObjectLibs ( "libjobControl.so" );
+        outputControl   timeStep;
+        outputInterval  1;
+
+        maxResiduals
+        {
+            p       1e-3;
+            U       1e-4;
+            T       1e-3;
+
+            // possibly check turbulence fields
+            "(k|epsilon|omega)" 1e-3;
+        }
+    }
+}
 
 // ************************************************************************* //
diff --git a/tutorials/compressible/rhoPorousSimpleFoam/angledDuctImplicit/system/fvSchemes b/tutorials/compressible/rhoPorousSimpleFoam/angledDuctImplicit/system/fvSchemes
index c03ccb3612..a432624360 100644
--- a/tutorials/compressible/rhoPorousSimpleFoam/angledDuctImplicit/system/fvSchemes
+++ b/tutorials/compressible/rhoPorousSimpleFoam/angledDuctImplicit/system/fvSchemes
@@ -34,7 +34,6 @@ divSchemes
     div(phi,h)      Gauss upwind;
     div(phi,epsilon) Gauss upwind;
     div(phi,k)      Gauss upwind;
-    div(phid,p)     Gauss linear;
 }
 
 laplacianSchemes
diff --git a/tutorials/compressible/rhoPorousSimpleFoam/angledDuctImplicit/system/fvSolution b/tutorials/compressible/rhoPorousSimpleFoam/angledDuctImplicit/system/fvSolution
index 19a717dfb3..e126f210c5 100644
--- a/tutorials/compressible/rhoPorousSimpleFoam/angledDuctImplicit/system/fvSolution
+++ b/tutorials/compressible/rhoPorousSimpleFoam/angledDuctImplicit/system/fvSolution
@@ -51,7 +51,7 @@ SIMPLE
 {
     nUCorrectors    2;
     nNonOrthogonalCorrectors 0;
-    pMin            pMin [ 1 -1 -2 0 0 0 0 ] 100;
+    pMin            pMin [ 1 -1 -2 0 0 ] 100;
 }
 
 relaxationFactors
diff --git a/tutorials/heatTransfer/buoyantBoussinesqPimpleFoam/hotRoom/0/T b/tutorials/heatTransfer/buoyantBoussinesqPimpleFoam/hotRoom/0/T
deleted file mode 100644
index 98cc191af7..0000000000
--- a/tutorials/heatTransfer/buoyantBoussinesqPimpleFoam/hotRoom/0/T
+++ /dev/null
@@ -1,443 +0,0 @@
-/*--------------------------------*- 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      T;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-dimensions      [0 0 0 1 0 0 0];
-
-internalField   uniform 300;
-
-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
-)
-;
-    }
-    ceiling
-    {
-        type            fixedValue;
-        value           uniform 300;
-    }
-    fixedWalls
-    {
-        type            zeroGradient;
-    }
-}
-
-// ************************************************************************* //
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/multiphase/interFoam/laminar/damBreak/0/alpha1 b/tutorials/heatTransfer/buoyantBoussinesqPimpleFoam/hotRoom/0/p_rgh
similarity index 71%
rename from tutorials/multiphase/interFoam/laminar/damBreak/0/alpha1
rename to tutorials/heatTransfer/buoyantBoussinesqPimpleFoam/hotRoom/0/p_rgh
index 057f3d1500..270f014558 100644
--- a/tutorials/multiphase/interFoam/laminar/damBreak/0/alpha1
+++ b/tutorials/heatTransfer/buoyantBoussinesqPimpleFoam/hotRoom/0/p_rgh
@@ -10,41 +10,35 @@ FoamFile
     version     2.0;
     format      ascii;
     class       volScalarField;
-    object      alpha;
+    object      p_rgh;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-dimensions      [0 0 0 0 0 0 0];
+dimensions      [0 2 -2 0 0 0 0];
 
 internalField   uniform 0;
 
 boundaryField
 {
-    leftWall
+    floor
     {
-        type            zeroGradient;
-    }
-
-    rightWall
-    {
-        type            zeroGradient;
-    }
-
-    lowerWall
-    {
-        type            zeroGradient;
-    }
-
-    atmosphere
-    {
-        type            inletOutlet;
-        inletValue      uniform 0;
+        type            buoyantPressure;
+        rho             rhok;
         value           uniform 0;
     }
 
-    defaultFaces
+    ceiling
     {
-        type            empty;
+        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/buoyantBoussinesqSimpleFoam/hotRoom/0/T b/tutorials/heatTransfer/buoyantPimpleFoam/hotRoom/0/p_rgh
similarity index 73%
rename from tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/hotRoom/0/T
rename to tutorials/heatTransfer/buoyantPimpleFoam/hotRoom/0/p_rgh
index 77b58a1a00..371d93c38a 100644
--- a/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/hotRoom/0/T
+++ b/tutorials/heatTransfer/buoyantPimpleFoam/hotRoom/0/p_rgh
@@ -10,29 +10,32 @@ FoamFile
     version     2.0;
     format      ascii;
     class       volScalarField;
-    object      T;
+    object      p_rgh;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-dimensions      [0 0 0 1 0 0 0];
+dimensions      [1 -1 -2 0 0 0 0];
 
-internalField   uniform 300;
+internalField   uniform 1e5;
 
 boundaryField
 {
     floor
     {
-        type            fixedValue;
-        value           uniform 300; 
+        type            buoyantPressure;
+        value           uniform 1e5;
     }
+
     ceiling
     {
-        type            fixedValue;
-        value           uniform 300;
+        type            buoyantPressure;
+        value           uniform 1e5;
     }
+
     fixedWalls
     {
-        type            zeroGradient;
+        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..a936e37b0f 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-3;
+            h       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..f43238ba65 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-3;
+            h       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..5dea594b3c 100644
--- a/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoom/system/controlDict
+++ b/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoom/system/controlDict
@@ -45,5 +45,26 @@ timePrecision   6;
 
 runTimeModifiable true;
 
+functions
+{
+    residualControl1
+    {
+        type            residualControl;
+        functionObjectLibs ( "libjobControl.so" );
+        outputControl   timeStep;
+        outputInterval  1;
+
+        maxResiduals
+        {
+            p_rgh   1e-2;
+            U       1e-3;
+            h       1e-3;
+            G       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..d66c31284a 100644
--- a/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoomFvDOM/system/controlDict
+++ b/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoomFvDOM/system/controlDict
@@ -45,5 +45,26 @@ timePrecision   6;
 
 runTimeModifiable true;
 
+functions
+{
+    residualControl1
+    {
+        type            residualControl;
+        functionObjectLibs ( "libjobControl.so" );
+        outputControl   timeStep;
+        outputInterval  1;
+
+        maxResiduals
+        {
+            p_rgh   1e-2;
+            U       1e-3;
+            h       1e-3;
+            G       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;
diff --git a/tutorials/incompressible/pimpleDyMFoam/movingCone/system/fvSolution b/tutorials/incompressible/pimpleDyMFoam/movingCone/system/fvSolution
index cc1c2297d1..6bc425e2f6 100644
--- a/tutorials/incompressible/pimpleDyMFoam/movingCone/system/fvSolution
+++ b/tutorials/incompressible/pimpleDyMFoam/movingCone/system/fvSolution
@@ -35,8 +35,7 @@ solvers
 
     pFinal
     {
-        solver          PCG;
-        preconditioner  DIC;
+        $p;
         tolerance       1e-06;
         relTol          0;
     }
@@ -46,6 +45,13 @@ solvers
         solver          PBiCG;
         preconditioner  DILU;
         tolerance       1e-05;
+        relTol          0.1;
+    }
+
+    UFinal
+    {
+        $U;
+        tolerance       1e-05;
         relTol          0;
     }
 
diff --git a/tutorials/incompressible/pimpleDyMFoam/wingMotion/wingMotion2D_pimpleDyMFoam/system/fvSolution b/tutorials/incompressible/pimpleDyMFoam/wingMotion/wingMotion2D_pimpleDyMFoam/system/fvSolution
index 2400b59b39..02becd019b 100644
--- a/tutorials/incompressible/pimpleDyMFoam/wingMotion/wingMotion2D_pimpleDyMFoam/system/fvSolution
+++ b/tutorials/incompressible/pimpleDyMFoam/wingMotion/wingMotion2D_pimpleDyMFoam/system/fvSolution
@@ -32,16 +32,9 @@ solvers
 
     pFinal
     {
-        solver           GAMG;
+        $p;
         tolerance        1e-7;
         relTol           0;
-        smoother         GaussSeidel;
-        nPreSweeps       0;
-        nPostSweeps      2;
-        cacheAgglomeration on;
-        agglomerator     faceAreaPair;
-        nCellsInCoarsestLevel 10;
-        mergeLevels      1;
     }
 
     "(U|k|omega)"
@@ -49,6 +42,13 @@ solvers
         solver          PBiCG;
         preconditioner  DILU;
         tolerance       1e-06;
+        relTol          0.1;
+    }
+
+    "(U|k|omega)Final"
+    {
+        $U;
+        tolerance       1e-06;
         relTol          0;
     }
 
diff --git a/tutorials/incompressible/pimpleFoam/t-junction-with-fan/system/fvSolution b/tutorials/incompressible/pimpleFoam/t-junction-with-fan/system/fvSolution
index 1b7fb06c4d..637611c63f 100644
--- a/tutorials/incompressible/pimpleFoam/t-junction-with-fan/system/fvSolution
+++ b/tutorials/incompressible/pimpleFoam/t-junction-with-fan/system/fvSolution
@@ -41,7 +41,7 @@ solvers
         mergeLevels     1;
     }
 
-    U
+    "(U|k|epsilon)"
     {
         solver          PBiCG;
         preconditioner  DILU;
@@ -49,26 +49,9 @@ solvers
         relTol          0.1;
     }
 
-    UFinal
+    "(U|k|epsilon)Final"
     {
-        solver          PBiCG;
-        preconditioner  DILU;
-        tolerance       1e-05;
-        relTol          0;
-    }
-
-    k
-    {
-        solver          PBiCG;
-        preconditioner  DILU;
-        tolerance       1e-05;
-        relTol          0;
-    }
-
-    epsilon
-    {
-        solver          PBiCG;
-        preconditioner  DILU;
+        $U;
         tolerance       1e-05;
         relTol          0;
     }
diff --git a/tutorials/incompressible/pimpleFoam/t-junction/system/fvSolution b/tutorials/incompressible/pimpleFoam/t-junction/system/fvSolution
index 1b7fb06c4d..637611c63f 100644
--- a/tutorials/incompressible/pimpleFoam/t-junction/system/fvSolution
+++ b/tutorials/incompressible/pimpleFoam/t-junction/system/fvSolution
@@ -41,7 +41,7 @@ solvers
         mergeLevels     1;
     }
 
-    U
+    "(U|k|epsilon)"
     {
         solver          PBiCG;
         preconditioner  DILU;
@@ -49,26 +49,9 @@ solvers
         relTol          0.1;
     }
 
-    UFinal
+    "(U|k|epsilon)Final"
     {
-        solver          PBiCG;
-        preconditioner  DILU;
-        tolerance       1e-05;
-        relTol          0;
-    }
-
-    k
-    {
-        solver          PBiCG;
-        preconditioner  DILU;
-        tolerance       1e-05;
-        relTol          0;
-    }
-
-    epsilon
-    {
-        solver          PBiCG;
-        preconditioner  DILU;
+        $U;
         tolerance       1e-05;
         relTol          0;
     }
diff --git a/tutorials/incompressible/porousSimpleFoam/angledDuctExplicit/Allrun b/tutorials/incompressible/porousSimpleFoam/angledDuctExplicit/Allrun
new file mode 120000
index 0000000000..09efeaf930
--- /dev/null
+++ b/tutorials/incompressible/porousSimpleFoam/angledDuctExplicit/Allrun
@@ -0,0 +1 @@
+../angledDuctImplicit/Allrun
\ No newline at end of file
diff --git a/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/Allrun b/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/Allrun
new file mode 100755
index 0000000000..689a832a04
--- /dev/null
+++ b/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/Allrun
@@ -0,0 +1,10 @@
+cd ${0%/*} || exit 1    # run from this directory                              
+
+m4 constant/polyMesh/blockMeshDict.m4 > constant/polyMesh/blockMeshDict
+
+# Source tutorial run functions
+. $WM_PROJECT_DIR/bin/tools/RunFunctions
+
+runApplication blockMesh
+runApplication porousSimpleFoam
+
diff --git a/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/constant/polyMesh/blockMeshDict.m4 b/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/constant/polyMesh/blockMeshDict.m4
new file mode 100644
index 0000000000..fb67da75da
--- /dev/null
+++ b/tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/constant/polyMesh/blockMeshDict.m4
@@ -0,0 +1,165 @@
+/*--------------------------------*- 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       dictionary;
+    object      blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// block definition for a porosity with an angled inlet/outlet
+// the porosity is not aligned with the main axes
+//
+dnl> -----------------------------------------------------------------
+dnl> <STANDARD DEFINTIONS>
+dnl>
+changecom(//)changequote([,]) dnl>
+define(calc, [esyscmd(perl -e 'print ($1)')]) dnl>
+define(VCOUNT, 0)  dnl>
+define(vlabel, [[// ]pt VCOUNT ($1) define($1, VCOUNT)define([VCOUNT], incr(VCOUNT))])  dnl>
+dnl>
+define(hex2D, hex ($1b $2b $3b $4b $1f $2f $3f $4f)) dnl>
+define(quad2D, ($1f $1b $2b $2f))  dnl>
+define(frontQuad, ($1f $2f $3f $4f)) dnl>
+define(backQuad, ($4b $3b $2b $1b)) dnl>
+dnl>
+dnl> </STANDARD DEFINTIONS>
+dnl> -----------------------------------------------------------------
+dnl>
+define(ncells, 20) dnl>
+define(ninlet, 15) dnl>
+define(nporo, 20) dnl>
+define(noutlet, 20) dnl>
+dnl>
+define(x0,0) dnl>
+define(y0,0) dnl>
+define(y0,0) dnl>
+define(Cos,0.7071067812) dnl>   == cos(45)
+define(Sin,0.7071067812) dnl>   == sin(45)
+dnl>
+define(width,50) dnl>
+define(zBack,calc(-width/2)) dnl>
+define(zFront,calc(width/2)) dnl>
+define(leninlet,150)dnl>
+define(lenporo,100)dnl>
+define(lenoutlet,100)dnl>
+dnl>
+define(xhyp,calc(Sin*width)) dnl>
+define(yhyp,calc(Cos*width)) dnl>
+define(xinlet,leninlet)dnl>
+define(xporo,calc(Cos*lenporo)) dnl>
+define(yporo,calc(Sin*lenporo)) dnl>
+define(xoutlet,calc(xporo + Cos*lenoutlet)) dnl>
+define(youtlet,calc(yporo + Sin*lenoutlet)) dnl>
+dnl>
+
+convertToMeters 0.001;
+
+vertices
+(
+    // inlet region
+    ( -xinlet  y0  zBack )  vlabel(in1b)
+    ( -xinlet yhyp  zBack ) vlabel(in2b)
+    ( -xinlet  y0  zFront )  vlabel(in1f)
+    ( -xinlet yhyp  zFront ) vlabel(in2f)
+
+    // join inlet->outlet
+    (  x0 y0  zBack )    vlabel(join1b)
+    ( -xhyp   yhyp  zBack ) vlabel(join2b)
+    (  x0 y0  zFront )    vlabel(join1f)
+    ( -xhyp   yhyp  zFront ) vlabel(join2f)
+
+    // porosity ends ->outlet
+    ( xporo yporo  zBack )  vlabel(poro1b)
+    ( calc(xporo - xhyp) calc(yporo + yhyp)  zBack )  vlabel(poro2b)
+    ( xporo yporo  zFront )  vlabel(poro1f)
+    ( calc(xporo - xhyp) calc(yporo + yhyp)  zFront )  vlabel(poro2f)
+
+    // outlet
+    ( xoutlet youtlet zBack ) vlabel(out1b)
+    ( calc(xoutlet - xhyp) calc(youtlet + yhyp) zBack ) vlabel(out2b)
+    ( xoutlet youtlet zFront ) vlabel(out1f)
+    ( calc(xoutlet - xhyp) calc(youtlet + yhyp) zFront ) vlabel(out2f)
+);
+
+blocks
+(
+    // inlet block
+    hex2D(in1, join1, join2, in2)
+    inlet ( ninlet ncells ncells ) simpleGrading (1 1 1)
+
+    // porosity block
+    hex2D(join1, poro1, poro2, join2)
+    porosity ( nporo ncells ncells ) simpleGrading (1 1 1)
+
+    // outlet block
+    hex2D(poro1, out1, out2, poro2)
+    outlet ( noutlet ncells ncells )  simpleGrading (1 1 1)
+);
+
+edges
+(
+);
+
+patches
+(
+    // is there no way of defining all my 'defaultFaces' to be 'wall'?
+    wall front
+    (
+    // inlet block
+    frontQuad(in1, join1, join2, in2)
+    // outlet block
+    frontQuad(poro1, out1, out2, poro2)
+    )
+
+    wall back
+    (
+    // inlet block
+    backQuad(in1, join1, join2, in2)
+    // outlet block
+    backQuad(poro1, out1, out2, poro2)
+    )
+
+    wall wall
+    (
+    // inlet block
+    quad2D(in1, join1)
+    quad2D(join2, in2)
+    // outlet block
+    quad2D(poro1, out1)
+    quad2D(out2, poro2)
+    )
+
+    wall porosityWall
+    (
+    // porosity block
+    frontQuad(join1, poro1, poro2, join2)
+    // porosity block
+    backQuad(join1, poro1, poro2, join2)
+    // porosity block
+    quad2D(join1, poro1)
+    quad2D(poro2, join2)
+    )
+
+    patch inlet
+    (
+    quad2D(in2, in1)
+    )
+
+    patch outlet
+    (
+    quad2D(out2, out1)
+    )
+);
+
+mergePatchPairs
+(
+);
+
+// ************************************************************************* //
diff --git a/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/coalCloud1Properties b/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/coalCloud1Properties
index 5624bb507b..8ea186c54b 100644
--- a/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/coalCloud1Properties
+++ b/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/coalCloud1Properties
@@ -33,7 +33,7 @@ PhaseChangeModel LiquidEvaporation;
 
 DevolatilisationModel ConstantRateDevolatilisation;
 
-SurfaceReactionModel COxidationDiffusionLimitedRate;
+SurfaceReactionModel COxidationKineticDiffusionLimitedRate;
 
 PostProcessingModel none;
 
@@ -167,10 +167,12 @@ ConstantRateDevolatilisationCoeffs
     volatileResidualCoeff 0.001;
 }
 
-COxidationDiffusionLimitedRateCoeffs
+COxidationKineticDiffusionLimitedRateCoeffs
 {
-    Sb              Sb [ 0 0 0 0 0 ] 1;
-    D               D [ 0 2 -1 0 0 ] 1.5e-05;
+    Sb              Sb [ 0 0 0 0 0 ] 2.667;
+    C1              C1 [ 0 0 0 0 0] 5.0E-12;
+    C2              C2 [ 0 0 0 0 0] 0.002;
+    E               E  [ 0 0 0 0 0] 7.9E+07;
 }
 
 
diff --git a/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/energySourceProperties b/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/energySourceProperties
index 4c28e875b3..ba68f2616f 100644
--- a/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/energySourceProperties
+++ b/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/constant/energySourceProperties
@@ -20,10 +20,10 @@ FoamFile
     {
         active          true;
         timeStart       0.15;
-        duration        0.1;
+        duration        0.2;
         selectionMode   cellSet;
         volumeMode      absolute;
-        fieldData       10000;
+        fieldData       15000;
         cellSet         ignitionCells;
     }
 );
diff --git a/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/system/fvSolution b/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/system/fvSolution
index ffb16df0c0..aa49989c5b 100644
--- a/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/system/fvSolution
+++ b/tutorials/lagrangian/coalChemistryFoam/simplifiedSiwek/system/fvSolution
@@ -22,13 +22,34 @@ solvers
         solver          PCG;
         preconditioner  DIC;
         tolerance       1e-05;
+        relTol          0.1;
+    }
+
+   "(rho|G)Final"
+    {
+        $rho;
+        tolerance       1e-05;
         relTol          0;
     }
 
     p
     {
-        solver          PBiCG;
-        preconditioner  DILU;
+        solver          GAMG;
+        tolerance       1e-06;
+        relTol          0.01;    
+        smoother        GaussSeidel;
+        nPreSweeps      0;             
+        nPostSweeps     2;             
+        nBottomSweeps   2;             
+        cacheAgglomeration false;      
+        nCellsInCoarsestLevel 10;      
+        agglomerator    faceAreaPair;  
+        mergeLevels     1;             
+    }
+
+    pFinal
+    {
+        $p;
         tolerance       1e-06;
         relTol          0;
     }
@@ -38,6 +59,13 @@ solvers
         solver          PBiCG;
         preconditioner  DILU;
         tolerance       1e-05;
+        relTol          0.1;
+    }
+
+    "(U|hs|R|k|epsilon)Final"
+    {
+        $U;
+        tolerance       1e-05;
         relTol          0;
     }
 
@@ -46,6 +74,14 @@ solvers
         solver          PBiCG;
         preconditioner  DILU;
         tolerance       1e-06;
+        relTol          0.1;
+    }
+
+    "(Yi|CO2|O2|N2|CH4|H2|H2O|CO)Final"
+    {
+        $Yi;
+        preconditioner  DILU;
+        tolerance       1e-06;
         relTol          0;
     }
 }
@@ -53,7 +89,7 @@ solvers
 PISO
 {
     transonic       yes;
-    nCorrectors     2;
+    nCorrectors     3;
     nNonOrthogonalCorrectors 0;
     momentumPredictor yes;
 }
diff --git a/tutorials/lagrangian/reactingParcelFilmFoam/hotBoxes/constant/surfaceFilmProperties b/tutorials/lagrangian/reactingParcelFilmFoam/hotBoxes/constant/surfaceFilmProperties
index 1045475fae..4099a1fff3 100644
--- a/tutorials/lagrangian/reactingParcelFilmFoam/hotBoxes/constant/surfaceFilmProperties
+++ b/tutorials/lagrangian/reactingParcelFilmFoam/hotBoxes/constant/surfaceFilmProperties
@@ -15,38 +15,18 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-//surfaceFilmModel kinematicSingleLayer;
 surfaceFilmModel thermoSingleLayer;
 
 filmRegionName  wallFilmRegion;
 
 active          true;
 
-kinematicSingleLayerCoeffs
-{
-    thermoModel constant;
-    rho0        rho0 [1 -3 0 0 0] 1000;
-    mu0         mu0 [1 -1 -1 0 0] 1e-3;
-    sigma0      sigma0 [1 0 -2 0 0] 0.07;
-
-    deltaStable deltaStable [0 1 0 0 0] 0.0005;
-    Cf          0.005;
-
-    injectionModel removeInjection; // none; // cloudInjection;
-}
-
 
 thermoSingleLayerCoeffs
 {
     thermoModel singleComponent; // constant
     liquid      H2O;
-//    thermoModel constant;
-//    rho0        rho0 [1 -3 0 0 0] 1000;
-//    mu0         mu0 [1 -1 -1 0 0] 1e-3;
-//    sigma0      sigma0 [1 0 -2 0 0] 0.07;
-//    cp0         cp0 [0 2 -2 -1 0] 4187;
 
-    deltaStable deltaStable [0 1 0 0 0] 0;
     Cf          0.005;
 
     upperSurfaceModels
@@ -75,7 +55,9 @@ thermoSingleLayerCoeffs
     {
         cloudName    reactingCloud1;
 
-        particlesPerParcel 100.0;
+        deltaStable         0.0;
+
+        particlesPerParcel  100.0;
 
         parcelPDF
         {
diff --git a/tutorials/lagrangian/reactingParcelFilmFoam/panel/constant/surfaceFilmProperties b/tutorials/lagrangian/reactingParcelFilmFoam/panel/constant/surfaceFilmProperties
index 9ff37a7d55..73c1c06c67 100644
--- a/tutorials/lagrangian/reactingParcelFilmFoam/panel/constant/surfaceFilmProperties
+++ b/tutorials/lagrangian/reactingParcelFilmFoam/panel/constant/surfaceFilmProperties
@@ -16,7 +16,6 @@ FoamFile
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 surfaceFilmModel kinematicSingleLayer;
-//surfaceFilmModel thermoSingleLayer;
 
 filmRegionName  wallFilmRegion;
 
@@ -24,30 +23,13 @@ active          true;
 
 kinematicSingleLayerCoeffs
 {
-    deltaStable deltaStable [0 1 0 0 0] 0;
+    thermoModel constant;
     Cf          0.005;
 
-    rho         rho [1 -3 0 0 0] 1000;
-    mu          mu [1 -1 -1 0 0] 1e-3;
-    sigma       sigma [1 0 -2 0 0] 0.07;
+    rho0        rho0 [1 -3 0 0 0] 1000;
+    mu0         mu0 [1 -1 -1 0 0] 1e-3;
+    sigma0      sigma0 [1 0 -2 0 0] 0.07;
 
     injectionModel none;
 }
 
-
-thermoSingleLayerCoeffs
-{
-    deltaStable deltaStable [0 1 0 0 0] 0;
-    Cf          0.005;
-    htcw        htcw [1 0 -3 -1 0] 1e-8; // 0.0;
-    htcs        htcs [1 0 -3 -1 0] 1e-8; // 0.0;
-
-    rho         rho [1 -3 0 0 0] 1000;
-    mu          mu [1 -1 -1 0 0] 1e-3;
-    sigma       sigma [1 0 -2 0 0] 0.07;
-    cp          cp [0 2 -2 -1 0] 4187;
-
-    injectionModel none;
-
-    phaseChangeModel none;
-}
diff --git a/tutorials/lagrangian/reactingParcelFilmFoam/panel/constant/thermophysicalProperties b/tutorials/lagrangian/reactingParcelFilmFoam/panel/constant/thermophysicalProperties
index acfbe00bce..e57087805d 100644
--- a/tutorials/lagrangian/reactingParcelFilmFoam/panel/constant/thermophysicalProperties
+++ b/tutorials/lagrangian/reactingParcelFilmFoam/panel/constant/thermophysicalProperties
@@ -23,17 +23,25 @@ foamChemistryFile "$FOAM_CASE/chemkin/foam.inp";
 
 foamChemistryThermoFile "$FOAM_CASE/chemkin/foam.dat";
 
-liquidComponents
-(
-    H2O
-);
-
-solidComponents
-();
-
 inertSpecie     N2;
 
-H2O             H2O defaultCoeffs;
+liquids
+{
+    liquidComponents
+    (
+        H2O
+    );
+
+    H2O             H2O defaultCoeffs;
+}
+
+solids
+{
+    solidComponents
+    ();
+}
+
+
 
 
 // ************************************************************************* //
diff --git a/tutorials/lagrangian/reactingParcelFilmFoam/panel/system/wallFilmRegion/fvSchemes b/tutorials/lagrangian/reactingParcelFilmFoam/panel/system/wallFilmRegion/fvSchemes
index b6a31afe89..0471d1d8f9 100644
--- a/tutorials/lagrangian/reactingParcelFilmFoam/panel/system/wallFilmRegion/fvSchemes
+++ b/tutorials/lagrangian/reactingParcelFilmFoam/panel/system/wallFilmRegion/fvSchemes
@@ -36,6 +36,7 @@ gradSchemes
 {
     default none;
     grad(pL) Gauss linear;
+    grad(sigmaf) Gauss linear;
     snGradCorr(deltaf) Gauss linear;
     snGradCorr(pp) Gauss linear;
     snGradCorr(pu) Gauss linear;
diff --git a/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/Allclean b/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/Allclean
new file mode 100755
index 0000000000..f3d1d6f142
--- /dev/null
+++ b/tutorials/multiphase/MRFMultiphaseInterFoam/mixerVessel2D/Allclean
@@ -0,0 +1,10 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # run from this directory
+
+# Source tutorial clean functions
+. $WM_PROJECT_DIR/bin/tools/CleanFunctions
+
+cleanCase
+rm 0/alphas > /dev/null 2>&1
+
+# ----------------------------------------------------------------- end-of-file
diff --git a/tutorials/multiphase/interFoam/laminar/Allclean b/tutorials/multiphase/interFoam/laminar/Allclean
index 448fc2d28d..909068ce6c 100755
--- a/tutorials/multiphase/interFoam/laminar/Allclean
+++ b/tutorials/multiphase/interFoam/laminar/Allclean
@@ -9,12 +9,7 @@ loseCases="damBreakFine"
 
 for case in $keepCases
 do
-    (cd $case && foamCleanTutorials)
-
-    if [ "$case" = "damBreak" ]
-    then
-        cp $case/0/alpha1.org $case/0/alpha1
-    fi
+    (cd $case && foamCleanTutorials && rm -rf 0/alpha1)
 done
 
 for case in $loseCases
diff --git a/tutorials/multiphase/interFoam/laminar/Allrun b/tutorials/multiphase/interFoam/laminar/Allrun
index 6a4557d324..50f30eed5f 100755
--- a/tutorials/multiphase/interFoam/laminar/Allrun
+++ b/tutorials/multiphase/interFoam/laminar/Allrun
@@ -31,6 +31,7 @@ setDamBreakFine ()
     cd damBreak || exit
 
     runApplication blockMesh
+    cp 0/alpha1.org 0/alpha1
     runApplication setFields
     runApplication $application
 )
@@ -43,10 +44,10 @@ cloneCase damBreak damBreakFine
 
     # Modify case
     setDamBreakFine
-    cp ../damBreak/0/alpha1.org 0/alpha1
 
     # And execute
     runApplication blockMesh
+    cp ../damBreak/0/alpha1.org 0/alpha1
     runApplication setFields
     runApplication decomposePar
     runParallel $application 4
@@ -58,6 +59,7 @@ cloneCase damBreak damBreakFine
     cd capillaryRise || exit
 
     runApplication blockMesh
+    cp 0/alpha1.org 0/alpha1
     runApplication setFields
     runApplication $application
 )