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

applications/solvers/combustion/fireFoam/fireFoam.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,8 +33,8 @@ Description
 #include "hsCombustionThermo.H"
 #include "turbulenceModel.H"
 #include "combustionModel.H"
-
 #include "radiationModel.H"
+#include "pimpleLoop.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -57,7 +57,7 @@ int main(int argc, char *argv[])
 
     while (runTime.run())
     {
-        #include "readPISOControls.H"
+        #include "readPIMPLEControls.H"
         #include "readTimeControls.H"
         #include "compressibleCourantNo.H"
         #include "setDeltaT.H"
@@ -68,14 +68,13 @@ int main(int argc, char *argv[])
         #include "rhoEqn.H"
 
         // --- Pressure-velocity PIMPLE corrector loop
-        for (int oCorr=0; oCorr<nOuterCorr; oCorr++)
+        for
+        (
+            pimpleLoop pimpleCorr(mesh, nOuterCorr);
+            pimpleCorr.loop();
+            pimpleCorr++
+        )
         {
-            bool finalIter = oCorr == nOuterCorr-1;
-            if (finalIter)
-            {
-                mesh.data::add("finalIteration", true);
-            }
-
             #include "UEqn.H"
 
             #include "ftEqn.H"
@@ -86,11 +85,6 @@ int main(int argc, char *argv[])
             {
                 #include "pEqn.H"
             }
-
-            if (finalIter)
-            {
-                mesh.data::remove("finalIteration");
-            }
         }
 
         turbulence->correct();

applications/solvers/combustion/fireFoam/pEqn.H

@@ -33,7 +33,7 @@ for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
             p_rgh.select
             (
                 (
-                    finalIter
+                    pimpleCorr.finalIter()
                  && corr == nCorr-1
                  && nonOrth == nNonOrthCorr
                 )

applications/solvers/combustion/rhoReactingFoam/pEqn.H

@@ -45,9 +45,9 @@
                     p.select
                     (
                         (
-                            finalIter
+                            pimpleCorr.finalIter()
-                            && corr == nCorr-1
+                         && corr == nCorr-1
-                            && nonOrth == nNonOrthCorr
+                         && nonOrth == nNonOrthCorr
                         )
                     )
                 )
@@ -89,9 +89,9 @@
                     p.select
                     (
                         (
-                            finalIter
+                            pimpleCorr.finalIter()
-                            && corr == nCorr-1
+                         && corr == nCorr-1
-                            && nonOrth == nNonOrthCorr
+                         && nonOrth == nNonOrthCorr
                         )
                     )
                 )

applications/solvers/combustion/rhoReactingFoam/rhoReactingFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2009-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2009-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -36,6 +36,7 @@ Description
 #include "rhoChemistryModel.H"
 #include "chemistrySolver.H"
 #include "multivariateScheme.H"
+#include "pimpleLoop.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -69,14 +70,14 @@ int main(int argc, char *argv[])
         #include "chemistry.H"
         #include "rhoEqn.H"
 
-        for (label oCorr=1; oCorr <= nOuterCorr; oCorr++)
+        // --- Pressure-velocity PIMPLE corrector loop
+        for
+        (
+            pimpleLoop pimpleCorr(mesh, nOuterCorr);
+            pimpleCorr.loop();
+            pimpleCorr++
+        )
         {
-            bool finalIter = oCorr == nOuterCorr-1;
-            if (finalIter)
-            {
-                mesh.data::add("finalIteration", true);
-            }
-
             #include "UEqn.H"
             #include "YEqn.H"
             #include "hsEqn.H"
@@ -86,11 +87,6 @@ int main(int argc, char *argv[])
             {
                 #include "pEqn.H"
             }
-
-            if (finalIter)
-            {
-                mesh.data::remove("finalIteration");
-            }
         }
 
         turbulence->correct();

applications/solvers/compressible/rhoPimpleFoam/Allwmake

