Merge remote branch 'opencfd/master' into olesenm

ReleaseNotes-dev

@@ -198,7 +198,7 @@
     {
         type            codedFixedValue;
         value           uniform 0;
-        redirectType    fixedValue10;
+        redirectType    ramped;
 
         code
         #{
@@ -277,6 +277,26 @@
        triSurfaceMesh).
     + =nearWallFields=: constructs field with on selected patches interpolated
        internal field for further postprocessing.
+    + =coded=: uses the dynamic code compilation from =#codeStream=
+    to provide an in-line functionObject. E.g.
+    #+BEGIN_SRC c++
+    functions
+    {
+        pAverage
+        {
+            functionObjectLibs ("libutilityFunctionObjects.so");
+            type coded;
+            redirectType average;
+            code
+            #{
+                const volScalarField& p = mesh().lookupObject<volScalarField>("p");
+                Info<<"p avg:" << average(p) << endl;
+            #};
+        }
+    }
+    #+END_SRC
+    See also [[./doc/changes/dynamicCode.org]]
+
 
 * New tutorials
   There is a large number of new tutorials for existing and new solvers in the

applications/solvers/combustion/fireFoam/fireFoam.C

@@ -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,7 +45,7 @@
                     p.select
                     (
                         (
-                            finalIter
+                            pimpleCorr.finalIter()
                          && corr == nCorr-1
                          && nonOrth == nNonOrthCorr
                         )
@@ -89,7 +89,7 @@
                     p.select
                     (
                         (
-                            finalIter
+                            pimpleCorr.finalIter()
                          && corr == nCorr-1
                          && nonOrth == nNonOrthCorr
                         )

applications/solvers/combustion/rhoReactingFoam/rhoReactingFoam.C

@@ -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,13 +37,11 @@ if (transonic)
             (
                 p.select
                 (
-                    (
+                    pimpleCorr.finalIter()
-                        finalIter
                  && corr == nCorr-1
                  && nonOrth == nNonOrthCorr
                 )
             )
-            )
         );
 
         if (nonOrth == nNonOrthCorr)
@@ -77,13 +75,11 @@ else
             (
                 p.select
                 (
-                    (
+                    pimpleCorr.finalIter()
-                        finalIter
                  && corr == nCorr-1
                  && nonOrth == nNonOrthCorr
                 )
             )
-            )
         );
 
         if (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

@@ -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

@@ -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,13 +27,11 @@ for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
         (
             p.select
             (
-                (
+                pimpleCorr.finalIter()
-                    finalIter
              && corr == nCorr-1
              && nonOrth == nNonOrthCorr
             )
         )
-        )
     );
 
     if (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

@@ -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

@@ -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,6 +87,17 @@ int main(int argc, char *argv[])
         #include "timeScales.H"
 
         #include "rhoEqn.H"
+
+        // --- Pressure-velocity PIMPLE corrector loop
+        for
+        (
+            pimpleLoop pimpleCorr(mesh, nOuterCorr);
+            pimpleCorr.loop();
+            pimpleCorr++
+        )
+        {
+            turbulence->correct();
+
             #include "UEqn.H"
             #include "YEqn.H"
             #include "hsEqn.H"
@@ -95,8 +107,7 @@ int main(int argc, char *argv[])
             {
                 #include "pEqn.H"
             }
-
+        }
-        turbulence->correct();
 
         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)
         {

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

@@ -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,13 +32,11 @@ if (transonic)
             (
                 p.select
                 (
-                    (
+                    pimpleCorr.finalIter()
-                        finalIter
                  && corr == nCorr-1
                  && nonOrth == nNonOrthCorr
                 )
             )
-            )
         );
 
         if (nonOrth == nNonOrthCorr)
@@ -73,13 +71,11 @@ else
             (
                 p.select
                 (
-                    (
+                    pimpleCorr.finalIter()
-                        finalIter
                  && corr == nCorr-1
                  && nonOrth == nNonOrthCorr
                 )
             )
-            )
         );
 
         if (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)
