diff --git a/applications/solvers/compressible/rhoPorousMRFPimpleFoam/UEqn.H b/applications/solvers/compressible/rhoPorousMRFPimpleFoam/UEqn.H
index ea80c2e7d..3e1b83048 100644
--- a/applications/solvers/compressible/rhoPorousMRFPimpleFoam/UEqn.H
+++ b/applications/solvers/compressible/rhoPorousMRFPimpleFoam/UEqn.H
@@ -7,37 +7,23 @@ tmp<fvVectorMatrix> UEqn
   + turbulence->divDevRhoReff(U)
 );
 
-if (oCorr == nOuterCorr-1)
-{
-    if (mesh.solutionDict().relax("UFinal"))
-    {
-        UEqn().relax(mesh.solutionDict().relaxationFactor("UFinal"));
-    }
-    else
-    {
-        UEqn().relax(1);
-    }
-}
-else
-{
-    UEqn().relax();
-}
+UEqn().relax
+(
+    mesh.solutionDict().relaxationFactor(U.select(pimple.finalIter()))
+);
 
 mrfZones.addCoriolis(rho, UEqn());
 pZones.addResistance(UEqn());
 
 volScalarField rUA = 1.0/UEqn().A();
 
-if (momentumPredictor)
+if (pimple.momentumPredictor())
 {
-    if (oCorr == nOuterCorr-1)
-    {
-        solve(UEqn() == -fvc::grad(p), mesh.solutionDict().solver("UFinal"));
-    }
-    else
-    {
-        solve(UEqn() == -fvc::grad(p));
-    }
+    solve
+    (
+        UEqn() == -fvc::grad(p),
+        mesh.solutionDict().solver((U.select(pimple.finalIter())))
+    );
 }
 else
 {
diff --git a/applications/solvers/compressible/rhoPorousMRFPimpleFoam/createFields.H b/applications/solvers/compressible/rhoPorousMRFPimpleFoam/createFields.H
index b9a86ef99..ee3beea43 100644
--- a/applications/solvers/compressible/rhoPorousMRFPimpleFoam/createFields.H
+++ b/applications/solvers/compressible/rhoPorousMRFPimpleFoam/createFields.H
@@ -68,3 +68,5 @@
 
     porousZones pZones(mesh);
     Switch pressureImplicitPorosity(false);
+
+    mesh.schemesDict().setFluxRequired(p.name());
diff --git a/applications/solvers/compressible/rhoPorousMRFPimpleFoam/pEqn.H b/applications/solvers/compressible/rhoPorousMRFPimpleFoam/pEqn.H
index b714af2df..920ba0e5d 100644
--- a/applications/solvers/compressible/rhoPorousMRFPimpleFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoPorousMRFPimpleFoam/pEqn.H
@@ -3,12 +3,12 @@ rho = thermo.rho();
 volScalarField rUA = 1.0/UEqn().A();
 U = rUA*UEqn().H();
 
-if (nCorr <= 1)
+if (pimple.nCorrPISO() <= 1)
 {
     UEqn.clear();
 }
 
-if (transonic)
+if (pimple.transonic())
 {
     surfaceScalarField phid
     (
@@ -21,7 +21,7 @@ if (transonic)
     );
     mrfZones.relativeFlux(fvc::interpolate(psi), phid);
 
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
@@ -30,21 +30,12 @@ if (transonic)
           - fvm::laplacian(rho*rUA, p)
         );
 
-        if
+        pEqn.solve
         (
-            oCorr == nOuterCorr-1
-            && corr == nCorr-1
-            && nonOrth == nNonOrthCorr
-        )
-        {
-            pEqn.solve(mesh.solutionDict().solver("pFinal"));
-        }
-        else
-        {
-            pEqn.solve();
-        }
+            mesh.solutionDict().solver(p.select(pimple.finalInnerIter()))
+        );
 
-        if (nonOrth == nNonOrthCorr)
+        if (pimple.finalNonOrthogonalIter())
         {
             phi == pEqn.flux();
         }
@@ -60,7 +51,7 @@ else
         );
     mrfZones.relativeFlux(fvc::interpolate(rho), phi);
 
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    while (pimple.correctNonOrthogonal())
     {
         // Pressure corrector
         fvScalarMatrix pEqn
@@ -70,21 +61,12 @@ else
           - fvm::laplacian(rho*rUA, p)
         );
 
-        if
+        pEqn.solve
         (
-            oCorr == nOuterCorr-1
-         && corr == nCorr-1
-         && nonOrth == nNonOrthCorr
-        )
-        {
-            pEqn.solve(mesh.solutionDict().solver("pFinal"));
-        }
-        else
-        {
-            pEqn.solve();
-        }
+            mesh.solutionDict().solver(p.select(pimple.finalInnerIter()))
+        );
 
-        if (nonOrth == nNonOrthCorr)
+        if (pimple.finalNonOrthogonalIter())
         {
             phi += pEqn.flux();
         }
diff --git a/applications/solvers/compressible/rhoPorousMRFPimpleFoam/rhoPorousMRFPimpleFoam.C b/applications/solvers/compressible/rhoPorousMRFPimpleFoam/rhoPorousMRFPimpleFoam.C
index a097479ae..94c3c077c 100644
--- a/applications/solvers/compressible/rhoPorousMRFPimpleFoam/rhoPorousMRFPimpleFoam.C
+++ b/applications/solvers/compressible/rhoPorousMRFPimpleFoam/rhoPorousMRFPimpleFoam.C
@@ -39,6 +39,7 @@ Description
 #include "bound.H"
 #include "MRFZones.H"
 #include "porousZones.H"
+#include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -47,6 +48,9 @@ int main(int argc, char *argv[])
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createMesh.H"
+
+    pimpleControl pimple(mesh);
+
     #include "createFields.H"
     #include "initContinuityErrs.H"
     #include "createTimeControls.H"
