diff --git a/applications/solvers/DNS/dnsFoam/createFields.H b/applications/solvers/DNS/dnsFoam/createFields.H
index 36fc513e0..34a2b963f 100644
--- a/applications/solvers/DNS/dnsFoam/createFields.H
+++ b/applications/solvers/DNS/dnsFoam/createFields.H
@@ -27,3 +27,5 @@
     );
 
 #   include "createPhi.H"
+
+    mesh.schemesDict().setFluxRequired(p.name());
diff --git a/applications/solvers/DNS/dnsFoam/dnsFoam.C b/applications/solvers/DNS/dnsFoam/dnsFoam.C
index 84e2f0fbb..6219abde1 100644
--- a/applications/solvers/DNS/dnsFoam/dnsFoam.C
+++ b/applications/solvers/DNS/dnsFoam/dnsFoam.C
@@ -35,6 +35,7 @@ Description
 #include "fft.H"
 #include "calcEk.H"
 #include "graph.H"
+#include "pisoControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -44,6 +45,9 @@ int main(int argc, char *argv[])
 
     #include "createTime.H"
     #include "createMeshNoClear.H"
+
+    pisoControl piso(mesh);
+
     #include "readTransportProperties.H"
     #include "createFields.H"
     #include "readTurbulenceProperties.H"
@@ -57,8 +61,6 @@ int main(int argc, char *argv[])
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        #include "readPISOControls.H"
-
         force.internalField() = ReImSum
         (
             fft::reverseTransform
@@ -83,7 +85,7 @@ int main(int argc, char *argv[])
 
         // --- PISO loop
 
-        for (int corr=1; corr<=1; corr++)
+        while (piso.correct())
         {
             volScalarField rUA = 1.0/UEqn.A();
 
diff --git a/applications/solvers/electromagnetics/mhdFoam/createFields.H b/applications/solvers/electromagnetics/mhdFoam/createFields.H
index b54cae44c..d5a0e62de 100644
--- a/applications/solvers/electromagnetics/mhdFoam/createFields.H
+++ b/applications/solvers/electromagnetics/mhdFoam/createFields.H
@@ -104,4 +104,7 @@
 
     label pRefCell = 0;
     scalar pRefValue = 0.0;
-    setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue);
+    setRefCell(p, piso.dict(), pRefCell, pRefValue);
+    mesh.schemesDict().setFluxRequired(p.name());
+
+    mesh.schemesDict().setFluxRequired(pB.name());
diff --git a/applications/solvers/electromagnetics/mhdFoam/mhdFoam.C b/applications/solvers/electromagnetics/mhdFoam/mhdFoam.C
index 44394f123..06058ce45 100644
--- a/applications/solvers/electromagnetics/mhdFoam/mhdFoam.C
+++ b/applications/solvers/electromagnetics/mhdFoam/mhdFoam.C
@@ -50,6 +50,7 @@ Description
 
 #include "fvCFD.H"
 #include "OSspecific.H"
+#include "pisoControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -59,6 +60,9 @@ int main(int argc, char *argv[])
 
 #   include "createTime.H"
 #   include "createMesh.H"
+
+    pisoControl piso(mesh);
+
 #   include "createFields.H"
 #   include "initContinuityErrs.H"
 
@@ -69,7 +73,6 @@ int main(int argc, char *argv[])
 
     while (runTime.loop())
     {
-#       include "readPISOControls.H"
 #       include "readBPISOControls.H"
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
@@ -91,7 +94,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();
 
@@ -100,7 +103,7 @@ int main(int argc, char *argv[])
                 phi = (fvc::interpolate(U) & mesh.Sf())
                     + fvc::ddtPhiCorr(rUA, U, phi);
 
-                for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+                while (piso.correctNonOrthogonal())
                 {
                     fvScalarMatrix pEqn
                     (
@@ -110,7 +113,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/equationReaderDemo/UEqn.H b/applications/solvers/equationReaderDemo/UEqn.H
index eb40b7665..e6fb6bf89 100644
--- a/applications/solvers/equationReaderDemo/UEqn.H
+++ b/applications/solvers/equationReaderDemo/UEqn.H
@@ -8,9 +8,4 @@
 
     UEqn().relax();
 
-    eqnResidual = solve
-    (
-        UEqn() == -fvc::grad(p)
-    ).initialResidual();
-
-    maxResidual = max(eqnResidual, maxResidual);
+    solve(UEqn() == -fvc::grad(p));
diff --git a/applications/solvers/equationReaderDemo/createFields.H b/applications/solvers/equationReaderDemo/createFields.H
index ab4249625..239a6938f 100644
--- a/applications/solvers/equationReaderDemo/createFields.H
+++ b/applications/solvers/equationReaderDemo/createFields.H
@@ -31,7 +31,8 @@
 
     label pRefCell = 0;
     scalar pRefValue = 0.0;
-    setRefCell(p, mesh.solutionDict().subDict("SIMPLE"), pRefCell, pRefValue);
+    setRefCell(p, simple.dict(), pRefCell, pRefValue);
+    mesh.schemesDict().setFluxRequired(p.name());
 
 
     singlePhaseTransportModel laminarTransport(U, phi);
diff --git a/applications/solvers/equationReaderDemo/equationReaderDemo.C b/applications/solvers/equationReaderDemo/equationReaderDemo.C
index 78cb4890b..d7efd06ab 100644
--- a/applications/solvers/equationReaderDemo/equationReaderDemo.C
+++ b/applications/solvers/equationReaderDemo/equationReaderDemo.C
@@ -36,6 +36,7 @@ Author
 #include "singlePhaseTransportModel.H"
 #include "RASModel.H"
 #include "IOEquationReader.H"
+#include "simpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -44,6 +45,9 @@ int main(int argc, char *argv[])
 #   include "setRootCase.H"
 #   include "createTime.H"
 #   include "createMesh.H"
+
+    simpleControl simple(mesh);
+
 #   include "createFields.H"
 #   include "createEquationReader.H"
 #   include "loadEquationData.H"
@@ -54,22 +58,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;
 
-    dictionary simple = mesh.solutionDict().subDict("SIMPLE");
-
-    int nNonOrthCorr =
-        simple.lookupOrDefault<int>("nNonOrthogonalCorrectors", 0);
-
-    bool momentumPredictor =
-        simple.lookupOrDefault<Switch>("momentumPredictor", true);
-
-#       include "initConvergenceCheck.H"
-
-        p.storePrevIter();
-
         // Pressure-velocity SIMPLE corrector
         {
 #           include "UEqn.H"
@@ -85,8 +77,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/equationReaderDemo/pEqn.H b/applications/solvers/equationReaderDemo/pEqn.H
index 42ebfc225..22921f652 100644
--- a/applications/solvers/equationReaderDemo/pEqn.H
+++ b/applications/solvers/equationReaderDemo/pEqn.H
@@ -7,7 +7,7 @@
     adjustPhi(phi, U, p);
 
     // Non-orthogonal pressure corrector loop
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    while (simple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
@@ -16,18 +16,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())
         {
             phi -= pEqn.flux();
         }
diff --git a/tutorials/DNS/dnsFoam/boxTurb16/constant/polyMesh/boundary b/tutorials/DNS/dnsFoam/boxTurb16/constant/polyMesh/boundary
index ef4afc88a..40c2ce70b 100644
--- a/tutorials/DNS/dnsFoam/boxTurb16/constant/polyMesh/boundary
+++ b/tutorials/DNS/dnsFoam/boxTurb16/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/DNS/dnsFoam/boxTurb16/system/fvSchemes b/tutorials/DNS/dnsFoam/boxTurb16/system/fvSchemes
index 128546641..07f21ed4c 100644
--- a/tutorials/DNS/dnsFoam/boxTurb16/system/fvSchemes
+++ b/tutorials/DNS/dnsFoam/boxTurb16/system/fvSchemes
@@ -50,10 +50,4 @@ snGradSchemes
     default         corrected;
 }
 
-fluxRequired
-{
-    default         no;
-    p;
-}
-
 // ************************************************************************* //
diff --git a/tutorials/electromagnetics/mhdFoam/hartmann/constant/polyMesh/boundary b/tutorials/electromagnetics/mhdFoam/hartmann/constant/polyMesh/boundary
index 769ac0f42..42bdf98f2 100644
--- a/tutorials/electromagnetics/mhdFoam/hartmann/constant/polyMesh/boundary
+++ b/tutorials/electromagnetics/mhdFoam/hartmann/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/electromagnetics/mhdFoam/hartmann/system/fvSchemes b/tutorials/electromagnetics/mhdFoam/hartmann/system/fvSchemes
index 03f2aba28..add92f8e6 100644
--- a/tutorials/electromagnetics/mhdFoam/hartmann/system/fvSchemes
+++ b/tutorials/electromagnetics/mhdFoam/hartmann/system/fvSchemes
@@ -56,11 +56,4 @@ snGradSchemes
     default         corrected;
 }
 
-fluxRequired
-{
-    default         no;
-    p;
-    pB;
-}
-
 // ************************************************************************* //
diff --git a/tutorials/equationReader/equationReaderDemo/pitzDaily/system/fvSchemes b/tutorials/equationReader/equationReaderDemo/pitzDaily/system/fvSchemes
index 7d773025d..e69ec2c98 100644
--- a/tutorials/equationReader/equationReaderDemo/pitzDaily/system/fvSchemes
+++ b/tutorials/equationReader/equationReaderDemo/pitzDaily/system/fvSchemes
@@ -35,8 +35,7 @@ divSchemes
     div(phi,R)      Gauss upwind;
     div(R)          Gauss linear;
     div(phi,nuTilda) Gauss upwind;
-    div((nuEff*dev(grad(U).T()))) Gauss linear; // backwards compatibility
-    div((nuEff*dev(T(grad(U))))) Gauss linear;
+    div((nuEff*dev2(T(grad(U))))) Gauss linear;
 }
 
 laplacianSchemes
@@ -61,10 +60,4 @@ snGradSchemes
     default         corrected;
 }
 
-fluxRequired
-{
-    default         no;
-    p;
-}
-
 // ************************************************************************* //
diff --git a/tutorials/multiSolver/multiSolverDemo/parTeeFitting2d/system/multiFvSchemes b/tutorials/multiSolver/multiSolverDemo/parTeeFitting2d/system/multiFvSchemes
index 137a5fc46..77740ea9f 100644
--- a/tutorials/multiSolver/multiSolverDemo/parTeeFitting2d/system/multiFvSchemes
+++ b/tutorials/multiSolver/multiSolverDemo/parTeeFitting2d/system/multiFvSchemes
@@ -53,12 +53,6 @@ multiSolver
         {
             default         corrected;
         }
-
-        fluxRequired
-        {
-            default         no;
-            p;
-        }
     }
 
     icoFoam2