Merge branch 'master' of /home/noisy3/OpenFOAM/OpenFOAM-dev

ReleaseNotes-dev

@@ -120,6 +120,8 @@
     + reuses tracking state in interpolation - improves consistency and speed
 *** *Updated* split cyclics into two separate patches.
     See [[./doc/changes/splitCyclic.txt]]
+  * *New* cyclics (and all coupled patches) have optional
+    matchTolerance entry to allow looser area matching tolerance.
   * *Updated* interpolation (volPointInterpolation) now works without the
     globalPointPatch. Moving mesh cases can now be run non-parallel and
     continued in parallel and reconstructed without any limitation.
@@ -264,8 +266,6 @@
     internal field instead of boundary field.
     + =ensight=: new output format for all sampledSurfaces.
   + Function objects:
-    + =residualControl=: new function object to allow users to terminate steady
-    state calculations when the defined residual levels are achieved
     + =abortCalculation=: watches for presence of the named file in the
     $FOAM_CASE directory and aborts the calculation if it is present
     + =timeActivatedFileUpdate=: performs a file copy/replacement once a
@@ -312,8 +312,3 @@
 * Other
   + compilable with =clang=
     In your prefs.sh set the WM_COMPILER to Clang
-  + dlclose error: upon exit of e.g. paraFoam you can get the error
-    Inconsistency detected by ld.so: dl-close.c: 731: _dl_close: Assertion
-    `map->l_init_called' failed!
-    This seems to happen with a mix of system libraries and a thirdParty
-    compiler.

applications/solvers/DNS/dnsFoam/createFields.H

@@ -26,4 +26,4 @@
         mesh
     );
 
-#   include "createPhi.H"
+    #include "createPhi.H"

applications/solvers/DNS/dnsFoam/dnsFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -110,7 +110,12 @@ int main(int argc, char *argv[])
 
         if (runTime.outputTime())
         {
-            calcEk(U, K).write(runTime.timePath()/"Ek", runTime.graphFormat());
+            calcEk(U, K).write
+            (
+                runTime.path()/"graphs"/runTime.timeName(),
+                "Ek",
+                runTime.graphFormat()
+            );
         }
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"

applications/solvers/basic/laplacianFoam/laplacianFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -30,6 +30,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
+#include "simpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -41,17 +42,17 @@ int main(int argc, char *argv[])
     #include "createMesh.H"
     #include "createFields.H"
 
+    simpleControl simple(mesh);
+
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nCalculating temperature distribution\n" << endl;
 
-    while (runTime.loop())
+    while (simple.loop())
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        #include "readSIMPLEControls.H"
+        for (int nonOrth=0; nonOrth<=simple.nNonOrthCorr(); nonOrth++)
-
-        for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
         {
             solve
             (

applications/solvers/basic/scalarTransportFoam/scalarTransportFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -30,32 +30,30 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
+#include "simpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 int main(int argc, char *argv[])
 {
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "createFields.H"
 
-#   include "setRootCase.H"
+    simpleControl simple(mesh);
 
-#   include "createTime.H"
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-#   include "createMesh.H"
-#   include "createFields.H"
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nCalculating scalar transport\n" << endl;
 
-#   include "CourantNo.H"
+    #include "CourantNo.H"
 
-    while (runTime.loop())
+    while (simple.loop())
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-#       include "readSIMPLEControls.H"
+        for (int nonOrth=0; nonOrth<=simple.nNonOrthCorr(); nonOrth++)
-
-        for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
         {
             solve
             (

applications/solvers/combustion/PDRFoam/PDRFoam.C

@@ -53,6 +53,18 @@ Description
     PDR (porosity/distributed resistance) modelling is included to handle
     regions containing blockages which cannot be resolved by the mesh.
 
+    The fields used by this solver are:
+
+    betav:  Volume porosity
+    Lobs:   Average diameter of obstacle in cell (m)
+    Aw:     Obstacle surface area per unit volume (1/m)
+    CR:     Drag tensor (1/m)
+    CT:     Turbulence generation parameter (1/m)
+    Nv:     Number of obstacles in cell per unit volume (m^-2)
+    nsv:    Tensor whose diagonal indicates the number to substract from
+            Nv to get the number of obstacles crossing the flow in each
+            direction.
+
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
@@ -64,6 +76,7 @@ Description
 #include "ignition.H"
 #include "Switch.H"
 #include "bound.H"
+#include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -81,6 +94,8 @@ int main(int argc, char *argv[])
     #include "compressibleCourantNo.H"
     #include "setInitialDeltaT.H"
 
+    pimpleControl pimple(mesh);
+
     scalar StCoNum = 0.0;
 
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -90,7 +105,6 @@ int main(int argc, char *argv[])
     while (runTime.run())
     {
         #include "readTimeControls.H"
-        #include "readPISOControls.H"
         #include "compressibleCourantNo.H"
         #include "setDeltaT.H"
 
@@ -98,26 +112,34 @@ int main(int argc, char *argv[])
         Info<< "\n\nTime = " << runTime.timeName() << endl;
 
         #include "rhoEqn.H"
-        #include "UEqn.H"
 
-        // --- PISO loop
+        // --- Pressure-velocity PIMPLE corrector loop
-        for (int corr=1; corr<=nCorr; corr++)
+        for (pimple.start(); pimple.loop(); pimple++)
         {
-            #include "bEqn.H"
+            #include "UEqn.H"
-            #include "ftEqn.H"
-            #include "huEqn.H"
-            #include "hEqn.H"
 
-            if (!ign.ignited())
+            // --- PISO loop
+            for (int corr=1; corr<=pimple.nCorr(); corr++)
             {
-                hu == h;
+                #include "bEqn.H"
+                #include "ftEqn.H"
+                #include "huEqn.H"
+                #include "hEqn.H"
+
+                if (!ign.ignited())
+                {
+                    hu == h;
+                }
+
+                #include "pEqn.H"
             }
 
-            #include "pEqn.H"
+            if (pimple.turbCorr())
+            {
+                turbulence->correct();
+            }
         }
 
-        turbulence->correct();
-
         runTime.write();
 
         Info<< "\nExecutionTime = "

applications/solvers/combustion/PDRFoam/UEqn.H

@@ -9,7 +9,7 @@
 
     volSymmTensorField invA(inv(I*UEqn.A() + drag->Dcu()));
 
-    if (momentumPredictor)
+    if (pimple.momentumPredictor())
     {
         U = invA & (UEqn.H() - betav*fvc::grad(p));
         U.correctBoundaryConditions();

applications/solvers/combustion/PDRFoam/XiModels/algebraic/algebraic.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2011 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

applications/solvers/combustion/PDRFoam/laminarFlameSpeed/SCOPE/SCOPELaminarFlameSpeed.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2011 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

applications/solvers/combustion/PDRFoam/pEqn.H

@@ -3,7 +3,7 @@ rho = thermo.rho();
 volScalarField rAU(1.0/UEqn.A());
 U = invA & UEqn.H();
 
-if (transonic)
+if (pimple.transonic())
 {
     surfaceScalarField phid
     (
@@ -15,7 +15,7 @@ if (transonic)
         )
     );
 
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
     {
         fvScalarMatrix pEqn
         (
@@ -24,9 +24,12 @@ if (transonic)
           - fvm::laplacian(rho*invA, p)
         );
 
-        pEqn.solve();
+        pEqn.solve
+        (
+            mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
+        );
 
-        if (nonOrth == nNonOrthCorr)
+        if (nonOrth == pimple.nNonOrthCorr())
         {
             phi == pEqn.flux();
         }
@@ -41,7 +44,7 @@ else
           + fvc::ddtPhiCorr(rAU, rho, U, phi)
         );
 
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
     {
         fvScalarMatrix pEqn
         (
@@ -50,9 +53,12 @@ else
           - fvm::laplacian(rho*invA, p)
         );
 
-        pEqn.solve();
+        pEqn.solve
+        (
+            mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
+        );
 
-        if (nonOrth == nNonOrthCorr)
+        if (nonOrth == pimple.nNonOrthCorr())
         {
             phi += pEqn.flux();
         }

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -55,6 +55,7 @@ Description
 #include "laminarFlameSpeed.H"
 #include "ignition.H"
 #include "Switch.H"
+#include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -72,6 +73,8 @@ int main(int argc, char *argv[])
     #include "compressibleCourantNo.H"
     #include "setInitialDeltaT.H"
 
+    pimpleControl pimple(mesh);
+
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
@@ -79,17 +82,17 @@ int main(int argc, char *argv[])
     while (runTime.run())
     {
         #include "readTimeControls.H"
-        #include "readPISOControls.H"
         #include "compressibleCourantNo.H"
         #include "setDeltaT.H"
 
         runTime++;
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
+        #include "rhoEqn.H"
+
         // --- Pressure-velocity PIMPLE corrector loop
-        for (int oCorr=0; oCorr<nOuterCorr; oCorr++)
+        for (pimple.start(); pimple.loop(); pimple++)
         {
-            #include "rhoEqn.H"
             #include "UEqn.H"
 
             #include "ftEqn.H"
@@ -103,12 +106,15 @@ int main(int argc, char *argv[])
             }
 
             // --- PISO loop
-            for (int corr=1; corr<=nCorr; corr++)
+            for (int corr=0; corr<pimple.nCorr(); corr++)
             {
                 #include "pEqn.H"
             }
 
-            turbulence->correct();
+            if (pimple.turbCorr())
+            {
+                turbulence->correct();
+            }
         }
 
         rho = thermo.rho();

applications/solvers/combustion/XiFoam/pEqn.H

@@ -3,7 +3,7 @@ rho = thermo.rho();
 volScalarField rAU(1.0/UEqn.A());
 U = rAU*UEqn.H();
 
-if (transonic)
+if (pimple.transonic())
 {
     surfaceScalarField phid
     (
@@ -15,7 +15,7 @@ if (transonic)
         )
     );
 
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
     {
         fvScalarMatrix pEqn
         (
@@ -24,9 +24,12 @@ if (transonic)
           - fvm::laplacian(rho*rAU, p)
         );
 
-        pEqn.solve();
+        pEqn.solve
+        (
+            mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
+        );
 
-        if (nonOrth == nNonOrthCorr)
+        if (nonOrth == pimple.nNonOrthCorr())
         {
             phi == pEqn.flux();
         }
@@ -41,7 +44,7 @@ else
           + fvc::ddtPhiCorr(rAU, rho, U, phi)
         );
 
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
     {
         fvScalarMatrix pEqn
         (
@@ -50,9 +53,12 @@ else
           - fvm::laplacian(rho*rAU, p)
         );
 
-        pEqn.solve();
+        pEqn.solve
+        (
+            mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
+        );
 
-        if (nonOrth == nNonOrthCorr)
+        if (nonOrth == pimple.nNonOrthCorr())
         {
             phi += pEqn.flux();
         }

applications/solvers/combustion/coldEngineFoam/coldEngineFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -35,6 +35,7 @@ Description
 #include "basicPsiThermo.H"
 #include "turbulenceModel.H"
 #include "OFstream.H"
+#include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -51,13 +52,14 @@ int main(int argc, char *argv[])
     #include "setInitialDeltaT.H"
     #include "startSummary.H"
 
+    pimpleControl pimple(mesh);
+
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
 
     while (runTime.run())
     {
-        #include "readPISOControls.H"
         #include "readEngineTimeControls.H"
         #include "compressibleCourantNo.H"
         #include "setDeltaT.H"
@@ -71,16 +73,23 @@ int main(int argc, char *argv[])
 
         #include "rhoEqn.H"
 
-        #include "UEqn.H"
+        // --- Pressure-velocity PIMPLE corrector loop
-
+        for (pimple.start(); pimple.loop(); pimple++)
-        // --- PISO loop
-        for (int corr=1; corr<=nCorr; corr++)
         {
-            #include "hEqn.H"
+            #include "UEqn.H"
-            #include "pEqn.H"
-        }
 
-        turbulence->correct();
+            // --- PISO loop
+            for (int corr=1; corr<=pimple.nCorr(); corr++)
+            {
+                #include "hEqn.H"
+                #include "pEqn.H"
+            }
+
+            if (pimple.turbCorr())
+            {
+                turbulence->correct();
+            }
+        }
 
         runTime.write();
 

applications/solvers/combustion/dieselEngineFoam/UEqn.H

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

applications/solvers/combustion/dieselEngineFoam/dieselEngineFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -43,6 +43,7 @@ Description
 #include "volPointInterpolation.H"
 #include "thermoPhysicsTypes.H"
 #include "mathematicalConstants.H"
+#include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -61,13 +62,14 @@ int main(int argc, char *argv[])
     #include "setInitialDeltaT.H"
     #include "startSummary.H"
 
+    pimpleControl pimple(mesh);
+
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
 
     while (runTime.run())
     {
-        #include "readPISOControls.H"
         #include "readEngineTimeControls.H"
         #include "compressibleCourantNo.H"
         #include "setDeltaT.H"
@@ -77,10 +79,6 @@ int main(int argc, char *argv[])
         Info<< "Crank angle = " << runTime.theta() << " CA-deg" << endl;
 
         mesh.move();
-        const_cast<volPointInterpolation&>
-        (
-            volPointInterpolation::New(mesh)
-        ).updateMesh();
 
         dieselSpray.evolve();
 
@@ -106,22 +104,26 @@ int main(int argc, char *argv[])
 
         chemistrySh = kappa*chemistry.Sh()();
 
-        #include "rhoEqn.H"
-        #include "UEqn.H"
 
-        for (label ocorr=1; ocorr <= nOuterCorr; ocorr++)
+        #include "rhoEqn.H"
+
+        for (pimple.start(); pimple.loop(); pimple++)
         {
+            #include "UEqn.H"
             #include "YEqn.H"
             #include "hsEqn.H"
 
             // --- PISO loop
-            for (int corr=1; corr<=nCorr; corr++)
+            for (int corr=1; corr<=pimple.nCorr(); corr++)
             {
                 #include "pEqn.H"
             }
-        }
 
-        turbulence->correct();
+            if (pimple.turbCorr())
+            {
+                turbulence->correct();
+            }
+        }
 
         #include "logSummary.H"
         #include "spraySummary.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++)
+    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
     {
         fvScalarMatrix pEqn
         (
@@ -23,9 +23,12 @@ if (transonic)
             Sevap
         );
 
-        pEqn.solve();
+        pEqn.solve
+        (
+            mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
+        );
 
-        if (nonOrth == nNonOrthCorr)
+        if (nonOrth == pimple.nNonOrthCorr())
         {
             phi == pEqn.flux();
         }
@@ -36,7 +39,7 @@ else
     phi = fvc::interpolate(rho)
          *((fvc::interpolate(U) & mesh.Sf()) - fvc::meshPhi(rho, U));
 
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
     {
         fvScalarMatrix pEqn
         (
@@ -47,9 +50,12 @@ else
             Sevap
         );
 
-        pEqn.solve();
+        pEqn.solve
+        (
+            mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
+        );
 
-        if (nonOrth == nNonOrthCorr)
+        if (nonOrth == pimple.nNonOrthCorr())
         {
             phi += pEqn.flux();
         }

applications/solvers/combustion/dieselFoam/dieselFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -41,6 +41,7 @@ Description
 #include "OFstream.H"
 #include "Switch.H"
 #include "mathematicalConstants.H"
+#include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -58,13 +59,14 @@ int main(int argc, char *argv[])
     #include "compressibleCourantNo.H"
     #include "setInitialDeltaT.H"
 
+    pimpleControl pimple(mesh);
+
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
 
     while (runTime.run())
     {
-        #include "readPISOControls.H"
         #include "compressibleCourantNo.H"
         #include "setDeltaT.H"
 
@@ -98,21 +100,25 @@ int main(int argc, char *argv[])
         chemistrySh = kappa*chemistry.Sh()();
 
         #include "rhoEqn.H"
-        #include "UEqn.H"
 
-        for (label ocorr=1; ocorr <= nOuterCorr; ocorr++)
+        // --- Pressure-velocity PIMPLE corrector loop
+        for (pimple.start(); pimple.loop(); pimple++)
         {
+            #include "UEqn.H"
             #include "YEqn.H"
             #include "hsEqn.H"
 
             // --- PISO loop
-            for (int corr=1; corr<=nCorr; corr++)
+            for (int corr=0; corr<pimple.nCorr(); corr++)
             {
                 #include "pEqn.H"
             }
-        }
 
-        turbulence->correct();
+            if (pimple.turbCorr())
+            {
+                turbulence->correct();
+            }
+        }
 
         #include "spraySummary.H"
 

applications/solvers/combustion/dieselFoam/pEqn.H

@@ -3,7 +3,7 @@ rho = thermo.rho();
 volScalarField rAU(1.0/UEqn.A());
 U = rAU*UEqn.H();
 
-if (transonic)
+if (pimple.transonic())
 {
     surfaceScalarField phid
     (
@@ -15,7 +15,7 @@ if (transonic)
         )
     );
 
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
     {
         fvScalarMatrix pEqn
         (
@@ -26,9 +26,12 @@ if (transonic)
             Sevap
         );
 
-        pEqn.solve();
+        pEqn.solve
+        (
+            mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
+        );
 
-        if (nonOrth == nNonOrthCorr)
+        if (nonOrth == pimple.nNonOrthCorr())
         {
             phi == pEqn.flux();
         }
@@ -43,7 +46,7 @@ else
           + fvc::ddtPhiCorr(rAU, rho, U, phi)
         );
 
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
     {
         fvScalarMatrix pEqn
         (
@@ -54,9 +57,12 @@ else
             Sevap
         );
 
-        pEqn.solve();
+        pEqn.solve
+        (
+            mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
+        );
 
-        if (nonOrth == nNonOrthCorr)
+        if (nonOrth == pimple.nNonOrthCorr())
         {
             phi += pEqn.flux();
         }

applications/solvers/combustion/engineFoam/UEqn.H

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

applications/solvers/combustion/engineFoam/engineFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -58,6 +58,7 @@ Description
 #include "Switch.H"
 #include "OFstream.H"
 #include "mathematicalConstants.H"
+#include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -75,13 +76,14 @@ int main(int argc, char *argv[])
     #include "setInitialDeltaT.H"
     #include "startSummary.H"
 
+    pimpleControl pimple(mesh);
+
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
 
     while (runTime.run())
     {
-        #include "readPISOControls.H"
         #include "readEngineTimeControls.H"
         #include "compressibleCourantNo.H"
         #include "setDeltaT.H"
@@ -94,11 +96,11 @@ int main(int argc, char *argv[])
 
         #include "rhoEqn.H"
 
-        #include "UEqn.H"
+        // --- Pressure-velocity PIMPLE corrector loop
-
+        for (pimple.start(); pimple.loop(); pimple++)
-        // --- PISO loop
-        for (int corr=1; corr<=nCorr; corr++)
         {
+            #include "UEqn.H"
+
             #include "ftEqn.H"
             #include "bEqn.H"
             #include "huEqn.H"
@@ -109,10 +111,17 @@ int main(int argc, char *argv[])
                 hu == h;
             }
 
-            #include "pEqn.H"
+            // --- PISO loop
-        }
+            for (int corr=1; corr<=pimple.nCorr(); corr++)
+            {
+                #include "pEqn.H"
+            }
 
-        turbulence->correct();
+            if (pimple.turbCorr())
+            {
+                turbulence->correct();
+            }
+        }
 
         #include "logSummary.H"
 

applications/solvers/combustion/engineFoam/pEqn.H

@@ -3,7 +3,7 @@ rho = thermo.rho();
 volScalarField rAU(1.0/UEqn.A());
 U = rAU*UEqn.H();
 
-if (transonic)
+if (pimple.transonic())
 {
     surfaceScalarField phid
     (
@@ -12,18 +12,21 @@ if (transonic)
        *((fvc::interpolate(U) & mesh.Sf()) - fvc::meshPhi(rho, U))
     );
 
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
     {
         fvScalarMatrix pEqn
         (
             fvm::ddt(psi, p)
-          + fvm::div(phid, p, "div(phid,p)")
+          + fvm::div(phid, p)
           - fvm::laplacian(rho*rAU, p)
         );
 
-        pEqn.solve();
+        pEqn.solve
+        (
+            mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
+        );
 
-        if (nonOrth == nNonOrthCorr)
+        if (nonOrth == pimple.nNonOrthCorr())
         {
             phi == pEqn.flux();
         }
@@ -34,7 +37,7 @@ else
     phi = fvc::interpolate(rho)
          *((fvc::interpolate(U) & mesh.Sf()) - fvc::meshPhi(rho, U));
 
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
     {
         fvScalarMatrix pEqn
         (
@@ -43,9 +46,12 @@ else
           - fvm::laplacian(rho*rAU, p)
         );
 
-        pEqn.solve();
+        pEqn.solve
+        (
+            mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
+        );
 
-        if (nonOrth == nNonOrthCorr)
+        if (nonOrth == pimple.nNonOrthCorr())
         {
             phi += pEqn.flux();
         }

applications/solvers/combustion/fireFoam/fireFoam.C

@@ -34,7 +34,7 @@ Description
 #include "turbulenceModel.H"
 #include "combustionModel.H"
 #include "radiationModel.H"
-#include "pimpleLoop.H"
+#include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -51,13 +51,14 @@ int main(int argc, char *argv[])
     #include "compressibleCourantNo.H"
     #include "setInitialDeltaT.H"
 
+    pimpleControl pimple(mesh);
+
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
 
     while (runTime.run())
     {
-        #include "readPIMPLEControls.H"
         #include "readTimeControls.H"
         #include "compressibleCourantNo.H"
         #include "setDeltaT.H"
@@ -68,12 +69,7 @@ int main(int argc, char *argv[])
         #include "rhoEqn.H"
 
         // --- Pressure-velocity PIMPLE corrector loop
-        for
+        for (pimple.start(); pimple.loop(); pimple++)
-        (
-            pimpleLoop pimpleCorr(mesh, nOuterCorr);
-            pimpleCorr.loop();
-            pimpleCorr++
-        )
         {
             #include "UEqn.H"
 
@@ -81,13 +77,16 @@ int main(int argc, char *argv[])
             #include "fuhsEqn.H"
 
             // --- PISO loop
-            for (int corr=0; corr<nCorr; corr++)
+            for (int corr=0; corr<pimple.nCorr(); corr++)
             {
                 #include "pEqn.H"
             }
-        }
 
-        turbulence->correct();
+            if (pimple.turbCorr())
+            {
+                turbulence->correct();
+            }
+        }
 
         rho = thermo.rho();
 

applications/solvers/combustion/fireFoam/pEqn.H

@@ -15,7 +15,7 @@ surfaceScalarField phiU
 
 phi = phiU - rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf();
 
-for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
 {
     surfaceScalarField rhorAUf(fvc::interpolate(rho*rAU));
 
@@ -28,20 +28,10 @@ for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
 
     p_rghEqn.solve
     (
-        mesh.solver
+        mesh.solver(p_rgh.select(pimple.finalInnerIter(corr, nonOrth)))
-        (
-            p_rgh.select
-            (
-                (
-                    pimpleCorr.finalIter()
-                 && corr == nCorr-1
-                 && nonOrth == nNonOrthCorr
-                )
-            )
-        )
     );
 
-    if (nonOrth == nNonOrthCorr)
+    if (nonOrth == pimple.nNonOrthCorr())
     {
         phi += p_rghEqn.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/chemistry.H

@@ -1,3 +1,4 @@
+if (chemistry.chemistry())
 {
     Info<< "Solving chemistry" << endl;
 
@@ -10,17 +11,29 @@
     // turbulent time scale
     if (turbulentReaction)
     {
-        volScalarField tk
+        tmp<volScalarField> tepsilon(turbulence->epsilon());
-        (
+        const volScalarField& epsilon = tepsilon();
-            Cmix*sqrt(turbulence->muEff()/rho/turbulence->epsilon())
+        tmp<volScalarField> tmuEff(turbulence->muEff());
-        );
+        const volScalarField& muEff = tmuEff();
-        volScalarField tc
+        tmp<volScalarField> ttc(chemistry.tc());
-        (
+        const volScalarField& tc = ttc();
-            chemistry.tc()
-        );
 
-        // Chalmers PaSR model
+        forAll(epsilon, i)
-        kappa = (runTime.deltaT() + tc)/(runTime.deltaT() + tc + tk);
+        {
+            if (epsilon[i] > 0)
+            {
+                // Chalmers PaSR model
+                scalar tk = Cmix.value()*Foam::sqrt(muEff[i]/rho[i]/epsilon[i]);
+                kappa[i] =
+                    (runTime.deltaTValue() + tc[i])
+                   /(runTime.deltaTValue() + tc[i] + tk);
+            }
+            else
+            {
+                // Return to laminar combustion
+                kappa[i] = 1.0;
+            }
+        }
     }
     else
     {

applications/solvers/combustion/reactingFoam/pEqn.H

@@ -3,7 +3,7 @@ rho = thermo.rho();
 volScalarField rAU(1.0/UEqn.A());
 U = rAU*UEqn.H();
 
-if (transonic)
+if (pimple.transonic())
 {
     surfaceScalarField phid
     (
@@ -15,7 +15,7 @@ if (transonic)
         )
     );
 
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
     {
         fvScalarMatrix pEqn
         (
@@ -24,9 +24,12 @@ if (transonic)
           - fvm::laplacian(rho*rAU, p)
         );
 
-        pEqn.solve();
+        pEqn.solve
+        (
+            mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
+        );
 
-        if (nonOrth == nNonOrthCorr)
+        if (nonOrth == pimple.nNonOrthCorr())
         {
             phi == pEqn.flux();
         }
@@ -41,7 +44,7 @@ else
           + fvc::ddtPhiCorr(rAU, rho, U, phi)
         );
 
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
     {
         fvScalarMatrix pEqn
         (
@@ -50,9 +53,12 @@ else
           - fvm::laplacian(rho*rAU, p)
         );
 
-        pEqn.solve();
+        pEqn.solve
+        (
+            mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
+        );
 
-        if (nonOrth == nNonOrthCorr)
+        if (nonOrth == pimple.nNonOrthCorr())
         {
             phi += pEqn.flux();
         }

applications/solvers/combustion/reactingFoam/reactingFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -35,6 +35,7 @@ Description
 #include "psiChemistryModel.H"
 #include "chemistrySolver.H"
 #include "multivariateScheme.H"
+#include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -51,6 +52,8 @@ int main(int argc, char *argv[])
     #include "compressibleCourantNo.H"
     #include "setInitialDeltaT.H"
 
+    pimpleControl pimple(mesh);
+
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
@@ -58,7 +61,6 @@ int main(int argc, char *argv[])
     while (runTime.run())
     {
         #include "readTimeControls.H"
-        #include "readPISOControls.H"
         #include "compressibleCourantNo.H"
         #include "setDeltaT.H"
 
@@ -68,20 +70,23 @@ int main(int argc, char *argv[])
         #include "chemistry.H"
         #include "rhoEqn.H"
 
-        for (label ocorr=1; ocorr <= nOuterCorr; ocorr++)
+        for (pimple.start(); pimple.loop(); pimple++)
         {
             #include "UEqn.H"
             #include "YEqn.H"
             #include "hsEqn.H"
 
             // --- PISO loop
-            for (int corr=1; corr<=nCorr; corr++)
+            for (int corr=0; corr<pimple.nCorr(); corr++)
             {
                 #include "pEqn.H"
             }
-        }
 
-        turbulence->correct();
+            if (pimple.turbCorr())
+            {
+                turbulence->correct();
+            }
+        }
 
         if (runTime.write())
         {

applications/solvers/combustion/rhoReactingFoam/Make/options

@@ -5,7 +5,9 @@ EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
     -I$(LIB_SRC)/ODE/lnInclude \
-    -I$(LIB_SRC)/finiteVolume/lnInclude
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(FOAM_SOLVERS)/combustion/reactingFoam
+
 
 EXE_LIBS = \
     -lcompressibleTurbulenceModel \

applications/solvers/combustion/rhoReactingFoam/UEqn.H

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

applications/solvers/combustion/rhoReactingFoam/chemistry.H

@@ -1,28 +0,0 @@
-{
-    Info<< "Solving chemistry" << endl;
-
-    chemistry.solve
-    (
-        runTime.value() - runTime.deltaTValue(),
-        runTime.deltaTValue()
-    );
-
-    // turbulent time scale
-    if (turbulentReaction)
-    {
-        volScalarField tk
-        (
-            Cmix*sqrt(turbulence->muEff()/rho/turbulence->epsilon())
-        );
-        volScalarField tc(chemistry.tc());
-
-        // Chalmers PaSR model
-        kappa = (runTime.deltaT() + tc)/(runTime.deltaT() + tc + tk);
-    }
-    else
-    {
-        kappa = 1.0;
-    }
-
-    chemistrySh = kappa*chemistry.Sh()();
-}

applications/solvers/combustion/rhoReactingFoam/pEqn.H

@@ -8,7 +8,7 @@
     volScalarField rAU(1.0/UEqn.A());
     U = rAU*UEqn.H();
 
-    if (transonic)
+    if (pimple.transonic())
     {
         surfaceScalarField phiv
         (
@@ -30,7 +30,7 @@
           + correction(fvm::ddt(psi, p) + fvm::div(phid, p))
         );
 
-        for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+        for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
         {
             fvScalarMatrix pEqn
             (
@@ -40,20 +40,10 @@
 
             pEqn.solve
             (
-                mesh.solver
+                mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
-                (
-                    p.select
-                    (
-                        (
-                            pimpleCorr.finalIter()
-                         && corr == nCorr-1
-                         && nonOrth == nNonOrthCorr
-                        )
-                    )
-                )
             );
 
-            if (nonOrth == nNonOrthCorr)
+            if (nonOrth == pimple.nNonOrthCorr())
             {
                 phi += pEqn.flux();
             }
@@ -74,7 +64,7 @@
           + fvc::div(phi)
         );
 
-        for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+        for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
         {
             fvScalarMatrix pEqn
             (
@@ -84,20 +74,10 @@
 
             pEqn.solve
             (
-                mesh.solver
+                mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
-                (
-                    p.select
-                    (
-                        (
-                            pimpleCorr.finalIter()
-                         && corr == nCorr-1
-                         && nonOrth == nNonOrthCorr
-                        )
-                    )
-                )
             );
 
-            if (nonOrth == nNonOrthCorr)
+            if (nonOrth == pimple.nNonOrthCorr())
             {
                 phi += pEqn.flux();
             }

applications/solvers/combustion/rhoReactingFoam/rhoReactingFoam.C

@@ -36,7 +36,7 @@ Description
 #include "rhoChemistryModel.H"
 #include "chemistrySolver.H"
 #include "multivariateScheme.H"
-#include "pimpleLoop.H"
+#include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -53,6 +53,8 @@ int main(int argc, char *argv[])
     #include "compressibleCourantNo.H"
     #include "setInitialDeltaT.H"
 
+    pimpleControl pimple(mesh);
+
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
@@ -60,7 +62,6 @@ int main(int argc, char *argv[])
     while (runTime.run())
     {
         #include "readTimeControls.H"
-        #include "readPISOControls.H"
         #include "compressibleCourantNo.H"
         #include "setDeltaT.H"
 
@@ -71,25 +72,23 @@ int main(int argc, char *argv[])
         #include "rhoEqn.H"
 
         // --- Pressure-velocity PIMPLE corrector loop
-        for
+        for (pimple.start(); pimple.loop(); pimple++)
-        (
-            pimpleLoop pimpleCorr(mesh, nOuterCorr);
-            pimpleCorr.loop();
-            pimpleCorr++
-        )
         {
             #include "UEqn.H"
             #include "YEqn.H"
             #include "hsEqn.H"
 
             // --- PISO loop
-            for (int corr=1; corr<=nCorr; corr++)
+            for (int corr=1; corr<=pimple.nCorr(); corr++)
             {
                 #include "pEqn.H"
             }
-        }
 
-        turbulence->correct();
+            if (pimple.turbCorr())
+            {
+                turbulence->correct();
+            }
+        }
 
         rho = thermo.rho();
 

applications/solvers/compressible/rhoPimpleFoam/UEqn.H

@@ -11,7 +11,7 @@ UEqn().relax();
 
 volScalarField rAU(1.0/UEqn().A());
 
-if (momentumPredictor)
+if (pimple.momentumPredictor())
 {
     solve(UEqn() == -fvc::grad(p));
 }

applications/solvers/compressible/rhoPimpleFoam/pEqn.H

@@ -5,12 +5,12 @@ rho.relax();
 
 U = rAU*UEqn().H();
 
-if (nCorr <= 1)
+if (pimple.nCorr() <= 1)
 {
     UEqn.clear();
 }
 
-if (transonic)
+if (pimple.transonic())
 {
     surfaceScalarField phid
     (
@@ -22,7 +22,7 @@ if (transonic)
         )
     );
 
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
     {
         fvScalarMatrix pEqn
         (
@@ -33,18 +33,10 @@ if (transonic)
 
         pEqn.solve
         (
-            mesh.solver
+            mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
-            (
-                p.select
-                (
-                    pimpleCorr.finalIter()
-                 && corr == nCorr-1
-                 && nonOrth == nNonOrthCorr
-                )
-            )
         );
 
-        if (nonOrth == nNonOrthCorr)
+        if (nonOrth == pimple.nNonOrthCorr())
         {
             phi == pEqn.flux();
         }
@@ -59,7 +51,7 @@ else
           + fvc::ddtPhiCorr(rAU, rho, U, phi)
         );
 
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
     {
         // Pressure corrector
         fvScalarMatrix pEqn
@@ -71,18 +63,10 @@ else
 
         pEqn.solve
         (
-            mesh.solver
+            mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
-            (
-                p.select
-                (
-                    pimpleCorr.finalIter()
-                 && corr == nCorr-1
-                 && nonOrth == nNonOrthCorr
-                )
-            )
         );
 
-        if (nonOrth == nNonOrthCorr)
+        if (nonOrth == pimple.nNonOrthCorr())
         {
             phi += pEqn.flux();
         }

applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C

@@ -37,7 +37,7 @@ Description
 #include "basicPsiThermo.H"
 #include "turbulenceModel.H"
 #include "bound.H"
-#include "pimpleLoop.H"
+#include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -49,12 +49,13 @@ int main(int argc, char *argv[])
     #include "createFields.H"
     #include "initContinuityErrs.H"
 
+    pimpleControl pimple(mesh);
+
     Info<< "\nStarting time loop\n" << endl;
 
     while (runTime.run())
     {
         #include "readTimeControls.H"
-        #include "readPIMPLEControls.H"
         #include "compressibleCourantNo.H"
         #include "setDeltaT.H"
 
@@ -65,14 +66,9 @@ int main(int argc, char *argv[])
         #include "rhoEqn.H"
 
         // --- Pressure-velocity PIMPLE corrector loop
-        for
+        for (pimple.start(); pimple.loop(); pimple++)
-        (
-            pimpleLoop pimpleCorr(mesh, nOuterCorr);
-            pimpleCorr.loop();
-            pimpleCorr++
-        )
         {
-            if (nOuterCorr != 1)
+            if (pimple.nOuterCorr() != 1)
             {
                 p.storePrevIter();
                 rho.storePrevIter();
@@ -82,12 +78,15 @@ int main(int argc, char *argv[])
             #include "hEqn.H"
 
             // --- PISO loop
-            for (int corr=0; corr<nCorr; corr++)
+            for (int corr=0; corr<pimple.nCorr(); corr++)
             {
                 #include "pEqn.H"
             }
 
-            turbulence->correct();
+            if (pimple.turbCorr())
+            {
+                turbulence->correct();
+            }
         }
 
         runTime.write();

applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFLTSPimpleFoam/rhoPorousMRFLTSPimpleFoam.C

@@ -40,7 +40,7 @@ Description
 #include "MRFZones.H"
 #include "porousZones.H"
 #include "fvcSmooth.H"
-#include "pimpleLoop.H"
+#include "pimpleControl.H"
 #include "bound.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -50,18 +50,21 @@ int main(int argc, char *argv[])
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createMesh.H"
-    #include "readPIMPLEControls.H"
+
+    pimpleControl pimple(mesh);
+
     #include "setInitialrDeltaT.H"
     #include "createFields.H"
     #include "createZones.H"
     #include "initContinuityErrs.H"
 
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
     Info<< "\nStarting time loop\n" << endl;
 
     while (runTime.run())
     {
         #include "readTimeControls.H"
-        #include "readPIMPLEControls.H"
         #include "compressibleCourantNo.H"
         #include "setDeltaT.H"
 
@@ -74,14 +77,9 @@ int main(int argc, char *argv[])
         #include "rhoEqn.H"
 
         // --- Pressure-velocity PIMPLE corrector loop
-        for
+        for (pimple.start(); pimple.loop(); pimple++)
-        (
-            pimpleLoop pimpleCorr(mesh, nOuterCorr);
-            pimpleCorr.loop();
-            pimpleCorr++
-        )
         {
-            if (nOuterCorr != 1)
+            if (pimple.nOuterCorr() != 1)
             {
                 p.storePrevIter();
                 rho.storePrevIter();
@@ -93,7 +91,7 @@ int main(int argc, char *argv[])
             #include "hEqn.H"
 
             // --- PISO loop
-            for (int corr=0; corr<nCorr; corr++)
+            for (int corr=0; corr<pimple.nCorr(); corr++)
             {
                 #include "pEqn.H"
             }

applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFLTSPimpleFoam/setInitialrDeltaT.H

@@ -1,6 +1,6 @@
 scalar maxDeltaT
 (
-    pimple.lookupOrDefault<scalar>("maxDeltaT", GREAT)
+    pimple.dict().lookupOrDefault<scalar>("maxDeltaT", GREAT)
 );
 
 volScalarField rDeltaT

applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFLTSPimpleFoam/setrDeltaT.H

@@ -1,22 +1,24 @@
 {
+    const dictionary& pimpleDict = pimple.dict();
+
     scalar maxCo
     (
-        pimple.lookupOrDefault<scalar>("maxCo", 0.8)
+        pimpleDict.lookupOrDefault<scalar>("maxCo", 0.8)
     );
 
     scalar rDeltaTSmoothingCoeff
     (
-        pimple.lookupOrDefault<scalar>("rDeltaTSmoothingCoeff", 0.02)
+        pimpleDict.lookupOrDefault<scalar>("rDeltaTSmoothingCoeff", 0.02)
     );
 
     scalar rDeltaTDampingCoeff
     (
-        pimple.lookupOrDefault<scalar>("rDeltaTDampingCoeff", 1.0)
+        pimpleDict.lookupOrDefault<scalar>("rDeltaTDampingCoeff", 1.0)
     );
 
     scalar maxDeltaT
     (
-        pimple.lookupOrDefault<scalar>("maxDeltaT", GREAT)
+        pimpleDict.lookupOrDefault<scalar>("maxDeltaT", GREAT)
     );
 
     volScalarField rDeltaT0("rDeltaT0", rDeltaT);
@@ -29,7 +31,7 @@
        /((2*maxCo)*mesh.V()*rho.dimensionedInternalField())
     );
 
-    if (transonic)
+    if (pimple.transonic())
     {
         surfaceScalarField phid
         (

applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/UEqn.H

@@ -15,7 +15,7 @@ pZones.addResistance(UEqn());
 
 volScalarField rAU(1.0/UEqn().A());
 
-if (momentumPredictor)
+if (pimple.momentumPredictor())
 {
     solve(UEqn() == -fvc::grad(p));
 }

applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/pEqn.H

@@ -6,12 +6,12 @@ rho.relax();
 volScalarField rAU(1.0/UEqn().A());
 U = rAU*UEqn().H();
 
-if (nCorr <= 1)
+if (pimple.nCorr() <= 1)
 {
     UEqn.clear();
 }
 
-if (transonic)
+if (pimple.transonic())
 {
     surfaceScalarField phid
     (
@@ -24,7 +24,7 @@ if (transonic)
     );
     mrfZones.relativeFlux(fvc::interpolate(psi), phid);
 
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
     {
         fvScalarMatrix pEqn
         (
@@ -35,20 +35,10 @@ if (transonic)
 
         pEqn.solve
         (
-            mesh.solver
+            mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
-            (
-                p.select
-                (
-                    (
-                       pimpleCorr.finalIter()
-                     && corr == nCorr-1
-                     && nonOrth == nNonOrthCorr
-                    )
-                )
-            )
         );
 
-        if (nonOrth == nNonOrthCorr)
+        if (nonOrth == pimple.nNonOrthCorr())
         {
             phi == pEqn.flux();
         }
@@ -64,7 +54,7 @@ else
         );
     mrfZones.relativeFlux(fvc::interpolate(rho), phi);
 
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
     {
         // Pressure corrector
         fvScalarMatrix pEqn
@@ -76,20 +66,10 @@ else
 
         pEqn.solve
         (
-            mesh.solver
+            mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
-            (
-                p.select
-                (
-                    (
-                        pimpleCorr.finalIter()
-                     && corr == nCorr-1
-                     && nonOrth == nNonOrthCorr
-                    )
-                )
-            )
         );
 
-        if (nonOrth == nNonOrthCorr)
+        if (nonOrth == pimple.nNonOrthCorr())
         {
             phi += pEqn.flux();
         }

applications/solvers/compressible/rhoPimpleFoam/rhoPorousMRFPimpleFoam/rhoPorousMRFPimpleFoam.C

@@ -39,7 +39,7 @@ Description
 #include "bound.H"
 #include "MRFZones.H"
 #include "porousZones.H"
-#include "pimpleLoop.H"
+#include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -52,6 +52,8 @@ int main(int argc, char *argv[])
     #include "createZones.H"
     #include "initContinuityErrs.H"
 
+    pimpleControl pimple(mesh);
+
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
@@ -59,7 +61,6 @@ int main(int argc, char *argv[])
     while (runTime.run())
     {
         #include "readTimeControls.H"
-        #include "readPIMPLEControls.H"
         #include "compressibleCourantNo.H"
         #include "setDeltaT.H"
 
@@ -70,14 +71,9 @@ int main(int argc, char *argv[])
         #include "rhoEqn.H"
 
         // --- Pressure-velocity PIMPLE corrector loop
-        for
+        for (pimple.start(); pimple.loop(); pimple++)
-        (
-            pimpleLoop pimpleCorr(mesh, nOuterCorr);
-            pimpleCorr.loop();
-            pimpleCorr++
-        )
         {
-            if (nOuterCorr != 1)
+            if (pimple.nOuterCorr() != 1)
             {
                 p.storePrevIter();
                 rho.storePrevIter();
@@ -87,12 +83,15 @@ int main(int argc, char *argv[])
             #include "hEqn.H"
 
             // --- PISO loop
-            for (int corr=0; corr<nCorr; corr++)
+            for (int corr=0; corr<pimple.nCorr(); corr++)
             {
                 #include "pEqn.H"
             }
 
-            turbulence->correct();
+            if (pimple.turbCorr())
+            {
+                turbulence->correct();
+            }
         }
 
         runTime.write();

applications/solvers/compressible/rhoSimpleFoam/pEqn.H

@@ -9,7 +9,7 @@ UEqn.clear();
 
 bool closedVolume = false;
 
-if (transonic)
+if (simple.transonic())
 {
     surfaceScalarField phid
     (
@@ -17,7 +17,7 @@ if (transonic)
         fvc::interpolate(psi)*(fvc::interpolate(U) & mesh.Sf())
     );
 
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    for (int nonOrth=0; nonOrth<=simple.nNonOrthCorr(); nonOrth++)
     {
         fvScalarMatrix pEqn
         (
@@ -32,7 +32,7 @@ if (transonic)
 
         pEqn.solve();
 
-        if (nonOrth == nNonOrthCorr)
+        if (nonOrth == simple.nNonOrthCorr())
         {
             phi == pEqn.flux();
         }
@@ -43,7 +43,7 @@ else
     phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf());
     closedVolume = adjustPhi(phi, U, p);
 
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    for (int nonOrth=0; nonOrth<=simple.nNonOrthCorr(); nonOrth++)
     {
         fvScalarMatrix pEqn
         (
@@ -54,7 +54,7 @@ else
 
         pEqn.solve();
 
-        if (nonOrth == nNonOrthCorr)
+        if (nonOrth == simple.nNonOrthCorr())
         {
             phi -= pEqn.flux();
         }

applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/pEqn.H

@@ -11,7 +11,7 @@ UEqn.clear();
 
 bool closedVolume = false;
 
-if (transonic)
+if (simple.transonic())
 {
     surfaceScalarField phid
     (
@@ -20,7 +20,7 @@ if (transonic)
     );
     mrfZones.relativeFlux(fvc::interpolate(psi), phid);
 
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    for (int nonOrth=0; nonOrth<=simple.nNonOrthCorr(); nonOrth++)
     {
         tmp<fvScalarMatrix> tpEqn;
 
@@ -37,7 +37,7 @@ if (transonic)
 
         tpEqn().solve();
 
-        if (nonOrth == nNonOrthCorr)
+        if (nonOrth == simple.nNonOrthCorr())
         {
             phi == tpEqn().flux();
         }
@@ -50,7 +50,7 @@ else
 
     closedVolume = adjustPhi(phi, U, p);
 
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    for (int nonOrth=0; nonOrth<=simple.nNonOrthCorr(); nonOrth++)
     {
         tmp<fvScalarMatrix> tpEqn;
 
@@ -67,7 +67,7 @@ else
 
         tpEqn().solve();
 
-        if (nonOrth == nNonOrthCorr)
+        if (nonOrth == simple.nNonOrthCorr())
         {
             phi -= tpEqn().flux();
         }

applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/rhoPorousMRFSimpleFoam.C

@@ -36,6 +36,7 @@ Description
 #include "RASModel.H"
 #include "MRFZones.H"
 #include "thermalPorousZones.H"
+#include "simpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -48,16 +49,16 @@ int main(int argc, char *argv[])
     #include "createZones.H"
     #include "initContinuityErrs.H"
 
+    simpleControl simple(mesh);
+
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
 
-    while (runTime.loop())
+    while (simple.loop())
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        #include "readSIMPLEControls.H"
-
         p.storePrevIter();
         rho.storePrevIter();
 

applications/solvers/compressible/rhoSimpleFoam/rhoSimpleFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -33,6 +33,7 @@ Description
 #include "fvCFD.H"
 #include "basicPsiThermo.H"
 #include "RASModel.H"
+#include "simpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -44,16 +45,16 @@ int main(int argc, char *argv[])
     #include "createFields.H"
     #include "initContinuityErrs.H"
 
+    simpleControl simple(mesh);
+
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
 
-    while (runTime.loop())
+    while (simple.loop())
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        #include "readSIMPLEControls.H"
-
         p.storePrevIter();
         rho.storePrevIter();
 

applications/solvers/compressible/rhoSimplecFoam/pEqn.H

@@ -12,9 +12,9 @@ UEqn.clear();
 
 bool closedVolume = false;
 
-if (transonic)
+if (simple.transonic())
 {
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    for (int nonOrth=0; nonOrth<=simple.nNonOrthCorr(); nonOrth++)
     {
         surfaceScalarField phid
         (
@@ -54,7 +54,7 @@ if (transonic)
             pEqn.solve();
         }
 
-        if (nonOrth == nNonOrthCorr)
+        if (nonOrth == simple.nNonOrthCorr())
         {
             phi == phic + pEqn.flux();
         }
@@ -62,7 +62,7 @@ if (transonic)
 }
 else
 {
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    for (int nonOrth=0; nonOrth<=simple.nNonOrthCorr(); nonOrth++)
     {
         phi = fvc::interpolate(rho*U) & mesh.Sf();
         closedVolume = adjustPhi(phi, U, p);
@@ -88,7 +88,7 @@ else
         }
 
 
-        if (nonOrth == nNonOrthCorr)
+        if (nonOrth == simple.nNonOrthCorr())
         {
             phi += pEqn.flux();
         }
@@ -119,7 +119,7 @@ rho = thermo.rho();
 rho = max(rho, rhoMin);
 rho = min(rho, rhoMax);
 
-if (!transonic)
+if (!simple.transonic())
 {
     rho.relax();
 }

applications/solvers/compressible/rhoSimplecFoam/rhoSimplecFoam.C

@@ -35,6 +35,7 @@ Description
 #include "RASModel.H"
 #include "mixedFvPatchFields.H"
 #include "bound.H"
+#include "simpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -46,20 +47,20 @@ int main(int argc, char *argv[])
     #include "createFields.H"
     #include "initContinuityErrs.H"
 
+    simpleControl simple(mesh);
+
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
 
-    while (runTime.loop())
+    while (simple.loop())
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        #include "readSIMPLEControls.H"
-
         p.storePrevIter();
         rho.storePrevIter();
 
-        if (!transonic)
+        if (!simple.transonic())
         {
             rho.storePrevIter();
         }

applications/solvers/discreteMethods/molecularDynamics/mdEquilibrationFoam/mdEquilibrationFoam.C

@@ -2,7 +2,7 @@
  =========                   |
  \\      /   F ield          | OpenFOAM: The Open Source CFD Toolbox
   \\    /    O peration      |
-   \\  /     A nd            | Copyright (C) 2008-2010 OpenCFD Ltd.
+   \\  /     A nd            | Copyright (C) 2008-2011 OpenCFD Ltd.
     \\/      M anipulation   |
 -------------------------------------------------------------------------------
 License
@@ -34,10 +34,11 @@ Description
 
 int main(int argc, char *argv[])
 {
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
 
-#   include "setRootCase.H"
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-#   include "createTime.H"
-#   include "createMesh.H"
 
     Info<< "\nReading field U\n" << endl;
     volVectorField U
@@ -57,9 +58,9 @@ int main(int argc, char *argv[])
 
     moleculeCloud molecules(mesh, pot);
 
-#   include "temperatureAndPressureVariables.H"
+    #include "temperatureAndPressureVariables.H"
 
-#   include "readmdEquilibrationDict.H"
+    #include "readmdEquilibrationDict.H"
 
     label nAveragingSteps = 0;
 
@@ -67,18 +68,17 @@ int main(int argc, char *argv[])
 
     while (runTime.loop())
     {
-
         nAveragingSteps++;
 
         Info<< "Time = " << runTime.timeName() << endl;
 
         molecules.evolve();
 
-#       include "meanMomentumEnergyAndNMols.H"
+        #include "meanMomentumEnergyAndNMols.H"
 
-#       include "temperatureAndPressure.H"
+        #include "temperatureAndPressure.H"
 
-#       include "temperatureEquilibration.H"
+        #include "temperatureEquilibration.H"
 
         runTime.write();
 
@@ -96,3 +96,6 @@ int main(int argc, char *argv[])
 
     return 0;
 }
+
+
+// ************************************************************************* //

applications/solvers/discreteMethods/molecularDynamics/mdFoam/mdFoam.C

@@ -2,7 +2,7 @@
  =========                   |
  \\      /   F ield          | OpenFOAM: The Open Source CFD Toolbox
   \\    /    O peration      |
-   \\  /     A nd            | Copyright (C) 2008-2010 OpenCFD Ltd.
+   \\  /     A nd            | Copyright (C) 2008-2011 OpenCFD Ltd.
     \\/      M anipulation   |
 -------------------------------------------------------------------------------
 License
@@ -34,10 +34,11 @@ Description
 
 int main(int argc, char *argv[])
 {
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
 
-#   include "setRootCase.H"
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-#   include "createTime.H"
-#   include "createMesh.H"
 
     Info<< "\nReading field U\n" << endl;
     volVectorField U
@@ -57,7 +58,7 @@ int main(int argc, char *argv[])
 
     moleculeCloud molecules(mesh, pot);
 
-#   include "temperatureAndPressureVariables.H"
+    #include "temperatureAndPressureVariables.H"
 
     label nAveragingSteps = 0;
 
@@ -65,16 +66,15 @@ int main(int argc, char *argv[])
 
     while (runTime.loop())
     {
-
         nAveragingSteps++;
 
         Info<< "Time = " << runTime.timeName() << endl;
 
         molecules.evolve();
 
-#       include "meanMomentumEnergyAndNMols.H"
+        #include "meanMomentumEnergyAndNMols.H"
 
-#       include "temperatureAndPressure.H"
+        #include "temperatureAndPressure.H"
 
         runTime.write();
 
@@ -92,3 +92,6 @@ int main(int argc, char *argv[])
 
     return 0;
 }
+
+
+// ************************************************************************* //

applications/solvers/electromagnetics/magneticFoam/magneticFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2010-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -38,6 +38,7 @@ Description
 #include "OSspecific.H"
 #include "magnet.H"
 #include "electromagneticConstants.H"
+#include "simpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -65,7 +66,8 @@ int main(int argc, char *argv[])
     #include "createTime.H"
     #include "createMesh.H"
     #include "createFields.H"
-    #include "readSIMPLEControls.H"
+
+    simpleControl simple(mesh);
 
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -73,7 +75,7 @@ int main(int argc, char *argv[])
 
     runTime++;
 
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    for (int nonOrth=0; nonOrth<=simple.nNonOrthCorr(); nonOrth++)
     {
         solve(fvm::laplacian(murf, psi) + fvc::div(murf*Mrf));
     }

applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/UEqn.H

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

applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/buoyantBoussinesqPimpleFoam.C

@@ -48,7 +48,7 @@ Description
 #include "fvCFD.H"
 #include "singlePhaseTransportModel.H"
 #include "RASModel.H"
-#include "pimpleLoop.H"
+#include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -64,6 +64,8 @@ int main(int argc, char *argv[])
     #include "CourantNo.H"
     #include "setInitialDeltaT.H"
 
+    pimpleControl pimple(mesh);
+
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
@@ -73,19 +75,13 @@ int main(int argc, char *argv[])
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
         #include "readTimeControls.H"
-        #include "readPIMPLEControls.H"
         #include "CourantNo.H"
         #include "setDeltaT.H"
 
         // --- Pressure-velocity PIMPLE corrector loop
-        for
+        for (pimple.start(); pimple.loop(); pimple++)
-        (
-            pimpleLoop pimpleCorr(mesh, nOuterCorr);
-            pimpleCorr.loop();
-            pimpleCorr++
-        )
         {
-            if (nOuterCorr != 1)
+            if (pimple.nOuterCorr() != 1)
             {
                 p_rgh.storePrevIter();
             }
@@ -94,12 +90,15 @@ int main(int argc, char *argv[])
             #include "TEqn.H"
 
             // --- PISO loop
-            for (int corr=0; corr<nCorr; corr++)
+            for (int corr=0; corr<pimple.nCorr(); corr++)
             {
                 #include "pEqn.H"
             }
 
-            turbulence->correct();
+            if (pimple.turbCorr())
+            {
+                turbulence->correct();
+            }
         }
 
         runTime.write();

applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/pEqn.H

@@ -10,7 +10,7 @@
     surfaceScalarField buoyancyPhi(rAUf*ghf*fvc::snGrad(rhok)*mesh.magSf());
     phi -= buoyancyPhi;
 
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
     {
         fvScalarMatrix p_rghEqn
         (
@@ -21,20 +21,10 @@
 
         p_rghEqn.solve
         (
-            mesh.solver
+            mesh.solver(p_rgh.select(pimple.finalInnerIter(corr, nonOrth)))
-            (
-                p_rgh.select
-                (
-                    (
-                        pimpleCorr.finalIter()
-                     && corr == nCorr-1
-                     && nonOrth == nNonOrthCorr
-                    )
-                )
-            )
         );
 
-        if (nonOrth == nNonOrthCorr)
+        if (nonOrth == pimple.nNonOrthCorr())
         {
             // Calculate the conservative fluxes
             phi -= p_rghEqn.flux();

applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/UEqn.H

@@ -8,7 +8,7 @@
 
     UEqn().relax();
 
-    if (momentumPredictor)
+    if (simple.momentumPredictor())
     {
         solve
         (

applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/buoyantBoussinesqSimpleFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -48,6 +48,7 @@ Description
 #include "fvCFD.H"
 #include "singlePhaseTransportModel.H"
 #include "RASModel.H"
+#include "simpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -60,16 +61,16 @@ int main(int argc, char *argv[])
     #include "createFields.H"
     #include "initContinuityErrs.H"
 
+    simpleControl simple(mesh);
+
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
 
-    while (runTime.loop())
+    while (simple.loop())
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        #include "readSIMPLEControls.H"
-
         p_rgh.storePrevIter();
 
         // Pressure-velocity SIMPLE corrector

applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pEqn.H

@@ -11,7 +11,7 @@
     surfaceScalarField buoyancyPhi(rAUf*ghf*fvc::snGrad(rhok)*mesh.magSf());
     phi -= buoyancyPhi;
 
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    for (int nonOrth=0; nonOrth<=simple.nNonOrthCorr(); nonOrth++)
     {
         fvScalarMatrix p_rghEqn
         (
@@ -22,7 +22,7 @@
 
         p_rghEqn.solve();
 
-        if (nonOrth == nNonOrthCorr)
+        if (nonOrth == simple.nNonOrthCorr())
         {
             // Calculate the conservative fluxes
             phi -= p_rghEqn.flux();

applications/solvers/heatTransfer/buoyantPimpleFoam/UEqn.H

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

applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.C

@@ -37,7 +37,7 @@ Description
 #include "basicRhoThermo.H"
 #include "turbulenceModel.H"
 #include "fixedGradientFvPatchFields.H"
-#include "pimpleLoop.H"
+#include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -53,6 +53,8 @@ int main(int argc, char *argv[])
     #include "compressibleCourantNo.H"
     #include "setInitialDeltaT.H"
 
+    pimpleControl pimple(mesh);
+
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
@@ -60,7 +62,6 @@ int main(int argc, char *argv[])
     while (runTime.run())
     {
         #include "readTimeControls.H"
-        #include "readPIMPLEControls.H"
         #include "compressibleCourantNo.H"
         #include "setDeltaT.H"
 
@@ -71,14 +72,9 @@ int main(int argc, char *argv[])
         #include "rhoEqn.H"
 
         // --- Pressure-velocity PIMPLE corrector loop
-        for
+        for (pimple.start(); pimple.loop(); pimple++)
-        (
-            pimpleLoop pimpleCorr(mesh, nOuterCorr);
-            pimpleCorr.loop();
-            pimpleCorr++
-        )
         {
-            if (nOuterCorr != 1)
+            if (pimple.nOuterCorr() != 1)
             {
                 p_rgh.storePrevIter();
             }
@@ -87,16 +83,19 @@ int main(int argc, char *argv[])
             #include "hEqn.H"
 
             // --- PISO loop
-            for (int corr=0; corr<nCorr; corr++)
+            for (int corr=0; corr<pimple.nCorr(); corr++)
             {
                 #include "pEqn.H"
             }
 
-            turbulence->correct();
+            if (pimple.turbCorr())
-
+            {
-            rho = thermo.rho();
+                turbulence->correct();
+            }
         }
 
+        rho = thermo.rho();
+
         runTime.write();
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"

applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H

@@ -25,7 +25,7 @@
       + fvc::div(phi)
     );
 
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
     {
         fvScalarMatrix p_rghEqn
         (
@@ -35,20 +35,10 @@
 
         p_rghEqn.solve
         (
-            mesh.solver
+            mesh.solver(p_rgh.select(pimple.finalInnerIter(corr, nonOrth)))
-            (
-                p_rgh.select
-                (
-                    (
-                        pimpleCorr.finalIter()
-                     && corr == nCorr-1
-                     && nonOrth == nNonOrthCorr
-                    )
-                )
-            )
         );
 
-        if (nonOrth == nNonOrthCorr)
+        if (nonOrth == pimple.nNonOrthCorr())
         {
             // Calculate the conservative fluxes
             phi += p_rghEqn.flux();

applications/solvers/heatTransfer/buoyantSimpleFoam/UEqn.H

@@ -8,7 +8,7 @@
 
     UEqn().relax();
 
-    if (momentumPredictor)
+    if (simple.momentumPredictor())
     {
         solve
         (

applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -33,6 +33,7 @@ Description
 #include "basicPsiThermo.H"
 #include "RASModel.H"
 #include "fixedGradientFvPatchFields.H"
+#include "simpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -45,16 +46,16 @@ int main(int argc, char *argv[])
     #include "createFields.H"
     #include "initContinuityErrs.H"
 
+    simpleControl simple(mesh);
+
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
 
-    while (runTime.loop())
+    while (simple.loop())
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        #include "readSIMPLEControls.H"
-
         p_rgh.storePrevIter();
         rho.storePrevIter();
 

applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H

@@ -14,7 +14,7 @@
     surfaceScalarField buoyancyPhi(rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
     phi -= buoyancyPhi;
 
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    for (int nonOrth=0; nonOrth<=simple.nNonOrthCorr(); nonOrth++)
     {
         fvScalarMatrix p_rghEqn
         (
@@ -24,7 +24,7 @@
         p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
         p_rghEqn.solve();
 
-        if (nonOrth == nNonOrthCorr)
+        if (nonOrth == simple.nNonOrthCorr())
         {
             // Calculate the conservative fluxes
             phi -= p_rghEqn.flux();

applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/buoyantSimpleRadiationFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -35,6 +35,7 @@ Description
 #include "RASModel.H"
 #include "fixedGradientFvPatchFields.H"
 #include "radiationModel.H"
+#include "simpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -48,16 +49,16 @@ int main(int argc, char *argv[])
     #include "createRadiationModel.H"
     #include "initContinuityErrs.H"
 
+    simpleControl simple(mesh);
+
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
 
-    while (runTime.loop())
+    while (simple.loop())
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        #include "readSIMPLEControls.H"
-
         p_rgh.storePrevIter();
         rho.storePrevIter();
 

applications/solvers/heatTransfer/chtMultiRegionFoam/coupledDerivedFvPatchFields/turbulentTemperatureRadCoupledMixed/turbulentTemperatureRadCoupledMixedFvPatchScalarField.C

@@ -8,10 +8,10 @@
 License
     This file is part of OpenFOAM.
 
-    OpenFOAM is free software; you can redistribute it and/or modify it
+    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by the
+    under the terms of the GNU General Public License as published by
-    Free Software Foundation; either version 2 of the License, or (at your
+    the Free Software Foundation, either version 3 of the License, or
-    option) any later version.
+    (at your option) any later version.
 
     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
@@ -19,8 +19,7 @@ License
     for more details.
 
     You should have received a copy of the GNU General Public License
-    along with OpenFOAM; if not, write to the Free Software Foundation,
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 \*---------------------------------------------------------------------------*/
 

applications/solvers/incompressible/adjointShapeOptimizationFoam/adjointShapeOptimizationFoam.C

@@ -48,6 +48,7 @@ Description
 #include "fvCFD.H"
 #include "singlePhaseTransportModel.H"
 #include "RASModel.H"
+#include "simpleControl.H"
 
 template<class Type>
 void zeroCells
@@ -74,16 +75,16 @@ int main(int argc, char *argv[])
     #include "createFields.H"
     #include "initContinuityErrs.H"
 
+    simpleControl simple(mesh);
+
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
 
-    while (runTime.loop())
+    while (simple.loop())
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        #include "readSIMPLEControls.H"
-
         p.storePrevIter();
 
         laminarTransport.lookup("lambda") >> lambda;
@@ -121,7 +122,7 @@ int main(int argc, char *argv[])
             adjustPhi(phi, U, p);
 
             // Non-orthogonal pressure corrector loop
-            for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+            for (int nonOrth=0; nonOrth<=simple.nNonOrthCorr(); nonOrth++)
             {
                 fvScalarMatrix pEqn
                 (
@@ -131,7 +132,7 @@ int main(int argc, char *argv[])
                 pEqn.setReference(pRefCell, pRefValue);
                 pEqn.solve();
 
-                if (nonOrth == nNonOrthCorr)
+                if (nonOrth == simple.nNonOrthCorr())
                 {
                     phi -= pEqn.flux();
                 }
@@ -184,7 +185,7 @@ int main(int argc, char *argv[])
             adjustPhi(phia, Ua, pa);
 
             // Non-orthogonal pressure corrector loop
-            for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+            for (int nonOrth=0; nonOrth<=simple.nNonOrthCorr(); nonOrth++)
             {
                 fvScalarMatrix paEqn
                 (
@@ -194,7 +195,7 @@ int main(int argc, char *argv[])
                 paEqn.setReference(paRefCell, paRefValue);
                 paEqn.solve();
 
-                if (nonOrth == nNonOrthCorr)
+                if (nonOrth == simple.nNonOrthCorr())
                 {
                     phia -= paEqn.flux();
                 }

applications/solvers/incompressible/pimpleFoam/UEqn.H

@@ -11,7 +11,7 @@ UEqn().relax();
 
 volScalarField rAU(1.0/UEqn().A());
 
-if (momentumPredictor)
+if (pimple.momentumPredictor())
 {
     solve(UEqn() == -fvc::grad(p));
 }

applications/solvers/incompressible/pimpleFoam/pEqn.H

@@ -1,6 +1,6 @@
 U = rAU*UEqn().H();
 
-if (nCorr <= 1)
+if (pimple.nCorr() <= 1)
 {
     UEqn.clear();
 }
@@ -11,7 +11,7 @@ phi = (fvc::interpolate(U) & mesh.Sf())
 adjustPhi(phi, U, p);
 
 // Non-orthogonal pressure corrector loop
-for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
 {
     // Pressure corrector
     fvScalarMatrix pEqn
@@ -23,18 +23,10 @@ for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
 
     pEqn.solve
     (
-        mesh.solver
+        mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
-        (
-            p.select
-            (
-                pimpleCorr.finalIter()
-             && corr == nCorr-1
-             && nonOrth == nNonOrthCorr
-            )
-        )
     );
 
-    if (nonOrth == nNonOrthCorr)
+    if (nonOrth == pimple.nNonOrthCorr())
     {
         phi -= pEqn.flux();
     }

applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/UEqn.H

@@ -11,7 +11,7 @@ UEqn().relax();
 
 rAU = 1.0/UEqn().A();
 
-if (momentumPredictor)
+if (pimple.momentumPredictor())
 {
     solve(UEqn() == -fvc::grad(p));
 }

applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/correctPhi.H

@@ -50,7 +50,7 @@
         pcorrTypes
     );
 
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
     {
         fvScalarMatrix pcorrEqn
         (
@@ -60,7 +60,7 @@
         pcorrEqn.setReference(pRefCell, pRefValue);
         pcorrEqn.solve();
 
-        if (nonOrth == nNonOrthCorr)
+        if (nonOrth == pimple.nNonOrthCorr())
         {
             phi -= pcorrEqn.flux();
         }

applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pEqn.H

@@ -1,6 +1,6 @@
 U = rAU*UEqn().H();
 
-if (nCorr <= 1)
+if (pimple.nCorr() <= 1)
 {
     UEqn.clear();
 }
@@ -14,7 +14,7 @@ if (p.needReference())
     fvc::makeAbsolute(phi, U);
 }
 
-for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
 {
     fvScalarMatrix pEqn
     (
@@ -25,20 +25,10 @@ for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
 
     pEqn.solve
     (
-        mesh.solver
+        mesh.solver(p.select(pimple.finalInnerIter(corr, nonOrth)))
-        (
-            p.select
-            (
-                (
-                    pimpleCorr.finalIter()
-                 && corr == nCorr-1
-                 && nonOrth == nNonOrthCorr
-                )
-            )
-        )
     );
 
-    if (nonOrth == nNonOrthCorr)
+    if (nonOrth == pimple.nNonOrthCorr())
     {
         phi -= pEqn.flux();
     }

applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pimpleDyMFoam.C

@@ -36,7 +36,7 @@ Description
 #include "singlePhaseTransportModel.H"
 #include "turbulenceModel.H"
 #include "dynamicFvMesh.H"
-#include "pimpleLoop.H"
+#include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -46,11 +46,12 @@ int main(int argc, char *argv[])
 
     #include "createTime.H"
     #include "createDynamicFvMesh.H"
-    #include "readPIMPLEControls.H"
     #include "initContinuityErrs.H"
     #include "createFields.H"
     #include "readTimeControls.H"
 
+    pimpleControl pimple(mesh);
+
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
@@ -85,14 +86,9 @@ int main(int argc, char *argv[])
         }
 
         // --- Pressure-velocity PIMPLE corrector loop
-        for
+        for (pimple.start(); pimple.loop(); pimple++)
-        (
-            pimpleLoop pimpleCorr(mesh, nOuterCorr);
-            pimpleCorr.loop();
-            pimpleCorr++
-        )
         {
-            if (nOuterCorr != 1)
+            if (pimple.nOuterCorr() != 1)
             {
                 p.storePrevIter();
             }
@@ -100,12 +96,15 @@ int main(int argc, char *argv[])
             #include "UEqn.H"
 
             // --- PISO loop
-            for (int corr=0; corr<nCorr; corr++)
+            for (int corr=0; corr<pimple.nCorr(); corr++)
             {
                 #include "pEqn.H"
             }
 
-            turbulence->correct();
+            if (pimple.turbCorr())
+            {
+                turbulence->correct();
+            }
         }
 
         runTime.write();

applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/readControls.H

@@ -1,9 +1,10 @@
     #include "readTimeControls.H"
-    #include "readPIMPLEControls.H"
+
+    const dictionary& pimpleDict = pimple.dict();
 
     const bool correctPhi =
-        pimple.lookupOrDefault("correctPhi", false);
+        pimpleDict.lookupOrDefault("correctPhi", false);
 
     const bool checkMeshCourantNo =
-        pimple.lookupOrDefault("checkMeshCourantNo", false);
+        pimpleDict.lookupOrDefault("checkMeshCourantNo", false);
 

applications/solvers/incompressible/pimpleFoam/pimpleFoam.C

@@ -35,7 +35,7 @@ Description
 #include "fvCFD.H"
 #include "singlePhaseTransportModel.H"
 #include "turbulenceModel.H"
-#include "pimpleLoop.H"
+#include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -47,12 +47,15 @@ int main(int argc, char *argv[])
     #include "createFields.H"
     #include "initContinuityErrs.H"
 
+    pimpleControl pimple(mesh);
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
     Info<< "\nStarting time loop\n" << endl;
 
     while (runTime.run())
     {
         #include "readTimeControls.H"
-        #include "readPIMPLEControls.H"
         #include "CourantNo.H"
         #include "setDeltaT.H"
 
@@ -61,14 +64,9 @@ int main(int argc, char *argv[])
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
         // --- Pressure-velocity PIMPLE corrector loop
-        for
+        for (pimple.start(); pimple.loop(); pimple++)
-        (
-            pimpleLoop pimpleCorr(mesh, nOuterCorr);
-            pimpleCorr.loop();
-            pimpleCorr++
-        )
         {
-            if (nOuterCorr != 1)
+            if (pimple.nOuterCorr() != 1)
             {
                 p.storePrevIter();
             }
@@ -76,12 +74,15 @@ int main(int argc, char *argv[])
             #include "UEqn.H"
 
             // --- PISO loop
-            for (int corr=0; corr<nCorr; corr++)
+            for (int corr=0; corr<pimple.nCorr(); corr++)
             {
                 #include "pEqn.H"
             }
 
-            turbulence->correct();
+            if (pimple.turbCorr())
+            {
+                turbulence->correct();
+            }
         }
 
         runTime.write();

applications/solvers/incompressible/shallowWaterFoam/shallowWaterFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -33,6 +33,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
+#include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -44,6 +45,8 @@ int main(int argc, char *argv[])
     #include "readGravitationalAcceleration.H"
     #include "createFields.H"
 
+    pimpleControl pimple(mesh);
+
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
@@ -52,10 +55,10 @@ int main(int argc, char *argv[])
     {
         Info<< "\n Time = " << runTime.timeName() << nl << endl;
 
-        #include "readPISOControls.H"
         #include "CourantNo.H"
 
-        for (int ucorr=0; ucorr<nOuterCorr; ucorr++)
+        // --- Pressure-velocity PIMPLE corrector loop
+        for (pimple.start(); pimple.loop(); pimple++)
         {
             surfaceScalarField phiv("phiv", phi/fvc::interpolate(h));
 
@@ -67,7 +70,7 @@ int main(int argc, char *argv[])
 
             hUEqn.relax();
 
-            if (momentumPredictor)
+            if (pimple.momentumPredictor())
             {
                 if (rotating)
                 {
@@ -87,7 +90,7 @@ int main(int argc, char *argv[])
             }
 
             // --- PISO loop
-            for (int corr=0; corr<nCorr; corr++)
+            for (int corr=0; corr<pimple.nCorr(); corr++)
             {
                 volScalarField rAU(1.0/hUEqn.A());
                 surfaceScalarField ghrAUf(magg*fvc::interpolate(h*rAU));
@@ -107,7 +110,7 @@ int main(int argc, char *argv[])
                     + fvc::ddtPhiCorr(rAU, h, hU, phi)
                     - phih0;
 
-                for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+                for (int nonOrth=0; nonOrth<=pimple.nNonOrthCorr(); nonOrth++)
                 {
                     fvScalarMatrix hEqn
                     (
@@ -116,16 +119,15 @@ int main(int argc, char *argv[])
                       - fvm::laplacian(ghrAUf, h)
                     );
 
-                    if (ucorr < nOuterCorr-1 || corr < nCorr-1)
+                    hEqn.solve
-                    {
+                    (
-                        hEqn.solve();
+                        mesh.solver
-                    }
+                        (
-                    else
+                            h.select(pimple.finalInnerIter(corr, nonOrth))
-                    {
+                        )
-                        hEqn.solve(mesh.solver(h.name() + "Final"));
+                    );
-                    }
 
-                    if (nonOrth == nNonOrthCorr)
+                    if (nonOrth == pimple.nNonOrthCorr())
                     {
                         phi += hEqn.flux();
                     }

applications/solvers/incompressible/simpleFoam/Allwmake

@@ -3,6 +3,10 @@ cd ${0%/*} || exit 1    # run from this directory
 set -x
 
 wmake
+
+wmake SRFSimpleFoam
+wmake MRFSimpleFoam
+wmake windSimpleFoam
 wmake porousSimpleFoam
 
 # ----------------------------------------------------------------- end-of-file

applications/solvers/incompressible/simpleFoam/MRFSimpleFoam/MRFSimpleFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -22,19 +22,19 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Application
-    rhoPisoTwinParcelFoam
+    MRFSimpleFoam
 
 Description
-    Transient solver for compressible, turbulent flow with two thermo-clouds.
+    Steady-state solver for incompressible, turbulent flow of non-Newtonian
+    fluids with MRF regions.
 
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "basicPsiThermo.H"
+#include "singlePhaseTransportModel.H"
-#include "turbulenceModel.H"
+#include "RASModel.H"
-
+#include "MRFZones.H"
-#include "basicThermoCloud.H"
+#include "simpleControl.H"
-#include "basicKinematicCloud.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -44,54 +44,32 @@ int main(int argc, char *argv[])
 
     #include "createTime.H"
     #include "createMesh.H"
-    #include "readGravitationalAcceleration.H"
     #include "createFields.H"
-    #include "createClouds.H"
-    #include "readPISOControls.H"
     #include "initContinuityErrs.H"
-    #include "readTimeControls.H"
+
-    #include "compressibleCourantNo.H"
+    MRFZones mrfZones(mesh);
-    #include "setInitialDeltaT.H"
+    mrfZones.correctBoundaryVelocity(U);
+
+    simpleControl simple(mesh);
 
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
 
-    while (runTime.run())
+    while (simple.loop())
     {
-        #include "readTimeControls.H"
-        #include "readPISOControls.H"
-        #include "compressibleCourantNo.H"
-        #include "setDeltaT.H"
-
-        runTime++;
-
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        thermoCloud1.evolve();
+        p.storePrevIter();
 
-        kinematicCloud1.evolve();
+        // --- Pressure-velocity SIMPLE corrector
-
-
-        #include "rhoEqn.H"
-
-        // --- PIMPLE loop
-        for (int ocorr=1; ocorr<=nOuterCorr; ocorr++)
         {
             #include "UEqn.H"
-
+            #include "pEqn.H"
-            // --- PISO loop
-            for (int corr=1; corr<=nCorr; corr++)
-            {
-                #include "hsEqn.H"
-                #include "pEqn.H"
-            }
         }
 
         turbulence->correct();
 
-        rho = thermo.rho();
-
         runTime.write();
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"

applications/solvers/incompressible/simpleFoam/MRFSimpleFoam/Make/files

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

applications/solvers/incompressible/simpleFoam/MRFSimpleFoam/Make/options

@@ -1,4 +1,5 @@
 EXE_INC = \
+    -I.. \
     -I$(LIB_SRC)/turbulenceModels \
     -I$(LIB_SRC)/turbulenceModels/incompressible/RAS/RASModel \
     -I$(LIB_SRC)/transportModels \

applications/solvers/incompressible/simpleFoam/MRFSimpleFoam/UEqn.H

@@ -0,0 +1,12 @@
+    // Momentum predictor
+
+    tmp<fvVectorMatrix> UEqn
+    (
+        fvm::div(phi, U)
+        + turbulence->divDevReff(U)
+    );
+    mrfZones.addCoriolis(UEqn());
+
+    UEqn().relax();
+
+    solve(UEqn() == -fvc::grad(p));

applications/solvers/incompressible/simpleFoam/MRFSimpleFoam/pEqn.H

@@ -0,0 +1,37 @@
+{
+    p.boundaryField().updateCoeffs();
+
+    volScalarField rAU(1.0/UEqn().A());
+    U = rAU*UEqn().H();
+    UEqn.clear();
+
+    phi = fvc::interpolate(U, "interpolate(HbyA)") & mesh.Sf();
+    mrfZones.relativeFlux(phi);
+    adjustPhi(phi, U, p);
+
+    // Non-orthogonal pressure corrector loop
+    for (int nonOrth=0; nonOrth<=simple.nNonOrthCorr(); nonOrth++)
+    {
+        fvScalarMatrix pEqn
+        (
+            fvm::laplacian(rAU, p) == fvc::div(phi)
+        );
+
+        pEqn.setReference(pRefCell, pRefValue);
+        pEqn.solve();
+
+        if (nonOrth == simple.nNonOrthCorr())
+        {
+            phi -= pEqn.flux();
+        }
+    }
+
+    #include "continuityErrs.H"
+
+    // Explicitly relax pressure for momentum corrector
+    p.relax();
+
+    // Momentum corrector
+    U -= rAU*fvc::grad(p);
+    U.correctBoundaryConditions();
+}

applications/solvers/incompressible/simpleFoam/SRFSimpleFoam/Make/files

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

applications/solvers/incompressible/simpleFoam/SRFSimpleFoam/Make/options

@@ -1,4 +1,5 @@
 EXE_INC = \
+    -I.. \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/turbulenceModels \
     -I$(LIB_SRC)/turbulenceModels/incompressible/RAS/RASModel \

applications/solvers/incompressible/simpleFoam/SRFSimpleFoam/SRFSimpleFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -22,7 +22,7 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Application
-    simpleSRFFoam
+    SRFSimpleFoam
 
 Description
     Steady-state solver for incompressible, turbulent flow of non-Newtonian
@@ -34,6 +34,7 @@ Description
 #include "singlePhaseTransportModel.H"
 #include "RASModel.H"
 #include "SRFModel.H"
+#include "simpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -45,20 +46,19 @@ int main(int argc, char *argv[])
     #include "createFields.H"
     #include "initContinuityErrs.H"
 
+    simpleControl simple(mesh);
+
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
 
-    while (runTime.loop())
+    while (simple.loop())
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        #include "readSIMPLEControls.H"
-        #include "initConvergenceCheck.H"
-
         p.storePrevIter();
 
-        // Pressure-velocity SIMPLE corrector
+        // --- Pressure-velocity SIMPLE corrector
         {
             #include "UrelEqn.H"
             #include "pEqn.H"
@@ -87,8 +87,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/simpleFoam/SRFSimpleFoam/UrelEqn.H

@@ -1,4 +1,5 @@
     // Relative momentum predictor
+
     tmp<fvVectorMatrix> UrelEqn
     (
         fvm::div(phi, Urel)
@@ -8,9 +9,4 @@
 
     UrelEqn().relax();
 
-    eqnResidual = solve
+    solve(UrelEqn() == -fvc::grad(p));
-    (
-        UrelEqn() == -fvc::grad(p)
-    ).initialResidual();
-
-    maxResidual = max(eqnResidual, maxResidual);

applications/solvers/incompressible/simpleFoam/SRFSimpleFoam/pEqn.H

@@ -1,34 +1,37 @@
 {
     p.boundaryField().updateCoeffs();
-    volScalarField AUrel = UrelEqn().A();
+
-    Urel = UrelEqn().H()/AUrel;
+    volScalarField rAUrel(1.0/UrelEqn().A());
+    Urel = rAUrel*UrelEqn().H();
     UrelEqn.clear();
-    phi = fvc::interpolate(Urel) & mesh.Sf();
+
+    phi = fvc::interpolate(Urel, "interpolate(HbyA)") & mesh.Sf();
     adjustPhi(phi, Urel, p);
 
     // Non-orthogonal pressure corrector loop
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    for (int nonOrth=0; nonOrth<=simple.nNonOrthCorr(); nonOrth++)
     {
         fvScalarMatrix pEqn
         (
-            fvm::laplacian(1.0/AUrel, p) == fvc::div(phi)
+            fvm::laplacian(rAUrel, p) == fvc::div(phi)
         );
 
         pEqn.setReference(pRefCell, pRefValue);
+
         pEqn.solve();
 
-        if (nonOrth == nNonOrthCorr)
+        if (nonOrth == simple.nNonOrthCorr())
         {
             phi -= pEqn.flux();
-       }
+        }
     }
 
-#   include "continuityErrs.H"
+    #include "continuityErrs.H"
 
     // Explicitly relax pressure for momentum corrector
     p.relax();
 
     // Momentum corrector
-    Urel -= fvc::grad(p)/AUrel;
+    Urel -= rAUrel*fvc::grad(p);
     Urel.correctBoundaryConditions();
 }

applications/solvers/incompressible/simpleFoam/UEqn.H

@@ -1,4 +1,4 @@
-    // Solve the Momentum equation
+    // Momentum predictor
 
     tmp<fvVectorMatrix> UEqn
     (

applications/solvers/incompressible/simpleFoam/createFields.H

@@ -26,14 +26,13 @@
         mesh
     );
 
-#   include "createPhi.H"
+    #include "createPhi.H"
 
 
     label pRefCell = 0;
     scalar pRefValue = 0.0;
     setRefCell(p, mesh.solutionDict().subDict("SIMPLE"), pRefCell, pRefValue);
 
-
     singlePhaseTransportModel laminarTransport(U, phi);
 
     autoPtr<incompressible::RASModel> turbulence

applications/solvers/incompressible/simpleFoam/pEqn.H

@@ -1,24 +1,26 @@
+{
     p.boundaryField().updateCoeffs();
 
-    volScalarField AU(UEqn().A());
+    volScalarField rAU(1.0/UEqn().A());
-    U = UEqn().H()/AU;
+    U = rAU*UEqn().H();
     UEqn.clear();
-    phi = fvc::interpolate(U) & mesh.Sf();
+
+    phi = fvc::interpolate(U, "interpolate(HbyA)") & mesh.Sf();
     adjustPhi(phi, U, p);
 
     // Non-orthogonal pressure corrector loop
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    for (int nonOrth=0; nonOrth<=simple.nNonOrthCorr(); nonOrth++)
     {
         fvScalarMatrix pEqn
         (
-            fvm::laplacian(1.0/AU, p) == fvc::div(phi)
+            fvm::laplacian(rAU, p) == fvc::div(phi)
         );
 
         pEqn.setReference(pRefCell, pRefValue);
 
         pEqn.solve();
 
-        if (nonOrth == nNonOrthCorr)
+        if (nonOrth == simple.nNonOrthCorr())
         {
             phi -= pEqn.flux();
         }
@@ -30,5 +32,6 @@
     p.relax();
 
     // Momentum corrector
-    U -= fvc::grad(p)/AU;
+    U -= rAU*fvc::grad(p);
     U.correctBoundaryConditions();
+}

applications/solvers/incompressible/simpleFoam/porousSimpleFoam/createPorousZones.H

@@ -6,13 +6,7 @@
     if (pZones.size())
     {
         // nUCorrectors for pressureImplicitPorosity
-        if (mesh.solutionDict().subDict("SIMPLE").found("nUCorrectors"))
+        nUCorr = simple.dict().lookupOrDefault<int>("nUCorrectors", 0);
-        {
-            nUCorr = readInt
-            (
-                mesh.solutionDict().subDict("SIMPLE").lookup("nUCorrectors")
-            );
-        }
 
         if (nUCorr > 0)
         {

applications/solvers/incompressible/simpleFoam/porousSimpleFoam/pEqn.H

@@ -11,7 +11,7 @@ UEqn.clear();
 phi = fvc::interpolate(U) & mesh.Sf();
 adjustPhi(phi, U, p);
 
-for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+for (int nonOrth=0; nonOrth<=simple.nNonOrthCorr(); nonOrth++)
 {
     tmp<fvScalarMatrix> tpEqn;
 
@@ -35,7 +35,7 @@ for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
         tpEqn().solve();
     }
 
-    if (nonOrth == nNonOrthCorr)
+    if (nonOrth == simple.nNonOrthCorr())
     {
         phi -= tpEqn().flux();
     }

applications/solvers/incompressible/simpleFoam/porousSimpleFoam/porousSimpleFoam.C

@@ -34,6 +34,7 @@ Description
 #include "singlePhaseTransportModel.H"
 #include "RASModel.H"
 #include "porousZones.H"
+#include "simpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -42,6 +43,9 @@ int main(int argc, char *argv[])
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createMesh.H"
+
+    simpleControl simple(mesh);
+
     #include "createFields.H"
     #include "createPorousZones.H"
     #include "initContinuityErrs.H"
@@ -50,12 +54,10 @@ int main(int argc, char *argv[])
 
     Info<< "\nStarting time loop\n" << endl;
 
-    while (runTime.loop())
+    while (simple.loop())
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        #include "readSIMPLEControls.H"
-
         p.storePrevIter();
 
         // Pressure-velocity SIMPLE corrector

applications/solvers/incompressible/simpleFoam/simpleFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -32,6 +32,7 @@ Description
 #include "fvCFD.H"
 #include "singlePhaseTransportModel.H"
 #include "RASModel.H"
+#include "simpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -43,19 +44,19 @@ int main(int argc, char *argv[])
     #include "createFields.H"
     #include "initContinuityErrs.H"
 
+    simpleControl simple(mesh);
+
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
 
-    while (runTime.loop())
+    while (simple.loop())
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        #include "readSIMPLEControls.H"
-
         p.storePrevIter();
 
-        // Pressure-velocity SIMPLE corrector
+        // --- Pressure-velocity SIMPLE corrector
         {
             #include "UEqn.H"
             #include "pEqn.H"

applications/solvers/incompressible/simpleFoam/windSimpleFoam/Make/files

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

applications/solvers/incompressible/simpleFoam/windSimpleFoam/Make/options

@@ -1,4 +1,5 @@
 EXE_INC = \
+    -I.. \
     -I$(LIB_SRC)/turbulenceModels \
     -I$(LIB_SRC)/turbulenceModels/incompressible/RAS/RASModel \
     -I$(LIB_SRC)/transportModels \

applications/solvers/incompressible/simpleFoam/windSimpleFoam/UEqn.H

@@ -11,9 +11,4 @@
     // Add resistance on the actuation disks
     actuationDisks.addSu(UEqn());
 
-    eqnResidual = solve
+    solve(UEqn() == -fvc::grad(p));
-    (
-        UEqn() == -fvc::grad(p)
-    ).initialResidual();
-
-    maxResidual = max(eqnResidual, maxResidual);

applications/solvers/incompressible/simpleFoam/windSimpleFoam/windSimpleFoam.C

@@ -34,6 +34,7 @@ Description
 #include "singlePhaseTransportModel.H"
 #include "RASModel.H"
 #include "IObasicSourceList.H"
+#include "simpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -45,17 +46,18 @@ int main(int argc, char *argv[])
     #include "createFields.H"
     #include "initContinuityErrs.H"
 
+    IObasicSourceList actuationDisks(mesh);
+
+    simpleControl simple(mesh);
+
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
 
-    while (runTime.loop())
+    while (simple.loop())
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        #include "readSIMPLEControls.H"
-        #include "initConvergenceCheck.H"
-
         p.storePrevIter();
 
         // Pressure-velocity SIMPLE corrector
@@ -71,9 +73,6 @@ int main(int argc, char *argv[])
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
             << "  ClockTime = " << runTime.elapsedClockTime() << " s"
             << nl << endl;
-
-         #include "convergenceCheck.H"
-
     }
 
     Info<< "End\n" << endl;

applications/solvers/lagrangian/LTSReactingParcelFoam/LTSReactingParcelFoam.C

@@ -46,7 +46,7 @@ Description
 #include "timeActivatedExplicitSource.H"
 #include "SLGThermo.H"
 #include "fvcSmooth.H"
-#include "pimpleLoop.H"
+#include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -57,6 +57,9 @@ int main(int argc, char *argv[])
     #include "createTime.H"
     #include "createMesh.H"
     #include "readGravitationalAcceleration.H"
+
+    pimpleControl pimple(mesh);
+
     #include "readTimeControls.H"
     #include "readAdditionalSolutionControls.H"
     #include "createFields.H"
@@ -72,7 +75,6 @@ int main(int argc, char *argv[])
 
     while (runTime.run())
     {
-        #include "readPIMPLEControls.H"
         #include "readChemistryProperties.H"
         #include "readAdditionalSolutionControls.H"
         #include "readTimeControls.H"
@@ -89,14 +91,9 @@ int main(int argc, char *argv[])
         #include "rhoEqn.H"
 
         // --- Pressure-velocity PIMPLE corrector loop
-        for
+        for (pimple.start(); pimple.loop(); pimple++)
-        (
-            pimpleLoop pimpleCorr(mesh, nOuterCorr);
-            pimpleCorr.loop();
-            pimpleCorr++
-        )
         {
-            if (nOuterCorr != 1)
+            if (pimple.nOuterCorr() != 1)
             {
                 p.storePrevIter();
             }
@@ -108,7 +105,7 @@ int main(int argc, char *argv[])
             #include "hsEqn.H"
 
             // --- PISO loop
-            for (int corr=0; corr<nCorr; corr++)
+            for (int corr=0; corr<pimple.nCorr(); corr++)
             {
                 #include "pEqn.H"
             }

applications/solvers/lagrangian/LTSReactingParcelFoam/Make/options

@@ -19,7 +19,8 @@ EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/radiationModels/lnInclude \
     -I$(LIB_SRC)/ODE/lnInclude \
     -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
-    -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude
+    -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude
 
 EXE_LIBS = \
     -lfiniteVolume \
@@ -42,4 +43,5 @@ EXE_LIBS = \
     -lradiationModels \
     -lODE \
     -lregionModels \
-    -lsurfaceFilmModels
+    -lsurfaceFilmModels \
+    -lsampling