@@ -58,7 +62,6 @@ int main(int argc, char *argv[])
     while (runTime.run())
     {
         #include "readTimeControls.H"
-        #include "readPIMPLEControls.H"
         #include "compressibleCourantNo.H"
         #include "setDeltaT.H"
 
@@ -66,7 +69,7 @@ int main(int argc, char *argv[])
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        if (nOuterCorr != 1)
+        if (!pimple.firstIter())
         {
             p.storePrevIter();
             rho.storePrevIter();
@@ -75,13 +78,13 @@ int main(int argc, char *argv[])
         #include "rhoEqn.H"
 
         // --- Pressure-velocity PIMPLE corrector loop
-        for (int oCorr=0; oCorr<nOuterCorr; oCorr++)
+        while (pimple.loop())
         {
             #include "UEqn.H"
             #include "hEqn.H"
 
             // --- PISO loop
-            for (int corr = 0; corr < nCorr; corr++)
+            while (pimple.correct())
             {
                 #include "pEqn.H"
             }
diff --git a/applications/solvers/compressible/rhoPorousSimpleFoam/UEqn.H b/applications/solvers/compressible/rhoPorousSimpleFoam/UEqn.H
index 546094428..fd79c8400 100644
--- a/applications/solvers/compressible/rhoPorousSimpleFoam/UEqn.H
+++ b/applications/solvers/compressible/rhoPorousSimpleFoam/UEqn.H
@@ -34,12 +34,7 @@
     {
         pZones.addResistance(UEqn());
 
-        eqnResidual = solve
-        (
-            UEqn() == -fvc::grad(p)
-        ). initialResidual();
-
-        maxResidual = max(eqnResidual, maxResidual);
+        solve(UEqn() == -fvc::grad(p));
 
         trAU = 1.0/UEqn().A();
         trAU().rename("rAU");
diff --git a/applications/solvers/compressible/rhoPorousSimpleFoam/convergenceCheck.H b/applications/solvers/compressible/rhoPorousSimpleFoam/convergenceCheck.H
deleted file mode 100644
index 895806319..000000000
--- a/applications/solvers/compressible/rhoPorousSimpleFoam/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/compressible/rhoPorousSimpleFoam/createFields.H b/applications/solvers/compressible/rhoPorousSimpleFoam/createFields.H
index 1177cba2a..0e1ba2436 100644
--- a/applications/solvers/compressible/rhoPorousSimpleFoam/createFields.H
+++ b/applications/solvers/compressible/rhoPorousSimpleFoam/createFields.H
@@ -44,6 +44,7 @@
     label pRefCell = 0;
     scalar pRefValue = 0.0;
     setRefCell(p, mesh.solutionDict().subDict("SIMPLE"), pRefCell, pRefValue);
+    mesh.schemesDict().setFluxRequired(p.name());
 
     dimensionedScalar pMin
     (
diff --git a/applications/solvers/compressible/rhoPorousSimpleFoam/hEqn.H b/applications/solvers/compressible/rhoPorousSimpleFoam/hEqn.H
index 605b8820d..324dcd642 100644
--- a/applications/solvers/compressible/rhoPorousSimpleFoam/hEqn.H
+++ b/applications/solvers/compressible/rhoPorousSimpleFoam/hEqn.H
@@ -11,8 +11,7 @@
 
     hEqn.relax();
 
-    eqnResidual = hEqn.solve().initialResidual();
-    maxResidual = max(eqnResidual, maxResidual);
+    hEqn.solve().initialResidual();
 
     thermo.correct();
 }
diff --git a/applications/solvers/compressible/rhoPorousSimpleFoam/initConvergenceCheck.H b/applications/solvers/compressible/rhoPorousSimpleFoam/initConvergenceCheck.H
deleted file mode 100644
index b56197f22..000000000
--- a/applications/solvers/compressible/rhoPorousSimpleFoam/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/compressible/rhoPorousSimpleFoam/pEqn.H b/applications/solvers/compressible/rhoPorousSimpleFoam/pEqn.H
index 4d8e010f7..ab81dac30 100644
--- a/applications/solvers/compressible/rhoPorousSimpleFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoPorousSimpleFoam/pEqn.H
@@ -12,7 +12,7 @@ UEqn.clear();
 phi = fvc::interpolate(rho*U) & mesh.Sf();
 bool closedVolume = adjustPhi(phi, U, p);
 
-for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+while (simple.correctNonOrthogonal())
 {
     tmp<fvScalarMatrix> tpEqn;
 
@@ -26,18 +26,10 @@ for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     }
 
     tpEqn().setReference(pRefCell, pRefValue);
-    // retain the residual from the first iteration
-    if (nonOrth == 0)
-    {
-        eqnResidual = tpEqn().solve().initialResidual();
-        maxResidual = max(eqnResidual, maxResidual);
-    }
-    else
-    {
-        tpEqn().solve();
-    }
 
-    if (nonOrth == nNonOrthCorr)
+    tpEqn().solve();
+
+    if (simple.finalNonOrthogonalIter())
     {
         phi -= tpEqn().flux();
     }
diff --git a/applications/solvers/compressible/rhoPorousSimpleFoam/rhoPorousSimpleFoam.C b/applications/solvers/compressible/rhoPorousSimpleFoam/rhoPorousSimpleFoam.C
index b72e7efd2..817c33ad8 100644
--- a/applications/solvers/compressible/rhoPorousSimpleFoam/rhoPorousSimpleFoam.C
+++ b/applications/solvers/compressible/rhoPorousSimpleFoam/rhoPorousSimpleFoam.C
@@ -34,6 +34,7 @@ Description
 #include "basicPsiThermo.H"
 #include "RASModel.H"
 #include "porousZones.H"
+#include "simpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -42,6 +43,9 @@ int main(int argc, char *argv[])
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createMesh.H"
+
+    simpleControl simple(mesh);
+
     #include "createFields.H"
     #include "initContinuityErrs.H"
 
@@ -53,9 +57,6 @@ int main(int argc, char *argv[])
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        #include "readSIMPLEControls.H"
-        #include "initConvergenceCheck.H"
-
         p.storePrevIter();
         rho.storePrevIter();
 
@@ -72,8 +73,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/compressible/rhoSimpleFoam/UEqn.H b/applications/solvers/compressible/rhoSimpleFoam/UEqn.H
index 401e1203d..21ec2646b 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/UEqn.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/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/compressible/rhoSimpleFoam/convergenceCheck.H b/applications/solvers/compressible/rhoSimpleFoam/convergenceCheck.H
deleted file mode 100644
index 895806319..000000000
--- a/applications/solvers/compressible/rhoSimpleFoam/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/compressible/rhoSimpleFoam/createFields.H b/applications/solvers/compressible/rhoSimpleFoam/createFields.H
index 122280cfa..38aaa9e5c 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/createFields.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/createFields.H
@@ -43,6 +43,7 @@
     label pRefCell = 0;
     scalar pRefValue = 0.0;
     setRefCell(p, mesh.solutionDict().subDict("SIMPLE"), pRefCell, pRefValue);
+    mesh.schemesDict().setFluxRequired(p.name());
 
     dimensionedScalar rhoMax
     (
diff --git a/applications/solvers/compressible/rhoSimpleFoam/hEqn.H b/applications/solvers/compressible/rhoSimpleFoam/hEqn.H
index 57395e977..e6c74938e 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/hEqn.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/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/compressible/rhoSimpleFoam/initConvergenceCheck.H b/applications/solvers/compressible/rhoSimpleFoam/initConvergenceCheck.H
deleted file mode 100644
index b56197f22..000000000
--- a/applications/solvers/compressible/rhoSimpleFoam/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/compressible/rhoSimpleFoam/pEqn.H b/applications/solvers/compressible/rhoSimpleFoam/pEqn.H
index 04af6df75..3865aeb9f 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/pEqn.H
+++ b/applications/solvers/compressible/rhoSimpleFoam/pEqn.H
@@ -9,7 +9,7 @@ UEqn.clear();
 
 bool closedVolume = false;
 
-if (transonic)
+if (simple.transonic())
 {
     surfaceScalarField phid
     (
@@ -17,7 +17,7 @@ if (transonic)
         fvc::interpolate(psi)*(fvc::interpolate(U) & mesh.Sf())
     );
 
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    while (simple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
@@ -26,22 +26,13 @@ if (transonic)
         );
 
         // Relax the pressure equation to ensure diagonal-dominance
-        pEqn.relax(mesh.solutionDict().relaxationFactor("pEqn"));
+        pEqn.relax();
 
         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();
         }
@@ -52,7 +43,7 @@ else
     phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf());
     closedVolume = adjustPhi(phi, U, p);
 
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    while (simple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
@@ -61,18 +52,9 @@ else
 
         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/applications/solvers/compressible/rhoSimpleFoam/rhoSimpleFoam.C b/applications/solvers/compressible/rhoSimpleFoam/rhoSimpleFoam.C
index 6c63d09c0..189bd4db2 100644
--- a/applications/solvers/compressible/rhoSimpleFoam/rhoSimpleFoam.C
+++ b/applications/solvers/compressible/rhoSimpleFoam/rhoSimpleFoam.C
@@ -33,6 +33,7 @@ Description
 #include "fvCFD.H"
 #include "basicPsiThermo.H"
 #include "RASModel.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 "createFields.H"
     #include "initContinuityErrs.H"
 
@@ -52,9 +56,6 @@ int main(int argc, char *argv[])
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-#       include "readSIMPLEControls.H"
-#       include "initConvergenceCheck.H"
-
         p.storePrevIter();
         rho.storePrevIter();
 
@@ -72,8 +73,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/compressible/rhoSonicFoam/createFields.H b/applications/solvers/compressible/rhoSonicFoam/createFields.H
index f0c0e4a22..aa121b96e 100644
--- a/applications/solvers/compressible/rhoSonicFoam/createFields.H
+++ b/applications/solvers/compressible/rhoSonicFoam/createFields.H
@@ -94,3 +94,5 @@
         rho*Cv*T + 0.5*rho*magSqr(rhoU/rho),
         T.boundaryField().types()
     );
+
+    mesh.schemesDict().setFluxRequired(p.name());
diff --git a/applications/solvers/compressible/rhopSonicFoam/createFields.H b/applications/solvers/compressible/rhopSonicFoam/createFields.H
index fed843271..c8e478e88 100644
--- a/applications/solvers/compressible/rhopSonicFoam/createFields.H
+++ b/applications/solvers/compressible/rhopSonicFoam/createFields.H
@@ -127,3 +127,4 @@
     fields.add(magRhoU);
     fields.add(H);
 
+    mesh.schemesDict().setFluxRequired(p.name());
diff --git a/applications/solvers/compressible/rhopSonicFoam/rhopSonicFoam.C b/applications/solvers/compressible/rhopSonicFoam/rhopSonicFoam.C
index 5d1e32170..58d362fec 100644
--- a/applications/solvers/compressible/rhopSonicFoam/rhopSonicFoam.C
+++ b/applications/solvers/compressible/rhopSonicFoam/rhopSonicFoam.C
@@ -36,6 +36,7 @@ Description
 #include "MUSCL.H"
 #include "LimitedScheme.H"
 #include "boundaryTypes.H"
+#include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -45,6 +46,9 @@ int main(int argc, char *argv[])
 #   include "setRootCase.H"
 #   include "createTime.H"
 #   include "createMesh.H"
+
+    pimpleControl pimple(mesh);
+
 #   include "readThermodynamicProperties.H"
 #   include "createFields.H"
 #   include "createTimeControls.H"
@@ -57,8 +61,11 @@ int main(int argc, char *argv[])
     {
         Info<< "Time = " << runTime.value() << nl << endl;
 
-#       include "readPISOControls.H"
-        scalar HbyAblend = readScalar(piso.lookup("HbyAblend"));
+        scalar HbyAblend =
+            readScalar
+            (
+                mesh.solutionDict().subDict("PIMPLE").lookup("HbyAblend")
+            );
 
 #       include "readTimeControls.H"
 
@@ -72,7 +79,7 @@ int main(int argc, char *argv[])
 
 #       include "setDeltaT.H"
 
-        for (int outerCorr = 0; outerCorr < nOuterCorr; outerCorr++)
+        while (pimple.loop())
         {
             magRhoU = mag(rhoU);
             H = (rhoE + p)/rho;
@@ -118,7 +125,7 @@ int main(int argc, char *argv[])
             psi = 1.0/(R*T);
             p = rho/psi;
 
-            for (int corr = 0; corr < nCorr; corr++)
+            while (pimple.correct())
             {
                 volScalarField rrhoUA = 1.0/rhoUEqn.A();
                 surfaceScalarField rrhoUAf("rrhoUAf", fvc::interpolate(rrhoUA));
diff --git a/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/constant/polyMesh/blockMeshDict b/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/constant/polyMesh/blockMeshDict
index afb2d3392..9664dc64d 100644
--- a/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/constant/polyMesh/blockMeshDict
+++ b/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/constant/polyMesh/blockMeshDict
@@ -15,148 +15,409 @@ FoamFile
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 // General macros to create 2D/extruded-2D meshes
 
+
+
+
+
+
+
+
+
+
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 convertToMeters 0.1;
 
 // Hub radius
 
+
 // Impeller-tip radius
 
+
 // Baffle-tip radius
 
+
 // Tank radius
 
+
 // MRF region radius
 
+
 // Thickness of 2D slab
 
+
 // Base z
 
+
 // Top z
 
+
 // Number of cells radially between hub and impeller tip
 
+
 // Number of cells radially in each of the two regions between
 // impeller and baffle tips
 
+
 // Number of cells radially between baffle tip and tank
 
+
 // Number of cells azimuthally in each of the 8 blocks
 
+
 // Number of cells in the thickness of the slab
 
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 vertices
 (
-    (0.2 0 0) // Vertex r0b = 0
-    (0.2 0 0) // Vertex r0sb = 1
-    (0.141421356364228 -0.141421356110391 0) // Vertex r1b = 2
-    (3.58979347393082e-10 -0.2 0) // Vertex r2b = 3
-    (3.58979347393082e-10 -0.2 0) // Vertex r2sb = 4
-    (-0.141421355856554 -0.141421356618065 0) // Vertex r3b = 5
-    (-0.2 7.17958694786164e-10 0) // Vertex r4b = 6
-    (-0.2 7.17958694786164e-10 0) // Vertex r4sb = 7
-    (-0.141421355856554 0.141421356618065 0) // Vertex r5b = 8
-    (3.58979347393082e-10 0.2 0) // Vertex r6b = 9
-    (3.58979347393082e-10 0.2 0) // Vertex r6sb = 10
-    (0.141421356364228 0.141421356110391 0) // Vertex r7b = 11
+    (0.2 0 0) // Vertex r0b = 0 
+    (0.2 0 0) // Vertex r0sb = 1 
+    (0.141421356364228 -0.141421356110391 0) // Vertex r1b = 2 
+    (3.58979347393082e-10 -0.2 0) // Vertex r2b = 3 
+    (3.58979347393082e-10 -0.2 0) // Vertex r2sb = 4 
+    (-0.141421355856554 -0.141421356618065 0) // Vertex r3b = 5 
+    (-0.2 7.17958694786164e-10 0) // Vertex r4b = 6 
+    (-0.2 7.17958694786164e-10 0) // Vertex r4sb = 7 
+    (-0.141421355856554 0.141421356618065 0) // Vertex r5b = 8 
+    (3.58979347393082e-10 0.2 0) // Vertex r6b = 9 
+    (3.58979347393082e-10 0.2 0) // Vertex r6sb = 10 
+    (0.141421356364228 0.141421356110391 0) // Vertex r7b = 11 
 
-    (0.5 0 0) // Vertex rb0b = 12
-    (0.353553390910569 -0.353553390275978 0) // Vertex rb1b = 13
-    (8.97448368482705e-10 -0.5 0) // Vertex rb2b = 14
-    (-0.353553389641386 -0.353553391545162 0) // Vertex rb3b = 15
-    (-0.5 1.79489673696541e-09 0) // Vertex rb4b = 16
-    (-0.353553389641386 0.353553391545162 0) // Vertex rb5b = 17
-    (8.97448368482705e-10 0.5 0) // Vertex rb6b = 18
-    (0.353553390910569 0.353553390275978 0) // Vertex rb7b = 19
+    (0.5 0 0) // Vertex rb0b = 12 
+    (0.353553390910569 -0.353553390275978 0) // Vertex rb1b = 13 
+    (8.97448368482705e-10 -0.5 0) // Vertex rb2b = 14 
+    (-0.353553389641386 -0.353553391545162 0) // Vertex rb3b = 15 
+    (-0.5 1.79489673696541e-09 0) // Vertex rb4b = 16 
+    (-0.353553389641386 0.353553391545162 0) // Vertex rb5b = 17 
+    (8.97448368482705e-10 0.5 0) // Vertex rb6b = 18 
+    (0.353553390910569 0.353553390275978 0) // Vertex rb7b = 19 
 
-    (0.6 0 0) // Vertex ri0b = 20
-    (0.424264069092683 -0.424264068331174 0) // Vertex ri1b = 21
-    (1.07693804217925e-09 -0.6 0) // Vertex ri2b = 22
-    (-0.424264067569663 -0.424264069854194 0) // Vertex ri3b = 23
-    (-0.6 2.15387608435849e-09 0) // Vertex ri4b = 24
-    (-0.424264067569663 0.424264069854194 0) // Vertex ri5b = 25
-    (1.07693804217925e-09 0.6 0) // Vertex ri6b = 26
-    (0.424264069092683 0.424264068331174 0) // Vertex ri7b = 27
+    (0.6 0 0) // Vertex ri0b = 20 
+    (0.424264069092683 -0.424264068331174 0) // Vertex ri1b = 21 
+    (1.07693804217925e-09 -0.6 0) // Vertex ri2b = 22 
+    (-0.424264067569663 -0.424264069854194 0) // Vertex ri3b = 23 
+    (-0.6 2.15387608435849e-09 0) // Vertex ri4b = 24 
+    (-0.424264067569663 0.424264069854194 0) // Vertex ri5b = 25 
+    (1.07693804217925e-09 0.6 0) // Vertex ri6b = 26 
+    (0.424264069092683 0.424264068331174 0) // Vertex ri7b = 27 
 
-    (0.7 0 0) // Vertex Rb0b = 28
-    (0.494974747274797 -0.494974746386369 0) // Vertex Rb1b = 29
-    (1.25642771587579e-09 -0.7 0) // Vertex Rb2b = 30
-    (-0.49497474549794 -0.494974748163226 0) // Vertex Rb3b = 31
-    (-0.7 2.51285543175157e-09 0) // Vertex Rb4b = 32
-    (-0.49497474549794 0.494974748163226 0) // Vertex Rb5b = 33
-    (1.25642771587579e-09 0.7 0) // Vertex Rb6b = 34
-    (0.494974747274797 0.494974746386369 0) // Vertex Rb7b = 35
+    (0.7 0 0) // Vertex Rb0b = 28 
+    (0.494974747274797 -0.494974746386369 0) // Vertex Rb1b = 29 
+    (1.25642771587579e-09 -0.7 0) // Vertex Rb2b = 30 
+    (-0.49497474549794 -0.494974748163226 0) // Vertex Rb3b = 31 
+    (-0.7 2.51285543175157e-09 0) // Vertex Rb4b = 32 
+    (-0.49497474549794 0.494974748163226 0) // Vertex Rb5b = 33 
+    (1.25642771587579e-09 0.7 0) // Vertex Rb6b = 34 
+    (0.494974747274797 0.494974746386369 0) // Vertex Rb7b = 35 
 
-    (1 0 0) // Vertex R0b = 36
-    (0.707106781821139 -0.707106780551956 0) // Vertex R1b = 37
-    (0.707106781821139 -0.707106780551956 0) // Vertex R1sb = 38
-    (1.79489673696541e-09 -1 0) // Vertex R2b = 39
-    (-0.707106779282772 -0.707106783090323 0) // Vertex R3b = 40
-    (-0.707106779282772 -0.707106783090323 0) // Vertex R3sb = 41
-    (-1 3.58979347393082e-09 0) // Vertex R4b = 42
-    (-0.707106779282772 0.707106783090323 0) // Vertex R5b = 43
-    (-0.707106779282772 0.707106783090323 0) // Vertex R5sb = 44
-    (1.79489673696541e-09 1 0) // Vertex R6b = 45
-    (0.707106781821139 0.707106780551956 0) // Vertex R7b = 46
-    (0.707106781821139 0.707106780551956 0) // Vertex R7sb = 47
+    (1 0 0) // Vertex R0b = 36 
+    (0.707106781821139 -0.707106780551956 0) // Vertex R1b = 37 
+    (0.707106781821139 -0.707106780551956 0) // Vertex R1sb = 38 
+    (1.79489673696541e-09 -1 0) // Vertex R2b = 39 
+    (-0.707106779282772 -0.707106783090323 0) // Vertex R3b = 40 
+    (-0.707106779282772 -0.707106783090323 0) // Vertex R3sb = 41 
+    (-1 3.58979347393082e-09 0) // Vertex R4b = 42 
+    (-0.707106779282772 0.707106783090323 0) // Vertex R5b = 43 
+    (-0.707106779282772 0.707106783090323 0) // Vertex R5sb = 44 
+    (1.79489673696541e-09 1 0) // Vertex R6b = 45 
+    (0.707106781821139 0.707106780551956 0) // Vertex R7b = 46 
+    (0.707106781821139 0.707106780551956 0) // Vertex R7sb = 47 
 
-    (0.2 0 0.1) // Vertex r0t = 48
-    (0.2 0 0.1) // Vertex r0st = 49
-    (0.141421356364228 -0.141421356110391 0.1) // Vertex r1t = 50
-    (3.58979347393082e-10 -0.2 0.1) // Vertex r2t = 51
-    (3.58979347393082e-10 -0.2 0.1) // Vertex r2st = 52
-    (-0.141421355856554 -0.141421356618065 0.1) // Vertex r3t = 53
-    (-0.2 7.17958694786164e-10 0.1) // Vertex r4t = 54
-    (-0.2 7.17958694786164e-10 0.1) // Vertex r4st = 55
-    (-0.141421355856554 0.141421356618065 0.1) // Vertex r5t = 56
-    (3.58979347393082e-10 0.2 0.1) // Vertex r6t = 57
-    (3.58979347393082e-10 0.2 0.1) // Vertex r6st = 58
-    (0.141421356364228 0.141421356110391 0.1) // Vertex r7t = 59
+    (0.2 0 0.1) // Vertex r0t = 48 
+    (0.2 0 0.1) // Vertex r0st = 49 
+    (0.141421356364228 -0.141421356110391 0.1) // Vertex r1t = 50 
+    (3.58979347393082e-10 -0.2 0.1) // Vertex r2t = 51 
+    (3.58979347393082e-10 -0.2 0.1) // Vertex r2st = 52 
+    (-0.141421355856554 -0.141421356618065 0.1) // Vertex r3t = 53 
+    (-0.2 7.17958694786164e-10 0.1) // Vertex r4t = 54 
+    (-0.2 7.17958694786164e-10 0.1) // Vertex r4st = 55 
+    (-0.141421355856554 0.141421356618065 0.1) // Vertex r5t = 56 
+    (3.58979347393082e-10 0.2 0.1) // Vertex r6t = 57 
+    (3.58979347393082e-10 0.2 0.1) // Vertex r6st = 58 
+    (0.141421356364228 0.141421356110391 0.1) // Vertex r7t = 59 
 
-    (0.5 0 0.1) // Vertex rb0t = 60
-    (0.353553390910569 -0.353553390275978 0.1) // Vertex rb1t = 61
-    (8.97448368482705e-10 -0.5 0.1) // Vertex rb2t = 62
-    (-0.353553389641386 -0.353553391545162 0.1) // Vertex rb3t = 63
-    (-0.5 1.79489673696541e-09 0.1) // Vertex rb4t = 64
-    (-0.353553389641386 0.353553391545162 0.1) // Vertex rb5t = 65
-    (8.97448368482705e-10 0.5 0.1) // Vertex rb6t = 66
-    (0.353553390910569 0.353553390275978 0.1) // Vertex rb7t = 67
+    (0.5 0 0.1) // Vertex rb0t = 60 
+    (0.353553390910569 -0.353553390275978 0.1) // Vertex rb1t = 61 
+    (8.97448368482705e-10 -0.5 0.1) // Vertex rb2t = 62 
+    (-0.353553389641386 -0.353553391545162 0.1) // Vertex rb3t = 63 
+    (-0.5 1.79489673696541e-09 0.1) // Vertex rb4t = 64 
+    (-0.353553389641386 0.353553391545162 0.1) // Vertex rb5t = 65 
+    (8.97448368482705e-10 0.5 0.1) // Vertex rb6t = 66 
+    (0.353553390910569 0.353553390275978 0.1) // Vertex rb7t = 67 
 
-    (0.6 0 0.1) // Vertex ri0t = 68
-    (0.424264069092683 -0.424264068331174 0.1) // Vertex ri1t = 69
-    (1.07693804217925e-09 -0.6 0.1) // Vertex ri2t = 70
-    (-0.424264067569663 -0.424264069854194 0.1) // Vertex ri3t = 71
-    (-0.6 2.15387608435849e-09 0.1) // Vertex ri4t = 72
-    (-0.424264067569663 0.424264069854194 0.1) // Vertex ri5t = 73
-    (1.07693804217925e-09 0.6 0.1) // Vertex ri6t = 74
-    (0.424264069092683 0.424264068331174 0.1) // Vertex ri7t = 75
+    (0.6 0 0.1) // Vertex ri0t = 68 
+    (0.424264069092683 -0.424264068331174 0.1) // Vertex ri1t = 69 
+    (1.07693804217925e-09 -0.6 0.1) // Vertex ri2t = 70 
+    (-0.424264067569663 -0.424264069854194 0.1) // Vertex ri3t = 71 
+    (-0.6 2.15387608435849e-09 0.1) // Vertex ri4t = 72 
+    (-0.424264067569663 0.424264069854194 0.1) // Vertex ri5t = 73 
+    (1.07693804217925e-09 0.6 0.1) // Vertex ri6t = 74 
+    (0.424264069092683 0.424264068331174 0.1) // Vertex ri7t = 75 
 
-    (0.7 0 0.1) // Vertex Rb0t = 76
-    (0.494974747274797 -0.494974746386369 0.1) // Vertex Rb1t = 77
-    (1.25642771587579e-09 -0.7 0.1) // Vertex Rb2t = 78
-    (-0.49497474549794 -0.494974748163226 0.1) // Vertex Rb3t = 79
-    (-0.7 2.51285543175157e-09 0.1) // Vertex Rb4t = 80
-    (-0.49497474549794 0.494974748163226 0.1) // Vertex Rb5t = 81
-    (1.25642771587579e-09 0.7 0.1) // Vertex Rb6t = 82
-    (0.494974747274797 0.494974746386369 0.1) // Vertex Rb7t = 83
+    (0.7 0 0.1) // Vertex Rb0t = 76 
+    (0.494974747274797 -0.494974746386369 0.1) // Vertex Rb1t = 77 
+    (1.25642771587579e-09 -0.7 0.1) // Vertex Rb2t = 78 
+    (-0.49497474549794 -0.494974748163226 0.1) // Vertex Rb3t = 79 
+    (-0.7 2.51285543175157e-09 0.1) // Vertex Rb4t = 80 
+    (-0.49497474549794 0.494974748163226 0.1) // Vertex Rb5t = 81 
+    (1.25642771587579e-09 0.7 0.1) // Vertex Rb6t = 82 
+    (0.494974747274797 0.494974746386369 0.1) // Vertex Rb7t = 83 
 
-    (1 0 0.1) // Vertex R0t = 84
-    (0.707106781821139 -0.707106780551956 0.1) // Vertex R1t = 85
-    (0.707106781821139 -0.707106780551956 0.1) // Vertex R1st = 86
-    (1.79489673696541e-09 -1 0.1) // Vertex R2t = 87
-    (-0.707106779282772 -0.707106783090323 0.1) // Vertex R3t = 88
-    (-0.707106779282772 -0.707106783090323 0.1) // Vertex R3st = 89
-    (-1 3.58979347393082e-09 0.1) // Vertex R4t = 90
-    (-0.707106779282772 0.707106783090323 0.1) // Vertex R5t = 91
-    (-0.707106779282772 0.707106783090323 0.1) // Vertex R5st = 92
-    (1.79489673696541e-09 1 0.1) // Vertex R6t = 93
-    (0.707106781821139 0.707106780551956 0.1) // Vertex R7t = 94
-    (0.707106781821139 0.707106780551956 0.1) // Vertex R7st = 95
+    (1 0 0.1) // Vertex R0t = 84 
+    (0.707106781821139 -0.707106780551956 0.1) // Vertex R1t = 85 
+    (0.707106781821139 -0.707106780551956 0.1) // Vertex R1st = 86 
+    (1.79489673696541e-09 -1 0.1) // Vertex R2t = 87 
+    (-0.707106779282772 -0.707106783090323 0.1) // Vertex R3t = 88 
+    (-0.707106779282772 -0.707106783090323 0.1) // Vertex R3st = 89 
+    (-1 3.58979347393082e-09 0.1) // Vertex R4t = 90 
+    (-0.707106779282772 0.707106783090323 0.1) // Vertex R5t = 91 
+    (-0.707106779282772 0.707106783090323 0.1) // Vertex R5st = 92 
+    (1.79489673696541e-09 1 0.1) // Vertex R6t = 93 
+    (0.707106781821139 0.707106780551956 0.1) // Vertex R7t = 94 
+    (0.707106781821139 0.707106780551956 0.1) // Vertex R7st = 95 
 );
 
 blocks
diff --git a/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/constant/polyMesh/boundary b/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/constant/polyMesh/boundary
index 02383fd43..c42c32d0c 100644
--- a/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/constant/polyMesh/boundary
+++ b/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/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/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/system/fvSchemes b/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/system/fvSchemes
index a8cc61ad4..427aa072d 100644
--- a/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/system/fvSchemes
+++ b/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/system/fvSchemes
@@ -32,7 +32,7 @@ divSchemes
     div(phi,h)      Gauss limitedLinear 1;
     div(phi,k)      Gauss limitedLinear 1;
     div(phi,epsilon) Gauss limitedLinear 1;
-    div((muEff*dev2(grad(U).T()))) Gauss linear;
+    div((muEff*dev2(T(grad(U))))) Gauss linear;
     div(phiU,p) Gauss linear;
 }
 
@@ -57,10 +57,4 @@ snGradSchemes
     default         corrected;
 }
 
-fluxRequired
-{
-    default         no;
-    p               ;
-}
-
 // ************************************************************************* //
diff --git a/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/system/fvSolution b/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/system/fvSolution
index a62d3ed4d..6d60a61a0 100644
--- a/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/system/fvSolution
+++ b/tutorials/compressible/rhoPorousMRFPimpleFoam/mixerVessel2D/system/fvSolution
@@ -112,7 +112,9 @@ PIMPLE
 relaxationFactors
 {
     U               1;
+    UFinal          1;
     h               1;
+    hFinal          1;
     k               1;
     epsilon         1;
 }
diff --git a/tutorials/compressible/rhoPorousSimpleFoam/angledDuctExplicit/system/fvSchemes b/tutorials/compressible/rhoPorousSimpleFoam/angledDuctExplicit/system/fvSchemes
index ec261c3ee..2981c1b64 100644
--- a/tutorials/compressible/rhoPorousSimpleFoam/angledDuctExplicit/system/fvSchemes
+++ b/tutorials/compressible/rhoPorousSimpleFoam/angledDuctExplicit/system/fvSchemes
@@ -30,7 +30,7 @@ gradSchemes
 divSchemes
 {
     div(phi,U)      Gauss upwind;
-    div((muEff*dev2(grad(U).T()))) Gauss linear;
+    div((muEff*dev2(T(grad(U))))) Gauss linear;
     div(phi,h)      Gauss upwind;
     div(phi,epsilon) Gauss upwind;
     div(phi,k)      Gauss upwind;
@@ -59,10 +59,4 @@ snGradSchemes
     default         corrected;
 }
 
-fluxRequired
-{
-    default         no;
-    p               ;
-}
-
 // ************************************************************************* //
diff --git a/tutorials/compressible/rhoPorousSimpleFoam/angledDuctImplicit/system/fvSchemes b/tutorials/compressible/rhoPorousSimpleFoam/angledDuctImplicit/system/fvSchemes
index 985fae22f..97d7e575a 100644
--- a/tutorials/compressible/rhoPorousSimpleFoam/angledDuctImplicit/system/fvSchemes
+++ b/tutorials/compressible/rhoPorousSimpleFoam/angledDuctImplicit/system/fvSchemes
@@ -30,7 +30,7 @@ gradSchemes
 divSchemes
 {
     div(phi,U)      Gauss upwind;
-    div((muEff*dev2(grad(U).T()))) Gauss linear;
+    div((muEff*dev2(T(grad(U))))) Gauss linear;
     div(phi,h)      Gauss upwind;
     div(phi,epsilon) Gauss upwind;
     div(phi,k)      Gauss upwind;
@@ -58,10 +58,4 @@ snGradSchemes
     default         corrected;
 }
 
-fluxRequired
-{
-    default         no;
-    p               ;
-}
-
 // ************************************************************************* //
diff --git a/tutorials/compressible/rhoPorousSimpleFoam/beaverJoseph/0/T b/tutorials/compressible/rhoPorousSimpleFoam/beaverJoseph/0/T
index a68c56a13..3bd0d273c 100644
--- a/tutorials/compressible/rhoPorousSimpleFoam/beaverJoseph/0/T
+++ b/tutorials/compressible/rhoPorousSimpleFoam/beaverJoseph/0/T
@@ -27,16 +27,20 @@ boundaryField
     }
     outlet
     {
-        type inletOutlet;
-        inletValue $internalField;
+        type            inletOutlet;
+        inletValue      $internalField;
     }
     lowerWall
     {
-        type zeroGradient;
+        type            zeroGradient;
     }
     upperWall
     {
-        type zeroGradient;
+        type            zeroGradient;
+    }
+    frontAndBack
+    {
+        type            empty;
     }
 }
 
diff --git a/tutorials/compressible/rhoPorousSimpleFoam/beaverJoseph/0/U b/tutorials/compressible/rhoPorousSimpleFoam/beaverJoseph/0/U
index a02bcd0b3..4e7644bbe 100644
--- a/tutorials/compressible/rhoPorousSimpleFoam/beaverJoseph/0/U
+++ b/tutorials/compressible/rhoPorousSimpleFoam/beaverJoseph/0/U
@@ -45,6 +45,10 @@ boundaryField
         type            fixedValue;
         value           uniform (0 0 0);
     }
+    frontAndBack
+    {
+        type            empty;
+    }
 }
 
 // ************************************************************************* //
diff --git a/tutorials/compressible/rhoPorousSimpleFoam/beaverJoseph/0/p b/tutorials/compressible/rhoPorousSimpleFoam/beaverJoseph/0/p
index 9fb308d2b..2983e7faa 100644
--- a/tutorials/compressible/rhoPorousSimpleFoam/beaverJoseph/0/p
+++ b/tutorials/compressible/rhoPorousSimpleFoam/beaverJoseph/0/p
@@ -35,7 +35,11 @@ boundaryField
     outlet
     {
         type            fixedValue;
-        value $internalField;
+        value           $internalField;
+    }
+    frontAndBack
+    {
+        type            empty;
     }
 }
 
