diff --git a/applications/solvers/combustion/PDRFoam/PDRFoam.C b/applications/solvers/combustion/PDRFoam/PDRFoam.C
index e3a5e8c1f..3d4311526 100644
--- a/applications/solvers/combustion/PDRFoam/PDRFoam.C
+++ b/applications/solvers/combustion/PDRFoam/PDRFoam.C
@@ -64,6 +64,7 @@ Description
 #include "ignition.H"
 #include "Switch.H"
 #include "bound.H"
+#include "pisoControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -73,6 +74,9 @@ int main(int argc, char *argv[])
 
 #   include "createTime.H"
 #   include "createMesh.H"
+
+    pisoControl piso(mesh);
+
 #   include "readCombustionProperties.H"
 #   include "readGravitationalAcceleration.H"
 #   include "createFields.H"
@@ -90,7 +94,6 @@ int main(int argc, char *argv[])
     while (runTime.run())
     {
 #       include "readTimeControls.H"
-#       include "readPISOControls.H"
 #       include "CourantNo.H"
 #       include "setDeltaT.H"
 
@@ -102,7 +105,7 @@ int main(int argc, char *argv[])
 #       include "UEqn.H"
 
         // --- PISO loop
-        for (int corr=1; corr<=nCorr; corr++)
+        while (piso.correct())
         {
 #           include "bEqn.H"
 #           include "ftEqn.H"
diff --git a/applications/solvers/combustion/PDRFoam/PDRFoamAutoRefine.C b/applications/solvers/combustion/PDRFoam/PDRFoamAutoRefine.C
index fee3bd322..263e6b90b 100644
--- a/applications/solvers/combustion/PDRFoam/PDRFoamAutoRefine.C
+++ b/applications/solvers/combustion/PDRFoam/PDRFoamAutoRefine.C
@@ -66,6 +66,7 @@ Description
 #include "Switch.H"
 #include "bound.H"
 #include "dynamicRefineFvMesh.H"
+#include "pisoControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -75,10 +76,12 @@ int main(int argc, char *argv[])
 
 #   include "createTime.H"
 #   include "createDynamicFvMesh.H"
+
+    pisoControl piso(mesh);
+
 #   include "readCombustionProperties.H"
 #   include "readGravitationalAcceleration.H"
 #   include "createFields.H"
-#   include "readPISOControls.H"
 #   include "initContinuityErrs.H"
 #   include "createTimeControls.H"
 #   include "setInitialDeltaT.H"
@@ -92,7 +95,6 @@ int main(int argc, char *argv[])
     while (runTime.run())
     {
 #       include "readTimeControls.H"
-#       include "readPISOControls.H"
 #       include "CourantNo.H"
 
 #       include "setDeltaT.H"
@@ -168,7 +170,7 @@ int main(int argc, char *argv[])
 #       include "UEqn.H"
 
         // --- PISO loop
-        for (int corr=1; corr<=nCorr; corr++)
+        while (piso.correct())
         {
 #           include "bEqn.H"
 #           include "ftEqn.H"
diff --git a/applications/solvers/combustion/PDRFoam/UEqn.H b/applications/solvers/combustion/PDRFoam/UEqn.H
index aea6ad9a4..af2ea9510 100644
--- a/applications/solvers/combustion/PDRFoam/UEqn.H
+++ b/applications/solvers/combustion/PDRFoam/UEqn.H
@@ -9,7 +9,7 @@
 
     volSymmTensorField invA = inv(I*UEqn.A() + drag->Dcu());
 
-    if (momentumPredictor)
+    if (piso.momentumPredictor())
     {
         U = invA & (UEqn.H() - betav*fvc::grad(p));
         U.correctBoundaryConditions();
diff --git a/applications/solvers/combustion/PDRFoam/createFields.H b/applications/solvers/combustion/PDRFoam/createFields.H
index 74a18ab6f..6d3a4bdf5 100644
--- a/applications/solvers/combustion/PDRFoam/createFields.H
+++ b/applications/solvers/combustion/PDRFoam/createFields.H
@@ -163,3 +163,5 @@
     fields.add(h);
     fields.add(hu);
     flameWrinkling->addXi(fields);
+
+    mesh.schemesDict().setFluxRequired(p.name());
diff --git a/applications/solvers/combustion/PDRFoam/pEqn.H b/applications/solvers/combustion/PDRFoam/pEqn.H
index 524c8eac7..250fdf2a0 100644
--- a/applications/solvers/combustion/PDRFoam/pEqn.H
+++ b/applications/solvers/combustion/PDRFoam/pEqn.H
@@ -3,7 +3,7 @@ rho = thermo.rho();
 volScalarField rUA = 1.0/UEqn.A();
 U = invA & UEqn.H();
 
-if (transonic)
+if (piso.transonic())
 {
     surfaceScalarField phid
     (
@@ -15,7 +15,7 @@ if (transonic)
         )
     );
 
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    while (piso.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
@@ -26,7 +26,7 @@ if (transonic)
 
         pEqn.solve();
 
-        if (nonOrth == nNonOrthCorr)
+        if (piso.finalNonOrthogonalIter())
         {
             phi == pEqn.flux();
         }
@@ -41,7 +41,7 @@ else
           + fvc::ddtPhiCorr(rUA, rho, U, phi)
         );
 
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    while (piso.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
@@ -52,7 +52,7 @@ else
 
         pEqn.solve();
 
-        if (nonOrth == nNonOrthCorr)
+        if (piso.finalNonOrthogonalIter())
         {
             phi += pEqn.flux();
         }
diff --git a/applications/solvers/combustion/rhoReactingFoam/UEqn.H b/applications/solvers/combustion/rhoReactingFoam/UEqn.H
index 9697c6e1e..1e626d75b 100644
--- a/applications/solvers/combustion/rhoReactingFoam/UEqn.H
+++ b/applications/solvers/combustion/rhoReactingFoam/UEqn.H
@@ -9,7 +9,7 @@
 
     UEqn.relax();
 
-    if (momentumPredictor)
+    if (pimple.momentumPredictor())
     {
         solve(UEqn == -fvc::grad(p));
     }
diff --git a/applications/solvers/combustion/rhoReactingFoam/createFields.H b/applications/solvers/combustion/rhoReactingFoam/createFields.H
index d4119b8c9..035d5a77b 100644
--- a/applications/solvers/combustion/rhoReactingFoam/createFields.H
+++ b/applications/solvers/combustion/rhoReactingFoam/createFields.H
@@ -97,3 +97,5 @@ DimensionedField<scalar, volMesh> chemistrySh
     mesh,
     dimensionedScalar("chemistrySh", dimEnergy/dimTime/dimVolume, 0.0)
 );
+
+mesh.schemesDict().setFluxRequired(p.name());
diff --git a/applications/solvers/combustion/rhoReactingFoam/pEqn.H b/applications/solvers/combustion/rhoReactingFoam/pEqn.H
index 32653a927..222e8c959 100644
--- a/applications/solvers/combustion/rhoReactingFoam/pEqn.H
+++ b/applications/solvers/combustion/rhoReactingFoam/pEqn.H
@@ -8,7 +8,7 @@
     volScalarField rUA = 1.0/UEqn.A();
     U = rUA*UEqn.H();
 
-    if (transonic)
+    if (pimple.transonic())
     {
         surfaceScalarField phiv =
             (fvc::interpolate(U) & mesh.Sf())
@@ -22,7 +22,7 @@
             fvc::interpolate(thermo.psi())*phiv
         );
 
-        for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+        while (pimple.correctNonOrthogonal())
         {
             fvScalarMatrix pEqn
             (
@@ -31,21 +31,12 @@
               - fvm::laplacian(rho*rUA, p)
             );
 
-            if
+            pEqn.solve
             (
-                ocorr == nOuterCorr
-             && corr == nCorr
-             && nonOrth == nNonOrthCorr
-            )
-            {
-                pEqn.solve(mesh.solutionDict().solver(p.name() + "Final"));
-            }
-            else
-            {
-                pEqn.solve();
-            }
+                mesh.solutionDict().solver(p.select(pimple.finalInnerIter()))
+            );
 
-            if (nonOrth == nNonOrthCorr)
+            if (pimple.finalNonOrthogonalIter())
             {
                 phi += pEqn.flux();
             }
@@ -60,7 +51,7 @@
               + fvc::ddtPhiCorr(rUA, rho, U, phi)
             );
 
-        for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+        while (pimple.correctNonOrthogonal())
         {
             fvScalarMatrix pEqn
             (
@@ -69,21 +60,12 @@
               - fvm::laplacian(rho*rUA, p)
             );
 
-            if
+            pEqn.solve
             (
-                ocorr == nOuterCorr
-             && corr == nCorr
-             && nonOrth == nNonOrthCorr
-            )
-            {
-                pEqn.solve(mesh.solutionDict().solver(p.name() + "Final"));
-            }
-            else
-            {
-                pEqn.solve();
-            }
+                mesh.solutionDict().solver(p.select(pimple.finalInnerIter()))
+            );
 
-            if (nonOrth == nNonOrthCorr)
+            if (pimple.finalNonOrthogonalIter())
             {
                 phi += pEqn.flux();
             }
diff --git a/applications/solvers/combustion/rhoReactingFoam/rhoReactingFoam.C b/applications/solvers/combustion/rhoReactingFoam/rhoReactingFoam.C
index dc729c5f3..645d3a2cc 100644
--- a/applications/solvers/combustion/rhoReactingFoam/rhoReactingFoam.C
+++ b/applications/solvers/combustion/rhoReactingFoam/rhoReactingFoam.C
@@ -36,6 +36,7 @@ Description
 #include "rhoChemistryModel.H"
 #include "chemistrySolver.H"
 #include "multivariateScheme.H"
+#include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -44,6 +45,9 @@ int main(int argc, char *argv[])
 #   include "setRootCase.H"
 #   include "createTime.H"
 #   include "createMesh.H"
+
+    pimpleControl pimple(mesh);
+
 #   include "readChemistryProperties.H"
 #   include "readGravitationalAcceleration.H"
 #   include "createFields.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"
 
@@ -69,14 +72,14 @@ int main(int argc, char *argv[])
 #       include "chemistry.H"
 #       include "rhoEqn.H"
 
-        for (label ocorr=1; ocorr <= nOuterCorr; ocorr++)
+        while (pimple.loop())
         {
 #           include "UEqn.H"
 #           include "YEqn.H"
 #           include "hsEqn.H"
 
             // --- PISO loop
-            for (int corr=1; corr<=nCorr; corr++)
+            while (pimple.correct())
             {
 #               include "pEqn.H"
             }
diff --git a/applications/solvers/engine/sonicTurbDyMEngineFoam/UEqn.H b/applications/solvers/engine/sonicTurbDyMEngineFoam/UEqn.H
index 21a4d5fab..07f87de6a 100644
--- a/applications/solvers/engine/sonicTurbDyMEngineFoam/UEqn.H
+++ b/applications/solvers/engine/sonicTurbDyMEngineFoam/UEqn.H
@@ -6,13 +6,9 @@
       + turbulence->divDevRhoReff(U)
     );
 
-    if (oCorr == nOuterCorr - 1)
-    {
-        UEqn.relax(1);
-    }
-    else
-    {
-        UEqn.relax();
-    }
+    UEqn.relax
+    (
+        mesh.solutionDict().relaxationFactor(U.select(pimple.finalIter()))
+    );
 
     solve(UEqn == -fvc::grad(p));
diff --git a/applications/solvers/engine/sonicTurbDyMEngineFoam/createFields.H b/applications/solvers/engine/sonicTurbDyMEngineFoam/createFields.H
index 030e65640..6acbbcbb6 100644
--- a/applications/solvers/engine/sonicTurbDyMEngineFoam/createFields.H
+++ b/applications/solvers/engine/sonicTurbDyMEngineFoam/createFields.H
@@ -128,3 +128,5 @@
         ),
         fvc::interpolate(rho) * fvc::meshPhi(U)
     );
+
+    mesh.schemesDict().setFluxRequired(p.name());
diff --git a/applications/solvers/engine/sonicTurbDyMEngineFoam/pEqn.H b/applications/solvers/engine/sonicTurbDyMEngineFoam/pEqn.H
index 05075d8d7..45d692a63 100644
--- a/applications/solvers/engine/sonicTurbDyMEngineFoam/pEqn.H
+++ b/applications/solvers/engine/sonicTurbDyMEngineFoam/pEqn.H
@@ -4,18 +4,13 @@
     rUA = 1.0/UEqn.A();
     U = rUA*UEqn.H();
 
-    if (nOuterCorr != 1)
-    {
-        p.storePrevIter();
-    }
-
-    if (transonic)
+    if (pimple.transonic())
     {
         surfaceScalarField phid =
             fvc::interpolate(thermo.psi())*
             ((fvc::interpolate(U) & mesh.Sf()) - fvc::meshPhi(rho, U));
 
-        for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+        while (pimple.correctNonOrthogonal())
         {
             fvScalarMatrix pEqn
             (
@@ -26,7 +21,7 @@
 
             pEqn.solve();
 
-            if (nonOrth == nNonOrthCorr)
+            if (pimple.finalNonOrthogonalIter())
             {
                 phi == pEqn.flux();
             }
@@ -37,7 +32,7 @@
         phi = fvc::interpolate(rho)*
             ((fvc::interpolate(U) & mesh.Sf()) - fvc::meshPhi(rho, U));
 
-        for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+        while (pimple.correctNonOrthogonal())
         {
             fvScalarMatrix pEqn
             (
@@ -48,7 +43,7 @@
 
             pEqn.solve();
 
-            if (nonOrth == nNonOrthCorr)
+            if (pimple.finalNonOrthogonalIter())
             {
                 phi += pEqn.flux();
             }
@@ -56,7 +51,7 @@
     }
 
     // Explicitly relax pressure except for last corrector
-    if (oCorr != nOuterCorr - 1)
+    if (!pimple.finalIter())
     {
         p.relax();
     }
diff --git a/applications/solvers/engine/sonicTurbDyMEngineFoam/sonicTurbDyMEngineFoam.C b/applications/solvers/engine/sonicTurbDyMEngineFoam/sonicTurbDyMEngineFoam.C
index 04225e83b..5b2b6e2f3 100644
--- a/applications/solvers/engine/sonicTurbDyMEngineFoam/sonicTurbDyMEngineFoam.C
+++ b/applications/solvers/engine/sonicTurbDyMEngineFoam/sonicTurbDyMEngineFoam.C
@@ -39,6 +39,7 @@ Description
 #include "turbulenceModel.H"
 #include "Switch.H"
 #include "OFstream.H"
+#include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -49,7 +50,9 @@ int main(int argc, char *argv[])
 
 #   include "createEngineTime.H"
 #   include "createEngineDynamicMesh.H"
-#   include "readPIMPLEControls.H"
+
+    pimpleControl pimple(mesh);
+
 #   include "createFields.H"
 #   include "initContinuityErrs.H"
 #   include "createTimeControls.H"
@@ -67,7 +70,6 @@ int main(int argc, char *argv[])
 
     while (runTime.run())
     {
-#       include "readPIMPLEControls.H"
 #       include "checkTotalVolume.H"
 #       include "readEngineTimeControls.H"
 #       include "compressibleCourantNo.H"
@@ -106,19 +108,18 @@ int main(int argc, char *argv[])
         }
 
         // Pressure-velocity corrector
-        int oCorr = 0;
-        do
+        while (pimple.loop())
         {
 #           include "rhoEqn.H"
 #           include "UEqn.H"
 
             // --- PISO loop
-            for (int corr = 1; corr <= nCorr; corr++)
+            while (pimple.correct())
             {
 #               include "pEqn.H"
 #               include "hEqn.H"
             }
-        } while (++oCorr < nOuterCorr);
+        }
 
         turbulence->correct();
 
diff --git a/applications/solvers/engine/turbDyMEngineFoam/UEqn.H b/applications/solvers/engine/turbDyMEngineFoam/UEqn.H
index d60edb253..eb6b24fce 100644
--- a/applications/solvers/engine/turbDyMEngineFoam/UEqn.H
+++ b/applications/solvers/engine/turbDyMEngineFoam/UEqn.H
@@ -5,7 +5,7 @@
       + turbulence->divDevReff(U)
     );
 
-    if (momentumPredictor)
+    if (piso.momentumPredictor())
     {
         solve(UEqn == -fvc::grad(p));
     }
diff --git a/applications/solvers/engine/turbDyMEngineFoam/correctPhi.H b/applications/solvers/engine/turbDyMEngineFoam/correctPhi.H
index 5a8dc6a33..42cac00cf 100644
--- a/applications/solvers/engine/turbDyMEngineFoam/correctPhi.H
+++ b/applications/solvers/engine/turbDyMEngineFoam/correctPhi.H
@@ -28,7 +28,7 @@
 
     mesh.schemesDict().setFluxRequired(pcorr.name());
 
-    for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    while (piso.correctNonOrthogonal())
     {
         fvScalarMatrix pcorrEqn
         (
@@ -38,7 +38,7 @@
         pcorrEqn.setReference(pRefCell, pRefValue);
         pcorrEqn.solve();
 
-        if (nonOrth == nNonOrthCorr)
+        if (piso.finalNonOrthogonalIter())
         {
             phi -= pcorrEqn.flux();
         }
diff --git a/applications/solvers/engine/turbDyMEngineFoam/createControls.H b/applications/solvers/engine/turbDyMEngineFoam/createControls.H
index d76028fda..04d4121db 100644
--- a/applications/solvers/engine/turbDyMEngineFoam/createControls.H
+++ b/applications/solvers/engine/turbDyMEngineFoam/createControls.H
@@ -1,8 +1,6 @@
 #   include "createTimeControls.H"
-#   include "readPISOControls.H"
 
-    bool checkMeshCourantNo = false;
-    if (piso.found("checkMeshCourantNo"))
-    {
-        checkMeshCourantNo = Switch(piso.lookup("checkMeshCourantNo"));
-    }
+bool checkMeshCourantNo
+(
+    piso.dict().lookupOrDefault("checkMeshCourantNo", false)
+);
diff --git a/applications/solvers/engine/turbDyMEngineFoam/createFields.H b/applications/solvers/engine/turbDyMEngineFoam/createFields.H
index 2db24767d..9cad74217 100644
--- a/applications/solvers/engine/turbDyMEngineFoam/createFields.H
+++ b/applications/solvers/engine/turbDyMEngineFoam/createFields.H
@@ -33,7 +33,8 @@
 
     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());
 
     scalar totalVolume = sum(mesh.V()).value();
 
diff --git a/applications/solvers/engine/turbDyMEngineFoam/readControls.H b/applications/solvers/engine/turbDyMEngineFoam/readControls.H
index 375b524f2..bf221077f 100644
--- a/applications/solvers/engine/turbDyMEngineFoam/readControls.H
+++ b/applications/solvers/engine/turbDyMEngineFoam/readControls.H
@@ -1,8 +1,3 @@
 #   include "readTimeControls.H"
-#   include "readPISOControls.H"
 
-    checkMeshCourantNo = false;
-    if (piso.found("checkMeshCourantNo"))
-    {
-        checkMeshCourantNo = Switch(piso.lookup("checkMeshCourantNo"));
-    }
+checkMeshCourantNo = piso.dict().lookupOrDefault("checkMeshCourantNo", false);
diff --git a/applications/solvers/engine/turbDyMEngineFoam/turbDyMEngineFoam.C b/applications/solvers/engine/turbDyMEngineFoam/turbDyMEngineFoam.C
index 780a16a7f..1a8d46463 100644
--- a/applications/solvers/engine/turbDyMEngineFoam/turbDyMEngineFoam.C
+++ b/applications/solvers/engine/turbDyMEngineFoam/turbDyMEngineFoam.C
@@ -35,6 +35,7 @@ Description
 #include "turbulenceModel.H"
 #include "dynamicFvMesh.H"
 #include "engineTime.H"
+#include "pisoControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -45,6 +46,9 @@ int main(int argc, char *argv[])
 
 #   include "createEngineTime.H"
 #   include "createDynamicFvMesh.H"
+
+    pisoControl piso(mesh);
+
 #   include "initContinuityErrs.H"
 #   include "createFields.H"
 #   include "createControls.H"
@@ -89,7 +93,7 @@ int main(int argc, char *argv[])
 #       include "UEqn.H"
 
         // --- PISO loop
-        for (int corr = 0; corr < nCorr; corr++)
+        while (piso.correct())
         {
             rUA = 1.0/UEqn.A();
 
@@ -99,7 +103,7 @@ int main(int argc, char *argv[])
 
             adjustPhi(phi, U, p);
 
-            for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+            while (piso.correctNonOrthogonal())
             {
                 fvScalarMatrix pEqn
                 (
@@ -108,16 +112,12 @@ int main(int argc, char *argv[])
 
                 pEqn.setReference(pRefCell, pRefValue);
 
-                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();
                 }