+    {
+        surfaceScalarField phiU("phiU", phi/fvc::interpolate(rho));
+
+        pWork() += fvc::div(phiU*fvc::interpolate(p)) - p*fvc::div(phiU);
+
+        if (pressureWorkTimeDerivative)
         {
             pWork() += fvc::ddt(p);
         }
-    if (eWork)
-    {
-        pWork() = -p*fvc::div(phi/fvc::interpolate(rho));
-    }
-    if (hWork)
-    {
-        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,23 +78,28 @@ 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();
+        }
 
         if (runTime.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

@@ -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

@@ -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();

applications/test/PtrList/Test-PtrList.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
@@ -56,7 +56,10 @@ public:
         Info<<"delete Scalar: " << data_ << endl;
     }
 
-    autoPtr<Scalar> clone() const;
+    autoPtr<Scalar> clone() const
+    {
+        return autoPtr<Scalar>(new Scalar(data_));
+    }
 
     friend Ostream& operator<<(Ostream& os, const Scalar& val)
     {

applications/test/dataEntry/Make/options

@@ -8,4 +8,5 @@ EXE_LIBS = \
     -lfiniteVolume \
     -llagrangianIntermediate \
     -lradiationModels \
-    -lthermophysicalFunctions
+    -lthermophysicalFunctions \
+    -lregionModels

applications/test/liquid/Make/options

@@ -1,5 +1,5 @@
 EXE_INC = \
-    -I$(LIB_SRC)/thermophysicalModels/liquids/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/properties/liquidProperties/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/thermophysicalFunctions/lnInclude
 
-EXE_LIBS = -lliquids -lthermophysicalFunctions
+EXE_LIBS = -lliquidProperties -lthermophysicalFunctions

applications/test/passiveParticle/Test-passiveParticle.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,7 +55,7 @@ int main(int argc, char *argv[])
         Pout<< "Starting particles:" << particles.size() << endl;
 
         Pout<< "Adding a particle." << endl;
-        particles.addParticle(new passiveParticle(particles, vector::zero, -1));
+        particles.addParticle(new passiveParticle(mesh, vector::zero, -1));
 
         forAllConstIter(passiveParticleCloud, particles, iter)
         {

bin/tools/CleanFunctions

@@ -52,7 +52,7 @@ cleanTimeDirectories()
         rm -rf ./${timeDir} ./-${timeDir} > /dev/null 2>&1
         zeros="0$zeros"
     done
-    rm -rf ./[1-9]* ./-[1-9]* ./log ./log.* ./log-* ./logSummary.* ./.fxLock ./*.xml ./ParaView* ./paraFoam* ./*.OpenFOAM > /dev/null 2>&1
+    rm -rf ./[1-9]* ./-[1-9]* ./log ./log.* ./log-* ./logSummary.* ./.fxLock ./*.xml ./ParaView* ./paraFoam* ./*.OpenFOAM ./.setSet > /dev/null 2>&1
 }
 
 

doc/changes/dynamicCode.org

@@ -8,7 +8,7 @@
 # Copyright (c) 2011 OpenCFD Ltd.
 
 * Dictionary preprocessing directive: =#codeStream=
-  This is a dictionary preprocessing directive ('=functionEntry=') which
+  This is a dictionary preprocessing directive (=functionEntry=) which
   provides a snippet of OpenFOAM C++ code which gets compiled and executed to
   provide the actual dictionary entry. The snippet gets provided as three
   sections of C++ code which just gets inserted into a template:
@@ -17,8 +17,11 @@
     =dict.lookup= to find current dictionary values.
   - optional =codeInclude= section: any #include statements to include OpenFOAM
     files.
-  - optional 'codeOptions' section: any extra compilation flags to be added to
+  - optional =codeOptions= section: any extra compilation flags to be added to
-    =EXE_INC= in =Make/options=
+    =EXE_INC= in =Make/options=. These usually are =-I= include directory
+    options.
+  - optional =codeLibs= section: any extra compilation flags to be added to
+    =LIB_LIBS= in =Make/options=.
 
   To ease inputting mulit-line code there is the =#{ #}= syntax. Anything in
   between these two delimiters becomes a string with all newlines, quotes etc
@@ -26,6 +29,7 @@
 
   Example: Look up dictionary entries and do some calculation
   #+BEGIN_SRC c++
+
     startTime       0;
     endTime         100;
     ..
@@ -43,30 +47,30 @@
 
 * Implementation
   - the =#codeStream= entry reads the dictionary following it, extracts the
-    =code=, =codeInclude=, =codeOptions= sections (these are just strings) and
+    =code=, =codeInclude=, =codeOptions=, =codeLibs= sections (these are just strings) and
     calculates the SHA1 checksum of the contents.
   - it copies a template file
-    =(~OpenFOAM/codeTemplates/dynamicCode/codeStreamTemplate.C)= or
+    =(etc/codeTemplates/dynamicCode/codeStreamTemplate.C)= or
     =($FOAM_CODE_TEMPLATES/codeStreamTemplate.C)=, substituting all
-    occurences of =code=, =codeInclude=, =codeOptions=.
+    occurences of =code=, =codeInclude=, =codeOptions=, =codeLibs=.
   - it writes library source files to =dynamicCode/<SHA1>= and compiles
     it using =wmake libso=.
   - the resulting library is generated under
     =dynamicCode/platforms/$WM_OPTIONS/lib= and is loaded (=dlopen=, =dlsym=)
-    and the function executed
+    and the function executed.
   - the function will have written its output into the Ostream which then gets
     used to construct the entry to replace the whole =#codeStream= section.
   - using the SHA1 means that same code will only be compiled and loaded once.
 
 * Boundary condition: =codedFixedValue=
-  This uses the code from codeStream to have an in-line specialised
+  This uses the same framework as codeStream to have an in-line specialised
-  =fixedValueFvPatchScalarField=. For now only for scalars:
+  =fixedValueFvPatchField=.
   #+BEGIN_SRC c++
   outlet
   {
       type            codedFixedValue;
       value           uniform 0;
-      redirectType    fixedValue10;
+      redirectType    ramp;
 
       code
       #{
@@ -75,13 +79,16 @@
   }
   #+END_SRC
   It by default always includes =fvCFD.H= and adds the =finiteVolume= library to
-  the include search path.
+  the include search path and the linked libraries. Any other libraries will
+  need
+  to be added using the =codeInclude=, =codeLibs=, =codeOptions= section or provided through
+  the =libs= entry in the =system/controlDict=.
 
   A special form is where the code is not supplied in-line but instead comes
   from the =codeDict= dictionary in the =system= directory. It should contain
-  a =fixedValue10= entry:
+  a =ramp= entry:
   #+BEGIN_SRC c++
-  fixedValue10
+  ramp
   {
       code
       #{
@@ -90,7 +97,49 @@
   }
   #+END_SRC
   The advantage of using this indirect way is that it supports
-  runTimeModifiable so any change of the code will be picked up next iteration.
+  =runTimeModifiable= so any change of the code will be picked up next iteration.
+
+* Function object: =coded=
+  This uses the same framework as codeStream to have an in-line specialised
+  =functionObject=.
+  #+BEGIN_SRC c++
+  functions
+  {
+      pAverage
+      {
+          functionObjectLibs ("libutilityFunctionObjects.so");
+          type coded;
+          redirectType average;
+          outputControl outputTime;
+          code
+          #{
+              const volScalarField& p = mesh().lookupObject<volScalarField>("p");
+              Info<<"p avg:" << average(p) << endl;
+          #};
+      }
+  }
+  #+END_SRC
+  This dynamic code framework uses the following entries
+  + =codeData=: declaration (in .H file) of local (null-constructable) data
+  + =codeInclude=: (.C file) usual include section
+  + =codeRead=: (.C file) executed upon dictionary read
+  + =codeExecute=: (.C file) executed upon functionObject execute
+  + =codeEnd=: (.C file) executed upon functionObject end
+  + =code=: (.C file) executed upon functionObject write. This is the usual place
+  for simple functionObject.
+  + =codeLibs=, =codeOptions=: usual
+
+  =coded= by default always includes =fvCFD.H= and adds the =finiteVolume= library to
+  the include search path and the linked libraries. Any other libraries will
+  need to be added explicitly (see =codeInclude=, =codeLibs=, =codeOptions= sections) or provided through
+  the =libs= entry in the =system/controlDict=.
+
+  =coded= is an =OutputFilter= type =functionObject= so supports the usual
+  + =region=: non-default region
+  + =enabled=: enable/disable
+  + =outputControl=: =timeStep= or =outputTime=
+  + =outputInterval=: in case of =timeStep=
+  entries.
 
 * Security
   Allowing the case to execute C++ code does introduce security risks.  A
@@ -128,6 +177,11 @@
       #{
           -I$(LIB_SRC)/finiteVolume/lnInclude
       #};
+
+      codeLibs
+      #{
+          -lfiniteVolume
+      #};
   };
   #+END_SRC
 
@@ -144,7 +198,7 @@
     prints 'uniform 12.34;'. Note the ';' at the end. It is advised to use the
     =writeEntry= as above to handle this and also e.g. binary streams (=codeStream=
     inherits the stream type from the dictionary)
-  + the =code=, =codeInclude=, =codeOptions= entries are just like any other
+  + the =code=, =codeInclude=, =codeOptions=, =codeLibs= entries are just like any other
     dictionary string entry so there has to be a ';' after the string
   + the =#codeStream= entry (itself a dictionary) has to end in a ';'
 
@@ -153,14 +207,37 @@
   Following applications read
   the field as a dictionary, not as an =IOdictionary=:
   - =foamFormatConvert=
-  - =changeDictionaryDict=
+  - =changeDictionary=
   - =foamUpgradeCyclics=
   These applications will usually switch off all '#' processing which
-  just preserves the entries as strings (including all formatting).
+  just preserves the entries as strings (including all
+  formatting). =changeDictionary= has the =-enableFunctionEntries= option for if
+  one does want to evaluate any preprocessing in the changeDictionaryDict.
 
 * Other
+  - paraFoam: paraview currently does not export symbols on loaded libraries
+    (more specific : it does not add 'RTLD_GLOBAL' to the dlopen flags) so
+    one will have to add the used additional libraries (libfiniteVolume,
+    lib..) either to the =codeLibs= linkage section (preferred) or to the 'libs' entry in system/controlDict to prevent getting
+    an error of the form
+
+        --> FOAM FATAL IO ERROR:
+        Failed loading library "libcodeStream_3cd388ceb070a2f8b0ae61782adbc21c5687ce6f.so"
+
+    By default =#codeStream= links in =libOpenFOAM= and =codedFixedValue= and =coded=
+    functionObject link in both =libOpenFOAM= and =libfiniteVolume=.
+
   - parallel running not tested a lot. What about distributed data
     (i.e. non-=NFS=) parallel?
-  - paraview has been patched so it will pass in RTLD_GLOBAL when loading
+
-    the OpenFOAM reader module. This is necessary for above dictionary
+  - codedFixedValue could be extended to provide local data however
-    processing to work.
+    in terms of complexity this is not really worthwhile.
+
+  - all templates come from (in order of preference)
+    =FOAM_TEMPLATE_DIR=
+    =~/.OpenFOAM/dev/codeTemplates/dynamicCode=
+    =etc/codeTemplates/dynamicCode=
+
+  - any generated C++ code will display line numbers relative to the original
+    dictionary (using the '#line' directive) to ease finding compilation
+    errors.

etc/codeTemplates/dynamicCode/FilterFunctionObjectTemplate.C

@@ -0,0 +1,65 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "FilterFunctionObjectTemplate.H"
+
+// * * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * //
+
+extern "C"
+{
+    // dynamicCode:
+    // SHA1 = ${SHA1sum}
+    //
+    // unique function name that can be checked if the correct library version
+    // has been loaded
+    void ${typeName}_${SHA1sum}(bool load)
+    {
+        if (load)
+        {
+            // code that can be explicitly executed after loading
+        }
+        else
+        {
+            // code that can be explicitly executed before unloading
+        }
+    }
+}
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineNamedTemplateTypeNameAndDebug(${typeName}FilterFunctionObject, 0);
+
+    //addToRunTimeSelectionTable
+    addRemovableToRunTimeSelectionTable
+    (
+        functionObject,
+        ${typeName}FilterFunctionObject,
+        dictionary
+    );
+}
+
+// ************************************************************************* //

etc/codeTemplates/dynamicCode/FilterFunctionObjectTemplate.H

@@ -0,0 +1,51 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+Description
+    FunctionObject wrapper around functionObjectTemplate to allow them
+    to be created via the functions entry within controlDict.
+
+SourceFiles
+    FilterFunctionObject.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef FilterFunctionObject_H
+#define FilterFunctionObject_H
+
+#include "functionObjectTemplate.H"
+#include "OutputFilterFunctionObject.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    typedef OutputFilterFunctionObject<${typeName}FunctionObject>
+        ${typeName}FilterFunctionObject;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

etc/codeTemplates/dynamicCode/IOfunctionObjectTemplate.H

@@ -0,0 +1,50 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+Typedef
+    Foam::IOfunctionObjectTemplate
+
+Description
+    Instance of the generic IOOutputFilter for ${typeName}FunctionObject.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef IOfunctionObjectTemplate_H
+#define IOfunctionObjectTemplate_H
+
+#include "functionObjectTemplate.H"
+#include "IOOutputFilter.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    typedef IOOutputFilter<${typeName}FunctionObject>
+    IO${typeName}FunctionObject;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

etc/codeTemplates/dynamicCode/functionObjectTemplate.C

@@ -0,0 +1,146 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "functionObjectTemplate.H"
+#include "Time.H"
+#include "fvCFD.H"
+
+//{{{ begin codeInclude
+${codeInclude}
+//}}} end codeInclude
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(${typeName}FunctionObject, 0);
+
+
+// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
+
+//{{{ begin localCode
+${localCode}
+//}}} end localCode
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+const objectRegistry& ${typeName}FunctionObject::obr() const
+{
+    return obr_;
+}
+
+
+const fvMesh& ${typeName}FunctionObject::mesh() const
+{
+    return refCast<const fvMesh>(obr_);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+${typeName}FunctionObject::${typeName}FunctionObject
+(
+    const word& name,
+    const objectRegistry& obr,
+    const dictionary& dict,
+    const bool
+)
+:
+    name_(name),
+    obr_(obr)
+{
+    read(dict);
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+${typeName}FunctionObject::~${typeName}FunctionObject()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void ${typeName}FunctionObject::read(const dictionary& dict)
+{
+    if (${verbose:-false})
+    {
+        Info<<"read ${typeName} sha1: ${SHA1sum}\n";
+    }
+
+//{{{ begin code
+    ${codeRead}
+//}}} end code
+}
+
+
+void ${typeName}FunctionObject::execute()
+{
+    if (${verbose:-false})
+    {
+        Info<<"execute ${typeName} sha1: ${SHA1sum}\n";
+    }
+
+//{{{ begin code
+    ${codeExecute}
+//}}} end code
+}
+
+
+void ${typeName}FunctionObject::end()
+{
+    if (${verbose:-false})
+    {
+        Info<<"end ${typeName} sha1: ${SHA1sum}\n";
+    }
+
+//{{{ begin code
+    ${codeEnd}
+//}}} end code
+}
+
+
+void ${typeName}FunctionObject::write()
+{
+    if (${verbose:-false})
+    {
+        Info<<"write ${typeName} sha1: ${SHA1sum}\n";
+    }
+
+//{{{ begin code
+    ${code}
+//}}} end code
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+
+// ************************************************************************* //

etc/codeTemplates/dynamicCode/functionObjectTemplate.H

@@ -0,0 +1,143 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+Description
+    Template for use with dynamic code generation of a
+    OutputFilter functionObject.
+
+SourceFiles
+    functionObjectTemplate.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef functionObjectTemplate_H
+#define functionObjectTemplate_H
+
+#include "stringList.H"
+#include "pointField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+class objectRegistry;
+class dictionary;
+class mapPolyMesh;
+class fvMesh;
+
+/*---------------------------------------------------------------------------*\
+                         A templated functionObject
+\*---------------------------------------------------------------------------*/
+
+class ${typeName}FunctionObject
+{
+    // Private data
+
+        //- Name of this set of system calls
+        word name_;
+
+        //- Registry
+        const objectRegistry& obr_;
+
+//{{{ begin codeData
+        ${codeData}
+//}}} end codeData
+
+
+    // Private Member Functions
+
+        const objectRegistry& obr() const;
+
+        const fvMesh& mesh() const;
+
+        //- Disallow default bitwise copy construct
+        ${typeName}FunctionObject(const ${typeName}FunctionObject&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const ${typeName}FunctionObject&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("${typeName}");
+
+
+    // Constructors
+
+        //- Construct for given objectRegistry and dictionary.
+        //  Allow the possibility to load fields from files
+        ${typeName}FunctionObject
+        (
+            const word& name,
+            const objectRegistry& unused,
+            const dictionary&,
+            const bool loadFromFilesUnused = false
+        );
+
+
+    //- Destructor
+    virtual ~${typeName}FunctionObject();
+
+
+    // Member Functions
+
+        //- Return name of the system call set
+        virtual const word& name() const
+        {
+            return name_;
+        }
+
+        //- Read the system calls
+        virtual void read(const dictionary&);
+
+        //- Execute the "executeCalls" at each time-step
+        virtual void execute();
+
+        //- Execute the "endCalls" at the final time-loop
+        virtual void end();
+
+        //- Write, execute the "writeCalls"
+        virtual void write();
+
+        //- Update for changes of mesh
+        virtual void updateMesh(const mapPolyMesh&)
+        {}
+
+        //- Update for changes of mesh
+        virtual void movePoints(const pointField&)
+        {}
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/OSspecific/POSIX/POSIX.C

@@ -37,6 +37,7 @@ Description
 #include "fileStat.H"
 #include "timer.H"
 #include "IFstream.H"
+#include "DynamicList.H"
 
 #include <fstream>
 #include <cstdlib>
@@ -52,6 +53,7 @@ Description
 #include <sys/socket.h>
 #include <netdb.h>
 #include <dlfcn.h>
+#include <link.h>
 
 #include <netinet/in.h>
 
@@ -1107,10 +1109,20 @@ void* Foam::dlOpen(const fileName& lib)
 {
     if (POSIX::debug)
     {
-        Info<< "dlOpen(const fileName&)"
+        std::cout<< "dlOpen(const fileName&)"
-            << " : dlopen of " << lib << endl;
+            << " : dlopen of " << lib << std::endl;
     }
-    return ::dlopen(lib.c_str(), RTLD_LAZY|RTLD_GLOBAL);
+    void* handle = ::dlopen(lib.c_str(), RTLD_LAZY|RTLD_GLOBAL);
+
+    if (POSIX::debug)
+    {
+        std::cout
+            << "dlOpen(const fileName&)"
+            << " : dlopen of " << lib
+            << " handle " << handle << std::endl;
+    }
+
+    return handle;
 }
 
 
@@ -1118,8 +1130,9 @@ bool Foam::dlClose(void* handle)
 {
     if (POSIX::debug)
     {
-        Info<< "dlClose(void*)"
+        std::cout
-            << " : dlclose" << endl;
+            << "dlClose(void*)"
+            << " : dlclose of handle " << handle << std::endl;
     }
     return ::dlclose(handle) == 0;
 }
@@ -1129,8 +1142,9 @@ void* Foam::dlSym(void* handle, const std::string& symbol)
 {
     if (POSIX::debug)
     {
-        Info<< "dlSym(void*, const std::string&)"
+        std::cout
-            << " : dlsym of " << symbol << endl;
+            << "dlSym(void*, const std::string&)"
+            << " : dlsym of " << symbol << std::endl;
     }
     // clear any old errors - see manpage dlopen
     (void) ::dlerror();
@@ -1158,8 +1172,9 @@ bool Foam::dlSymFound(void* handle, const std::string& symbol)
     {
         if (POSIX::debug)
         {
-            Info<< "dlSymFound(void*, const std::string&)"
+            std::cout
-                << " : dlsym of " << symbol << endl;
+                << "dlSymFound(void*, const std::string&)"
+                << " : dlsym of " << symbol << std::endl;
         }
 
         // clear any old errors - see manpage dlopen
@@ -1178,4 +1193,32 @@ bool Foam::dlSymFound(void* handle, const std::string& symbol)
 }
 
 
+static int collectLibsCallback
+(
+    struct dl_phdr_info *info,
+    size_t size,
+    void *data
+)
+{
+    Foam::DynamicList<Foam::fileName>* ptr =
+        reinterpret_cast<Foam::DynamicList<Foam::fileName>*>(data);
+    ptr->append(info->dlpi_name);
+    return 0;
+}
+
+
+Foam::fileNameList Foam::dlLoaded()
+{
+    DynamicList<fileName> libs;
+    dl_iterate_phdr(collectLibsCallback, &libs);
+    if (POSIX::debug)
+    {
+        std::cout
+            << "dlLoaded()"
+            << " : determined loaded libraries :" << libs.size() << std::endl;
+    }
+    return libs;
+}
+
+
 // ************************************************************************* //

src/OpenFOAM/db/dictionary/functionEntries/codeStream/codeStream.C

@@ -136,7 +136,9 @@ bool Foam::functionEntries::codeStream::execute
                 (
                     "EXE_INC = -g \\\n"
                   + context.options()
-                  + "\n\nLIB_LIBS ="
+                  + "\n\nLIB_LIBS = \\\n"
+                  + "    -lOpenFOAM \\\n"
+                  + context.libs()
                 );
 
                 if (!dynCode.copyOrCreateFiles(true))
@@ -172,6 +174,8 @@ bool Foam::functionEntries::codeStream::execute
                 "functionEntries::codeStream::execute(..)",
                 parentDict
             )   << "Failed loading library " << libPath << nl
+                << "Did you add all libraries to the 'libs' entry"
+                << " in system/controlDict?"
                 << exit(FatalIOError);
         }
 

src/OpenFOAM/db/dynamicLibrary/dynamicCode/dynamicCode.C

@@ -31,8 +31,6 @@ License
 #include "OFstream.H"
 #include "OSspecific.H"
 #include "dictionary.H"
-#include "dlLibraryTable.H"
-
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -524,7 +522,7 @@ bool Foam::dynamicCode::copyOrCreateFiles(const bool verbose) const
 
 bool Foam::dynamicCode::wmakeLibso() const
 {
-    const Foam::string wmakeCmd("wmake -s libso " + this->codeRelPath());
+    const Foam::string wmakeCmd("wmake -s libso " + this->codePath());
     Info<< "Invoking " << wmakeCmd << endl;
 
     if (Foam::system(wmakeCmd))
@@ -557,22 +555,4 @@ bool Foam::dynamicCode::upToDate(const dynamicCodeContext& context) const
 }
 
 
-// bool Foam::dynamicCode::openLibrary() const
-// {
-//     return dlLibraryTable::openLibrary(this->libPath(), false);
-// }
-//
-//
-// bool Foam::dynamicCode::closeLibrary() const
-// {
-//     return dlLibraryTable::closeLibrary(this->libPath(), false);
-// }
-//
-//
-// void* Foam::dynamicCode::findLibrary() const
-// {
-//     return dlLibraryTable::findLibrary(this->libPath());
-// }
-
-
 // ************************************************************************* //

src/OpenFOAM/db/dynamicLibrary/dynamicCode/dynamicCode.H

@@ -37,10 +37,8 @@ SourceFiles
 #define dynamicCode_H
 
 #include "Tuple2.H"
-#include "SHA1Digest.H"
 #include "HashTable.H"
 #include "DynamicList.H"
-#include "dlLibraryTable.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -51,6 +49,7 @@ namespace Foam
 class dynamicCodeContext;
 class ISstream;
 class OSstream;
+class SHA1Digest;
 
 /*---------------------------------------------------------------------------*\
                         Class dynamicCode Declaration
@@ -283,16 +282,6 @@ public:
         //- Compile a libso
         bool wmakeLibso() const;
 
-//        //- Open the libPath() library
-//        bool openLibrary() const;
-//
-//        //- Close the libPath() library
-//        bool closeLibrary() const;
-//
-//        //- Find the handle of the libPath() library
-//        void* findLibrary() const;
-
-
 };
 
 

src/OpenFOAM/db/dynamicLibrary/dynamicCode/dynamicCodeContext.C

@@ -34,44 +34,78 @@ License
 Foam::dynamicCodeContext::dynamicCodeContext(const dictionary& dict)
 :
     dict_(dict),
-    code_(stringOps::trim(dict["code"])),
+    code_(),
     localCode_(),
     include_(),
-    options_()
+    options_(),
+    libs_()
 {
     // expand dictionary entries
+
+    {
+        const entry& codeEntry = dict.lookupEntry("code", false, false);
+        code_ = stringOps::trim(codeEntry.stream());
         stringOps::inplaceExpand(code_, dict);
+        addLineDirective(code_, codeEntry.startLineNumber(), dict.name());
+    }
 
     // note: removes any leading/trailing whitespace
     // - necessary for compilation options, convenient for includes
     // and body.
 
     // optional
-    if (dict.found("localCode"))
+    const entry* includePtr = dict.lookupEntryPtr
+    (
+        "codeInclude",
+        false,
+        false
+    );
+    if (includePtr)
     {
-        localCode_ = stringOps::trim(dict["localCode"]);
+        include_ = stringOps::trim(includePtr->stream());
-        stringOps::inplaceExpand(localCode_, dict);
-    }
-
-    // optional
-    if (dict.found("codeInclude"))
-    {
-        include_ = stringOps::trim(dict["codeInclude"]);
         stringOps::inplaceExpand(include_, dict);
+        addLineDirective(include_, includePtr->startLineNumber(), dict.name());
     }
 
     // optional
-    if (dict.found("codeOptions"))
+    const entry* optionsPtr = dict.lookupEntryPtr
+    (
+        "codeOptions",
+        false,
+        false
+    );
+    if (optionsPtr)
     {
-        options_ = stringOps::trim(dict["codeOptions"]);
+        options_ = stringOps::trim(optionsPtr->stream());
         stringOps::inplaceExpand(options_, dict);
     }
 
+    // optional
+    const entry* libsPtr = dict.lookupEntryPtr("codeLibs", false, false);
+    if (libsPtr)
+    {
+        libs_ = stringOps::trim(libsPtr->stream());
+        stringOps::inplaceExpand(libs_, dict);
+    }
+
     // calculate SHA1 digest from include, options, localCode, code
     OSHA1stream os;
-    os  << include_ << options_ << localCode_ << code_;
+    os  << include_ << options_ << libs_ << localCode_ << code_;
     sha1_ = os.digest();
 }
 
 
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+void Foam::dynamicCodeContext::addLineDirective
+(
+    string& code,
+    const label lineNum,
+    const fileName& name
+)
+{
+    code = "#line " + Foam::name(lineNum) + " \"" + name + "\"\n" + code;
+}
+
+
 // ************************************************************************* //

src/OpenFOAM/db/dynamicLibrary/dynamicCode/dynamicCodeContext.H

@@ -51,6 +51,7 @@ namespace Foam
 class dynamicCodeContext
 {
     // Private data
+
         //- The parent dictionary context
         const dictionary& dict_;
 
@@ -66,6 +67,9 @@ class dynamicCodeContext
         //- Optional "codeOptions" entry
         string options_;
 
+        //- Optional "codeLib" entry
+        string libs_;
+
         //- Calculated SHA1Digest
         SHA1Digest sha1_;
 
@@ -96,6 +100,12 @@ public:
             return options_;
         }
 
+        //- Return the code-libs
+        const string& libs() const
+        {
+            return libs_;
+        }
+
         //- Return the code
         const string& code() const
         {
@@ -114,6 +124,13 @@ public:
             return sha1_;
         }
 
+        //- Helper: add #line directive
+        static void addLineDirective
+        (
+            string&,
+            const label lineNum,
+            const fileName& name
+        );
 
 };
 

src/OpenFOAM/db/functionObjects/OutputFilterFunctionObject/OutputFilterFunctionObject.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
@@ -39,6 +39,7 @@ void Foam::OutputFilterFunctionObject<OutputFilter>::readDict()
     dict_.readIfPresent("storeFilter", storeFilter_);
 }
 
+
 template<class OutputFilter>
 void Foam::OutputFilterFunctionObject<OutputFilter>::allocateFilter()
 {
@@ -68,12 +69,14 @@ void Foam::OutputFilterFunctionObject<OutputFilter>::allocateFilter()
     }
 }
 
+
 template<class OutputFilter>
 void Foam::OutputFilterFunctionObject<OutputFilter>::destroyFilter()
 {
     ptr_.reset();
 }
 
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class OutputFilter>

src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedFieldIO.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
@@ -51,7 +51,7 @@ void Foam::DimensionedField<Type, GeoMesh>::readIfPresent
 {
     if
     (
-        (this->headerOk() && this->readOpt() == IOobject::READ_IF_PRESENT)
+        (this->readOpt() == IOobject::READ_IF_PRESENT && this->headerOk())
      || this->readOpt() == IOobject::MUST_READ
      || this->readOpt() == IOobject::MUST_READ_IF_MODIFIED
     )

src/OpenFOAM/include/OSspecific.H

@@ -196,6 +196,9 @@ void* dlSym(void* handle, const std::string& symbol);
 //- Report if symbol in a dlopened library could be found
 bool dlSymFound(void* handle, const std::string& symbol);
 
+//- Return all loaded libraries
+fileNameList dlLoaded();
+
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

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/finiteVolume/fields/fvPatchFields/derived/codedFixedValue/codedFixedValueFvPatchField.C

@@ -37,6 +37,9 @@ License
 #include "stringOps.H"
 #include "IOdictionary.H"
 
+#include <dlfcn.h>
+#include <link.h>
+
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 template<class Type>
@@ -238,10 +241,6 @@ void Foam::codedFixedValueFvPatchField<Type>::createLibrary
         // Write files for new library
         if (!dynCode.upToDate(context))
         {
-            Info<< "Using dynamicCode for patch " << this->patch().name()
-                << " on field " << this->dimensionedInternalField().name()
-                << endl;
-
             // filter with this context
             dynCode.reset(context);
 
@@ -268,9 +267,12 @@ void Foam::codedFixedValueFvPatchField<Type>::createLibrary
             dynCode.setMakeOptions
             (
                 "EXE_INC = -g \\\n"
-                "-I$(LIB_SRC)/finiteVolume/lnInclude\\\n"
+                "-I$(LIB_SRC)/finiteVolume/lnInclude \\\n"
               + context.options()
-              + "\n\nLIB_LIBS = "
+              + "\n\nLIB_LIBS = \\\n"
+              + "    -lOpenFOAM \\\n"
+              + "    -lfiniteVolume \\\n"
+              + context.libs()
             );
 
             if (!dynCode.copyOrCreateFiles(true))
@@ -337,6 +339,12 @@ void Foam::codedFixedValueFvPatchField<Type>::updateLibrary() const
         return;
     }
 
+    Info<< "Using dynamicCode for patch " << this->patch().name()
+        << " on field " << this->dimensionedInternalField().name() << nl
+        << "at line " << codeDict.startLineNumber()
+        << " in " << codeDict.name() << endl;
+
+
     // remove instantiation of fvPatchField provided by library
     redirectPatchFieldPtr_.clear();
 

src/finiteVolume/fields/fvPatchFields/derived/codedFixedValue/codedFixedValueFvPatchField.H

@@ -80,6 +80,9 @@ SourceFiles
 
 #include "fixedValueFvPatchFields.H"
 
+#include <dlfcn.h>
+#include <link.h>
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -119,6 +122,13 @@ class codedFixedValueFvPatchField
         //- Global loader/unloader function type
         typedef void (*loaderFunctionType)(bool);
 
+        static int collectLibsCallback
+        (
+            struct dl_phdr_info *info,
+            size_t size,
+            void *data
+        );
+
         //- Load specified library and execute globalFuncName(true)
         static void* loadLibrary
         (

src/finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValueFvPatchField.C

@@ -29,8 +29,8 @@ License
 #include "triSurface.H"
 #include "vector2D.H"
 #include "OFstream.H"
-#include "long.H"
 #include "AverageIOField.H"
+#include "Random.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -50,6 +50,7 @@ timeVaryingMappedFixedValueFvPatchField
     fixedValueFvPatchField<Type>(p, iF),
     fieldTableName_(iF.name()),
     setAverage_(false),
+    perturb_(0),
     referenceCS_(NULL),
     nearestVertex_(0),
     nearestVertexWeight_(0),
@@ -76,6 +77,7 @@ timeVaryingMappedFixedValueFvPatchField
     fixedValueFvPatchField<Type>(ptf, p, iF, mapper),
     fieldTableName_(ptf.fieldTableName_),
     setAverage_(ptf.setAverage_),
+    perturb_(ptf.perturb_),
     referenceCS_(NULL),
     nearestVertex_(0),
     nearestVertexWeight_(0),
@@ -101,6 +103,7 @@ timeVaryingMappedFixedValueFvPatchField
     fixedValueFvPatchField<Type>(p, iF),
     fieldTableName_(iF.name()),
     setAverage_(readBool(dict.lookup("setAverage"))),
+    perturb_(dict.lookupOrDefault("perturb", 1E-5)),
     referenceCS_(NULL),
     nearestVertex_(0),
     nearestVertexWeight_(0),
@@ -135,6 +138,7 @@ timeVaryingMappedFixedValueFvPatchField
     fixedValueFvPatchField<Type>(ptf),
     fieldTableName_(ptf.fieldTableName_),
     setAverage_(ptf.setAverage_),
+    perturb_(ptf.perturb_),
     referenceCS_(ptf.referenceCS_),
     nearestVertex_(ptf.nearestVertex_),
     nearestVertexWeight_(ptf.nearestVertexWeight_),
@@ -160,6 +164,7 @@ timeVaryingMappedFixedValueFvPatchField
     fixedValueFvPatchField<Type>(ptf, iF),
     fieldTableName_(ptf.fieldTableName_),
     setAverage_(ptf.setAverage_),
+    perturb_(ptf.perturb_),
     referenceCS_(ptf.referenceCS_),
     nearestVertex_(ptf.nearestVertex_),
     nearestVertexWeight_(ptf.nearestVertexWeight_),
@@ -332,7 +337,27 @@ void timeVaryingMappedFixedValueFvPatchField<Type>::readSamplePoints()
     (
         referenceCS().localPosition(samplePoints)
     );
-    const vectorField& localVertices = tlocalVertices();
+    vectorField& localVertices = tlocalVertices();
+
+    const boundBox bb(localVertices, true);
+    const point bbMid(bb.midpoint());
+
+    if (debug)
+    {
+        Info<< "timeVaryingMappedFixedValueFvPatchField :"
+            << " Perturbing points with " << perturb_
+            << " fraction of a random position inside " << bb
+            << " to break any ties on regular meshes."
+            << nl << endl;
+    }
+
+    Random rndGen(123456);
+    forAll(localVertices, i)
+    {
+        localVertices[i] +=
+            perturb_
+           *(rndGen.position(bb.min(), bb.max())-bbMid);
+    }
 
     // Determine triangulation
     List<vector2D> localVertices2D(localVertices.size());
@@ -342,33 +367,16 @@ void timeVaryingMappedFixedValueFvPatchField<Type>::readSamplePoints()
         localVertices2D[i][1] = localVertices[i][1];
     }
 
-    tmp<pointField> localFaceCentres
+    triSurface s(triSurfaceTools::delaunay2D(localVertices2D));
+
+    tmp<pointField> tlocalFaceCentres
     (
         referenceCS().localPosition
         (
             this->patch().patch().faceCentres()
         )
     );
-
+    const pointField& localFaceCentres = tlocalFaceCentres();
-    if (debug)
-    {
-        OFstream str
-        (
-            this->db().time().path()/this->patch().name()
-          + "_localFaceCentres.obj"
-        );
-        Pout<< "readSamplePoints :"
-            << " Dumping face centres to " << str.name() << endl;
-
-        forAll(localFaceCentres(), i)
-        {
-            const point& p = localFaceCentres()[i];
-            str<< "v " << p.x() << ' ' << p.y() << ' ' << p.z() << nl;
-        }
-    }
-
-
-    triSurface s(triSurfaceTools::delaunay2D(localVertices2D));
 
     if (debug)
     {
@@ -380,9 +388,9 @@ void timeVaryingMappedFixedValueFvPatchField<Type>::readSamplePoints()
         Pout<< "readSamplePoints :"
             << " Dumping face centres to " << str.name() << endl;
 
-        forAll(localFaceCentres(), i)
+        forAll(localFaceCentres, i)
         {
-            const point& p = localFaceCentres()[i];
+            const point& p = localFaceCentres[i];
             str<< "v " << p.x() << ' ' << p.y() << ' ' << p.z() << nl;
         }
     }
@@ -760,6 +768,7 @@ void timeVaryingMappedFixedValueFvPatchField<Type>::write(Ostream& os) const
 {
     fvPatchField<Type>::write(os);
     os.writeKeyword("setAverage") << setAverage_ << token::END_STATEMENT << nl;
+    os.writeKeyword("peturb") << perturb_ << token::END_STATEMENT << nl;
 
     if (fieldTableName_ != this->dimensionedInternalField().name())
     {

src/finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValueFvPatchField.H

@@ -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
@@ -27,6 +27,36 @@ Class
 Description
     Foam::timeVaryingMappedFixedValueFvPatchField
 
+    Interpolates from a set of supplied points in space and time. Supplied
+    data in constant/boundaryData/<patchname>:
+        - points : pointField with locations
+        - ddd    : supplied values at time ddd
+    Points need to be more or less on a plane since get triangulated in 2D.
+
+    At startup this bc does the triangulation and determines linear
+    interpolation (triangle it is in and weights to the 3 vertices)
+    for every face centre. Interpolates linearly inbetween times.
+
+    @verbatim
+        inlet
+        {
+            type            timeVaryingMappedFixedValue;
+
+            // Maintain average to that of the supplied values
+            setAverage      false;
+
+            // Optional: change perturbation (default 1E-5) to avoid any ties
+            // in triangulating regular geometries.
+            //perturb       0.0;
+
+            // Optional: use name instead of patchname for location of data
+            //fieldTableName samples;
+        }
+    @endverbatim
+
+    Switch on debug flag to have it dump the triangulation (in transformed
+    space) and transform face centres.
+
 SourceFiles
     timeVaryingMappedFixedValueFvPatchField.C
 
@@ -62,6 +92,9 @@ class timeVaryingMappedFixedValueFvPatchField
         //- If true adjust the mapped field to maintain average value
         bool setAverage_;
 
+        //- Fraction of perturbation (fraction of bounding box) to add
+        scalar perturb_;
+
         //- Coordinate system
         autoPtr<coordinateSystem> referenceCS_;
 

src/lagrangian/coalCombustion/coalParcel/makeCoalParcelSubmodels.C

@@ -50,27 +50,23 @@ License
 
 namespace Foam
 {
-    typedef coalCloud::cloudType coalCloud_R;
-    typedef coalCloud_R::cloudType coalCloud_T;
-    typedef coalCloud_T::cloudType coalCloud_K;
-
     // Kinematic sub-models
-    makeThermoParcelForces(coalCloud_K);
+    makeThermoParcelForces(coalCloud);
-    makeParcelDispersionModels(coalCloud_K);
+    makeParcelDispersionModels(coalCloud);
-    makeReactingMultiphaseParcelInjectionModels(coalCloud_K);
+    makeReactingMultiphaseParcelInjectionModels(coalCloud);
-    makeParcelPatchInteractionModels(coalCloud_K);
+    makeParcelPatchInteractionModels(coalCloud);
-    makeParcelPostProcessingModels(coalCloud_K);
+    makeParcelPostProcessingModels(coalCloud);
 
     // Thermo sub-models
-    makeParcelHeatTransferModels(coalCloud_T);
+    makeParcelHeatTransferModels(coalCloud);
 
     // Reacting sub-models
-    makeReactingMultiphaseParcelCompositionModels(coalCloud_R);
+    makeReactingMultiphaseParcelCompositionModels(coalCloud);
-    makeReactingParcelPhaseChangeModels(coalCloud_R);
+    makeReactingParcelPhaseChangeModels(coalCloud);
 
     // Reacting multiphase sub-models
     makeReactingMultiphaseParcelDevolatilisationModels(coalCloud);
-    makeReactingParcelSurfaceFilmModels(coalCloud_K);
+    makeReactingParcelSurfaceFilmModels(coalCloud);
     makeCoalParcelSurfaceReactionModels(coalCloud);
 }
 

src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationDiffusionLimitedRate/COxidationDiffusionLimitedRate.C

@@ -157,12 +157,15 @@ Foam::scalar Foam::COxidationDiffusionLimitedRate<CloudType>::calculate
     dMassSRCarrier[O2GlobalId_] -= dmO2;
     dMassSRCarrier[CO2GlobalId_] += dmCO2;
 
-    const scalar HC = thermo.solids().properties()[CsLocalId_].H(T);
+    const scalar HsC = thermo.solids().properties()[CsLocalId_].Hs(T);
     const scalar HCO2 = thermo.carrier().H(CO2GlobalId_, T);
-    const scalar HO2 = thermo.carrier().H(O2GlobalId_, T);
+
+    // carrier enthalpy transfer handled by change in composition
+    // const scalar HsO2 = thermo.carrier().Hs(O2GlobalId_, T);
+    // dhsTrans -= dmO2*HsO2;
 
     // Heat of reaction [J]
-    return dmC*HC + dmO2*HO2 - dmCO2*HCO2;
+    return dmC*HsC - dmCO2*HCO2;
 }
 
 

src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationKineticDiffusionLimitedRate/COxidationKineticDiffusionLimitedRate.C

@@ -160,12 +160,15 @@ Foam::scalar Foam::COxidationKineticDiffusionLimitedRate<CloudType>::calculate
     dMassSRCarrier[O2GlobalId_] -= dmO2;
     dMassSRCarrier[CO2GlobalId_] += dmCO2;
 
-    const scalar HC = thermo.solids().properties()[CsLocalId_].H(T);
+    const scalar HsC = thermo.solids().properties()[CsLocalId_].Hs(T);
     const scalar HCO2 = thermo.carrier().H(CO2GlobalId_, T);
-    const scalar HO2 = thermo.carrier().H(O2GlobalId_, T);
+
+    // carrier enthalpy transfer handled by change in composition
+    // const scalar HsO2 = thermo.carrier().Hs(O2GlobalId_, T);
+    // dhsTrans -= dmO2*HsO2;
 
     // Heat of reaction [J]
-    return dmC*HC + dmO2*HO2 - dmCO2*HCO2;
+    return dmC*HsC - dmCO2*HCO2;
 }
 
 

src/lagrangian/coalCombustion/submodels/surfaceReactionModel/COxidationMurphyShaddix/COxidationMurphyShaddix.C

@@ -220,12 +220,15 @@ Foam::scalar Foam::COxidationMurphyShaddix<CloudType>::calculate
     // Add to particle mass transfer
     dMassSolid[CsLocalId_] += dOmega*WC_;
 
-    const scalar HC = thermo.solids().properties()[CsLocalId_].H(T);
+    const scalar HsC = thermo.solids().properties()[CsLocalId_].Hs(T);
     const scalar HCO2 = thermo.carrier().H(CO2GlobalId_, T);
-    const scalar HO2 = thermo.carrier().H(O2GlobalId_, T);
+
+    // carrier enthalpy transfer handled by change in composition
+    // const scalar HsO2 = thermo.carrier().Hs(O2GlobalId_, T);
+    // dhsTrans -= dmO2*HsO2;
 
     // Heat of reaction
-    return dOmega*(WC_*HC + WO2_*HO2 - (WC_ + WO2_)*HCO2);
+    return dOmega*(WC_*HsC - (WC_ + WO2_)*HCO2);
 }
 
 

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

@@ -165,6 +165,11 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
     const label cellI
 )
 {
+    typedef typename TrackData::cloudType::reactingCloudType reactingCloudType;
+    const CompositionModel<reactingCloudType>& composition =
+        td.cloud().composition();
+
+
     // Define local properties at beginning of timestep
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     const scalar np0 = this->nParticle_;
@@ -178,9 +183,9 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
     const scalar pc = this->pc_;
 
     const scalarField& YMix = this->Y_;
-    const label idG = td.cloud().composition().idGas();
+    const label idG = composition.idGas();
-    const label idL = td.cloud().composition().idLiquid();
+    const label idL = composition.idLiquid();
-    const label idS = td.cloud().composition().idSolid();
+    const label idS = composition.idSolid();
 
 
     // Calc surface values
@@ -221,7 +226,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
     scalar NCpW = 0.0;
 
     // Surface concentrations of emitted species
-    scalarField Cs(td.cloud().composition().carrier().species().size(), 0.0);
+    scalarField Cs(composition.carrier().species().size(), 0.0);
 
     // Calc mass and enthalpy transfer due to phase change
     this->calcPhaseChange
@@ -271,10 +276,6 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
         Cs
     );
 
-    // Correct surface values due to emitted species
-    this->correctSurfaceValues(td, cellI, Ts, Cs, rhos, mus, Prs, kappas);
-    Res = this->Re(U0, d0, rhos, mus);
-
 
     // Surface reactions
     // ~~~~~~~~~~~~~~~~~
@@ -283,13 +284,9 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
     scalarField dMassSRGas(YGas_.size(), 0.0);
     scalarField dMassSRLiquid(YLiquid_.size(), 0.0);
     scalarField dMassSRSolid(YSolid_.size(), 0.0);
-    scalarField dMassSRCarrier
+    scalarField dMassSRCarrier(composition.carrier().species().size(), 0.0);
-    (
-        td.cloud().composition().carrier().species().size(),
-        0.0
-    );
 
-    // Clac mass and enthalpy transfer due to surface reactions
+    // Calc mass and enthalpy transfer due to surface reactions
     calcSurfaceReactions
     (
         td,
@@ -313,6 +310,11 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
     );
 
 
+    // Correct surface values due to emitted species
+    this->correctSurfaceValues(td, cellI, Ts, Cs, rhos, mus, Prs, kappas);
+    Res = this->Re(U0, d0, rhos, mus);
+
+
     // Update component mass fractions
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
@@ -324,11 +326,13 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
         updateMassFractions(mass0, dMassGas, dMassLiquid, dMassSolid);
 
 
+
+
     // Heat transfer
     // ~~~~~~~~~~~~~
 
     // Calculate new particle temperature
-    scalar Cuh = 0.0;
+    scalar Sph = 0.0;
     scalar T1 =
         this->calcHeatTransfer
         (
@@ -345,7 +349,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
             NCpW,
             Sh,
             dhsTrans,
-            Cuh
+            Sph
         );
 
 
@@ -380,25 +384,33 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
         // Transfer mass lost from particle to carrier mass source
         forAll(YGas_, i)
         {
-            label gid = td.cloud().composition().localToGlobalCarrierId(GAS, i);
+            label gid = composition.localToGlobalCarrierId(GAS, i);
             td.cloud().rhoTrans(gid)[cellI] += np0*dMassGas[i];
+//            td.cloud().hsTrans()[cellI] +=
+//                np0*dMassGas[i]*composition.carrier().Hs(gid, T0);
         }
         forAll(YLiquid_, i)
         {
-            label gid = td.cloud().composition().localToGlobalCarrierId(LIQ, i);
+            label gid = composition.localToGlobalCarrierId(LIQ, i);
             td.cloud().rhoTrans(gid)[cellI] += np0*dMassLiquid[i];
+//            td.cloud().hsTrans()[cellI] +=
+//                np0*dMassLiquid[i]*composition.carrier().Hs(gid, T0);
         }
 /*
         // No mapping between solid components and carrier phase
         forAll(YSolid_, i)
         {
-            label gid = td.cloud().composition().localToGlobalCarrierId(SLD, i);
+            label gid = composition.localToGlobalCarrierId(SLD, i);
             td.cloud().rhoTrans(gid)[cellI] += np0*dMassSolid[i];
+//            td.cloud().hsTrans()[cellI] +=
+//                np0*dMassSolid[i]*composition.carrier().Hs(gid, T0);
         }
 */
         forAll(dMassSRCarrier, i)
         {
             td.cloud().rhoTrans(i)[cellI] += np0*dMassSRCarrier[i];
+//            td.cloud().hsTrans()[cellI] +=
+//                np0*dMassSRCarrier[i]*composition.carrier().Hs(i, T0);
         }
 
         // Update momentum transfer
@@ -411,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;
     }
 
 
@@ -427,14 +439,12 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
             // Absorb parcel into carrier phase
             forAll(YGas_, i)
             {
-                label gid =
+                label gid = composition.localToGlobalCarrierId(GAS, i);
-                    td.cloud().composition().localToGlobalCarrierId(GAS, i);
                 td.cloud().rhoTrans(gid)[cellI] += np0*mass1*YMix[GAS]*YGas_[i];
             }
             forAll(YLiquid_, i)
             {
-                label gid =
+                label gid = composition.localToGlobalCarrierId(LIQ, i);
-                    td.cloud().composition().localToGlobalCarrierId(LIQ, i);
                 td.cloud().rhoTrans(gid)[cellI] +=
                     np0*mass1*YMix[LIQ]*YLiquid_[i];
             }
@@ -442,8 +452,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
             // No mapping between solid components and carrier phase
             forAll(YSolid_, i)
             {
-                label gid =
+                label gid = composition.localToGlobalCarrierId(SLD, i);
-                    td.cloud().composition().localToGlobalCarrierId(SLD, i);
                 td.cloud().rhoTrans(gid)[cellI] +=
                     np0*mass1*YMix[SLD]*YSolid_[i];
             }
@@ -508,6 +517,11 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcDevolatilisation
         return;
     }
 
+    typedef typename TrackData::cloudType::reactingCloudType reactingCloudType;
+    const CompositionModel<reactingCloudType>& composition =
+        td.cloud().composition();
+
+
     // Total mass of volatiles evolved
     td.cloud().devolatilisation().calculate
     (
@@ -535,10 +549,9 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcDevolatilisation
         // Note: hardcoded gaseous diffusivities for now
         // TODO: add to carrier thermo
         const scalar beta = sqr(cbrt(15.0) + cbrt(15.0));
-        const label id =
+        const label id = composition.localToGlobalCarrierId(GAS, i);
-            td.cloud().composition().localToGlobalCarrierId(GAS, i);
+        const scalar Cp = composition.carrier().Cp(id, Ts);
-        const scalar Cp = td.cloud().composition().carrier().Cp(id, Ts);
+        const scalar W = composition.carrier().W(id);
-        const scalar W = td.cloud().composition().carrier().W(id);
         const scalar Ni = dMassDV[i]/(this->areaS(d)*dt*W);
 
         // Dab calc'd using API vapour mass diffusivity function

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

@@ -252,6 +252,11 @@ void Foam::ReactingParcel<ParcelType>::calc
     const label cellI
 )
 {
+    typedef typename TrackData::cloudType::reactingCloudType reactingCloudType;
+    const CompositionModel<reactingCloudType>& composition =
+        td.cloud().composition();
+
+
     // Define local properties at beginning of time step
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     const scalar np0 = this->nParticle_;
@@ -301,7 +306,7 @@ void Foam::ReactingParcel<ParcelType>::calc
     scalar NCpW = 0.0;
 
     // Surface concentrations of emitted species
-    scalarField Cs(td.cloud().composition().carrier().species().size(), 0.0);
+    scalarField Cs(composition.carrier().species().size(), 0.0);
 
     // Calc mass and enthalpy transfer due to phase change
     calcPhaseChange
@@ -337,7 +342,7 @@ void Foam::ReactingParcel<ParcelType>::calc
     // ~~~~~~~~~~~~~
 
     // Calculate new particle temperature
-    scalar Cuh = 0.0;
+    scalar Sph = 0.0;
     scalar T1 =
         this->calcHeatTransfer
         (
@@ -354,7 +359,7 @@ void Foam::ReactingParcel<ParcelType>::calc
             NCpW,
             Sh,
             dhsTrans,
-            Cuh
+            Sph
         );
 
 
@@ -389,8 +394,10 @@ void Foam::ReactingParcel<ParcelType>::calc
         // Transfer mass lost from particle to carrier mass source
         forAll(dMassPC, i)
         {
-            label gid = td.cloud().composition().localToGlobalCarrierId(0, i);
+            label gid = composition.localToGlobalCarrierId(0, i);
             td.cloud().rhoTrans(gid)[cellI] += np0*dMassPC[i];
+//            td.cloud().hsTrans()[cellI] +=
+//                np0*dMassPC[i]*composition.carrier().Hs(gid, T0);
         }
 
         // Update momentum transfer
@@ -403,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;
     }
 
 
@@ -418,13 +425,12 @@ void Foam::ReactingParcel<ParcelType>::calc
             // Absorb parcel into carrier phase
             forAll(Y_, i)
             {
-                label gid =
+                label gid = composition.localToGlobalCarrierId(0, i);
-                    td.cloud().composition().localToGlobalCarrierId(0, i);
                 td.cloud().rhoTrans(gid)[cellI] += np0*mass1*Y_[i];
             }
             td.cloud().UTrans()[cellI] += np0*mass1*U1;
             td.cloud().hsTrans()[cellI] +=
-                np0*mass1*td.cloud().composition().H(0, Y_, pc_, T1);
+                np0*mass1*composition.H(0, Y_, pc_, T1);
         }
     }
 