diff --git a/tutorials/compressible/rhoPorousSimpleFoam/beaverJoseph/system/fvSchemes b/tutorials/compressible/rhoPorousSimpleFoam/beaverJoseph/system/fvSchemes
index 985fae22f..97d7e575a 100644
--- a/tutorials/compressible/rhoPorousSimpleFoam/beaverJoseph/system/fvSchemes
+++ b/tutorials/compressible/rhoPorousSimpleFoam/beaverJoseph/system/fvSchemes
@@ -30,7 +30,7 @@ gradSchemes
 divSchemes
 {
     div(phi,U)      Gauss upwind;
-    div((muEff*dev2(grad(U).T()))) Gauss linear;
+    div((muEff*dev2(T(grad(U))))) Gauss linear;
     div(phi,h)      Gauss upwind;
     div(phi,epsilon) Gauss upwind;
     div(phi,k)      Gauss upwind;
@@ -58,10 +58,4 @@ snGradSchemes
     default         corrected;
 }
 
-fluxRequired
-{
-    default         no;
-    p               ;
-}
-
 // ************************************************************************* //
diff --git a/tutorials/compressible/rhoPorousSimpleFoam/porousPlug/0/T b/tutorials/compressible/rhoPorousSimpleFoam/porousPlug/0/T
index a68c56a13..4f2711da2 100644
--- a/tutorials/compressible/rhoPorousSimpleFoam/porousPlug/0/T
+++ b/tutorials/compressible/rhoPorousSimpleFoam/porousPlug/0/T
@@ -38,6 +38,10 @@ boundaryField
     {
         type zeroGradient;
     }
