Merge remote branch 'OpenCFD/master' into olesenm

2025-11-28 03:28:01 +00:00 · 2010-11-08 08:07:29 +01:00
parent d5acd22a63 9c350ab501
commit 8b59f1d0d5
674 changed files with 19530 additions and 5336 deletions
--- a/applications/solvers/combustion/reactingFoam/chemistry.H
+++ b/applications/solvers/combustion/reactingFoam/chemistry.H
@ -11,7 +11,7 @@
    if (turbulentReaction)
    {
        volScalarField tk =
-                Cmix*sqrt(turbulence->muEff()/rho/turbulence->epsilon());
+            Cmix*sqrt(turbulence->muEff()/rho/turbulence->epsilon());
        volScalarField tc = chemistry.tc();

        // Chalmers PaSR model
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C
@ -31,7 +31,7 @@ Description
 \*---------------------------------------------------------------------------*/

 #include "fvCFD.H"
-#include "basicPsiThermo.H"
+#include "basicRhoThermo.H"
 #include "turbulenceModel.H"
 #include "fixedGradientFvPatchFields.H"
 #include "regionProperties.H"
@ -121,7 +121,7 @@ int main(int argc, char *argv[])
            << nl << endl;
    }

-    Info << "End\n" << endl;
+    Info<< "End\n" << endl;

    return 0;
 }
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/chtMultiRegionSimpleFoam.C
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/chtMultiRegionSimpleFoam.C
@ -30,7 +30,7 @@ Description
 \*---------------------------------------------------------------------------*/

 #include "fvCFD.H"
-#include "basicPsiThermo.H"
+#include "basicRhoThermo.H"
 #include "turbulenceModel.H"
 #include "fixedGradientFvPatchFields.H"
 #include "regionProperties.H"
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/createFluidFields.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/createFluidFields.H
@ -1,5 +1,5 @@
    // Initialise fluid field pointer lists
-    PtrList<basicPsiThermo> thermoFluid(fluidRegions.size());
+    PtrList<basicRhoThermo> thermoFluid(fluidRegions.size());
    PtrList<volScalarField> rhoFluid(fluidRegions.size());
    PtrList<volScalarField> KFluid(fluidRegions.size());
    PtrList<volVectorField> UFluid(fluidRegions.size());
@ -28,7 +28,7 @@
        thermoFluid.set
        (
            i,
-            basicPsiThermo::New(fluidRegions[i]).ptr()
+            basicRhoThermo::New(fluidRegions[i]).ptr()
        );

        Info<< "    Adding to rhoFluid\n" << endl;
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H
@ -1,5 +1,4 @@
 {
-    // From buoyantSimpleFoam
    rho = thermo.rho();
    rho = max(rho, rhoMin[i]);
    rho = min(rho, rhoMax[i]);
@ -13,6 +12,8 @@

    phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf());
    bool closedVolume = adjustPhi(phi, U, p_rgh);
+    dimensionedScalar compressibility = fvc::domainIntegrate(psi);
+    bool compressible = (compressibility.value() > SMALL);

    surfaceScalarField buoyancyPhi = rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf();
    phi -= buoyancyPhi;
@ -25,7 +26,11 @@
            fvm::laplacian(rhorAUf, p_rgh) == fvc::div(phi)
        );

-        p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
+        p_rghEqn.setReference
+        (
+            pRefCell,
+            compressible ? getRefCellValue(p_rgh, pRefCell) : pRefValue
+        );

        p_rghEqn.solve();

@ -50,10 +55,10 @@

    // For closed-volume cases adjust the pressure level
    // to obey overall mass continuity
-    if (closedVolume)
+    if (closedVolume && compressible)
    {
-        p += (initialMass - fvc::domainIntegrate(psi*p))
-            /fvc::domainIntegrate(psi);
+        p += (initialMass - fvc::domainIntegrate(thermo.rho()))
+            /compressibility;
        p_rgh = p - rho*gh;
    }

--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/setRegionFluidFields.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/setRegionFluidFields.H
@ -1,6 +1,6 @@
    const fvMesh& mesh = fluidRegions[i];

-    basicPsiThermo& thermo = thermoFluid[i];
+    basicRhoThermo& thermo = thermoFluid[i];
    volScalarField& rho = rhoFluid[i];
    volScalarField& K = KFluid[i];
    volVectorField& U = UFluid[i];
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleCourantNo.C
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleCourantNo.C
@ -34,16 +34,13 @@ Foam::scalar Foam::compressibleCourantNo
    const surfaceScalarField& phi
 )
 {
-    scalar CoNum = 0.0;
-    scalar meanCoNum = 0.0;
-
    scalarField sumPhi =
        fvc::surfaceSum(mag(phi))().internalField()
       /rho.internalField();

-    CoNum = 0.5*gMax(sumPhi/mesh.V().field())*runTime.deltaTValue();
+    scalar CoNum = 0.5*gMax(sumPhi/mesh.V().field())*runTime.deltaTValue();

-    meanCoNum =
+    scalar meanCoNum =
        0.5*(gSum(sumPhi)/gSum(mesh.V().field()))*runTime.deltaTValue();

    Info<< "Region: " << mesh.name() << " Courant Number mean: " << meanCoNum
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H
@ -1,5 +1,5 @@
    // Initialise fluid field pointer lists
-    PtrList<basicPsiThermo> thermoFluid(fluidRegions.size());
+    PtrList<basicRhoThermo> thermoFluid(fluidRegions.size());
    PtrList<volScalarField> rhoFluid(fluidRegions.size());
    PtrList<volScalarField> KFluid(fluidRegions.size());
    PtrList<volVectorField> UFluid(fluidRegions.size());
@ -23,7 +23,7 @@
        thermoFluid.set
        (
            i,
-            basicPsiThermo::New(fluidRegions[i]).ptr()
+            basicRhoThermo::New(fluidRegions[i]).ptr()
        );

        Info<< "    Adding to rhoFluid\n" << endl;
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H
@ -1,5 +1,7 @@
 {
    bool closedVolume = p_rgh.needReference();
+    dimensionedScalar compressibility = fvc::domainIntegrate(psi);
+    bool compressible = (compressibility.value() > SMALL);

    rho = thermo.rho();

@ -19,34 +21,48 @@

    phi = phiU - rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf();

-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
    {
-        fvScalarMatrix p_rghEqn
+        fvScalarMatrix p_rghDDtEqn
        (
-            fvm::ddt(psi, p_rgh) + fvc::ddt(psi, rho)*gh
+            fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
          + fvc::div(phi)
-          - fvm::laplacian(rhorAUf, p_rgh)
        );

-        p_rghEqn.solve
-        (
-            mesh.solver
+        // Thermodynamic density needs to be updated by psi*d(p) after the
+        // pressure solution - done in 2 parts. Part 1:
+        thermo.rho() -= psi*p_rgh;
+
+        for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+        {
+            fvScalarMatrix p_rghEqn
            (
-                p_rgh.select
+                p_rghDDtEqn
+              - fvm::laplacian(rhorAUf, p_rgh)
+            );
+
+            p_rghEqn.solve
+            (
+                mesh.solver
                (
+                    p_rgh.select
                    (
-                        oCorr == nOuterCorr-1
-                     && corr == nCorr-1
-                     && nonOrth == nNonOrthCorr
+                        (
+                           oCorr == nOuterCorr-1
+                        && corr == nCorr-1
+                        && nonOrth == nNonOrthCorr
+                        )
                    )
                )
-            )
-        );
+            );

-        if (nonOrth == nNonOrthCorr)
-        {
-            phi += p_rghEqn.flux();
+            if (nonOrth == nNonOrthCorr)
+            {
+                phi += p_rghEqn.flux();
+            }
        }
+
+        // Second part of thermodynamic density update
+        thermo.rho() += psi*p_rgh;
    }

    // Correct velocity field
@ -66,10 +82,10 @@

    // For closed-volume cases adjust the pressure and density levels
    // to obey overall mass continuity
-    if (closedVolume)
+    if (closedVolume && compressible)
    {
-        p += (initialMass - fvc::domainIntegrate(psi*p))
-            /fvc::domainIntegrate(psi);
+        p += (initialMass - fvc::domainIntegrate(thermo.rho()))
+            /compressibility;
        rho = thermo.rho();
        p_rgh = p - rho*gh;
    }
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H
@ -1,6 +1,6 @@
    fvMesh& mesh = fluidRegions[i];

-    basicPsiThermo& thermo = thermoFluid[i];
+    basicRhoThermo& thermo = thermoFluid[i];
    volScalarField& rho = rhoFluid[i];
    volScalarField& K = KFluid[i];
    volVectorField& U = UFluid[i];
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/readSolidMultiRegionPISOControls.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/readSolidMultiRegionPISOControls.H
@ -1,5 +0,0 @@
-    const dictionary& piso = solidRegions[i].solutionDict().subDict("PISO");
-
-    const int nNonOrthCorr =
-        piso.lookupOrDefault<int>("nNonOrthogonalCorrectors", 0);
-
--- a/applications/solvers/lagrangian/coalChemistryFoam/UEqn.H
+++ b/applications/solvers/lagrangian/coalChemistryFoam/UEqn.H
@ -5,8 +5,8 @@
      + turbulence->divDevRhoReff(U)
     ==
        rho.dimensionedInternalField()*g
-      + coalParcels.SU()
-      + limestoneParcels.SU()
+      + coalParcels.SU(U)
+      + limestoneParcels.SU(U)
    );

    UEqn.relax();
--- a/applications/solvers/lagrangian/coalChemistryFoam/YEqn.H
+++ b/applications/solvers/lagrangian/coalChemistryFoam/YEqn.H
@ -25,7 +25,7 @@ tmp<fv::convectionScheme<scalar> > mvConvection
              + mvConvection->fvmDiv(phi, Yi)
              - fvm::laplacian(turbulence->muEff(), Yi)
              ==
-                coalParcels.Srho(i)
+                coalParcels.SYi(i, Yi)
              + kappa*chemistry.RR(i)().dimensionedInternalField()
            );

