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

applications/solvers/combustion/dieselEngineFoam/Make/options

@@ -16,6 +16,7 @@ EXE_INC = \
     -I$(LIB_SRC)/ODE/lnInclude \
     -I$(LIB_SRC)/engine/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/combustionModels/lnInclude
 
 EXE_LIBS = \
     -lengine \
@@ -35,4 +36,5 @@ EXE_LIBS = \
     -llaminarFlameSpeedModels \
     -lchemistryModel \
     -lODE \
-    -ldistributionModels
+    -ldistributionModels \
+    -lcombustionModels

applications/solvers/combustion/dieselEngineFoam/YEqn.H

@@ -10,7 +10,8 @@ tmp<fv::convectionScheme<scalar> > mvConvection
 );
 
 {
-
+    combustion->correct();
+    dQ = combustion->dQ();
     label inertIndex = -1;
     volScalarField Yt(0.0*Y[0]);
 
@@ -20,17 +21,19 @@ tmp<fv::convectionScheme<scalar> > mvConvection
         {
             volScalarField& Yi = Y[i];
 
-            solve
+            fvScalarMatrix YiEqn
             (
                 fvm::ddt(rho, Yi)
               + mvConvection->fvmDiv(phi, Yi)
               - fvm::laplacian(turbulence->muEff(), Yi)
               ==
                 dieselSpray.evaporationSource(i)
-              + kappa*chemistry.RR(i),
+              + combustion->R(Yi)
-                mesh.solver("Yi")
             );
 
+            YiEqn.relax();
+            YiEqn.solve(mesh.solver("Yi"));
+
             Yi.max(0.0);
             Yt += Yi;
         }

applications/solvers/combustion/dieselEngineFoam/createFields.H

@@ -1,10 +1,14 @@
-Info<< nl << "Reading thermophysicalProperties" << endl;
+Info<< "Creating combustion model\n" << endl;
 
-autoPtr<psiChemistryModel> pChemistry
+autoPtr<combustionModels::psiChemistryCombustionModel> combustion
 (
-    psiChemistryModel::New(mesh)
+    combustionModels::psiChemistryCombustionModel::New
+    (
+        mesh
+    )
 );
-psiChemistryModel& chemistry = pChemistry();
+
+psiChemistryModel& chemistry = combustion->pChemistry();
 
 hsCombustionThermo& thermo = chemistry.thermo();
 
@@ -55,20 +59,6 @@ volScalarField& hs = thermo.hs();
 
 #include "compressibleCreatePhi.H"
 
-volScalarField kappa
-(
-    IOobject
-    (
-        "kappa",
-        runTime.timeName(),
-        mesh,
-        IOobject::NO_READ,
-        IOobject::AUTO_WRITE
-    ),
-    mesh,
-    dimensionedScalar("zero", dimless, 0.0)
-);
-
 Info << "Creating turbulence model.\n" << nl;
 autoPtr<compressible::turbulenceModel> turbulence
 (
@@ -81,6 +71,9 @@ autoPtr<compressible::turbulenceModel> turbulence
     )
 );
 
+// Set the turbulence into the combustion model
+combustion->setTurbulence(turbulence());
+
 Info<< "Creating field DpDt\n" << endl;
 volScalarField DpDt
 (
@@ -96,16 +89,16 @@ forAll(Y, i)
 }
 fields.add(hs);
 
-DimensionedField<scalar, volMesh> chemistrySh
+volScalarField dQ
 (
     IOobject
     (
-        "chemistry::Sh",
+        "dQ",
         runTime.timeName(),
         mesh,
         IOobject::NO_READ,
-        IOobject::NO_WRITE
+        IOobject::AUTO_WRITE
     ),
     mesh,
-    dimensionedScalar("chemistrySh", dimEnergy/dimTime/dimVolume, 0.0)
+    dimensionedScalar("dQ", dimEnergy/dimTime, 0.0)
 );

applications/solvers/combustion/dieselEngineFoam/dieselEngineFoam.C

@@ -32,8 +32,8 @@ Description
 #include "fvCFD.H"
 #include "engineTime.H"
 #include "engineMesh.H"
-#include "hCombustionThermo.H"
 #include "turbulenceModel.H"
+#include "psiChemistryCombustionModel.H"
 #include "spray.H"
 #include "psiChemistryModel.H"
 #include "chemistrySolver.H"
@@ -54,7 +54,6 @@ int main(int argc, char *argv[])
     #include "createEngineMesh.H"
     #include "createFields.H"
     #include "readGravitationalAcceleration.H"
-    #include "readCombustionProperties.H"
     #include "createSpray.H"
     #include "initContinuityErrs.H"
     #include "readEngineTimeControls.H"
@@ -82,29 +81,6 @@ int main(int argc, char *argv[])
 
         dieselSpray.evolve();
 
-        Info<< "Solving chemistry" << endl;
-
-        chemistry.solve
-        (
-            runTime.value() - runTime.deltaTValue(),
-            runTime.deltaTValue()
-        );
-
-        // turbulent time scale
-        {
-            volScalarField tk
-            (
-                Cmix*sqrt(turbulence->muEff()/rho/turbulence->epsilon())
-            );
-            volScalarField tc(chemistry.tc());
-
-            // Chalmers PaSR model
-            kappa = (runTime.deltaT() + tc)/(runTime.deltaT() + tc + tk);
-        }
-
-        chemistrySh = kappa*chemistry.Sh()();
-
-
         #include "rhoEqn.H"
 
         for (pimple.start(); pimple.loop(); pimple++)
@@ -130,10 +106,7 @@ int main(int argc, char *argv[])
 
         rho = thermo.rho();
 
-        if (runTime.write())
+        runTime.write();
-        {
-            chemistry.dQ()().write();
-        }
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
             << "  ClockTime = " << runTime.elapsedClockTime() << " s"

applications/solvers/combustion/dieselEngineFoam/hsEqn.H

@@ -1,14 +1,21 @@
 {
-    solve
+    fvScalarMatrix hsEqn
     (
         fvm::ddt(rho, hs)
       + mvConvection->fvmDiv(phi, hs)
       - fvm::laplacian(turbulence->alphaEff(), hs)
      ==
         DpDt
-     + dieselSpray.heatTransferSource()().dimensionedInternalField()
+      + combustion->Sh()
-     + chemistrySh
+      + dieselSpray.heatTransferSource()()
+
     );
 
+    hsEqn.relax();
+    hsEqn.solve();
+
     thermo.correct();
+
+    Info<< "min/max(T) = "
+        << min(T).value() << ", " << max(T).value() << endl;
 }

applications/solvers/combustion/dieselEngineFoam/readCombustionProperties.H

@@ -1,18 +0,0 @@
-Info<< "Reading combustion properties\n" << endl;
-
-IOdictionary combustionProperties
-(
-    IOobject
-    (
-        "combustionProperties",
-        runTime.constant(),
-        mesh,
-        IOobject::MUST_READ_IF_MODIFIED,
-        IOobject::NO_WRITE
-    )
-);
-
-dimensionedScalar Cmix
-(
-    combustionProperties.lookup("Cmix")
-);

applications/solvers/combustion/dieselFoam/Make/options

@@ -14,7 +14,8 @@ EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/laminarFlameSpeed/lnInclude \
     -I$(LIB_SRC)/../applications/solvers/reactionThermo/XiFoam \
     -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
-    -I$(LIB_SRC)/ODE/lnInclude
+    -I$(LIB_SRC)/ODE/lnInclude \
+    -I$(LIB_SRC)/combustionModels/lnInclude
 
 EXE_LIBS = \
     -lcompressibleTurbulenceModel \
@@ -33,4 +34,5 @@ EXE_LIBS = \
     -lchemistryModel \
     -lODE \
     -ldistributionModels \
-    -lfiniteVolume
+    -lfiniteVolume \
+    -lcombustionModels

applications/solvers/combustion/dieselFoam/dieselFoam.C

@@ -30,12 +30,11 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "hCombustionThermo.H"
+#include "psiChemistryCombustionModel.H"
 #include "turbulenceModel.H"
-#include "spray.H"
 #include "psiChemistryModel.H"
 #include "chemistrySolver.H"
-
+#include "spray.H"
 #include "multivariateScheme.H"
 #include "IFstream.H"
 #include "OFstream.H"
@@ -52,7 +51,6 @@ int main(int argc, char *argv[])
     #include "createMesh.H"
     #include "createFields.H"
     #include "readGravitationalAcceleration.H"
-    #include "readCombustionProperties.H"
     #include "createSpray.H"
     #include "initContinuityErrs.H"
     #include "readTimeControls.H"
@@ -79,26 +77,6 @@ int main(int argc, char *argv[])
 
         Info<< "Solving chemistry" << endl;
 
-        chemistry.solve
-        (
-            runTime.value() - runTime.deltaTValue(),
-            runTime.deltaTValue()
-        );
-
-        // turbulent time scale
-        {
-            volScalarField tk
-            (
-                Cmix*sqrt(turbulence->muEff()/rho/turbulence->epsilon())
-            );
-            volScalarField tc(chemistry.tc());
-
-            // Chalmers PaSR model
-            kappa = (runTime.deltaT() + tc)/(runTime.deltaT()+tc+tk);
-        }
-
-        chemistrySh = kappa*chemistry.Sh()();
-
         #include "rhoEqn.H"
 
         // --- Pressure-velocity PIMPLE corrector loop
@@ -124,10 +102,7 @@ int main(int argc, char *argv[])
 
         rho = thermo.rho();
 
-        if (runTime.write())
+        runTime.write();
-        {
-            chemistry.dQ()().write();
-        }
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
             << "  ClockTime = " << runTime.elapsedClockTime() << " s"

applications/solvers/combustion/fireFoam/YhsEqn.H

@@ -20,7 +20,6 @@ tmp<fv::convectionScheme<scalar> > mvConvection
         if (Y[i].name() != inertSpecie)
         {
             volScalarField& Yi = Y[i];
-            fvScalarMatrix R(combustion->R(Yi));
 
             fvScalarMatrix YiEqn
             (
@@ -30,7 +29,7 @@ tmp<fv::convectionScheme<scalar> > mvConvection
               ==
                 parcels.SYi(i, Yi)
               + surfaceFilm.Srho(i)
-              + R
+              + combustion->R(Yi)
             );
 
             YiEqn.relax();
@@ -55,7 +54,7 @@ tmp<fv::convectionScheme<scalar> > mvConvection
       - fvm::laplacian(turbulence->alphaEff(), hs)
      ==
         DpDt
-      + dQ
+      + combustion->Sh()
       + radiation->Shs(thermo)
       + parcels.Sh(hs)
       + surfaceFilm.Sh()
@@ -66,5 +65,6 @@ tmp<fv::convectionScheme<scalar> > mvConvection
 
     thermo.correct();
 
-    Info<< "min/max(T) = " << min(T).value() << ", " << max(T).value() << endl;
+    Info<< "min/max(T) = "
+        << min(T).value() << ", " << max(T).value() << endl;
 }

applications/solvers/combustion/fireFoam/createFields.H

@@ -1,10 +1,16 @@
+    Info<< "Creating combustion model\n" << endl;
+
+    autoPtr<combustionModels::psiCombustionModel> combustion
+    (
+        combustionModels::psiCombustionModel::New
+        (
+            mesh
+        )
+    );
+
     Info<< "Reading thermophysical properties\n" << endl;
 
-    autoPtr<hsCombustionThermo> pThermo
+    hsCombustionThermo& thermo = combustion->thermo();
-    (
-        hsCombustionThermo::New(mesh)
-    );
-    hsCombustionThermo& thermo = pThermo();
 
     SLGThermo slgThermo(mesh, thermo);
 
@@ -60,30 +66,8 @@
         )
     );
 
-    IOdictionary combustionProperties
+    // Set the turbulence into the combustion model
-    (
+    combustion->setTurbulence(turbulence());
-        IOobject
-        (
-            "combustionProperties",
-            runTime.constant(),
-            mesh,
-            IOobject::MUST_READ_IF_MODIFIED,
-            IOobject::NO_WRITE
-        )
-    );
-
-    Info<< "Creating combustion model\n" << endl;
-    autoPtr<combustionModel> combustion
-    (
-        combustionModel::combustionModel::New
-        (
-            combustionProperties,
-            thermo,
-            turbulence(),
-            phi,
-            rho
-        )
-    );
 
     volScalarField dQ
     (
@@ -96,7 +80,7 @@
             IOobject::AUTO_WRITE
         ),
         mesh,
-        dimensionedScalar("dQ", dimMass/pow3(dimTime)/dimLength, 0.0)
+        dimensionedScalar("dQ", dimEnergy/dimTime, 0.0)
     );
 
     Info<< "Creating field DpDt\n" << endl;

applications/solvers/combustion/fireFoam/fireFoam.C

@@ -38,9 +38,8 @@ Description
 #include "pyrolysisModel.H"
 #include "radiationModel.H"
 #include "SLGThermo.H"
-#include "hsCombustionThermo.H"
 #include "solidChemistryModel.H"
-#include "combustionModel.H"
+#include "psiCombustionModel.H"
 #include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -51,7 +50,6 @@ int main(int argc, char *argv[])
 
     #include "createTime.H"
     #include "createMesh.H"
-    #include "readChemistryProperties.H"
     #include "readGravitationalAcceleration.H"
     #include "createFields.H"
     #include "createClouds.H"

applications/solvers/combustion/fireFoam/readChemistryProperties.H

@@ -1,23 +0,0 @@
-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/combustion/reactingFoam/Make/options

@@ -5,7 +5,8 @@ EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
     -I$(LIB_SRC)/ODE/lnInclude \
-    -I$(LIB_SRC)/finiteVolume/lnInclude
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/combustionModels/lnInclude
 
 EXE_LIBS = \
     -lcompressibleTurbulenceModel \
@@ -16,4 +17,5 @@ EXE_LIBS = \
     -lbasicThermophysicalModels \
     -lchemistryModel \
     -lODE \
-    -lfiniteVolume
+    -lfiniteVolume \
+    -lcombustionModels

applications/solvers/combustion/reactingFoam/YEqn.H

@@ -10,6 +10,8 @@ tmp<fv::convectionScheme<scalar> > mvConvection
 );
 
 {
+    combustion->correct();
+    dQ = combustion->dQ();
     label inertIndex = -1;
     volScalarField Yt(0.0*Y[0]);
 
@@ -19,16 +21,18 @@ tmp<fv::convectionScheme<scalar> > mvConvection
         {
             volScalarField& Yi = Y[i];
 
-            solve
+            fvScalarMatrix YiEqn
             (
                 fvm::ddt(rho, Yi)
               + mvConvection->fvmDiv(phi, Yi)
               - fvm::laplacian(turbulence->muEff(), Yi)
              ==
-                kappa*chemistry.RR(i),
+                combustion->R(Yi)
-                mesh.solver("Yi")
             );
 
+            YiEqn.relax();
+            YiEqn.solve(mesh.solver("Yi"));
+
             Yi.max(0.0);
             Yt += Yi;
         }

applications/solvers/combustion/reactingFoam/chemistry.H

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

applications/solvers/combustion/reactingFoam/createFields.H

@@ -1,9 +1,14 @@
-Info<< nl << "Reading thermophysicalProperties" << endl;
+Info<< "Creating combustion model\n" << endl;
-autoPtr<psiChemistryModel> pChemistry
+
+autoPtr<combustionModels::psiChemistryCombustionModel> combustion
 (
-    psiChemistryModel::New(mesh)
+    combustionModels::psiChemistryCombustionModel::New
+    (
+        mesh
+    )
 );
-psiChemistryModel& chemistry = pChemistry();
+
+psiChemistryModel& chemistry = combustion->pChemistry();
 
 hsCombustionThermo& thermo = chemistry.thermo();
 
@@ -45,20 +50,6 @@ const volScalarField& T = thermo.T();
 
 #include "compressibleCreatePhi.H"
 
-volScalarField kappa
-(
-    IOobject
-    (
-        "kappa",
-        runTime.timeName(),
-        mesh,
-        IOobject::NO_READ,
-        IOobject::AUTO_WRITE
-    ),
-    mesh,
-    dimensionedScalar("zero", dimless, 0.0)
-);
-
 Info << "Creating turbulence model.\n" << nl;
 autoPtr<compressible::turbulenceModel> turbulence
 (
@@ -71,6 +62,9 @@ autoPtr<compressible::turbulenceModel> turbulence
     )
 );
 
