diff --git a/applications/solvers/heatTransfer/boussinesqBuoyantFoam/boussinesqBuoyantFoam.C b/applications/solvers/heatTransfer/boussinesqBuoyantFoam/boussinesqBuoyantFoam.C
index fbb31ef8d..a2daaf199 100644
--- a/applications/solvers/heatTransfer/boussinesqBuoyantFoam/boussinesqBuoyantFoam.C
+++ b/applications/solvers/heatTransfer/boussinesqBuoyantFoam/boussinesqBuoyantFoam.C
@@ -31,6 +31,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
+#include "pisoControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -41,6 +42,9 @@ int main(int argc, char *argv[])
 
 #   include "createTime.H"
 #   include "createMesh.H"
+
+    pisoControl piso(mesh);
+
 #   include "readTransportProperties.H"
 #   include "readGravitationalAcceleration.H"
 #   include "createFields.H"
@@ -50,11 +54,10 @@ int main(int argc, char *argv[])
 
     Info<< "\nStarting time loop\n" << endl;
 
-    for (runTime++; !runTime.end(); runTime++)
+    while (runTime.loop())
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-#       include "readPISOControls.H"
 #       include "CourantNo.H"
 
         fvVectorMatrix UEqn
@@ -70,7 +73,7 @@ int main(int argc, char *argv[])
 
         // --- PISO loop
 
-        for (int corr = 0; corr < nCorr; corr++)
+        while (piso.correct())
         {
             volScalarField rUA = 1.0/UEqn.A();
 
@@ -80,7 +83,7 @@ int main(int argc, char *argv[])
 
             adjustPhi(phi, U, p);
 
-            for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+            while (piso.correctNonOrthogonal())
             {
                 fvScalarMatrix pEqn
                 (
@@ -90,7 +93,7 @@ int main(int argc, char *argv[])
                 pEqn.setReference(pRefCell, pRefValue);
                 pEqn.solve();
 
-                if (nonOrth == nNonOrthCorr)
+                if (piso.finalNonOrthogonalIter())
                 {
                     phi -= pEqn.flux();
                 }
diff --git a/applications/solvers/heatTransfer/boussinesqBuoyantFoam/createFields.H b/applications/solvers/heatTransfer/boussinesqBuoyantFoam/createFields.H
index ce61861ee..bfac7811b 100644
--- a/applications/solvers/heatTransfer/boussinesqBuoyantFoam/createFields.H
+++ b/applications/solvers/heatTransfer/boussinesqBuoyantFoam/createFields.H
@@ -49,6 +49,7 @@
     label pRefCell = 0;
     scalar pRefValue = 0.0;
     setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue);
+    mesh.schemesDict().setFluxRequired(p.name());
 
     Info<< "Reading/calculating field rho\n" << endl;
     volScalarField rho
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqPisoFoam/UEqn.H b/applications/solvers/heatTransfer/buoyantBoussinesqPisoFoam/UEqn.H
index 35387f4ff..fa68de296 100644
--- a/applications/solvers/heatTransfer/buoyantBoussinesqPisoFoam/UEqn.H
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqPisoFoam/UEqn.H
@@ -9,7 +9,7 @@
 
     UEqn.relax();
 
-    if (momentumPredictor)
+    if (piso.momentumPredictor())
     {
         solve
         (
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqPisoFoam/buoyantBoussinesqPisoFoam.C b/applications/solvers/heatTransfer/buoyantBoussinesqPisoFoam/buoyantBoussinesqPisoFoam.C
index 73547b6a3..e91ead7cd 100644
--- a/applications/solvers/heatTransfer/buoyantBoussinesqPisoFoam/buoyantBoussinesqPisoFoam.C
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqPisoFoam/buoyantBoussinesqPisoFoam.C
@@ -48,6 +48,7 @@ Description
 #include "fvCFD.H"
 #include "singlePhaseTransportModel.H"
 #include "turbulenceModel.H"
+#include "pisoControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -56,6 +57,9 @@ int main(int argc, char *argv[])
 #   include "setRootCase.H"
 #   include "createTime.H"
 #   include "createMesh.H"
+
+    pisoControl piso(mesh);
+
 #   include "readGravitationalAcceleration.H"
 #   include "createFields.H"
 #   include "initContinuityErrs.H"
@@ -72,7 +76,6 @@ int main(int argc, char *argv[])
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
 #       include "readTimeControls.H"
-#       include "readPISOControls.H"
 #       include "CourantNo.H"
 #       include "setDeltaT.H"
 
@@ -80,7 +83,7 @@ int main(int argc, char *argv[])
 #       include "TEqn.H"
 
         // --- PISO loop
-        for (int corr = 0; corr < nCorr; corr++)
+        while (piso.correct())
         {
 #           include "pEqn.H"
         }
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqPisoFoam/createFields.H b/applications/solvers/heatTransfer/buoyantBoussinesqPisoFoam/createFields.H
index b2fa0af92..697bc6aba 100644
--- a/applications/solvers/heatTransfer/buoyantBoussinesqPisoFoam/createFields.H
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqPisoFoam/createFields.H
@@ -61,6 +61,7 @@
         pRefCell,
         pRefValue
     );
+    mesh.schemesDict().setFluxRequired(p.name());
 
 
     // Kinematic density for buoyancy force
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqPisoFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantBoussinesqPisoFoam/pEqn.H
index 02d86340c..2601d0eb6 100644
--- a/applications/solvers/heatTransfer/buoyantBoussinesqPisoFoam/pEqn.H
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqPisoFoam/pEqn.H
@@ -12,23 +12,16 @@
 
     phi = phiU + rUAf*fvc::interpolate(rhok)*(g & mesh.Sf());
 
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    while (piso.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
             fvm::laplacian(rUAf, p) == fvc::div(phi)
         );
 
-        if (corr == nCorr - 1 && nonOrth == nNonOrthCorr)
-        {
-            pEqn.solve(mesh.solutionDict().solver(p.name() + "Final"));
-        }
-        else
-        {
-            pEqn.solve(mesh.solutionDict().solver(p.name()));
-        }
+        pEqn.solve(mesh.solutionDict().solver(p.select(piso.finalInnerIter())));
 
-        if (nonOrth == nNonOrthCorr)
+        if (piso.finalNonOrthogonalIter())
         {
             phi -= pEqn.flux();
         }
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/TEqn.H b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/TEqn.H
index 7452a0495..3bdeea045 100644
--- a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/TEqn.H
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/TEqn.H
@@ -14,8 +14,7 @@
 
     TEqn.relax();
 
-    eqnResidual = TEqn.solve().initialResidual();
-    maxResidual = max(eqnResidual, maxResidual);
+    TEqn.solve();
 
     rhok = 1.0 - beta*(T - TRef);
 }
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/UEqn.H b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/UEqn.H
index bcf99e166..292f44b4b 100644
--- a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/UEqn.H
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/UEqn.H
@@ -8,7 +8,7 @@
 
     UEqn().relax();
 
-    eqnResidual = solve
+    solve
     (
         UEqn()
       ==
@@ -19,6 +19,4 @@
               - fvc::snGrad(p)*mesh.magSf()
             )
         )
-    ).initialResidual();
-
-    maxResidual = max(eqnResidual, maxResidual);
+    );
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/buoyantBoussinesqSimpleFoam.C b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/buoyantBoussinesqSimpleFoam.C
index acfd8ff48..341f0a3a4 100644
--- a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/buoyantBoussinesqSimpleFoam.C
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/buoyantBoussinesqSimpleFoam.C
@@ -48,6 +48,7 @@ Description
 #include "fvCFD.H"
 #include "singlePhaseTransportModel.H"
 #include "RASModel.H"
