Merge branch 'clouds'

applications/solvers/lagrangian/coalChemistryFoam/UEqn.H

@@ -5,8 +5,8 @@
       + turbulence->divDevRhoReff(U)
      ==
         rho.dimensionedInternalField()*g
-      + coalParcels.SU()
+      + coalParcels.SU(U)
-      + limestoneParcels.SU()
+      + limestoneParcels.SU(U)
     );
 
     UEqn.relax();

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

applications/solvers/lagrangian/coalChemistryFoam/hsEqn.H

@@ -6,8 +6,8 @@
       - fvm::laplacian(turbulence->alphaEff(), hs)
      ==
         DpDt
-     +  coalParcels.Sh()
+     +  coalParcels.Sh(hs)
-     +  limestoneParcels.Sh()
+     +  limestoneParcels.Sh(hs)
      +  enthalpySource.Su()
      +  radiation->Shs(thermo)
      +  chemistrySh

applications/solvers/lagrangian/coalChemistryFoam/rhoEqn.H

@@ -35,7 +35,7 @@ Description
         fvm::ddt(rho)
       + fvc::div(phi)
       ==
-        coalParcels.Srho()
+        coalParcels.Srho(rho)
     );
 }
 

applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/UEqn.H

@@ -6,7 +6,7 @@
       + turbulence->divDevRhoReff(U)
      ==
         rho.dimensionedInternalField()*g
-      + parcels.SU()
+      + parcels.SU(U)
       + momentumSource.Su()
     );
 

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")

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

applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/rhoEqn.H

@@ -35,7 +35,7 @@ Description
         fvm::ddt(rho)
       + fvc::div(phi)
       ==
-        parcels.Srho()
+        parcels.Srho(rho)
       + massSource.SuTot()
     );
 

applications/solvers/lagrangian/reactingParcelFilmFoam/UEqn.H

@@ -5,7 +5,7 @@
       + turbulence->divDevRhoReff(U)
      ==
         rho.dimensionedInternalField()*g
-      + parcels.SU()
+      + parcels.SU(U)
     );
 
     UEqn.relax();

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")

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

applications/solvers/lagrangian/reactingParcelFilmFoam/rhoEqn.H

@@ -35,7 +35,7 @@ Description
         fvm::ddt(rho)
       + fvc::div(phi)
       ==
-        parcels.Srho()
+        parcels.Srho(rho)
       + surfaceFilm.Srho()
     );
 }

applications/solvers/lagrangian/reactingParcelFoam/UEqn.H

@@ -5,7 +5,7 @@
       + turbulence->divDevRhoReff(U)
      ==
         rho.dimensionedInternalField()*g
-      + parcels.SU()
+      + parcels.SU(U)
     );
 
     UEqn.relax();

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")
             );

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

applications/solvers/lagrangian/reactingParcelFoam/rhoEqn.H

@@ -35,7 +35,7 @@ Description
         fvm::ddt(rho)
       + fvc::div(phi)
       ==
-        parcels.Srho()
+        parcels.Srho(rho)
     );
 }
 

applications/solvers/lagrangian/steadyReactingParcelFoam/Make/files

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

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

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));
+    }
+

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);
+}

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()();
+}

applications/solvers/lagrangian/steadyReactingParcelFoam/createClouds.H

@@ -0,0 +1,9 @@
+Info<< "\nConstructing reacting cloud" << endl;
+basicReactingMultiphaseCloud parcels
+(
+    "reactingCloud1",
+    rho,
+    U,
+    g,
+    slgThermo
+);

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"
+);

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"

applications/solvers/lagrangian/steadyReactingParcelFoam/createPorousZones.H

@@ -0,0 +1,3 @@
+    Info<< "Creating porous zones" << nl << endl;
+
+    porousZones pZones(mesh);

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;
+}

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"
+}

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

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;
+}

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)
+);
+
+
+// ************************************************************************* //

applications/solvers/lagrangian/steadyReactingParcelFoam/rhoEqn.H

@@ -21,34 +21,28 @@ License
     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.
+
 \*---------------------------------------------------------------------------*/
 
-template<class CloudType>
-const Foam::dictionary& Foam::PostProcessingModel<CloudType>::dict() const
 {
-    return dict_;
+    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;
 }
 
-
-template<class CloudType>
-const CloudType& Foam::PostProcessingModel<CloudType>::owner() const
-{
-    return owner_;
-}
-
-
-template<class CloudType>
-CloudType& Foam::PostProcessingModel<CloudType>::owner()
-{
-    return owner_;
-}
-
-
-template<class CloudType>
-const Foam::dictionary& Foam::PostProcessingModel<CloudType>::coeffDict() const
-{
-    return coeffDict_;
-}
-
-
 // ************************************************************************* //

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));
+}

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);
+}
+
+
+// ************************************************************************* //

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;
+}
+
+
+// ************************************************************************* //

applications/utilities/postProcessing/lagrangian/steadyParticleTracks/Make/files

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

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

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"));

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 );
+
+// ************************************************************************* //
+

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;
+}
+
+
+// ************************************************************************* //

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
+
+// ************************************************************************* //

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
+
+// ************************************************************************* //

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

src/OpenFOAM/Make/files

@@ -61,7 +61,8 @@ sha1 = primitives/hashes/SHA1
 $(sha1)/SHA1.C
 $(sha1)/SHA1Digest.C
 
-primitives/random/Random.C
+primitives/random/Random/Random.C
+primitives/random/cachedRandom/cachedRandom.H
 
 containers/HashTables/HashTable/HashTableCore.C
 containers/HashTables/StaticHashTable/StaticHashTableCore.C

src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedField.C

@@ -81,14 +81,20 @@ DimensionedField<Type, GeoMesh>::DimensionedField
 (
     const IOobject& io,
     const Mesh& mesh,
-    const dimensionSet& dims
+    const dimensionSet& dims,
+    const bool checkIOFlags
 )
 :
     regIOobject(io),
     Field<Type>(GeoMesh::size(mesh)),
     mesh_(mesh),
     dimensions_(dims)
-{}
+{
+    if (checkIOFlags)
+    {
+        readIfPresent();
+    }
+}
 
 
 template<class Type, class GeoMesh>
@@ -96,14 +102,20 @@ DimensionedField<Type, GeoMesh>::DimensionedField
 (
     const IOobject& io,
     const Mesh& mesh,
-    const dimensioned<Type>& dt
+    const dimensioned<Type>& dt,
+    const bool checkIOFlags
 )
 :
     regIOobject(io),
     Field<Type>(GeoMesh::size(mesh), dt.value()),
     mesh_(mesh),
     dimensions_(dt.dimensions())
-{}
+{
+    if (checkIOFlags)
+    {
+        readIfPresent();
+    }
+}
 
 
 template<class Type, class GeoMesh>

src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedField.H

@@ -94,6 +94,11 @@ private:
         dimensionSet dimensions_;
 
 
+    // Private Member Functions
+
+        void readIfPresent(const word& fieldDictEntry = "value");
+
+
 public:
 
     //- Runtime type information
@@ -122,7 +127,8 @@ public:
         (
             const IOobject&,
             const Mesh& mesh,
-            const dimensionSet&
+            const dimensionSet&,
+            const bool checkIOFlags = true
         );
 
         //- Construct from components
@@ -130,7 +136,8 @@ public:
         (
             const IOobject&,
             const Mesh& mesh,
-            const dimensioned<Type>&
+            const dimensioned<Type>&,
+            const bool checkIOFlags = true
         );
 
         //- Construct from Istream
@@ -141,12 +148,6 @@ public:
             const word& fieldDictEntry="value"
         );
 
-        void readField
-        (
-            const dictionary& fieldDict,
-            const word& fieldDictEntry="value"
-        );
-
         //- Construct as copy
         DimensionedField
         (
@@ -222,6 +223,12 @@ public:
 
     // Member Functions
 
+        void readField
+        (
+            const dictionary& fieldDict,
+            const word& fieldDictEntry = "value"
+        );
+
         //- Return mesh
         inline const Mesh& mesh() const;
 

src/OpenFOAM/fields/DimensionedFields/DimensionedField/DimensionedFieldIO.C

@@ -47,6 +47,21 @@ void DimensionedField<Type, GeoMesh>::readField
 }
 
 