+// Set the turbulence into the combustion model
+combustion->setTurbulence(turbulence());
+
 Info<< "Creating field DpDt\n" << endl;
 volScalarField DpDt
 (
@@ -85,16 +79,16 @@ forAll(Y, i)
 }
 fields.add(hs);
 
-DimensionedField<scalar, volMesh> chemistrySh
+volScalarField dQ
 (
     IOobject
     (
-        "chemistry::Sh",
+        "dQ",
         runTime.timeName(),
         mesh,
         IOobject::NO_READ,
-        IOobject::NO_WRITE
+        IOobject::AUTO_WRITE
     ),
     mesh,
-    dimensionedScalar("chemistrySh", dimEnergy/dimTime/dimVolume, 0.0)
+    dimensionedScalar("dQ", dimEnergy/dimTime, 0.0)
 );

applications/solvers/combustion/reactingFoam/hsEqn.H

@@ -7,7 +7,7 @@
 //      - fvm::laplacian(turbulence->muEff(), hs)  // unit lewis no.
      ==
         DpDt
-      + chemistrySh
+      + combustion->Sh()
     );
 
     hsEqn.relax();

applications/solvers/combustion/reactingFoam/reactingFoam.C

@@ -30,10 +30,8 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "hCombustionThermo.H"
 #include "turbulenceModel.H"
-#include "psiChemistryModel.H"
+#include "psiChemistryCombustionModel.H"
-#include "chemistrySolver.H"
 #include "multivariateScheme.H"
 #include "pimpleControl.H"
 
@@ -44,7 +42,6 @@ int main(int argc, char *argv[])
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createMesh.H"
-    #include "readChemistryProperties.H"
     #include "readGravitationalAcceleration.H"
     #include "createFields.H"
     #include "initContinuityErrs.H"
@@ -67,7 +64,6 @@ int main(int argc, char *argv[])
         runTime++;
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        #include "chemistry.H"
         #include "rhoEqn.H"
 
         for (pimple.start(); pimple.loop(); pimple++)
@@ -88,11 +84,6 @@ int main(int argc, char *argv[])
             }
         }
 
-        if (runTime.write())
-        {
-            chemistry.dQ()().write();
-        }
-
         runTime.write();
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"

applications/solvers/combustion/reactingFoam/readChemistryProperties.H

@@ -1,23 +0,0 @@
-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/combustion/rhoReactingFoam/Make/options

@@ -6,7 +6,8 @@ EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
     -I$(LIB_SRC)/ODE/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(FOAM_SOLVERS)/combustion/reactingFoam
+    -I$(FOAM_SOLVERS)/combustion/reactingFoam \
+    -I$(LIB_SRC)/combustionModels/lnInclude
 
 
 EXE_LIBS = \
@@ -18,4 +19,5 @@ EXE_LIBS = \
     -lbasicThermophysicalModels \
     -lchemistryModel \
     -lODE \
-    -lfiniteVolume
+    -lfiniteVolume \
+    -lcombustionModels

applications/solvers/combustion/rhoReactingFoam/YEqn.H

@@ -10,6 +10,8 @@ tmp<fv::convectionScheme<scalar> > mvConvection
 );
 
 {
+    combustion->correct();
+    dQ = combustion->dQ();
     label inertIndex = -1;
     volScalarField Yt(0.0*Y[0]);
 
@@ -19,14 +21,13 @@ tmp<fv::convectionScheme<scalar> > mvConvection
         {
             volScalarField& Yi = Y[i];
 
-            solve
+            fvScalarMatrix YiEqn
             (
                 fvm::ddt(rho, Yi)
               + mvConvection->fvmDiv(phi, Yi)
               - fvm::laplacian(turbulence->muEff(), Yi)
              ==
-                kappa*chemistry.RR(i),
+                combustion->R(Yi)
-                mesh.solver("Yi")
             );
 
             Yi.max(0.0);

applications/solvers/combustion/rhoReactingFoam/createFields.H

@@ -1,9 +1,14 @@
-Info<< nl << "Reading thermophysicalProperties" << endl;
+Info<< "Creating combustion model\n" << endl;
-autoPtr<rhoChemistryModel> pChemistry
+
+autoPtr<combustionModels::rhoChemistryCombustionModel> combustion
 (
-    rhoChemistryModel::New(mesh)
+    combustionModels::rhoChemistryCombustionModel::New
+    (
+        mesh
+    )
 );
-rhoChemistryModel& chemistry = pChemistry();
+
+rhoChemistryModel& chemistry = combustion->pChemistry();
 
 hsReactionThermo& thermo = chemistry.thermo();
 
@@ -46,19 +51,6 @@ const volScalarField& T = thermo.T();
 
 #include "compressibleCreatePhi.H"
 
-volScalarField kappa
-(
-    IOobject
-    (
-        "kappa",
-        runTime.timeName(),
-        mesh,
-        IOobject::NO_READ,
-        IOobject::AUTO_WRITE
-    ),
-    mesh,
-    dimensionedScalar("zero", dimless, 0.0)
-);
 
 Info << "Creating turbulence model.\n" << nl;
 autoPtr<compressible::turbulenceModel> turbulence
@@ -72,6 +64,9 @@ autoPtr<compressible::turbulenceModel> turbulence
     )
 );
 
+// Set the turbulence into the combustion model
+combustion->setTurbulence(turbulence());
+
 Info<< "Creating field DpDt\n" << endl;
 volScalarField DpDt
 (
@@ -86,16 +81,16 @@ forAll(Y, i)
 }
 fields.add(hs);
 
-DimensionedField<scalar, volMesh> chemistrySh
+volScalarField dQ
 (
     IOobject
     (
-        "chemistry::Sh",
+        "dQ",
         runTime.timeName(),
         mesh,
         IOobject::NO_READ,
-        IOobject::NO_WRITE
+        IOobject::AUTO_WRITE
     ),
     mesh,
-    dimensionedScalar("chemistrySh", dimEnergy/dimTime/dimVolume, 0.0)
+    dimensionedScalar("dQ", dimEnergy/dimTime, 0.0)
 );

applications/solvers/combustion/rhoReactingFoam/hsEqn.H

@@ -4,9 +4,10 @@
         fvm::ddt(rho, hs)
       + mvConvection->fvmDiv(phi, hs)
       - fvm::laplacian(turbulence->alphaEff(), hs)
+//    - fvm::laplacian(turbulence->muEff(), hs)  // unit lewis no.
      ==
         DpDt
-      + chemistrySh
+      + combustion->Sh()
     );
 
     hsEqn.relax();
@@ -14,6 +15,6 @@
 
     thermo.correct();
 
-    Info<< "T gas min/max   = " << min(T).value() << ", "
+   Info<< "min/max(T) = "
-        << max(T).value() << endl;
+        << min(T).value() << ", " << max(T).value() << endl;
 }

applications/solvers/combustion/rhoReactingFoam/readChemistryProperties.H

@@ -1,23 +0,0 @@
-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/combustion/rhoReactingFoam/rhoReactingFoam.C

@@ -31,10 +31,8 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "hReactionThermo.H"
+#include "rhoChemistryCombustionModel.H"
 #include "turbulenceModel.H"
-#include "rhoChemistryModel.H"
-#include "chemistrySolver.H"
 #include "multivariateScheme.H"
 #include "pimpleControl.H"
 
@@ -45,7 +43,6 @@ int main(int argc, char *argv[])
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createMesh.H"
-    #include "readChemistryProperties.H"
     #include "readGravitationalAcceleration.H"
     #include "createFields.H"
     #include "initContinuityErrs.H"
@@ -68,7 +65,6 @@ int main(int argc, char *argv[])
         runTime++;
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        #include "chemistry.H"
         #include "rhoEqn.H"
 
         // --- Pressure-velocity PIMPLE corrector loop
@@ -92,10 +88,7 @@ int main(int argc, char *argv[])
 
         rho = thermo.rho();
 
-        if (runTime.write())
+        runTime.write();
-        {
-            chemistry.dQ()().write();
-        }
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
             << "  ClockTime = " << runTime.elapsedClockTime() << " s"

applications/solvers/lagrangian/LTSReactingParcelFoam/LTSReactingParcelFoam.C

@@ -36,11 +36,9 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "hReactionThermo.H"
 #include "turbulenceModel.H"
 #include "basicReactingMultiphaseCloud.H"
-#include "rhoChemistryModel.H"
+#include "rhoChemistryCombustionModel.H"
-#include "chemistrySolver.H"
 #include "radiationModel.H"
 #include "porousZones.H"
 #include "timeActivatedExplicitSource.H"
@@ -75,7 +73,6 @@ int main(int argc, char *argv[])
 
     while (runTime.run())
     {
-        #include "readChemistryProperties.H"
         #include "readAdditionalSolutionControls.H"
         #include "readTimeControls.H"
 
@@ -85,7 +82,6 @@ int main(int argc, char *argv[])
 
         parcels.evolve();
 
-        #include "chemistry.H"
         #include "timeScales.H"
 
         #include "rhoEqn.H"
@@ -111,10 +107,7 @@ int main(int argc, char *argv[])
             }
         }
 
-        if (runTime.write())
+        runTime.write();
-        {
-            chemistry.dQ()().write();
-        }
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
             << "  ClockTime = " << runTime.elapsedClockTime() << " s"

applications/solvers/lagrangian/LTSReactingParcelFoam/Make/options

@@ -20,7 +20,8 @@ EXE_INC = \
     -I$(LIB_SRC)/ODE/lnInclude \
     -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
     -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
-    -I$(LIB_SRC)/sampling/lnInclude
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/combustionModels/lnInclude
 
 EXE_LIBS = \
     -lfiniteVolume \
@@ -44,4 +45,5 @@ EXE_LIBS = \
     -lODE \
     -lregionModels \
     -lsurfaceFilmModels \
-    -lsampling
+    -lsampling \
+    -lcombustionModels

applications/solvers/lagrangian/LTSReactingParcelFoam/YEqn.H

@@ -9,6 +9,9 @@ tmp<fv::convectionScheme<scalar> > mvConvection
     )
 );
 
+combustion->correct();
+dQ = combustion->dQ();
+
 if (solveSpecies)
 {
     label inertIndex = -1;
@@ -19,6 +22,7 @@ if (solveSpecies)
         if (Y[i].name() != inertSpecie)
         {
             volScalarField& Yi = Y[i];
+
             solve
             (
                 fvm::ddt(rho, Yi)
@@ -26,7 +30,7 @@ if (solveSpecies)
               - fvm::laplacian(turbulence->muEff(), Yi)
              ==
                 parcels.SYi(i, Yi)
-              + kappa*chemistry.RR(i)().dimensionedInternalField()
+              + combustion->R(Yi)
               + massSource.Su(i),
                 mesh.solver("Yi")
             );

applications/solvers/lagrangian/LTSReactingParcelFoam/createFields.H

@@ -1,10 +1,14 @@
-    Info<< "Reading thermophysical properties\n" << endl;
+    Info<< "Creating combustion model\n" << endl;
 
-    autoPtr<rhoChemistryModel> pChemistry
+    autoPtr<combustionModels::rhoChemistryCombustionModel> combustion
     (
-        rhoChemistryModel::New(mesh)
+        combustionModels::rhoChemistryCombustionModel::New
+        (
+            mesh
+        )
     );
-    rhoChemistryModel& chemistry = pChemistry();
+
+    rhoChemistryModel& chemistry = combustion->pChemistry();
 
     hsReactionThermo& thermo = chemistry.thermo();
 
@@ -57,20 +61,6 @@
 
     #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("PIMPLE").lookup("rhoMax")
@@ -93,6 +83,9 @@
         )
     );
 
+    // Set the turbulence into the combustion model
+    combustion->setTurbulence(turbulence());
+
     Info<< "Creating multi-variate interpolation scheme\n" << endl;
     multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
 
@@ -102,20 +95,21 @@
     }
     fields.add(hs);
 
-    DimensionedField<scalar, volMesh> chemistrySh
+    volScalarField dQ
     (
         IOobject
         (
-            "chemistry::Sh",
+            "dQ",
             runTime.timeName(),
             mesh,
             IOobject::NO_READ,
-            IOobject::NO_WRITE
+            IOobject::AUTO_WRITE
         ),
         mesh,
-        dimensionedScalar("chemistrySh", dimEnergy/dimTime/dimVolume, 0.0)
+        dimensionedScalar("dQ", dimEnergy/dimTime, 0.0)
     );
 
+
     volScalarField rDeltaT
     (
         IOobject

applications/solvers/lagrangian/LTSReactingParcelFoam/hsEqn.H

@@ -9,7 +9,7 @@
       + parcels.Sh(hs)
       + radiation->Shs(thermo)
       + energySource.Su()
-      + chemistrySh
+      + combustion->Sh()
     );
 
     hsEqn.solve();

applications/solvers/lagrangian/LTSReactingParcelFoam/readChemistryProperties.H

@@ -1,23 +0,0 @@
-// 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/LTSReactingParcelFoam/timeScales.H

@@ -77,7 +77,7 @@ Info<< "Time scales min/max:" << endl;
                 DpDt
               + parcels.hsTrans()/(mesh.V()*runTime.deltaT())
               + energySource.Su()
-              + chemistrySh
+              + combustion->Sh()()
             )
            /rho
         );

applications/solvers/lagrangian/coalChemistryFoam/Make/options

@@ -21,6 +21,7 @@ EXE_INC = \
     -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
     -I$(LIB_SRC)/ODE/lnInclude \
     -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/combustionModels/lnInclude \
     -I$(FOAM_SOLVERS)/combustion/reactingFoam
 
 
@@ -47,4 +48,5 @@ EXE_LIBS = \
     -lregionModels \
     -lsurfaceFilmModels \
     -lODE \
-    -lsampling
+    -lsampling \
+    -lcombustionModels

applications/solvers/lagrangian/coalChemistryFoam/YEqn.H