@@ -434,7 +440,7 @@ void Foam::ReactingParcel<ParcelType>::calc
 
     else
     {
-        this->Cp_ = td.cloud().composition().Cp(0, Y_, pc_, T1);
+        this->Cp_ = composition.Cp(0, Y_, pc_, T1);
         this->T_ = T1;
         this->U_ = U1;
 
@@ -484,6 +490,11 @@ void Foam::ReactingParcel<ParcelType>::calcPhaseChange
         return;
     }
 
+    typedef typename TrackData::cloudType::reactingCloudType reactingCloudType;
+    const CompositionModel<reactingCloudType>& composition =
+        td.cloud().composition();
+
+
     // Calculate mass transfer due to phase change
     td.cloud().phaseChange().calculate
     (
@@ -511,19 +522,18 @@ void Foam::ReactingParcel<ParcelType>::calcPhaseChange
 
     forAll(YComponents, i)
     {
-        const label idc =
+        const label idc = composition.localToGlobalCarrierId(idPhase, i);
-            td.cloud().composition().localToGlobalCarrierId(idPhase, i);
+        const label idl = composition.globalIds(idPhase)[i];
-        const label idl = td.cloud().composition().globalIds(idPhase)[i];
 
         const scalar dh = td.cloud().phaseChange().dh(idc, idl, pc_, T);
         Sh -= dMassPC[i]*dh/dt;
 
         // Update particle surface thermo properties
         const scalar Dab =
-            td.cloud().composition().liquids().properties()[idl].D(pc_, Ts, Wc);
+            composition.liquids().properties()[idl].D(pc_, Ts, Wc);
 
-        const scalar Cp = td.cloud().composition().carrier().Cp(idc, Ts);
+        const scalar Cp = composition.carrier().Cp(idc, Ts);
-        const scalar W = td.cloud().composition().carrier().W(idc);
+        const scalar W = composition.carrier().W(idc);
         const scalar Ni = dMassPC[i]/(this->areaS(d)*dt*W);
 
         // Molar flux of species coming from the particle (kmol/m^2/s)

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

@@ -30,7 +30,7 @@ License
 
 template<class ParcelType>
 Foam::string Foam::ReactingParcel<ParcelType>::propHeader =
-    ThermoParcel<ParcelType>::propHeader
+    ParcelType::propHeader
   + " mass0"
   + " nPhases(Y1..YN)";
 

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

@@ -208,8 +208,8 @@ void Foam::ThermoParcel<ParcelType>::calc
     // Sum Ni*Cpi*Wi of emission species
     scalar NCpW = 0.0;
 
-    // Calculate new particle velocity
+    // Calculate new particle temperature
-    scalar Cuh = 0.0;
+    scalar Sph = 0.0;
     scalar T1 =
         this->calcHeatTransfer
         (
@@ -226,7 +226,7 @@ void Foam::ThermoParcel<ParcelType>::calc
             NCpW,
             Sh,
             dhsTrans,
-            Cuh
+            Sph
         );
 
 
@@ -267,7 +267,7 @@ void Foam::ThermoParcel<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;
     }
 
     // Set new particle properties
@@ -294,7 +294,7 @@ Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
     const scalar NCpW,
     const scalar Sh,
     scalar& dhsTrans,
-    scalar& Cuh
+    scalar& Sph
 )
 {
     if (!td.cloud().heatTransfer().active())
@@ -317,6 +317,7 @@ Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
 
     htc = max(htc, ROOTVSMALL);
     const scalar As = this->areaS(d);
+
     scalar ap = Tc_ + Sh/As/htc;
     scalar bp = 6.0*(Sh/As + htc*(Tc_ - T));
     if (td.cloud().radiation())
@@ -337,9 +338,9 @@ Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
 
     scalar Tnew = max(Tres.value(), td.cloud().constProps().TMin());
 
-    dhsTrans += dt*htc*As*(0.5*(T + Tnew) - Tc_);
+    Sph = dt*htc*As;
 
-    Cuh = dt*bp;
+    dhsTrans += Sph*(0.5*(T + Tnew) - Tc_);
 
     return Tnew;
 }