@@ -3,6 +3,7 @@ cd ${0%/*} || exit 1    # run from this directory
 set -x
 
 wmake
+wmake rhoLTSPimpleFoam
 wmake rhoPorousMRFPimpleFoam
 
 # ----------------------------------------------------------------- end-of-file

applications/solvers/compressible/rhoPimpleFoam/pEqn.H

@@ -37,11 +37,9 @@ if (transonic)
             (
                 p.select
                 (
-                    (
+                    pimpleCorr.finalIter()
-                        finalIter
+                 && corr == nCorr-1
-                     && corr == nCorr-1
+                 && nonOrth == nNonOrthCorr
-                     && nonOrth == nNonOrthCorr
-                    )
                 )
             )
         );
@@ -77,11 +75,9 @@ else
             (
                 p.select
                 (
-                    (
+                    pimpleCorr.finalIter()
-                        finalIter
+                 && corr == nCorr-1
-                     && corr == nCorr-1
+                 && nonOrth == nNonOrthCorr
-                     && nonOrth == nNonOrthCorr
-                    )
                 )
             )
         );

applications/solvers/compressible/rhoPimpleFoam/rhoLTSPimpleFoam/Make/files

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

applications/solvers/compressible/rhoPimpleFoam/rhoLTSPimpleFoam/Make/options

@@ -0,0 +1,15 @@
+EXE_INC = \
+    -I.. \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
+    -I$(LIB_SRC)/finiteVolume/cfdTools \
+    -I$(LIB_SRC)/finiteVolume/lnInclude
+
+EXE_LIBS = \
+    -lbasicThermophysicalModels \
+    -lspecie \
+    -lcompressibleTurbulenceModel \
+    -lcompressibleRASModels \
+    -lcompressibleLESModels \
+    -lfiniteVolume \
+    -lmeshTools

applications/solvers/compressible/rhoPimpleFoam/rhoLTSPimpleFoam/rhoLTSPimpleFoam.C

@@ -0,0 +1,112 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2011 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    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 Free Software Foundation, either version 3 of the License, or
+    (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
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    rhoLTSPimpleFoam
+
+Description
+    Transient solver for laminar or turbulent flow of compressible fluids
+    for HVAC and similar applications.
+
+    Uses the flexible PIMPLE (PISO-SIMPLE) solution for time-resolved and
+    pseudo-transient simulations with support for local time-stepping for
+    efficient steady-state solution.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "basicPsiThermo.H"
+#include "turbulenceModel.H"
+#include "fvcSmooth.H"
+#include "pimpleLoop.H"
+#include "bound.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "readPIMPLEControls.H"
+    #include "setInitialrDeltaT.H"
+    #include "createFields.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"
+
+        runTime++;
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        #include "setrDeltaT.H"
+
+        #include "rhoEqn.H"
+
+        // --- Pressure-velocity PIMPLE corrector loop
+        for
+        (
+            pimpleLoop pimpleCorr(mesh, nOuterCorr);
+            pimpleCorr.loop();
+            pimpleCorr++
+        )
+        {
+            if (nOuterCorr != 1)
+            {
+                p.storePrevIter();
+                rho.storePrevIter();
+            }
+
+            turbulence->correct();
+
+            #include "UEqn.H"
+            #include "hEqn.H"
+
+            // --- PISO loop
+            for (int corr=0; corr<nCorr; corr++)
+            {
+                #include "pEqn.H"
+            }
+        }
+
+        runTime.write();
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

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

@@ -0,0 +1,19 @@
+scalar maxDeltaT
+(
+    pimple.lookupOrDefault<scalar>("maxDeltaT", GREAT)
+);
+
+volScalarField rDeltaT
+(
+    IOobject
+    (
+        "rDeltaT",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh,
+    1/dimensionedScalar("maxDeltaT", dimTime, maxDeltaT),
+    zeroGradientFvPatchScalarField::typeName
+);

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

@@ -0,0 +1,79 @@
+{
+    scalar maxCo
+    (
+        pimple.lookupOrDefault<scalar>("maxCo", 0.8)
+    );
+
+    scalar rDeltaTSmoothingCoeff
+    (
+        pimple.lookupOrDefault<scalar>("rDeltaTSmoothingCoeff", 0.02)
+    );
+
+    scalar rDeltaTDampingCoeff
+    (
+        pimple.lookupOrDefault<scalar>("rDeltaTDampingCoeff", 1.0)
+    );
+
+    scalar maxDeltaT
+    (
+        pimple.lookupOrDefault<scalar>("maxDeltaT", GREAT)
+    );
+
+    volScalarField rDeltaT0 = rDeltaT;
+
+    // Set the reciprocal time-step from the local Courant number
+    rDeltaT.dimensionedInternalField() = max
+    (
+        1/dimensionedScalar("maxDeltaT", dimTime, maxDeltaT),
+        fvc::surfaceSum(mag(phi))().dimensionedInternalField()
+       /((2*maxCo)*mesh.V()*rho.dimensionedInternalField())
+    );
+
+    if (transonic)
+    {
+        surfaceScalarField phid
+        (
+            "phid",
+            fvc::interpolate(psi)*(fvc::interpolate(U) & mesh.Sf())
+        );
+
+        rDeltaT.dimensionedInternalField() = max
+        (
+            rDeltaT.dimensionedInternalField(),
+            fvc::surfaceSum(mag(phid))().dimensionedInternalField()
+            /((2*maxCo)*mesh.V()*psi.dimensionedInternalField())
+        );
+    }
+
+    // Update tho boundary values of the reciprocal time-step
+    rDeltaT.correctBoundaryConditions();
+
+    Info<< "Flow time scale min/max = "
+        << gMin(1/rDeltaT.internalField())
+        << ", " << gMax(1/rDeltaT.internalField()) << endl;
+
+    if (rDeltaTSmoothingCoeff < 1.0)
+    {
+        fvc::smooth(rDeltaT, rDeltaTSmoothingCoeff);
+    }
+
+    Info<< "Smoothed flow time scale min/max = "
+        << gMin(1/rDeltaT.internalField())
+        << ", " << gMax(1/rDeltaT.internalField()) << endl;
+
+    // Limit rate of change of time scale
+    // - reduce as much as required
+    // - only increase at a fraction of old time scale
+    if
+    (
+        rDeltaTDampingCoeff < 1.0
+     && runTime.timeIndex() > runTime.startTimeIndex() + 1
+    )
+    {
+        rDeltaT = rDeltaT0*max(rDeltaT/rDeltaT0, 1.0 - rDeltaTDampingCoeff);
+
+        Info<< "Damped flow time scale min/max = "
+            << gMin(1/rDeltaT.internalField())
+            << ", " << gMax(1/rDeltaT.internalField()) << endl;
+    }
+}

applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C

@@ -37,6 +37,7 @@ Description
 #include "basicPsiThermo.H"
 #include "turbulenceModel.H"
 #include "bound.H"
+#include "pimpleLoop.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -64,14 +65,13 @@ int main(int argc, char *argv[])
         #include "rhoEqn.H"
 
         // --- Pressure-velocity PIMPLE corrector loop
-        for (int oCorr=0; oCorr<nOuterCorr; oCorr++)
+        for
+        (
+            pimpleLoop pimpleCorr(mesh, nOuterCorr);
+            pimpleCorr.loop();
+            pimpleCorr++
+        )
         {
-            bool finalIter = oCorr == nOuterCorr-1;
-            if (finalIter)
-            {
-                mesh.data::add("finalIteration", true);
-            }
-
             if (nOuterCorr != 1)
             {
                 p.storePrevIter();
@@ -88,11 +88,6 @@ int main(int argc, char *argv[])
             }
 
             turbulence->correct();
-
-            if (finalIter)
-            {
-                mesh.data::remove("finalIteration");
-            }
         }
 
         runTime.write();

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

@@ -40,7 +40,7 @@ if (transonic)
                 p.select
                 (
                     (
-                        finalIter
+                       pimpleCorr.finalIter()
                      && corr == nCorr-1
                      && nonOrth == nNonOrthCorr
                     )
@@ -81,7 +81,7 @@ else
                 p.select
                 (
                     (
-                        finalIter
+                        pimpleCorr.finalIter()
                      && corr == nCorr-1
                      && nonOrth == nNonOrthCorr
                     )

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

@@ -39,6 +39,7 @@ Description
 #include "bound.H"
 #include "MRFZones.H"
 #include "porousZones.H"
+#include "pimpleLoop.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -69,14 +70,13 @@ int main(int argc, char *argv[])
         #include "rhoEqn.H"
 
         // --- Pressure-velocity PIMPLE corrector loop
-        for (int oCorr=0; oCorr<nOuterCorr; oCorr++)
+        for
+        (
+            pimpleLoop pimpleCorr(mesh, nOuterCorr);
+            pimpleCorr.loop();
+            pimpleCorr++
+        )
         {
-            bool finalIter = oCorr == nOuterCorr-1;
-            if (finalIter)
-            {
-                mesh.data::add("finalIteration", true);
-            }
-
             if (nOuterCorr != 1)
             {
                 p.storePrevIter();
@@ -93,11 +93,6 @@ int main(int argc, char *argv[])
             }
 
             turbulence->correct();
-
-            if (finalIter)
-            {
-                mesh.data::remove("finalIteration");
-            }
         }
 
         runTime.write();

applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/TEqn.H

@@ -12,7 +12,7 @@
     );
 
     TEqn.relax();
-    TEqn.solve(mesh.solver(T.select(finalIter)));
+    TEqn.solve();
 
     rhok = 1.0 - beta*(T - TRef);
 }

applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/UEqn.H

@@ -21,7 +21,6 @@
                   - ghf*fvc::snGrad(rhok)
                   - fvc::snGrad(p_rgh)
                 )*mesh.magSf()
-            ),
+            )
-            mesh.solver(U.select(finalIter))
         );
     }

applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/buoyantBoussinesqPimpleFoam.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 "pimpleLoop.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -77,14 +78,13 @@ int main(int argc, char *argv[])
         #include "setDeltaT.H"
 
         // --- Pressure-velocity PIMPLE corrector loop
-        for (int oCorr=0; oCorr<nOuterCorr; oCorr++)
+        for
+        (
+            pimpleLoop pimpleCorr(mesh, nOuterCorr);
+            pimpleCorr.loop();
+            pimpleCorr++
+        )
         {
-            bool finalIter = oCorr == nOuterCorr-1;
-            if (finalIter)
-            {
-                mesh.data::add("finalIteration", true);
-            }
-
             if (nOuterCorr != 1)
             {
                 p_rgh.storePrevIter();
@@ -100,11 +100,6 @@ int main(int argc, char *argv[])
             }
 
             turbulence->correct();
-
-            if (finalIter)
-            {
-                mesh.data::remove("finalIteration");
-            }
         }
 
         runTime.write();

applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/pEqn.H

@@ -26,7 +26,7 @@
                 p_rgh.select
                 (
                     (
-                        finalIter
+                        pimpleCorr.finalIter()
                      && corr == nCorr-1
                      && nonOrth == nNonOrthCorr
                     )

applications/solvers/heatTransfer/buoyantPimpleFoam/UEqn.H

@@ -21,7 +21,6 @@
                   - ghf*fvc::snGrad(rho)
                   - fvc::snGrad(p_rgh)
                 )*mesh.magSf()
-            ),
+            )
-            mesh.solver(U.select(finalIter))
         );
     }

applications/solvers/heatTransfer/buoyantPimpleFoam/buoyantPimpleFoam.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
@@ -37,6 +37,7 @@ Description
 #include "basicRhoThermo.H"
 #include "turbulenceModel.H"
 #include "fixedGradientFvPatchFields.H"
+#include "pimpleLoop.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -70,14 +71,13 @@ int main(int argc, char *argv[])
         #include "rhoEqn.H"
 
         // --- Pressure-velocity PIMPLE corrector loop
-        for (int oCorr=0; oCorr<nOuterCorr; oCorr++)
+        for
+        (
+            pimpleLoop pimpleCorr(mesh, nOuterCorr);
+            pimpleCorr.loop();
+            pimpleCorr++
+        )
         {
-            bool finalIter = oCorr == nOuterCorr-1;
-            if (finalIter)
-            {
-                mesh.data::add("finalIteration", true);
-            }
-
             if (nOuterCorr != 1)
             {
                 p_rgh.storePrevIter();
@@ -95,11 +95,6 @@ int main(int argc, char *argv[])
             turbulence->correct();
 
             rho = thermo.rho();
-
-            if (finalIter)
-            {
-                mesh.data::remove("finalIteration");
-            }
         }
 
         runTime.write();

applications/solvers/heatTransfer/buoyantPimpleFoam/hEqn.H

@@ -9,7 +9,7 @@
     );
 
     hEqn.relax();
-    hEqn.solve(mesh.solver(h.select(finalIter)));
+    hEqn.solve();
 
     thermo.correct();
 }

applications/solvers/heatTransfer/buoyantPimpleFoam/pEqn.H

@@ -40,7 +40,7 @@
                 p_rgh.select
                 (
                     (
-                        finalIter
+                        pimpleCorr.finalIter()
                      && corr == nCorr-1
                      && nonOrth == nNonOrthCorr
                     )

applications/solvers/incompressible/pimpleFoam/pEqn.H

@@ -27,11 +27,9 @@ for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
         (
             p.select
             (
-                (
+                pimpleCorr.finalIter()
-                    finalIter
+             && corr == nCorr-1
-                 && corr == nCorr-1
+             && nonOrth == nNonOrthCorr
-                 && nonOrth == nNonOrthCorr
-                )
             )
         )
     );

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

@@ -30,7 +30,7 @@ for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
             p.select
             (
                 (
-                    finalIter
+                    pimpleCorr.finalIter()
                  && corr == nCorr-1
                  && nonOrth == nNonOrthCorr
                 )

applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pimpleDyMFoam.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
@@ -36,6 +36,7 @@ Description
 #include "singlePhaseTransportModel.H"
 #include "turbulenceModel.H"
 #include "dynamicFvMesh.H"
+#include "pimpleLoop.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -83,15 +84,14 @@ int main(int argc, char *argv[])
             #include "meshCourantNo.H"
         }
 
-        // --- PIMPLE loop
+        // --- Pressure-velocity PIMPLE corrector loop
-        for (int oCorr=0; oCorr<nOuterCorr; oCorr++)
+        for
+        (
+            pimpleLoop pimpleCorr(mesh, nOuterCorr);
+            pimpleCorr.loop();
+            pimpleCorr++
+        )
         {
-            bool finalIter = oCorr == nOuterCorr-1;
-            if (finalIter)
-            {
-                mesh.data::add("finalIteration", true);
-            }
-
             if (nOuterCorr != 1)
             {
                 p.storePrevIter();
@@ -106,11 +106,6 @@ int main(int argc, char *argv[])
             }
 
             turbulence->correct();
-
-            if (finalIter)
-            {
-                mesh.data::remove("finalIteration");
-            }
         }
 
         runTime.write();

applications/solvers/incompressible/pimpleFoam/pimpleFoam.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 "fvCFD.H"
 #include "singlePhaseTransportModel.H"
 #include "turbulenceModel.H"
+#include "pimpleLoop.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -60,14 +61,13 @@ int main(int argc, char *argv[])
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
         // --- Pressure-velocity PIMPLE corrector loop
-        for (int oCorr=0; oCorr<nOuterCorr; oCorr++)
+        for
+        (
+            pimpleLoop pimpleCorr(mesh, nOuterCorr);
+            pimpleCorr.loop();
+            pimpleCorr++
+        )
         {
-            bool finalIter = oCorr == nOuterCorr-1;
-            if (finalIter)
-            {
-                mesh.data::add("finalIteration", true);
-            }
-
             if (nOuterCorr != 1)
             {
                 p.storePrevIter();
@@ -82,11 +82,6 @@ int main(int argc, char *argv[])
             }
 
             turbulence->correct();
-
-            if (finalIter)
-            {
-                mesh.data::remove("finalIteration");
-            }
         }
 
         runTime.write();

applications/solvers/lagrangian/LTSReactingParcelFoam/LTSReactingParcelFoam.C

@@ -46,6 +46,7 @@ Description
 #include "timeActivatedExplicitSource.H"
 #include "SLGThermo.H"
 #include "fvcSmooth.H"
+#include "pimpleLoop.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -71,7 +72,7 @@ int main(int argc, char *argv[])
 
     while (runTime.run())
     {
-        #include "readPISOControls.H"
+        #include "readPIMPLEControls.H"
         #include "readChemistryProperties.H"
         #include "readAdditionalSolutionControls.H"
         #include "readTimeControls.H"
@@ -86,17 +87,32 @@ int main(int argc, char *argv[])
         #include "timeScales.H"
 
         #include "rhoEqn.H"
-        #include "UEqn.H"
-        #include "YEqn.H"
-        #include "hsEqn.H"
 
-        // --- PISO loop
+        // --- Pressure-velocity PIMPLE corrector loop
-        for (int corr=0; corr<nCorr; corr++)
+        for
+        (
+            pimpleLoop pimpleCorr(mesh, nOuterCorr);
+            pimpleCorr.loop();
+            pimpleCorr++
+        )
         {
-            #include "pEqn.H"
+            if (nOuterCorr != 1)
-        }
+            {
+                p.storePrevIter();
+            }
 
-        turbulence->correct();
+            turbulence->correct();
+
+            #include "UEqn.H"
+            #include "YEqn.H"
+            #include "hsEqn.H"
+
+            // --- PISO loop
+            for (int corr=0; corr<nCorr; corr++)
+            {
+                #include "pEqn.H"
+            }
+        }
 
         if (runTime.write())
         {

applications/solvers/lagrangian/LTSReactingParcelFoam/YEqn.H

@@ -1,4 +1,3 @@
-
 tmp<fv::convectionScheme<scalar> > mvConvection
 (
     fv::convectionScheme<scalar>::New
@@ -10,7 +9,6 @@ tmp<fv::convectionScheme<scalar> > mvConvection
     )
 );
 
-
 if (solveSpecies)
 {
     label inertIndex = -1;

applications/solvers/lagrangian/LTSReactingParcelFoam/createFields.H

@@ -73,12 +73,12 @@
 
     dimensionedScalar rhoMax
     (
-        mesh.solutionDict().subDict("PISO").lookup("rhoMax")
+        mesh.solutionDict().subDict("PIMPLE").lookup("rhoMax")
     );
 
     dimensionedScalar rhoMin
     (
-        mesh.solutionDict().subDict("PISO").lookup("rhoMin")
+        mesh.solutionDict().subDict("PIMPLE").lookup("rhoMin")
     );
 
     Info<< "Creating turbulence model\n" << endl;
@@ -116,11 +116,11 @@
         dimensionedScalar("chemistrySh", dimEnergy/dimTime/dimVolume, 0.0)
     );
 
-    volScalarField invTauFlow
+    volScalarField invTau
     (
         IOobject
         (
-            "invTauFlow",
+            "invTau",
             runTime.timeName(),
             mesh,
             IOobject::READ_IF_PRESENT,

applications/solvers/lagrangian/LTSReactingParcelFoam/pEqn.H

@@ -40,7 +40,18 @@
           - fvm::laplacian(rho*rAU, p)
         );
 
-        pEqn.solve();
+        pEqn.solve
+        (
+            mesh.solver
+            (
+                p.select
+                (
+                    pimpleCorr.finalIter()
+                 && corr == nCorr-1
+                 && nonOrth == nNonOrthCorr
+                )
+            )
+        );
 
         if (nonOrth == nNonOrthCorr)
         {
@@ -48,7 +59,7 @@
         }
     }
 
-    Info<< "p min/max = " << min(p).value() << ", " << max(p).value() << endl;
+    p.relax();
 
     // Second part of thermodynamic density update
     thermo.rho() += psi*p;
@@ -64,4 +75,6 @@
     rho = min(rho, rhoMax);
 
     #include "setPressureWork.H"
+
+    Info<< "p min/max = " << min(p).value() << ", " << max(p).value() << endl;
 }

applications/solvers/lagrangian/LTSReactingParcelFoam/readAdditionalSolutionControls.H

@@ -1,7 +1,9 @@
 dictionary additional = mesh.solutionDict().subDict("additional");
 
-bool eWork = additional.lookupOrDefault("eWork", true);
+// pressure work term for enthalpy equation
-bool hWork = additional.lookupOrDefault("hWork", true);
+bool pressureWork = additional.lookupOrDefault("pressureWork", true);
+bool pressureWorkTimeDerivative =
+    additional.lookupOrDefault("pressureWorkTimeDerivative", true);
 
 // flag to activate solve transport for each specie (Y vector)
 bool solveSpecies = additional.lookupOrDefault("solveSpecies", true);

applications/solvers/lagrangian/LTSReactingParcelFoam/setPressureWork.H

@@ -1,10 +1,13 @@
 DpDt == dimensionedScalar("zero", DpDt.dimensions(), 0.0);
 
-if (eWork)
+if (pressureWork)
 {
-    DpDt += -p*fvc::div(phi/fvc::interpolate(rho));
+    surfaceScalarField phiU("phiU", phi/fvc::interpolate(rho));
-}
+
-if (hWork)
+    DpDt += fvc::div(phiU*fvc::interpolate(p)) - p*fvc::div(phiU);
-{
+
-    DpDt += fvc::div(phi/fvc::interpolate(rho)*fvc::interpolate(p));
+    if (pressureWorkTimeDerivative)
+    {
+        DpDt += fvc::ddt(p);
+    }
 }

applications/solvers/lagrangian/LTSReactingParcelFoam/timeScales.H

@@ -28,129 +28,40 @@ Info<< "Time scales min/max:" << endl;
 
 {
     // Cache old time scale field
-    tmp<volScalarField> tinvTauFlow0
+    tmp<volScalarField> tinvTau0
     (
         new volScalarField
         (
             IOobject
             (
-                "invTauFlow0",
+                "invTau0",
                 runTime.timeName(),
                 mesh,
                 IOobject::NO_READ,
                 IOobject::NO_WRITE,
                 false
             ),
-            invTauFlow
+            invTau
         )
     );
-    const volScalarField& invTauFlow0 = tinvTauFlow0();
+    const volScalarField& invTau0 = tinvTau0();
 
 
     // Flow time scale
     // ~~~~~~~~~~~~~~~
     {
-        invTauFlow =
+        invTau =
             fvc::surfaceSum
             (
                 mag(phi)*mesh.deltaCoeffs()/(maxCo*mesh.magSf())
             )
            /rho;
 
-        invTauFlow.max(1.0/maxDeltaT);
+        invTau.max(1.0/maxDeltaT);
 
         Info<< "    Flow        = "
-            << gMin(1/invTauFlow.internalField()) << ", "
+            << gMin(1/invTau.internalField()) << ", "
-            << gMax(1/invTauFlow.internalField()) << endl;
+            << gMax(1/invTau.internalField()) << endl;
-    }
-
-
-    // Mass source time scale
-    // ~~~~~~~~~~~~~~~~~~~~~~
-
-    {
-        scalarField tau
-        (
-            runTime.deltaTValue()*mag(parcels.Srho() + massSource.SuTot())
-        );
-
-        tau = alphaTauRho*rho/(tau + ROOTVSMALL);
-
-        Info<< "    Density     = "
-            << min(maxDeltaT, gMin(tau)) << ", "
-            << min(maxDeltaT, gMax(tau)) << endl;
-
-        invTauFlow.internalField() = max(invTauFlow.internalField(), 1/tau);
-    }
-
-
-    // Momentum source time scale
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-    {
-/*
-        // Method 1 - mag(U) limit using 'small' nominal velocity
-        scalarField tau
-        (
-            runTime.deltaTValue()
-           *mag
-            (
-                rho.dimensionedInternalField()*g
-              + parcels.UTrans()/(mesh.V()*runTime.deltaT())
-              + momentumSource.Su()
-            )
-           /rho
-        );
-
-        const scalar nomMagU(dimensionedScalar("1", dimVelocity, 1));
-        tau = alphaTauU*(nomMagU + mag(U))/(tau + ROOTVSMALL);
-*/
-/*
-        // Method 2 - based on fluxes and Co-like limit
-        volVectorField UEqnRhs
-        (
-            IOobject
-            (
-                "UEqnRhs",
-                runTime.timeName(),
-                mesh,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE,
-                false
-            ),
-            mesh,
-            dimensionedVector("zero", dimDensity*dimAcceleration, vector::zero),
-            zeroGradientFvPatchVectorField::typeName
-        );
-
-        UEqnRhs.internalField() =
-            rho.dimensionedInternalField()*g
-          + parcels.UTrans()/(mesh.V()*runTime.deltaT())
-          + momentumSource.Su();
-
-        surfaceScalarField phiSU
-        (
-            "phiSU",
-            fvc::interpolate(runTime.deltaT()*UEqnRhs) & mesh.Sf()
-        );
-
-        scalarField tau
-        (
-            alphaTauU*rho
-           /fvc::surfaceSum
-            (
-                mag(phi + phiSU)*mesh.deltaCoeffs()/mesh.magSf()
-              + dimensionedScalar("SMALL", dimDensity/dimTime, ROOTVSMALL)
-            )
-        );
-
-*/
-/*
-        Info<< "    Momentum    = " << min(maxDeltaT, gMin(tau)) << ", "
-            << min(maxDeltaT, gMax(tau)) << endl;
-
-        invTauFlow.internalField() = max(invTauFlow.internalField(), 1/tau);
-*/
     }
 
 