+    defaultFaces
+    {
+        type            empty;
+    }
 }
 
 // ************************************************************************* //
diff --git a/tutorials/compressible/rhoPorousSimpleFoam/porousPlug/0/U b/tutorials/compressible/rhoPorousSimpleFoam/porousPlug/0/U
index df82aea72..60530c711 100644
--- a/tutorials/compressible/rhoPorousSimpleFoam/porousPlug/0/U
+++ b/tutorials/compressible/rhoPorousSimpleFoam/porousPlug/0/U
@@ -43,6 +43,10 @@ boundaryField
         type            fixedValue;
         value           uniform (0 0 0);
     }
+    defaultFaces
+    {
+        type            empty;
+    }
 }
 
 // ************************************************************************* //
diff --git a/tutorials/compressible/rhoPorousSimpleFoam/porousPlug/0/p b/tutorials/compressible/rhoPorousSimpleFoam/porousPlug/0/p
index 9fb308d2b..e152e8bc0 100644
--- a/tutorials/compressible/rhoPorousSimpleFoam/porousPlug/0/p
+++ b/tutorials/compressible/rhoPorousSimpleFoam/porousPlug/0/p
@@ -37,6 +37,10 @@ boundaryField
         type            fixedValue;
         value $internalField;
     }
+    defaultFaces
+    {
+        type            empty;
+    }
 }
 
 // ************************************************************************* //
