diff --git a/applications/solvers/coupled/blockCoupledScalarTransportFoam/blockCoupledScalarTransportFoam.C b/applications/solvers/coupled/blockCoupledScalarTransportFoam/blockCoupledScalarTransportFoam.C
index 413ee4e5a..66a486307 100644
--- a/applications/solvers/coupled/blockCoupledScalarTransportFoam/blockCoupledScalarTransportFoam.C
+++ b/applications/solvers/coupled/blockCoupledScalarTransportFoam/blockCoupledScalarTransportFoam.C
@@ -39,7 +39,7 @@ Description
 #include "foamTime.H"
 #include "fvMesh.H"
 #include "fvBlockMatrix.H"
-
+#include "simpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -49,6 +49,9 @@ int main(int argc, char *argv[])
 #   include "setRootCase.H"
 #   include "createTime.H"
 #   include "createMesh.H"
+
+    simpleControl simple(mesh);
+
 #   include "createFields.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -57,13 +60,11 @@ int main(int argc, char *argv[])
 
 #   include "CourantNo.H"
 
-    for (runTime++; !runTime.end(); runTime++)
+    while (simple.loop())
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-#       include "readSIMPLEControls.H"
-
-        for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+        while (simple.correctNonOrthogonal())
         {
             fvScalarMatrix TEqn
             (
diff --git a/applications/solvers/coupled/conjugateHeatFoam/UEqn.H b/applications/solvers/coupled/conjugateHeatFoam/UEqn.H
index df6f90ac0..10aeeb0c6 100644
--- a/applications/solvers/coupled/conjugateHeatFoam/UEqn.H
+++ b/applications/solvers/coupled/conjugateHeatFoam/UEqn.H
@@ -9,7 +9,7 @@
 
     UEqn.relax();
 
-    if (momentumPredictor)
+    if (piso.momentumPredictor())
     {
         solve
         (
diff --git a/applications/solvers/coupled/conjugateHeatFoam/conjugateHeatFoam.C b/applications/solvers/coupled/conjugateHeatFoam/conjugateHeatFoam.C
index 89bfc93d2..5b00fd9da 100644
--- a/applications/solvers/coupled/conjugateHeatFoam/conjugateHeatFoam.C
+++ b/applications/solvers/coupled/conjugateHeatFoam/conjugateHeatFoam.C
@@ -38,6 +38,7 @@ Description
 #include "thermalModel.H"
 #include "singlePhaseTransportModel.H"
 #include "RASModel.H"
+#include "pisoControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -47,6 +48,9 @@ int main(int argc, char *argv[])
 #   include "createTime.H"
 #   include "createFluidMesh.H"
 #   include "createSolidMesh.H"
+
+    pisoControl piso(mesh);
+
 #   include "readGravitationalAcceleration.H"
 #   include "createFields.H"
 #   include "createSolidFields.H"
@@ -64,7 +68,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"
 
@@ -75,7 +78,7 @@ int main(int argc, char *argv[])
 
         p_rgh.storePrevIter();
 
-        for (int corr = 0; corr < nCorr; corr++)
+        while (piso.correct())
         {
 #           include "pEqn.H"
         }
diff --git a/applications/solvers/coupled/conjugateHeatFoam/createFields.H b/applications/solvers/coupled/conjugateHeatFoam/createFields.H
index 720ef6deb..fc5504c68 100644
--- a/applications/solvers/coupled/conjugateHeatFoam/createFields.H
+++ b/applications/solvers/coupled/conjugateHeatFoam/createFields.H
@@ -71,6 +71,7 @@
         pRefCell,
         pRefValue
     );
+    mesh.schemesDict().setFluxRequired(p_rgh.name());
 
     autoPtr<incompressible::RASModel> turbulence
     (
diff --git a/applications/solvers/coupled/conjugateHeatFoam/pEqn.H b/applications/solvers/coupled/conjugateHeatFoam/pEqn.H
index 57e9a156d..e85345491 100644
--- a/applications/solvers/coupled/conjugateHeatFoam/pEqn.H
+++ b/applications/solvers/coupled/conjugateHeatFoam/pEqn.H
@@ -10,7 +10,7 @@
     surfaceScalarField buoyancyPhi(rUAf*ghf*fvc::snGrad(rhok)*mesh.magSf());
     phi -= buoyancyPhi;
 
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    while (piso.correctNonOrthogonal())
     {
         fvScalarMatrix p_rghEqn
         (
@@ -19,16 +19,15 @@
 
         p_rghEqn.setReference(pRefCell, pRefValue);
 
-        if (corr == nCorr-1 && nonOrth == nNonOrthCorr)
-        {
-            p_rghEqn.solve(mesh.solutionDict().solver(p_rgh.name() + "Final"));
-        }
-        else
-        {
-            p_rghEqn.solve(mesh.solutionDict().solver(p_rgh.name()));
-        }
+        p_rghEqn.solve
+        (
+            mesh.solutionDict().solver
+            (
+                p_rgh.select(piso.finalInnerIter())
+            )
+        );
 
-        if (nonOrth == nNonOrthCorr)
+        if (piso.finalNonOrthogonalIter())
         {
             // Calculate the conservative fluxes
             phi -= p_rghEqn.flux();
diff --git a/applications/solvers/coupled/conjugateHeatSimpleFoam/UEqn.H b/applications/solvers/coupled/conjugateHeatSimpleFoam/UEqn.H
index 59633497a..32510f030 100644
--- a/applications/solvers/coupled/conjugateHeatSimpleFoam/UEqn.H
+++ b/applications/solvers/coupled/conjugateHeatSimpleFoam/UEqn.H
@@ -9,7 +9,7 @@
 
     UEqn().relax();
 
-    eqnResidual = solve
+    solve
     (
         UEqn()
       ==
@@ -22,6 +22,4 @@
                 )*mesh.magSf()
             )
         )
-    ).initialResidual();
-
-    maxResidual = max(eqnResidual, maxResidual);
+    );
diff --git a/applications/solvers/coupled/conjugateHeatSimpleFoam/conjugateHeatSimpleFoam.C b/applications/solvers/coupled/conjugateHeatSimpleFoam/conjugateHeatSimpleFoam.C
index affeaecf9..2ba10da9b 100644
--- a/applications/solvers/coupled/conjugateHeatSimpleFoam/conjugateHeatSimpleFoam.C
+++ b/applications/solvers/coupled/conjugateHeatSimpleFoam/conjugateHeatSimpleFoam.C
@@ -38,6 +38,7 @@ Description
 #include "thermalModel.H"
 #include "singlePhaseTransportModel.H"
 #include "RASModel.H"
+#include "simpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -47,6 +48,9 @@ int main(int argc, char *argv[])
 #   include "createTime.H"
 #   include "createFluidMesh.H"
 #   include "createSolidMesh.H"
+
+    simpleControl simple(mesh);
+
 #   include "readGravitationalAcceleration.H"
 #   include "createFields.H"
 #   include "createSolidFields.H"
@@ -56,13 +60,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"
-
         // Detach patches
 #       include "detachPatches.H"
 
diff --git a/applications/solvers/coupled/conjugateHeatSimpleFoam/createFields.H b/applications/solvers/coupled/conjugateHeatSimpleFoam/createFields.H
index 159d180cc..a60a1d76b 100644
--- a/applications/solvers/coupled/conjugateHeatSimpleFoam/createFields.H
+++ b/applications/solvers/coupled/conjugateHeatSimpleFoam/createFields.H
@@ -71,6 +71,7 @@
         pRefCell,
         pRefValue
     );
+    mesh.schemesDict().setFluxRequired(p_rgh.name());
 
     autoPtr<incompressible::RASModel> turbulence
     (
diff --git a/applications/solvers/coupled/conjugateHeatSimpleFoam/initConvergenceCheck.H b/applications/solvers/coupled/conjugateHeatSimpleFoam/initConvergenceCheck.H
deleted file mode 100644
index b56197f22..000000000
--- a/applications/solvers/coupled/conjugateHeatSimpleFoam/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/coupled/conjugateHeatSimpleFoam/pEqn.H b/applications/solvers/coupled/conjugateHeatSimpleFoam/pEqn.H
index 559b052c6..dddf81e5b 100644
--- a/applications/solvers/coupled/conjugateHeatSimpleFoam/pEqn.H
+++ b/applications/solvers/coupled/conjugateHeatSimpleFoam/pEqn.H
@@ -11,7 +11,7 @@
     surfaceScalarField buoyancyPhi(rUAf*ghf*fvc::snGrad(rhok)*mesh.magSf());
     phi -= buoyancyPhi;
 
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    while (simple.correctNonOrthogonal())
     {
         fvScalarMatrix p_rghEqn
         (
@@ -20,18 +20,9 @@
 
         p_rghEqn.setReference(pRefCell, pRefValue);
 
-        // retain the residual from the first iteration
-        if (nonOrth == 0)
-        {
-            eqnResidual = p_rghEqn.solve().initialResidual();
-            maxResidual = max(eqnResidual, maxResidual);
-        }
-        else
-        {
-            p_rghEqn.solve();
-        }
+        p_rghEqn.solve();
 
-        if (nonOrth == nNonOrthCorr)
+        if (simple.finalNonOrthogonalIter())
         {
             // Calculate the conservative fluxes
             phi -= p_rghEqn.flux();
diff --git a/tutorials/coupled/blockCoupledScalarTransportFoam/blockCoupledSwirlTest/constant/polyMesh/boundary b/tutorials/coupled/blockCoupledScalarTransportFoam/blockCoupledSwirlTest/constant/polyMesh/boundary
index cf8b0f726..8c24ffa42 100644
--- a/tutorials/coupled/blockCoupledScalarTransportFoam/blockCoupledSwirlTest/constant/polyMesh/boundary
+++ b/tutorials/coupled/blockCoupledScalarTransportFoam/blockCoupledSwirlTest/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/coupled/blockCoupledScalarTransportFoam/blockCoupledSwirlTest/system/fvSchemes b/tutorials/coupled/blockCoupledScalarTransportFoam/blockCoupledSwirlTest/system/fvSchemes
index 4f86f5b77..2834dfa86 100644
--- a/tutorials/coupled/blockCoupledScalarTransportFoam/blockCoupledSwirlTest/system/fvSchemes
+++ b/tutorials/coupled/blockCoupledScalarTransportFoam/blockCoupledSwirlTest/system/fvSchemes
@@ -48,9 +48,4 @@ snGradSchemes
     default         corrected;
 }
 
-fluxRequired
-{
-    default         no;
-}
-
 // ************************************************************************* //
diff --git a/tutorials/coupled/conjugateHeatFoam/conjugateCavity/0/solid/materials b/tutorials/coupled/conjugateHeatFoam/conjugateCavity/0/solid/materials
index 8d96601d3..019e71299 100644
--- a/tutorials/coupled/conjugateHeatFoam/conjugateCavity/0/solid/materials
+++ b/tutorials/coupled/conjugateHeatFoam/conjugateCavity/0/solid/materials
@@ -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
 {
@@ -39,4 +39,5 @@ boundaryField
     }
 }
 
+
 // ************************************************************************* //
diff --git a/tutorials/coupled/conjugateHeatFoam/conjugateCavity/system/fvSchemes b/tutorials/coupled/conjugateHeatFoam/conjugateCavity/system/fvSchemes
index 414342c91..f7c603ec2 100644
--- a/tutorials/coupled/conjugateHeatFoam/conjugateCavity/system/fvSchemes
+++ b/tutorials/coupled/conjugateHeatFoam/conjugateCavity/system/fvSchemes
@@ -36,7 +36,7 @@ divSchemes
     // div(phi,R)      Gauss upwind;
     // div(R)          Gauss linear;
     // div(phi,nuTilda) Gauss upwind;
-    div((nuEff*dev(grad(U).T()))) Gauss linear;
+    div((nuEff*dev2(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes
@@ -62,10 +62,4 @@ snGradSchemes
     default         corrected;
 }
 
-fluxRequired
-{
-    default         no;
-    p_rgh;
-}
-
 // ************************************************************************* //
diff --git a/tutorials/coupled/conjugateHeatFoam/conjugateCavity/system/solid/fvSchemes b/tutorials/coupled/conjugateHeatFoam/conjugateCavity/system/solid/fvSchemes
index 8ae0a8c2d..fa14cec57 100644
--- a/tutorials/coupled/conjugateHeatFoam/conjugateCavity/system/solid/fvSchemes
+++ b/tutorials/coupled/conjugateHeatFoam/conjugateCavity/system/solid/fvSchemes
@@ -46,9 +46,4 @@ snGradSchemes
     default         uncorrected;
 }
 
-fluxRequired
-{
-    default         no;
-}
-
 // ************************************************************************* //
diff --git a/tutorials/coupled/conjugateHeatSimpleFoam/conjugateCavity/0/solid/materials b/tutorials/coupled/conjugateHeatSimpleFoam/conjugateCavity/0/solid/materials
index 8d96601d3..019e71299 100644
--- a/tutorials/coupled/conjugateHeatSimpleFoam/conjugateCavity/0/solid/materials
+++ b/tutorials/coupled/conjugateHeatSimpleFoam/conjugateCavity/0/solid/materials
@@ -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
 {
@@ -39,4 +39,5 @@ boundaryField
     }
 }
 
+
 // ************************************************************************* //
diff --git a/tutorials/coupled/conjugateHeatSimpleFoam/conjugateCavity/system/fvSchemes b/tutorials/coupled/conjugateHeatSimpleFoam/conjugateCavity/system/fvSchemes
index f18d60481..84fcec45f 100644
--- a/tutorials/coupled/conjugateHeatSimpleFoam/conjugateCavity/system/fvSchemes
+++ b/tutorials/coupled/conjugateHeatSimpleFoam/conjugateCavity/system/fvSchemes
@@ -36,7 +36,7 @@ divSchemes
     // div(phi,R)      Gauss upwind;
     // div(R)          Gauss linear;
     // div(phi,nuTilda) Gauss upwind;
-    div((nuEff*dev(grad(U).T()))) Gauss linear;
+    div((nuEff*dev2(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes
@@ -62,10 +62,4 @@ snGradSchemes
     default         corrected;
 }
 
-fluxRequired
-{
-    default         no;
-    p_rgh;
-}
-
 // ************************************************************************* //
diff --git a/tutorials/coupled/conjugateHeatSimpleFoam/conjugateCavity/system/solid/fvSchemes b/tutorials/coupled/conjugateHeatSimpleFoam/conjugateCavity/system/solid/fvSchemes
index 8ae0a8c2d..fa14cec57 100644
--- a/tutorials/coupled/conjugateHeatSimpleFoam/conjugateCavity/system/solid/fvSchemes
+++ b/tutorials/coupled/conjugateHeatSimpleFoam/conjugateCavity/system/solid/fvSchemes
@@ -46,9 +46,4 @@ snGradSchemes
     default         uncorrected;
 }
 
-fluxRequired
-{
-    default         no;
-}
-
 // ************************************************************************* //