@@ -176,42 +87,7 @@ Info<< "Time scales min/max:" << endl;
         Info<< "    Temperature = " << min(maxDeltaT, gMin(tau)) << ", "
             << min(maxDeltaT, gMax(tau)) << endl;
 
-        invTauFlow.internalField() = max(invTauFlow.internalField(), 1/tau);
+        invTau.internalField() = max(invTau.internalField(), 1/tau);
-    }
-
-
-    // Specie source time scale
-    // ~~~~~~~~~~~~~~~~~~~~~~~~
-
-    {
-        scalarField tau(mesh.nCells(), ROOTVGREAT);
-        forAll(Y, fieldI)
-        {
-            const volScalarField& Yi = Y[fieldI];
-            const scalarField deltaYi
-            (
-                runTime.deltaTValue()
-               *mag
-                (
-                    kappa*chemistry.RR(fieldI)()
-                  + massSource.Su(fieldI)
-                  + parcels.Srho(fieldI)
-                )
-               /rho
-            );
-            tau =
-                min
-                (
-                    tau,
-                    alphaTauSpecie
-                   /(deltaYi/(Yi + specieMaxUnbound) + ROOTVSMALL)
-                );
-        }
-
-        Info<< "    Specie      = " << min(maxDeltaT, gMin(tau)) << ", "
-            << min(maxDeltaT, gMax(tau)) << endl;
-
-        invTauFlow.internalField() = max(invTauFlow.internalField(), 1/tau);
     }
 
 