+template<class Type, class GeoMesh>
+void DimensionedField<Type, GeoMesh>::readIfPresent(const word& fieldDictEntry)
+{
+    if
+    (
+        (this->headerOk() && this->readOpt() == IOobject::READ_IF_PRESENT)
+     || this->readOpt() == IOobject::MUST_READ
+     || this->readOpt() == IOobject::MUST_READ_IF_MODIFIED
+    )
+    {
+        readField(dictionary(readStream(typeName)), fieldDictEntry);
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class Type, class GeoMesh>

src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C

@@ -209,7 +209,7 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
     const word& patchFieldType
 )
 :
-    DimensionedField<Type, GeoMesh>(io, mesh, ds),
+    DimensionedField<Type, GeoMesh>(io, mesh, ds, false),
     timeIndex_(this->time().timeIndex()),
     field0Ptr_(NULL),
     fieldPrevIterPtr_(NULL),
@@ -240,7 +240,7 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
     const wordList& actualPatchTypes
 )
 :
-    DimensionedField<Type, GeoMesh>(io, mesh, ds),
+    DimensionedField<Type, GeoMesh>(io, mesh, ds, false),
     timeIndex_(this->time().timeIndex()),
     field0Ptr_(NULL),
     fieldPrevIterPtr_(NULL),
@@ -267,7 +267,7 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
     const word& patchFieldType
 )
 :
-    DimensionedField<Type, GeoMesh>(io, mesh, dt),
+    DimensionedField<Type, GeoMesh>(io, mesh, dt, false),
     timeIndex_(this->time().timeIndex()),
     field0Ptr_(NULL),
     fieldPrevIterPtr_(NULL),
@@ -297,7 +297,7 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
     const wordList& actualPatchTypes
 )
 :
-    DimensionedField<Type, GeoMesh>(io, mesh, dt),
+    DimensionedField<Type, GeoMesh>(io, mesh, dt, false),
     timeIndex_(this->time().timeIndex()),
     field0Ptr_(NULL),
     fieldPrevIterPtr_(NULL),
@@ -351,7 +351,7 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
     const Mesh& mesh
 )
 :
-    DimensionedField<Type, GeoMesh>(io, mesh, dimless),
+    DimensionedField<Type, GeoMesh>(io, mesh, dimless, false),
     timeIndex_(this->time().timeIndex()),
     field0Ptr_(NULL),
     fieldPrevIterPtr_(NULL),
@@ -392,7 +392,7 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
     Istream& is
 )
 :
-    DimensionedField<Type, GeoMesh>(io, mesh, dimless),
+    DimensionedField<Type, GeoMesh>(io, mesh, dimless, false),
     timeIndex_(this->time().timeIndex()),
     field0Ptr_(NULL),
     fieldPrevIterPtr_(NULL),
@@ -431,7 +431,7 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
     const dictionary& dict
 )
 :
-    DimensionedField<Type, GeoMesh>(io, mesh, dimless),
+    DimensionedField<Type, GeoMesh>(io, mesh, dimless, false),
     timeIndex_(this->time().timeIndex()),
     field0Ptr_(NULL),
     fieldPrevIterPtr_(NULL),

src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixOperations.C

@@ -202,9 +202,12 @@ void Foam::lduMatrix::operator+=(const lduMatrix& A)
     }
     else
     {
-        FatalErrorIn("lduMatrix::operator+=(const lduMatrix& A)")
+        if (debug > 1)
+        {
+            WarningIn("lduMatrix::operator+=(const lduMatrix& A)")
                 << "Unknown matrix type combination"
-            << abort(FatalError);
+                << endl;
+        }
     }
 }
 
@@ -270,9 +273,12 @@ void Foam::lduMatrix::operator-=(const lduMatrix& A)
     }
     else
     {
-        FatalErrorIn("lduMatrix::operator-=(const lduMatrix& A)")
+        if (debug > 1)
+        {
+            WarningIn("lduMatrix::operator-=(const lduMatrix& A)")
                 << "Unknown matrix type combination"
-            << abort(FatalError);
+                << endl;
+        }
     }
 }
 

src/OpenFOAM/primitives/random/cachedRandom/cachedRandom.C