+#include "simpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -56,6 +57,9 @@ int main(int argc, char *argv[])
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createMesh.H"
+
+    simpleControl simple(mesh);
+
     #include "readGravitationalAcceleration.H"
     #include "createFields.H"
     #include "initContinuityErrs.H"
@@ -64,15 +68,10 @@ int main(int argc, char *argv[])
 
     Info<< "\nStarting time loop\n" << endl;
 
-    while (runTime.loop())
+    while (simple.loop())
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        #include "readSIMPLEControls.H"
-        #include "initConvergenceCheck.H"
-
-        p.storePrevIter();
-
         // Pressure-velocity SIMPLE corrector
         {
             #include "UEqn.H"
@@ -87,8 +86,6 @@ int main(int argc, char *argv[])
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
             << "  ClockTime = " << runTime.elapsedClockTime() << " s"
             << nl << endl;
-
-        #include "convergenceCheck.H"
     }
 
     Info<< "End\n" << endl;
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/convergenceCheck.H b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/convergenceCheck.H
deleted file mode 100644
index 895806319..000000000
--- a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/convergenceCheck.H
+++ /dev/null
@@ -1,9 +0,0 @@
-// check convergence
-
-if (maxResidual < convergenceCriterion)
-{
-    Info<< "reached convergence criterion: " << convergenceCriterion << endl;
-    runTime.writeAndEnd();
-    Info<< "latestTime = " << runTime.timeName() << endl;
-}
-
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/createFields.H b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/createFields.H
index 0e668e9d9..663863637 100644
--- a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/createFields.H
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/createFields.H
@@ -64,7 +64,7 @@
         pRefCell,
         pRefValue
     );