@@ -11,6 +11,8 @@ tmp<fv::convectionScheme<scalar> > mvConvection
 
 
 {
+    combustion->correct();
+    dQ = combustion->dQ();
     label inertIndex = -1;
     volScalarField Yt(0.0*Y[0]);
 
@@ -19,16 +21,20 @@ tmp<fv::convectionScheme<scalar> > mvConvection
         if (Y[i].name() != inertSpecie)
         {
             volScalarField& Yi = Y[i];
-            solve
+
+            fvScalarMatrix YiEqn
             (
                 fvm::ddt(rho, Yi)
               + mvConvection->fvmDiv(phi, Yi)
               - fvm::laplacian(turbulence->muEff(), Yi)
               ==
                 coalParcels.SYi(i, Yi)
-              + kappa*chemistry.RR(i)().dimensionedInternalField()
+              + combustion->R(Yi)
             );
 
+            YiEqn.relax();
+            YiEqn.solve(mesh.solver("Yi"));
+
             Yi.max(0.0);
             Yt += Yi;
         }

applications/solvers/lagrangian/coalChemistryFoam/coalChemistryFoam.C

@@ -36,12 +36,10 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "hCombustionThermo.H"
 #include "turbulenceModel.H"
 #include "basicThermoCloud.H"
 #include "coalCloud.H"
-#include "psiChemistryModel.H"
+#include "psiChemistryCombustionModel.H"
-#include "chemistrySolver.H"
 #include "timeActivatedExplicitSource.H"
 #include "radiationModel.H"
 #include "SLGThermo.H"
@@ -55,7 +53,6 @@ int main(int argc, char *argv[])
 
     #include "createTime.H"
     #include "createMesh.H"
-    #include "readChemistryProperties.H"
     #include "readGravitationalAcceleration.H"
     #include "createFields.H"
     #include "createClouds.H"
@@ -88,7 +85,6 @@ int main(int argc, char *argv[])
 
         limestoneParcels.evolve();
 
-        #include "chemistry.H"
         #include "rhoEqn.H"
 
         // --- Pressure-velocity PIMPLE corrector loop
@@ -112,10 +108,7 @@ int main(int argc, char *argv[])
 
         rho = thermo.rho();
 
-        if (runTime.write())
+        runTime.write();
-        {
-            chemistry.dQ()().write();
-        }
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
             << "  ClockTime = " << runTime.elapsedClockTime() << " s"

applications/solvers/lagrangian/coalChemistryFoam/createFields.H

@@ -1,10 +1,14 @@
-    Info<< "Reading thermophysical properties\n" << endl;
+    Info<< "Creating combustion model\n" << endl;
 
-    autoPtr<psiChemistryModel> pChemistry
+    autoPtr<combustionModels::psiChemistryCombustionModel> combustion
     (
-        psiChemistryModel::New(mesh)
+        combustionModels::psiChemistryCombustionModel::New
+        (
+            mesh
+        )
     );
-    psiChemistryModel& chemistry = pChemistry();
+
+    psiChemistryModel& chemistry = combustion->pChemistry();
 
     hsCombustionThermo& thermo = chemistry.thermo();
 
@@ -96,20 +100,6 @@
 
     #include "compressibleCreatePhi.H"
 
-    DimensionedField<scalar, volMesh> kappa
-    (
-        IOobject
-        (
-            "kappa",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh,
-        dimensionedScalar("zero", dimless, 0.0)
-    );
-
     Info<< "Creating turbulence model\n" << endl;
     autoPtr<compressible::turbulenceModel> turbulence
     (
@@ -122,6 +112,9 @@
         )
     );
 
+    // Set the turbulence into the combustion model
+    combustion->setTurbulence(turbulence());
+
     Info<< "Creating field DpDt\n" << endl;
     volScalarField DpDt
     (
@@ -138,16 +131,16 @@
         "hs"
     );
 
-    DimensionedField<scalar, volMesh> chemistrySh
+    volScalarField dQ
     (
         IOobject
         (
-            "chemistry::Sh",
+            "dQ",
             runTime.timeName(),
             mesh,
             IOobject::NO_READ,
-            IOobject::NO_WRITE
+            IOobject::AUTO_WRITE
         ),
         mesh,
-        dimensionedScalar("chemistrySh", dimEnergy/dimTime/dimVolume, 0.0)
+        dimensionedScalar("dQ", dimEnergy/dimTime, 0.0)
     );

applications/solvers/lagrangian/coalChemistryFoam/hsEqn.H

@@ -6,15 +6,14 @@
       - fvm::laplacian(turbulence->alphaEff(), hs)
      ==
         DpDt
+     +  combustion->Sh()
      +  coalParcels.Sh(hs)
      +  limestoneParcels.Sh(hs)
      +  enthalpySource.Su()
      +  radiation->Shs(thermo)
-     +  chemistrySh
     );
 
     hsEqn.relax();
-
     hsEqn.solve();
 
     thermo.correct();

applications/solvers/lagrangian/coalChemistryFoam/readChemistryProperties.H

@@ -1,23 +0,0 @@
-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/porousExplicitSourceReactingParcelFoam/Make/options

@@ -21,6 +21,7 @@ EXE_INC = \
     -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
     -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
     -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/combustionModels/lnInclude \
     -I$(FOAM_SOLVERS)/combustion/reactingFoam
 
 
@@ -46,4 +47,5 @@ EXE_LIBS = \
     -lODE \
     -lregionModels \
     -lsurfaceFilmModels \
-    -lsampling
+    -lsampling \
+    -lcombustionModels

applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/YEqn.H

@@ -10,6 +10,8 @@ tmp<fv::convectionScheme<scalar> > mvConvection
     )
 );
 
+combustion->correct();
+dQ = combustion->dQ();
 
 if (solveSpecies)
 {
@@ -21,6 +23,7 @@ if (solveSpecies)
         if (Y[i].name() != inertSpecie)
         {
             volScalarField& Yi = Y[i];
+
             solve
             (
                 fvm::ddt(rho, Yi)
@@ -28,7 +31,7 @@ if (solveSpecies)
               - fvm::laplacian(turbulence->muEff(), Yi)
               ==
                 parcels.SYi(i, Yi)
-              + kappa*chemistry.RR(i)().dimensionedInternalField()
+              + combustion->R(Yi)
               + massSource.Su(i),
                 mesh.solver("Yi")
             );

applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/createFields.H

@@ -1,10 +1,14 @@
-    Info<< "Reading thermophysical properties\n" << endl;
+    Info<< "Creating combustion model\n" << endl;
 
-    autoPtr<rhoChemistryModel> pChemistry
+    autoPtr<combustionModels::rhoChemistryCombustionModel> combustion
     (
-        rhoChemistryModel::New(mesh)
+        combustionModels::rhoChemistryCombustionModel::New
+        (
+            mesh
+        )
     );
-    rhoChemistryModel& chemistry = pChemistry();
+
+    rhoChemistryModel& chemistry = combustion->pChemistry();
 
     hsReactionThermo& thermo = chemistry.thermo();
 
@@ -57,20 +61,6 @@
 
     #include "compressibleCreatePhi.H"
 
-    DimensionedField<scalar, volMesh> kappa
-    (
-        IOobject
-        (
-            "kappa",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh,
-        dimensionedScalar("zero", dimless, 0.0)
-    );
-
     Info<< "Creating turbulence model\n" << endl;
     autoPtr<compressible::turbulenceModel> turbulence
     (
@@ -83,6 +73,9 @@
         )
     );
 
+     // Set the turbulence into the combustion model
+    combustion->setTurbulence(turbulence());
+
     Info<< "Creating multi-variate interpolation scheme\n" << endl;
     multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
 
@@ -92,16 +85,16 @@
     }
     fields.add(hs);
 
-    DimensionedField<scalar, volMesh> chemistrySh
+    volScalarField dQ
     (
         IOobject
         (
-            "chemistry::Sh",
+            "dQ",
             runTime.timeName(),
             mesh,
             IOobject::NO_READ,
-            IOobject::NO_WRITE
+            IOobject::AUTO_WRITE
         ),
         mesh,
-        dimensionedScalar("chemistrySh", dimEnergy/dimTime/dimVolume, 0.0)
+        dimensionedScalar("dQ", dimEnergy/dimTime, 0.0)
     );

applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/hsEqn.H

@@ -39,7 +39,7 @@
           + parcels.Sh(hs)
           + radiation->Shs(thermo)
           + energySource.Su()
-          + chemistrySh
+          + combustion->Sh()
         );
 
         thermo.correct();

applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/porousExplicitSourceReactingParcelFoam.C

@@ -40,11 +40,9 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "hReactionThermo.H"
 #include "turbulenceModel.H"
 #include "basicReactingMultiphaseCloud.H"
-#include "rhoChemistryModel.H"
+#include "rhoChemistryCombustionModel.H"
-#include "chemistrySolver.H"
 #include "radiationModel.H"
 #include "porousZones.H"
 #include "timeActivatedExplicitSource.H"
@@ -59,7 +57,6 @@ int main(int argc, char *argv[])
 
     #include "createTime.H"
     #include "createMesh.H"
-    #include "readChemistryProperties.H"
     #include "readGravitationalAcceleration.H"
     #include "createFields.H"
     #include "createRadiationModel.H"
@@ -90,7 +87,6 @@ int main(int argc, char *argv[])
 
         parcels.evolve();
 
-        #include "chemistry.H"
         #include "rhoEqn.H"
 
         // --- Pressure-velocity PIMPLE corrector loop
@@ -114,10 +110,7 @@ int main(int argc, char *argv[])
 
         rho = thermo.rho();
 
-        if (runTime.write())
+        runTime.write();
-        {
-            chemistry.dQ()().write();
-        }
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
             << "  ClockTime = " << runTime.elapsedClockTime() << " s"

applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/readChemistryProperties.H

@@ -1,23 +0,0 @@
-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/reactingParcelFilmFoam/Make/options

@@ -20,6 +20,7 @@ EXE_INC = \
     -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
     -I$(LIB_SRC)/ODE/lnInclude \
     -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/combustionModels/lnInclude \
     -I$(FOAM_SOLVERS)/combustion/reactingFoam
 
 
@@ -41,4 +42,5 @@ EXE_LIBS = \
     -lsurfaceFilmModels \
     -llagrangianIntermediate \
     -lODE \
-    -lsampling
+    -lsampling \
+    -lcombustionModels

applications/solvers/lagrangian/reactingParcelFilmFoam/YEqn.H

@@ -11,6 +11,8 @@ tmp<fv::convectionScheme<scalar> > mvConvection
 
 
 {
+    combustion->correct();
+    dQ = combustion->dQ();
     label inertIndex = -1;
     volScalarField Yt(0.0*Y[0]);
 
@@ -19,6 +21,7 @@ tmp<fv::convectionScheme<scalar> > mvConvection
         if (Y[i].name() != inertSpecie)
         {
             volScalarField& Yi = Y[i];
+
             solve
             (
                 fvm::ddt(rho, Yi)
@@ -27,7 +30,7 @@ tmp<fv::convectionScheme<scalar> > mvConvection
               ==
                 parcels.SYi(i, Yi)
               + surfaceFilm.Srho(i)
-              + kappa*chemistry.RR(i)().dimensionedInternalField(),
+              + combustion->R(Yi),
                 mesh.solver("Yi")
             );
 

applications/solvers/lagrangian/reactingParcelFilmFoam/createFields.H

@@ -1,10 +1,14 @@
-    Info<< "Reading thermophysical properties\n" << endl;
+    Info<< "Creating combustion model\n" << endl;
 
-    autoPtr<psiChemistryModel> pChemistry
+    autoPtr<combustionModels::psiChemistryCombustionModel> combustion
     (
-        psiChemistryModel::New(mesh)
+        combustionModels::psiChemistryCombustionModel::New
+        (
+            mesh
+        )
     );
-    psiChemistryModel& chemistry = pChemistry();
+
+    psiChemistryModel& chemistry = combustion->pChemistry();
 
     hsCombustionThermo& thermo = chemistry.thermo();
 
@@ -50,21 +54,6 @@
 
     #include "compressibleCreatePhi.H"
 
-    Info<< "Creating field kappa\n" << endl;
-    DimensionedField<scalar, volMesh> kappa
-    (
-        IOobject
-        (
-            "kappa",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh,
-        dimensionedScalar("zero", dimless, 0.0)
-    );
-
     Info<< "Creating turbulence model\n" << endl;
     autoPtr<compressible::turbulenceModel> turbulence
     (
@@ -77,6 +66,9 @@
         )
     );
 
+    // Set the turbulence into the combustion model
+    combustion->setTurbulence(turbulence());
+
     Info<< "Creating field DpDt\n" << endl;
     volScalarField DpDt
     (
@@ -131,16 +123,16 @@
         additionalControlsDict.lookup("solvePrimaryRegion")
     );
 
-    DimensionedField<scalar, volMesh> chemistrySh
+    volScalarField dQ
     (
         IOobject
         (
-            "chemistry::Sh",
+            "dQ",
             runTime.timeName(),
             mesh,
             IOobject::NO_READ,
-            IOobject::NO_WRITE
+            IOobject::AUTO_WRITE
         ),
         mesh,
-        dimensionedScalar("chemistrySh", dimEnergy/dimTime/dimVolume, 0.0)
+        dimensionedScalar("dQ", dimEnergy/dimTime, 0.0)
     );

applications/solvers/lagrangian/reactingParcelFilmFoam/hsEqn.H

@@ -9,7 +9,7 @@
       + parcels.Sh(hs)
       + surfaceFilm.Sh()
       + radiation->Shs(thermo)
-      + chemistrySh
+      + combustion->Sh()
     );
 
     hsEqn.relax();

applications/solvers/lagrangian/reactingParcelFilmFoam/reactingParcelFilmFoam.C

@@ -31,12 +31,10 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "hCombustionThermo.H"
 #include "turbulenceModel.H"
 #include "basicReactingCloud.H"
 #include "surfaceFilmModel.H"
-#include "psiChemistryModel.H"
+#include "psiChemistryCombustionModel.H"
-#include "chemistrySolver.H"
 #include "radiationModel.H"
 #include "SLGThermo.H"
 #include "pimpleControl.H"
@@ -49,7 +47,6 @@ int main(int argc, char *argv[])
 
     #include "createTime.H"
     #include "createMesh.H"
-    #include "readChemistryProperties.H"
     #include "readGravitationalAcceleration.H"
     #include "createFields.H"
     #include "createClouds.H"
@@ -83,7 +80,6 @@ int main(int argc, char *argv[])
 
         if (solvePrimaryRegion)
         {
-            #include "chemistry.H"
             #include "rhoEqn.H"
 
             // --- PIMPLE loop
@@ -107,10 +103,7 @@ int main(int argc, char *argv[])
 
             rho = thermo.rho();
 
-            if (runTime.write())
+            runTime.write();
-            {
-                chemistry.dQ()().write();
-            }
         }
         else
         {

applications/solvers/lagrangian/reactingParcelFilmFoam/readChemistryProperties.H

@@ -1,23 +0,0 @@
-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/reactingParcelFoam/Make/options

@@ -20,6 +20,7 @@ EXE_INC = \
     -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
     -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
     -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/combustionModels/lnInclude \
     -I$(FOAM_SOLVERS)/combustion/reactingFoam
 
 
@@ -45,4 +46,5 @@ EXE_LIBS = \
     -lODE \
     -lregionModels \
     -lsurfaceFilmModels \
-    -lsampling
+    -lsampling \
+    -lcombustionModels

applications/solvers/lagrangian/reactingParcelFoam/YEqn.H

@@ -11,6 +11,8 @@ tmp<fv::convectionScheme<scalar> > mvConvection
 
 
 {
+    combustion->correct();
+    dQ = combustion->dQ();
     label inertIndex = -1;
     volScalarField Yt(0.0*Y[0]);
 
@@ -19,17 +21,20 @@ tmp<fv::convectionScheme<scalar> > mvConvection
         if (Y[i].name() != inertSpecie)
         {
             volScalarField& Yi = Y[i];
-            solve
+
+            fvScalarMatrix YiEqn
             (
                 fvm::ddt(rho, Yi)
               + mvConvection->fvmDiv(phi, Yi)
               - fvm::laplacian(turbulence->muEff(), Yi)
               ==
                 parcels.SYi(i, Yi)
-              + kappa*chemistry.RR(i)().dimensionedInternalField(),
+              + combustion->R(Yi)
-                mesh.solver("Yi")
             );
 
+            YiEqn.relax();
+            YiEqn.solve(mesh.solver("Yi"));
+
             Yi.max(0.0);
             Yt += Yi;
         }

applications/solvers/lagrangian/reactingParcelFoam/createFields.H

@@ -1,10 +1,14 @@
-    Info<< "Reading thermophysical properties\n" << endl;
+    Info<< "Creating combustion model\n" << endl;
 
-    autoPtr<psiChemistryModel> pChemistry
+    autoPtr<combustionModels::psiChemistryCombustionModel> combustion
     (
-        psiChemistryModel::New(mesh)
+        combustionModels::psiChemistryCombustionModel::New
+        (
+            mesh
+        )
     );
-    psiChemistryModel& chemistry = pChemistry();
+
+    psiChemistryModel& chemistry = combustion->pChemistry();
 
     hsCombustionThermo& thermo = chemistry.thermo();
 
@@ -57,20 +61,6 @@
 
     #include "compressibleCreatePhi.H"
 
-    DimensionedField<scalar, volMesh> kappa
-    (
-        IOobject
-        (
-            "kappa",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh,
-        dimensionedScalar("zero", dimless, 0.0)
-    );
-
     Info<< "Creating turbulence model\n" << endl;
     autoPtr<compressible::turbulenceModel> turbulence
     (
@@ -83,6 +73,9 @@
         )
     );
 
+    // Set the turbulence into the combustion model
+    combustion->setTurbulence(turbulence());
+
     Info<< "Creating field DpDt\n" << endl;
     volScalarField DpDt
     (
@@ -98,16 +91,16 @@
     }
     fields.add(hs);
 
-    DimensionedField<scalar, volMesh> chemistrySh
+    volScalarField dQ
     (
         IOobject
         (
-            "chemistry::Sh",
+            "dQ",
             runTime.timeName(),
             mesh,
             IOobject::NO_READ,
-            IOobject::NO_WRITE
+            IOobject::AUTO_WRITE
         ),
         mesh,
-        dimensionedScalar("chemistry::Sh", dimEnergy/dimTime/dimVolume, 0.0)
+        dimensionedScalar("dQ", dimEnergy/dimTime, 0.0)
     );

applications/solvers/lagrangian/reactingParcelFoam/hsEqn.H

@@ -8,7 +8,7 @@
         DpDt
      +  parcels.Sh(hs)
      +  radiation->Shs(thermo)
-     +  chemistrySh
+     +  combustion->Sh()
     );
 
     hEqn.relax();

applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C

@@ -31,11 +31,9 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "hCombustionThermo.H"
 #include "turbulenceModel.H"
 #include "basicReactingCloud.H"
-#include "psiChemistryModel.H"
+#include "psiChemistryCombustionModel.H"
-#include "chemistrySolver.H"
 #include "radiationModel.H"
 #include "SLGThermo.H"
 #include "pimpleControl.H"
@@ -48,7 +46,6 @@ int main(int argc, char *argv[])
 
     #include "createTime.H"
     #include "createMesh.H"
-    #include "readChemistryProperties.H"
     #include "readGravitationalAcceleration.H"
     #include "createFields.H"
     #include "createClouds.H"
@@ -76,7 +73,6 @@ int main(int argc, char *argv[])
 
         parcels.evolve();
 
-        #include "chemistry.H"
         #include "rhoEqn.H"
 
         // --- Pressure-velocity PIMPLE corrector loop
@@ -100,10 +96,7 @@ int main(int argc, char *argv[])
 
         rho = thermo.rho();
 
-        if (runTime.write())
+        runTime.write();
-        {
-            chemistry.dQ()().write();
-        }
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
             << "  ClockTime = " << runTime.elapsedClockTime() << " s"

applications/solvers/lagrangian/reactingParcelFoam/readChemistryProperties.H

@@ -1,23 +0,0 @@
-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/sprayFoam/Make/options

@@ -21,6 +21,7 @@ EXE_INC = \
     -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
     -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
     -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/combustionModels/lnInclude \
     -I$(FOAM_SOLVERS)/lagrangian/reactingParcelFoam
 
 
@@ -47,4 +48,5 @@ EXE_LIBS = \
     -lODE \
     -lregionModels \
     -lsurfaceFilmModels \
-    -lsampling
+    -lsampling \
+    -lcombustionModels

applications/solvers/lagrangian/sprayFoam/chemistry.H

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

applications/solvers/lagrangian/sprayFoam/sprayFoam.C

@@ -31,11 +31,9 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "hCombustionThermo.H"
 #include "turbulenceModel.H"
 #include "basicSprayCloud.H"
-#include "psiChemistryModel.H"
+#include "psiChemistryCombustionModel.H"
-#include "chemistrySolver.H"
 #include "radiationModel.H"
 #include "SLGThermo.H"
 #include "pimpleControl.H"
@@ -48,7 +46,6 @@ int main(int argc, char *argv[])
 
     #include "createTime.H"
     #include "createMesh.H"
-    #include "readChemistryProperties.H"
     #include "readGravitationalAcceleration.H"
     #include "createFields.H"
     #include "createClouds.H"
@@ -76,7 +73,6 @@ int main(int argc, char *argv[])
 
         parcels.evolve();
 
-        #include "chemistry.H"
         #include "rhoEqn.H"
 
         // --- Pressure-velocity PIMPLE corrector loop

applications/solvers/multiphase/interMixingFoam/incompressibleThreePhaseMixture/threePhaseMixture.C

@@ -21,9 +21,6 @@ License
     You should have received a copy of the GNU General Public License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
-Class
-    threePhaseMixture
-
 \*---------------------------------------------------------------------------*/
 
 #include "threePhaseMixture.H"

applications/utilities/mesh/manipulation/polyDualMesh/meshDualiser.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -21,9 +21,6 @@ License
     You should have received a copy of the GNU General Public License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
-Class
-    meshDualiser
-
 \*---------------------------------------------------------------------------*/
 
 #include "meshDualiser.H"

applications/utilities/mesh/manipulation/setSet/writeFuns.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -22,18 +22,18 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::writeFunctions
+    Foam::writeFuns
 
 Description
     Various functions for collecting and writing binary data.
 
 SourceFiles
-    writeFunctions.C
+    writeFuns.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef writeFunctions_H
+#ifndef writeFuns_H
-#define writeFunctions_H
+#define writeFuns_H
 
 #include "labelList.H"
 #include "floatScalar.H"

applications/utilities/postProcessing/dataConversion/foamToVTK/writeFuns.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -22,18 +22,18 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::writeFunctions
+    Foam::writeFuns
 
 Description
     Various functions for collecting and writing binary data.
 
 SourceFiles
-    writeFunctions.C
+    writeFuns.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef writeFunctions_H
+#ifndef writeFuns_H
-#define writeFunctions_H
+#define writeFuns_H
 
 #include "floatScalar.H"
 #include "DynamicList.H"

bin/paraFoam

@@ -70,6 +70,8 @@ export LC_ALL=C
 # reader extension
 extension=OpenFOAM
 
+requirePV=1
+
 # parse options
 while [ "$#" -gt 0 ]
 do
@@ -97,10 +99,12 @@ do
         ;;
     -touch)
         optTouch=true
+        requirePV=0
         shift
         ;;
     -touchAll)
         optTouch=all