src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.H

@@ -243,7 +243,7 @@ protected:
             const scalar NCpW,         // Sum of N*Cp*W of emission species
             const scalar Sh,           // explicit particle enthalpy source
             scalar& dhsTrans,          // sensible enthalpy transfer to carrier
-            scalar& Cuh                // linearised heat transfer coefficient
+            scalar& Sph                // linearised heat transfer coefficient
         );
 
 

src/lagrangian/intermediate/parcels/derived/basicKinematicCollidingParcel/makeBasicKinematicCollidingParcelSubmodels.C

@@ -39,16 +39,13 @@ License
 namespace Foam
 {
     // Kinematic sub-models
-    typedef basicKinematicCollidingCloud::cloudType
+    makeParcelForces(basicKinematicCollidingCloud);
-        basicKinematicCollidingCloud_K;
+    makeParcelDispersionModels(basicKinematicCollidingCloud);
-
+    makeParcelInjectionModels(basicKinematicCollidingCloud);
-    makeParcelForces(basicKinematicCollidingCloud_K);
-    makeParcelDispersionModels(basicKinematicCollidingCloud_K);
-    makeParcelInjectionModels(basicKinematicCollidingCloud_K);
     makeParcelCollisionModels(basicKinematicCollidingCloud);
-    makeParcelPatchInteractionModels(basicKinematicCollidingCloud_K);
+    makeParcelPatchInteractionModels(basicKinematicCollidingCloud);
-    makeParcelPostProcessingModels(basicKinematicCollidingCloud_K);
+    makeParcelPostProcessingModels(basicKinematicCollidingCloud);
-    makeParcelSurfaceFilmModels(basicKinematicCollidingCloud_K);
+    makeParcelSurfaceFilmModels(basicKinematicCollidingCloud);
 }
 
 