-
+    mesh.schemesDict().setFluxRequired(p.name());
 
     // Kinematic density for buoyancy force
     volScalarField rhok
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/initConvergenceCheck.H b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/initConvergenceCheck.H
deleted file mode 100644
index b56197f22..000000000
--- a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/initConvergenceCheck.H
+++ /dev/null
@@ -1,7 +0,0 @@
-// initialize values for convergence checks
-
-    scalar eqnResidual = 1, maxResidual = 0;
-    scalar convergenceCriterion = 0;
-
-    simple.readIfPresent("convergence", convergenceCriterion);
-
diff --git a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pEqn.H
index 6408c7ccc..256a336f3 100644
--- a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pEqn.H
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pEqn.H
@@ -12,7 +12,7 @@
         rUAf*fvc::interpolate(rhok)*(g & mesh.Sf());
     phi += buoyancyPhi;
 
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    while (simple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
@@ -21,18 +21,9 @@
 
         pEqn.setReference(pRefCell, pRefValue);
 
-        // retain the residual from the first iteration
-        if (nonOrth == 0)
-        {
-            eqnResidual = pEqn.solve().initialResidual();
-            maxResidual = max(eqnResidual, maxResidual);
-        }
-        else
-        {
-            pEqn.solve();
-        }
+        pEqn.solve();
 
-        if (nonOrth == nNonOrthCorr)
+        if (simple.finalNonOrthogonalIter())
         {
             // Calculate the conservative fluxes
             phi -= pEqn.flux();
diff --git a/applications/solvers/heatTransfer/buoyantPisoFoam/UEqn.H b/applications/solvers/heatTransfer/buoyantPisoFoam/UEqn.H
index d4878d063..230f18661 100644
--- a/applications/solvers/heatTransfer/buoyantPisoFoam/UEqn.H
+++ b/applications/solvers/heatTransfer/buoyantPisoFoam/UEqn.H
@@ -9,7 +9,7 @@
 
     UEqn.relax();
 
-    if (momentumPredictor)
+    if (piso.momentumPredictor())
     {
         solve
         (
diff --git a/applications/solvers/heatTransfer/buoyantPisoFoam/buoyantPisoFoam.C b/applications/solvers/heatTransfer/buoyantPisoFoam/buoyantPisoFoam.C
index 2a277a9bf..05dd6a70f 100644
--- a/applications/solvers/heatTransfer/buoyantPisoFoam/buoyantPisoFoam.C
+++ b/applications/solvers/heatTransfer/buoyantPisoFoam/buoyantPisoFoam.C
@@ -37,6 +37,7 @@ Description
 #include "basicRhoThermo.H"
 #include "turbulenceModel.H"
 #include "fixedGradientFvPatchFields.H"
+#include "pisoControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -45,6 +46,9 @@ int main(int argc, char *argv[])
 #   include "setRootCase.H"
 #   include "createTime.H"
 #   include "createMesh.H"
+
+    pisoControl piso(mesh);
+
 #   include "readGravitationalAcceleration.H"
 #   include "createFields.H"
 #   include "initContinuityErrs.H"
@@ -59,7 +63,6 @@ int main(int argc, char *argv[])
     while (runTime.run())
     {
 #       include "readTimeControls.H"
-#       include "readPISOControls.H"
 #       include "compressibleCourantNo.H"
 #       include "setDeltaT.H"
 
@@ -75,7 +78,7 @@ int main(int argc, char *argv[])
 
         // --- PISO loop
 
-        for (int corr = 0; corr < nCorr; corr++)
+        while (piso.correct())
         {
 #           include "pEqn.H"
         }
diff --git a/applications/solvers/heatTransfer/buoyantPisoFoam/createFields.H b/applications/solvers/heatTransfer/buoyantPisoFoam/createFields.H
index b8ac5595e..616d1c170 100644
--- a/applications/solvers/heatTransfer/buoyantPisoFoam/createFields.H
+++ b/applications/solvers/heatTransfer/buoyantPisoFoam/createFields.H
@@ -65,3 +65,5 @@
     dimensionedScalar initialMass = fvc::domainIntegrate(rho);
 
     dimensionedScalar totalVolume = sum(mesh.V());
+
+    mesh.schemesDict().setFluxRequired(p.name());
diff --git a/applications/solvers/heatTransfer/buoyantPisoFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantPisoFoam/pEqn.H
index f78fb9483..64edba028 100644
--- a/applications/solvers/heatTransfer/buoyantPisoFoam/pEqn.H
+++ b/applications/solvers/heatTransfer/buoyantPisoFoam/pEqn.H
@@ -23,7 +23,7 @@
 
     phi = phiU + rhorUAf*fvc::interpolate(rho)*(g & mesh.Sf());
 
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    while (piso.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
@@ -32,16 +32,9 @@
           - fvm::laplacian(rhorUAf, p)
         );
 
-        if (corr == nCorr - 1 && nonOrth == nNonOrthCorr)
-        {
-            pEqn.solve(mesh.solutionDict().solver(p.name() + "Final"));
-        }
-        else
-        {
-            pEqn.solve(mesh.solutionDict().solver(p.name()));
-        }
+        pEqn.solve(mesh.solutionDict().solver(p.select(piso.finalInnerIter())));
 
-        if (nonOrth == nNonOrthCorr)
+        if (piso.finalNonOrthogonalIter())
         {
             phi += pEqn.flux();
         }
diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/UEqn.H b/applications/solvers/heatTransfer/buoyantSimpleFoam/UEqn.H
index 0170d6e38..560988360 100644
--- a/applications/solvers/heatTransfer/buoyantSimpleFoam/UEqn.H
+++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/UEqn.H
@@ -8,7 +8,7 @@
 
     UEqn().relax();
 
-    eqnResidual = solve
+    solve
     (
         UEqn()
      ==
@@ -17,6 +17,4 @@
             fvc::interpolate(rho)*(g & mesh.Sf())
           - fvc::snGrad(p)*mesh.magSf()
         )
-    ).initialResidual();
-
-    maxResidual = max(eqnResidual, maxResidual);
+    );
diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C b/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C
index 913f3634b..9ff38ec08 100644
--- a/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C
+++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C
@@ -33,6 +33,7 @@ Description
 #include "basicPsiThermo.H"
 #include "RASModel.H"
 #include "fixedGradientFvPatchFields.H"
+#include "simpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -41,6 +42,9 @@ int main(int argc, char *argv[])
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createMesh.H"
+
+    simpleControl simple(mesh);
+
     #include "readGravitationalAcceleration.H"
     #include "createFields.H"
     #include "initContinuityErrs.H"
@@ -49,16 +53,10 @@ int main(int argc, char *argv[])
 
     Info<< "\nStarting time loop\n" << endl;
 
-    while (runTime.loop())
+    while (simple.loop())
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        #include "readSIMPLEControls.H"
-        #include "initConvergenceCheck.H"
-
-        p.storePrevIter();
-        rho.storePrevIter();
-
         // Pressure-velocity SIMPLE corrector
         {
             #include "UEqn.H"
@@ -73,8 +71,6 @@ int main(int argc, char *argv[])
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
             << "  ClockTime = " << runTime.elapsedClockTime() << " s"
             << nl << endl;
-
-        #include "convergenceCheck.H"
     }
 
     Info<< "End\n" << endl;
diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/convergenceCheck.H b/applications/solvers/heatTransfer/buoyantSimpleFoam/convergenceCheck.H
deleted file mode 100644
index 895806319..000000000
--- a/applications/solvers/heatTransfer/buoyantSimpleFoam/convergenceCheck.H
+++ /dev/null
@@ -1,9 +0,0 @@
-// check convergence
-
-if (maxResidual < convergenceCriterion)
-{
-    Info<< "reached convergence criterion: " << convergenceCriterion << endl;
-    runTime.writeAndEnd();
-    Info<< "latestTime = " << runTime.timeName() << endl;
-}
-
diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H b/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H
index b304ace10..a428c20ee 100644
--- a/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H
+++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H
@@ -67,3 +67,5 @@
 
 
     dimensionedScalar initialMass = fvc::domainIntegrate(rho);
+
+    mesh.schemesDict().setFluxRequired(p.name());
diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/hEqn.H b/applications/solvers/heatTransfer/buoyantSimpleFoam/hEqn.H
index 0c507ce3d..b91a17ee7 100644
--- a/applications/solvers/heatTransfer/buoyantSimpleFoam/hEqn.H
+++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/hEqn.H
@@ -11,8 +11,7 @@
 
     hEqn.relax();
 
-    eqnResidual = hEqn.solve().initialResidual();
-    maxResidual = max(eqnResidual, maxResidual);
+    hEqn.solve();
 
     thermo.correct();
 }
diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/initConvergenceCheck.H b/applications/solvers/heatTransfer/buoyantSimpleFoam/initConvergenceCheck.H
deleted file mode 100644
index c920b6708..000000000
--- a/applications/solvers/heatTransfer/buoyantSimpleFoam/initConvergenceCheck.H
+++ /dev/null
@@ -1,7 +0,0 @@
-    // initialize values for convergence checks
-
-    scalar eqnResidual = 1, maxResidual = 0;
-    scalar convergenceCriterion = 0;
-
-    simple.readIfPresent("convergence", convergenceCriterion);
-
diff --git a/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H b/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H
index 11768ea2f..0321e0748 100644
--- a/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H
+++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H
@@ -14,7 +14,7 @@
         rhorUAf*fvc::interpolate(rho)*(g & mesh.Sf());
     phi += buoyancyPhi;
 
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    while (simple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
@@ -23,18 +23,9 @@
 
         pEqn.setReference(pRefCell, pRefValue);
 
-        // retain the residual from the first iteration
-        if (nonOrth == 0)
-        {
-            eqnResidual = pEqn.solve().initialResidual();
-            maxResidual = max(eqnResidual, maxResidual);
-        }
-        else
-        {
-            pEqn.solve();
-        }
+        pEqn.solve();
 
-        if (nonOrth == nNonOrthCorr)
+        if (simple.finalNonOrthogonalIter())
         {
             // For closed-volume cases adjust the pressure and density levels
             // to obey overall mass continuity
diff --git a/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/buoyantSimpleRadiationFoam.C b/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/buoyantSimpleRadiationFoam.C
index f5ee87a62..3f451f99c 100644
--- a/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/buoyantSimpleRadiationFoam.C
+++ b/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/buoyantSimpleRadiationFoam.C
@@ -35,6 +35,7 @@ Description
 #include "RASModel.H"
 #include "fixedGradientFvPatchFields.H"
 #include "radiationModel.H"
+#include "simpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -43,6 +44,9 @@ int main(int argc, char *argv[])
 #   include "setRootCase.H"
 #   include "createTime.H"
 #   include "createMesh.H"
+
+    simpleControl simple(mesh);
+
 #   include "readGravitationalAcceleration.H"
 #   include "createFields.H"
 #   include "createRadiationModel.H"
@@ -52,16 +56,10 @@ int main(int argc, char *argv[])
 
     Info<< "\nStarting time loop\n" << endl;
 
-    while (runTime.loop())
+    while (simple.loop())
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-#       include "readSIMPLEControls.H"
-#       include "initConvergenceCheck.H"
-
-        p.storePrevIter();
-        rho.storePrevIter();
-
         // Pressure-velocity SIMPLE corrector
         {
 #           include "UEqn.H"
@@ -76,8 +74,6 @@ int main(int argc, char *argv[])
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
             << "  ClockTime = " << runTime.elapsedClockTime() << " s"
             << nl << endl;
-
-#       include "convergenceCheck.H"
     }
 
     Info<< "End\n" << endl;
diff --git a/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/hEqn.H b/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/hEqn.H
index 24b17645d..57de44e4a 100644
--- a/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/hEqn.H
+++ b/applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/hEqn.H
@@ -12,8 +12,7 @@
 
     hEqn.relax();
 
-    eqnResidual = hEqn.solve().initialResidual();
-    maxResidual = max(eqnResidual, maxResidual);
+    hEqn.solve();
 
     thermo.correct();
 
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H
index 49bffaa4d..5e584cce7 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H
@@ -15,3 +15,5 @@
     volScalarField& h = thermo.h();
 
     const dimensionedScalar massIni("massIni", dimMass, initialMassFluid[i]);
+
+    mesh.schemesDict().setFluxRequired(p.name());
diff --git a/applications/solvers/heatTransfer/chtMultiRegionSimpleFoam/fluid/setRegionFluidFields.H b/applications/solvers/heatTransfer/chtMultiRegionSimpleFoam/fluid/setRegionFluidFields.H
index 656fb42e5..4c1b4d921 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionSimpleFoam/fluid/setRegionFluidFields.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionSimpleFoam/fluid/setRegionFluidFields.H
@@ -22,3 +22,5 @@
 
     const label pRefCell = pRefCellFluid[i];
     const scalar pRefValue = pRefValueFluid[i];
+
+    mesh.schemesDict().setFluxRequired(p.name());
diff --git a/tutorials/heatTransfer/boussinesqBuoyantFoam/heatedCavity/constant/polyMesh/boundary b/tutorials/heatTransfer/boussinesqBuoyantFoam/heatedCavity/constant/polyMesh/boundary
index 70a873217..8c1c27580 100644
--- a/tutorials/heatTransfer/boussinesqBuoyantFoam/heatedCavity/constant/polyMesh/boundary
+++ b/tutorials/heatTransfer/boussinesqBuoyantFoam/heatedCavity/constant/polyMesh/boundary
@@ -1,9 +1,9 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
 | \\      /  F ield         | foam-extend: Open Source CFD                    |
-|  \\    /   O peration     | Version:     3.2                                |
+|  \\    /   O peration     | Version:     4.0                                |
 |   \\  /    A nd           | Web:         http://www.foam-extend.org         |
-|    \\/     M anipulation  |                                                 |
+|    \\/     M anipulation  | For copyright notice see file Copyright         |
 \*---------------------------------------------------------------------------*/
 FoamFile
 {
@@ -44,7 +44,6 @@ FoamFile
     frontAndBack
     {
         type            empty;
-        physicalType    empty;
         nFaces          800;
         startFace       840;
     }
diff --git a/tutorials/heatTransfer/boussinesqBuoyantFoam/heatedCavity/system/fvSchemes b/tutorials/heatTransfer/boussinesqBuoyantFoam/heatedCavity/system/fvSchemes
index 46bdf14a2..96ec38f0e 100644
--- a/tutorials/heatTransfer/boussinesqBuoyantFoam/heatedCavity/system/fvSchemes
+++ b/tutorials/heatTransfer/boussinesqBuoyantFoam/heatedCavity/system/fvSchemes
@@ -52,10 +52,4 @@ snGradSchemes
     default         corrected;
 }
 
-fluxRequired
-{
-    default         no;
-    p;
-}
-
 // ************************************************************************* //
diff --git a/tutorials/heatTransfer/boussinesqBuoyantFoam/heatedCavity/system/fvSolution b/tutorials/heatTransfer/boussinesqBuoyantFoam/heatedCavity/system/fvSolution
index c91fea86e..e231e3689 100644
--- a/tutorials/heatTransfer/boussinesqBuoyantFoam/heatedCavity/system/fvSolution
+++ b/tutorials/heatTransfer/boussinesqBuoyantFoam/heatedCavity/system/fvSolution
@@ -16,9 +16,29 @@ FoamFile
 
 solvers
 {
-    p               ICCG 1e-06 0;
-    U               BICCG 1e-05 0;
-    T               BICCG 1e-05 0;
+    p
+    {
+        solver          PCG;
+        preconditioner  DIC;
+        tolerance       1e-06;
+        relTol          0;
+    }
+
+    U
+    {
+        solver          PBiCG;
+        preconditioner  DILU;
+        tolerance       1e-05;
+        relTol          0;
+    }
+
+    T
+    {
+        solver          PBiCG;
+        preconditioner  DILU;
+        tolerance       1e-05;
+        relTol          0;
+    }
 }
 
 PISO
diff --git a/tutorials/heatTransfer/buoyantBoussinesqPisoFoam/hotRoom/Allclean b/tutorials/heatTransfer/buoyantBoussinesqPisoFoam/hotRoom/Allclean
index 83fcec83e..93433055a 100755
--- a/tutorials/heatTransfer/buoyantBoussinesqPisoFoam/hotRoom/Allclean
+++ b/tutorials/heatTransfer/buoyantBoussinesqPisoFoam/hotRoom/Allclean
@@ -4,4 +4,5 @@
 . $WM_PROJECT_DIR/bin/tools/CleanFunctions
 
 cleanCase
+rm -f constant/polyMesh/boundary
 cp 0/T.org 0/T
diff --git a/tutorials/heatTransfer/buoyantBoussinesqPisoFoam/hotRoom/system/fvSchemes b/tutorials/heatTransfer/buoyantBoussinesqPisoFoam/hotRoom/system/fvSchemes
index 83620fd7b..1bafd5c74 100644
--- a/tutorials/heatTransfer/buoyantBoussinesqPisoFoam/hotRoom/system/fvSchemes
+++ b/tutorials/heatTransfer/buoyantBoussinesqPisoFoam/hotRoom/system/fvSchemes
@@ -34,7 +34,7 @@ divSchemes
     div(phi,epsilon) Gauss upwind;
     div(phi,R)      Gauss upwind;
     div(R)          Gauss linear;
-    div((nuEff*dev(grad(U).T()))) Gauss linear;
+    div((nuEff*dev2(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes
@@ -58,10 +58,4 @@ snGradSchemes
     default         corrected;
 }
 
-fluxRequired
-{
-    default         no;
-    p               ;
-}
-
 // ************************************************************************* //
diff --git a/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/hotRoom/Allclean b/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/hotRoom/Allclean
index 83fcec83e..93433055a 100755
--- a/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/hotRoom/Allclean
+++ b/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/hotRoom/Allclean
@@ -4,4 +4,5 @@
 . $WM_PROJECT_DIR/bin/tools/CleanFunctions
 
 cleanCase
+rm -f constant/polyMesh/boundary
 cp 0/T.org 0/T
diff --git a/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/hotRoom/system/fvSchemes b/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/hotRoom/system/fvSchemes
index 120dc9023..f34199123 100644
--- a/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/hotRoom/system/fvSchemes
+++ b/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/hotRoom/system/fvSchemes
@@ -34,7 +34,7 @@ divSchemes
     div(phi,epsilon) Gauss upwind;
     div(phi,R)      Gauss upwind;
     div(R)          Gauss linear;
-    div((nuEff*dev(grad(U).T()))) Gauss linear;
+    div((nuEff*dev2(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes
@@ -58,10 +58,4 @@ snGradSchemes
     default         corrected;
 }
 
-fluxRequired
-{
-    default         no;
-    p               ;
-}
-
 // ************************************************************************* //
diff --git a/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/iglooWithFridges/constant/polyMesh/boundary b/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/iglooWithFridges/constant/polyMesh/boundary
index 2c17efa04..2d9d01dd7 100644
--- a/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/iglooWithFridges/constant/polyMesh/boundary
+++ b/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/iglooWithFridges/constant/polyMesh/boundary
@@ -1,9 +1,9 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
 | \\      /  F ield         | foam-extend: Open Source CFD                    |
-|  \\    /   O peration     | Version:     3.2                                |
+|  \\    /   O peration     | Version:     4.0                                |
 |   \\  /    A nd           | Web:         http://www.foam-extend.org         |
-|    \\/     M anipulation  |                                                 |
+|    \\/     M anipulation  | For copyright notice see file Copyright         |
 \*---------------------------------------------------------------------------*/
 FoamFile
 {
@@ -15,43 +15,61 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-6
+9
 (
     maxY
     {
         type            empty;
-        nFaces          400;
-        startFace       22800;
+        nFaces          0;
+        startFace       59065;
     }
     minX
     {
         type            empty;
-        nFaces          400;
-        startFace       23200;
+        nFaces          0;
+        startFace       59065;
     }
     maxX
     {
         type            empty;
-        nFaces          400;
-        startFace       23600;
+        nFaces          0;
+        startFace       59065;
     }
     minY
     {
         type            empty;
-        nFaces          400;
-        startFace       24000;
+        nFaces          0;
+        startFace       59065;
     }
     ground
     {
         type            wall;
-        nFaces          400;
-        startFace       24400;
+        nFaces          590;
+        startFace       59065;
     }
     maxZ
     {
         type            empty;
-        nFaces          400;
-        startFace       24800;
+        nFaces          0;
+        startFace       59655;
+    }
+    igloo_region0
+    {
+        type            wall;
+        nFaces          2260;
+        startFace       59655;
+    }
+    twoFridgeFreezers_seal_0
+    {
+        type            wall;
+        nFaces          1344;
+        startFace       61915;
+    }
+    twoFridgeFreezers_herring_1
+    {
+        type            wall;
+        nFaces          1116;
+        startFace       63259;
     }
 )
 
diff --git a/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/iglooWithFridges/system/fvSchemes b/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/iglooWithFridges/system/fvSchemes
index 120dc9023..f34199123 100644
--- a/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/iglooWithFridges/system/fvSchemes
+++ b/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/iglooWithFridges/system/fvSchemes
@@ -34,7 +34,7 @@ divSchemes
     div(phi,epsilon) Gauss upwind;
     div(phi,R)      Gauss upwind;
     div(R)          Gauss linear;
-    div((nuEff*dev(grad(U).T()))) Gauss linear;
+    div((nuEff*dev2(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes
@@ -58,10 +58,4 @@ snGradSchemes
     default         corrected;
 }
 
-fluxRequired
-{
-    default         no;
-    p               ;
-}
-
 // ************************************************************************* //
diff --git a/tutorials/heatTransfer/buoyantPisoFoam/hotRoom/system/fvSchemes b/tutorials/heatTransfer/buoyantPisoFoam/hotRoom/system/fvSchemes
index fbd5c70cf..9832e1b35 100644
--- a/tutorials/heatTransfer/buoyantPisoFoam/hotRoom/system/fvSchemes
+++ b/tutorials/heatTransfer/buoyantPisoFoam/hotRoom/system/fvSchemes
@@ -35,7 +35,7 @@ divSchemes
     div(phi,R)      Gauss upwind;
     div(phiU,p)     Gauss linear;
     div(R)          Gauss linear;
-    div((muEff*dev2(grad(U).T()))) Gauss linear;
+    div((muEff*dev2(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes
@@ -59,10 +59,4 @@ snGradSchemes
     default         corrected;
 }
 
-fluxRequired
-{
-    default         no;
-    p               ;
-}
-
 // ************************************************************************* //
diff --git a/tutorials/heatTransfer/buoyantSimpleFoam/hotRoom/0/epsilon b/tutorials/heatTransfer/buoyantSimpleFoam/hotRoom/0/epsilon
index 8f61fdf4f..31bf6d6e0 100644
--- a/tutorials/heatTransfer/buoyantSimpleFoam/hotRoom/0/epsilon
+++ b/tutorials/heatTransfer/buoyantSimpleFoam/hotRoom/0/epsilon
@@ -1,9 +1,9 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
 | \\      /  F ield         | foam-extend: Open Source CFD                    |
-|  \\    /   O peration     | Version:     3.2                                |
+|  \\    /   O peration     | Version:     4.0                                |
 |   \\  /    A nd           | Web:         http://www.foam-extend.org         |
-|    \\/     M anipulation  |                                                 |
+|    \\/     M anipulation  | For copyright notice see file Copyright         |
 \*---------------------------------------------------------------------------*/
 FoamFile
 {
@@ -24,27 +24,31 @@ boundaryField
     floor
     {
         type            compressible::epsilonWallFunction;
+        refValue        uniform 0;
+        value           uniform 0.01;
         Cmu             0.09;
         kappa           0.41;
         E               9.8;
-        value           uniform 0.01;
     }
     ceiling
     {
         type            compressible::epsilonWallFunction;
+        refValue        uniform 0;
+        value           uniform 0.01;
         Cmu             0.09;
         kappa           0.41;
         E               9.8;
-        value           uniform 0.01;
     }
     fixedWalls
     {
         type            compressible::epsilonWallFunction;
+        refValue        uniform 0;
+        value           uniform 0.01;
         Cmu             0.09;
         kappa           0.41;
         E               9.8;
-        value           uniform 0.01;
     }
 }
 
+
 // ************************************************************************* //
diff --git a/tutorials/heatTransfer/buoyantSimpleFoam/hotRoom/0/k b/tutorials/heatTransfer/buoyantSimpleFoam/hotRoom/0/k
index 22b76a2b4..bf5caa993 100644
--- a/tutorials/heatTransfer/buoyantSimpleFoam/hotRoom/0/k
+++ b/tutorials/heatTransfer/buoyantSimpleFoam/hotRoom/0/k
@@ -1,9 +1,9 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
 | \\      /  F ield         | foam-extend: Open Source CFD                    |
-|  \\    /   O peration     | Version:     3.2                                |
+|  \\    /   O peration     | Version:     4.0                                |
 |   \\  /    A nd           | Web:         http://www.foam-extend.org         |
-|    \\/     M anipulation  |                                                 |
+|    \\/     M anipulation  | For copyright notice see file Copyright         |
 \*---------------------------------------------------------------------------*/
 FoamFile
 {
@@ -38,4 +38,5 @@ boundaryField
     }
 }
 
+
 // ************************************************************************* //
diff --git a/tutorials/heatTransfer/buoyantSimpleFoam/hotRoom/system/fvSchemes b/tutorials/heatTransfer/buoyantSimpleFoam/hotRoom/system/fvSchemes
index 096978d29..9f1a8cd96 100644
--- a/tutorials/heatTransfer/buoyantSimpleFoam/hotRoom/system/fvSchemes
+++ b/tutorials/heatTransfer/buoyantSimpleFoam/hotRoom/system/fvSchemes
@@ -33,7 +33,7 @@ divSchemes
     div(phi,epsilon) Gauss upwind;
     div(phi,R)      Gauss upwind;
     div(R)          Gauss linear;
-    div((muEff*dev2(grad(U).T()))) Gauss linear;
+    div((muEff*dev2(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes
@@ -57,10 +57,4 @@ snGradSchemes
     default         corrected;
 }
 
-fluxRequired
-{
-    default         no;
-    p;
-}
-
 // ************************************************************************* //
diff --git a/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoom/0/epsilon b/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoom/0/epsilon
index d60a224ed..7a8fa872c 100644
--- a/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoom/0/epsilon
+++ b/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoom/0/epsilon
@@ -1,9 +1,9 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
 | \\      /  F ield         | foam-extend: Open Source CFD                    |
-|  \\    /   O peration     | Version:     3.2                                |
+|  \\    /   O peration     | Version:     4.0                                |
 |   \\  /    A nd           | Web:         http://www.foam-extend.org         |
-|    \\/     M anipulation  |                                                 |
+|    \\/     M anipulation  | For copyright notice see file Copyright         |
 \*---------------------------------------------------------------------------*/
 FoamFile
 {
@@ -24,35 +24,40 @@ boundaryField
     box
     {
         type            compressible::epsilonWallFunction;
+        refValue        uniform 0;
+        value           uniform 0.01;
         Cmu             0.09;
         kappa           0.41;
         E               9.8;
-        value           uniform 0.01;
     }
     floor
     {
         type            compressible::epsilonWallFunction;
+        refValue        uniform 0;
+        value           uniform 0.01;
         Cmu             0.09;
         kappa           0.41;
         E               9.8;
-        value           uniform 0.01;
     }
     ceiling
     {
         type            compressible::epsilonWallFunction;
+        refValue        uniform 0;
+        value           uniform 0.01;
         Cmu             0.09;
         kappa           0.41;
         E               9.8;
-        value           uniform 0.01;
     }
     fixedWalls
     {
         type            compressible::epsilonWallFunction;
+        refValue        uniform 0;
+        value           uniform 0.01;
         Cmu             0.09;
         kappa           0.41;
         E               9.8;
-        value           uniform 0.01;
     }
 }
 
+
 // ************************************************************************* //
diff --git a/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoom/0/k b/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoom/0/k
index 8c2bd1711..a52bacc9b 100644
--- a/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoom/0/k
+++ b/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoom/0/k
@@ -1,9 +1,9 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
 | \\      /  F ield         | foam-extend: Open Source CFD                    |
-|  \\    /   O peration     | Version:     3.2                                |
+|  \\    /   O peration     | Version:     4.0                                |
 |   \\  /    A nd           | Web:         http://www.foam-extend.org         |
-|    \\/     M anipulation  |                                                 |
+|    \\/     M anipulation  | For copyright notice see file Copyright         |
 \*---------------------------------------------------------------------------*/
 FoamFile
 {
@@ -43,4 +43,5 @@ boundaryField
     }
 }
 
+
 // ************************************************************************* //
diff --git a/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoom/system/fvSchemes b/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoom/system/fvSchemes
index fe8e1b62c..5c7b4c2dc 100644
--- a/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoom/system/fvSchemes
+++ b/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoom/system/fvSchemes
@@ -33,7 +33,7 @@ divSchemes
     div(phi,epsilon) Gauss upwind;
     div(phi,R)      Gauss upwind;
     div(R)          Gauss linear;
-    div((muEff*dev2(grad(U).T()))) Gauss linear;
+    div((muEff*dev2(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes
@@ -58,10 +58,4 @@ snGradSchemes
     default         corrected;
 }
 
-fluxRequired
-{
-    default         no;
-    p;
-}
-
 // ************************************************************************* //
diff --git a/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoomFvDOM/system/fvSchemes b/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoomFvDOM/system/fvSchemes
index 5b4dda112..4560d748a 100644
--- a/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoomFvDOM/system/fvSchemes
+++ b/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoomFvDOM/system/fvSchemes
@@ -35,7 +35,7 @@ divSchemes
     div(phi,R)      Gauss upwind;
     div(R)          Gauss linear;
     div(Ji,Ii_h)    Gauss linearUpwind Gauss linear; //Gauss upwind;
-    div((muEff*dev2(grad(U).T()))) Gauss linear;
+    div((muEff*dev2(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes
@@ -60,10 +60,4 @@ snGradSchemes
     default         corrected;
 }
 
-fluxRequired
-{
-    default         no;
-    p;
-}
-
 // ************************************************************************* //
diff --git a/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoomFvDOMSpecularCeiling/system/fvSchemes b/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoomFvDOMSpecularCeiling/system/fvSchemes
index 5b4dda112..4560d748a 100644
--- a/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoomFvDOMSpecularCeiling/system/fvSchemes
+++ b/tutorials/heatTransfer/buoyantSimpleRadiationFoam/hotRadiationRoomFvDOMSpecularCeiling/system/fvSchemes
@@ -35,7 +35,7 @@ divSchemes
     div(phi,R)      Gauss upwind;
     div(R)          Gauss linear;
     div(Ji,Ii_h)    Gauss linearUpwind Gauss linear; //Gauss upwind;
-    div((muEff*dev2(grad(U).T()))) Gauss linear;
+    div((muEff*dev2(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes
@@ -60,10 +60,4 @@ snGradSchemes
     default         corrected;
 }
 
-fluxRequired
-{
-    default         no;
-    p;
-}
-
 // ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/bottomAir/fvSchemes b/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/bottomAir/fvSchemes
index cbdb5c43e..cf84871b7 100644
--- a/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/bottomAir/fvSchemes
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/bottomAir/fvSchemes
@@ -36,7 +36,7 @@ divSchemes
     div(phi,epsilon) Gauss upwind;
     div(phi,R)      Gauss upwind;
     div(R)          Gauss linear;
-    div((muEff*dev2(grad(U).T()))) Gauss linear;
+    div((muEff*dev2(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes
@@ -61,10 +61,4 @@ snGradSchemes
     default         limited 0.333;
 }
 
-fluxRequired
-{
-    default         no;
-    p;
-}
-
 // ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/fvSchemes b/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/fvSchemes
index dffa496fd..947c673e9 100644
--- a/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/fvSchemes
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/fvSchemes
@@ -36,7 +36,7 @@ divSchemes
     div(phi,epsilon) Gauss upwind;
     div(phi,R)      Gauss upwind;
     div(R)          Gauss linear;
-    div((muEff*dev2(grad(U).T()))) Gauss linear;
+    div((muEff*dev2(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes
@@ -60,10 +60,4 @@ snGradSchemes
     default         limited 0.333;
 }
 
-fluxRequired
-{
-    default         no;
-    p;
-}
-
 // ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/heater/fvSchemes b/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/heater/fvSchemes
index c45fbbe2f..ced621f2d 100644
--- a/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/heater/fvSchemes
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/multiRegionHeater/system/heater/fvSchemes
@@ -45,9 +45,4 @@ snGradSchemes
     default         limited 0.333;
 }
 
-fluxRequired
-{
-    default         no;
-}
-
 // ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/constant/polyMesh/boundary b/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/constant/polyMesh/boundary
index a5a1a2810..2c11571c2 100644
--- a/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/constant/polyMesh/boundary
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/constant/polyMesh/boundary
@@ -1,9 +1,9 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
 | \\      /  F ield         | foam-extend: Open Source CFD                    |
-|  \\    /   O peration     | Version:     3.2                                |
+|  \\    /   O peration     | Version:     4.0                                |
 |   \\  /    A nd           | Web:         http://www.foam-extend.org         |
-|    \\/     M anipulation  |                                                 |
+|    \\/     M anipulation  | For copyright notice see file Copyright         |
 \*---------------------------------------------------------------------------*/
 FoamFile
 {
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/bottomAir/fvSchemes b/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/bottomAir/fvSchemes
index 994d37fa3..8819b4e9a 100644
--- a/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/bottomAir/fvSchemes
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/bottomAir/fvSchemes
@@ -36,7 +36,7 @@ divSchemes
     div(phi,epsilon) Gauss upwind;
     div(phi,R)      Gauss upwind;
     div(R)          Gauss linear;
-    div((muEff*dev2(grad(U).T()))) Gauss linear;
+    div((muEff*dev2(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes
@@ -60,10 +60,4 @@ snGradSchemes
     default         limited 0.333;
 }
 
-fluxRequired
-{
-    default         no;
-    p;
-}
-
 // ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/fvSchemes b/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/fvSchemes
index c4d92e41f..0e3d160c6 100644
--- a/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/fvSchemes
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/fvSchemes
@@ -38,8 +38,4 @@ snGradSchemes
 {
 }
 
-fluxRequired
-{
-}
-
 // ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/heater/fvSchemes b/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/heater/fvSchemes
index c45fbbe2f..ced621f2d 100644
--- a/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/heater/fvSchemes
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/heater/fvSchemes
@@ -45,9 +45,4 @@ snGradSchemes
     default         limited 0.333;
 }
 
-fluxRequired
-{
-    default         no;
-}
-
 // ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/leftSolid/fvSchemes b/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/leftSolid/fvSchemes
index c45fbbe2f..ced621f2d 100644
--- a/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/leftSolid/fvSchemes
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/leftSolid/fvSchemes
@@ -45,9 +45,4 @@ snGradSchemes
     default         limited 0.333;
 }
 
-fluxRequired
-{
-    default         no;
-}
-
 // ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/rightSolid/fvSchemes b/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/rightSolid/fvSchemes
index c45fbbe2f..ced621f2d 100644
--- a/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/rightSolid/fvSchemes
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/rightSolid/fvSchemes
@@ -45,9 +45,4 @@ snGradSchemes
     default         limited 0.333;
 }
 
-fluxRequired
-{
-    default         no;
-}
-
 // ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/topAir/fvSchemes b/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/topAir/fvSchemes
index 994d37fa3..8819b4e9a 100644
--- a/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/topAir/fvSchemes
+++ b/tutorials/heatTransfer/chtMultiRegionFoam/snappyMultiRegionHeater/system/topAir/fvSchemes
@@ -36,7 +36,7 @@ divSchemes
     div(phi,epsilon) Gauss upwind;
     div(phi,R)      Gauss upwind;
     div(R)          Gauss linear;
-    div((muEff*dev2(grad(U).T()))) Gauss linear;
+    div((muEff*dev2(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes
@@ -60,10 +60,4 @@ snGradSchemes
     default         limited 0.333;
 }
 
-fluxRequired
-{
-    default         no;
-    p;
-}
-
 // ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/fvSchemes b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/fvSchemes
index c4d92e41f..0e3d160c6 100644
--- a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/fvSchemes
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/fvSchemes
@@ -38,8 +38,4 @@ snGradSchemes
 {
 }
 
-fluxRequired
-{
-}
-
 // ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/heater/fvSchemes b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/heater/fvSchemes
index 57b67d2f5..de7b6e7c8 100644
--- a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/heater/fvSchemes
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/heater/fvSchemes
@@ -45,9 +45,4 @@ snGradSchemes
     default         limited 0.333;
 }
 
-fluxRequired
-{
-    default         no;
-}
-
 // ************************************************************************* //
diff --git a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/topAir/fvSchemes b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/topAir/fvSchemes
index 9e51abb1e..c4b9fd122 100644
--- a/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/topAir/fvSchemes
+++ b/tutorials/heatTransfer/chtMultiRegionSimpleFoam/multiRegionHeater/system/topAir/fvSchemes
@@ -34,7 +34,7 @@ divSchemes
     div(phi,epsilon) Gauss upwind;
     div(phi,R)      Gauss upwind;
     div(R)          Gauss linear;
-    div((muEff*dev2(grad(U).T()))) Gauss linear;
+    div((muEff*dev2(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes
@@ -58,10 +58,4 @@ snGradSchemes
     default         limited 0.333;
 }
 
-fluxRequired
-{
-    default         no;
-    p;
-}
-
 // ************************************************************************* //