Merged vanilla compatibility branch. Author: Henrik Rusche, FSB team. Merge: Henrik Rusche

Conflicts:
	applications/solvers/incompressible/channelFoam/channelFoam.C
	applications/solvers/incompressible/icoFoam/icoFoam.C
	applications/solvers/incompressible/nonNewtonianIcoFoam/nonNewtonianIcoFoam.C
	applications/solvers/incompressible/pimpleDyMFoam/UEqn.H
	applications/solvers/incompressible/pimpleDyMFoam/pimpleDyMFoam.C
	applications/solvers/incompressible/pimpleFoam/UEqn.H
	applications/solvers/incompressible/pimpleFoam/pEqn.H
	applications/solvers/incompressible/simpleFoam/UEqn.H
	applications/solvers/incompressible/simpleFoam/pEqn.H
	src/foam/fields/GeometricFields/GeometricTensorNFields/GeometricTensorNFields.H
	src/turbulenceModels/compressible/RAS/LRR/LRR.C
	src/turbulenceModels/compressible/RAS/LaunderSharmaKE/LaunderSharmaKE.C
	src/turbulenceModels/compressible/turbulenceModel/laminar/laminar.C
	src/turbulenceModels/incompressible/LES/GenEddyVisc/GenEddyVisc.C
	src/turbulenceModels/incompressible/LES/Smagorinsky2/Smagorinsky2.C
	src/turbulenceModels/incompressible/LES/SpalartAllmaras/SpalartAllmaras.C
	src/turbulenceModels/incompressible/LES/kOmegaSSTSAS/kOmegaSSTSAS.C
	src/turbulenceModels/incompressible/LES/laminar/laminar.C
	src/turbulenceModels/incompressible/RAS/LamBremhorstKE/LamBremhorstKE.C
	src/turbulenceModels/incompressible/RAS/LaunderSharmaKE/LaunderSharmaKE.C
	src/turbulenceModels/incompressible/RAS/LienCubicKE/LienCubicKE.C
	src/turbulenceModels/incompressible/RAS/NonlinearKEShih/NonlinearKEShih.C
	src/turbulenceModels/incompressible/RAS/RNGkEpsilon/RNGkEpsilon.C
	src/turbulenceModels/incompressible/RAS/SpalartAllmaras/SpalartAllmaras.C
	src/turbulenceModels/incompressible/RAS/coupledKEpsilon/coupledKEpsilon.C
	src/turbulenceModels/incompressible/RAS/kEpsilon/kEpsilon.C
	src/turbulenceModels/incompressible/RAS/kOmegaSST/kOmegaSST.C
	src/turbulenceModels/incompressible/RAS/laminar/laminar.C
	src/turbulenceModels/incompressible/RAS/qZeta/qZeta.C
	src/turbulenceModels/incompressible/RAS/realizableKE/realizableKE.C
	src/turbulenceModels/incompressible/turbulenceModel/laminar/laminar.C
	tutorials/basic/laplacianFoam/Case1.1_mixingPlane_sector_AB_60deg/system/fvSchemes
	tutorials/compressible/dbnsFoam/forwardStep/system/fvSchemes
	tutorials/compressible/dbnsFoam/shockTube/system/fvSchemes
	tutorials/compressible/dbnsTurbFoam/naca0012/system/fvSchemes
	tutorials/compressible/rhoPimpleFoam/angledDuct/system/fvSchemes
	tutorials/immersedBoundary/backwardStepShortTurbulentSimpleIbFoam/system/fvSchemes
	tutorials/immersedBoundary/pitzDailyBodyFitted/system/fvSchemes
	tutorials/immersedBoundary/pitzDailyLaminarSimpleIbFoam/system/fvSchemes
	tutorials/immersedBoundary/pitzDailyTurbulentSimpleIbFoam/system/fvSchemes
	tutorials/immersedBoundary/porousBumpChannelBodyFitted/system/fvSchemes
	tutorials/immersedBoundary/porousBumpChannelIb/system/fvSchemes
	tutorials/incompressible/MRFSimpleFoam/mixerVessel2D/system/fvSchemes
	tutorials/incompressible/boundaryFoam/boundaryLaunderSharma/system/fvSchemes
	tutorials/incompressible/boundaryFoam/boundaryTwoWallsFlowSolution/system/fvSchemes
	tutorials/incompressible/boundaryFoam/boundaryWallFunctions/system/fvSchemes
	tutorials/incompressible/channelFoam/channel395/system/fvSchemes
	tutorials/incompressible/icoFoam/elbow/constant/polyMesh/boundary
	tutorials/incompressible/pimpleDyMFoam/axialTurbine/system/fvSchemes
	tutorials/incompressible/pimpleDyMFoam/movingCone/system/fvSchemes
	tutorials/incompressible/pimpleDyMFoam/movingCone/system/fvSolution
	tutorials/incompressible/pimpleDyMFoam/movingCylinders/system/fvSchemes
	tutorials/incompressible/pimpleDyMFoam/wingMotion/wingMotion2D_pimpleDyMFoam/constant/polyMesh/boundary
	tutorials/incompressible/pimpleDyMFoam/wingMotion/wingMotion2D_pimpleDyMFoam/system/fvSchemes
	tutorials/incompressible/pimpleDyMFoam/wingMotion/wingMotion2D_simpleFoam/constant/polyMesh/boundary
	tutorials/incompressible/pimpleDyMFoam/wingMotion/wingMotion2D_simpleFoam/system/fvSchemes
	tutorials/incompressible/pimpleDyMFoam/wingMotion/wingMotion_snappyHexMesh/constant/polyMesh/boundary
	tutorials/incompressible/pimpleDyMFoam/wingMotion/wingMotion_snappyHexMesh/system/fvSchemes
	tutorials/incompressible/pimpleFoam/t-junction/system/fvSchemes
	tutorials/incompressible/pisoFoam/les/pitzDaily/system/fvSchemes
	tutorials/incompressible/pisoFoam/les/pitzDailyDirectMapped/system/fvSchemes
	tutorials/incompressible/pisoFoam/ras/cavity/system/fvSchemes
	tutorials/incompressible/porousSimpleFoam/angledDuctExplicit/system/fvSchemes
	tutorials/incompressible/porousSimpleFoam/angledDuctImplicit/system/fvSchemes
	tutorials/incompressible/porousSimpleFoam/beaverJoseph/system/fvSchemes
	tutorials/incompressible/porousSimpleFoam/porousPlug/system/fvSchemes
	tutorials/incompressible/simpleFoam/mixingPlaneAxial/system/fvSchemes
	tutorials/incompressible/simpleFoam/mixingPlaneDomADomB/system/fvSchemes
	tutorials/incompressible/simpleFoam/motorBike/system/fvSchemes
	tutorials/incompressible/simpleFoam/pitzDaily/system/fvSchemes
	tutorials/incompressible/simpleFoam/pitzDaily3Blocks/system/fvSchemes
	tutorials/incompressible/simpleFoam/pitzDailyExptInlet/system/fvSchemes
	tutorials/incompressible/simpleSRFFoam/mixer/system/fvSchemes

Fixing compilation of DBNS

applications/solvers/basic/laplacianFoam/laplacianFoam.C

@@ -30,7 +30,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-
+#include "simpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -41,6 +41,9 @@ int main(int argc, char *argv[])
 
 #   include "createTime.H"
 #   include "createMesh.H"
+
+    simpleControl simple(mesh);
+
 #   include "createFields.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -51,9 +54,7 @@ int main(int argc, char *argv[])
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-#       include "readSIMPLEControls.H"
+        while (simple.correctNonOrthogonal())
-
-        for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
         {
             solve
             (

applications/solvers/basic/potentialDyMFoam/correctPhi.H

@@ -34,6 +34,8 @@
         fvc::makeAbsolute(phi, U);
     }
 
+    mesh.schemesDict().setFluxRequired(pcorr.name());
+
     for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
         fvScalarMatrix pcorrEqn

applications/solvers/basic/potentialDyMFoam/createFields.H

@@ -48,3 +48,4 @@
     label pRefCell = 0;
     scalar pRefValue = 0.0;
     setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue);
+    mesh.schemesDict().setFluxRequired(p.name());

applications/solvers/basic/potentialDyMFoam/potentialDyMFoam.C

@@ -34,6 +34,7 @@ Author
 
 #include "fvCFD.H"
 #include "dynamicFvMesh.H"
+#include "pisoControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -45,6 +46,9 @@ int main(int argc, char *argv[])
 #   include "setRootCase.H"
 #   include "createTime.H"
 #   include "createDynamicFvMesh.H"
+
+    pisoControl piso(mesh);
+
 #   include "createFields.H"
 #   include "initTotalVolume.H"
 