--- a/applications/solvers/lagrangian/coalChemistryFoam/hsEqn.H
+++ b/applications/solvers/lagrangian/coalChemistryFoam/hsEqn.H
@ -6,8 +6,8 @@
      - fvm::laplacian(turbulence->alphaEff(), hs)
     ==
        DpDt
-     +  coalParcels.Sh()
-     +  limestoneParcels.Sh()
+     +  coalParcels.Sh(hs)
+     +  limestoneParcels.Sh(hs)
     +  enthalpySource.Su()
     +  radiation->Shs(thermo)
     +  chemistrySh
--- a/applications/solvers/lagrangian/coalChemistryFoam/rhoEqn.H
+++ b/applications/solvers/lagrangian/coalChemistryFoam/rhoEqn.H
@ -35,7 +35,7 @@ Description
        fvm::ddt(rho)
      + fvc::div(phi)
      ==
-        coalParcels.Srho()
+        coalParcels.Srho(rho)
    );
 }

--- a/applications/solvers/lagrangian/incompressibleUncoupledKinematicParcelDyMFoam/incompressibleUncoupledKinematicParcelDyMFoam.C
+++ b/applications/solvers/lagrangian/incompressibleUncoupledKinematicParcelDyMFoam/incompressibleUncoupledKinematicParcelDyMFoam.C
@ -51,7 +51,7 @@ int main(int argc, char *argv[])

    #include "setRootCase.H"
    #include "createTime.H"
-    #   include "createDynamicFvMesh.H"
+    #include "createDynamicFvMesh.H"

    #include "readGravitationalAcceleration.H"
    #include "createFields.H"
--- a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/UEqn.H
+++ b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/UEqn.H
@ -6,7 +6,7 @@
      + turbulence->divDevRhoReff(U)
     ==
        rho.dimensionedInternalField()*g
-      + parcels.SU()
+      + parcels.SU(U)
      + momentumSource.Su()
    );

--- a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/YEqn.H
+++ b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/YEqn.H
@ -26,7 +26,7 @@ tmp<fv::convectionScheme<scalar> > mvConvection
              + mvConvection->fvmDiv(phi, Yi)
              - fvm::laplacian(turbulence->muEff(), Yi)
              ==
-                parcels.Srho(i)
+                parcels.SYi(i, Yi)
              + kappa*chemistry.RR(i)().dimensionedInternalField()
              + massSource.Su(i),
                mesh.solver("Yi")
--- a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/hsEqn.H
+++ b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/hsEqn.H
@ -37,7 +37,7 @@
          - fvm::laplacian(turbulence->alphaEff(), hs)
         ==
            pWork()
-          + parcels.Sh()
+          + parcels.Sh(hs)
          + radiation->Shs(thermo)
          + energySource.Su()
          + chemistrySh
--- a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/rhoEqn.H
+++ b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/rhoEqn.H
@ -35,7 +35,7 @@ Description
        fvm::ddt(rho)
      + fvc::div(phi)
      ==
-        parcels.Srho()
+        parcels.Srho(rho)
      + massSource.SuTot()
    );

--- a/applications/solvers/lagrangian/reactingParcelFilmFoam/UEqn.H
+++ b/applications/solvers/lagrangian/reactingParcelFilmFoam/UEqn.H
@ -5,7 +5,7 @@
      + turbulence->divDevRhoReff(U)
     ==
        rho.dimensionedInternalField()*g
-      + parcels.SU()
+      + parcels.SU(U)
    );

    UEqn.relax();
--- a/applications/solvers/lagrangian/reactingParcelFilmFoam/YEqn.H
+++ b/applications/solvers/lagrangian/reactingParcelFilmFoam/YEqn.H
@ -25,7 +25,7 @@ tmp<fv::convectionScheme<scalar> > mvConvection
              + mvConvection->fvmDiv(phi, Yi)
              - fvm::laplacian(turbulence->muEff(), Yi)
              ==
-                parcels.Srho(i)
+                parcels.SYi(i, Yi)
              + surfaceFilm.Srho(i)
              + kappa*chemistry.RR(i)().dimensionedInternalField(),
                mesh.solver("Yi")