diff --git a/tutorials/compressible/rhoPorousSimpleFoam/porousPlug/system/fvSchemes b/tutorials/compressible/rhoPorousSimpleFoam/porousPlug/system/fvSchemes
index 985fae22f..97d7e575a 100644
--- a/tutorials/compressible/rhoPorousSimpleFoam/porousPlug/system/fvSchemes
+++ b/tutorials/compressible/rhoPorousSimpleFoam/porousPlug/system/fvSchemes
@@ -30,7 +30,7 @@ gradSchemes
 divSchemes
 {
     div(phi,U)      Gauss upwind;
-    div((muEff*dev2(grad(U).T()))) Gauss linear;
+    div((muEff*dev2(T(grad(U))))) Gauss linear;
     div(phi,h)      Gauss upwind;
     div(phi,epsilon) Gauss upwind;
     div(phi,k)      Gauss upwind;
@@ -58,10 +58,4 @@ snGradSchemes
     default         corrected;
 }
 
-fluxRequired
-{
-    default         no;
-    p               ;
-}
-
 // ************************************************************************* //
diff --git a/tutorials/compressible/rhoSonicFoam/forwardStep/system/fvSchemes b/tutorials/compressible/rhoSonicFoam/forwardStep/system/fvSchemes
index bbbf51c95..181d65a72 100644
--- a/tutorials/compressible/rhoSonicFoam/forwardStep/system/fvSchemes
+++ b/tutorials/compressible/rhoSonicFoam/forwardStep/system/fvSchemes
@@ -51,10 +51,4 @@ snGradSchemes
     default         corrected;
 }
 