@@ -0,0 +1,159 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "cachedRandom.H"
+#include <cstdlib>
+
+#if INT_MAX    != 2147483647
+#    error "INT_MAX    != 2147483647"
+#    error "The random number generator may not work!"
+#endif
+
+// * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * //
+
+Foam::scalar Foam::cachedRandom::scalar01()
+{
+    if (sampleI_ < 0)
+    {
+        return drand48();
+    }
+
+    if (sampleI_ == samples_.size() - 1)
+    {
+        scalar s = samples_[sampleI_];
+        sampleI_ = 0;
+        return s;
+    }
+    else
+    {
+        scalar s = samples_[sampleI_];
+        sampleI_++;
+        return s;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::cachedRandom::cachedRandom(const label seed, const label count)
+:
+    seed_(1),
+    samples_(0),
+    sampleI_(-1)
+{
+    if (seed > 1)
+    {
+        seed_ = seed;
+    }
+
+    // Samples will be cached if count > 0
+    if (count > 0)
+    {
+        samples_.setSize(count);
+        sampleI_ = 0;
+    }
+
+    // Initialise samples
+    srand48(seed_);
+    forAll(samples_, i)
+    {
+        samples_[i] = drand48();
+    }
+}
+
+
+Foam::cachedRandom::cachedRandom(const cachedRandom& cr, const bool reset)
+:
+    seed_(cr.seed_),
+    samples_(cr.samples_),
+    sampleI_(cr.sampleI_)
+{
+    if (sampleI_ == -1)
+    {
+        WarningIn
+        (
+            "Foam::cachedRandom::cachedRandom(const cachedRandom& cr)"
+        )   << "Copy constructor called, but samples not being cached. "
+            << "This may lead to non-repeatable behaviour" << endl;
+
+        srand48(seed_);
+    }
+    else if (reset)
+    {
+        sampleI_ = 0;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::cachedRandom::~cachedRandom()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<>
+Foam::label Foam::cachedRandom::sample01()
+{
+    return round(scalar01());
+}
+
+
+template<>
+Foam::scalar Foam::cachedRandom::sample01()
+{
+    return scalar01();
+}
+
+
+template<>
+Foam::label Foam::cachedRandom::position(const label& start, const label& end)
+{
+    return start + round(scalar01()*(end - start));
+}
+
+
+template<>
+Foam::scalar Foam::cachedRandom::position
+(
+    const scalar& start,
+    const scalar& end
+)
+{
+    return start + scalar01()*(end - start);
+}
+
+
+void Foam::cachedRandom::operator=(const cachedRandom& cr)
+{
+    seed_ = cr.seed_;
+    samples_ = cr.samples_;
+    sampleI_ = cr.sampleI_;
+}
+
+
+// ************************************************************************* //

src/OpenFOAM/primitives/random/cachedRandom/cachedRandom.H

@@ -0,0 +1,182 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::cachedRandom
+
+Description
+    Random number generator.
+
+    Pre-computes and caches samples on construction, so that when sample01()
+    is called, the function simply returns the next (pre-computed) sample. On
+    reaching the last sample, the sample sequence is repeated.
+
+    Constructed using a seed and sample count. If the supplied count is
+    negative, no caching is performed, and a new sample is generated on each
+    call to sample01().
+
+    Note: the copy constructor cannot be used if count = -1.
+
+SourceFiles
+    cachedRandomI.H
+    cachedRandom.C
+    cachedRandomTemplates.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef cachedRandom_H
+#define cachedRandom_H
+
+#include "scalarList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+class cachedRandom;
+
+/*---------------------------------------------------------------------------*\
+                       Class cachedRandom Declaration
+\*---------------------------------------------------------------------------*/
+
+class cachedRandom
+{
+private:
+
+    // Private data
+
+        //- Initial random number seed
+        label seed_;
+
+        //- List of scalar samples
+        scalarList samples_;
+
+        //- Current sample marker
+        label sampleI_;
+
+
+    // Private Member Functions
+
+        //- Returns the current sample
+        scalar scalar01();
+
+
+public:
+
+
+    // Constructors
+
+        //- Construct given seed and sample count
+        cachedRandom(const label seed, const label count);
+
+        //- Copy constructor with optional reset of sampleI
+        cachedRandom(const cachedRandom& cr, const bool reset = false);
+
+
+    // Destructor
+    ~cachedRandom();
+
+
+    // Member functions
+
+        // Access
+
+            //- Return const access to the initial random number seed
+            inline label seed() const;
+
+            //- Return const access to the list of samples
+            inline const scalarList& samples() const;
+
+            //- Return the current sample marker
+            inline label sampleI() const;
+
+
+        // Manipulation
+
+            //- Return non-const access to the sample marker
+            inline label& sampleI();
+
+
+        // Evaluation
+
+            //- Return a sample whose components lie in the range 0-1
+            template<class Type>
+            Type sample01();
+
+            //- Return a sample between start and end
+            template<class Type>
+            Type position(const Type& start, const Type& end);
+
+            //- Randomise value in the range 0-1
+            template<class Type>
+            void randomise01(Type& value);
+
+
+        // Operators
+
+            //- Assignment operator
+            void operator=(const cachedRandom& cr);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// Template specialisations
+
+template<>
+label cachedRandom::sample01<label>();
+
+template<>
+scalar cachedRandom::sample01<scalar>();
+
+template<>
+label cachedRandom::position<label>(const label& start, const label& end);
+
+template<>
+scalar cachedRandom::position<scalar>
+(
+    const scalar& start,
+    const scalar& end
+);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "cachedRandomI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "cachedRandomTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/OpenFOAM/primitives/random/cachedRandom/cachedRandomI.H

@@ -0,0 +1,55 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "cachedRandom.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+inline Foam::label Foam::cachedRandom::seed() const
+{
+    return seed_;
+}
+
+
+inline const Foam::scalarList& Foam::cachedRandom::samples() const
+{
+    return samples_;
+}
+
+
+inline Foam::label Foam::cachedRandom::sampleI() const
+{
+    return sampleI_;
+}
+
+
+inline Foam::label& Foam::cachedRandom::sampleI()
+{
+    return sampleI_;
+}
+
+
+// ************************************************************************* //

src/OpenFOAM/primitives/random/cachedRandom/cachedRandomTemplates.C

@@ -0,0 +1,65 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "cachedRandom.H"
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type>
+Type Foam::cachedRandom::sample01()
+{
+    Type value;
+    for (direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
+    {
+        value.component(cmpt) = scalar01();
+    }
+
+    return value;
+}
+
+
+template<class Type>
+Type Foam::cachedRandom::position(const Type& start, const Type& end)
+{
+    Type value(start);
+    for (direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
+    {
+        value.component(cmpt) +=
+            scalar01()*(end.component(cmpt) - start.component(cmpt));
+    }
+
+    return value;
+}
+
+
+template<class Type>
+void Foam::cachedRandom::randomise01(Type& value)
+{
+    value = sample01<Type>();
+}
+
+
+// ************************************************************************* //

src/lagrangian/basic/Cloud/Cloud.C

@@ -384,9 +384,19 @@ void Foam::Cloud<ParticleType>::deleteParticle(ParticleType& p)
 }
 
 
+template<class ParticleType>
+void Foam::Cloud<ParticleType>::cloudReset(const Cloud<ParticleType>& c)
+{
+    // Reset particle cound and particles only
+    // - not changing the cloud object registry or reference to the polyMesh
+    particleCount_ = 0;
+    IDLList<ParticleType>::operator=(c);
+}
+
+
 template<class ParticleType>
 template<class TrackingData>
-void Foam::Cloud<ParticleType>::move(TrackingData& td)
+void Foam::Cloud<ParticleType>::move(TrackingData& td, const scalar trackTime)
 {
     const polyBoundaryMesh& pbm = pMesh().boundaryMesh();
     const globalMeshData& pData = polyMesh_.globalData();
@@ -444,7 +454,7 @@ void Foam::Cloud<ParticleType>::move(TrackingData& td)
             ParticleType& p = pIter();
 
             // Move the particle
-            bool keepParticle = p.move(td);
+            bool keepParticle = p.move(td, trackTime);
 
             // If the particle is to be kept
             // (i.e. it hasn't passed through an inlet or outlet)

src/lagrangian/basic/Cloud/Cloud.H

@@ -318,10 +318,13 @@ public:
             //- Remove particle from cloud and delete
             void deleteParticle(ParticleType&);
 
+            //- Reset the particles
+            void cloudReset(const Cloud<ParticleType>& c);
+
             //- Move the particles
             //  passing the TrackingData to the track function
             template<class TrackingData>
-            void move(TrackingData& td);
+            void move(TrackingData& td, const scalar trackTime);
 
             //- Remap the cells of particles corresponding to the
             //  mesh topology change

src/lagrangian/basic/InteractionLists/InteractionLists.C

@@ -979,7 +979,10 @@ void Foam::InteractionLists<ParticleType>::fillReferredParticleCloud()
 
             forAllConstIter(typename IDLList<ParticleType>, refCell, iter)
             {
-                cloud_.addParticle(iter().clone().ptr());
+                cloud_.addParticle
+                (
+                    static_cast<ParticleType*>(iter().clone().ptr())
+                );
             }
         }
     }
@@ -1076,6 +1079,29 @@ void Foam::InteractionLists<ParticleType>::writeReferredWallFaces() const
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
+template<class ParticleType>
+Foam::InteractionLists<ParticleType>::InteractionLists(const polyMesh& mesh)
+:
+    mesh_(mesh),
+    cloud_(mesh_, "NULL_Cloud", IDLList<ParticleType>()),
+    writeCloud_(false),
+    cellMapPtr_(),
+    wallFaceMapPtr_(),
+    globalTransforms_(mesh_),
+    maxDistance_(0.0),
+    dil_(),
+    dwfil_(),
+    ril_(),
+    rilInverse_(),
+    cellIndexAndTransformToDistribute_(),
+    wallFaceIndexAndTransformToDistribute_(),
+    referredWallFaces_(),
+    UName_("unknown_UName"),
+    referredWallData_(),
+    referredParticles_()
+{}
+
+
 template<class ParticleType>
 Foam::InteractionLists<ParticleType>::InteractionLists
 (
@@ -1106,6 +1132,7 @@ Foam::InteractionLists<ParticleType>::InteractionLists
     buildInteractionLists();
 }
 
+
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 template<class ParticleType>

src/lagrangian/basic/InteractionLists/InteractionLists.H

@@ -216,6 +216,9 @@ public:
 
     // Constructors
 
+        //- Construct null from mesh
+        InteractionLists(const polyMesh& mesh);
+
         //- Construct and call function to create all information from
         //  the mesh
         InteractionLists

src/lagrangian/basic/Particle/Particle.C

@@ -198,6 +198,25 @@ Foam::Particle<ParticleType>::Particle(const Particle<ParticleType>& p)
 {}
 
 
+template<class ParticleType>
+Foam::Particle<ParticleType>::Particle
+(
+    const Particle<ParticleType>& p,
+    const Cloud<ParticleType>& c
+)
+:
+    cloud_(c),
+    position_(p.position_),
+    cellI_(p.cellI_),
+    faceI_(p.faceI_),
+    stepFraction_(p.stepFraction_),
+    tetFaceI_(p.tetFaceI_),
+    tetPtI_(p.tetPtI_),
+    origProc_(p.origProc_),
+    origId_(p.origId_)
+{}
+
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class ParticleType>
@@ -480,6 +499,9 @@ Foam::scalar Foam::Particle<ParticleType>::trackToFace
 
     } while (faceI_ < 0);
 
+    ParticleType& p = static_cast<ParticleType&>(*this);
+    p.hitFace(td);
+
     if (cloud_.internalFace(faceI_))
     {
         // Change tet ownership because a tri face has been crossed,
@@ -508,8 +530,6 @@ Foam::scalar Foam::Particle<ParticleType>::trackToFace
     }
     else
     {
-        ParticleType& p = static_cast<ParticleType&>(*this);
-
         label origFaceI = faceI_;
         label patchI = patch(faceI_);
 
@@ -622,9 +642,7 @@ Foam::scalar Foam::Particle<ParticleType>::trackToFace
                 vector nHat = wallTri.normal();
                 nHat /= mag(nHat);
 
-                const ParticleType& p = static_cast<const ParticleType&>(*this);
+                const scalar r = p.wallImpactDistance(nHat);
-
-                scalar r = p.wallImpactDistance(nHat);
 
                 // Removing (approximately) the wallTri normal
                 // component of the existing correction, to avoid the
@@ -675,6 +693,12 @@ void Foam::Particle<ParticleType>::transformProperties(const vector&)
 {}
 
 
+template<class ParticleType>
+template<class TrackData>
+void Foam::Particle<ParticleType>::hitFace(TrackData&)
+{}
+
+
 template<class ParticleType>
 template<class TrackData>
 bool Foam::Particle<ParticleType>::hitPatch

src/lagrangian/basic/Particle/Particle.H

@@ -219,6 +219,10 @@ protected:
 
     // Patch interactions
 
+        //- Overridable function to handle the particle hitting a face
+        template<class TrackData>
+        void hitFace(TrackData& td);
+
         //- Overridable function to handle the particle hitting a
         //  patch.  Executed before other patch-hitting functions.
         //  trackFraction is passed in to allow mesh motion to
@@ -361,10 +365,28 @@ public:
         //- Construct as a copy
         Particle(const Particle& p);
 
+        //- Construct as a copy
+        Particle(const Particle& p, const Cloud<ParticleType>& c);
+
         //- Construct a clone
-        autoPtr<ParticleType> clone() const
+        virtual autoPtr<Particle<ParticleType> > clone() const
         {
-            return autoPtr<Particle>(new Particle(*this));
+            return autoPtr<Particle<ParticleType> >
+            (
+                new Particle<ParticleType>(*this)
+            );
+        }
+
+        //- Construct a clone
+        virtual autoPtr<Particle<ParticleType> > clone
+        (
+            const Cloud<ParticleType>& c
+        ) const
+        {
+            return autoPtr<Particle<ParticleType> >
+            (
+                new Particle<ParticleType>(*this, c)
+            );
         }
 
 

src/lagrangian/basic/indexedParticle/indexedParticle.H

@@ -108,9 +108,12 @@ public:
         {}
 
         //- Construct and return a clone
-        autoPtr<indexedParticle> clone() const
+        virtual autoPtr<Particle<indexedParticle> > clone() const
         {
-            return autoPtr<indexedParticle>(new indexedParticle(*this));
+            return autoPtr<Particle<indexedParticle> >
+            (
+                new indexedParticle(*this)
+            );
         }
 
 

src/lagrangian/basic/passiveParticle/passiveParticle.H

@@ -99,9 +99,12 @@ public:
         {}
 
         //- Construct and return a clone
-        autoPtr<passiveParticle> clone() const
+        virtual autoPtr<Particle<passiveParticle> > clone() const
         {
-            return autoPtr<passiveParticle>(new passiveParticle(*this));
+            return autoPtr<Particle<passiveParticle> >
+            (
+                new passiveParticle(*this)
+            );
         }
 };
 

src/lagrangian/coalCombustion/coalParcel/coalParcel.C

@@ -111,6 +111,16 @@ Foam::coalParcel::coalParcel(const coalParcel& p)
 }
 
 
+Foam::coalParcel::coalParcel
+(
+    const coalParcel& p,
+    const ReactingMultiphaseCloud<coalParcel>& c
+)
+:
+    ReactingMultiphaseParcel<coalParcel>(p, c)
+{}
+
+
 // * * * * * * * * * * * * * * * *  Destructors  * * * * * * * * * * * * * * //
 
 Foam::coalParcel::~coalParcel()

src/lagrangian/coalCombustion/coalParcel/coalParcel.H

@@ -103,10 +103,33 @@ public:
         //- Construct as a copy
         coalParcel(const coalParcel& p);
 
-        //- Construct and return a clone
+        //- Construct as a copy
-        autoPtr<coalParcel> clone() const
+        coalParcel
+        (
+            const coalParcel& p,
+            const ReactingMultiphaseCloud<coalParcel>& c
+        );
+
+        //- Construct and return a (basic particle) clone
+        virtual autoPtr<Particle<coalParcel> > clone() const
         {
-            return autoPtr<coalParcel>(new coalParcel(*this));
+            return autoPtr<Particle<coalParcel> >(new coalParcel(*this));
+        }
+
+        //- Construct and return a (basic particle) clone
+        virtual autoPtr<Particle<coalParcel> > clone
+        (
+            const Cloud<coalParcel>& c
+        ) const
+        {
+            return autoPtr<Particle<coalParcel> >
+            (
+                new coalParcel
+                (
+                    *this,
+                    static_cast<const ReactingMultiphaseCloud<coalParcel>&>(c)
+                )
+            );
         }
 
 

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

@@ -26,6 +26,8 @@ License
 #include "COxidationDiffusionLimitedRate.H"
 #include "mathematicalConstants.H"
 
+using namespace Foam::constant;
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class CloudType>
@@ -36,8 +38,8 @@ Foam::COxidationDiffusionLimitedRate<CloudType>::COxidationDiffusionLimitedRate
 )
 :
     SurfaceReactionModel<CloudType>(dict, owner, typeName),
-    Sb_(dimensionedScalar(this->coeffDict().lookup("Sb")).value()),
+    Sb_(readScalar(this->coeffDict().lookup("Sb"))),
-    D_(dimensionedScalar(this->coeffDict().lookup("D")).value()),
+    D_(readScalar(this->coeffDict().lookup("D"))),
     CsLocalId_(-1),
     O2GlobalId_(owner.composition().globalCarrierId("O2")),
     CO2GlobalId_(owner.composition().globalCarrierId("CO2")),
@@ -67,9 +69,31 @@ Foam::COxidationDiffusionLimitedRate<CloudType>::COxidationDiffusionLimitedRate
         )   << "Stoichiometry of reaction, Sb, must be greater than zero" << nl
             << exit(FatalError);
     }