--- a/applications/solvers/lagrangian/reactingParcelFilmFoam/hsEqn.H
+++ b/applications/solvers/lagrangian/reactingParcelFilmFoam/hsEqn.H
@ -6,7 +6,7 @@
      - fvm::laplacian(turbulence->alphaEff(), hs)
     ==
        DpDt
-      + parcels.Sh()
+      + parcels.Sh(hs)
      + surfaceFilm.Sh()
      + radiation->Shs(thermo)
      + chemistrySh
--- a/applications/solvers/lagrangian/reactingParcelFilmFoam/rhoEqn.H
+++ b/applications/solvers/lagrangian/reactingParcelFilmFoam/rhoEqn.H
@ -35,7 +35,7 @@ Description
        fvm::ddt(rho)
      + fvc::div(phi)
      ==
-        parcels.Srho()
+        parcels.Srho(rho)
      + surfaceFilm.Srho()
    );
 }
--- a/applications/solvers/lagrangian/reactingParcelFoam/UEqn.H
+++ b/applications/solvers/lagrangian/reactingParcelFoam/UEqn.H
@ -5,7 +5,7 @@
      + turbulence->divDevRhoReff(U)
     ==
        rho.dimensionedInternalField()*g
-      + parcels.SU()
+      + parcels.SU(U)
    );

    UEqn.relax();
--- a/applications/solvers/lagrangian/reactingParcelFoam/YEqn.H
+++ b/applications/solvers/lagrangian/reactingParcelFoam/YEqn.H
@ -25,7 +25,7 @@ tmp<fv::convectionScheme<scalar> > mvConvection
              + mvConvection->fvmDiv(phi, Yi)
              - fvm::laplacian(turbulence->muEff(), Yi)
              ==
-                parcels.Srho(i)
+                parcels.SYi(i, Yi)
              + kappa*chemistry.RR(i)().dimensionedInternalField(),
                mesh.solver("Yi")
            );
--- a/applications/solvers/lagrangian/reactingParcelFoam/hsEqn.H
+++ b/applications/solvers/lagrangian/reactingParcelFoam/hsEqn.H
@ -6,7 +6,7 @@
      - fvm::laplacian(turbulence->alphaEff(), hs)
     ==
        DpDt
-     +  parcels.Sh()
+     +  parcels.Sh(hs)
     +  radiation->Shs(thermo)
     +  chemistrySh
    );
--- a/applications/solvers/lagrangian/reactingParcelFoam/rhoEqn.H
+++ b/applications/solvers/lagrangian/reactingParcelFoam/rhoEqn.H
@ -35,7 +35,7 @@ Description
        fvm::ddt(rho)
      + fvc::div(phi)
      ==
-        parcels.Srho()
+        parcels.Srho(rho)
    );
 }