@@ -219,21 +95,21 @@ Info<< "Time scales min/max:" << endl;
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     // - reduce as much as required for flow, but limit source contributions
     const dimensionedScalar deltaTRamp("deltaTRamp", dimless, 1/(1 + 0.2));
-    invTauFlow = max(invTauFlow, invTauFlow0*deltaTRamp);
+    invTau = max(invTau, invTau0*deltaTRamp);
-    tinvTauFlow0.clear();
+    tinvTau0.clear();
 
 
     // Limit the largest time scale
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    invTauFlow.max(1/maxDeltaT);
+    invTau.max(1/maxDeltaT);
 
 
     // Spatially smooth the time scale field
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    fvc::smooth(invTauFlow, alphaTauSmooth);
+    fvc::smooth(invTau, alphaTauSmooth);
 
-    Info<< "    Overall     = " << min(1/invTauFlow).value()
+    Info<< "    Overall     = " << min(1/invTau).value()
-        << ", " << max(1/invTauFlow).value() << nl << endl;
+        << ", " << max(1/invTau).value() << nl << endl;
 }
 
 

applications/solvers/lagrangian/coalChemistryFoam/UEqn.H

@@ -13,5 +13,5 @@
 
     if (momentumPredictor)
     {
-        solve(UEqn == -fvc::grad(p), mesh.solver(U.select(finalIter)));
+        solve(UEqn == -fvc::grad(p));
     }

applications/solvers/lagrangian/coalChemistryFoam/coalChemistryFoam.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
@@ -45,6 +45,7 @@ Description
 #include "timeActivatedExplicitSource.H"
 #include "radiationModel.H"
 #include "SLGThermo.H"
+#include "pimpleLoop.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -89,15 +90,14 @@ int main(int argc, char *argv[])
         #include "chemistry.H"
         #include "rhoEqn.H"
 
-        // --- PIMPLE loop
+        // --- Pressure-velocity PIMPLE corrector loop
-        for (int oCorr=0; oCorr<nOuterCorr; oCorr++)
+        for
+        (
+            pimpleLoop pimpleCorr(mesh, nOuterCorr);
+            pimpleCorr.loop();
+            pimpleCorr++
+        )
         {
-            bool finalIter = oCorr == nOuterCorr - 1;
-            if (finalIter)
-            {
-                mesh.data::add("finalIteration", true);
-            }
-
             #include "UEqn.H"
             #include "YEqn.H"
             #include "hsEqn.H"
@@ -107,11 +107,6 @@ int main(int argc, char *argv[])
             {
                 #include "pEqn.H"
             }
-
-            if (finalIter)
-            {
-                mesh.data::remove("finalIteration");
-            }
         }
 
         turbulence->correct();

applications/solvers/lagrangian/coalChemistryFoam/hsEqn.H

@@ -15,7 +15,7 @@
 
     hsEqn.relax();
 
-    hsEqn.solve(mesh.solver(hs.select(finalIter)));
+    hsEqn.solve();
 
     thermo.correct();
 

applications/solvers/lagrangian/coalChemistryFoam/pEqn.H

@@ -32,11 +32,9 @@ if (transonic)
             (
                 p.select
                 (
-                    (
+                    pimpleCorr.finalIter()
-                        finalIter
+                 && corr == nCorr-1
-                     && corr == nCorr-1
+                 && nonOrth == nNonOrthCorr
-                     && nonOrth == nNonOrthCorr
-                    )
                 )
             )
         );
@@ -73,11 +71,9 @@ else
             (
                 p.select
                 (
-                    (
+                    pimpleCorr.finalIter()
-                        finalIter
+                 && corr == nCorr-1
-                     && corr == nCorr-1
+                 && nonOrth == nNonOrthCorr
-                     && nonOrth == nNonOrthCorr
-                    )
                 )
             )
         );

applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/YEqn.H

@@ -11,6 +11,7 @@ tmp<fv::convectionScheme<scalar> > mvConvection
 );
 
 
+if (solveSpecies)
 {
     label inertIndex = -1;
     volScalarField Yt(0.0*Y[0]);

applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/hsEqn.H

@@ -16,17 +16,16 @@
         )
     );
 
-    if (dpdt)
+    if (pressureWork)
     {
-        pWork() += fvc::ddt(p);
+        surfaceScalarField phiU("phiU", phi/fvc::interpolate(rho));
-    }
+
-    if (eWork)
+        pWork() += fvc::div(phiU*fvc::interpolate(p)) - p*fvc::div(phiU);
-    {
+
-        pWork() = -p*fvc::div(phi/fvc::interpolate(rho));
+        if (pressureWorkTimeDerivative)
-    }
+        {
-    if (hWork)
+            pWork() += fvc::ddt(p);
-    {
+        }
-        pWork() += fvc::div(phi/fvc::interpolate(rho)*fvc::interpolate(p));
     }
 
     {

applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/pEqn.H

@@ -40,14 +40,17 @@
           - fvm::laplacian(rho*rAU, p)
         );
 