+
+    const scalar YCloc = owner.composition().Y0(idSolid)[CsLocalId_];
+    const scalar YSolidTot = owner.composition().YMixture0()[idSolid];
+    Info<< "    C(s): particle mass fraction = " << YCloc*YSolidTot << endl;
 }
 
 
+template<class CloudType>
+Foam::COxidationDiffusionLimitedRate<CloudType>::COxidationDiffusionLimitedRate
+(
+    const COxidationDiffusionLimitedRate<CloudType>& srm
+)
+:
+    SurfaceReactionModel<CloudType>(srm),
+    Sb_(srm.Sb_),
+    D_(srm.D_),
+    CsLocalId_(srm.CsLocalId_),
+    O2GlobalId_(srm.O2GlobalId_),
+    CO2GlobalId_(srm.CO2GlobalId_),
+    WC_(srm.WC_),
+    WO2_(srm.WO2_),
+    HcCO2_(srm.HcCO2_)
+{}
+
+
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 template<class CloudType>
@@ -80,13 +104,6 @@ Foam::COxidationDiffusionLimitedRate<CloudType>::
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-template<class CloudType>
-bool Foam::COxidationDiffusionLimitedRate<CloudType>::active() const
-{
-    return true;
-}
-
-
 template<class CloudType>
 Foam::scalar Foam::COxidationDiffusionLimitedRate<CloudType>::calculate
 (
@@ -119,12 +136,13 @@ Foam::scalar Foam::COxidationDiffusionLimitedRate<CloudType>::calculate
         return 0.0;
     }
 
+    const SLGThermo& thermo = this->owner().thermo();
+
     // Local mass fraction of O2 in the carrier phase
-    const scalar YO2 = this->owner().thermo().carrier().Y(O2GlobalId_)[cellI];
+    const scalar YO2 = thermo.carrier().Y(O2GlobalId_)[cellI];
 
     // Change in C mass [kg]
-    scalar dmC =
+    scalar dmC = 4.0*mathematical::pi*d*D_*YO2*Tc*rhoc/(Sb_*(T + Tc))*dt;
-        4.0*constant::mathematical::pi*d*D_*YO2*Tc*rhoc/(Sb_*(T + Tc))*dt;
 
     // Limit mass transfer by availability of C
     dmC = min(mass*fComb, dmC);
@@ -142,8 +160,12 @@ Foam::scalar Foam::COxidationDiffusionLimitedRate<CloudType>::calculate
     dMassSRCarrier[O2GlobalId_] -= dmO2;
     dMassSRCarrier[CO2GlobalId_] += dmCO2;
 
+    const scalar HC = thermo.solids().properties()[CsLocalId_].H(T);
+    const scalar HCO2 = thermo.carrier().H(CO2GlobalId_, T);
+    const scalar HO2 = thermo.carrier().H(O2GlobalId_, T);
+
     // Heat of reaction [J]
-    return -HcCO2_*dmCO2;
+    return dmC*HC + dmO2*HO2 - dmCO2*HCO2;
 }
 
 

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

@@ -105,6 +105,21 @@ public:
             CloudType& owner
         );
 
+        //- Construct copy
+        COxidationDiffusionLimitedRate
+        (
+            const COxidationDiffusionLimitedRate<CloudType>& srm
+        );
+
+        //- Construct and return a clone
+        virtual autoPtr<SurfaceReactionModel<CloudType> > clone() const
+        {
+            return autoPtr<SurfaceReactionModel<CloudType> >
+            (
+                new COxidationDiffusionLimitedRate<CloudType>(*this)
+            );
+        }
+
 
     //- Destructor
     virtual ~COxidationDiffusionLimitedRate();
@@ -112,9 +127,6 @@ public:
 
     // Member Functions
 