+        requirePV=0
         shift
         ;;
     --)
@@ -120,7 +124,7 @@ done
 #
 # check that reader module has been built
 #
-if [ ! -f $PV_PLUGIN_PATH/libPV3FoamReader_SM.so ]
+if [ $requirePV -eq 1 -a ! -f $PV_PLUGIN_PATH/libPV3FoamReader_SM.so ]
 then
     cat<< BUILDREADER
 

etc/config/settings.csh

@@ -209,10 +209,10 @@ case ThirdParty:
         breaksw
     case Gcc46:
     case Gcc46++0x:
-        set gcc_version=gcc-4.6.0
+        set gcc_version=gcc-4.6.1
-        set gmp_version=gmp-5.0.1
+        set gmp_version=gmp-5.0.2
-        set mpfr_version=mpfr-2.4.2
+        set mpfr_version=mpfr-3.0.1
-        set mpc_version=mpc-0.8.1
+        set mpc_version=mpc-0.9
         breaksw
     case Gcc45:
     case Gcc45++0x:

etc/config/settings.sh

@@ -228,10 +228,11 @@ OpenFOAM | ThirdParty)
         mpfr_version=mpfr-2.4.2
         ;;
     Gcc46 | Gcc46++0x)
-        gcc_version=gcc-4.6.0
+        gcc_version=gcc-4.6.1
-        gmp_version=gmp-5.0.1
+        gmp_version=gmp-5.0.2
-        mpfr_version=mpfr-2.4.2
+        mpfr_version=mpfr-3.0.1
-        mpc_version=mpc-0.8.1
+        mpc_version=mpc-0.9
+        gmpPACKAGE=gmp-5.0.2
         ;;
     Gcc45 | Gcc45++0x)
         gcc_version=gcc-4.5.2

src/OSspecific/POSIX/fileMonitor.C

@@ -21,9 +21,6 @@ License
     You should have received a copy of the GNU General Public License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
-Class
-    fileMonitor
-
 \*----------------------------------------------------------------------------*/
 
 #include "fileMonitor.H"

src/OpenFOAM/db/error/StaticAssert.H

@@ -22,7 +22,7 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::StaticAssertFailed
+    Foam::StaticAssertionFailed
 
 Description
     Macros and classes to provide static (compile-time) assertions.

src/OpenFOAM/fields/GeometricFields/geometricOneField/geometricOneField.H

@@ -22,7 +22,7 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::GeometricField
+    Foam::geometricOneField
 
 Description
     A class representing the concept of a GeometricField of 1 used to avoid

src/OpenFOAM/fields/GeometricFields/geometricZeroField/geometricZeroField.H

@@ -22,7 +22,7 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::GeometricField
+    Foam::geometricZeroField
 
 Description
     A class representing the concept of a GeometricField of 1 used to avoid

src/OpenFOAM/interpolations/interpolationTable/tableReaders/tableReaders.H

@@ -21,9 +21,6 @@ License
     You should have received a copy of the GNU General Public License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
-Class
-    Foam::tableReader
-
 SourceFiles
     tableReaders.C
 

src/OpenFOAM/meshes/polyMesh/zones/zone/zone.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -22,7 +22,7 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::cellZone
+    Foam::zone
 
 Description
     Base class for zones

src/combustionModels/FSD/FSD.C