-fluxRequired
-{
-    default         no;
-    p               ;
-}
-
 // ************************************************************************* //
diff --git a/tutorials/compressible/rhoSonicFoam/shockTube/0/T b/tutorials/compressible/rhoSonicFoam/shockTube/0/T
index 88a96c6a3..c73accee3 100644
--- a/tutorials/compressible/rhoSonicFoam/shockTube/0/T
+++ b/tutorials/compressible/rhoSonicFoam/shockTube/0/T
@@ -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
 {
@@ -17,7 +17,7 @@ FoamFile
 
 dimensions      [0 0 0 1 0 0 0];
 
-internalField   nonuniform List<scalar>
+internalField   nonuniform List<scalar> 
 1000
 (
 348.432
@@ -1035,4 +1035,5 @@ boundaryField
     }
 }
 
+
 // ************************************************************************* //
diff --git a/tutorials/compressible/rhoSonicFoam/shockTube/0/U b/tutorials/compressible/rhoSonicFoam/shockTube/0/U
index dfd6a3b28..0aef90a6e 100644
--- a/tutorials/compressible/rhoSonicFoam/shockTube/0/U
+++ b/tutorials/compressible/rhoSonicFoam/shockTube/0/U
@@ -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
 {
@@ -31,4 +31,5 @@ boundaryField
     }
 }
 