--- a/applications/solvers/lagrangian/steadyReactingParcelFoam/Make/files
+++ b/applications/solvers/lagrangian/steadyReactingParcelFoam/Make/files
@ -0,0 +1,3 @@
+steadyReactingParcelFoam.C
+
+EXE = $(FOAM_APPBIN)/steadyReactingParcelFoam
--- a/applications/solvers/lagrangian/steadyReactingParcelFoam/Make/options
+++ b/applications/solvers/lagrangian/steadyReactingParcelFoam/Make/options
@ -0,0 +1,43 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I${LIB_SRC}/meshTools/lnInclude \
+    -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
+    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
+    -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
+    -I$(LIB_SRC)/lagrangian/coalCombustion/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/pdfs/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/liquids/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/liquidMixture/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/solids/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/solidMixture/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/thermophysicalFunctions/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
+    -I$(LIB_SRC)/ODE/lnInclude \
+    -I$(LIB_SRC)/surfaceFilmModels/lnInclude
+
+EXE_LIBS = \
+    -lfiniteVolume \
+    -lmeshTools \
+    -lcompressibleTurbulenceModel \
+    -lcompressibleRASModels \
+    -lcompressibleLESModels \
+    -llagrangian \
+    -llagrangianIntermediate \
+    -lspecie \
+    -lbasicThermophysicalModels \
+    -lliquids \
+    -lliquidMixture \
+    -lsolids \
+    -lsolidMixture \
+    -lthermophysicalFunctions \
+    -lreactionThermophysicalModels \
+    -lSLGThermo \
+    -lchemistryModel \
+    -lradiation \
+    -lODE \
+    -lsurfaceFilmModels
--- a/applications/solvers/lagrangian/steadyReactingParcelFoam/UEqn.H
+++ b/applications/solvers/lagrangian/steadyReactingParcelFoam/UEqn.H
@ -0,0 +1,19 @@
+    fvVectorMatrix UEqn
+    (
+//        fvm::ddt(rho, U)
+        pZones.ddt(rho, U)
+      + fvm::div(phi, U)
+      + turbulence->divDevRhoReff(U)
+     ==
+        rho.dimensionedInternalField()*g
+      + parcels.SU(U)
+      + momentumSource.Su()
+    );
+
+    pZones.addResistance(UEqn);
+
+    if (momentumPredictor)
+    {
+        solve(UEqn == -fvc::grad(p));
+    }
+
--- a/applications/solvers/lagrangian/steadyReactingParcelFoam/YEqn.H
+++ b/applications/solvers/lagrangian/steadyReactingParcelFoam/YEqn.H
@ -0,0 +1,47 @@
+
+tmp<fv::convectionScheme<scalar> > mvConvection
+(
+    fv::convectionScheme<scalar>::New
+    (
+        mesh,
+        fields,
+        phi,
+        mesh.divScheme("div(phi,Yi_h)")
+    )
+);
+
+
+if (solveSpecies)
+{
+    label inertIndex = -1;
+    volScalarField Yt = 0.0*Y[0];
+
+    forAll(Y, i)
+    {
+        if (Y[i].name() != inertSpecie)
+        {
+            volScalarField& Yi = Y[i];
+            solve
+            (
+                fvm::ddt(rho, Yi)
+              + mvConvection->fvmDiv(phi, Yi)
+              - fvm::laplacian(turbulence->muEff(), Yi)
+              ==
+                parcels.SYi(i, Yi)
+              + kappa*chemistry.RR(i)().dimensionedInternalField()
+              + massSource.Su(i),
+                mesh.solver("Yi")
+            );
+
+            Yi.max(0.0);
+            Yt += Yi;
+        }
+        else
+        {
+            inertIndex = i;
+        }
+    }
+
+    Y[inertIndex] = scalar(1) - Yt;
+    Y[inertIndex].max(0.0);
+}
--- a/applications/solvers/lagrangian/steadyReactingParcelFoam/chemistry.H
+++ b/applications/solvers/lagrangian/steadyReactingParcelFoam/chemistry.H
@ -0,0 +1,28 @@
+if (chemistry.chemistry())
+{
+    Info<< "Solving chemistry" << endl;
+
+    // update reaction rates
+    chemistry.calculate();
+
+    // turbulent time scale
+    if (turbulentReaction)
+    {
+        typedef DimensionedField<scalar, volMesh> dsfType;
+
+        const dimensionedScalar e0("e0", sqr(dimLength)/pow3(dimTime), SMALL);
+
+        const dsfType tk =
+            Cmix*sqrt(turbulence->muEff()/rho/(turbulence->epsilon() + e0));
+
+        const dsfType tc = chemistry.tc()().dimensionedInternalField();
+
+        kappa = tc/(tc + tk);
+    }
+    else
+    {
+        kappa = 1.0;
+    }
+
+    chemistrySh = kappa*chemistry.Sh()();
+}
--- a/applications/solvers/lagrangian/steadyReactingParcelFoam/createClouds.H
+++ b/applications/solvers/lagrangian/steadyReactingParcelFoam/createClouds.H
@ -0,0 +1,9 @@
+Info<< "\nConstructing reacting cloud" << endl;
+basicReactingMultiphaseCloud parcels
+(
+    "reactingCloud1",
+    rho,
+    U,
+    g,
+    slgThermo
+);
--- a/applications/solvers/lagrangian/steadyReactingParcelFoam/createExplicitSources.H
+++ b/applications/solvers/lagrangian/steadyReactingParcelFoam/createExplicitSources.H
@ -0,0 +1,27 @@
+Info<< "Creating mass source\n" << endl;
+scalarTimeActivatedExplicitSourceList massSource
+(
+    "mass",
+    mesh,
+    dimMass/dimTime/dimVolume,
+    composition.species()
+);
+
+
+Info<< "Creating momentum source\n" << endl;
+vectorTimeActivatedExplicitSourceList momentumSource
+(
+    "momentum",
+    mesh,
+    dimMass*dimVelocity/dimTime/dimVolume,
+    "U"
+);
+
+Info<< "Creating energy source\n" << endl;
+scalarTimeActivatedExplicitSourceList energySource
+(
+    "energy",
+    mesh,
+    dimEnergy/dimTime/dimVolume,
+    "h"
+);
--- a/applications/solvers/lagrangian/steadyReactingParcelFoam/createFields.H
+++ b/applications/solvers/lagrangian/steadyReactingParcelFoam/createFields.H
@ -0,0 +1,149 @@
+    Info<< "Reading thermophysical properties\n" << endl;
+
+    autoPtr<rhoChemistryModel> pChemistry
+    (
+        rhoChemistryModel::New(mesh)
+    );
+    rhoChemistryModel& chemistry = pChemistry();
+
+    hsReactionThermo& thermo = chemistry.thermo();
+
+    SLGThermo slgThermo(mesh, thermo);
+
+    basicMultiComponentMixture& composition = thermo.composition();
+    PtrList<volScalarField>& Y = composition.Y();
+
+    const word inertSpecie(thermo.lookup("inertSpecie"));
+
+    if (!composition.contains(inertSpecie))
+    {
+        FatalErrorIn(args.executable())
+            << "Specified inert specie '" << inertSpecie << "' not found in "
+            << "species list. Available species:" << composition.species()
+            << exit(FatalError);
+    }
+
+    volScalarField& p = thermo.p();
+    volScalarField& hs = thermo.hs();
+    const volScalarField& T = thermo.T();
+    const volScalarField& psi = thermo.psi();
+
+    volScalarField rho
+    (
+        IOobject
+        (
+            "rho",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        thermo.rho()
+    );
+
+    Info<< "\nReading field U\n" << endl;
+    volVectorField U
+    (
+        IOobject
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    #include "compressibleCreatePhi.H"
+
+    DimensionedField<scalar, volMesh> kappa
+    (
+        IOobject
+        (
+            "kappa",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("zero", dimless, 0.0)
+    );
+
+    dimensionedScalar rhoMax
+    (
+        mesh.solutionDict().subDict("SIMPLE").lookup("rhoMax")
+    );
+
+    dimensionedScalar rhoMin
+    (
+        mesh.solutionDict().subDict("SIMPLE").lookup("rhoMin")
+    );
+
+    Info<< "Creating turbulence model\n" << endl;
+    autoPtr<compressible::turbulenceModel> turbulence
+    (
+        compressible::turbulenceModel::New
+        (
+            rho,
+            U,
+            phi,
+            thermo
+        )
+    );
+
+    Info<< "Creating multi-variate interpolation scheme\n" << endl;
+    multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
+
+    forAll(Y, i)
+    {
+        fields.add(Y[i]);
+    }
+    fields.add(hs);
+
+    DimensionedField<scalar, volMesh> chemistrySh
+    (
+        IOobject
+        (
+            "chemistry::Sh",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("chemistrySh", dimEnergy/dimTime/dimVolume, 0.0)
+    );
+
+    volScalarField invTauFlow
+    (
+        IOobject
+        (
+            "invTauFlow",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("one", dimless/dimTime, 1),
+        zeroGradientFvPatchScalarField::typeName
+    );
+
+    Info<< "Creating field DpDt\n" << endl;
+    volScalarField DpDt
+    (
+        IOobject
+        (
+            "DpDt",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("zero", dimPressure/dimTime, 0.0)
+    );
+
+    #include "setPressureWork.H"
--- a/applications/solvers/lagrangian/steadyReactingParcelFoam/createPorousZones.H
+++ b/applications/solvers/lagrangian/steadyReactingParcelFoam/createPorousZones.H
@ -0,0 +1,3 @@
+    Info<< "Creating porous zones" << nl << endl;
+
+    porousZones pZones(mesh);
--- a/applications/solvers/lagrangian/steadyReactingParcelFoam/hsEqn.H
+++ b/applications/solvers/lagrangian/steadyReactingParcelFoam/hsEqn.H
@ -0,0 +1,23 @@
+{
+    fvScalarMatrix hsEqn
+    (
+        fvm::ddt(rho, hs)
+      + mvConvection->fvmDiv(phi, hs)
+      - fvm::laplacian(turbulence->alphaEff(), hs)
+     ==
+        DpDt
+      + parcels.Sh(hs)
+      + radiation->Shs(thermo)
+      + energySource.Su()
+      + chemistrySh
+    );
+
+    hsEqn.solve();
+
+    thermo.correct();
+
+    radiation->correct();
+
+    Info<< "T gas min/max   = " << min(T).value() << ", "
+        << max(T).value() << endl;
+}
--- a/applications/solvers/lagrangian/steadyReactingParcelFoam/pEqn.H
+++ b/applications/solvers/lagrangian/steadyReactingParcelFoam/pEqn.H
@ -0,0 +1,65 @@
+{
+    rho = thermo.rho();
+
+    // Thermodynamic density needs to be updated by psi*d(p) after the
+    // pressure solution - done in 2 parts. Part 1:
+    thermo.rho() -= psi*p;
+
+    volScalarField rAU = 1.0/UEqn.A();
+    U = rAU*UEqn.H();
+
+    if (pZones.size() > 0)
+    {
+        // ddtPhiCorr not well defined for cases with porosity
+        phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf());
+    }
+    else
+    {
+        phi =
+            fvc::interpolate(rho)
+           *(
+                (fvc::interpolate(U) & mesh.Sf())
+              + fvc::ddtPhiCorr(rAU, rho, U, phi)
+            );
+    }
+
+    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    {
+        fvScalarMatrix pEqn
+        (
+            fvc::ddt(rho) + psi*correction(fvm::ddt(p))
+          + fvc::div(phi)
+          - fvm::laplacian(rho*rAU, p)
+        ==
+            parcels.Srho()
+          + massSource.SuTot()
+        );
+
+        pEqn.solve();
+
+        if (nonOrth == nNonOrthCorr)
+        {
+            phi += pEqn.flux();
+        }
+    }
+
+    // Explicitly relax pressure for momentum corrector
+    p.relax();
+
+    Info<< "p min/max = " << min(p).value() << ", " << max(p).value() << endl;
+
+    // Second part of thermodynamic density update
+    thermo.rho() += psi*p;
+
+    #include "rhoEqn.H" // NOTE: flux and time scales now inconsistent
+    #include "compressibleContinuityErrs.H"
+
+    U -= rAU*fvc::grad(p);
+    U.correctBoundaryConditions();
+
+    rho = thermo.rho();
+    rho = max(rho, rhoMin);
+    rho = min(rho, rhoMax);
+
+    #include "setPressureWork.H"
+}
--- a/applications/solvers/lagrangian/steadyReactingParcelFoam/readAdditionalSolutionControls.H
+++ b/applications/solvers/lagrangian/steadyReactingParcelFoam/readAdditionalSolutionControls.H
@ -0,0 +1,7 @@
+dictionary additional = mesh.solutionDict().subDict("additional");
+
+bool eWork = additional.lookupOrDefault("eWork", true);
+bool hWork = additional.lookupOrDefault("hWork", true);
+
+// flag to activate solve transport for each specie (Y vector)
+bool solveSpecies = additional.lookupOrDefault("solveSpecies", true);
--- a/applications/solvers/lagrangian/steadyReactingParcelFoam/readChemistryProperties.H
+++ b/applications/solvers/lagrangian/steadyReactingParcelFoam/readChemistryProperties.H
@ -0,0 +1,23 @@
+// Info<< "Reading chemistry properties\n" << endl;
+
+IOdictionary chemistryProperties
+(
+    IOobject
+    (
+        "chemistryProperties",
+        runTime.constant(),
+        mesh,
+        IOobject::MUST_READ_IF_MODIFIED,
+        IOobject::NO_WRITE,
+        false
+    )
+);
+
+Switch turbulentReaction(chemistryProperties.lookup("turbulentReaction"));
+
+dimensionedScalar Cmix("Cmix", dimless, 1.0);
+
+if (turbulentReaction)
+{
+    chemistryProperties.lookup("Cmix") >> Cmix;
+}
--- a/applications/solvers/lagrangian/steadyReactingParcelFoam/readTimeControls.H
+++ b/applications/solvers/lagrangian/steadyReactingParcelFoam/readTimeControls.H
@ -0,0 +1,70 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+// Maximum flow Courant number
+scalar maxCo(readScalar(runTime.controlDict().lookup("maxCo")));
+
+// Maximum time scale
+scalar maxDeltaT = readScalar(runTime.controlDict().lookup("maxDeltaT"));
+
+// Smoothing parameter (0-1) when smoothing iterations > 0
+scalar alphaTauSmooth
+(
+    runTime.controlDict().lookupOrDefault("alphaTauSmooth", 0.1)
+);
+
+// Maximum change in cell density per iteration (relative to previous value)
+scalar alphaTauRho
+(
+    runTime.controlDict().lookupOrDefault("alphaTauRho", 0.05)
+);
+
+// Maximum change in cell velocity per iteration (relative to previous value)
+scalar alphaTauU
+(
+    runTime.controlDict().lookupOrDefault("alphaTauU", 0.05)
+);
+
+// Maximum change in cell temperature per iteration (relative to previous value)
+scalar alphaTauTemp
+(
+    runTime.controlDict().lookupOrDefault("alphaTauTemp", 0.05)
+);
+
+// Max specie mass fraction that can be consumed/gained per chemistry
+// integration step
+scalar alphaTauSpecie
+(
+    runTime.controlDict().lookupOrDefault("alphaTauSpecie", 0.05)
+);
+
+// Maximum unboundedness allowed (fraction of 1)
+scalar specieMaxUnbound
+(
+    runTime.controlDict().lookupOrDefault("specieMaxUnbound", 0.01)
+);
+
+
+// ************************************************************************* //
--- a/applications/solvers/lagrangian/steadyReactingParcelFoam/rhoEqn.H
+++ b/applications/solvers/lagrangian/steadyReactingParcelFoam/rhoEqn.H
@ -0,0 +1,48 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2008-2010 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/>.
+
+Global
+    rhoEqn
+
+Description
+    Solve the continuity for density.
+
+\*---------------------------------------------------------------------------*/
+
+{
+    fvScalarMatrix rhoEqn
+    (
+        fvm::ddt(rho)
+      + fvc::div(phi)
+      ==
+        parcels.Srho(rho)
+      + massSource.SuTot()
+    );
+
+    rhoEqn.solve();
+
+    Info<< "rho min/max = " << min(rho).value() << ", " << max(rho).value()
+        << endl;
+}
+
+// ************************************************************************* //
--- a/applications/solvers/lagrangian/steadyReactingParcelFoam/setPressureWork.H
+++ b/applications/solvers/lagrangian/steadyReactingParcelFoam/setPressureWork.H
@ -0,0 +1,10 @@
+DpDt == dimensionedScalar("zero", DpDt.dimensions(), 0.0);
+
+if (eWork)
+{
+    DpDt += -p*fvc::div(phi/fvc::interpolate(rho));
+}
+if (hWork)
+{
+    DpDt += fvc::div(phi/fvc::interpolate(rho)*fvc::interpolate(p));
+}
--- a/applications/solvers/lagrangian/steadyReactingParcelFoam/steadyReactingParcelFoam.C
+++ b/applications/solvers/lagrangian/steadyReactingParcelFoam/steadyReactingParcelFoam.C
@ -0,0 +1,123 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2008-2010 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
+    porousExplicitSourceReactingParcelFoam
+
+Description
+    Transient PISO solver for compressible, laminar or turbulent flow with
+    reacting multiphase Lagrangian parcels for porous media, including explicit
+    sources for mass, momentum and energy
+
+    The solver includes:
+    - reacting multiphase parcel cloud
+    - porous media
+    - mass, momentum and energy sources
+
+    Note: ddtPhiCorr not used here when porous zones are active
+    - not well defined for porous calculations
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "hReactionThermo.H"
+#include "turbulenceModel.H"
+#include "basicReactingMultiphaseCloud.H"
+#include "rhoChemistryModel.H"
+#include "chemistrySolver.H"
+#include "radiationModel.H"
+#include "porousZones.H"
+#include "timeActivatedExplicitSource.H"
+#include "SLGThermo.H"
+#include "fvcSmooth.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "readGravitationalAcceleration.H"
+    #include "readTimeControls.H"
+    #include "readAdditionalSolutionControls.H"
+    #include "createFields.H"
+    #include "createRadiationModel.H"
+    #include "createClouds.H"
+    #include "createExplicitSources.H"
+    #include "createPorousZones.H"
+    #include "initContinuityErrs.H"
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.run())
+    {
+        #include "readSIMPLEControls.H"
+        #include "readChemistryProperties.H"
+        #include "readAdditionalSolutionControls.H"
+        #include "readTimeControls.H"
+
+        runTime++;
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        p.storePrevIter();
+
+        // --- Pressure-velocity corrector
+        {
+            parcels.evolve();
+
+            #include "chemistry.H"
+
+            #include "timeScales.H"
+
+            #include "rhoEqn.H"
+            #include "UEqn.H"
+            #include "YEqn.H"
+            #include "hsEqn.H"
+
+            #include "pEqn.H"
+
+            turbulence->correct();
+        }
+
+        if (runTime.write())
+        {
+            chemistry.dQ()().write();
+        }
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return(0);
+}
+
+
+// ************************************************************************* //
--- a/applications/solvers/lagrangian/steadyReactingParcelFoam/timeScales.H
+++ b/applications/solvers/lagrangian/steadyReactingParcelFoam/timeScales.H
@ -0,0 +1,228 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+Info<< "Time scales min/max:" << endl;
+
+
+{
+    // Cache old time scale field
+    tmp<volScalarField> tinvTauFlow0
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "invTauFlow0",
+                runTime.timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            invTauFlow
+        )
+    );
+    const volScalarField& invTauFlow0 = tinvTauFlow0();
+
+
+    // Flow time scale
+    // ~~~~~~~~~~~~~~~
+    {
+        invTauFlow =
+            fvc::surfaceSum
+            (
+                mag(phi)*mesh.deltaCoeffs()/(maxCo*mesh.magSf())
+            )
+           /rho;
+
+        invTauFlow.max(1.0/maxDeltaT);
+
+        Info<< "    Flow        = " << gMin(1/invTauFlow.internalField()) << ", "
+            << gMax(1/invTauFlow.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);
+*/
+    }
+
+
+    // Temperature source time scale
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    {
+        scalarField tau =
+            runTime.deltaTValue()
+           *mag
+            (
+                DpDt
+              + parcels.hsTrans()/(mesh.V()*runTime.deltaT())
+              + energySource.Su()
+              + chemistrySh
+            )
+           /rho;
+
+        tau = alphaTauTemp*thermo.Cp()*T/(tau + ROOTVSMALL);
+
+        Info<< "    Temperature = " << min(maxDeltaT, gMin(tau)) << ", "
+            << min(maxDeltaT, gMax(tau)) << endl;
+
+        invTauFlow.internalField() = max(invTauFlow.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);
+    }
+
+
+    // Limit rate of change of time scale
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    // - 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);
+    tinvTauFlow0.clear();
+
+
+    // Limit the largest time scale
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    invTauFlow.max(1/maxDeltaT);
+
+
+    // Spatially smooth the time scale field
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    fvc::smooth(invTauFlow, alphaTauSmooth);
+
+    Info<< "    Overall     = " << min(1/invTauFlow).value()
+        << ", " << max(1/invTauFlow).value() << nl << endl;
+}
+
+
+// ************************************************************************* //
--- a/applications/solvers/multiphase/interFoam/LTSInterFoam/setrDeltaT.H
+++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/setrDeltaT.H
@ -103,7 +103,7 @@
        fvc::sweep(rDeltaT, alpha1, nAlphaSweepIter, alphaSpreadDiff);
    }

-    Info<< "Flow time scale min/max = "
+    Info<< "Smoothed flow time scale min/max = "
        << gMin(1/rDeltaT.internalField())
        << ", " << gMax(1/rDeltaT.internalField()) << endl;

@ -116,13 +116,12 @@
     && runTime.timeIndex() > runTime.startTimeIndex() + 1
    )
    {
-        Info<< "Damping rDeltaT" << endl;
        rDeltaT = rDeltaT0*max(rDeltaT/rDeltaT0, 1.0 - rDeltaTDampingCoeff);
-    }