-        //- Flag to indicate whether model activates surface reaction model
-        virtual bool active() const;
-
         //- Update surface reactions
         virtual scalar calculate
         (

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

@@ -42,10 +42,9 @@ COxidationKineticDiffusionLimitedRate
         owner,
         typeName
     ),
-    Sb_(dimensionedScalar(this->coeffDict().lookup("Sb")).value()),
+    C1_(readScalar(this->coeffDict().lookup("C1"))),
-    C1_(dimensionedScalar(this->coeffDict().lookup("C1")).value()),
+    C2_(readScalar(this->coeffDict().lookup("C2"))),
-    C2_(dimensionedScalar(this->coeffDict().lookup("C2")).value()),
+    E_(readScalar(this->coeffDict().lookup("E"))),
-    E_(dimensionedScalar(this->coeffDict().lookup("E")).value()),
     CsLocalId_(-1),
     O2GlobalId_(owner.composition().globalCarrierId("O2")),
     CO2GlobalId_(owner.composition().globalCarrierId("CO2")),
@@ -63,19 +62,29 @@ COxidationKineticDiffusionLimitedRate
     WC_ = WCO2 - WO2_;
     HcCO2_ = owner.thermo().carrier().Hc(CO2GlobalId_);
 
-    if (Sb_ < 0)
+    const scalar YCloc = owner.composition().Y0(idSolid)[CsLocalId_];
-    {
+    const scalar YSolidTot = owner.composition().YMixture0()[idSolid];
-        FatalErrorIn
+    Info<< "    C(s): particle mass fraction = " << YCloc*YSolidTot << endl;
+}
+
+
+template<class CloudType>
+Foam::COxidationKineticDiffusionLimitedRate<CloudType>::
+COxidationKineticDiffusionLimitedRate
 (
-            "COxidationKineticDiffusionLimitedRate"
+    const COxidationKineticDiffusionLimitedRate<CloudType>& srm
-            "("
+)
-                "const dictionary&, "
+:
-                "CloudType&"
+    SurfaceReactionModel<CloudType>(srm),
-            ")"
+    C1_(srm.C1_),
-        )   << "Stoichiometry of reaction, Sb, must be greater than zero" << nl
+    C2_(srm.C2_),
-            << exit(FatalError);
+    E_(srm.E_),
-    }
+    CsLocalId_(srm.CsLocalId_),
-}
+    O2GlobalId_(srm.O2GlobalId_),
+    CO2GlobalId_(srm.CO2GlobalId_),
+    WC_(srm.WC_),
+    WO2_(srm.WO2_)
+{}
 
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
@@ -88,13 +97,6 @@ Foam::COxidationKineticDiffusionLimitedRate<CloudType>::
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-template<class CloudType>
-bool Foam::COxidationKineticDiffusionLimitedRate<CloudType>::active() const
-{
-    return true;
-}
-
-
 template<class CloudType>
 Foam::scalar Foam::COxidationKineticDiffusionLimitedRate<CloudType>::calculate
 (
@@ -127,8 +129,10 @@ Foam::scalar Foam::COxidationKineticDiffusionLimitedRate<CloudType>::calculate
         return 0.0;
     }
 
+    const SLGThermo& thermo = this->owner().thermo();
+
     // Local mass fraction of O2 in the carrier phase
-    const scalar YO2 = this->owner().thermo().carrier().Y(O2GlobalId_)[cellI];
+    const scalar YO2 = thermo.carrier().Y(O2GlobalId_)[cellI];
 
     // Diffusion rate coefficient
     const scalar D0 = C1_/d*pow(0.5*(T + Tc), 0.75);
@@ -145,21 +149,28 @@ Foam::scalar Foam::COxidationKineticDiffusionLimitedRate<CloudType>::calculate
     // Limit mass transfer by availability of C
     dmC = min(mass*fComb, dmC);
 
+    // Molar consumption
+    const scalar dOmega = dmC/WC_;
+
     // Change in O2 mass [kg]
-    const scalar dmO2 = dmC/WC_*Sb_*WO2_;
+    const scalar dmO2 = dOmega*WO2_;
 
     // Mass of newly created CO2 [kg]
-    const scalar dmCO2 = dmC + dmO2;
+    const scalar dmCO2 = dOmega*(WC_ + WO2_);
 
     // Update local particle C mass
-    dMassSolid[CsLocalId_] += dmC;
+    dMassSolid[CsLocalId_] += dOmega*WC_;
 
     // Update carrier O2 and CO2 mass
     dMassSRCarrier[O2GlobalId_] -= dmO2;
     dMassSRCarrier[CO2GlobalId_] += dmCO2;
 
+    const scalar HC = thermo.solids().properties()[CsLocalId_].H(T);
+    const scalar HCO2 = thermo.carrier().H(CO2GlobalId_, T);
+    const scalar HO2 = thermo.carrier().H(O2GlobalId_, T);
+
     // Heat of reaction [J]
-    return -HcCO2_*dmCO2;
+    return dOmega*(WC_*HC + WO2_*HO2 - (WC_ + WO2_)*HCO2);
 }
 
 

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

@@ -61,9 +61,6 @@ class COxidationKineticDiffusionLimitedRate
 
         // Model constants
 
-            //- Stoichiometry of reaction
-            const scalar Sb_;
-
             //- Mass diffusion limited rate constant, C1
             const scalar C1_;
 
@@ -113,6 +110,21 @@ public:
             CloudType& owner
         );
 
+        //- Construct copy
+        COxidationKineticDiffusionLimitedRate
+        (
+            const COxidationKineticDiffusionLimitedRate<CloudType>& srm
+        );
+
+        //- Construct and return a clone
+        virtual autoPtr<SurfaceReactionModel<CloudType> > clone() const
+        {
+            return autoPtr<SurfaceReactionModel<CloudType> >
+            (
+                new COxidationKineticDiffusionLimitedRate<CloudType>(*this)
+            );
+        }
+
 
     //- Destructor
     virtual ~COxidationKineticDiffusionLimitedRate();
@@ -120,9 +132,6 @@ public:
 
     // Member Functions
 
-        //- Flag to indicate whether model activates surface reaction model
-        virtual bool active() const;
-
         //- Update surface reactions
         virtual scalar calculate
         (

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

@@ -50,14 +50,14 @@ Foam::COxidationMurphyShaddix<CloudType>::COxidationMurphyShaddix
         owner,
         typeName
     ),
-    D0_(dimensionedScalar(this->coeffDict().lookup("D0")).value()),
+    D0_(readScalar(this->coeffDict().lookup("D0"))),
-    rho0_(dimensionedScalar(this->coeffDict().lookup("rho0")).value()),
+    rho0_(readScalar(this->coeffDict().lookup("rho0"))),
-    T0_(dimensionedScalar(this->coeffDict().lookup("T0")).value()),
+    T0_(readScalar(this->coeffDict().lookup("T0"))),
-    Dn_(dimensionedScalar(this->coeffDict().lookup("Dn")).value()),
+    Dn_(readScalar(this->coeffDict().lookup("Dn"))),
-    A_(dimensionedScalar(this->coeffDict().lookup("A")).value()),
+    A_(readScalar(this->coeffDict().lookup("A"))),
-    E_(dimensionedScalar(this->coeffDict().lookup("E")).value()),
+    E_(readScalar(this->coeffDict().lookup("E"))),
-    n_(dimensionedScalar(this->coeffDict().lookup("n")).value()),
+    n_(readScalar(this->coeffDict().lookup("n"))),
-    WVol_(dimensionedScalar(this->coeffDict().lookup("WVol")).value()),
+    WVol_(readScalar(this->coeffDict().lookup("WVol"))),
     CsLocalId_(-1),
     O2GlobalId_(owner.composition().globalCarrierId("O2")),
     CO2GlobalId_(owner.composition().globalCarrierId("CO2")),
@@ -72,9 +72,36 @@ Foam::COxidationMurphyShaddix<CloudType>::COxidationMurphyShaddix
     WO2_ = owner.thermo().carrier().W(O2GlobalId_);
     const scalar WCO2 = owner.thermo().carrier().W(CO2GlobalId_);
     WC_ = WCO2 - WO2_;
+
+    const scalar YCloc = owner.composition().Y0(idSolid)[CsLocalId_];
+    const scalar YSolidTot = owner.composition().YMixture0()[idSolid];
+    Info<< "    C(s): particle mass fraction = " << YCloc*YSolidTot << endl;
 }
 
 