+
 // ************************************************************************* //
diff --git a/tutorials/compressible/rhoSonicFoam/shockTube/0/p b/tutorials/compressible/rhoSonicFoam/shockTube/0/p
index ccc492aa0..0b19ee1b3 100644
--- a/tutorials/compressible/rhoSonicFoam/shockTube/0/p
+++ b/tutorials/compressible/rhoSonicFoam/shockTube/0/p
@@ -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
 {
@@ -17,7 +17,7 @@ FoamFile
 
 dimensions      [1 -1 -2 0 0 0 0];
 
-internalField   nonuniform List<scalar>
+internalField   nonuniform List<scalar> 
 1000
 (
 100000
@@ -1035,4 +1035,5 @@ boundaryField
     }
 }
 
+
 // ************************************************************************* //
diff --git a/tutorials/compressible/rhoSonicFoam/shockTube/system/fvSchemes b/tutorials/compressible/rhoSonicFoam/shockTube/system/fvSchemes
index bbbf51c95..181d65a72 100644
--- a/tutorials/compressible/rhoSonicFoam/shockTube/system/fvSchemes
+++ b/tutorials/compressible/rhoSonicFoam/shockTube/system/fvSchemes
@@ -51,10 +51,4 @@ snGradSchemes
     default         corrected;
 }
 
-fluxRequired
-{
-    default         no;
-    p               ;
-}
-
 // ************************************************************************* //