@@ -54,7 +58,6 @@ int main(int argc, char *argv[])
 
     while (runTime.loop())
     {
-#       include "readPISOControls.H"
 #       include "checkTotalVolume.H"
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
@@ -75,7 +78,7 @@ int main(int argc, char *argv[])
         // Solve potential flow equations
         adjustPhi(phi, U, p);
 
-        for (int nonOrth = 0; nonOrth <= nNonOrthCorr; nonOrth++)
+        while (piso.correctNonOrthogonal())
         {
             p.storePrevIter();
 
@@ -98,7 +101,7 @@ int main(int argc, char *argv[])
             pEqn.setReference(pRefCell, pRefValue);
             pEqn.solve();
 
-            if (nonOrth == nNonOrthCorr)
+            if (piso.finalNonOrthogonalIter())
             {
                 phi -= pEqn.flux();
             }

applications/solvers/basic/potentialFoam/createFields.H

@@ -50,3 +50,4 @@
     label pRefCell = 0;
     scalar pRefValue = 0.0;
     setRefCell(p, mesh.solutionDict().subDict("SIMPLE"), pRefCell, pRefValue);
+    mesh.schemesDict().setFluxRequired(p.name());

applications/solvers/basic/potentialFoam/potentialFoam.C

@@ -31,6 +31,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
+#include "simpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -43,8 +44,10 @@ int main(int argc, char *argv[])
 
 #   include "createTime.H"
 #   include "createMesh.H"
+
+    simpleControl simple(mesh);
+
 #   include "createFields.H"
-#   include "readSIMPLEControls.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -52,7 +55,7 @@ int main(int argc, char *argv[])
 
     adjustPhi(phi, U, p);
 
-    for (int nonOrth = 0; nonOrth <= nNonOrthCorr; nonOrth++)
+    while (simple.correctNonOrthogonal())
     {
         p.storePrevIter();
 
@@ -75,7 +78,7 @@ int main(int argc, char *argv[])
         pEqn.setReference(pRefCell, pRefValue);
         pEqn.solve();
 
-        if (nonOrth == nNonOrthCorr)
+        if (simple.finalNonOrthogonalIter())
         {
             phi -= pEqn.flux();
         }

applications/solvers/basic/scalarTransportFoam/scalarTransportFoam.C

@@ -30,6 +30,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
+#include "simpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -39,6 +40,9 @@ int main(int argc, char *argv[])
 
 #   include "createTime.H"
 #   include "createMesh.H"
+
+    simpleControl simple(mesh);
+
 #   include "createFields.H"
 
 
@@ -52,9 +56,7 @@ int main(int argc, char *argv[])
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-#       include "readSIMPLEControls.H"
+        while (simple.correctNonOrthogonal())
-
-        for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
         {
             solve
             (

applications/solvers/combustion/PDRFoam/PDRFoam.C

@@ -77,7 +77,7 @@ int main(int argc, char *argv[])
 #   include "readGravitationalAcceleration.H"
 #   include "createFields.H"
 #   include "initContinuityErrs.H"
-#   include "readTimeControls.H"
+#   include "createTimeControls.H"
 #   include "CourantNo.H"
 #   include "setInitialDeltaT.H"
 

applications/solvers/combustion/PDRFoam/PDRFoamAutoRefine.C

@@ -80,7 +80,7 @@ int main(int argc, char *argv[])
 #   include "createFields.H"
 #   include "readPISOControls.H"
 #   include "initContinuityErrs.H"
-#   include "readTimeControls.H"
+#   include "createTimeControls.H"
 #   include "setInitialDeltaT.H"
 
     scalar StCoNum = 0.0;

applications/solvers/combustion/XiFoam/UEqn.H

@@ -9,7 +9,7 @@
 
     UEqn.relax();
 
-    if (momentumPredictor)
+    if (pimple.momentumPredictor())
     {
         solve(UEqn == -fvc::grad(p));
     }

applications/solvers/combustion/XiFoam/XiFoam.C

@@ -55,6 +55,7 @@ Description
 #include "laminarFlameSpeed.H"
 #include "ignition.H"
 #include "Switch.H"
+#include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -64,11 +65,14 @@ int main(int argc, char *argv[])
 
 #   include "createTime.H"
 #   include "createMesh.H"
+
+    pimpleControl pimple(mesh);
+
 #   include "readCombustionProperties.H"
 #   include "readGravitationalAcceleration.H"
 #   include "createFields.H"
 #   include "initContinuityErrs.H"
-#   include "readTimeControls.H"
+#   include "createTimeControls.H"
 #   include "compressibleCourantNo.H"
 #   include "setInitialDeltaT.H"
 
@@ -79,7 +83,6 @@ int main(int argc, char *argv[])
     while (runTime.run())
     {
 #       include "readTimeControls.H"
-#       include "readPISOControls.H"
 #       include "compressibleCourantNo.H"
 #       include "setDeltaT.H"
 
@@ -87,7 +90,7 @@ int main(int argc, char *argv[])
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
         // --- Pressure-velocity PIMPLE corrector loop
-        for (int oCorr=0; oCorr<nOuterCorr; oCorr++)
+        while (pimple.loop())
         {
 #           include "rhoEqn.H"
 #           include "UEqn.H"
@@ -103,7 +106,7 @@ int main(int argc, char *argv[])
             }
 
             // --- PISO loop
-            for (int corr=1; corr<=nCorr; corr++)
+            while (pimple.correct())
             {
 #              include "pEqn.H"
             }

applications/solvers/combustion/XiFoam/createFields.H

@@ -132,3 +132,5 @@
     fields.add(b);
     fields.add(h);
     fields.add(hu);
+
+    mesh.schemesDict().setFluxRequired(p.name());

applications/solvers/combustion/XiFoam/pEqn.H

@@ -3,7 +3,7 @@ rho = thermo.rho();
 volScalarField rUA = 1.0/UEqn.A();
 U = rUA*UEqn.H();
 
-if (transonic)
+if (pimple.transonic())
 {
     surfaceScalarField phid
     (
@@ -15,7 +15,7 @@ if (transonic)
         )
     );
 
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
@@ -26,7 +26,7 @@ if (transonic)
 
         pEqn.solve();
 
-        if (nonOrth == nNonOrthCorr)
+        if (pimple.finalNonOrthogonalIter())
         {
             phi == pEqn.flux();
         }
@@ -41,7 +41,7 @@ else
           + fvc::ddtPhiCorr(rUA, rho, U, phi)
         );
 
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
@@ -52,7 +52,7 @@ else
 
         pEqn.solve();
 
-        if (nonOrth == nNonOrthCorr)
+        if (pimple.finalNonOrthogonalIter())
         {
             phi += pEqn.flux();
         }

applications/solvers/combustion/coldEngineFoam/coldEngineFoam.C

@@ -35,6 +35,7 @@ Description
 #include "basicPsiThermo.H"
 #include "turbulenceModel.H"
 #include "OFstream.H"
+#include "pisoControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -44,8 +45,12 @@ int main(int argc, char *argv[])
 
 #   include "createEngineTime.H"
 #   include "createEngineMesh.H"
+
+    pisoControl piso(mesh);
+
 #   include "createFields.H"
 #   include "initContinuityErrs.H"
+#   include "createTimeControls.H"
 #   include "readEngineTimeControls.H"
 #   include "compressibleCourantNo.H"
 #   include "setInitialDeltaT.H"
@@ -57,7 +62,6 @@ int main(int argc, char *argv[])
 
     while (runTime.run())
     {
-#       include "readPISOControls.H"
 #       include "readEngineTimeControls.H"
 #       include "compressibleCourantNo.H"
 #       include "setDeltaT.H"
@@ -74,7 +78,7 @@ int main(int argc, char *argv[])
 #       include "UEqn.H"
 
         // --- PISO loop
-        for (int corr=1; corr<=nCorr; corr++)
+        while (piso.correct())
         {
 #           include "hEqn.H"
 #           include "pEqn.H"

applications/solvers/combustion/coldEngineFoam/createFields.H

@@ -57,3 +57,5 @@
     Info<< "Creating field DpDt\n" << endl;
     volScalarField DpDt =
         fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
+
+    mesh.schemesDict().setFluxRequired(p.name());

applications/solvers/combustion/dieselEngineFoam/UEqn.H

@@ -10,7 +10,7 @@
 
     UEqn.relax();
 
-    if (momentumPredictor)
+    if (pimple.momentumPredictor())
     {
         solve(UEqn == -fvc::grad(p));
     }

applications/solvers/combustion/dieselEngineFoam/createFields.H

@@ -107,3 +107,5 @@ DimensionedField<scalar, volMesh> chemistrySh
     mesh,
     dimensionedScalar("chemistrySh", dimEnergy/dimTime/dimVolume, 0.0)
 );
+
+mesh.schemesDict().setFluxRequired(p.name());

applications/solvers/combustion/dieselEngineFoam/dieselEngineFoam.C

@@ -42,6 +42,7 @@ Description
 #include "OFstream.H"
 #include "volPointInterpolation.H"
 #include "thermoPhysicsTypes.H"
+#include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -50,11 +51,15 @@ int main(int argc, char *argv[])
 #   include "setRootCase.H"
 #   include "createEngineTime.H"
 #   include "createEngineMesh.H"
+
+    pimpleControl pimple(mesh);
+
 #   include "createFields.H"
 #   include "readGravitationalAcceleration.H"
 #   include "readCombustionProperties.H"
 #   include "createSpray.H"
 #   include "initContinuityErrs.H"
+#   include "createTimeControls.H"
 #   include "readEngineTimeControls.H"
 #   include "compressibleCourantNo.H"
 #   include "setInitialDeltaT.H"
@@ -66,7 +71,6 @@ int main(int argc, char *argv[])
 
     while (runTime.run())
     {
-#       include "readPISOControls.H"
 #       include "readEngineTimeControls.H"
 #       include "compressibleCourantNo.H"
 #       include "setDeltaT.H"
@@ -102,13 +106,13 @@ int main(int argc, char *argv[])
 #       include "rhoEqn.H"
 #       include "UEqn.H"
 
-        for (label ocorr=1; ocorr <= nOuterCorr; ocorr++)
+        while (pimple.loop())
         {
 #           include "YEqn.H"
 #           include "hsEqn.H"
 
             // --- PISO loop
-            for (int corr=1; corr<=nCorr; corr++)
+            while (pimple.correct())
             {
 #               include "pEqn.H"
             }

applications/solvers/combustion/dieselEngineFoam/pEqn.H

@@ -3,7 +3,7 @@ rho = thermo.rho();
 volScalarField A = UEqn.A();
 U = UEqn.H()/A;
 
-if (transonic)
+if (pimple.transonic())
 {
     surfaceScalarField phid
     (
@@ -12,7 +12,7 @@ if (transonic)
        *((fvc::interpolate(U) & mesh.Sf()) - fvc::meshPhi(rho, U))
     );
 
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
@@ -25,7 +25,7 @@ if (transonic)
 
         pEqn.solve();
 
-        if (nonOrth == nNonOrthCorr)
+        if (pimple.finalNonOrthogonalIter())
         {
             phi == pEqn.flux();
         }
@@ -36,7 +36,7 @@ else
     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
         (
@@ -49,7 +49,7 @@ else
 
         pEqn.solve();
 
-        if (nonOrth == nNonOrthCorr)
+        if (pimple.finalNonOrthogonalIter())
         {
             phi += pEqn.flux();
         }

applications/solvers/combustion/dieselFoam/dieselFoam.C

@@ -40,6 +40,7 @@ Description
 #include "IFstream.H"
 #include "OFstream.H"
 #include "Switch.H"
+#include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -48,12 +49,15 @@ int main(int argc, char *argv[])
 #   include "setRootCase.H"
 #   include "createTime.H"
 #   include "createMesh.H"
+
+    pimpleControl pimple(mesh);
+
 #   include "createFields.H"
 #   include "readGravitationalAcceleration.H"
 #   include "readCombustionProperties.H"
 #   include "createSpray.H"
 #   include "initContinuityErrs.H"
-#   include "readTimeControls.H"
+#   include "createTimeControls.H"
 #   include "compressibleCourantNo.H"
 #   include "setInitialDeltaT.H"
 
@@ -63,7 +67,6 @@ int main(int argc, char *argv[])
 
     while (runTime.run())
     {
-#       include "readPISOControls.H"
 #       include "compressibleCourantNo.H"
 #       include "setDeltaT.H"
 
@@ -97,13 +100,13 @@ int main(int argc, char *argv[])
 #       include "rhoEqn.H"
 #       include "UEqn.H"
 
-        for (label ocorr=1; ocorr <= nOuterCorr; ocorr++)
+        while (pimple.loop())
         {
 #           include "YEqn.H"
 #           include "hsEqn.H"
 
             // --- PISO loop
-            for (int corr=1; corr<=nCorr; corr++)
+            while (pimple.correct())
             {
 #               include "pEqn.H"
             }

applications/solvers/combustion/dieselFoam/pEqn.H

@@ -3,7 +3,7 @@ rho = thermo.rho();
 volScalarField rUA = 1.0/UEqn.A();
 U = rUA*UEqn.H();
 
-if (transonic)
+if (pimple.transonic())
 {
     surfaceScalarField phid
     (
@@ -15,7 +15,7 @@ if (transonic)
         )
     );
 
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
@@ -28,7 +28,7 @@ if (transonic)
 
         pEqn.solve();
 
-        if (nonOrth == nNonOrthCorr)
+        if (pimple.finalNonOrthogonalIter())
         {
             phi == pEqn.flux();
         }
@@ -43,7 +43,7 @@ else
           + fvc::ddtPhiCorr(rUA, rho, U, phi)
         );
 
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
@@ -56,7 +56,7 @@ else
 
         pEqn.solve();
 
-        if (nonOrth == nNonOrthCorr)
+        if (pimple.finalNonOrthogonalIter())
         {
             phi += pEqn.flux();
         }

applications/solvers/combustion/engineFoam/UEqn.H

@@ -5,7 +5,7 @@
       + turbulence->divDevRhoReff(U)
     );
 
-    if (momentumPredictor)
+    if (piso.momentumPredictor())
     {
         solve(UEqn == -fvc::grad(p));
     }

applications/solvers/combustion/engineFoam/engineFoam.C

@@ -57,6 +57,7 @@ Description
 #include "ignition.H"
 #include "Switch.H"
 #include "OFstream.H"
+#include "pisoControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -66,9 +67,13 @@ int main(int argc, char *argv[])
 
 #   include "createEngineTime.H"
 #   include "createEngineMesh.H"
+
+    pisoControl piso(mesh);
+
 #   include "readCombustionProperties.H"
 #   include "createFields.H"
 #   include "initContinuityErrs.H"
+#   include "createTimeControls.H"
 #   include "readEngineTimeControls.H"
 #   include "compressibleCourantNo.H"
 #   include "setInitialDeltaT.H"
@@ -80,7 +85,6 @@ int main(int argc, char *argv[])
 
     while (runTime.run())
     {
-#       include "readPISOControls.H"
 #       include "readEngineTimeControls.H"
 #       include "compressibleCourantNo.H"
 #       include "setDeltaT.H"
@@ -96,7 +100,7 @@ int main(int argc, char *argv[])
 #       include "UEqn.H"
 
         // --- PISO loop
-        for (int corr=1; corr<=nCorr; corr++)
+        while (piso.correct()) 
         {
 #           include "ftEqn.H"
 #           include "bEqn.H"

applications/solvers/combustion/engineFoam/pEqn.H

@@ -3,7 +3,7 @@ rho = thermo.rho();
 volScalarField rUA = 1.0/UEqn.A();
 U = rUA*UEqn.H();
 
-if (transonic)
+if (piso.transonic())
 {
     surfaceScalarField phid
     (
@@ -12,7 +12,7 @@ if (transonic)
        *((fvc::interpolate(U) & mesh.Sf()) - fvc::meshPhi(rho, U))
     );
 
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    while (piso.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
@@ -23,7 +23,7 @@ if (transonic)
 
         pEqn.solve();
 
-        if (nonOrth == nNonOrthCorr)
+        if (piso.finalNonOrthogonalIter())
         {
             phi == pEqn.flux();
         }
@@ -34,7 +34,7 @@ else
     phi = fvc::interpolate(rho)
          *((fvc::interpolate(U) & mesh.Sf()) - fvc::meshPhi(rho, U));
 
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    while (piso.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
@@ -45,7 +45,7 @@ else
 
         pEqn.solve();
 
-        if (nonOrth == nNonOrthCorr)
+        if (piso.finalNonOrthogonalIter())
         {
             phi += pEqn.flux();
         }

applications/solvers/combustion/fireFoam/UEqn.H

@@ -7,8 +7,6 @@ fvVectorMatrix UEqn
 
 UEqn.relax();
 
-if (oCorr == nOuterCorr - 1)
-{
 solve
 (
     UEqn
@@ -18,19 +16,5 @@ if (oCorr == nOuterCorr - 1)
         fvc::interpolate(rho)*(g & mesh.Sf())
       - fvc::snGrad(p)*mesh.magSf()
     ),
-        mesh.solutionDict().solver("UFinal")
+    mesh.solutionDict().solver((U.select(pimple.finalIter())))
 );
-}
-else
-{
-    solve
-    (
-        UEqn
-      ==
-        fvc::reconstruct
-        (
-            fvc::interpolate(rho)*(g & mesh.Sf())
-          - fvc::snGrad(p)*mesh.magSf()
-        )
-    );
-}

applications/solvers/combustion/fireFoam/createFields.H

@@ -128,3 +128,5 @@ if (composition.contains("fu"))
 }
 
 fields.add(hs);
+
+mesh.schemesDict().setFluxRequired(p.name());

applications/solvers/combustion/fireFoam/fireFoam.C

@@ -33,6 +33,7 @@ Description
 #include "hsCombustionThermo.H"
 #include "turbulenceModel.H"
 #include "combustionModel.H"
+#include "pimpleControl.H"
 
 #include "radiationModel.H"
 
@@ -43,11 +44,14 @@ int main(int argc, char *argv[])
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createMesh.H"
+
+    pimpleControl pimple(mesh);
+
     #include "readGravitationalAcceleration.H"
     #include "initContinuityErrs.H"
     #include "createFields.H"
     #include "createRadiationModel.H"
-    #include "readTimeControls.H"
+    #include "createTimeControls.H"
     #include "compressibleCourantNo.H"
     #include "setInitialDeltaT.H"
 
@@ -57,7 +61,6 @@ int main(int argc, char *argv[])
 
     while (runTime.run())
     {
-        #include "readPISOControls.H"
         #include "readTimeControls.H"
         #include "compressibleCourantNo.H"
         #include "setDeltaT.H"
@@ -68,7 +71,7 @@ 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"
 
@@ -76,7 +79,7 @@ int main(int argc, char *argv[])
             #include "fuhsEqn.H"
 
             // --- PISO loop
-            for (int corr = 0; corr < nCorr; corr++)
+            while (pimple.correct())
             {
                 #include "pEqn.H"
             }

applications/solvers/combustion/fireFoam/pEqn.H

@@ -17,7 +17,7 @@ surfaceScalarField phiU
 
 phi = phiU + rhorUAf*fvc::interpolate(rho)*(g & mesh.Sf());
 
-for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+while (pimple.correctNonOrthogonal())
 {
        surfaceScalarField rhorUAf = fvc::interpolate(rho*rUA);
 
@@ -30,16 +30,12 @@ for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
 
         closedVolume = p.needReference();
 
-        if (corr == nCorr-1 && nonOrth == nNonOrthCorr)
+        pEqn.solve
-        {
+        (
-            pEqn.solve(mesh.solutionDict().solver(p.name() + "Final"));
+            mesh.solutionDict().solver(p.select(pimple.finalInnerIter()))
-        }
+        );
-        else
-        {
-            pEqn.solve(mesh.solutionDict().solver(p.name()));
-        }
 
-        if (nonOrth == nNonOrthCorr)
+        if (pimple.finalNonOrthogonalIter())
         {
             phi += pEqn.flux();
         }

applications/solvers/combustion/reactingFoam/UEqn.H

@@ -9,7 +9,7 @@
 
     UEqn.relax();
 
-    if (momentumPredictor)
+    if (pimple.momentumPredictor())
     {
         solve(UEqn == -fvc::grad(p));
     }

applications/solvers/combustion/reactingFoam/createFields.H

@@ -96,3 +96,5 @@ DimensionedField<scalar, volMesh> chemistrySh
     mesh,
     dimensionedScalar("chemistrySh", dimEnergy/dimTime/dimVolume, 0.0)
 );
+
+mesh.schemesDict().setFluxRequired(p.name());

applications/solvers/combustion/reactingFoam/pEqn.H

@@ -3,7 +3,7 @@ rho = thermo.rho();
 volScalarField rUA = 1.0/UEqn.A();
 U = rUA*UEqn.H();
 
-if (transonic)
+if (pimple.transonic())
 {
     surfaceScalarField phid
     (
@@ -15,7 +15,7 @@ if (transonic)
         )
     );
 
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
@@ -26,7 +26,7 @@ if (transonic)
 
         pEqn.solve();
 
-        if (nonOrth == nNonOrthCorr)
+        if (pimple.finalNonOrthogonalIter())
         {
             phi == pEqn.flux();
         }
@@ -41,7 +41,7 @@ else
           + fvc::ddtPhiCorr(rUA, rho, U, phi)
         );
 
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
@@ -52,7 +52,7 @@ else
 
         pEqn.solve();
 
-        if (nonOrth == nNonOrthCorr)
+        if (pimple.finalNonOrthogonalIter())
         {
             phi += pEqn.flux();
         }

applications/solvers/combustion/reactingFoam/reactingFoam.C

@@ -35,6 +35,7 @@ Description
 #include "psiChemistryModel.H"
 #include "chemistrySolver.H"
 #include "multivariateScheme.H"
+#include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -43,11 +44,14 @@ 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"
 #   include "initContinuityErrs.H"
-#   include "readTimeControls.H"
+#   include "createTimeControls.H"
 #   include "compressibleCourantNo.H"
 #   include "setInitialDeltaT.H"
 
@@ -58,7 +62,6 @@ int main(int argc, char *argv[])
     while (runTime.run())
     {
 #       include "readTimeControls.H"
-#       include "readPISOControls.H"
 #       include "compressibleCourantNo.H"
 #       include "setDeltaT.H"
 
@@ -68,14 +71,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"
             }

applications/solvers/combustion/rhoReactingFoam/rhoReactingFoam.C

@@ -48,7 +48,7 @@ int main(int argc, char *argv[])
 #   include "readGravitationalAcceleration.H"
 #   include "createFields.H"
 #   include "initContinuityErrs.H"
-#   include "readTimeControls.H"
+#   include "createTimeControls.H"
 #   include "compressibleCourantNo.H"
 #   include "setInitialDeltaT.H"
 

applications/solvers/compressible/dbnsFoam/dbnsFoam.C

@@ -54,6 +54,7 @@ int main(int argc, char *argv[])
 #   include "createTime.H"
 #   include "createMesh.H"
 #   include "createFields.H"
+#   include "createTimeControls.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

applications/solvers/compressible/dbnsTurbFoam/dbnsTurbFoam.C

@@ -55,6 +55,7 @@ int main(int argc, char *argv[])
 #   include "createTime.H"
 #   include "createMesh.H"
 #   include "createFields.H"
+#   include "createTimeControls.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

applications/solvers/compressible/rhoCentralFoam/rhoCentralFoam.C

@@ -45,7 +45,7 @@ int main(int argc, char *argv[])
 #   include "createMesh.H"
 #   include "createFields.H"
 #   include "readThermophysicalProperties.H"
-#   include "readTimeControls.H"
+#   include "createTimeControls.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

applications/solvers/compressible/rhoPimpleFoam/UEqn.H

@@ -7,34 +7,20 @@ tmp<fvVectorMatrix> UEqn
   + turbulence->divDevRhoReff(U)
 );
 
-if (oCorr == nOuterCorr - 1)
+UEqn().relax
-{
+(
-    if (mesh.solutionDict().relax("UFinal"))
+    mesh.solutionDict().relaxationFactor(U.select(pimple.finalIter()))
-    {
+);
-        UEqn().relax(mesh.solutionDict().relaxationFactor("UFinal"));
-    }
-    else
-    {
-        UEqn().relax(1);
-    }
-}
-else
-{
-    UEqn().relax();
-}
 
 volScalarField rUA = 1.0/UEqn().A();
 
-if (momentumPredictor)
+if (pimple.momentumPredictor())
 {
-    if (oCorr == nOuterCorr-1)
+    solve
-    {
+    (
-        solve(UEqn() == -fvc::grad(p), mesh.solutionDict().solver("UFinal"));
+        UEqn() == -fvc::grad(p),
-    }
+        mesh.solutionDict().solver((U.select(pimple.finalIter())))
-    else
+    );
-    {
-        solve(UEqn() == -fvc::grad(p));
-    }
 }
 else
 {

applications/solvers/compressible/rhoPimpleFoam/createFields.H

@@ -62,3 +62,5 @@
     Info<< "Creating field DpDt\n" << endl;
     volScalarField DpDt =
         fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
+
+    mesh.schemesDict().setFluxRequired(p.name());

applications/solvers/compressible/rhoPimpleFoam/hEqn.H

@@ -8,16 +8,11 @@
         DpDt
     );
 
-    if (oCorr == nOuterCorr-1)
+    hEqn.relax
-    {
+    (
-        hEqn.relax();
+        mesh.solutionDict().relaxationFactor(h.select(pimple.finalIter()))
-        hEqn.solve(mesh.solutionDict().solver("hFinal"));
+    );
-    }
+    hEqn.solve(mesh.solutionDict().solver((h.select(pimple.finalIter()))));
-    else
-    {
-        hEqn.relax();
-        hEqn.solve();
-    }
 
     thermo.correct();
 }

applications/solvers/compressible/rhoPimpleFoam/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
     (
@@ -20,7 +20,7 @@ if (transonic)
         )
     );
 
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
@@ -29,21 +29,12 @@ if (transonic)
           - fvm::laplacian(rho*rUA, p)
         );
 
-        if
+        pEqn.solve
         (
-            oCorr == nOuterCorr-1
+            mesh.solutionDict().solver(p.select(pimple.finalInnerIter()))
-            && corr == nCorr-1
+        );
-            && nonOrth == nNonOrthCorr
-        )
-        {
-            pEqn.solve(mesh.solutionDict().solver("pFinal"));
-        }
-        else
-        {
-            pEqn.solve();
-        }
 
-        if (nonOrth == nNonOrthCorr)
+        if (pimple.finalNonOrthogonalIter())
         {
             phi == pEqn.flux();
         }
@@ -58,7 +49,7 @@ else
         //+ fvc::ddtPhiCorr(rUA, rho, U, phi)
         );
 
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    while (pimple.correctNonOrthogonal())
     {
         // Pressure corrector
         fvScalarMatrix pEqn
@@ -68,21 +59,12 @@ else
           - fvm::laplacian(rho*rUA, p)
         );
 
-        if
+        pEqn.solve
         (
-            oCorr == nOuterCorr-1
+            mesh.solutionDict().solver(p.select(pimple.finalInnerIter()))
-         && corr == nCorr-1
+        );
-         && nonOrth == nNonOrthCorr
-        )
-        {
-            pEqn.solve(mesh.solutionDict().solver("pFinal"));
-        }
-        else
-        {
-            pEqn.solve();
-        }
 
-        if (nonOrth == nNonOrthCorr)
+        if (pimple.finalNonOrthogonalIter())
         {
             phi += pEqn.flux();
         }

applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C

@@ -37,6 +37,7 @@ Description
 #include "basicPsiThermo.H"
 #include "turbulenceModel.H"
 #include "bound.H"
+#include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -45,8 +46,12 @@ 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"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -55,7 +60,6 @@ int main(int argc, char *argv[])
     while (runTime.run())
     {
         #include "readTimeControls.H"
-        #include "readPIMPLEControls.H"
         #include "compressibleCourantNo.H"
         #include "setDeltaT.H"
 
@@ -63,7 +67,7 @@ int main(int argc, char *argv[])
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        if (nOuterCorr != 1)
+        if (!pimple.firstIter())
         {
             p.storePrevIter();
             rho.storePrevIter();
@@ -72,13 +76,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"
             }

applications/solvers/compressible/rhoPisoFoam/UEqn.H

@@ -5,7 +5,7 @@
       + turbulence->divDevRhoReff(U)
     );
 
-    if (momentumPredictor)
+    if (piso.momentumPredictor())
     {
         solve(UEqn == -fvc::grad(p));
     }

applications/solvers/compressible/rhoPisoFoam/createFields.H

@@ -56,3 +56,5 @@
     Info<< "Creating field DpDt\n" << endl;
     volScalarField DpDt =
         fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
+
+    mesh.schemesDict().setFluxRequired(p.name());

applications/solvers/compressible/rhoPisoFoam/pEqn.H

@@ -3,7 +3,7 @@ rho = thermo.rho();
 volScalarField rUA = 1.0/UEqn.A();
 U = rUA*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();
         }

applications/solvers/compressible/rhoPisoFoam/rhoPisoFoam.C

@@ -32,6 +32,7 @@ Description
 #include "fvCFD.H"
 #include "basicPsiThermo.H"
 #include "turbulenceModel.H"
+#include "pisoControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -41,9 +42,12 @@ int main(int argc, char *argv[])
 
     #include "createTime.H"
     #include "createMesh.H"
+
+    pisoControl piso(mesh);
+
     #include "createFields.H"
     #include "initContinuityErrs.H"
-    #include "readTimeControls.H"
+    #include "createTimeControls.H"
     #include "compressibleCourantNo.H"
     #include "setInitialDeltaT.H"
 
@@ -55,7 +59,6 @@ int main(int argc, char *argv[])
     while (runTime.run())
     {
         #include "readTimeControls.H"
-        #include "readPISOControls.H"
         #include "compressibleCourantNo.H"
         #include "setDeltaT.H"
 
@@ -67,7 +70,7 @@ int main(int argc, char *argv[])
         #include "UEqn.H"
 
         // --- PISO loop
-        for (int corr = 1; corr <= nCorr; corr++)
+        while (piso.correct())
         {
             #include "hEqn.H"
             #include "pEqn.H"

applications/solvers/compressible/rhoPorousMRFPimpleFoam/UEqn.H

@@ -7,37 +7,23 @@ tmp<fvVectorMatrix> UEqn
   + turbulence->divDevRhoReff(U)
 );
 
-if (oCorr == nOuterCorr-1)
+UEqn().relax
-{
+(
-    if (mesh.solutionDict().relax("UFinal"))
+    mesh.solutionDict().relaxationFactor(U.select(pimple.finalIter()))
-    {
+);
-        UEqn().relax(mesh.solutionDict().relaxationFactor("UFinal"));
-    }
-    else
-    {
-        UEqn().relax(1);
-    }
-}
-else
-{
-    UEqn().relax();
-}
 
 mrfZones.addCoriolis(rho, UEqn());
 pZones.addResistance(UEqn());
 
 volScalarField rUA = 1.0/UEqn().A();
 
-if (momentumPredictor)
+if (pimple.momentumPredictor())
 {
-    if (oCorr == nOuterCorr-1)
+    solve
-    {
+    (
-        solve(UEqn() == -fvc::grad(p), mesh.solutionDict().solver("UFinal"));
+        UEqn() == -fvc::grad(p),
-    }
+        mesh.solutionDict().solver((U.select(pimple.finalIter())))
-    else
+    );
-    {
-        solve(UEqn() == -fvc::grad(p));
-    }
 }
 else
 {

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
+            mesh.solutionDict().solver(p.select(pimple.finalInnerIter()))
-            && corr == nCorr-1
+        );
-            && nonOrth == nNonOrthCorr
-        )
-        {
-            pEqn.solve(mesh.solutionDict().solver("pFinal"));
-        }
-        else
-        {
-            pEqn.solve();
-        }
 
-        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
+            mesh.solutionDict().solver(p.select(pimple.finalInnerIter()))
-         && corr == nCorr-1
+        );
-         && nonOrth == nNonOrthCorr
-        )
-        {
-            pEqn.solve(mesh.solutionDict().solver("pFinal"));
-        }
-        else
-        {
-            pEqn.solve();
-        }
 
-        if (nonOrth == nNonOrthCorr)
+        if (pimple.finalNonOrthogonalIter())
         {
             phi += pEqn.flux();
         }

applications/solvers/compressible/rhoPorousMRFPimpleFoam/rhoPorousMRFPimpleFoam.C

@@ -39,6 +39,7 @@ Description
 #include "bound.H"
 #include "MRFZones.H"
 #include "porousZones.H"
+#include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -47,8 +48,12 @@ 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"
 
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -57,7 +62,6 @@ int main(int argc, char *argv[])
     while (runTime.run())
     {
         #include "readTimeControls.H"
-        #include "readPIMPLEControls.H"
         #include "compressibleCourantNo.H"
         #include "setDeltaT.H"
 
@@ -65,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();
@@ -74,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"
             }

applications/solvers/compressible/rhopSonicFoam/rhopSonicFoam.C

@@ -47,6 +47,7 @@ int main(int argc, char *argv[])
 #   include "createMesh.H"
 #   include "readThermodynamicProperties.H"
 #   include "createFields.H"
+#   include "createTimeControls.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

applications/solvers/compressible/sonicDyMFoam/createControls.H

@@ -1,4 +1,4 @@
-#   include "readTimeControls.H"
+#   include "createTimeControls.H"
 #   include "readPIMPLEControls.H"
 
     bool correctPhi = false;

applications/solvers/compressible/sonicDyMFoam/readControls.H

@@ -1,13 +1,13 @@
 #   include "readTimeControls.H"
 #   include "readPIMPLEControls.H"
 
-    bool correctPhi = false;
+    correctPhi = false;
     if (pimple.found("correctPhi"))
     {
         correctPhi = Switch(pimple.lookup("correctPhi"));
     }
 
-    bool checkMeshCourantNo = false;
+    checkMeshCourantNo = false;
     if (pimple.found("checkMeshCourantNo"))
     {
         checkMeshCourantNo = Switch(pimple.lookup("checkMeshCourantNo"));

applications/solvers/compressible/sonicDyMFoam/sonicDyMFoam.C

@@ -58,6 +58,7 @@ int main(int argc, char *argv[])
 #   include "createDynamicFvMesh.H"
 #   include "createFields.H"
 #   include "initContinuityErrs.H"
+#   include "createControls.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

applications/solvers/compressible/sonicFoam/sonicFoam.C

@@ -52,6 +52,7 @@ int main(int argc, char *argv[])
 #   include "createMesh.H"
 #   include "createFields.H"
 #   include "initContinuityErrs.H"
+#   include "createTimeControls.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

applications/solvers/coupled/conjugateHeatFoam/conjugateHeatFoam.C

@@ -51,7 +51,7 @@ int main(int argc, char *argv[])
 #   include "createFields.H"
 #   include "createSolidFields.H"
 #   include "initContinuityErrs.H"
-#   include "readTimeControls.H"
+#   include "createTimeControls.H"
 #   include "CourantNo.H"
 #   include "setInitialDeltaT.H"
 

applications/solvers/engine/icoDyMEngineFoam/correctPhi.H

@@ -24,9 +24,11 @@
         pcorrTypes
     );
 
-    phi = (fvc::interpolate(U) & mesh.Sf());
+    phi = fvc::interpolate(U) & mesh.Sf();
 
-    for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    mesh.schemesDict().setFluxRequired(pcorr.name());
+
+    while (piso.correctNonOrthogonal())
     {
         fvScalarMatrix pcorrEqn
         (
@@ -36,7 +38,7 @@
         pcorrEqn.setReference(pRefCell, pRefValue);
         pcorrEqn.solve();
 
-        if (nonOrth == nNonOrthCorr)
+        if (piso.finalNonOrthogonalIter())
         {
             phi -= pcorrEqn.flux();
         }

applications/solvers/engine/icoDyMEngineFoam/createControls.H

@@ -0,0 +1,6 @@
+#   include "createTimeControls.H"
+
+    bool checkMeshCourantNo
+    (
+        piso.dict().lookupOrDefault("checkMeshCourantNo", false)
+    );

applications/solvers/engine/icoDyMEngineFoam/createFields.H

@@ -53,6 +53,7 @@
     label pRefCell = 0;
     scalar pRefValue = 0.0;
     setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue);
+    mesh.schemesDict().setFluxRequired(p.name());
 
     scalar totalVolume = sum(mesh.V()).value();
 
@@ -69,16 +70,9 @@
         zeroGradientFvPatchScalarField::typeName
     );
 
+    bool momentumPredictor = piso.momentumPredictor();
 
     {
-        const dictionary& piso = mesh.solutionDict().subDict("PISO");
-
-        bool momentumPredictor = true;
-        if (piso.found("momentumPredictor"))
-        {
-            momentumPredictor = Switch(piso.lookup("momentumPredictor"));
-        }
-
         bool momentumPredictorSave = momentumPredictor;
         momentumPredictor = false;
 #       include "UEqn.H"

applications/solvers/engine/icoDyMEngineFoam/icoDyMEngineFoam.C

@@ -36,6 +36,7 @@ Author
 #include "fvCFD.H"
 #include "dynamicFvMesh.H"
 #include "engineTime.H"
+#include "pisoControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -45,8 +46,12 @@ int main(int argc, char *argv[])
 #   include "setRootCase.H"
 #   include "createEngineTime.H"
 #   include "createDynamicFvMesh.H"
+
+    pisoControl piso(mesh);
+
 #   include "initContinuityErrs.H"
 #   include "createFields.H"
+#   include "createControls.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -89,7 +94,7 @@ int main(int argc, char *argv[])
 
         // --- PISO loop
 
-        for (int corr = 0; corr < nCorr; corr++)
+        while (piso.correct())
         {
             rUA = 1.0/UEqn.A();
 
@@ -99,7 +104,8 @@ int main(int argc, char *argv[])
 
             adjustPhi(phi, U, p);
 
-            for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+            // Non-orthogonal pressure corrector loop
+            while (piso.correctNonOrthogonal())
             {
                 fvScalarMatrix pEqn
                 (
@@ -107,17 +113,12 @@ int main(int argc, char *argv[])
                 );
 
                 pEqn.setReference(pRefCell, pRefValue);
+                pEqn.solve
+                (
+                    mesh.solutionDict().solver(p.select(piso.finalInnerIter()))
+                );
 
-                if (corr == nCorr - 1 && nonOrth == nNonOrthCorr)
+                if (piso.finalNonOrthogonalIter())
-                {
-                    pEqn.solve(mesh.solutionDict().solver(p.name() + "Final"));
-                }
-                else
-                {
-                    pEqn.solve(mesh.solutionDict().solver(p.name()));
-                }
-
-                if (nonOrth == nNonOrthCorr)
                 {
                     phi -= pEqn.flux();
                 }

applications/solvers/engine/icoDyMEngineFoam/readControls.H

@@ -1,8 +1,4 @@
 #   include "readTimeControls.H"
-#   include "readPISOControls.H"
 
-    bool checkMeshCourantNo = false;
+    checkMeshCourantNo =
-    if (piso.found("checkMeshCourantNo"))
+        piso.dict().lookupOrDefault("checkMeshCourantNo", false);
-    {
-        checkMeshCourantNo = Switch(piso.lookup("checkMeshCourantNo"));
-    }

applications/solvers/engine/sonicTurbDyMEngineFoam/sonicTurbDyMEngineFoam.C

@@ -52,6 +52,7 @@ int main(int argc, char *argv[])
 #   include "readPIMPLEControls.H"
 #   include "createFields.H"
 #   include "initContinuityErrs.H"
+#   include "createTimeControls.H"
 #   include "readEngineTimeControls.H"
 #   include "CourantNo.H"
 #   include "setInitialDeltaT.H"

applications/solvers/engine/turbDyMEngineFoam/correctPhi.H

@@ -24,7 +24,9 @@
         pcorrTypes
     );
 
-    phi = (fvc::interpolate(U) & mesh.Sf());
+    phi = fvc::interpolate(U) & mesh.Sf();
+
+    mesh.schemesDict().setFluxRequired(pcorr.name());
 
     for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {

applications/solvers/engine/turbDyMEngineFoam/createControls.H

@@ -0,0 +1,8 @@
+#   include "createTimeControls.H"
+#   include "readPISOControls.H"
+
+    bool checkMeshCourantNo = false;
+    if (piso.found("checkMeshCourantNo"))
+    {
+        checkMeshCourantNo = Switch(piso.lookup("checkMeshCourantNo"));
+    }

applications/solvers/engine/turbDyMEngineFoam/readControls.H

@@ -1,7 +1,7 @@
 #   include "readTimeControls.H"
 #   include "readPISOControls.H"
 
-    bool checkMeshCourantNo = false;
+    checkMeshCourantNo = false;
     if (piso.found("checkMeshCourantNo"))
     {
         checkMeshCourantNo = Switch(piso.lookup("checkMeshCourantNo"));

applications/solvers/engine/turbDyMEngineFoam/turbDyMEngineFoam.C

@@ -47,6 +47,7 @@ int main(int argc, char *argv[])
 #   include "createDynamicFvMesh.H"
 #   include "initContinuityErrs.H"
 #   include "createFields.H"
+#   include "createControls.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

applications/solvers/finiteArea/liquidFilmFoam/liquidFilmFoam.C

@@ -49,6 +49,7 @@ int main(int argc, char *argv[])
 #   include "readTransportProperties.H"
 #   include "createFaFields.H"
 #   include "createFvFields.H"
+#   include "createTimeControls.H"
 
     Info << "\nStarting time loop\n" << endl;
 

applications/solvers/heatTransfer/buoyantBoussinesqPisoFoam/buoyantBoussinesqPisoFoam.C

@@ -59,7 +59,7 @@ int main(int argc, char *argv[])
 #   include "readGravitationalAcceleration.H"
 #   include "createFields.H"
 #   include "initContinuityErrs.H"
-#   include "readTimeControls.H"
+#   include "createTimeControls.H"
 #   include "CourantNo.H"
 #   include "setInitialDeltaT.H"
 

applications/solvers/heatTransfer/buoyantPisoFoam/buoyantPisoFoam.C

@@ -48,7 +48,7 @@ int main(int argc, char *argv[])
 #   include "readGravitationalAcceleration.H"
 #   include "createFields.H"
 #   include "initContinuityErrs.H"
-#   include "readTimeControls.H"
+#   include "createTimeControls.H"
 #   include "compressibleCourantNo.H"
 #   include "setInitialDeltaT.H"
 

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C

@@ -55,7 +55,7 @@ int main(int argc, char *argv[])
 
     #include "initContinuityErrs.H"
 
-    #include "readTimeControls.H"
+    #include "createTimeControls.H"
     #include "readSolidTimeControls.H"
 
 

applications/solvers/immersedBoundary/icoDyMIbFoam/createFields.H

@@ -17,6 +17,7 @@
         transportProperties.lookup("nu")
     );
 
+
     Info<< "Reading field p\n" << endl;
     volScalarField p
     (
@@ -31,6 +32,7 @@
         mesh
     );
 
+
     Info<< "Reading field U\n" << endl;
     volVectorField U
     (
@@ -47,6 +49,8 @@
 
 #   include "createPhi.H"
 
+
     label pRefCell = 0;
     scalar pRefValue = 0.0;
     setRefCell(p, mesh.solutionDict().subDict("PIMPLE"), pRefCell, pRefValue);
+    mesh.schemesDict().setFluxRequired(p.name());

applications/solvers/immersedBoundary/icoDyMIbFoam/icoDyMIbFoam.C

@@ -38,6 +38,7 @@ Author
 #include "dynamicFvMesh.H"
 #include "immersedBoundaryFvPatch.H"
 #include "immersedBoundaryAdjustPhi.H"
+#include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -46,6 +47,9 @@ int main(int argc, char *argv[])
 #   include "setRootCase.H"
 #   include "createTime.H"
 #   include "createDynamicFvMesh.H"
+
+    pimpleControl pimple(mesh);
+
 #   include "createFields.H"
 #   include "initContinuityErrs.H"
 
@@ -68,12 +72,10 @@ int main(int argc, char *argv[])
         // Make the fluxes relative to the mesh motion
         fvc::makeRelative(phi, U);
 
-#       include "readPIMPLEControls.H"
 #       include "CourantNo.H"
 
         // Pressure-velocity corrector
-        int oCorr = 0;
+        while (pimple.loop())
-        do
         {
             fvVectorMatrix UEqn
             (
@@ -82,10 +84,13 @@ int main(int argc, char *argv[])
               - fvm::laplacian(nu, U)
             );
 
+            if (pimple.momentumPredictor())
+            {
                 solve(UEqn == -fvc::grad(p));
+            }
 
             // --- PISO loop
-            for (int corr = 0; corr < nCorr; corr++)
+            while (pimple.correct())
             {
                 volScalarField rUA = 1.0/UEqn.A();
 
@@ -99,7 +104,8 @@ int main(int argc, char *argv[])
                 immersedBoundaryAdjustPhi(phi, U);
                 adjustPhi(phi, U, p);
 
-                for (int nonOrth = 0; nonOrth <= nNonOrthCorr; nonOrth++)
+                // Non-orthogonal pressure corrector loop
+                while (pimple.correctNonOrthogonal())
                 {
                     fvScalarMatrix pEqn
                     (
@@ -107,9 +113,15 @@ int main(int argc, char *argv[])
                     );
 
                     pEqn.setReference(pRefCell, pRefValue);
-                    pEqn.solve();
+                    pEqn.solve
+                    (
+                        mesh.solutionDict().solver
+                        (
+                            p.select(pimple.finalInnerIter())
+                        )
+                    );
 
-                    if (nonOrth == nNonOrthCorr)
+                    if (pimple.finalNonOrthogonalIter())
                     {
                         phi -= pEqn.flux();
                     }
@@ -123,7 +135,7 @@ int main(int argc, char *argv[])
                 U -= rUA*fvc::grad(p);
                 U.correctBoundaryConditions();
             }
-        } while (++oCorr < nOuterCorr);
+        }
 
         runTime.write();
 

applications/solvers/immersedBoundary/icoIbFoam/createFields.H

@@ -17,6 +17,7 @@
         transportProperties.lookup("nu")
     );
 
+
     Info<< "Reading field p\n" << endl;
     volScalarField p
     (
@@ -46,10 +47,10 @@
         mesh
     );
 
-
 #   include "createPhi.H"
 
 
     label pRefCell = 0;
     scalar pRefValue = 0.0;
     setRefCell(p, mesh.solutionDict().subDict("PIMPLE"), pRefCell, pRefValue);
+    mesh.schemesDict().setFluxRequired(p.name());

applications/solvers/immersedBoundary/icoIbFoam/icoIbFoam.C

@@ -36,6 +36,7 @@ Author
 #include "fvCFD.H"
 #include "immersedBoundaryFvPatch.H"
 #include "immersedBoundaryAdjustPhi.H"
+#include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -45,6 +46,9 @@ int main(int argc, char *argv[])
 
 #   include "createTime.H"
 #   include "createMesh.H"
+
+    pimpleControl pimple(mesh);
+
 #   include "createIbMasks.H"
 #   include "createFields.H"
 #   include "initContinuityErrs.H"
@@ -57,12 +61,10 @@ int main(int argc, char *argv[])
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-#       include "readPIMPLEControls.H"
 #       include "CourantNo.H"
 
         // Pressure-velocity corrector
-        int oCorr = 0;
+        while (pimple.loop())
-        do
         {
             fvVectorMatrix UEqn
             (
@@ -71,10 +73,13 @@ int main(int argc, char *argv[])
               - fvm::laplacian(nu, U)
             );
 
+            if (pimple.momentumPredictor())
+            {
                 solve(UEqn == -fvc::grad(p));
+            }
 
             // --- PISO loop
-            for (int corr = 0; corr < nCorr; corr++)
+            while (pimple.correct())
             {
                 volScalarField rUA = 1.0/UEqn.A();
 
@@ -88,7 +93,8 @@ int main(int argc, char *argv[])
                 immersedBoundaryAdjustPhi(phi, U);
                 adjustPhi(phi, U, p);
 
-                for (int nonOrth = 0; nonOrth <= nNonOrthCorr; nonOrth++)
+                // Non-orthogonal pressure corrector loop
+                while (pimple.correctNonOrthogonal())
                 {
                     fvScalarMatrix pEqn
                     (
@@ -96,9 +102,15 @@ int main(int argc, char *argv[])
                     );
 
                     pEqn.setReference(pRefCell, pRefValue);
-                    pEqn.solve();
+                    pEqn.solve
+                    (
+                        mesh.solutionDict().solver
+                        (
+                            p.select(pimple.finalInnerIter())
+                        )
+                    );
 
-                    if (nonOrth == nNonOrthCorr)
+                    if (pimple.finalNonOrthogonalIter())
                     {
                         phi -= pEqn.flux();
                     }
@@ -109,7 +121,7 @@ int main(int argc, char *argv[])
                 U -= rUA*fvc::grad(p);
                 U.correctBoundaryConditions();
             }
-        } while (++oCorr < nOuterCorr);
+        }
 
         runTime.write();
 

applications/solvers/immersedBoundary/interIbFoam/UEqn.H

@@ -16,7 +16,7 @@
 
     UEqn.relax();
 
-    if (momentumPredictor)
+    if (pimple.momentumPredictor())
     {
         solve
         (

applications/solvers/immersedBoundary/interIbFoam/correctPhi.H

@@ -36,7 +36,9 @@
     immersedBoundaryAdjustPhi(phi, U);
     adjustPhi(phi, U, pcorr);
 
-    for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    mesh.schemesDict().setFluxRequired(pcorr.name());
+
+    while(pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pcorrEqn
         (
@@ -46,7 +48,7 @@
         pcorrEqn.setReference(pdRefCell, pdRefValue);
         pcorrEqn.solve();
 
-        if (nonOrth == nNonOrthCorr)
+        if (pimple.finalNonOrthogonalIter())
         {
             phi -= pcorrEqn.flux();
         }

applications/solvers/immersedBoundary/interIbFoam/createFields.H

@@ -128,6 +128,8 @@
         );
     }
 
+    mesh.schemesDict().setFluxRequired(pd.name());
+
     // Construct interface from alpha1 distribution
     interfaceProperties interface(alpha1, U, twoPhaseProperties);
 

applications/solvers/immersedBoundary/interIbFoam/interIbFoam.C

@@ -45,6 +45,7 @@ Author
 #include "interfaceProperties.H"
 #include "twoPhaseMixture.H"
 #include "turbulenceModel.H"
+#include "pimpleControl.H"
 
 #include "immersedBoundaryFvPatch.H"
 #include "immersedBoundaryAdjustPhi.H"
@@ -56,12 +57,14 @@ int main(int argc, char *argv[])
 #   include "setRootCase.H"
 #   include "createTime.H"
 #   include "createMesh.H"
+
+    pimpleControl pimple(mesh);
+
 #   include "readGravitationalAcceleration.H"
-#   include "readPIMPLEControls.H"
 #   include "initContinuityErrs.H"
 #   include "createFields.H"
 #   include "createIbMasks.H"
-#   include "readTimeControls.H"
+#   include "createTimeControls.H"
 #   include "correctPhi.H"
 #   include "CourantNo.H"
 #   include "setInitialDeltaT.H"
@@ -72,7 +75,6 @@ int main(int argc, char *argv[])
 
     while (runTime.run())
     {
-#       include "readPIMPLEControls.H"
 #       include "readTimeControls.H"
 #       include "immersedBoundaryCourantNo.H"
 #       include "setDeltaT.H"
@@ -82,8 +84,7 @@ int main(int argc, char *argv[])
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
         // Pressure-velocity corrector
-        int oCorr = 0;
+        while (pimple.loop())
-        do
         {
             twoPhaseProperties.correct();
 
@@ -92,7 +93,7 @@ int main(int argc, char *argv[])
 #           include "UEqn.H"
 
             // --- PISO loop
-            for (int corr = 0; corr < nCorr; corr++)
+            while (pimple.correct())
             {
 #               include "pEqn.H"
             }
@@ -117,7 +118,7 @@ int main(int argc, char *argv[])
             }
 
             turbulence->correct();
-        } while (++oCorr < nOuterCorr);
+        }
 
         runTime.write();
 

applications/solvers/immersedBoundary/interIbFoam/limitU.H

@@ -1,5 +1,6 @@
 {
-    scalar limitMagU = readScalar(pimple.lookup("limitMagU"));
+    scalar limitMagU =
+        readScalar(mesh.solutionDict().subDict("PIMPLE").lookup("limitMagU"));
 
     volScalarField magU(mag(U));
 

applications/solvers/immersedBoundary/interIbFoam/pEqn.H

@@ -1,5 +1,5 @@
 {
-    if (nOuterCorr != 1)
+    if (!pimple.firstIter())
     {
         pd.storePrevIter();
     }
@@ -29,7 +29,7 @@
         )*rUAf*mesh.magSf();
 
 
-    for(int nonOrth = 0; nonOrth <= nNonOrthCorr; nonOrth++)
+    while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pdEqn
         (
@@ -38,23 +38,20 @@
 
         pdEqn.setReference(pdRefCell, pdRefValue);
 
-        if (corr == nCorr - 1 && nonOrth == nNonOrthCorr)
+        pdEqn.solve
-        {
+        (
-            pdEqn.solve(mesh.solutionDict().solver(pd.name() + "Final"));
+            mesh.solutionDict().solver(pd.select(pimple.finalInnerIter()))
-        }
+        );
-        else
-        {
-            pdEqn.solve(mesh.solutionDict().solver(pd.name()));
-        }
 
-        if (nonOrth == nNonOrthCorr)
+
+        if (pimple.finalNonOrthogonalIter())
         {
             phi -= pdEqn.flux();
         }
     }
 
     // Explicitly relax pressure except for last corrector
-    if (oCorr != nOuterCorr - 1)
+    if (!pimple.finalIter())
     {
         pd.relax();
     }

applications/solvers/immersedBoundary/porousSimpleIbFoam/UEqn.H

@@ -9,9 +9,4 @@
 
     pZones.addResistance(UEqn());
 
-    eqnResidual = solve
+    solve(UEqn() == -fvc::grad(p));
-    (
-        UEqn() == -fvc::grad(p)
-    ).initialResidual();
-
-    maxResidual = max(eqnResidual, maxResidual);

applications/solvers/immersedBoundary/porousSimpleIbFoam/convergenceCheck.H

@@ -1,9 +0,0 @@
-// check convergence
-
-if (maxResidual < convergenceCriterion)
-{
-    Info<< "reached convergence criterion: " << convergenceCriterion << endl;
-    runTime.writeAndEnd();
-    Info<< "latestTime = " << runTime.timeName() << endl;
-}
-

applications/solvers/immersedBoundary/porousSimpleIbFoam/createFields.H

@@ -32,7 +32,7 @@
     label pRefCell = 0;
     scalar pRefValue = 0.0;
     setRefCell(p, mesh.solutionDict().subDict("SIMPLE"), pRefCell, pRefValue);
-
+    mesh.schemesDict().setFluxRequired(p.name());
 
     singlePhaseTransportModel laminarTransport(U, phi);
 

applications/solvers/immersedBoundary/porousSimpleIbFoam/initConvergenceCheck.H

@@ -1,7 +0,0 @@
-// initialize values for convergence checks
-
-    scalar eqnResidual = 1, maxResidual = 0;
-    scalar convergenceCriterion = 0;
-
-    simple.readIfPresent("convergence", convergenceCriterion);
-

applications/solvers/immersedBoundary/porousSimpleIbFoam/pEqn.H

@@ -13,7 +13,7 @@
     adjustPhi(phi, U, p);
 
     // Non-orthogonal pressure corrector loop
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    while (simple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
@@ -22,18 +22,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();
-        }
 
-        if (nonOrth == nNonOrthCorr)
+        if (simple.finalNonOrthogonalIter())
         {
             phi -= pEqn.flux();
         }

applications/solvers/immersedBoundary/porousSimpleIbFoam/porousSimpleIbFoam.C

@@ -37,6 +37,7 @@ Author
 #include "singlePhaseTransportModel.H"
 #include "RASModel.H"
 #include "porousZones.H"
+#include "simpleControl.H"
 
 #include "immersedBoundaryFvPatch.H"
 #include "immersedBoundaryAdjustPhi.H"
@@ -48,6 +49,9 @@ int main(int argc, char *argv[])
 #   include "setRootCase.H"
 #   include "createTime.H"
 #   include "createMesh.H"
+
+    simpleControl simple(mesh);
+
 #   include "createIbMasks.H"
 #   include "createFields.H"
 #   include "initContinuityErrs.H"
@@ -60,9 +64,6 @@ int main(int argc, char *argv[])
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-#       include "readSIMPLEControls.H"
-#       include "initConvergenceCheck.H"
-
         p.storePrevIter();
 
         // Pressure-velocity SIMPLE corrector
@@ -78,8 +79,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;

applications/solvers/immersedBoundary/potentialIbFoam/Make/options

@@ -1,6 +1,8 @@
 EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/dynamicMesh/dynamicFvMesh/lnInclude \
+    -I$(LIB_SRC)/dynamicMesh/dynamicMesh/lnInclude \
     -I$(LIB_SRC)/immersedBoundary/immersedBoundary/lnInclude
 
 EXE_LIBS = \
@@ -9,5 +11,6 @@ EXE_LIBS = \
     -lsurfMesh \
     -lsampling \
     -ldynamicMesh \
+    -ldynamicFvMesh \
     -limmersedBoundary \
     -llduSolvers

applications/solvers/immersedBoundary/potentialIbFoam/createFields.H

@@ -47,3 +47,4 @@
     label pRefCell = 0;
     scalar pRefValue = 0.0;
     setRefCell(p, mesh.solutionDict().subDict("SIMPLE"), pRefCell, pRefValue);
+    mesh.schemesDict().setFluxRequired(p.name());

applications/solvers/immersedBoundary/potentialIbFoam/potentialIbFoam.C

@@ -32,6 +32,7 @@ Description
 #include "fvCFD.H"
 #include "immersedBoundaryFvPatch.H"
 #include "immersedBoundaryAdjustPhi.H"
+#include "simpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -43,16 +44,18 @@ int main(int argc, char *argv[])
 
 #   include "createTime.H"
 #   include "createMesh.H"
+
+    simpleControl simple(mesh);
+
 #   include "createIbMasks.H"
 #   include "createFields.H"
-#   include "readSIMPLEControls.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< nl << "Calculating potential flow" << endl;
 
     // Do correctors over the complete set
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    while (simple.correctNonOrthogonal())
     {
         phi = faceIbMask*(linearInterpolate(U) & mesh.Sf());
 
@@ -86,7 +89,7 @@ int main(int argc, char *argv[])
         // Correct the flux
         phi -= pEqn.flux();
 
-        if (nonOrth != nNonOrthCorr)
+        if (!simple.finalNonOrthogonalIter())
         {
             p.relax();
         }

applications/solvers/immersedBoundary/simpleIbFoam/UEqn.H

@@ -7,9 +7,4 @@
 
     UEqn().relax();
 
-    eqnResidual = solve
+    solve(UEqn() == -fvc::grad(p));
-    (
-        UEqn() == -fvc::grad(p)
-    ).initialResidual();
-
-    maxResidual = max(eqnResidual, maxResidual);

applications/solvers/immersedBoundary/simpleIbFoam/convergenceCheck.H

@@ -1,9 +0,0 @@
-// check convergence
-
-if (maxResidual < convergenceCriterion)
-{
-    Info<< "reached convergence criterion: " << convergenceCriterion << endl;
-    runTime.writeAndEnd();
-    Info<< "latestTime = " << runTime.timeName() << endl;
-}
-

applications/solvers/immersedBoundary/simpleIbFoam/createFields.H

@@ -32,7 +32,7 @@
     label pRefCell = 0;
     scalar pRefValue = 0.0;
     setRefCell(p, mesh.solutionDict().subDict("SIMPLE"), pRefCell, pRefValue);
-
+    mesh.schemesDict().setFluxRequired(p.name());
 
     singlePhaseTransportModel laminarTransport(U, phi);
 

applications/solvers/immersedBoundary/simpleIbFoam/initConvergenceCheck.H

@@ -1,7 +0,0 @@
-// initialize values for convergence checks
-
-    scalar eqnResidual = 1, maxResidual = 0;
-    scalar convergenceCriterion = 0;
-
-    simple.readIfPresent("convergence", convergenceCriterion);
-

applications/solvers/immersedBoundary/simpleIbFoam/pEqn.H

@@ -13,7 +13,7 @@
     adjustPhi(phi, U, p);
 
     // Non-orthogonal pressure corrector loop
-    for (int nonOrth = 0; nonOrth <= nNonOrthCorr; nonOrth++)
+    while (simple.correctNonOrthogonal())
     {
         fvScalarMatrix pEqn
         (
@@ -22,18 +22,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();
-        }
 
-        if (nonOrth == nNonOrthCorr)
+        if (simple.finalNonOrthogonalIter())
         {
             phi -= pEqn.flux();
         }

applications/solvers/immersedBoundary/simpleIbFoam/simpleIbFoam.C

@@ -36,6 +36,7 @@ Author
 #include "fvCFD.H"
 #include "singlePhaseTransportModel.H"
 #include "RASModel.H"
+#include "simpleControl.H"
 
 #include "immersedBoundaryFvPatch.H"
 #include "immersedBoundaryAdjustPhi.H"
@@ -47,6 +48,9 @@ int main(int argc, char *argv[])
 #   include "setRootCase.H"
 #   include "createTime.H"
 #   include "createMesh.H"
+
+    simpleControl simple(mesh);
+
 #   include "createIbMasks.H"
 #   include "createFields.H"
 #   include "initContinuityErrs.H"
@@ -59,9 +63,6 @@ int main(int argc, char *argv[])
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-#       include "readSIMPLEControls.H"
-#       include "initConvergenceCheck.H"
-
         p.storePrevIter();
 
         // Pressure-velocity SIMPLE corrector
@@ -77,8 +78,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;

applications/solvers/incompressible/MRFSimpleFoam/MRFSimpleFoam.C

@@ -34,6 +34,7 @@ Description
 #include "incompressible/singlePhaseTransportModel/singlePhaseTransportModel.H"
 #include "incompressible/RAS/RASModel/RASModel.H"
 #include "MRFZones.H"
+#include "simpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -44,6 +45,9 @@ int main(int argc, char *argv[])
 
 #   include "createTime.H"
 #   include "createMesh.H"
+
+    simpleControl simple(mesh);
+
 #   include "createFields.H"
 #   include "initContinuityErrs.H"
 
@@ -55,8 +59,6 @@ int main(int argc, char *argv[])
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-#       include "readSIMPLEControls.H"
-
         p.storePrevIter();
 
         // Pressure-velocity SIMPLE corrector
@@ -83,7 +85,7 @@ int main(int argc, char *argv[])
             adjustPhi(phi, U, p);
 
             // Non-orthogonal pressure corrector loop
-            for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+            while (simple.correctNonOrthogonal())
             {
                 fvScalarMatrix pEqn
                 (
@@ -93,7 +95,7 @@ int main(int argc, char *argv[])
                 pEqn.setReference(pRefCell, pRefValue);
                 pEqn.solve();
 
-                if (nonOrth == nNonOrthCorr)
+                if (simple.finalNonOrthogonalIter())
                 {
                     phi -= pEqn.flux();
                 }

applications/solvers/incompressible/MRFSimpleFoam/createFields.H

@@ -32,7 +32,7 @@
     label pRefCell = 0;
     scalar pRefValue = 0.0;
     setRefCell(p, mesh.solutionDict().subDict("SIMPLE"), pRefCell, pRefValue);
-
+    mesh.schemesDict().setFluxRequired(p.name());
 
     singlePhaseTransportModel laminarTransport(U, phi);
 

applications/solvers/incompressible/channelFoam/channelFoam.C

@@ -39,6 +39,7 @@ Author
 #include "IFstream.H"
 #include "OFstream.H"
 #include "Random.H"
+#include "pisoControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -47,6 +48,9 @@ int main(int argc, char *argv[])
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createMesh.H"
+
+    pisoControl piso(mesh);
+
     #include "readTransportProperties.H"
     #include "createFields.H"
     #include "initContinuityErrs.H"
@@ -58,8 +62,6 @@ int main(int argc, char *argv[])
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        #include "readPISOControls.H"
-
         #include "CourantNo.H"
 
         sgsModel->correct();
@@ -76,7 +78,7 @@ int main(int argc, char *argv[])
         // Time derivative matrix
         fvVectorMatrix ddtUEqn(fvm::ddt(U));
 
-        if (momentumPredictor)
+        if (piso.momentumPredictor())
         {
             solve(ddtUEqn + HUEqn == -fvc::grad(p));
         }
@@ -87,14 +89,14 @@ int main(int argc, char *argv[])
 
         // --- PISO loop
 
-        for (int corr = 0; corr < nCorr; corr++)
+        while (piso.correct())
         {
             U = HUEqn.H()/aU;
             phi = (fvc::interpolate(U) & mesh.Sf());
 
             adjustPhi(phi, U, p);
 
-            for (int nonOrth = 0; nonOrth <= nNonOrthCorr; nonOrth++)
+            while (piso.correctNonOrthogonal())
             {
                 fvScalarMatrix pEqn
                 (
@@ -102,17 +104,12 @@ int main(int argc, char *argv[])
                 );
 
                 pEqn.setReference(pRefCell, pRefValue);
+                pEqn.solve
+                (
+                    mesh.solutionDict().solver(p.select(piso.finalInnerIter()))
+                );
 
-                if (corr == nCorr-1 && nonOrth == nNonOrthCorr)
+                if (piso.finalNonOrthogonalIter())
-                {
-                    pEqn.solve(mesh.solutionDict().solver(p.name() + "Final"));
-                }
-                else
-                {
-                    pEqn.solve(mesh.solutionDict().solver(p.name()));
-                }
-
-                if (nonOrth == nNonOrthCorr)
                 {
                     phi -= pEqn.flux();
                 }

applications/solvers/incompressible/channelFoam/createFields.H

@@ -33,7 +33,7 @@
     label pRefCell = 0;
     scalar pRefValue = 0.0;
     setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue);
-
+    mesh.schemesDict().setFluxRequired(p.name());
 
     singlePhaseTransportModel laminarTransport(U, phi);
 

applications/solvers/incompressible/consistentIcoFoam/Make/files

@@ -1,3 +0,0 @@
-consistentIcoFoam.C
-
-EXE = $(FOAM_APPBIN)/consistentIcoFoam

applications/solvers/incompressible/consistentIcoFoam/Make/options

@@ -1,6 +0,0 @@
-EXE_INC = \
-    -I$(LIB_SRC)/finiteVolume/lnInclude
-
-EXE_LIBS = \
-    -lfiniteVolume \
-    -llduSolvers

applications/solvers/incompressible/consistentIcoFoam/consistentIcoFoam.C

@@ -1,125 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | foam-extend: Open Source CFD
-   \\    /   O peration     | Version:     3.2
-    \\  /    A nd           | Web:         http://www.foam-extend.org
-     \\/     M anipulation  | For copyright notice see file Copyright
--------------------------------------------------------------------------------
-License
-    This file is part of foam-extend.
-
-    foam-extend is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by the
-    Free Software Foundation, either version 3 of the License, or (at your
-    option) any later version.
-
-    foam-extend is distributed in the hope that it will be useful, but
-    WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
-    General Public License for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
-
-Application
-    consistentIcoFoam
-
-Description
-    Transient solver for incompressible, laminar flow of Newtonian fluids.
-    Consistent formulation without time-step and relaxation dependence by Jasak
-
-Author
-    Hrvoje Jasak, Wikki Ltd.  All rights reserved
-
-\*---------------------------------------------------------------------------*/
-
-#include "fvCFD.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-int main(int argc, char *argv[])
-{
-#   include "setRootCase.H"
-
-#   include "createTime.H"
-#   include "createMesh.H"
-#   include "createFields.H"
-#   include "initContinuityErrs.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-    Info<< "\nStarting time loop\n" << endl;
-
-    while (runTime.loop())
-    {
-        Info<< "Time = " << runTime.timeName() << nl << endl;
-
-#       include "readPISOControls.H"
-#       include "CourantNo.H"
-
-        // Convection-diffusion matrix
-        fvVectorMatrix HUEqn
-        (
-            fvm::div(phi, U)
-          - fvm::laplacian(nu, U)
-        );
-
-        // ddt matrix
-        fvVectorMatrix ddtUEqn(fvm::ddt(U));
-
-        solve(ddtUEqn + HUEqn == -fvc::grad(p));
-
-        // Prepare clean Ap without time derivative contribution
-        // HJ, 26/Oct/2015
-        volScalarField aU = HUEqn.A();
-
-        // --- PISO loop
-
-        for (int corr = 0; corr < nCorr; corr++)
-        {
-            U = HUEqn.H()/aU;
-            phi = (fvc::interpolate(U) & mesh.Sf());
-
-            adjustPhi(phi, U, p);
-
-            for (int nonOrth = 0; nonOrth <= nNonOrthCorr; nonOrth++)
-            {
-                fvScalarMatrix pEqn
-                (
-                    fvm::laplacian(1/aU, p) == fvc::div(phi)
-                );
-
-                pEqn.setReference(pRefCell, pRefValue);
-                pEqn.solve();
-
-                if (nonOrth == nNonOrthCorr)
-                {
-                    phi -= pEqn.flux();
-                }
-            }
-
-#           include "continuityErrs.H"
-
-            // Note: cannot call H(U) here because the velocity is not complete
-            // HJ, 22/Jan/2016
-            U = 1.0/(aU + ddtUEqn.A())*
-                (
-                    U*aU - fvc::grad(p) + ddtUEqn.H()
-                );
-            U.correctBoundaryConditions();
-        }
-
-        runTime.write();
-
-        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
-            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
-            << nl << endl;
-    }
-
-    Info<< "End\n" << endl;
-
-    return 0;
-}
-
-
-// ************************************************************************* //