+template<class CloudType>
+Foam::COxidationMurphyShaddix<CloudType>::COxidationMurphyShaddix
+(
+    const COxidationMurphyShaddix<CloudType>& srm
+)
+:
+    SurfaceReactionModel<CloudType>(srm),
+    D0_(srm.D0_),
+    rho0_(srm.rho0_),
+    T0_(srm.T0_),
+    Dn_(srm.Dn_),
+    A_(srm.A_),
+    E_(srm.E_),
+    n_(srm.n_),
+    WVol_(srm.WVol_),
+    CsLocalId_(srm.CsLocalId_),
+    O2GlobalId_(srm.O2GlobalId_),
+    CO2GlobalId_(srm.CO2GlobalId_),
+    WC_(srm.WC_),
+    WO2_(srm.WO2_)
+{}
+
+
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 template<class CloudType>
@@ -84,13 +111,6 @@ Foam::COxidationMurphyShaddix<CloudType>::~COxidationMurphyShaddix()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-template<class CloudType>
-bool Foam::COxidationMurphyShaddix<CloudType>::active() const
-{
-    return true;
-}
-
-
 template<class CloudType>
 Foam::scalar Foam::COxidationMurphyShaddix<CloudType>::calculate
 (
@@ -123,9 +143,10 @@ Foam::scalar Foam::COxidationMurphyShaddix<CloudType>::calculate
         return 0.0;
     }
 
+    const SLGThermo& thermo = this->owner().thermo();
+
     // Cell carrier phase O2 species density [kg/m^3]
-    const scalar rhoO2 =
+    const scalar rhoO2 = rhoc*thermo.carrier().Y(O2GlobalId_)[cellI];
-        rhoc*this->owner().thermo().carrier().Y(O2GlobalId_)[cellI];
 
     if (rhoO2 < SMALL)
     {
@@ -204,11 +225,9 @@ Foam::scalar Foam::COxidationMurphyShaddix<CloudType>::calculate
     // Add to particle mass transfer
     dMassSolid[CsLocalId_] += dOmega*WC_;
 
-    const scalar HC =
+    const scalar HC = thermo.solids().properties()[CsLocalId_].H(T);
-        this->owner().composition().solids().properties()[CsLocalId_].Hf()
+    const scalar HCO2 = thermo.carrier().H(CO2GlobalId_, T);
-      + this->owner().composition().solids().properties()[CsLocalId_].cp()*T;
+    const scalar HO2 = thermo.carrier().H(O2GlobalId_, T);
-    const scalar HCO2 = this->owner().thermo().carrier().H(CO2GlobalId_, T);
-    const scalar HO2 = this->owner().thermo().carrier().H(O2GlobalId_, T);
 
     // Heat of reaction
     return dOmega*(WC_*HC + WO2_*HO2 - (WC_ + WO2_)*HCO2);

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

@@ -131,6 +131,21 @@ public:
             CloudType& owner
         );
 
+        //- Construct copy
+        COxidationMurphyShaddix
+        (
+            const COxidationMurphyShaddix<CloudType>& srm
+        );
+
+        //- Construct and return a clone
+        virtual autoPtr<SurfaceReactionModel<CloudType> > clone() const
+        {
+            return autoPtr<SurfaceReactionModel<CloudType> >
+            (
+                new COxidationMurphyShaddix<CloudType>(*this)
+            );
+        }
+
 
     //- Destructor
     virtual ~COxidationMurphyShaddix();
@@ -138,9 +153,6 @@ public:
 
     // Member Functions
 