diff --git a/tutorials/compressible/rhopSonicFoam/shockTube/0/T b/tutorials/compressible/rhopSonicFoam/shockTube/0/T
index bdd6097ef..f3385c002 100644
--- a/tutorials/compressible/rhopSonicFoam/shockTube/0/T
+++ b/tutorials/compressible/rhopSonicFoam/shockTube/0/T
@@ -1,22 +1,127 @@
 /*--------------------------------*- 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
 {
     version     2.0;
     format      ascii;
     class       volScalarField;
+    location    "0";
     object      T;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 dimensions      [0 0 0 1 0 0 0];
 
-internalField   uniform 1;
+internalField   nonuniform List<scalar> 
+100
+(
+348.432
+348.432
+348.432
+348.432
+348.432
+348.432
+348.432
+348.432
+348.432
+348.432
+348.432
+348.432
+348.432
+348.432
+348.432
+348.432
+348.432
+348.432
+348.432
+348.432
+348.432
+348.432
+348.432
+348.432
+348.432
+348.432
+348.432
+348.432
+348.432
+348.432
+348.432
+348.432
+348.432
+348.432
+348.432
+348.432
+348.432
+348.432
+348.432
+348.432
+348.432
+348.432
+348.432
+348.432
+348.432
+348.432
+348.432
+348.432
+348.432
+348.432
+278.746
+278.746
+278.746
+278.746
+278.746
+278.746
+278.746
+278.746
+278.746
+278.746
+278.746
+278.746
+278.746
+278.746
+278.746
+278.746
+278.746
+278.746
+278.746
+278.746
+278.746
+278.746
+278.746
+278.746
+278.746
+278.746
+278.746
+278.746
+278.746
+278.746
+278.746
+278.746
+278.746
+278.746
+278.746
+278.746
+278.746
+278.746
+278.746
+278.746
+278.746
+278.746
+278.746
+278.746
+278.746
+278.746
+278.746
+278.746
+278.746
+278.746
+)
+;
 
 boundaryField
 {
@@ -24,11 +129,11 @@ boundaryField
     {
         type            zeroGradient;
     }
-
     empty
     {
         type            empty;
     }
 }
 
+
 // ************************************************************************* //
diff --git a/tutorials/compressible/rhopSonicFoam/shockTube/0/U b/tutorials/compressible/rhopSonicFoam/shockTube/0/U
index 12c79ed59..0aef90a6e 100644
--- a/tutorials/compressible/rhopSonicFoam/shockTube/0/U
+++ b/tutorials/compressible/rhopSonicFoam/shockTube/0/U
@@ -1,15 +1,16 @@
 /*--------------------------------*- 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
 {
     version     2.0;
     format      ascii;
     class       volVectorField;
+    location    "0";
     object      U;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -24,11 +25,11 @@ boundaryField
     {
         type            zeroGradient;
     }
-
     empty
     {
         type            empty;
     }
 }
 
+
 // ************************************************************************* //
diff --git a/tutorials/compressible/rhopSonicFoam/shockTube/0/p b/tutorials/compressible/rhopSonicFoam/shockTube/0/p
index 96d0f410e..438517d73 100644
--- a/tutorials/compressible/rhopSonicFoam/shockTube/0/p
+++ b/tutorials/compressible/rhopSonicFoam/shockTube/0/p
@@ -1,22 +1,127 @@
 /*--------------------------------*- 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
 {
     version     2.0;
     format      ascii;
     class       volScalarField;
+    location    "0";
     object      p;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 dimensions      [1 -1 -2 0 0 0 0];
 
-internalField   uniform 0;
+internalField   nonuniform List<scalar> 
+100
+(
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+100000
+10000
+10000
+10000
+10000
+10000
+10000
+10000
+10000
+10000
+10000
+10000
+10000
+10000
+10000
+10000
+10000
+10000
+10000
+10000
+10000
+10000
+10000
+10000
+10000
+10000
+10000
+10000
+10000
+10000
+10000
+10000
+10000
+10000
+10000
+10000
+10000
+10000
+10000
+10000
+10000
+10000
+10000
+10000
+10000
+10000
+10000
+10000
+10000
+10000
+10000
+)
+;
 
 boundaryField
 {
@@ -24,11 +129,11 @@ boundaryField
     {
         type            zeroGradient;
     }
-
     empty
     {
         type            empty;
     }
 }
 
+
 // ************************************************************************* //
diff --git a/tutorials/compressible/rhopSonicFoam/shockTube/constant/polyMesh/boundary b/tutorials/compressible/rhopSonicFoam/shockTube/constant/polyMesh/boundary
index c81fa42aa..53a52f02b 100644
--- a/tutorials/compressible/rhopSonicFoam/shockTube/constant/polyMesh/boundary
+++ b/tutorials/compressible/rhopSonicFoam/shockTube/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/compressible/rhopSonicFoam/shockTube/system/fvSchemes b/tutorials/compressible/rhopSonicFoam/shockTube/system/fvSchemes
index c7dee39fb..902ca1719 100644
--- a/tutorials/compressible/rhopSonicFoam/shockTube/system/fvSchemes
+++ b/tutorials/compressible/rhopSonicFoam/shockTube/system/fvSchemes
@@ -46,10 +46,4 @@ snGradSchemes
     default         corrected;
 }
 
-fluxRequired
-{
-    default         no;
-    p               ;
-}
-
 // ************************************************************************* //
diff --git a/tutorials/compressible/rhopSonicFoam/shockTube/system/fvSolution b/tutorials/compressible/rhopSonicFoam/shockTube/system/fvSolution
index 3bb92efb1..6a581ece1 100644
--- a/tutorials/compressible/rhopSonicFoam/shockTube/system/fvSolution
+++ b/tutorials/compressible/rhopSonicFoam/shockTube/system/fvSolution
@@ -58,7 +58,7 @@ solvers
     }
 }
 
-PISO
+PIMPLE
 {
     nOuterCorrectors 3;
     nCorrectors     1;
diff --git a/tutorials/compressible/rhopSonicFoam/wedge15Ma5/constant/polyMesh/boundary b/tutorials/compressible/rhopSonicFoam/wedge15Ma5/constant/polyMesh/boundary
index b7d6204fe..bec7824fb 100644
--- a/tutorials/compressible/rhopSonicFoam/wedge15Ma5/constant/polyMesh/boundary
+++ b/tutorials/compressible/rhopSonicFoam/wedge15Ma5/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/compressible/rhopSonicFoam/wedge15Ma5/system/fvSchemes b/tutorials/compressible/rhopSonicFoam/wedge15Ma5/system/fvSchemes
index c92521cc5..902ca1719 100644
--- a/tutorials/compressible/rhopSonicFoam/wedge15Ma5/system/fvSchemes
+++ b/tutorials/compressible/rhopSonicFoam/wedge15Ma5/system/fvSchemes
@@ -46,10 +46,4 @@ snGradSchemes
     default         corrected;
 }
 
-fluxRequired
-{
-    default         yes;
-    p               ;
-}
-
 // ************************************************************************* //
diff --git a/tutorials/compressible/rhopSonicFoam/wedge15Ma5/system/fvSolution b/tutorials/compressible/rhopSonicFoam/wedge15Ma5/system/fvSolution
index f7cc03d3a..dc452d140 100644
--- a/tutorials/compressible/rhopSonicFoam/wedge15Ma5/system/fvSolution
+++ b/tutorials/compressible/rhopSonicFoam/wedge15Ma5/system/fvSolution
@@ -58,7 +58,7 @@ solvers
     }
 }
 
-PISO
+PIMPLE
 {
     nOuterCorrectors 2;
     nCorrectors     2;