-        if (corr == nCorr-1 && nonOrth == nNonOrthCorr)
+        pEqn.solve
-        {
+        (
-            pEqn.solve(mesh.solver("pFinal"));
+            mesh.solver
-        }
+            (
-        else
+                p.select
-        {
+                (
-            pEqn.solve();
+                    corr == nCorr-1
-        }
+                 && nonOrth == nNonOrthCorr
+                )
+            )
+        );
 
         if (nonOrth == nNonOrthCorr)
         {

applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/readAdditionalSolutionControls.H

@@ -1,5 +1,9 @@
 dictionary additional = mesh.solutionDict().subDict("additional");
 
-bool dpdt = additional.lookupOrDefault("dpdt", true);
+// pressure work term for enthalpy equation
-bool eWork = additional.lookupOrDefault("eWork", true);
+bool pressureWork = additional.lookupOrDefault("pressureWork", true);
-bool hWork = additional.lookupOrDefault("hWork", true);
+bool pressureWorkTimeDerivative =
+    additional.lookupOrDefault("pressureWorkTimeDerivative", true);
+
+// flag to activate solve transport for each specie (Y vector)
+bool solveSpecies = additional.lookupOrDefault("solveSpecies", true);

applications/solvers/lagrangian/reactingParcelFoam/pEqn.H

@@ -26,7 +26,18 @@ if (transonic)
             parcels.Srho()
         );
 
-        pEqn.solve();
+        pEqn.solve
+        (
+            mesh.solver
+            (
+                p.select
+                (
+                    pimpleCorr.finalIter()
+                 && corr == nCorr-1
+                 && nonOrth == nNonOrthCorr
+                )
+            )
+        );
 
         if (nonOrth == nNonOrthCorr)
         {
@@ -54,7 +65,18 @@ else
             parcels.Srho()
         );
 
-        pEqn.solve();
+        pEqn.solve
+        (
+            mesh.solver
+            (
+                p.select
+                (
+                    pimpleCorr.finalIter()
+                 && corr == nCorr-1
+                 && nonOrth == nNonOrthCorr
+                )
+            )
+        );
 
         if (nonOrth == nNonOrthCorr)
         {
@@ -69,4 +91,5 @@ else
 U -= rAU*fvc::grad(p);
 U.correctBoundaryConditions();
 
+
 DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);

applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.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
@@ -38,6 +38,7 @@ Description
 #include "chemistrySolver.H"
 #include "radiationModel.H"
 #include "SLGThermo.H"
+#include "pimpleLoop.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -77,24 +78,29 @@ int main(int argc, char *argv[])
         #include "chemistry.H"
         #include "rhoEqn.H"
 
-        // --- PIMPLE loop
+        // --- Pressure-velocity PIMPLE corrector loop
-        for (int ocorr=1; ocorr<=nOuterCorr; ocorr++)
+        for
+        (
+            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=0; corr<nCorr; corr++)
             {
-                #include "hsEqn.H"
                 #include "pEqn.H"
             }
+
+            turbulence->correct();
+
+            rho = thermo.rho();
         }
 