-    Info<< "Flow time scale min/max = "
-        << gMin(1/rDeltaT.internalField())
-        << ", " << gMax(1/rDeltaT.internalField()) << endl;
+        Info<< "Damped flow time scale min/max = "
+            << gMin(1/rDeltaT.internalField())
+            << ", " << gMax(1/rDeltaT.internalField()) << endl;
+    }

    label nAlphaSubCycles
    (
--- a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict
+++ b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict
@ -1,7 +1,7 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
 | \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  1.6                                   |
+|  \\    /   O peration     | Version:  1.7.1                                 |
 |   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/
@ -106,7 +106,8 @@ castellatedMeshControls
    // Explicit feature edge refinement
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

-    // Specifies a level for any cell intersected by its edges.
+    // Specifies a level for any cell intersected by explicitly provided
+    // edges.
    // This is a featureEdgeMesh, read from constant/triSurface for now.
    features
    (
--- a/applications/utilities/postProcessing/lagrangian/steadyParticleTracks/Make/files
+++ b/applications/utilities/postProcessing/lagrangian/steadyParticleTracks/Make/files
@ -0,0 +1,3 @@
+steadyParticleTracks.C
+
+EXE = $(FOAM_APPBIN)/steadyParticleTracks
--- a/applications/utilities/postProcessing/lagrangian/steadyParticleTracks/Make/options
+++ b/applications/utilities/postProcessing/lagrangian/steadyParticleTracks/Make/options
@ -0,0 +1,9 @@
+EXE_INC = \
+    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude
+
+EXE_LIBS = \
+    -llagrangian \
+    -lmeshTools \
+    -lfiniteVolume
--- a/applications/utilities/postProcessing/lagrangian/steadyParticleTracks/createFields.H
+++ b/applications/utilities/postProcessing/lagrangian/steadyParticleTracks/createFields.H
@ -0,0 +1,16 @@
+word dictName(args.optionLookupOrDefault<word>("dict", "particleTrackDict"));
+
+IOdictionary propsDict
+(
+    IOobject
+    (
+        dictName,
+        runTime.constant(),
+        mesh,
+        IOobject::MUST_READ_IF_MODIFIED
+    )
+);
+
+word cloudName(propsDict.lookup("cloudName"));
+
+List<word> userFields(propsDict.lookup("fields"));
--- a/applications/utilities/postProcessing/lagrangian/steadyParticleTracks/particleTrackDict
+++ b/applications/utilities/postProcessing/lagrangian/steadyParticleTracks/particleTrackDict
@ -0,0 +1,23 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      particleTrackDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+cloudName       reactingCloud1Tracks;
+
+fields ( d U T );
+
+// ************************************************************************* //
+
--- a/applications/utilities/postProcessing/lagrangian/steadyParticleTracks/steadyParticleTracks.C
+++ b/applications/utilities/postProcessing/lagrangian/steadyParticleTracks/steadyParticleTracks.C
@ -0,0 +1,326 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2008-2010 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
+    steadyParticleTracks
+
+Description
+    Generates a VTK file of particle tracks for cases that were computed using
+    a steady-state cloud
+    NOTE: case must be re-constructed (if running in parallel) before use
+
+\*---------------------------------------------------------------------------*/
+
+#include "argList.H"
+#include "Cloud.H"
+#include "IOdictionary.H"
+#include "fvMesh.H"
+#include "Time.H"
+#include "timeSelector.H"
+#include "OFstream.H"
+#include "passiveParticleCloud.H"
+
+#include "SortableList.H"
+#include "IOobjectList.H"
+#include "PtrList.H"
+#include "Field.H"
+#include "steadyParticleTracksTemplates.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+using namespace Foam;
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+label validateFields
+(
+    const List<word>& userFields,
+    const IOobjectList& cloudObjs
+)
+{
+    List<bool> ok(userFields.size(), false);
+
+    forAll(userFields, i)
+    {
+        ok[i] = ok[i] || fieldOk<label>(cloudObjs, userFields[i]);
+        ok[i] = ok[i] || fieldOk<scalar>(cloudObjs, userFields[i]);
+        ok[i] = ok[i] || fieldOk<vector>(cloudObjs, userFields[i]);
+        ok[i] = ok[i] || fieldOk<sphericalTensor>(cloudObjs, userFields[i]);
+        ok[i] = ok[i] || fieldOk<symmTensor>(cloudObjs, userFields[i]);
+        ok[i] = ok[i] || fieldOk<tensor>(cloudObjs, userFields[i]);
+    }
+
+    label nOk = 0;
+    forAll(ok, i)
+    {
+        if (ok[i])
+        {
+            nOk++;
+        }
+        else
+        {
+            Info << "\n*** Warning: user specified field '" << userFields[i]
+                 << "' unavailable" << endl;
+        }
+    }
+
+    return nOk;
+}
+
+
+template<>
+void writeVTK(OFstream& os, const label& value)
+{
+    os  << value;
+}
+
+
+template<>
+void writeVTK(OFstream& os, const scalar& value)
+{
+    os  << value;
+}
+
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    argList::noParallel();
+    timeSelector::addOptions();
+    #include "addRegionOption.H"
+    argList::validOptions.insert("dict", "");
+
+    #include "setRootCase.H"
+
+    #include "createTime.H"
+    instantList timeDirs = timeSelector::select0(runTime, args);
+    #include "createNamedMesh.H"
+    #include "createFields.H"
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    fileName vtkPath(runTime.path()/"VTK");
+    mkDir(vtkPath);
+
+    typedef HashTable<label, labelPair, labelPair::Hash<> > trackTableType;
+
+    forAll(timeDirs, timeI)
+    {
+        runTime.setTime(timeDirs[timeI], timeI);
+        Info<< "Time = " << runTime.timeName() << endl;
+
+        fileName vtkTimePath(runTime.path()/"VTK"/runTime.timeName());
+        mkDir(vtkTimePath);
+
+        Info<< "    Reading particle positions" << endl;
+
+        PtrList<passiveParticle> particles(0);
+
+        // transfer particles to (more convenient) list
+        {
+            passiveParticleCloud ppc(mesh, cloudName);
+            Info<< "\n    Read " << returnReduce(ppc.size(), sumOp<label>())
+                << " particles" << endl;
+
+            particles.setSize(ppc.size());
+
+            label i = 0;
+            forAllIter(passiveParticleCloud, ppc, iter)
+            {
+                particles.set(i++, ppc.remove(&iter()));
+            }
+
+            // myCloud should now be empty
+        }
+
+        List<label> particleToTrack(particles.size());
+        label nTracks = 0;
+
+        {
+            trackTableType trackTable;
+            forAll(particles, i)
+            {
+                const label origProc = particles[i].origProc();
+                const label origId = particles[i].origId();
+
+                const trackTableType::const_iterator& iter =
+                    trackTable.find(labelPair(origProc, origId));
+
+                if (iter == trackTable.end())
+                {
+                    particleToTrack[i] = nTracks;
+                    trackTable.insert(labelPair(origProc, origId), nTracks);
+                    nTracks++;
+                }
+                else
+                {
+                    particleToTrack[i] = iter();
+                }
+            }
+        }
+
+        if (nTracks == 0)
+        {
+            Info<< "\n    No track data" << endl;
+        }
+        else
+        {
+            Info<< "\n    Generating " << nTracks << " tracks" << endl;
+
+            // determine length of each track
+            labelList trackLengths(nTracks, 0);
+            forAll(particleToTrack, i)
+            {
+                const label trackI = particleToTrack[i];
+                trackLengths[trackI]++;
+            }
+
+            // particle "age" property used to sort the tracks
+            List<SortableList<scalar> > agePerTrack(nTracks);
+
+            forAll(trackLengths, i)
+            {
+                const label length = trackLengths[i];
+                agePerTrack[i].setSize(length);
+            }
+
+            // store the particle age per track
+            IOobjectList cloudObjs
+            (
+                mesh,
+                runTime.timeName(),
+                cloud::prefix/cloudName
+            );
+
+            // TODO: gather age across all procs
+            {
+                tmp<scalarField> tage =
+                    readParticleField<scalar>("age", cloudObjs);
+                const scalarField& age = tage();
+                List<label> trackSamples(nTracks, 0);
+                forAll(particleToTrack, i)
+                {
+                    const label trackI = particleToTrack[i];
+                    const label sampleI = trackSamples[trackI];
+                    agePerTrack[trackI][sampleI] = age[i];
+                    trackSamples[trackI]++;
+                }
+                tage.clear();
+            }
+
+
+            if (Pstream::master())
+            {
+                OFstream os(vtkTimePath/"particleTracks.vtk");
+
+                Info<< "\n    Writing particle tracks to " << os.name() << endl;
+
+                label nPoints = sum(trackLengths);
+
+                os  << "# vtk DataFile Version 2.0" << nl
+                    << "particleTracks" << nl
+                    << "ASCII" << nl
+                    << "DATASET POLYDATA" << nl
+                    << "POINTS " << nPoints << " float" << nl;
+
+                Info<< "\n    Writing points" << endl;
+
+                {
+                    label offset = 0;
+                    forAll(agePerTrack, i)
+                    {
+                        agePerTrack[i].sort();
+                        const labelList& ids = agePerTrack[i].indices();
+
+                        forAll(ids, j)
+                        {
+                            const label localId = offset + ids[j];
+                            const vector& pos = particles[localId].position();
+                            os  << pos.x() << ' ' << pos.y() << ' ' << pos.z()
+                                << nl;
+                        }
+
+                        offset += trackLengths[i];
+                    }
+                }
+
+
+                // write track (line) connectivity to file
+
+                Info<< "\n    Writing track lines" << endl;
+                os  << "\nLINES " << nTracks << ' ' << nPoints + nTracks << nl;
+
+                // Write ids of track points to file
+                {
+                    label globalPtI = 0;
+                    forAll(agePerTrack, i)
+                    {
+                        os  << agePerTrack[i].size() << nl;
+
+                        forAll(agePerTrack[i], j)
+                        {
+                            os  << ' ' << globalPtI++;
+                            if (((j + 1) % 10 == 0) && (j != 0))
+                            {
+                                os  << nl;
+                            }
+                        }
+
+                        os  << nl;
+                    }
+                }
+
+
+                const label nFields = validateFields(userFields, cloudObjs);
+
+                os  << "POINT_DATA " << nPoints << nl
+                    << "FIELD attributes " << nFields << nl;
+
+                Info<< "\n    Processing fields" << nl << endl;
+
+                processFields<label>(os, agePerTrack, userFields, cloudObjs);
+                processFields<scalar>(os, agePerTrack, userFields, cloudObjs);
+                processFields<vector>(os, agePerTrack, userFields, cloudObjs);
+                processFields<sphericalTensor>
+                    (os, agePerTrack, userFields, cloudObjs);
+                processFields<symmTensor>
+                    (os, agePerTrack, userFields, cloudObjs);
+                processFields<tensor>(os, agePerTrack, userFields, cloudObjs);
+
+            }
+        }
+        Info<< endl;
+    }
+
+    Info<< "\ndone" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //
--- a/applications/utilities/postProcessing/lagrangian/steadyParticleTracks/steadyParticleTracksTemplates.C
+++ b/applications/utilities/postProcessing/lagrangian/steadyParticleTracks/steadyParticleTracksTemplates.C
@ -0,0 +1,187 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 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 "steadyParticleTracksTemplates.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class Type>
+bool fieldOk(const IOobjectList& cloudObjs, const word& name)
+{
+    IOobjectList objects(cloudObjs.lookupClass(IOField<Type>::typeName));
+
+    return (objects.lookup(name) != NULL);
+}
+
+
+template<class Type>
+tmp<Field<Type> > readParticleField
+(
+    const word& name,
+    const IOobjectList cloudObjs
+)
+{
+    IOobjectList objects(cloudObjs.lookupClass(IOField<Type>::typeName));
+
+    const IOobject* obj = objects.lookup(name);
+    if (obj != NULL)
+    {
+        IOField<Type> newField(*obj);
+        return tmp<Field<Type> >(new Field<Type>(newField.xfer()));
+    }
+
+    Info<< "error: cloud field name " << name << " not found" << endl;
+
+    return Field<Type>::null();
+}
+
+
+template<class Type>
+PtrList<List<Type> > readFields
+(
+    PtrList<List<Type> >& values,
+    const List<word>& fields,
+    const IOobjectList& cloudObjs
+)
+{
+    IOobjectList objects(cloudObjs.lookupClass(IOField<Type>::typeName));
+
+    label fieldI = 0;
+    forAllConstIter(IOobjectList, objects, iter)
+    {
+        const IOobject& obj = *iter();
+        forAll(fields, j)
+        {
+            if (obj.name() == fields[j])
+            {
+                Info<< "        reading field " << obj.name() << endl;
+                IOField<Type> newField(obj);
+                values.set(fieldI++, new List<Type>(newField.xfer()));
+                break;
+            }
+        }
+    }
+
+    return values;
+}
+
+
+template<class Type>
+void writeVTK(OFstream& os, const Type& value)
+{
+    os  << value.component(0);
+    for (label i=1; i<pTraits<Type>::nComponents; i++)
+    {
+        os  << ' ' << value.component(i);
+    }
+}
+
+
+template<class Type>
+void writeVTKFields
+(
+    OFstream& os,
+    const PtrList<List<Type> >& values,
+    const List<SortableList<scalar> >& agePerTrack,
+    const List<word>& fieldNames
+)
+{
+    label step = max(floor(8/pTraits<Type>::nComponents), 1);
+
+    forAll(values, fieldI)
+    {
+        Info<< "        writing field " << fieldNames[fieldI] << endl;
+        os  << nl << fieldNames[fieldI] << ' ' << pTraits<Type>::nComponents << ' '
+            << values[fieldI].size() << " float" << nl;
+        label offset = 0;
+        forAll(agePerTrack, trackI)
+        {
+            const List<label> ids = agePerTrack[trackI].indices() + offset;
+
+            List<Type> data(UIndirectList<Type>(values[fieldI], ids));
+            label nData = data.size() - 1;
+            forAll(data, i)
+            {
+                writeVTK<Type>(os, data[i]);
+                if (((i + 1) % step == 0) || (i == nData))
+                {
+                    os  << nl;
+                }
+                else
+                {
+                    os  << ' ';
+                }
+            }
+            offset += ids.size();
+        }
+    }
+}
+
+
+template<class Type>
+void processFields
+(
+    OFstream& os,
+    const List<SortableList<scalar> >& agePerTrack,
+    const List<word>& userFieldNames,
+    const IOobjectList& cloudObjs
+)
+{
+    IOobjectList objects(cloudObjs.lookupClass(IOField<Type>::typeName));
+
+    if (objects.size())
+    {
+        DynamicList<word> fieldNames(objects.size());
+        forAll(userFieldNames, i)
+        {
+            IOobject* obj = objects.lookup(userFieldNames[i]);
+            if (obj != NULL)
+            {
+                fieldNames.append(obj->name());
+            }
+        }
+        fieldNames.shrink();
+
+        PtrList<List<Type> > values(fieldNames.size());
+        readFields<Type>(values, fieldNames, cloudObjs);
+
+        writeVTKFields<Type>
+        (
+            os,
+            values,
+            agePerTrack,
+            fieldNames.xfer()
+        );
+    }
+}
+
+} // End namespace Foam
+
+// ************************************************************************* //
--- a/applications/utilities/postProcessing/lagrangian/steadyParticleTracks/steadyParticleTracksTemplates.H
+++ b/applications/utilities/postProcessing/lagrangian/steadyParticleTracks/steadyParticleTracksTemplates.H
@ -0,0 +1,93 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef steadyParticleTracksTemplates_H
+#define steadyParticleTracksTemplates_H
+
+#include "OFstream.H"
+#include "SortableList.H"
+#include "IOobjectList.H"
+#include "PtrList.H"
+#include "Field.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    bool fieldOk(const IOobjectList& cloudObjs, const word& name);
+
+    template<class Type>
+    tmp<Field<Type> > readParticleField
+    (
+        const word& name,
+        const IOobjectList cloudObjs
+    );
+
+    template<class Type>
+    PtrList<List<Type> > readFields
+    (
+        PtrList<List<Type> >& values,
+        const List<word>& fields,
+        const IOobjectList& cloudObjs
+    );
+
+    template<class Type>
+    void writeVTK(OFstream& os, const Type& value);
+
+    template<class Type>
+    void writeVTKFields
+    (
+        OFstream& os,
+        const PtrList<List<Type> >& values,
+        const List<SortableList<scalar> >& agePerTrack,
+        const List<word>& fieldNames
+    );
+
+    void processFields
+    (
+        OFstream& os,
+        const List<SortableList<scalar> >& agePerTrack,
+        const List<word>& userFieldNames,
+        const IOobjectList& cloudObjs
+    );
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "steadyParticleTracksTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
--- a/applications/utilities/postProcessing/miscellaneous/pdfPlot/createFields.H
+++ b/applications/utilities/postProcessing/miscellaneous/pdfPlot/createFields.H
@ -20,7 +20,7 @@
    const fileName pdfPath = runTime.path()/"pdf";
    mkDir(pdfPath);

-    Random rndGen(label(0));
+    cachedRandom rndGen(label(0), -1);

    autoPtr<pdfs::pdf> p(pdfs::pdf::New(pdfDictionary, rndGen));