src/lagrangian/intermediate/parcels/derived/basicReactingMultiphaseParcel/makeBasicReactingMultiphaseParcelSubmodels.C

@@ -48,29 +48,22 @@ License
 
 namespace Foam
 {
-    typedef basicReactingMultiphaseCloud::cloudType
-        basicReactingMultiphaseCloud_R;
-    typedef basicReactingMultiphaseCloud_R::cloudType
-        basicReactingMultiphaseCloud_T;
-    typedef basicReactingMultiphaseCloud_T::cloudType
-        basicReactingMultiphaseCloud_K;
-
     // Kinematic sub-models
-    makeThermoParcelForces(basicReactingMultiphaseCloud_K);
+    makeThermoParcelForces(basicReactingMultiphaseCloud);
-    makeParcelDispersionModels(basicReactingMultiphaseCloud_K);
+    makeParcelDispersionModels(basicReactingMultiphaseCloud);
-    makeReactingMultiphaseParcelInjectionModels(basicReactingMultiphaseCloud_K);
+    makeReactingMultiphaseParcelInjectionModels(basicReactingMultiphaseCloud);
-    makeParcelPatchInteractionModels(basicReactingMultiphaseCloud_K);
+    makeParcelPatchInteractionModels(basicReactingMultiphaseCloud);
-    makeParcelPostProcessingModels(basicReactingMultiphaseCloud_K);
+    makeParcelPostProcessingModels(basicReactingMultiphaseCloud);
 
     // Thermo sub-models
-    makeParcelHeatTransferModels(basicReactingMultiphaseCloud_T);
+    makeParcelHeatTransferModels(basicReactingMultiphaseCloud);
 
     // Reacting sub-models
     makeReactingMultiphaseParcelCompositionModels
     (
-        basicReactingMultiphaseCloud_R
+        basicReactingMultiphaseCloud
     );
-    makeReactingParcelPhaseChangeModels(basicReactingMultiphaseCloud_R);
+    makeReactingParcelPhaseChangeModels(basicReactingMultiphaseCloud);
 
     // Reacting multiphase sub-models
     makeReactingMultiphaseParcelDevolatilisationModels
@@ -79,7 +72,7 @@ namespace Foam
     );
     makeReactingParcelSurfaceFilmModels
     (
-        basicReactingMultiphaseCloud_K
+        basicReactingMultiphaseCloud
     );
     makeReactingMultiphaseParcelSurfaceReactionModels
     (

src/lagrangian/intermediate/parcels/derived/basicReactingParcel/makeBasicReactingParcelSubmodels.C

@@ -44,23 +44,23 @@ License
 
 namespace Foam
 {
-    typedef basicReactingCloud::cloudType basicReactingCloud_T;
+    typedef basicReactingCloud::thermoCloudType thermoCloudType2;
-    typedef basicReactingCloud_T::cloudType basicReactingCloud_K;
+    typedef basicReactingCloud::kinematicCloudType kinematicCloudType2;
 
     // Kinematic sub-models
-    makeThermoParcelForces(basicReactingCloud_K);
+    makeThermoParcelForces(basicReactingCloud);
-    makeParcelDispersionModels(basicReactingCloud_K);
+    makeParcelDispersionModels(basicReactingCloud);
-    makeReactingParcelInjectionModels(basicReactingCloud_K);
+    makeReactingParcelInjectionModels(basicReactingCloud);
-    makeParcelPatchInteractionModels(basicReactingCloud_K);
+    makeParcelPatchInteractionModels(basicReactingCloud);
-    makeParcelPostProcessingModels(basicReactingCloud_K);
+    makeParcelPostProcessingModels(basicReactingCloud);
 
     // Thermo sub-models
-    makeParcelHeatTransferModels(basicReactingCloud_T);
+    makeParcelHeatTransferModels(basicReactingCloud);
 
     // Reacting sub-models
     makeReactingParcelCompositionModels(basicReactingCloud);
     makeReactingParcelPhaseChangeModels(basicReactingCloud);
-    makeReactingParcelSurfaceFilmModels(basicReactingCloud_K);
+    makeReactingParcelSurfaceFilmModels(basicReactingCloud);
 }
 
 

src/lagrangian/intermediate/parcels/derived/basicThermoParcel/makeBasicThermoParcelSubmodels.C

@@ -41,18 +41,15 @@ License
 namespace Foam
 {
     // Kinematic sub-models
-    typedef basicThermoCloud::cloudType basicThermoCloud_K;
+    makeThermoParcelForces(basicThermoCloud);
-
+    makeParcelDispersionModels(basicThermoCloud);
-    // Kinematic sub-models
+    makeParcelInjectionModels(basicThermoCloud);
-    makeThermoParcelForces(basicThermoCloud_K);
+    makeParcelPatchInteractionModels(basicThermoCloud);
-    makeParcelDispersionModels(basicThermoCloud_K);
+    makeParcelPostProcessingModels(basicThermoCloud);
-    makeParcelInjectionModels(basicThermoCloud_K);
-    makeParcelPatchInteractionModels(basicThermoCloud_K);
-    makeParcelPostProcessingModels(basicThermoCloud_K);
 
     // Thermo sub-models
     makeParcelHeatTransferModels(basicThermoCloud);
-    makeParcelSurfaceFilmModels(basicThermoCloud_K);
+    makeParcelSurfaceFilmModels(basicThermoCloud);
 }
 
 

src/lagrangian/intermediate/parcels/include/makeParcelDispersionModels.H

@@ -38,7 +38,12 @@ License
                                                                               \
     makeDispersionModel(CloudType);                                           \
                                                                               \
-    defineNamedTemplateTypeNameAndDebug(DispersionRASModel<CloudType>, 0);    \
+    typedef CloudType::kinematicCloudType kinematicCloudType;                 \
+    defineNamedTemplateTypeNameAndDebug                                       \
+    (                                                                         \
+        DispersionRASModel<kinematicCloudType>,                               \
+        0                                                                     \
+    );                                                                        \
                                                                               \
     makeDispersionModelType(NoDispersion, CloudType);                         \
     makeDispersionModelType(GradientDispersionRAS, CloudType);                \

src/lagrangian/intermediate/submodels/ForceTypes/ParticleForceList/ParticleForceList.C

@@ -24,6 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "ParticleForceList.H"
+#include "entry.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -57,26 +58,46 @@ Foam::ParticleForceList<CloudType>::ParticleForceList
 {
     if (readFields)
     {
-        const wordList activeForces(dict.lookup("activeForces"));
+        wordList modelNames(dict.toc());
-
-        wordHashSet models;
-        models.insert(activeForces);
 
         Info<< "Constructing particle forces" << endl;
-        if (models.size() > 0)
+
+        if (modelNames.size() > 0)
         {
-            this->setSize(models.size());
+            this->setSize(modelNames.size());
 
             label i = 0;
-            forAllConstIter(wordHashSet, models, iter)
+            forAllConstIter(IDLList<entry>, dict, iter)
+            {
+                const word& model = iter().keyword();
+                if (iter().isDict())
                 {
-                const word& model = iter.key();
                     this->set
                     (
-                    i,
+                        i++,
-                    ParticleForce<CloudType>::New(owner, mesh, dict, model)
+                        ParticleForce<CloudType>::New
+                        (
+                            owner,
+                            mesh,
+                            iter().dict(),
+                            model
+                        )
                     );
-                i++;
+                }
+                else
+                {
+                    this->set
+                    (
+                        i++,
+                        ParticleForce<CloudType>::New
+                        (
+                            owner,
+                            mesh,
+                            dict,
+                            model
+                        )
+                    );
+                }
             }
         }
         else

src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/CollisionModel/CollisionModel.H

@@ -141,18 +141,27 @@ public:
 
 #define makeCollisionModel(CloudType)                                         \
                                                                               \
-    defineNamedTemplateTypeNameAndDebug(CollisionModel<CloudType>, 0);        \
+    typedef CloudType::collidingCloudType collidingCloudType;                 \
-                                                                              \
+    defineNamedTemplateTypeNameAndDebug                                       \
-    defineTemplateRunTimeSelectionTable(CollisionModel<CloudType>, dictionary);
+    (                                                                         \
+        CollisionModel<collidingCloudType>,                                   \
+        0                                                                     \
+    );                                                                        \
+    defineTemplateRunTimeSelectionTable                                       \
+    (                                                                         \
+        CollisionModel<collidingCloudType>,                                   \
+        dictionary                                                            \
+    );
 
 
 #define makeCollisionModelType(SS, CloudType)                                 \
                                                                               \
-    defineNamedTemplateTypeNameAndDebug(SS<CloudType>, 0);                    \
+    typedef CloudType::collidingCloudType collidingCloudType;                 \
+    defineNamedTemplateTypeNameAndDebug(SS<collidingCloudType>, 0);           \
                                                                               \
-    CollisionModel<CloudType>::                                               \
+    CollisionModel<collidingCloudType>::                                      \
-        adddictionaryConstructorToTable<SS<CloudType> >                       \
+        adddictionaryConstructorToTable<SS<collidingCloudType> >              \
-            add##SS##CloudType##ConstructorToTable_;
+            add##SS##CloudType##collidingCloudType##ConstructorToTable_;
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/lagrangian/intermediate/submodels/Kinematic/DispersionModel/DispersionModel/DispersionModel.H

@@ -132,22 +132,23 @@ public:
 
 #define makeDispersionModel(CloudType)                                        \
                                                                               \
-    defineTemplateTypeNameAndDebug(DispersionModel<CloudType>, 0);            \
+    typedef CloudType::kinematicCloudType kinematicCloudType;                 \
-                                                                              \
+    defineTemplateTypeNameAndDebug(DispersionModel<kinematicCloudType>, 0);   \
     defineTemplateRunTimeSelectionTable                                       \
     (                                                                         \
-        DispersionModel<CloudType>,                                           \
+        DispersionModel<kinematicCloudType>,                                  \
         dictionary                                                            \
     );
 
 
 #define makeDispersionModelType(SS, CloudType)                                \
                                                                               \
-    defineNamedTemplateTypeNameAndDebug(SS<CloudType>, 0);                    \
+    typedef CloudType::kinematicCloudType kinematicCloudType;                 \
+    defineNamedTemplateTypeNameAndDebug(SS<kinematicCloudType>, 0);           \
                                                                               \
-    DispersionModel<CloudType>::                                              \
+    DispersionModel<kinematicCloudType>::                                     \
-        adddictionaryConstructorToTable<SS<CloudType> >                       \
+        adddictionaryConstructorToTable<SS<kinematicCloudType> >              \
-            add##SS##CloudType##ConstructorToTable_;
+            add##SS##CloudType##kinematicCloudType##ConstructorToTable_;
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.H

@@ -328,17 +328,27 @@ public:
 
 #define makeInjectionModel(CloudType)                                         \
                                                                               \
-    defineNamedTemplateTypeNameAndDebug(InjectionModel<CloudType>, 0);        \
+    typedef CloudType::kinematicCloudType kinematicCloudType;                 \
-                                                                              \
+    defineNamedTemplateTypeNameAndDebug                                       \
-    defineTemplateRunTimeSelectionTable(InjectionModel<CloudType>, dictionary);
+    (                                                                         \
+        InjectionModel<kinematicCloudType>,                                   \
+        0                                                                     \
+    );                                                                        \
+    defineTemplateRunTimeSelectionTable                                       \
+    (                                                                         \
+        InjectionModel<kinematicCloudType>,                                   \
+        dictionary                                                            \
+    );
 
 
 #define makeInjectionModelType(SS, CloudType)                                 \
                                                                               \
-    defineNamedTemplateTypeNameAndDebug(SS<CloudType>, 0);                    \
+    typedef CloudType::kinematicCloudType kinematicCloudType;                 \
+    defineNamedTemplateTypeNameAndDebug(SS<kinematicCloudType>, 0);           \
                                                                               \
-    InjectionModel<CloudType>::adddictionaryConstructorToTable<SS<CloudType> >\
+    InjectionModel<kinematicCloudType>::                                      \
-        add##SS##CloudType##ConstructorToTable_;
+        adddictionaryConstructorToTable<SS<kinematicCloudType> >              \
+        add##SS##CloudType##kinematicCloudType##ConstructorToTable_;
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/NonSphereDrag/NonSphereDragForce.C

@@ -41,11 +41,10 @@ Foam::NonSphereDragForce<CloudType>::NonSphereDragForce
 (
     CloudType& owner,
     const fvMesh& mesh,
-    const dictionary& dict,
+    const dictionary& dict
-    const word& forceType
 )
 :
-    ParticleForce<CloudType>(owner, mesh, dict, forceType),
+    ParticleForce<CloudType>(owner, mesh, dict, typeName, true),
     phi_(readScalar(this->coeffs().lookup("phi"))),
     a_(exp(2.3288 - 6.4581*phi_ + 2.4486*sqr(phi_))),
     b_(0.0964 + 0.5565*phi_),

src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/NonSphereDrag/NonSphereDragForce.H

@@ -114,8 +114,7 @@ public:
         (
             CloudType& owner,
             const fvMesh& mesh,
-            const dictionary& dict,
+            const dictionary& dict
-            const word& forceType
         );
 
         //- Construct copy

src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/SphereDrag/SphereDragForce.C

@@ -48,11 +48,10 @@ Foam::SphereDragForce<CloudType>::SphereDragForce
 (
     CloudType& owner,
     const fvMesh& mesh,
-    const dictionary& dict,
+    const dictionary& dict
-    const word& forceType
 )
 :
-    ParticleForce<CloudType>(owner, mesh, dict)
+    ParticleForce<CloudType>(owner, mesh, dict, typeName, false)
 {}
 
 

src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/SphereDrag/SphereDragForce.H

@@ -66,8 +66,7 @@ public:
         (
             CloudType& owner,
             const fvMesh& mesh,
-            const dictionary& dict,
+            const dictionary& dict
-            const word& forceType
         );
 
         //- Construct copy

src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Gravity/GravityForce.C

@@ -32,11 +32,10 @@ Foam::GravityForce<CloudType>::GravityForce
 (
     CloudType& owner,
     const fvMesh& mesh,
-    const dictionary& dict,
+    const dictionary& dict
-    const word& forceType
 )
 :
-    ParticleForce<CloudType>(owner, mesh, dict),
+    ParticleForce<CloudType>(owner, mesh, dict, typeName, false),
     g_(owner.g().value())
 {}
 

src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/Gravity/GravityForce.H

@@ -73,8 +73,7 @@ public:
         (
             CloudType& owner,
             const fvMesh& mesh,
-            const dictionary& dict,
+            const dictionary& dict
-            const word& forceType
         );
 
         //- Construct copy

src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/NonInertialFrame/NonInertialFrameForce.C

@@ -33,11 +33,10 @@ Foam::NonInertialFrameForce<CloudType>::NonInertialFrameForce
 (
     CloudType& owner,
     const fvMesh& mesh,
-    const dictionary& dict,
+    const dictionary& dict
-    const word& forceType
 )
 :
-    ParticleForce<CloudType>(owner, mesh, dict),
+    ParticleForce<CloudType>(owner, mesh, dict, typeName, true),
     WName_
     (
         this->coeffs().template lookupOrDefault<word>

src/lagrangian/intermediate/submodels/Kinematic/ParticleForces/NonInertialFrame/NonInertialFrameForce.H

@@ -94,8 +94,7 @@ public:
         (
             CloudType& owner,
             const fvMesh& mesh,
-            const dictionary& dict,
+            const dictionary& dict
-            const word& forceType
         );
 
         //- Construct copy