-        turbulence->correct();
-
-        rho = thermo.rho();
-
         if (runTime.write())
         {
             chemistry.dQ()().write();

applications/solvers/multiphase/interFoam/LTSInterFoam/setrDeltaT.H

@@ -11,7 +11,7 @@
 
     scalar rDeltaTSmoothingCoeff
     (
-        piso.lookupOrDefault<scalar>("rDeltaTSmoothingCoeff", 0.1)
+        piso.lookupOrDefault<scalar>("rDeltaTSmoothingCoeff", 0.02)
     );
 
     label nAlphaSpreadIter

applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeFoam.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
@@ -46,6 +46,7 @@ Description
 #include "interfaceProperties.H"
 #include "phaseChangeTwoPhaseMixture.H"
 #include "turbulenceModel.H"
+#include "pimpleLoop.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -83,10 +84,13 @@ int main(int argc, char *argv[])
         turbulence->correct();
 
         // --- Pressure-velocity PIMPLE corrector loop
-        for (int oCorr=0; oCorr<nOuterCorr; oCorr++)
+        for
+        (
+            pimpleLoop pimpleCorr(mesh, nOuterCorr);
+            pimpleCorr.loop();
+            pimpleCorr++
+        )
         {
-            bool finalIter = oCorr == nOuterCorr-1;
-
             #include "UEqn.H"
 
             // --- PISO loop

applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H

@@ -39,7 +39,7 @@
             (
                 p_rgh.select
                 (
-                    finalIter
+                    pimpleCorr.finalIter()
                  && corr == nCorr-1
                  && nonOrth == nNonOrthCorr
                 )

applications/solvers/multiphase/twoLiquidMixingFoam/pEqn.H

@@ -29,7 +29,7 @@
             (
                 p_rgh.select
                 (
-                    finalIter
+                    pimpleCorr.finalIter()
                  && corr == nCorr-1
                  && nonOrth == nNonOrthCorr
                 )

applications/solvers/multiphase/twoLiquidMixingFoam/twoLiquidMixingFoam.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
@@ -34,6 +34,7 @@ Description
 #include "fvCFD.H"
 #include "twoPhaseMixture.H"
 #include "turbulenceModel.H"
+#include "pimpleLoop.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -66,14 +67,13 @@ int main(int argc, char *argv[])
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
         // --- Pressure-velocity PIMPLE corrector loop
-        for (int oCorr=0; oCorr<nOuterCorr; oCorr++)
+        for
+        (
+            pimpleLoop pimpleCorr(mesh, nOuterCorr);
+            pimpleCorr.loop();
+            pimpleCorr++
+        )
         {
-            bool finalIter = oCorr == nOuterCorr-1;
-            if (finalIter)
-            {
-                mesh.data::add("finalIteration", true);
-            }
-
             twoPhaseProperties.correct();
 
             #include "alphaEqn.H"
@@ -87,11 +87,6 @@ int main(int argc, char *argv[])
             }
 
             turbulence->correct();
-
-            if (finalIter)
-            {
-                mesh.data::remove("finalIteration");
-            }
         }
 
         runTime.write();

src/finiteVolume/cfdTools/general/pimpleLoop/pimpleLoop.H

@@ -0,0 +1,128 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    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 Free Software Foundation, either version 3 of the License, or
+    (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
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::pimpleLoop
+
+Description
+    PIMPLE loop class to formalise the iteration and automate the handling
+    of the "finalIteration" mesh data entry.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef pimpleLoop_H
+#define pimpleLoop_H
+
+#include "fvMesh.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class pimpleLoop Declaration
+\*---------------------------------------------------------------------------*/
+
+class pimpleLoop
+{
+    // Private data
+
+        //- Reference to the mesh
+        fvMesh& mesh_;
+
+        //- Number of PIMPLE correctors
+        const int nCorr_;
+
+        //- Current PIMPLE corrector
+        int corr_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        pimpleLoop(const pimpleLoop&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const pimpleLoop&);
+
+
+public:
+
+    // Constructors
+
+        //- Construct from components
+        pimpleLoop(fvMesh& mesh, const int nCorr)
+        :
+            mesh_(mesh),
+            nCorr_(nCorr),
+            corr_(0)
+        {}
+
+
+    //- Destructor
+    ~pimpleLoop()
+    {}
+
+
+    // Member Functions
+
+        bool loop()
+        {
+            if (finalIter())
+            {
+                mesh_.data::add("finalIteration", true);
+            }
+
+            return corr_ < nCorr_;
+        }
+
+        bool finalIter() const
+        {
+            return corr_ == nCorr_-1;
+        }
+
+
+    // Member Operators
+
+        void operator++(int)
+        {
+            if (finalIter())
+            {
+                mesh_.data::remove("finalIteration");
+            }
+
+            corr_++;
+        }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C

@@ -332,7 +332,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
     // ~~~~~~~~~~~~~
 
     // Calculate new particle temperature
-    scalar Cuh = 0.0;
+    scalar Sph = 0.0;
     scalar T1 =
         this->calcHeatTransfer
         (
@@ -349,7 +349,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
             NCpW,
             Sh,
             dhsTrans,
-            Cuh
+            Sph
         );
 
 
@@ -423,7 +423,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
         td.cloud().hsTrans()[cellI] += np0*dhsTrans;
 
         // Update sensible enthalpy coefficient
-        td.cloud().hsCoeff()[cellI] += np0*Cuh*this->areaS();
+        td.cloud().hsCoeff()[cellI] += np0*Sph;
     }
 
 

src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C

@@ -342,7 +342,7 @@ void Foam::ReactingParcel<ParcelType>::calc
     // ~~~~~~~~~~~~~
 
     // Calculate new particle temperature
-    scalar Cuh = 0.0;
+    scalar Sph = 0.0;
     scalar T1 =
         this->calcHeatTransfer
         (
@@ -359,7 +359,7 @@ void Foam::ReactingParcel<ParcelType>::calc
             NCpW,
             Sh,
             dhsTrans,
-            Cuh
+            Sph
         );
 
 
@@ -410,7 +410,7 @@ void Foam::ReactingParcel<ParcelType>::calc
         td.cloud().hsTrans()[cellI] += np0*dhsTrans;
 
         // Update sensible enthalpy coefficient
-        td.cloud().hsCoeff()[cellI] += np0*Cuh*this->areaS();
+        td.cloud().hsCoeff()[cellI] += np0*Sph;
     }
 
 

tutorials/combustion/fireFoam/les/smallPoolFire2D/system/fvSolution

@@ -92,7 +92,7 @@ solvers
 
 }
 
-PISO
+PIMPLE
 {
     momentumPredictor yes;
     nOuterCorrectors  1;

tutorials/combustion/fireFoam/les/smallPoolFire3D/system/fvSolution

@@ -83,7 +83,7 @@ solvers
 
 }
 
-PISO
+PIMPLE
 {
     momentumPredictor no;
     nOuterCorrectors  1;

tutorials/compressible/rhoLTSPimpleFoam/angledDuct/0/T

@@ -0,0 +1,54 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      T;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 1 0 0 0];
+
+internalField   uniform 293;
+
+boundaryField
+{
+
+    front
+    {
+        type            zeroGradient;
+    }
+    back
+    {
+        type            zeroGradient;
+    }
+    wall
+    {
+        type            zeroGradient;
+    }
+    porosityWall
+    {
+        type            zeroGradient;
+    }
+
+    inlet
+    {
+        type            fixedValue;
+        value           $internalField;
+    }
+    outlet
+    {
+        type            inletOutlet;
+        value           $internalField;
+        inletValue      $internalField;
+    }
+}
+
+// ************************************************************************* //

tutorials/compressible/rhoLTSPimpleFoam/angledDuct/0/U

@@ -0,0 +1,57 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volVectorField;
+    object      U;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (47.8 0 0);
+
+boundaryField
+{
+    front
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+    back
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+    wall
+    {
+        type            fixedValue;
+        value           uniform (0 0 0);
+    }
+    porosityWall
+    {
+        type            slip;
+        value           uniform (0 0 0);
+    }
+    inlet
+    {
+        type            flowRateInletVelocity;
+        flowRate        0.1;
+        value           uniform (0 0 0);
+    }
+    outlet
+    {
+        type            fluxCorrectedVelocity; //inletOutlet;
+        value           uniform (0 0 0);
+        inletValue      uniform (0 0 0);
+    }
+}
+
+// ************************************************************************* //

tutorials/compressible/rhoLTSPimpleFoam/angledDuct/0/alphat

@@ -0,0 +1,57 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      alphat;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -1 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    front
+    {
+        type            alphatWallFunction;
+        value           uniform 0;
+    }
+    back
+    {
+        type            alphatWallFunction;
+        value           uniform 0;
+    }
+    wall
+    {
+        type            alphatWallFunction;
+        value           uniform 0;
+    }
+    porosityWall
+    {
+        type            alphatWallFunction;
+        value           uniform 0;
+    }
+    inlet
+    {
+        type            calculated;
+        value           uniform 0;
+    }
+    outlet
+    {
+        type            calculated;
+        value           uniform 0;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/compressible/rhoLTSPimpleFoam/angledDuct/0/epsilon

@@ -0,0 +1,59 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      epsilon;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -3 0 0 0 0];
+
+internalField   uniform 200;
+
+boundaryField
+{
+    front
+    {
+        type            compressible::epsilonWallFunction;
+        value           uniform 200;
+    }
+    back
+    {
+        type            compressible::epsilonWallFunction;
+        value           uniform 200;
+    }
+    wall
+    {
+        type            compressible::epsilonWallFunction;
+        value           uniform 200;
+    }
+    porosityWall
+    {
+        type            compressible::epsilonWallFunction;
+        value           uniform 200;
+    }
+    inlet
+    {
+        type            compressible::turbulentMixingLengthDissipationRateInlet;
+        mixingLength    0.005;
+        value           uniform 200;
+    }
+    outlet
+    {
+        type            inletOutlet;
+        inletValue      uniform 200;
+        value           uniform 200;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/compressible/rhoLTSPimpleFoam/angledDuct/0/k

@@ -0,0 +1,59 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      k;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -2 0 0 0 0];
+
+internalField   uniform 1;
+
+boundaryField
+{
+    front
+    {
+        type            compressible::kqRWallFunction;
+        value           uniform 1;
+    }
+    back
+    {
+        type            compressible::kqRWallFunction;
+        value           uniform 1;
+    }
+    wall
+    {
+        type            compressible::kqRWallFunction;
+        value           uniform 1;
+    }
+    porosityWall
+    {
+        type            compressible::kqRWallFunction;
+        value           uniform 1;
+    }
+    inlet
+    {
+        type            turbulentIntensityKineticEnergyInlet;
+        intensity       0.05;
+        value           uniform 1;
+    }
+    outlet
+    {
+        type            inletOutlet;
+        inletValue      uniform 1;
+        value           uniform 1;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/compressible/rhoLTSPimpleFoam/angledDuct/0/mut

@@ -0,0 +1,57 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      mut;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -1 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    front
+    {
+        type            mutkWallFunction;
+        value           uniform 0;
+    }
+    back
+    {
+        type            mutkWallFunction;
+        value           uniform 0;
+    }
+    wall
+    {
+        type            mutkWallFunction;
+        value           uniform 0;
+    }
+    porosityWall
+    {
+        type            mutkWallFunction;
+        value           uniform 0;
+    }
+    inlet
+    {
+        type            calculated;
+        value           uniform 0;
+    }
+    outlet
+    {
+        type            calculated;
+        value           uniform 0;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/compressible/rhoLTSPimpleFoam/angledDuct/0/p

@@ -0,0 +1,51 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      p;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -2 0 0 0 0];
+
+internalField   uniform 1.0e5;
+
+boundaryField
+{
+    front
+    {
+        type            zeroGradient;
+    }
+    back
+    {
+        type            zeroGradient;
+    }
+    wall
+    {
+        type            zeroGradient;
+    }
+    porosityWall
+    {
+        type            zeroGradient;
+    }
+
+    inlet
+    {
+        type            zeroGradient;
+    }
+    outlet
+    {
+        type            fixedValue;
+        value           $internalField;
+    }
+}
+
+// ************************************************************************* //

tutorials/compressible/rhoLTSPimpleFoam/angledDuct/Allrun

@@ -0,0 +1,9 @@
+cd ${0%/*} || exit 1    # run from this directory
+
+m4 constant/polyMesh/blockMeshDict.m4 > constant/polyMesh/blockMeshDict
+
+# Source tutorial run functions
+. $WM_PROJECT_DIR/bin/tools/RunFunctions
+
+runApplication blockMesh
+runApplication `getApplication`

tutorials/compressible/rhoLTSPimpleFoam/angledDuct/constant/RASProperties

@@ -0,0 +1,25 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      RASProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+RASModel        kEpsilon;
+
+turbulence      on;
+
+printCoeffs     on;
+
+
+// ************************************************************************* //

tutorials/compressible/rhoLTSPimpleFoam/angledDuct/constant/polyMesh/blockMeshDict.m4

@@ -0,0 +1,165 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    `format'      ascii;
+    class       dictionary;
+    object      blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// block definition for a porosity with an angled inlet/outlet
+// the porosity is not aligned with the main axes
+//
+dnl> -----------------------------------------------------------------
+dnl> <STANDARD DEFINTIONS>
+dnl>
+changecom(//)changequote([,]) dnl>
+define(calc, [esyscmd(perl -e 'print ($1)')]) dnl>
+define(VCOUNT, 0)  dnl>
+define(vlabel, [[// ]pt VCOUNT ($1) define($1, VCOUNT)define([VCOUNT], incr(VCOUNT))])  dnl>
+dnl>
+define(hex2D, hex ($1b $2b $3b $4b $1f $2f $3f $4f)) dnl>
+define(quad2D, ($1f $1b $2b $2f))  dnl>
+define(frontQuad, ($1f $2f $3f $4f)) dnl>
+define(backQuad, ($4b $3b $2b $1b)) dnl>
+dnl>
+dnl> </STANDARD DEFINTIONS>
+dnl> -----------------------------------------------------------------
+dnl>
+define(ncells, 20) dnl>
+define(ninlet, 15) dnl>
+define(nporo, 20) dnl>
+define(noutlet, 20) dnl>
+dnl>
+define(x0,0) dnl>
+define(y0,0) dnl>
+define(y0,0) dnl>
+define(Cos,0.7071067812) dnl>   == cos(45)
+define(Sin,0.7071067812) dnl>   == sin(45)
+dnl>
+define(width,50) dnl>
+define(zBack,calc(-width/2)) dnl>
+define(zFront,calc(width/2)) dnl>
+define(leninlet,150)dnl>
+define(lenporo,100)dnl>
+define(lenoutlet,100)dnl>
+dnl>
+define(xhyp,calc(Sin*width)) dnl>
+define(yhyp,calc(Cos*width)) dnl>
+define(xinlet,leninlet)dnl>
+define(xporo,calc(Cos*lenporo)) dnl>
+define(yporo,calc(Sin*lenporo)) dnl>
+define(xoutlet,calc(xporo + Cos*lenoutlet)) dnl>
+define(youtlet,calc(yporo + Sin*lenoutlet)) dnl>
+dnl>
+
+convertToMeters 0.001;
+
+vertices
+(
+    // inlet region
+    ( -xinlet  y0  zBack )  vlabel(in1b)
+    ( -xinlet yhyp  zBack ) vlabel(in2b)
+    ( -xinlet  y0  zFront )  vlabel(in1f)
+    ( -xinlet yhyp  zFront ) vlabel(in2f)
+
+    // join inlet->outlet
+    (  x0 y0  zBack )    vlabel(join1b)
+    ( -xhyp   yhyp  zBack ) vlabel(join2b)
+    (  x0 y0  zFront )    vlabel(join1f)
+    ( -xhyp   yhyp  zFront ) vlabel(join2f)
+
+    // porosity ends ->outlet
+    ( xporo yporo  zBack )  vlabel(poro1b)
+    ( calc(xporo - xhyp) calc(yporo + yhyp)  zBack )  vlabel(poro2b)
+    ( xporo yporo  zFront )  vlabel(poro1f)
+    ( calc(xporo - xhyp) calc(yporo + yhyp)  zFront )  vlabel(poro2f)
+
+    // outlet
+    ( xoutlet youtlet zBack ) vlabel(out1b)
+    ( calc(xoutlet - xhyp) calc(youtlet + yhyp) zBack ) vlabel(out2b)
+    ( xoutlet youtlet zFront ) vlabel(out1f)
+    ( calc(xoutlet - xhyp) calc(youtlet + yhyp) zFront ) vlabel(out2f)
+);
+
+blocks
+(
+    // inlet block
+    hex2D(in1, join1, join2, in2)
+    inlet ( ninlet ncells ncells ) simpleGrading (1 1 1)
+
+    // porosity block
+    hex2D(join1, poro1, poro2, join2)
+    porosity ( nporo ncells ncells ) simpleGrading (1 1 1)
+
+    // outlet block
+    hex2D(poro1, out1, out2, poro2)
+    outlet ( noutlet ncells ncells )  simpleGrading (1 1 1)
+);
+
+edges
+(
+);
+
+patches
+(
+    // is there no way of defining all my 'defaultFaces' to be 'wall'?
+    wall front
+    (
+    // inlet block
+    frontQuad(in1, join1, join2, in2)
+    // outlet block
+    frontQuad(poro1, out1, out2, poro2)
+    )
+
+    wall back
+    (
+    // inlet block
+    backQuad(in1, join1, join2, in2)
+    // outlet block
+    backQuad(poro1, out1, out2, poro2)
+    )
+
+    wall wall
+    (
+    // inlet block
+    quad2D(in1, join1)
+    quad2D(join2, in2)
+    // outlet block
+    quad2D(poro1, out1)
+    quad2D(out2, poro2)
+    )
+
+    wall porosityWall
+    (
+    // porosity block
+    frontQuad(join1, poro1, poro2, join2)
+    // porosity block
+    backQuad(join1, poro1, poro2, join2)
+    // porosity block
+    quad2D(join1, poro1)
+    quad2D(poro2, join2)
+    )
+
+    patch inlet
+    (
+    quad2D(in2, in1)
+    )
+
+    patch outlet
+    (
+    quad2D(out2, out1)
+    )
+);
+
+mergePatchPairs
+(
+);
+
+// ************************************************************************* //

tutorials/compressible/rhoLTSPimpleFoam/angledDuct/constant/polyMesh/boundary

@@ -0,0 +1,58 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       polyBoundaryMesh;
+    location    "constant/polyMesh";
+    object      boundary;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+6
+(
+    front
+    {
+        type            wall;
+        nFaces          700;
+        startFace       63400;
+    }
+    back
+    {
+        type            wall;
+        nFaces          700;
+        startFace       64100;
+    }
+    wall
+    {
+        type            wall;
+        nFaces          1400;
+        startFace       64800;
+    }
+    porosityWall
+    {
+        type            wall;
+        nFaces          1600;
+        startFace       66200;
+    }
+    inlet
+    {
+        type            patch;
+        nFaces          400;
+        startFace       67800;
+    }
+    outlet
+    {
+        type            patch;
+        nFaces          400;
+        startFace       68200;
+    }
+)
+
+// ************************************************************************* //

tutorials/compressible/rhoLTSPimpleFoam/angledDuct/constant/porousZones

@@ -0,0 +1,36 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      porousZones;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+1
+(
+    porosity
+    {
+        coordinateSystem
+        {
+            e1  (0.70710678 0.70710678 0);
+            e2  (0 0 1);
+        }
+
+        Darcy
+        {
+            d   d [0 -2 0 0 0 0 0] (5e7 -1000 -1000);
+            f   f [0 -1 0 0 0 0 0] (0 0 0);
+        }
+    }
+)
+
+// ************************************************************************* //

tutorials/compressible/rhoLTSPimpleFoam/angledDuct/constant/thermophysicalProperties

@@ -0,0 +1,31 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      thermophysicalProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+thermoType      hPsiThermo<pureMixture<sutherlandTransport<specieThermo<hConstThermo<perfectGas>>>>>;
+
+mixture
+{
+    nMoles      1;
+    molWeight   28.9;
+    Cp          1007;
+    Hf          0;
+    As          1.4792e-06;
+    Ts          116;
+}
+
+
+// ************************************************************************* //

tutorials/compressible/rhoLTSPimpleFoam/angledDuct/constant/turbulenceProperties

@@ -0,0 +1,21 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      turbulenceProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+simulationType  RASModel;
+
+
+// ************************************************************************* //

tutorials/compressible/rhoLTSPimpleFoam/angledDuct/system/controlDict

@@ -0,0 +1,49 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+application     rhoLTSPimpleFoam;
+
+startFrom       startTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         500;
+
+deltaT          1;
+
+writeControl    adjustableRunTime;
+
+writeInterval   50;
+
+purgeWrite      0;
+
+writeFormat     ascii;
+
+writePrecision  6;
+
+writeCompression off;
+
+timeFormat      general;
+
+timePrecision   6;
+
+runTimeModifiable true;
+
+
+// ************************************************************************* //

tutorials/compressible/rhoLTSPimpleFoam/angledDuct/system/fvSchemes

@@ -0,0 +1,76 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+    default         localEuler rDeltaT;
+}
+
+gradSchemes
+{
+    default         Gauss linear;
+    grad(p)         Gauss linear;
+}
+
+divSchemes
+{
+    default         none;
+    div(phi,U)      Gauss upwind;
+    div(phid,p)     Gauss upwind;
+    div(phiU,p)     Gauss linear;
+    div(phi,h)      Gauss upwind;
+    div(phi,k)      Gauss upwind;
+    div(phi,epsilon) Gauss upwind;
+    div(phi,R)      Gauss upwind;
+    div(phi,omega)  Gauss upwind;
+    div((rho*R))    Gauss linear;
+    div(R)          Gauss linear;
+    div(U)          Gauss linear;
+    div((muEff*dev2(T(grad(U))))) Gauss linear;
+}
+
+laplacianSchemes
+{
+    default         none;
+    laplacian(muEff,U) Gauss linear corrected;
+    laplacian(mut,U) Gauss linear corrected;
+    laplacian(DkEff,k) Gauss linear corrected;
+    laplacian(DepsilonEff,epsilon) Gauss linear corrected;
+    laplacian(DREff,R) Gauss linear corrected;
+    laplacian(DomegaEff,omega) Gauss linear corrected;
+    laplacian((rho*(1|A(U))),p) Gauss linear corrected;
+    laplacian(alphaEff,h) Gauss linear corrected;
+}
+
+interpolationSchemes
+{
+    default         linear;
+}
+
+snGradSchemes
+{
+    default         corrected;
+}
+
+fluxRequired
+{
+    default         no;
+    p               ;
+}
+
+
+// ************************************************************************* //

tutorials/compressible/rhoLTSPimpleFoam/angledDuct/system/fvSolution

@@ -0,0 +1,68 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+    p
+    {
+        solver          PCG;
+        preconditioner  DIC;
+        tolerance       1e-06;
+        relTol          0.1;
+    }
+
+    pFinal
+    {
+        $p;
+        tolerance       1e-06;
+        relTol          0.1;
+    }
+
+    "(rho|U|h|k|epsilon|omega)"
+    {
+        solver          PBiCG;
+        preconditioner  DILU;
+        tolerance       1e-05;
+        relTol          0.1;
+    }
+
+    "(rho|U|h|k|epsilon|omega)Final"
+    {
+        $U;
+        tolerance       1e-05;
+        relTol          0.1;
+    }
+
+}
+
+PIMPLE
+{
+    momentumPredictor yes;
+    nOuterCorrectors  1;
+    nCorrectors       1;
+    nNonOrthogonalCorrectors 0;
+    rhoMin            rhoMin [ 1 -3 0 0 0 ] 0.5;
+    rhoMax            rhoMax [ 1 -3 0 0 0 ] 2.0;
+
+    maxCo             0.8;
+    rDeltaTSmoothingCoeff 0.02;
+    rDeltaTDampingCoeff 1;
+    maxDeltaT         1;
+}
+
+
+// ************************************************************************* //

tutorials/lagrangian/LTSReactingParcelFoam/counterFlowFlame2D/system/controlDict

@@ -15,7 +15,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-application     reactingFoam;
+application     LTSReactingParcelFoam;
 
 startFrom       latestTime;
 
@@ -51,8 +51,6 @@ maxCo           0.5;
 
 maxDeltaT       1;
 
-alphaTauSpecie  1e10;
-
 alphaTauTemp    0.005;
 
 // ************************************************************************* //

tutorials/lagrangian/LTSReactingParcelFoam/counterFlowFlame2D/system/fvSchemes

@@ -17,7 +17,7 @@ FoamFile
 
 ddtSchemes
 {
-    default         localEuler invTauFlow;
+    default         localEuler invTau;
 }
 
 gradSchemes

tutorials/lagrangian/LTSReactingParcelFoam/counterFlowFlame2D/system/fvSolution

@@ -24,6 +24,13 @@ solvers
         tolerance       0;
         relTol          0.1;
     }
+    rhoFinal
+    {
+        solver          PCG;
+        preconditioner  DIC;
+        tolerance       1e-4;
+        relTol          0;
+    }
     p
     {
         solver           GAMG;
@@ -44,6 +51,14 @@ solvers
         maxIter          50;
     };
 
+    pFinal
+    {
+        $p;
+
+        tolerance        1e-4;
+        relTol           0;
+    };
+
     "(U|Yi|hs|k|epsilon)"
     {
         solver          smoothSolver;
@@ -51,9 +66,17 @@ solvers
         tolerance       0;
         relTol          0.1;
     }
+
+    "(U|Yi|hs|k|epsilon)Final"
+    {
+        solver          smoothSolver;
+        smoother        GaussSeidel;
+        tolerance       1e-4;
+        relTol          0;
+    }
 }
 
-PISO
+PIMPLE
 {
     nCorrectors         2;
     nNonOrthogonalCorrectors 0;
@@ -64,8 +87,8 @@ PISO
 
 additional
 {
-    eWork           true;
+    pressureWork    true;
-    hWork           true;
+    pressureWorkTimeDerivative true;
     solveSpecies    true;
 }
 

tutorials/lagrangian/LTSReactingParcelFoam/verticalChannel/constant/reactingCloud1Properties

@@ -122,7 +122,7 @@ subModels
         massFlowRate    0.8e-03;
         parcelBasisType mass;
         patchName       inletCentral;
-        parcelsPerSecond 400;
+        parcelsPerSecond 100;
         duration        1; // NOTE: set to 1 for steady state
         U0              (0 40 0);
         flowRateProfile constant 1;

tutorials/lagrangian/LTSReactingParcelFoam/verticalChannel/system/controlDict

@@ -45,17 +45,13 @@ timePrecision   6;
 
 runTimeModifiable yes;
 
-maxCo               5;
+maxCo           5;
 
-alphaTauSmooth      1e10;
+alphaTauSmooth  1;
 
-alphaTauTemp        1e10;
+alphaTauTemp    1;
 
-alphaTauRho         1e10;
+maxDeltaT       1;
-
-alphaTauSpecie      1e10;
-
-maxDeltaT           1e10;
 
 functions
 {

tutorials/lagrangian/LTSReactingParcelFoam/verticalChannel/system/fvSchemes

@@ -17,7 +17,7 @@ FoamFile
 
 ddtSchemes
 {
-    default         localEuler invTauFlow;
+    default         localEuler invTau;
 }
 
 gradSchemes

tutorials/lagrangian/LTSReactingParcelFoam/verticalChannel/system/fvSolution

@@ -1,14 +1,14 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
 | \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  dev                                   |
+|  \\    /   O peration     | Version:  1.6                                   |
 |   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/
 FoamFile
 {
     version     2.0;
-    format      ascii;
+    format      binary;
     class       dictionary;
     location    "system";
     object      fvSolution;
@@ -17,11 +17,21 @@ FoamFile
 
 solvers
 {
+    "(rho|.*Final)"
+    {
+        solver          PCG;
+        preconditioner  DIC;
+        tolerance       0;
+        relTol          0.1;
+    }
+
+//    "(p|.*Final)"
     p
     {
         solver           GAMG;
+
         tolerance        0;
-        relTol           0.1;
+        relTol           0.05;
 
         smoother         DICGaussSeidel;
         nPreSweeps       0;
@@ -33,47 +43,39 @@ solvers
         agglomerator     faceAreaPair;
         mergeLevels      1;
 
-        maxIter          20;
+        maxIter          50;
     };
-    "(rho|G)"
+
+    pFinal
     {
-        solver           PCG;
+        $p;
-        preconditioner   DIC;
+        tolerance        1e-4;
-        tolerance        0;
+    }
-        relTol           0.1;
+
-        maxIter          20;
+    "(U|Yi|hs|k|omega|.*Final)"
-    };
-    "(Yi|hs)"
     {
-        solver           smoothSolver;
+        solver          smoothSolver;
-        smoother         GaussSeidel;
+        smoother        GaussSeidel;
-        tolerance        0;
+        tolerance       0;
-        relTol           0.1;
+        relTol          0.1;
-        maxIter          20;
+        maxIter         20;
-    };
+    }
-    "(U|k|omega)"
-    {
-        solver           smoothSolver;
-        smoother         GaussSeidel;
-        tolerance        0;
-        relTol           0.1;
-        maxIter          20;
-    };
 }
 
-PISO
+PIMPLE
 {
-    nCorrectors          2;
+    nOuterCorrectors    1;
-    nNonOrthogonalCorrectors 0; // 10;
+    nCorrectors         2;
+    nNonOrthogonalCorrectors 0;
     momentumPredictor yes;
-    rhoMin rhoMin [ 1 -3 0 0 0 ] 0; // 0.5;
+    rhoMin rhoMin [1 -3 0 0 0] 0.1;
-    rhoMax rhoMax [ 1 -3 0 0 0 ] 100; // 1.5;
+    rhoMax rhoMax [1 -3 0 0 0] 1.5;
 }
 
 additional
 {
-    eWork           true;
+    pressureWork    true;
-    hWork           true;
+    pressureWorkTimeDerivative false; // true;
     solveSpecies    true;
 }
 

tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/filter/system/fvSolution

@@ -74,9 +74,8 @@ PISO
 
 additional
 {
-    dpdt            true;
+    pressureWork    true;
-    hWork           true;
+    pressureWorkTimeDerivative true;
-    eWork           true;
 }
 
 relaxationFactors

tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/parcelInBox/system/fvSolution

@@ -74,9 +74,8 @@ PISO
 
 additional
 {
-    dpdt            false;
+    pressureWork    true;
-    eWork           true;
+    pressureWorkTimeDerivative true;
-    hWork           true;
 }
 
 

tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/system/fvSolution

@@ -82,9 +82,8 @@ SIMPLE
 
 additional
 {
-    dpdt            false;
+    pressureWork    false;
-    eWork           true;
+    pressureWorkTimeDerivative false;
-    hWork           true;
 }
 
 

tutorials/lagrangian/reactingParcelFilmFoam/evaporationTest/system/controlDict

@@ -15,7 +15,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-application     reactingParcelFoam;
+application     reactingParcelFilmFoam;
 
 startFrom       latestTime;
 

tutorials/lagrangian/reactingParcelFilmFoam/hotBoxes/system/controlDict

@@ -15,7 +15,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-application     reactingParcelFoam;
+application     reactingParcelFilmFoam;
 
 startFrom       latestTime; // startTime;
 

tutorials/lagrangian/reactingParcelFilmFoam/multipleBoxes/system/controlDict

@@ -15,7 +15,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-application     reactingParcelFoam;
+application     reactingParcelFilmFoam;
 
 startFrom       startTime;
 

tutorials/lagrangian/reactingParcelFilmFoam/panel/system/controlDict

@@ -15,7 +15,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-application     reactingParcelFoam;
+application     reactingParcelFilmFoam;
 
 startFrom       startTime;