-        //- Flag to indicate whether model activates surface reaction model
-        virtual bool active() const;
-
         //- Update surface reactions
         virtual scalar calculate
         (

src/lagrangian/dieselSpray/injector/commonRailInjector/commonRailInjector.C

@@ -25,7 +25,6 @@ License
 
 #include "commonRailInjector.H"
 #include "addToRunTimeSelectionTable.H"
-#include "Random.H"
 #include "mathematicalConstants.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -175,13 +174,13 @@ Foam::commonRailInjector::~commonRailInjector()
 
 void Foam::commonRailInjector::setTangentialVectors()
 {
-    Random rndGen(label(0));
+    cachedRandom rndGen(label(0), -1);
     scalar magV = 0.0;
     vector tangent;
 
     while (magV < SMALL)
     {
-        vector testThis = rndGen.vector01();
+        vector testThis = rndGen.sample01<vector>();
 
         tangent = testThis - (testThis & direction_)*direction_;
         magV = mag(tangent);
@@ -223,7 +222,7 @@ Foam::vector Foam::commonRailInjector::position
     const vector& axisOfSymmetry,
     const vector& axisOfWedge,
     const vector& axisOfWedgeNormal,
-    Random& rndGen
+    cachedRandom& rndGen
 ) const
 {
     if (twoD)
@@ -241,8 +240,8 @@ Foam::vector Foam::commonRailInjector::position
     else
     {
         // otherwise, disc injection
-        scalar iRadius = d_*rndGen.scalar01();
+        scalar iRadius = d_*rndGen.sample01<scalar>();
-        scalar iAngle = constant::mathematical::twoPi*rndGen.scalar01();
+        scalar iAngle = constant::mathematical::twoPi*rndGen.sample01<scalar>();
 
         return
         (

src/lagrangian/dieselSpray/injector/commonRailInjector/commonRailInjector.H

@@ -138,7 +138,7 @@ public:
             const vector& axisOfSymmetry,
             const vector& axisOfWedge,
             const vector& axisOfWedgeNormal,
-            Random& rndGen
+            cachedRandom& rndGen
         ) const;
 
         //- Return the number of holes

src/lagrangian/dieselSpray/injector/definedInjector/definedInjector.C

@@ -25,7 +25,6 @@ License
 
 #include "definedInjector.H"
 #include "addToRunTimeSelectionTable.H"
-#include "Random.H"
 #include "mathematicalConstants.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -158,13 +157,13 @@ Foam::definedInjector::~definedInjector()
 
 void Foam::definedInjector::setTangentialVectors()
 {
-    Random rndGen(label(0));
+    cachedRandom rndGen(label(0), -1);
     scalar magV = 0.0;
     vector tangent;
 
     while (magV < SMALL)
     {
-        vector testThis = rndGen.vector01();
+        vector testThis = rndGen.sample01<vector>();
 
         tangent = testThis - (testThis & direction_)*direction_;
         magV = mag(tangent);
@@ -204,7 +203,7 @@ Foam::vector Foam::definedInjector::position
     const vector& axisOfSymmetry,
     const vector& axisOfWedge,
     const vector& axisOfWedgeNormal,
-    Random& rndGen
+    cachedRandom& rndGen
 ) const
 {
     if (twoD)
@@ -222,8 +221,8 @@ Foam::vector Foam::definedInjector::position
     else
     {
         // otherwise, disc injection
-        scalar iRadius = d_*rndGen.scalar01();
+        scalar iRadius = d_*rndGen.sample01<scalar>();
-        scalar iAngle = constant::mathematical::twoPi*rndGen.scalar01();
+        scalar iAngle = constant::mathematical::twoPi*rndGen.sample01<scalar>();
 
         return
         (

src/lagrangian/dieselSpray/injector/definedInjector/definedInjector.H

@@ -139,7 +139,7 @@ public:
             const vector& axisOfSymmetry,
             const vector& axisOfWedge,
             const vector& axisOfWedgeNormal,
-            Random& rndGen
+            cachedRandom& rndGen
         ) const;
 
         //- Return the number of holes

src/lagrangian/dieselSpray/injector/injectorType/injectorType.H

@@ -43,7 +43,7 @@ SourceFiles
 #include "runTimeSelectionTables.H"
 #include "vector.H"
 #include "scalarField.H"
-#include "Random.H"
+#include "cachedRandom.H"
 #include "liquidMixture.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -120,7 +120,7 @@ public:
             const vector& axisOfSymmetry,
             const vector& axisOfWedge,
             const vector& axisOfWedgeNormal,
-            Random& rndGen
+            cachedRandom& rndGen
         ) const = 0;
 
         //- Return the number of holes

src/lagrangian/dieselSpray/injector/multiHoleInjector/multiHoleInjector.C

@@ -25,7 +25,6 @@ License
 
 #include "multiHoleInjector.H"
 #include "addToRunTimeSelectionTable.H"
-#include "Random.H"
 #include "unitConversion.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -194,7 +193,7 @@ void Foam::multiHoleInjector::setTangentialVectors()
         position_[i] = centerPosition_ + 0.5*nozzleTipDiameter_*dp;
     }
 
-    Random rndGen(label(0));
+    cachedRandom rndGen(label(0), -1);
 
     for (label i=0; i<nHoles_; i++)
     {
@@ -202,7 +201,7 @@ void Foam::multiHoleInjector::setTangentialVectors()
         scalar magV = 0;
         while (magV < SMALL)
         {
-            vector testThis = rndGen.vector01();
+            vector testThis = rndGen.sample01<vector>();
 
             tangent = testThis - (testThis & direction_[i])*direction_[i];
             magV = mag(tangent);
@@ -244,7 +243,7 @@ Foam::vector Foam::multiHoleInjector::position
     const vector& axisOfSymmetry,
     const vector& axisOfWedge,
     const vector& axisOfWedgeNormal,
-    Random& rndGen
+    cachedRandom& rndGen
 ) const
 {
     if (twoD)
@@ -262,8 +261,8 @@ Foam::vector Foam::multiHoleInjector::position
     else
     {
         // otherwise, disc injection
-        scalar iRadius = d_*rndGen.scalar01();
+        scalar iRadius = d_*rndGen.sample01<scalar>();
-        scalar iAngle = constant::mathematical::twoPi*rndGen.scalar01();
+        scalar iAngle = constant::mathematical::twoPi*rndGen.sample01<scalar>();
 
         return
         (

src/lagrangian/dieselSpray/injector/multiHoleInjector/multiHoleInjector.H

@@ -148,7 +148,7 @@ public:
             const vector& axisOfSymmetry,
             const vector& axisOfWedge,
             const vector& axisOfWedgeNormal,
-            Random& rndGen
+            cachedRandom& rndGen
         ) const;
 
         //- Return the number of holes

src/lagrangian/dieselSpray/injector/swirlInjector/swirlInjector.C

@@ -25,7 +25,6 @@ License
 
 #include "swirlInjector.H"
 #include "addToRunTimeSelectionTable.H"
-#include "Random.H"
 #include "mathematicalConstants.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -179,13 +178,13 @@ Foam::swirlInjector::~swirlInjector()
 
 void Foam::swirlInjector::setTangentialVectors()
 {
-    Random rndGen(label(0));
+    cachedRandom rndGen(label(0), -1);
     scalar magV = 0.0;
     vector tangent;
 
     while (magV < SMALL)
     {
-        vector testThis = rndGen.vector01();
+        vector testThis = rndGen.sample01<vector>();
 
         tangent = testThis - (testThis & direction_)*direction_;
         magV = mag(tangent);
@@ -226,7 +225,7 @@ Foam::vector Foam::swirlInjector::position
     const vector& axisOfSymmetry,
     const vector& axisOfWedge,
     const vector& axisOfWedgeNormal,
-    Random& rndGen
+    cachedRandom& rndGen
 ) const
 {
     if (twoD)
@@ -244,8 +243,8 @@ Foam::vector Foam::swirlInjector::position
     else
     {
         // otherwise, disc injection
-        scalar iRadius = d_*rndGen.scalar01();
+        scalar iRadius = d_*rndGen.sample01<scalar>();
-        scalar iAngle = constant::mathematical::twoPi*rndGen.scalar01();
+        scalar iAngle = constant::mathematical::twoPi*rndGen.sample01<scalar>();
 
         return
         (

src/lagrangian/dieselSpray/injector/swirlInjector/swirlInjector.H

@@ -157,7 +157,7 @@ public:
             const vector& axisOfSymmetry,
             const vector& axisOfWedge,
             const vector& axisOfWedgeNormal,
-            Random& rndGen
+            cachedRandom& rndGen
         ) const;
 
         //- Return the number of holes

src/lagrangian/dieselSpray/injector/unitInjector/unitInjector.C

@@ -25,7 +25,6 @@ License
 
 #include "unitInjector.H"
 #include "addToRunTimeSelectionTable.H"
-#include "Random.H"
 #include "mathematicalConstants.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -154,13 +153,13 @@ Foam::unitInjector::~unitInjector()
 
 void Foam::unitInjector::setTangentialVectors()
 {
-    Random rndGen(label(0));
+    cachedRandom rndGen(label(0), -1);
     scalar magV = 0.0;
     vector tangent;
 
     while (magV < SMALL)
     {
-        vector testThis = rndGen.vector01();
+        vector testThis = rndGen.sample01<vector>();
 
         tangent = testThis - (testThis & direction_)*direction_;
         magV = mag(tangent);
@@ -199,7 +198,7 @@ Foam::vector Foam::unitInjector::position
     const vector& axisOfSymmetry,
     const vector& axisOfWedge,
     const vector& axisOfWedgeNormal,
-    Random& rndGen
+    cachedRandom& rndGen
 ) const
 {
     if (twoD)
@@ -217,8 +216,8 @@ Foam::vector Foam::unitInjector::position
     else
     {
         // otherwise, disc injection
-        scalar iRadius = d_*rndGen.scalar01();
+        scalar iRadius = d_*rndGen.sample01<scalar>();
-        scalar iAngle = constant::mathematical::twoPi*rndGen.scalar01();
+        scalar iAngle = constant::mathematical::twoPi*rndGen.sample01<scalar>();
 
         return
         (

src/lagrangian/dieselSpray/injector/unitInjector/unitInjector.H

@@ -137,7 +137,7 @@ public:
             const vector& axisOfSymmetry,
             const vector& axisOfWedge,
             const vector& axisOfWedgeNormal,
-            Random& rndGen
+            cachedRandom& rndGen
         ) const;
 
         //- Return the number of holes

src/lagrangian/dieselSpray/parcel/parcel.C

@@ -94,14 +94,13 @@ Foam::parcel::parcel
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-bool Foam::parcel::move(spray& sDB)
+bool Foam::parcel::move(spray& sDB, const scalar trackTime)
 {
     const polyMesh& mesh = cloud().pMesh();
     const polyBoundaryMesh& pbMesh = mesh.boundaryMesh();
 
     const liquidMixture& fuels = sDB.fuels();
 
-    scalar deltaT = sDB.runTime().deltaTValue();
     label Nf = fuels.components().size();
     label Ns = sDB.composition().Y().size();
 
@@ -156,12 +155,12 @@ bool Foam::parcel::move(spray& sDB)
         pg,
         Yfg,
         m()*fuels.Y(X()),
-        deltaT
+        trackTime
     );
 
 
     // set the end-time for the track
-    scalar tEnd = (1.0 - stepFraction())*deltaT;
+    scalar tEnd = (1.0 - stepFraction())*trackTime;
 
     // set the maximum time step for this parcel
     scalar dtMax = min
@@ -215,7 +214,7 @@ bool Foam::parcel::move(spray& sDB)
             tEnd -= dt;
 
             // Set the current time-step fraction.
-            stepFraction() = 1.0 - tEnd/deltaT;
+            stepFraction() = 1.0 - tEnd/trackTime;
 
             if (onBoundary()) // hit face
             {
@@ -521,7 +520,7 @@ void Foam::parcel::updateParcelProperties
 
         scalar Taverage = TDroplet + (Tg - TDroplet)/3.0;
         // for a liquid Cl \approx Cp
-        scalar liquidcL = sDB.fuels().cp(pg, TDroplet, X());
+        scalar liquidcL = sDB.fuels().Cp(pg, TDroplet, X());
 
         cpMix = 0.0;
         for (label i=0; i<Ns; i++)

src/lagrangian/dieselSpray/parcel/parcel.H

@@ -349,7 +349,7 @@ public:
 
         // Parcel operations
 
-            bool move(spray& sprayData);
+            bool move(spray& sprayData, const scalar trackTime);
 
             //- Transform the position and physical properties of the particle
             //  according to the given transformation tensor

src/lagrangian/dieselSpray/parcel/setRelaxationTimes.C

@@ -65,7 +65,7 @@ void Foam::parcel::setRelaxationTimes
     // calculate mixture properties
     scalar W = 0.0;
     scalar kMixture = 0.0;
-    scalar cpMixture = 0.0;
+    scalar CpMixture = 0.0;
     scalar muf = 0.0;
 
     for (label i=0; i<Ns; i++)
@@ -74,7 +74,7 @@ void Foam::parcel::setRelaxationTimes
         W += Y/sDB.gasProperties()[i].W();
         // Using mass-fractions to average...
         kMixture += Y*sDB.gasProperties()[i].kappa(Tf);
-        cpMixture += Y*sDB.gasProperties()[i].Cp(Tf);
+        CpMixture += Y*sDB.gasProperties()[i].Cp(Tf);
         muf += Y*sDB.gasProperties()[i].mu(Tf);
     }
     W = 1.0/W;
@@ -98,7 +98,7 @@ void Foam::parcel::setRelaxationTimes
     scalar nuf = muf/rho;
 
     scalar liquidDensity = fuels.rho(pressure, T(), X());
-    scalar liquidcL = fuels.cp(pressure, T(), X());
+    scalar liquidcL = fuels.Cp(pressure, T(), X());
     scalar heatOfVapour = fuels.hl(pressure, T(), X());
 
     // calculate the partial rho of the fuel vapour
@@ -130,7 +130,7 @@ void Foam::parcel::setRelaxationTimes
     }
 
     scalar Reynolds = Re(Up, nuf);
-    scalar Prandtl = Pr(cpMixture, muf, kMixture);
+    scalar Prandtl = Pr(CpMixture, muf, kMixture);
 
     // calculate the characteritic times
 
@@ -280,7 +280,7 @@ void Foam::parcel::setRelaxationTimes
                 tauBoiling[i] = sDB.evaporation().boilingTime
                 (
                     fuels.properties()[i].rho(pressure, Td),
-                    fuels.properties()[i].cp(pressure, Td),
+                    fuels.properties()[i].Cp(pressure, Td),
                     heatOfVapour,
                     kMixture,
                     Nusselt,
@@ -292,7 +292,7 @@ void Foam::parcel::setRelaxationTimes
                     tBoilingSurface,
                     vapourSurfaceEnthalpy,
                     vapourFarEnthalpy,
-                    cpMixture,
+                    CpMixture,
                     temperature,
                     kLiquid
                 );

src/lagrangian/dieselSpray/spray/spray.C

@@ -68,7 +68,7 @@ Foam::spray::spray
     runTime_(U.time()),
     time0_(runTime_.value()),
     mesh_(U.mesh()),
-    rndGen_(label(0)),
+    rndGen_(label(0), -1),
     g_(g.value()),
 
     U_(U),

src/lagrangian/dieselSpray/spray/spray.H

@@ -39,7 +39,7 @@ Description
 #include "liquid.H"
 #include "autoPtr.H"
 #include "liquidMixture.H"
-#include "Random.H"
+#include "cachedRandom.H"
 #include "thermoPhysicsTypes.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -76,7 +76,7 @@ class spray
             const fvMesh& mesh_;
 
             //- Random number generator
-            Random rndGen_;
+            cachedRandom rndGen_;
 
         //- Acceleration due to gravity
         const vector& g_;
@@ -250,7 +250,7 @@ public:
             inline tmp<volScalarField> evaporationSource(const label i) const;
             inline tmp<volScalarField> heatTransferSource() const;
 
-            inline Random& rndGen();
+            inline cachedRandom& rndGen();
             inline label subCycles() const;
             inline const vector& g() const;
 

src/lagrangian/dieselSpray/spray/sprayI.H

@@ -221,7 +221,7 @@ inline tmp<volScalarField> spray::heatTransferSource() const
 }
 
 
-inline Random& spray::rndGen()
+inline cachedRandom& spray::rndGen()
 {
     return rndGen_;
 }

src/lagrangian/dieselSpray/spray/sprayInject.C

@@ -159,7 +159,8 @@ void Foam::spray::inject()
                             (runTime_.deltaTValue() - dt)
                            /runTime_.deltaTValue();
 
-                        bool keepParcel = pPtr->move(*this);
+                        bool keepParcel =
+                            pPtr->move(*this, runTime_.deltaTValue());
 
                         if (keepParcel)
                         {

src/lagrangian/dieselSpray/spray/sprayOps.C

@@ -76,7 +76,7 @@ void Foam::spray::move()
         srhos_[i] = 0.0;
     }
 
-    Cloud<parcel>::move(*this);
+    Cloud<parcel>::move(*this, runTime_.deltaTValue());
 }
 
 

src/lagrangian/dieselSpray/spraySubModels/atomizationModel/LISA/LISA.C

@@ -54,7 +54,6 @@ Foam::LISA::LISA
 :
     atomizationModel(dict, sm),
     coeffsDict_(dict.subDict(typeName + "Coeffs")),
-    rndGen_(sm.rndGen()),
     Cl_(readScalar(coeffsDict_.lookup("Cl"))),
     cTau_(readScalar(coeffsDict_.lookup("cTau"))),
     Q_(readScalar(coeffsDict_.lookup("Q"))),
@@ -361,8 +360,8 @@ void Foam::LISA::atomizeParcel
 
         do
         {
-            x = minValue + range*rndGen_.scalar01();
+            x = minValue + range*rndGen_.sample01<scalar>();
-            y = rndGen_.scalar01();
+            y = rndGen_.sample01<scalar>();
 
             scalar xx = pow(x/dD, nExp);
 

src/lagrangian/dieselSpray/spraySubModels/atomizationModel/LISA/LISA.H

@@ -68,7 +68,6 @@ private:
     // Private data
 
         dictionary coeffsDict_;
-        Random& rndGen_;
         scalar Cl_;
         scalar cTau_;
         scalar Q_;

src/lagrangian/dieselSpray/spraySubModels/atomizationModel/atomizationModel/atomizationModel.H

@@ -55,7 +55,7 @@ protected:
 
         const dictionary& dict_;
         spray& spray_;
-        Random& rndGen_;
+        cachedRandom& rndGen_;
 
 
 public:

src/lagrangian/dieselSpray/spraySubModels/atomizationModel/blobsSheetAtomization/blobsSheetAtomization.C

@@ -56,8 +56,7 @@ Foam::blobsSheetAtomization::blobsSheetAtomization
     atomizationModel(dict, sm),
     coeffsDict_(dict.subDict(typeName + "Coeffs")),
     B_(readScalar(coeffsDict_.lookup("B"))),
-    angle_(readScalar(coeffsDict_.lookup("angle"))),
+    angle_(readScalar(coeffsDict_.lookup("angle")))
-    rndGen_(sm.rndGen())
 {}
 
 

src/lagrangian/dieselSpray/spraySubModels/atomizationModel/blobsSheetAtomization/blobsSheetAtomization.H

@@ -67,7 +67,6 @@ private:
         dictionary coeffsDict_;
         scalar B_;
         scalar angle_;
-        Random& rndGen_;
 
 
 public:

src/lagrangian/dieselSpray/spraySubModels/breakupModel/SHF/SHF.C

@@ -53,7 +53,6 @@ Foam::SHF::SHF
     breakupModel(dict, sm),
     coeffsDict_(dict.subDict(typeName + "Coeffs")),
     g_(sm.g()),
-    rndGen_(sm.rndGen()),
     weCorrCoeff_(readScalar(coeffsDict_.lookup("weCorrCoeff"))),
     weBuCrit_(readScalar(coeffsDict_.lookup("weBuCrit"))),
     weBuBag_(readScalar(coeffsDict_.lookup("weBuBag"))),
@@ -186,9 +185,9 @@ void Foam::SHF::breakupParcel
 
             do
             {
-                x = cDmaxBM_*rndGen_.scalar01();
+                x = cDmaxBM_*rndGen_.sample01<scalar>();
                 d = sqr(x)*d05;
-                y = rndGen_.scalar01();
+                y = rndGen_.sample01<scalar>();
 
                 px =
                     x
@@ -217,9 +216,9 @@ void Foam::SHF::breakupParcel
             do
             {
 
-                x = cDmaxS_*rndGen_.scalar01();
+                x = cDmaxS_*rndGen_.sample01<scalar>();
                 d = sqr(x)*d05;
-                y = rndGen_.scalar01();
+                y = rndGen_.sample01<scalar>();
 
                 px =
                     x

src/lagrangian/dieselSpray/spraySubModels/breakupModel/SHF/SHF.H

@@ -66,8 +66,6 @@ private:
         // reference to gravity
         const vector& g_;
 
-        Random& rndGen_;
-
         // model constants
 
         scalar weCorrCoeff_;

src/lagrangian/dieselSpray/spraySubModels/breakupModel/TAB/TAB.C

@@ -175,7 +175,7 @@ void Foam::TAB::breakupParcel
 
                 label n = 0;
                 bool found = false;
-                scalar random = rndGen_.scalar01();
+                scalar random = rndGen_.sample01<scalar>();
                 while (!found && (n<99))
                 {
                     if (rrd_[n] > random)

src/lagrangian/dieselSpray/spraySubModels/breakupModel/breakupModel/breakupModel.H

@@ -56,7 +56,7 @@ protected:
         const dictionary& dict_;
 
         spray& spray_;
-        Random& rndGen_;
+        cachedRandom& rndGen_;
 
         Switch includeOscillation_;