@@ -0,0 +1,362 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+
+    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 "FSD.H"
+#include "addToRunTimeSelectionTable.H"
+#include "LESModel.H"
+
+namespace Foam
+{
+namespace combustionModels
+{
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class CombThermoType, class ThermoType>
+FSD<CombThermoType, ThermoType>::FSD
+(
+    const word& modelType, const fvMesh& mesh
+)
+:
+    singleStepCombustion<CombThermoType, ThermoType>(modelType, mesh),
+    reactionRateFlameArea_
+    (
+        reactionRateFlameArea::New
+        (
+            this->coeffs(),
+            this->mesh(),
+            *this
+        )
+    ),
+    ft_
+    (
+        IOobject
+        (
+            "ft",
+            this->mesh().time().timeName(),
+            this->mesh(),
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        this->mesh(),
+        dimensionedScalar("zero", dimless, 0.0)
+    ),
+    YFuelFuelStream_
+    (
+        dimensionedScalar("YFuelStream", dimless, 1.0)
+    ),
+    YO2OxiStream_
+    (
+        dimensionedScalar("YOxiStream", dimless, 0.23)
+    ),
+    Cv_(readScalar(this->coeffs().lookup("Cv"))),
+    C_(5.0),
+    ftMin_(0.0),
+    ftMax_(1.0),
+    ftDim_(300),
+    ftVarMin_(readScalar(this->coeffs().lookup("ftVarMin")))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * * //
+
+template<class CombThermoType, class ThermoType>
+FSD<CombThermoType, ThermoType>::~FSD()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+template<class CombThermoType, class ThermoType>
+void FSD<CombThermoType, ThermoType>::calculateSourceNorm()
+{
+    this->singleMixture_.fresCorrect();
+
+    const label fuelI = this->singleMixture_.fuelIndex();
+
+    const volScalarField& YFuel = this->thermo_->composition().Y()[fuelI];
+
+    const volScalarField& YO2 = this->thermo_->composition().Y("O2");
+
+    const dimensionedScalar s = this->singleMixture_.s();
+
+    ft_ =
+        (s*YFuel - (YO2 - YO2OxiStream_))/(s*YFuelFuelStream_ + YO2OxiStream_);
+
+
+    volVectorField nft = fvc::grad(ft_);
+
+    volScalarField mgft = mag(nft);
+
+    surfaceVectorField SfHat = this->mesh().Sf()/this->mesh().magSf();
+
+    volScalarField cAux = scalar(1) - ft_;
+
+    dimensionedScalar dMgft = 1.0e-3*
+        (ft_*cAux*mgft)().weightedAverage(this->mesh().V())
+       /((ft_*cAux)().weightedAverage(this->mesh().V()) + SMALL)
+      + dimensionedScalar("ddMgft", mgft.dimensions(), SMALL);
+
+    mgft += dMgft;
+
+    nft /= mgft;
+
+    const volVectorField& U = YO2.db().lookupObject<volVectorField>("U");
+
+    const volScalarField sigma =
+        (nft & nft)*fvc::div(U) - (nft & fvc::grad(U) & nft);
+
+    reactionRateFlameArea_->correct(sigma);
+
+    const volScalarField& omegaFuel = reactionRateFlameArea_->omega();
+
+
+    const scalar ftStoich =
+        YO2OxiStream_.value()
+       /(
+            s.value()*YFuelFuelStream_.value() + YO2OxiStream_.value()
+        );
+
+    tmp<volScalarField> tPc
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "Pc",
+                U.time().timeName(),
+                U.db(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            U.mesh(),
+            dimensionedScalar("Pc", dimless, 0)
+        )
+    );
+
+    volScalarField& pc = tPc();
+
+    tmp<volScalarField> tomegaFuel
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "omegaFuelBar",
+                U.time().timeName(),
+                U.db(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            U.mesh(),
+            dimensionedScalar
+            (
+                "omegaFuelBar",
+                omegaFuel.dimensions(),
+                0
+            )
+        )
+    );
+
+    volScalarField& omegaFuelBar = tomegaFuel();
+
+    // Calculation of the mixture fraction variance (ftVar)
+    const compressible::LESModel& lesModel =
+        YO2.db().lookupObject<compressible::LESModel>("LESProperties");
+
+    const volScalarField& delta = lesModel.delta();
+    const volScalarField ftVar = Cv_*sqr(delta)*sqr(mgft);
+
+    // Thickened flame (average flame thickness for counterflow configuration
+    // is 1.5 mm)
+
+    volScalarField  deltaF =
+        lesModel.delta()/dimensionedScalar("flame",dimLength, 1.5e-3);
+
+    // Linear correlation between delta and flame thickness
+    volScalarField omegaF = max(deltaF*(4.0/3.0) + (2.0/3.0), 1.0);
+
+    scalar deltaFt = 1.0/ftDim_;
+
+    forAll(ft_, cellI)
+    {
+        if(ft_[cellI] > ftMin_ && ft_[cellI] < ftMax_)
+        {
+            scalar ftCell = ft_[cellI];
+
+            if(ftVar[cellI] > ftVarMin_) //sub-grid beta pdf of ft_
+            {
+                scalar ftVarc = ftVar[cellI];
+                scalar a =
+                    max(ftCell*(ftCell*(1.0 - ftCell)/ftVarc - 1.0), 0.0);
+                scalar b = max(a/ftCell - a, 0.0);
+
+                for(int i=1; i<ftDim_; i++)
+                {
+                    scalar ft = i*deltaFt;
+                    pc[cellI] += pow(ft, a-1.0)*pow(1.0 - ft, b - 1.0)*deltaFt;
+                }
+
+                for(int i=1; i<ftDim_; i++)
+                {
+                    scalar ft = i*deltaFt;
+                    omegaFuelBar[cellI] +=
+                        omegaFuel[cellI]/omegaF[cellI]
+                       *exp
+                       (
+                        -sqr(ft - ftStoich)
+                        /(2.0*sqr(0.01*omegaF[cellI]))
+                       )
+                       *pow(ft, a - 1.0)
+                       *pow(1.0 - ft, b - 1.0)
+                       *deltaFt;
+                }
+                omegaFuelBar[cellI] /= max(pc[cellI], 1e-4);
+            }
+            else
+            {
+                omegaFuelBar[cellI] =
+                    (omegaFuel[cellI]/omegaF[cellI])
+                   *exp
+                   (
+                    -sqr(ftCell - ftStoich)/(2.0*sqr(0.01*omegaF[cellI]))
+                   );
+            }
+
+        }
+        else
+        {
+                omegaFuelBar[cellI] = 0.0;
+        }
+    }
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+//  Combustion progress variable (c).
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    List<label> productsIndex(2, -1);
+    {
+        label i = 0;
+        forAll (this->singleMixture_.specieProd(), specieI)
+        {
+            if (this->singleMixture_.specieProd()[specieI] < 0)
+            {
+                productsIndex[i] = specieI;
+                i++;
+            }
+        }
+    }
+
+
+    // Flamelet probability of the progress c based on IFC (reuse pc)
+    scalar YprodTotal = 0;
+    forAll (productsIndex, j)
+    {
+        YprodTotal += this->singleMixture_.Yprod0()[productsIndex[j]];
+    }
+
+    forAll(ft_, cellI)
+    {
+        if(ft_[cellI] < ftStoich)
+        {
+            pc[cellI] = ft_[cellI]*(YprodTotal/ftStoich);
+        }
+        else
+        {
+            pc[cellI] = (1.0 - ft_[cellI])*(YprodTotal/(1.0 - ftStoich));
+        }
+    }
+
+    tmp<volScalarField> tproducts
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "products",
+                U.time().timeName(),
+                U.db(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            U.mesh(),
+            dimensionedScalar("products", dimless, 0)
+        )
+    );
+
+    volScalarField& products = tproducts();
+
+    forAll (productsIndex, j)
+    {
+        label specieI = productsIndex[j];
+        const volScalarField& Yp = this->thermo_->composition().Y()[specieI];
+        products += Yp;
+    }
+
+    volScalarField c = max(scalar(1.0) - products/max(pc, 1e-5), 0.0);
+
+    pc = min(C_*c, scalar(1.0));
+
+    const volScalarField fres = this->singleMixture_.fres(fuelI);
+
+    this->wFuel_ == mgft*pc*omegaFuelBar;
+}
+
+
+template<class CombThermoType, class ThermoType>
+void FSD<CombThermoType, ThermoType>::correct()
+{
+    this->wFuel_ ==
+        dimensionedScalar("zero", dimMass/pow3(dimLength)/dimTime, 0.0);
+
+    if (this->active())
+    {
+        calculateSourceNorm();
+    }
+}
+
+
+template<class CombThermoType, class ThermoType>
+bool FSD<CombThermoType, ThermoType>::read()
+{
+    if (singleStepCombustion<CombThermoType, ThermoType>::read())
+    {
+        this->coeffs().lookup("Cv") >> Cv_ ;
+        this->coeffs().lookup("ftVarMin") >> ftVarMin_;
+        reactionRateFlameArea_->read(this->coeffs());
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace combustionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/combustionModels/FSD/FSD.H

@@ -0,0 +1,173 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::combustionModels::FSD
+
+Description
+
+    Flame Surface Dennsity (FDS) combustion model.
+
+    The fuel source term is given by mgft*pc*omegaFuelBar.
+
+    where:
+          mgft: filtered flame area.
+          pc:   probability of the combustion progress.
+          omegaFuelBar: filtered consumption speed per unit of flame area.
+
+    pc is considered from the IFC solution.
+    omegaFuelBar is calculated solving a relaxation equation which tends to
+    omegaEq. This omegaEq is obtained from the flamelet solution for
+    different strain rates and fit using a expential distribution.
+
+    The spacial distribution of the consumption speed (omega) is obtained also
+    from a strained flamelet solution and it is assumed to have a guassian
+    distribution.
+
+    If the grid resolution is not enough to resolve the flame, the consumption
+    speed distribution is linearly thickened conserving the overall heat
+    release.
+
+    If the turbulent fluctuation of the mixture fraction at the sub-grid level
+    is large (>1E-04) then a beta pdf is used for filtering.
+
+    At the moment the flame area combustion model is only fit to work in a LES
+    frame work. In RAS the subgrid fluctiuation has to be solved by an extra
+    transport equation.
+
+SourceFiles
+    FSD.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef FSD_H
+#define FSD_H
+
+#include "singleStepCombustion.H"
+#include "reactionRateFlameArea.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace combustionModels
+{
+
+/*---------------------------------------------------------------------------*\
+                          Class FSD Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CombThermoType, class ThermoType>
+class FSD
+:
+    public singleStepCombustion <CombThermoType, ThermoType>
+{
+    // Private data
+
+        //- Auto pointer to consumption speed per unit of flame area model
+        autoPtr<reactionRateFlameArea> reactionRateFlameArea_;
+
+        //- Mixture fraction
+        volScalarField ft_;
+
+        //- Fuel mass concentration on the fuel stream
+        dimensionedScalar YFuelFuelStream_;
+
+        //- Oxygen mass concentration on the oxydizer stream
+        dimensionedScalar YO2OxiStream_;
+
+        //- Similarity constant for the sub-grid ft fluctuations
+        scalar Cv_;
+
+        //- Model constant
+        scalar C_;
+
+        //- Lower flammability limit
+        scalar ftMin_;
+
+        //- Upper flammability limit
+        scalar ftMax_;
+
+        //- Dimension of the ft space. Used to integrate the beta-pdf
+        scalar ftDim_;
+
+        //- Minimum mixture freaction variance to calculate pdf
+        scalar ftVarMin_;
+
+
+    // Private Member Functions
+
+        //- Calculate the normalised fuel source term
+        void calculateSourceNorm();
+
+        //- Disallow copy construct
+        FSD(const FSD&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const FSD&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("FSD");
+
+
+    // Constructors
+
+        //- Construct from components
+        FSD(const word& modelType, const fvMesh& mesh);
+
+
+    // Destructor
+        virtual ~FSD();
+
+
+    // Evolution
+
+        //- Correct combustion rate
+        virtual void correct();
+
+
+    // I-O
+
+        //- Update properties
+        virtual bool read();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace combustionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "FSD.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/combustionModels/FSD/FSDs.C

@@ -0,0 +1,49 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "makeCombustionTypes.H"
+
+#include "thermoPhysicsTypes.H"
+#include "psiCombustionModel.H"
+#include "FSD.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace combustionModels
+{
+    makeCombustionTypesThermo
+    (
+        FSD,
+        psiCombustionModel,
+        gasThermoPhysics
+    );
+
+}
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/combustionModels/FSD/reactionRateFlameAreaModels/consumptionSpeed/consumptionSpeed.C

@@ -0,0 +1,137 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 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 "consumptionSpeed.H"
+
+/* * * * * * * * * * * * * * * private static data * * * * * * * * * * * * * */
+
+namespace Foam
+{
+    defineTypeNameAndDebug(consumptionSpeed, 0);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::consumptionSpeed::consumptionSpeed
+(
+    const dictionary& dict
+)
+:   omega0_(readScalar(dict.lookup("omega0"))),
+    eta_(readScalar(dict.lookup("eta"))),
+    sigmaExt_(readScalar(dict.lookup("sigmaExt"))),
+    omegaMin_(readScalar(dict.lookup("omegaMin")))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::consumptionSpeed::~consumptionSpeed()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::scalar Foam::consumptionSpeed::omega0Sigma
+(
+    scalar sigma,
+    scalar a
+) const
+{
+        scalar omega0 = 0.0;
+        if (sigma < sigmaExt_)
+        {
+            omega0 =
+            max
+            (
+                a*omega0_*(1.0 - exp(eta_*max(sigma, 0.0))),
+                omegaMin_
+            ) ;
+        }
+        return omega0;
+}
+
+
+Foam::tmp<Foam::volScalarField> Foam::consumptionSpeed::omega0Sigma
+(
+    const volScalarField& sigma
+)
+{
+    tmp<volScalarField> tomega0
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "omega0",
+                sigma.time().timeName(),
+                sigma.db(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            sigma.mesh(),
+            dimensionedScalar
+            (
+                "omega0",
+                dimensionSet(1, -2, -1, 0, 0, 0, 0),
+                0
+            )
+        )
+    );
+
+    volScalarField& omega0 = tomega0();
+
+    forAll(omega0, celli)
+    {
+        omega0[celli] = omega0Sigma(sigma[celli], 1.0);
+    }
+
+    forAll(omega0.boundaryField(), patchi)
+    {
+        forAll(omega0.boundaryField()[patchi], facei)
+        {
+            omega0.boundaryField()[patchi][facei] =
+                omega0Sigma
+                (
+                    sigma.boundaryField()[patchi][facei],
+                    1.0
+                );
+        }
+    }
+
+    return tomega0;
+}
+
+
+void  Foam::consumptionSpeed::read(const dictionary& dict)
+{
+    dict.lookup("omega0") >> omega0_ ;
+    dict.lookup("eta") >> eta_ ;
+    dict.lookup("sigmaExt") >> sigmaExt_;
+    dict.lookup("omegaMin") >> omegaMin_;
+}
+
+// ************************************************************************* //

src/combustionModels/FSD/reactionRateFlameAreaModels/consumptionSpeed/consumptionSpeed.H

@@ -0,0 +1,139 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::reactionRateFlameAreaModels::consumptionSpeed
+
+Description
+    Correlation function for laminar consumption speed obtained from flamelet
+    solution at increasing strain rates.
+
+SourceFiles
+    consumptionSpeed.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef consumptionSpeed_H
+#define consumptionSpeed_H
+
+#include "IOdictionary.H"
+#include "volFields.H"
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class consumptionSpeed Declaration
+\*---------------------------------------------------------------------------*/
+
+class consumptionSpeed
+{
+    // Private Data
+
+
+        //- Maximum consumption speed
+        scalar omega0_;
+
+        //- Exponential factor
+        scalar eta_;
+
+        //- Extinction strain
+        scalar sigmaExt_;
+
+        //- Minimum consumption speed
+        scalar omegaMin_;
+
+
+    // Private member functions
+
+        //- Return consumption rate
+        scalar omega0Sigma(scalar sigma, scalar a) const;
+
+        //- Disallow copy construct
+        consumptionSpeed(const consumptionSpeed&);
+
+         //- Disallow default bitwise assignment
+        void operator=(const consumptionSpeed&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("consumptionSpeed");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        consumptionSpeed(const dictionary& dict);
+
+
+    //- Destructor
+        virtual ~consumptionSpeed();
+
+
+    // Member functions
+
+        //- Return speed consumption rate temp
+        tmp<volScalarField> omega0Sigma(const volScalarField& sigma);
+
+
+        // Access functions
+
+            scalar omega0() const
+            {
+                return omega0_;
+            }
+
+            scalar eta() const
+            {
+                return eta_;
+            }
+
+            scalar sigmaExt() const
+            {
+                return sigmaExt_;
+            }
+
+            scalar omegaMin() const
+            {
+                return omegaMin_;
+            }
+
+
+    // I-O
+
+        //- Update properties
+        void read(const dictionary& dict);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/combustionModels/FSD/reactionRateFlameAreaModels/reactionRateFlameArea/reactionRateFlameArea.C

@@ -0,0 +1,107 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+
+    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 "reactionRateFlameArea.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(reactionRateFlameArea, 0);
+    defineRunTimeSelectionTable(reactionRateFlameArea, dictionary);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::reactionRateFlameArea::reactionRateFlameArea
+(
+    const dictionary& dict,
+    const fvMesh& mesh,
+    const combustionModel& combModel
+)
+:
+    coeffDict_(dictionary::null),
+    mesh_(mesh),
+    combModel_(combModel),
+    fuel_(dict.lookup("fuel")),
+    omega_
+    (
+        IOobject
+        (
+            "omega",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh_
+    )
+{
+}
+
+
+Foam::reactionRateFlameArea::reactionRateFlameArea
+(
+    const word& modelType,
+    const dictionary& dict,
+    const fvMesh& mesh,
+    const combustionModel& combModel
+)
+:
+    coeffDict_(dict.subDict(modelType + "Coeffs")),
+    mesh_(mesh),
+    combModel_(combModel),
+    fuel_(dict.lookup("fuel")),
+    omega_
+    (
+        IOobject
+        (
+            "omega",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh_
+    )
+{
+}
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::reactionRateFlameArea::~reactionRateFlameArea()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+bool Foam::reactionRateFlameArea::read(const dictionary& dict)
+{
+    dict.lookup("fuel") >> fuel_;
+    return true;
+
+}
+
+// ************************************************************************* //

src/combustionModels/FSD/reactionRateFlameAreaModels/reactionRateFlameArea/reactionRateFlameArea.H

@@ -0,0 +1,176 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+
+Class
+    Foam::reactionRateFlameArea
+
+Description
+    Abstract class for reaction rate per flame area unit
+
+SourceFiles
+    reactionRateFlameArea.C
+    reactionRateFlameAreaNew.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef reactionRateFlameArea_H
+#define reactionRateFlameArea_H
+
+#include "runTimeSelectionTables.H"
+#include "dictionary.H"
+#include "autoPtr.H"
+#include "volFields.H"
+#include "combustionModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+class fvMesh;
+
+/*---------------------------------------------------------------------------*\
+                           Class reactionRateFlameArea Declaration
+\*---------------------------------------------------------------------------*/
+
+class reactionRateFlameArea
+{
+
+protected:
+
+    // Protected data
+
+        //- Dictionary
+        dictionary coeffDict_;
+
+        //- Mesh reference
+        const fvMesh& mesh_;
+
+        //- Combstion model owner
+        const combustionModel& combModel_;
+
+        //- Fuel name
+        word fuel_;
+
+        //- Fuel consumption rate per unit of flame area
+        volScalarField omega_;
+
+
+private:
+
+    // Private member functions
+
+         //- Disallow copy construct
+        reactionRateFlameArea(const reactionRateFlameArea&);
+
+         //- Disallow default bitwise assignment
+        void operator=(const reactionRateFlameArea&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("reactionRateFlameArea");
+
+
+    // Declare run-time constructor selection table
+
+        declareRunTimeSelectionTable
+        (
+            autoPtr,
+            reactionRateFlameArea,
+            dictionary,
+            (
+                const word modelType,
+                const dictionary& dict,
+                const fvMesh& mesh,
+                const combustionModel& combModel
+            ),
+            (modelType, dict, mesh, combModel)
+        );
+
+
+    // Constructors
+
+        //- Construct from dictionary and hsCombustionThermo
+        reactionRateFlameArea
+        (
+            const dictionary& dict,
+            const fvMesh& mesh,
+            const combustionModel& combModel
+        );
+
+        //- Construct from components
+        reactionRateFlameArea
+        (
+            const word& modelType,
+            const dictionary& dict,
+            const fvMesh& mesh,
+            const combustionModel& combModel
+        );
+
+
+    // Selector
+
+        static autoPtr<reactionRateFlameArea> New
+        (
+            const dictionary& dict,
+            const fvMesh& mesh,
+            const combustionModel& combModel
+        );
+
+
+    // Destructor
+
+        virtual ~reactionRateFlameArea();
+
+
+    // Member functions
+
+        //- Access functions
+
+            //- Return omega
+            const volScalarField& omega() const
+            {
+                return omega_;
+            }
+
+
+        //- Correct omega
+        virtual void correct(const volScalarField& sigma) = 0;
+
+        //- Update from dictionary
+        virtual bool read(const dictionary& dictProperties);
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/combustionModels/FSD/reactionRateFlameAreaModels/reactionRateFlameArea/reactionRateFlameAreaNew.C

@@ -0,0 +1,69 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+
+    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 "reactionRateFlameArea.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::reactionRateFlameArea> Foam::reactionRateFlameArea::New
+(
+    const dictionary& dict,
+    const fvMesh& mesh,
+    const combustionModel& combModel
+)
+{
+    word reactionRateFlameAreaType
+    (
+        dict.lookup("reactionRateFlameArea")
+    );
+
+    Info<< "Selecting reaction rate flame area correlation "
+        << reactionRateFlameAreaType << endl;
+
+    dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(reactionRateFlameAreaType);
+
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    {
+        FatalIOErrorIn
+        (
+            "reactionRateFlameArea::New(const hsCombustionThermo&)",
+            dict
+        )   << "Unknown reactionRateFlameArea type "
+            << reactionRateFlameAreaType << endl << endl
+            << "Valid reaction rate flame area types are :" << endl
+            << dictionaryConstructorTablePtr_->toc()
+            << exit(FatalIOError);
+    }
+
+    const label tempOpen = reactionRateFlameAreaType.find('<');
+
+    const word className = reactionRateFlameAreaType(0, tempOpen);
+
+    return autoPtr<reactionRateFlameArea>
+        (cstrIter()(className, dict, mesh, combModel));
+}
+
+
+// ************************************************************************* //

src/combustionModels/FSD/reactionRateFlameAreaModels/relaxation/relaxation.C

@@ -0,0 +1,155 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+
+    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 "relaxation.H"
+#include "addToRunTimeSelectionTable.H"
+#include "fvm.H"
+#include "LESModel.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace reactionRateFlameAreaModels
+{
+    defineTypeNameAndDebug(relaxation, 0);
+    addToRunTimeSelectionTable
+    (
+        reactionRateFlameArea,
+        relaxation,
+        dictionary
+    );
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::reactionRateFlameAreaModels::relaxation::relaxation
+(
+    const word modelType,
+    const dictionary& dict,
+    const fvMesh& mesh,
+    const combustionModel& combModel
+)
+:
+    reactionRateFlameArea(modelType, dict, mesh, combModel),
+    correlation_(dict.subDict(typeName + "Coeffs").subDict(fuel_)),
+    C_(readScalar(dict.subDict(typeName + "Coeffs").lookup("C"))),
+    alpha_(readScalar(dict.subDict(typeName + "Coeffs").lookup("alpha")))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::reactionRateFlameAreaModels::relaxation::~relaxation()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+void Foam::reactionRateFlameAreaModels::relaxation::correct
+(
+    const volScalarField& sigma
+)
+{
+
+    dimensionedScalar omega0
+    (
+        "omega0",
+        dimensionSet(1, -2, -1, 0, 0, 0, 0),
+        correlation_.omega0()
+    );
+
+    dimensionedScalar sigmaExt
+    (
+        "sigmaExt",
+        dimensionSet(0, 0, -1, 0, 0, 0, 0),
+        correlation_.sigmaExt()
+    );
+
+    dimensionedScalar omegaMin
+    (
+        "omegaMin",
+        omega0.dimensions(),
+        1e-4
+    );
+
+    const compressible::LESModel& lesModel =
+        omega_.db().lookupObject<compressible::LESModel>("LESProperties");
+
+    // Total strain : resolved and sub-grid (just LES for now)
+    const volScalarField sigmaTotal =
+        sigma + alpha_*lesModel.epsilon()/(lesModel.k() + lesModel.kMin());
+
+    const volScalarField omegaInf = correlation_.omega0Sigma(sigmaTotal);
+
+    dimensionedScalar sigma0("sigma0", sigma.dimensions(), 0.0);
+
+    const volScalarField tau = C_*mag(sigmaTotal);
+
+    volScalarField Rc =
+        (tau*omegaInf*(omega0 - omegaInf) + sqr(omegaMin)*sigmaExt)
+       /(sqr(omega0 - omegaInf) + sqr(omegaMin));
+
+    const volScalarField rho(combModel_.rho());
+    const surfaceScalarField phi(combModel_.phi());
+
+    solve
+    (
+         fvm::ddt(rho, omega_)
+       + fvm::div(phi, omega_, "div(phi,omega)")
+      ==
+         rho*Rc*omega0
+       - fvm::SuSp(rho*(tau + Rc), omega_)
+    );
+
+    omega_.min(omega0);
+    omega_.max(0.0);
+}
+
+
+bool  Foam::reactionRateFlameAreaModels::relaxation::read
+(
+    const dictionary& dict
+)
+{
+    if (reactionRateFlameArea::read(dict))
+    {
+        coeffDict_ = dict.subDict(typeName + "Coeffs");
+        coeffDict_.lookup("C") >> C_;
+        coeffDict_.lookup("alpha") >> alpha_;
+        correlation_.read
+        (
+            coeffDict_.subDict(fuel_)
+        );
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+// ************************************************************************* //

src/combustionModels/FSD/reactionRateFlameAreaModels/relaxation/relaxation.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2009-2011 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2011-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -21,72 +21,97 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::combustionModel::noCombustion
+    Foam::reactionRateFlameAreaModels::relaxation
 
 Description
-    Dummy combustion model for 'none' option
+    Consumption rate per unit of flame area obtained from a relaxation equation
 
 SourceFiles
-    noCombustion.C
+    relaxation.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef noCombustion_H
+#ifndef relaxation_H
-#define noCombustion_H
+#define relaxation_H
 
-#include "combustionModel.H"
+#include "reactionRateFlameArea.H"
+#include "consumptionSpeed.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
-namespace combustionModels
+namespace reactionRateFlameAreaModels
 {
 
 /*---------------------------------------------------------------------------*\
-                      Class noCombustion Declaration
+                           Class relaxation Declaration
 \*---------------------------------------------------------------------------*/
 
-class noCombustion
+class relaxation
 :
-    public combustionModel
+    public reactionRateFlameArea
 {
+    // Private Data
+
+        //- Correlation
+        consumptionSpeed correlation_;
+
+        //- Proportionality constant for time scale in the relaxation Eq.
+        scalar C_;
+
+        //- Proportionality constant for sub-grid strain
+        scalar alpha_;
+
+
      // Private Member Functions
 
         //- Disallow copy construct
-        noCombustion(const noCombustion&);
+        relaxation(const relaxation&);
 
         //- Disallow default bitwise assignment
-        void operator=(const noCombustion&);
+        void operator=(const relaxation&);
 
 
 public:
 
     //- Runtime type information
-    TypeName("none");
+    TypeName("relaxation");
 
 
     // Constructors
 
-        //- Construct from components
+        //- Construct from dictionary and hsCombustionThermo
-        noCombustion
+        relaxation
         (
-            const dictionary& combustionProperties,
+            const word modelType,
-            hsCombustionThermo& thermo,
+            const dictionary& dictCoeffs,
-            const compressible::turbulenceModel& turbulence,
+            const fvMesh& mesh,
-            const surfaceScalarField& phi,
+            const combustionModel& combModel
-            const volScalarField& rho
         );
 
 
-    //- Destructor
+    // Destructor
-    virtual ~noCombustion();
+
+        virtual ~relaxation();
+
+
+    // Member functions
+
+        //- Correct omega
+        virtual void correct(const volScalarField& sigma);
+
+
+    // I-O
+
+        //- Update properties from given dictionary
+        virtual bool read(const dictionary& dictProperties);
 };
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-} // End namespace combustionModels
+} // End reactionRateFlameAreaModels
 } // End namespace Foam
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/combustionModels/Make/files

@@ -1,9 +1,27 @@
 combustionModel/combustionModel.C
-combustionModel/combustionModelNew.C
 
-infinitelyFastChemistry/infinitelyFastChemistry.C
+psiCombustionModel/psiCombustionModel.C
+psiCombustionModel/psiCombustionModelNew.C
 
-noCombustion/noCombustion.C
+rhoCombustionModel/rhoCombustionModel.C
+rhoCombustionModel/rhoCombustionModelNew.C
+
+infinitelyFastChemistry/infinitelyFastChemistrys.C
+
+psiChemistryCombustionModel/psiChemistryCombustionModel.C
+psiChemistryCombustionModel/psiChemistryCombustionModelNew.C
+
+rhoChemistryCombustionModel/rhoChemistryCombustionModel.C
+rhoChemistryCombustionModel/rhoChemistryCombustionModelNew.C
+
+PaSR/PaSRs.C
+
+FSD/reactionRateFlameAreaModels/consumptionSpeed/consumptionSpeed.C
+FSD/reactionRateFlameAreaModels/reactionRateFlameArea/reactionRateFlameArea.C
+FSD/reactionRateFlameAreaModels/reactionRateFlameArea/reactionRateFlameAreaNew.C
+FSD/reactionRateFlameAreaModels/relaxation/relaxation.C
+
+FSD/FSDs.C
 
 LIB = $(FOAM_LIBBIN)/libcombustionModels
 

src/combustionModels/Make/options

@@ -2,8 +2,14 @@ EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
     -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
+    -I$(LIB_SRC)/turbulenceModels/ \
+    -I$(LIB_SRC)/turbulenceModels/LES/LESdeltas/lnInclude \
+    -I$(LIB_SRC)/turbulenceModels/LES/LESfilters/lnInclude \
+    -I$(LIB_SRC)/turbulenceModels/compressible/LES/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude
 
 LIB_LIBS = \
-    -lfiniteVolume
+    -lfiniteVolume \
+    -lchemistryModel

src/combustionModels/PaSR/PaSR.C

@@ -0,0 +1,253 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+
+    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 "PaSR.H"
+#include "fvmSup.H"
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class CombThermoType>
+Foam::combustionModels::PaSR<CombThermoType>::PaSR
+(
+    const word& modelType,
+    const fvMesh& mesh
+)
+:
+    CombThermoType(modelType, mesh),
+    Cmix_(this->coeffs().lookup("Cmix")),
+    turbulentReaction_(this->coeffs().lookup("turbulentReaction")),
+    kappa_
+    (
+        IOobject
+        (
+            "kappa",
+            mesh.time().timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE,
+            false
+        ),
+        mesh,
+        dimensionedScalar("kappa", dimless, 0.0)
+    ),
+    useReactionRate_(this->coeffs().lookupOrDefault("useReactionRate", false))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * * //
+
+template<class CombThermoType>
+Foam::combustionModels::PaSR<CombThermoType>::~PaSR()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+
+template<class CombThermoType>
+Foam::tmp<Foam::volScalarField>
+Foam::combustionModels::PaSR<CombThermoType>::tc() const
+{
+    return this->pChemistry_->tc();
+}
+
+
+template<class CombThermoType>
+void Foam::combustionModels::PaSR<CombThermoType>::correct()
+{
+    if (this->active())
+    {
+        if (!useReactionRate_)
+        {
+            this->pChemistry_->solve
+            (
+                this->mesh().time().value()-this->mesh().time().deltaTValue(),
+                this->mesh().time().deltaTValue()
+            );
+        }
+        else
+        {
+            this->pChemistry_->calculate();
+        }
+
+        if (turbulentReaction_)
+        {
+            tmp<volScalarField> tepsilon(this->turbulence().epsilon());
+            const volScalarField& epsilon = tepsilon();
+            tmp<volScalarField> tmuEff(this->turbulence().muEff());
+            const volScalarField& muEff = tmuEff();
+            tmp<volScalarField> ttc(tc());
+            const volScalarField& tc = ttc();
+            forAll(epsilon, i)
+            {
+                if (epsilon[i] > 0)
+                {
+                    const dimensionedScalar e0
+                    (
+                        "e0",
+                        sqr(dimLength)/pow3(dimTime), SMALL
+                    );
+
+                    scalar tk =
+                        Cmix_.value()
+                       *Foam::sqrt
+                       (
+                            muEff[i]/this->rho()()[i]/(epsilon[i] + e0.value())
+                       );
+
+                    // Chalmers PaSR model
+                    if (!useReactionRate_)
+                    {
+                        kappa_[i] =
+                            ( this->mesh().time().deltaTValue() + tc[i])
+                        /( this->mesh().time().deltaTValue() + tc[i] + tk);
+                    }
+                    else
+                    {
+                        kappa_[i] = tc[i]/(tc[i] + tk);
+                    }
+                }
+                else
+                {
+                    // Return to laminar combustion
+                    kappa_[i] = 1.0;
+                }
+            }
+        }
+        else
+        {
+            kappa_ = 1.0;
+        }
+    }
+}
+
+
+template<class CombThermoType>
+Foam::tmp<Foam::fvScalarMatrix>
+Foam::combustionModels::PaSR<CombThermoType>::R(const volScalarField& Y) const
+{
+
+    tmp<fvScalarMatrix> tSu
+    (
+        new fvScalarMatrix(Y, dimMass/dimTime)
+    );
+
+    fvScalarMatrix& Su = tSu();
+
+    if (this->active())
+    {
+        const label specieI = this->thermo().composition().species()[Y.name()];
+
+        Su += kappa_*this->pChemistry_->RR(specieI);
+    }
+
+    return tSu;
+}
+
+
+template<class CombThermoType>
+Foam::tmp<Foam::volScalarField>
+Foam::combustionModels::PaSR<CombThermoType>::dQ() const
+{
+    tmp<volScalarField> tdQ
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "dQ",
+                this->mesh().time().timeName(),
+                this->mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            this->mesh(),
+            dimensionedScalar("dQ", dimEnergy/dimTime, 0.0),
+            zeroGradientFvPatchScalarField::typeName
+        )
+    );
+
+    if (this->active())
+    {
+        volScalarField& dQ = tdQ();
+        dQ = kappa_*this->pChemistry_->dQ();
+    }
+
+    return tdQ;
+}
+
+
+template<class CombThermoType>
+Foam::tmp<Foam::volScalarField>
+Foam::combustionModels::PaSR<CombThermoType>::Sh() const
+{
+    tmp<volScalarField> tSh
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "Sh",
+                this->mesh().time().timeName(),
+                this->mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            this->mesh(),
+            dimensionedScalar("zero", dimEnergy/dimTime/dimVolume, 0.0),
+            zeroGradientFvPatchScalarField::typeName
+        )
+    );
+
+    if (this->active())
+    {
+        scalarField& Sh = tSh();
+        Sh = kappa_*this->pChemistry_->Sh();
+    }
+
+    return tSh;
+}
+
+
+template<class CombThermoType>
+bool Foam::combustionModels::PaSR<CombThermoType>::read()
+{
+    if (CombThermoType::read())
+    {
+        this->coeffs().lookup("Cmix") >> Cmix_;
+        this->coeffs().lookup("turbulentReaction") >> turbulentReaction_;
+        this->coeffs().lookup("useReactionRate") >> useReactionRate_;
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+// ************************************************************************* //

src/combustionModels/PaSR/PaSR.H

@@ -0,0 +1,140 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::combustionModels::PaSR
+
+Description
+    Simple infinitely fast chemistry combustion model based on the principle
+    mixed is burnt. Additional parameter C is used to distribute the heat
+    release rate.in time
+
+SourceFiles
+    PaSR.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef PaSR_H
+#define PaSR_H
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace combustionModels
+{
+
+/*---------------------------------------------------------------------------*\
+                            Class PaSR Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CombThermoType>
+class PaSR
+:
+    public CombThermoType
+{
+    // Private data
+
+        //- Mixing constant
+        dimensionedScalar Cmix_;
+
+        //- Turbulent reaction switch
+        Switch turbulentReaction_;
+
+        //- Mixing parameter
+        volScalarField kappa_;
+
+        //- Use reaction rate
+        bool useReactionRate_;
+
+
+    // Private Member Functions
+
+        //- Return the chemical time scale
+        tmp<volScalarField> tc() const;
+
+        //- Disallow copy construct
+        PaSR(const PaSR&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const PaSR&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("PaSR");
+
+
+    // Constructors
+
+        //- Construct from components
+        PaSR
+        (
+            const word& modelType,
+            const fvMesh& mesh
+        );
+
+
+    //- Destructor
+    virtual ~PaSR();
+
+
+    // Member Functions
+
+        // Evolution
+
+            //- Correct combustion rate
+            virtual void correct();
+
+            //- Fuel consumption rate matrix.
+            virtual tmp<fvScalarMatrix> R(const volScalarField& Y) const;
+
+            //- Heat release rate calculated from fuel consumption rate matrix
+            virtual tmp<volScalarField> dQ() const;
+
+            //-  Return source for enthalpy equation [kg/m/s3]
+            virtual tmp<volScalarField> Sh() const;
+
+    // I-O
+
+            //- Update properties from given dictionary
+            virtual bool read();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace combustionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "PaSR.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/combustionModels/PaSR/PaSRs.C

@@ -0,0 +1,53 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "makeCombustionTypes.H"
+
+#include "psiChemistryCombustionModel.H"
+#include "rhoChemistryCombustionModel.H"
+#include "PaSR.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace combustionModels
+{
+    makeCombustionTypes
+    (
+        PaSR,
+        psiChemistryCombustionModel
+    );
+
+    makeCombustionTypes
+    (
+        PaSR,
+        rhoChemistryCombustionModel
+    );
+}
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/combustionModels/combustionModel/combustionModel.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2009-2011 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2011-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -23,84 +23,71 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "combustionModel.H"
-#include "surfaceFields.H"
-#include "fvScalarMatrix.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
     defineTypeNameAndDebug(combustionModel, 0);
-    defineRunTimeSelectionTable(combustionModel, dictionary);
+}
-};
-
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::combustionModel::combustionModel
-(
-    const dictionary& combustionProps,
-    hsCombustionThermo& thermo,
-    const compressible::turbulenceModel& turbulence,
-    const surfaceScalarField& phi,
-    const volScalarField& rho
-)
-:
-    coeffs_(dictionary::null),
-    thermo_(thermo),
-    turbulence_(turbulence),
-    mesh_(phi.mesh()),
-    phi_(phi),
-    rho_(rho)
-{}
-
-
 Foam::combustionModel::combustionModel
 (
     const word& modelType,
-    const dictionary& combustionProps,
+    const fvMesh& mesh
-    hsCombustionThermo& thermo,
-    const compressible::turbulenceModel& turbulence,
-    const surfaceScalarField& phi,
-    const volScalarField& rho
 )
 :
-    coeffs_(combustionProps.subDict(modelType + "Coeffs")),
+    IOdictionary
-    thermo_(thermo),
+    (
-    turbulence_(turbulence),
+        IOobject
-    mesh_(phi.mesh()),
+        (
-    phi_(phi),
+            "combustionProperties",
-    rho_(rho)
+            mesh.time().constant(),
+            mesh,
+            IOobject::MUST_READ_IF_MODIFIED,
+            IOobject::NO_WRITE
+        )
+    ),
+    turbulencePtr_(),
+    mesh_(mesh),
+    active_(lookupOrDefault<Switch>("active", true)),
+    coeffs_(subDict(modelType + "Coeffs")),
+    modelType_(modelType)
 {}
 
 
 // * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * * //
 
 Foam::combustionModel::~combustionModel()
-{}
+{
+    if (turbulencePtr_)
+    {
+        turbulencePtr_ = 0;
+    }
+}
 
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
-void Foam::combustionModel::correct()
+
+bool Foam::combustionModel::read()
 {
-    // do nothing
+    if (regIOobject::read())
+    {
+        this->lookup("active") >> active_;
+        coeffs_ = subDict(modelType_ + "Coeffs");
+        return true;
+    }
+    else
+    {
+        return false;
+    }
 }
 
 
-Foam::tmp<Foam::fvScalarMatrix> Foam::combustionModel::R
+Foam::tmp<Foam::volScalarField> Foam::combustionModel::Sh() const
-(
-    volScalarField& Y
-) const
-{
-    return tmp<fvScalarMatrix>
-    (
-        new fvScalarMatrix(Y, dimMass/dimTime*Y.dimensions())
-    );
-}
-
-
-Foam::tmp<Foam::volScalarField> Foam::combustionModel::dQ() const
 {
     return tmp<Foam::volScalarField>
     (
@@ -108,7 +95,7 @@ Foam::tmp<Foam::volScalarField> Foam::combustionModel::dQ() const
         (
             IOobject
             (
-                "dQ",
+                "Sh",
                 mesh_.time().timeName(),
                 mesh_,
                 IOobject::NO_READ,
@@ -121,33 +108,4 @@ Foam::tmp<Foam::volScalarField> Foam::combustionModel::dQ() const
 }
 
 
-Foam::tmp<Foam::volScalarField> Foam::combustionModel::wFuelNorm() const
-{
-    return tmp<Foam::volScalarField>
-    (
-        new volScalarField
-        (
-            IOobject
-            (
-                "wFuelNorm",
-                mesh_.time().timeName(),
-                mesh_,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            mesh_,
-            dimensionedScalar("zero", dimMass/dimTime/pow3(dimLength), 0.0)
-        )
-    );
-}
-
-
-bool Foam::combustionModel::read(const dictionary& combustionProps)
-{
-    coeffs_ = combustionProps.subDict(type() + "Coeffs");
-
-    return true;
-}
-
-
 // ************************************************************************* //

src/combustionModels/combustionModel/combustionModel.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2009-2011 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2011-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -24,8 +24,7 @@ Class
     Foam::combustionModel
 
 Description
-    Base class for all non-premixed combustion models based on single step
+    Base class for combustion models
-    chemistry
 
 SourceFiles
     combustionModel.C
@@ -36,9 +35,7 @@ SourceFiles
 #define combustionModel_H
 
 #include "IOdictionary.H"
-#include "hsCombustionThermo.H"
 #include "turbulenceModel.H"
-#include "runTimeSelectionTables.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -50,29 +47,28 @@ namespace Foam
 \*---------------------------------------------------------------------------*/
 
 class combustionModel
+:
+    public IOdictionary
 {
 
 protected:
 
     // Protected data
 
-        //- Dictionary of coefficients for the particular model
-        dictionary coeffs_;
-
-        //- Reference to the thermodynamics
-        hsCombustionThermo& thermo_;
-
         //- Reference to the turbulence model
-        const compressible::turbulenceModel& turbulence_;
+        compressible::turbulenceModel* turbulencePtr_;
 
         //- Reference to the mesh database
         const fvMesh& mesh_;
 
-        //- Reference to mass-flux field
+        //- Active
-        const surfaceScalarField& phi_;
+        Switch active_;
 
-        //- Reference to the density field
+        //- Dictionary of the model
-        const volScalarField& rho_;
+        dictionary coeffs_;
+
+        //- Model name
+        const word modelType_;
 
 
 private:
@@ -92,65 +88,11 @@ public:
     TypeName("combustionModel");
 
 
-    // Declare run-time constructor selection table
-
-        declareRunTimeSelectionTable
-        (
-            autoPtr,
-            combustionModel,
-            dictionary,
-            (
-                const dictionary& combustionProperties,
-                hsCombustionThermo& thermo,
-                const compressible::turbulenceModel& turbulence,
-                const surfaceScalarField& phi,
-                const volScalarField& rho
-            ),
-            (
-                combustionProperties,
-                thermo,
-                turbulence,
-                phi,
-                rho
-            )
-        );
-
-
-    // Selectors
-
-        //- Return a reference to the selected combustion model
-        static autoPtr<combustionModel> New
-        (
-            const dictionary& combustionProperties,
-            hsCombustionThermo& thermo,
-            const compressible::turbulenceModel& turbulence,
-            const surfaceScalarField& phi,
-            const volScalarField& rho
-        );
-
-
     // Constructors
 
-        //- Construct null from components
-        combustionModel
-        (
-            const dictionary& combustionProps,
-            hsCombustionThermo& thermo,
-            const compressible::turbulenceModel& turbulence,
-            const surfaceScalarField& phi,
-            const volScalarField& rho
-        );
 
         //- Construct from components
-        combustionModel
+        combustionModel(const word& modelType, const fvMesh& mesh);
-        (
-            const word& modelType,
-            const dictionary& combustionProperties,
-            hsCombustionThermo& thermo,
-            const compressible::turbulenceModel& turbulence,
-            const surfaceScalarField& phi,
-            const volScalarField& rho
-        );
 
 
     //- Destructor
@@ -159,34 +101,53 @@ public:
 
     // Member Functions
 
-        // Access functions
+        // Access
 
-            //- Access combustion dictionary
+
-            const dictionary coeffs() const
+            //- Return const access to the mesh database
-            {
+            inline const fvMesh& mesh() const;
-                return coeffs_;
+
-            }
+            //- Return const access to phi
+            inline const surfaceScalarField& phi() const;
+
+            //- Return const access to rho
+            virtual tmp<volScalarField> rho() const = 0;
+
+            //- Return access to turbulence
+            inline const compressible::turbulenceModel& turbulence() const;
+
+            //- Set turbulence
+            inline void setTurbulence
+            (
+                compressible::turbulenceModel& turbModel
+            );
+
+            //- Is combustion active?
+            inline const Switch& active() const;
+
+            //- Return const dictionary of the model
+            inline const dictionary& coeffs() const;
 
 
     // Evolution
 
         //- Correct combustion rate
-        virtual void correct();
+        virtual void correct() = 0;
 
         //- Fuel consumption rate matrix, i.e. source term for fuel equation
-        virtual tmp<fvScalarMatrix> R(volScalarField& Y) const;
+        virtual tmp<fvScalarMatrix> R(const volScalarField& Y) const = 0;
 
         //- Heat release rate calculated from fuel consumption rate matrix
-        virtual tmp<volScalarField> dQ() const;
+        virtual tmp<volScalarField> dQ() const = 0;
 
-        //- Return normalised consumption rate of (fu - fres)
+        //-  Return source for enthalpy equation [kg/m/s3]
-        virtual tmp<Foam::volScalarField> wFuelNorm() const;
+        virtual tmp<volScalarField> Sh() const;
 
 
     // I-O
 
         //- Update properties from given dictionary
-        virtual bool read(const dictionary& combustionProps);
+        virtual bool read();
 };
 
 
@@ -196,6 +157,10 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+#include "combustionModelI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
 #endif
 
 // ************************************************************************* //

src/combustionModels/combustionModel/combustionModelI.H

@@ -0,0 +1,99 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+inline const Foam::fvMesh& Foam::combustionModel::mesh() const
+{
+    return mesh_;
+}
+
+
+inline const Foam::surfaceScalarField& Foam::combustionModel::phi() const
+{
+    if (turbulencePtr_)
+    {
+        return turbulencePtr_->phi();
+    }
+    else
+    {
+        FatalErrorIn
+        (
+            "const Foam::compressible::turbulenceModel& "
+            "Foam::combustionModel::turbulence() const "
+        )   << "turbulencePtr_ is empty. Please use "
+            << "combustionModel::setTurbulence "
+            << "(compressible::turbulenceModel& )"
+            << abort(FatalError);
+
+        return turbulencePtr_->phi();
+    }
+}
+
+
+inline const Foam::compressible::turbulenceModel&
+Foam::combustionModel::turbulence() const
+{
+    if (turbulencePtr_)
+    {
+        return *turbulencePtr_;
+    }
+    else
+    {
+        FatalErrorIn
+        (
+            "const Foam::compressible::turbulenceModel& "
+            "Foam::combustionModel::turbulence() const "
+        )   << "turbulencePtr_ is empty. Please use "
+            << "combustionModel::setTurbulence "
+            << "(compressible::turbulenceModel& )"
+            << abort(FatalError);
+
+        return *turbulencePtr_;
+    }
+}
+
+
+inline const Foam::Switch& Foam::combustionModel::active() const
+{
+    return active_;
+}
+
+
+inline void Foam::combustionModel::setTurbulence
+(
+    compressible::turbulenceModel& turbModel
+)
+{
+    turbulencePtr_ = &turbModel;
+}
+
+
+inline const Foam::dictionary& Foam::combustionModel::coeffs() const
+{
+    return coeffs_;
+}
+
+// ************************************************************************* //

src/combustionModels/combustionModel/makeCombustionTypes.H

@@ -0,0 +1,76 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef makeCombustionTypes_H
+#define makeCombustionTypes_H
+
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#define makeCombustionTypesThermo(CombModel, Comb, Thermo)                    \
+                                                                              \
+    typedef CombModel<Comb, Thermo> CombModel##Comb##Thermo;                  \
+                                                                              \
+                                                                              \
+    defineTemplateTypeNameAndDebugWithName                                    \
+    (                                                                         \
+        CombModel##Comb##Thermo,                                              \
+        #CombModel"<"#Comb","#Thermo">",                                      \
+        0                                                                     \
+    );                                                                        \
+                                                                              \
+                                                                              \
+    addToRunTimeSelectionTable                                                \
+    (                                                                         \
+        Comb,                                                                 \
+        CombModel##Comb##Thermo,                                              \
+        dictionary                                                            \
+    );
+
+#define makeCombustionTypes(CombModel, CombThermoType)                        \
+                                                                              \
+    typedef CombModel<CombThermoType>                                         \
+        CombModel##CombThermoType;                                            \
+                                                                              \
+    defineTemplateTypeNameAndDebugWithName                                    \
+    (                                                                         \
+        CombModel##CombThermoType,                                            \
+        #CombModel"<"#CombThermoType">",                                      \
+        0                                                                     \
+    );                                                                        \
+                                                                              \
+    addToRunTimeSelectionTable                                                \
+    (                                                                         \
+        CombThermoType,                                                       \
+        CombModel##CombThermoType,                                            \
+        dictionary                                                            \
+    );
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/combustionModels/infinitelyFastChemistry/infinitelyFastChemistry.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2009-2011 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2011-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -23,131 +23,80 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "infinitelyFastChemistry.H"
-#include "addToRunTimeSelectionTable.H"
-#include "fvmSup.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
 namespace combustionModels
 {
-    defineTypeNameAndDebug(infinitelyFastChemistry, 0);
-    addToRunTimeSelectionTable
-    (
-        combustionModel,
-        infinitelyFastChemistry,
-        dictionary
-    );
-};
-};
-
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::combustionModels::infinitelyFastChemistry::infinitelyFastChemistry
+template<class CombThermoType, class ThermoType>
+infinitelyFastChemistry<CombThermoType, ThermoType>::infinitelyFastChemistry
 (
-    const dictionary& combustionProps,
+    const word& modelType, const fvMesh& mesh
-    hsCombustionThermo& thermo,
-    const compressible::turbulenceModel& turbulence,
-    const surfaceScalarField& phi,
-    const volScalarField& rho
 )
 :
-    combustionModel(typeName, combustionProps, thermo, turbulence, phi, rho),
+    singleStepCombustion<CombThermoType, ThermoType>(modelType, mesh),
-    C_(readScalar(coeffs_.lookup("C"))),
+    C_(readScalar(this->coeffs().lookup("C")))
-    singleMixture_
-    (
-        dynamic_cast<singleStepReactingMixture<gasThermoPhysics>&>(thermo)
-    ),
-    wFuelNorm_
-    (
-        IOobject
-        (
-            "wFuelNorm",
-            mesh_.time().timeName(),
-            mesh_,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        mesh_,
-        dimensionedScalar("zero", dimMass/pow3(dimLength)/dimTime, 0.0)
-    )
 {}
 
 
 // * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * * //
 
-Foam::combustionModels::infinitelyFastChemistry::~infinitelyFastChemistry()
+template<class CombThermoType, class ThermoType>
+infinitelyFastChemistry<CombThermoType, ThermoType>::~infinitelyFastChemistry()
 {}
 
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
-void Foam::combustionModels::infinitelyFastChemistry::correct()
+template<class CombThermoType, class ThermoType>
+void infinitelyFastChemistry<CombThermoType, ThermoType>::correct()
 {
-    singleMixture_.fresCorrect();
+    this->wFuel_ ==
+        dimensionedScalar("zero", dimMass/pow3(dimLength)/dimTime, 0.0);
 
-    const label fuelI = singleMixture_.fuelIndex();
+    if (this->active())
-
-    const volScalarField& YFuel = thermo_.composition().Y()[fuelI];
-
-    const dimensionedScalar s = singleMixture_.s();
-
-    if (thermo_.composition().contains("O2"))
     {
-        const volScalarField& YO2 = thermo_.composition().Y("O2");
+        this->singleMixture_.fresCorrect();
-        wFuelNorm_ == rho_/(mesh_.time().deltaT()*C_)*min(YFuel, YO2/s.value());
+
+        const label fuelI = this->singleMixture_.fuelIndex();
+
+        const volScalarField& YFuel = this->thermo_->composition().Y()[fuelI];
+
+        const dimensionedScalar s = this->singleMixture_.s();
+
+        if (this->thermo_->composition().contains("O2"))
+        {
+            const volScalarField& YO2 = this->thermo_->composition().Y("O2");
+
+            this->wFuel_ ==
+                 this->rho()/(this->mesh().time().deltaT()*C_)
+                *min(YFuel, YO2/s.value());
+        }
     }
 }
 
 
-Foam::tmp<Foam::fvScalarMatrix>
+template<class CombThermoType, class ThermoType>
-Foam::combustionModels::infinitelyFastChemistry::R(volScalarField& Y) const
+bool infinitelyFastChemistry<CombThermoType, ThermoType>::read()
 {
-    const label specieI = thermo_.composition().species()[Y.name()];
+    if (singleStepCombustion<CombThermoType, ThermoType>::read())
-
+    {
-    const label fNorm = singleMixture_.specieProd()[specieI];
+        this->coeffs().lookup("C") >> C_ ;
-
-    const volScalarField fres(singleMixture_.fres(specieI));
-
-    const volScalarField wSpecie
-    (
-        wFuelNorm_*singleMixture_.specieStoichCoeffs()[specieI]
-      / max(fNorm*(Y - fres), scalar(0.001))
-    );
-
-    return -fNorm*wSpecie*fres + fNorm*fvm::Sp(wSpecie, Y);
-}
-
-
-Foam::tmp<Foam::volScalarField>
-Foam::combustionModels::infinitelyFastChemistry::dQ() const
-{
-    const label fuelI = singleMixture_.fuelIndex();
-    volScalarField& YFuel = thermo_.composition().Y(fuelI);
-
-    return -singleMixture_.qFuel()*(R(YFuel) & YFuel);
-}
-
-
-Foam::tmp<Foam::volScalarField>
-Foam::combustionModels::infinitelyFastChemistry::wFuelNorm() const
-{
-    return wFuelNorm_;
-}
-
-
-bool Foam::combustionModels::infinitelyFastChemistry::read
-(
-    const dictionary& combustionProps
-)
-{
-    combustionModel::read(combustionProps);
-    coeffs_.lookup("C") >> C_ ;
-
         return true;
+    }
+    else
+    {
+        return false;
+    }
 }
 
 
-// ************************************************************************* //
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace combustionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/combustionModels/infinitelyFastChemistry/infinitelyFastChemistry.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2009-2011 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2011-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -35,10 +35,7 @@ SourceFiles
 
 #ifndef infinitelyFastChemistry_H
 #define infinitelyFastChemistry_H
-
+#include "singleStepCombustion.H"
-#include "combustionModel.H"
-#include "singleStepReactingMixture.H"
-#include "thermoPhysicsTypes.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -51,21 +48,16 @@ namespace combustionModels
                   Class infinitelyFastChemistry Declaration
 \*---------------------------------------------------------------------------*/
 
+template<class CombThermoType, class ThermoType>
 class infinitelyFastChemistry
 :
-    public combustionModel
+    public singleStepCombustion <CombThermoType, ThermoType>
 {
     // Private data
 
         //- Model constant
         scalar C_;
 
-        //- Reference to singleStepReactingMixture mixture
-        singleStepReactingMixture<gasThermoPhysics>& singleMixture_;
-
-        //- Normalised consumption rate of (fu - fres)
-        volScalarField wFuelNorm_;
-
 
     // Private Member Functions
 
@@ -87,11 +79,7 @@ public:
         //- Construct from components
         infinitelyFastChemistry
         (
-            const dictionary& combustionProps,
+            const word& modelType, const fvMesh& mesh
-            hsCombustionThermo& thermo,
-            const compressible::turbulenceModel& turbulence,
-            const surfaceScalarField& phi,
-            const volScalarField& rho
         );
 
 
@@ -106,20 +94,11 @@ public:
             //- Correct combustion rate
             virtual void correct();
 
-            //- Fuel consumption rate matrix, i.e. source term for fuel equation
-            virtual tmp<fvScalarMatrix> R(volScalarField& Y) const;
-
-            //- Heat release rate calculated from fuel consumption rate matrix
-            virtual tmp<volScalarField> dQ() const;
-
-            //- Return normalised consumption rate of (fu - fres)
-            virtual tmp<volScalarField> wFuelNorm() const;
-
 
     // I-O
 
-            //- Update properties from given dictionary
+            //- Update properties
-            virtual bool read(const dictionary& combustionProperties);
+            virtual bool read();
 };
 
 
@@ -128,6 +107,13 @@ public:
 } // End namespace combustionModels
 } // End namespace Foam
 
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "infinitelyFastChemistry.C"
+#endif
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #endif

src/combustionModels/infinitelyFastChemistry/infinitelyFastChemistrys.C

@@ -0,0 +1,65 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "makeCombustionTypes.H"
+
+#include "thermoPhysicsTypes.H"
+#include "psiCombustionModel.H"
+#include "rhoCombustionModel.H"
+#include "infinitelyFastChemistry.H"
+#include "singleStepCombustion.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace combustionModels
+{
+
+    makeCombustionTypesThermo
+    (
+        infinitelyFastChemistry,
+        psiCombustionModel,
+        gasThermoPhysics
+    );
+
+    makeCombustionTypesThermo
+    (
+        infinitelyFastChemistry,
+        psiCombustionModel,
+        constGasThermoPhysics
+    );
+
+    makeCombustionTypesThermo
+    (
+        infinitelyFastChemistry,
+        rhoCombustionModel,
+        gasThermoPhysics
+    );
+}
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/combustionModels/psiChemistryCombustionModel/psiChemistryCombustionModel.C

@@ -0,0 +1,76 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "psiChemistryCombustionModel.H"
+
+/* * * * * * * * * * * * * * * private static data * * * * * * * * * * * * * */
+
+namespace Foam
+{
+namespace combustionModels
+{
+    defineTypeNameAndDebug(psiChemistryCombustionModel, 0);
+    defineRunTimeSelectionTable(psiChemistryCombustionModel, dictionary);
+
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+psiChemistryCombustionModel::psiChemistryCombustionModel
+(
+    const word& modelType,
+    const fvMesh& mesh
+)
+:
+    combustionModel(modelType, mesh),
+    pChemistry_(psiChemistryModel::New(mesh))
+{}
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+psiChemistryCombustionModel::~psiChemistryCombustionModel()
+{}
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+
+bool psiChemistryCombustionModel::read()
+{
+    if (combustionModel::read())
+    {
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+} // End namespace combustionModels
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/combustionModels/psiChemistryCombustionModel/psiChemistryCombustionModel.H

@@ -0,0 +1,160 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::psiChemistryCombustionModel
+
+Description
+    Combustion models for compressibility-based thermodynamics
+
+SourceFiles
+    psiChemistryCombustionModelI.H
+    psiChemistryCombustionModel.C
+    psiChemistryCombustionModelNew.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef psiChemistryCombustionModel_H
+#define psiChemistryCombustionModel_H
+
+#include "combustionModel.H"
+#include "autoPtr.H"
+#include "runTimeSelectionTables.H"
+#include "psiChemistryModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace combustionModels
+{
+
+/*---------------------------------------------------------------------------*\
+                     class psiChemistryCombustionModel Declaration
+\*---------------------------------------------------------------------------*/
+
+class psiChemistryCombustionModel
+:
+    public combustionModel
+{
+
+    // Private Member Functions
+
+        //- Construct as copy (not implemented)
+        psiChemistryCombustionModel(const psiChemistryCombustionModel&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const psiChemistryCombustionModel&);
+
+
+protected:
+
+    // Protected data
+
+        //- Auto pointer to psiChemistry
+        autoPtr<psiChemistryModel> pChemistry_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("psiChemistryCombustionModel");
+
+
+    //- Declare run-time constructor selection tables
+    declareRunTimeSelectionTable
+    (
+        autoPtr,
+        psiChemistryCombustionModel,
+        dictionary,
+        (
+            const word& modelType,
+            const fvMesh& mesh
+        ),
+        (modelType, mesh)
+    );
+
+
+    // Constructors
+
+
+        //- Construct from components and thermo
+        psiChemistryCombustionModel
+        (
+            const word& modelType,
+            const fvMesh& mesh
+        );
+
+
+
+    //- Selector
+    static autoPtr<psiChemistryCombustionModel> New
+    (
+        const fvMesh& mesh
+    );
+
+
+    //- Destructor
+    virtual ~psiChemistryCombustionModel();
+
+
+    // Member Functions
+
+
+        //- Return access to the thermo package
+        inline psiChemistryModel& pChemistry();
+
+        //- Return const access to the thermo package
+        inline const psiChemistryModel& pChemistry() const;
+
+        //- Return const access to rho
+        inline tmp<volScalarField> rho() const;
+
+        //- Return const access to rho
+        inline const hsCombustionThermo& thermo() const;
+
+        //- Return non const access to rho
+        inline hsCombustionThermo& thermo();
+
+
+    // I-O
+
+        //- Update properties from given dictionary
+        virtual bool read();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+} // End namespace combustionModels
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "psiChemistryCombustionModelI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/combustionModels/psiChemistryCombustionModel/psiChemistryCombustionModelI.H

@@ -0,0 +1,59 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+inline Foam::psiChemistryModel&
+Foam::combustionModels::psiChemistryCombustionModel::pChemistry()
+{
+    return pChemistry_();
+}
+
+inline const Foam::psiChemistryModel&
+Foam::combustionModels::psiChemistryCombustionModel::
+pChemistry() const
+{
+    return pChemistry_();
+}
+
+inline Foam::tmp<Foam::volScalarField>
+Foam::combustionModels::psiChemistryCombustionModel::rho() const
+{
+    return pChemistry_->thermo().rho();
+}
+
+inline const Foam::hsCombustionThermo&
+Foam::combustionModels::psiChemistryCombustionModel::thermo() const
+{
+    return pChemistry_->thermo();
+}
+
+inline Foam::hsCombustionThermo&
+Foam::combustionModels::psiChemistryCombustionModel::thermo()
+{
+    return pChemistry_->thermo();
+}
+
+// ************************************************************************* //

src/combustionModels/psiChemistryCombustionModel/psiChemistryCombustionModelNew.C

@@ -0,0 +1,78 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-2011 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "psiChemistryCombustionModel.H"
+
+// * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::combustionModels::psiChemistryCombustionModel>
+Foam::combustionModels::psiChemistryCombustionModel::New
+(
+    const fvMesh& mesh
+)
+{
+    const word combModelName
+    (
+        IOdictionary
+        (
+            IOobject
+            (
+                "combustionProperties",
+                mesh.time().constant(),
+                mesh,
+                IOobject::MUST_READ,
+                IOobject::NO_WRITE,
+                false
+            )
+        ).lookup("combustionModel")
+    );
+
+    Info<< "Selecting combustion model " << combModelName << endl;
+
+    dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(combModelName);
+
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    {
+        FatalErrorIn
+        (
+            "psiChemistryCombustionModel::New"
+        )   << "Unknown psiChemistryCombustionModel type "
+            << combModelName << endl << endl
+            << "Valid  combustionModels are : " << endl
+            << dictionaryConstructorTablePtr_->toc()
+            << exit(FatalError);
+    }
+
+    const label tempOpen = combModelName.find('<');
+
+    const word className = combModelName(0, tempOpen);
+
+    return autoPtr<psiChemistryCombustionModel>
+        (cstrIter()(className, mesh));
+}
+
+
+// ************************************************************************* //

src/combustionModels/psiCombustionModel/psiCombustionModel.C

@@ -0,0 +1,78 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "psiCombustionModel.H"
+
+/* * * * * * * * * * * * * * * private static data * * * * * * * * * * * * * */
+
+namespace Foam
+{
+namespace combustionModels
+{
+    defineTypeNameAndDebug(psiCombustionModel, 0);
+    defineRunTimeSelectionTable(psiCombustionModel, dictionary);
+}
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+
+Foam::combustionModels::psiCombustionModel::psiCombustionModel
+(
+    const word& modelType,
+    const fvMesh& mesh
+)
+:
+    combustionModel(modelType, mesh),
+    thermo_(hsCombustionThermo::New(mesh))
+{}
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::combustionModels::psiCombustionModel::~psiCombustionModel()
+{}
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+
+bool Foam::combustionModels::psiCombustionModel::read()
+{
+    if (combustionModel::read())
+    {
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+Foam::tmp<Foam::volScalarField>
+Foam::combustionModels::psiCombustionModel::rho() const
+{
+    return thermo_->rho();
+}
+
+// ************************************************************************* //

src/combustionModels/psiCombustionModel/psiCombustionModel.H

@@ -0,0 +1,155 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::psiCombustionModel
+
+Description
+    Combustion models for compressibility-based thermodynamics
+
+SourceFiles
+    psiCombustionModelI.H
+    psiCombustionModel.C
+    psiCombustionModelNew.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef psiCombustionModel_H
+#define psiCombustionModel_H
+
+#include "combustionModel.H"
+#include "autoPtr.H"
+#include "runTimeSelectionTables.H"
+#include "hsCombustionThermo.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace combustionModels
+{
+
+/*---------------------------------------------------------------------------*\
+                     class psiCombustionModel Declaration
+\*---------------------------------------------------------------------------*/
+
+class psiCombustionModel
+:
+    public combustionModel
+{
+
+    // Private Member Functions
+
+        //- Construct as copy (not implemented)
+        psiCombustionModel(const psiCombustionModel&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const psiCombustionModel&);
+
+
+protected:
+
+    // Protected data
+
+
+        //- Thermo package
+        autoPtr<hsCombustionThermo> thermo_;
+
+
+
+
+public:
+
+    //- Runtime type information
+    TypeName("psiCombustionModel");
+
+
+    //- Declare run-time constructor selection tables
+    declareRunTimeSelectionTable
+    (
+        autoPtr,
+        psiCombustionModel,
+        dictionary,
+        (
+            const word& modelType,
+            const fvMesh& mesh
+        ),
+        (modelType, mesh)
+    );
+
+
+    // Constructors
+
+
+        //- Construct from components and thermo
+        psiCombustionModel
+        (
+            const word& modelType,
+            const fvMesh& mesh
+        );
+
+
+
+    //- Selector
+    static autoPtr<psiCombustionModel> New
+    (
+        const fvMesh& mesh
+    );
+
+
+    //- Destructor
+    virtual ~psiCombustionModel();
+
+
+    // Member Functions
+
+        //- Return access to the thermo package
+        inline hsCombustionThermo& thermo();
+
+        //- Return const access to the thermo package
+        inline const hsCombustionThermo& thermo() const;
+
+        //- Return tmp of rho
+        virtual tmp<volScalarField> rho() const;
+
+    // I-O
+
+        //- Update properties from given dictionary
+        virtual bool read();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+} // End namespace combustionModels
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "psiCombustionModelI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //