Merge commit 'OpenCFD/master' into olesenm

applications/solvers/combustion/dieselEngineFoam/createFields.H

@@ -6,7 +6,7 @@ autoPtr<psiChemistryModel> pChemistry
 );
 psiChemistryModel& chemistry = pChemistry();
 
-hCombustionThermo& thermo = chemistry.thermo();
+hsCombustionThermo& thermo = chemistry.thermo();
 
 basicMultiComponentMixture& composition = thermo.composition();
 PtrList<volScalarField>& Y = composition.Y();
@@ -50,7 +50,7 @@ volVectorField U
 volScalarField& p = thermo.p();
 const volScalarField& psi = thermo.psi();
 const volScalarField& T = thermo.T();
-volScalarField& h = thermo.h();
+volScalarField& hs = thermo.hs();
 
 
 #include "compressibleCreatePhi.H"
@@ -92,4 +92,18 @@ forAll(Y, i)
 {
     fields.add(Y[i]);
 }
-fields.add(h);
+fields.add(hs);
+
+DimensionedField<scalar, volMesh> chemistrySh
+(
+    IOobject
+    (
+        "chemistry::Sh",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("chemistrySh", dimEnergy/dimTime/dimVolume, 0.0)
+);

applications/solvers/combustion/dieselEngineFoam/dieselEngineFoam.C

@@ -23,7 +23,7 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Application
-    dieselFoam
+    dieselEngineFoam
 
 Description
     Solver for diesel engine spray and combustion.
@@ -103,13 +103,15 @@ int main(int argc, char *argv[])
             kappa = (runTime.deltaT() + tc)/(runTime.deltaT() + tc + tk);
         }
 
+        chemistrySh = kappa*chemistry.Sh()();
+
         #include "rhoEqn.H"
         #include "UEqn.H"
 
         for (label ocorr=1; ocorr <= nOuterCorr; ocorr++)
         {
             #include "YEqn.H"
-            #include "hEqn.H"
+            #include "hsEqn.H"
 
             // --- PISO loop
             for (int corr=1; corr<=nCorr; corr++)

applications/solvers/combustion/dieselEngineFoam/hEqn.H

@@ -1,13 +0,0 @@
-{
-    solve
-    (
-        fvm::ddt(rho, h)
-      + mvConvection->fvmDiv(phi, h)
-      - fvm::laplacian(turbulence->alphaEff(), h)
-     ==
-       DpDt
-     + dieselSpray.heatTransferSource()
-    );
-
-    thermo.correct();
-}

applications/solvers/combustion/dieselEngineFoam/hsEqn.H

@@ -0,0 +1,14 @@
+{
+    solve
+    (
+        fvm::ddt(rho, hs)
+      + mvConvection->fvmDiv(phi, hs)
+      - fvm::laplacian(turbulence->alphaEff(), hs)
+     ==
+       DpDt
+     + dieselSpray.heatTransferSource()().dimensionedInternalField()
+     + chemistrySh
+    );
+
+    thermo.correct();
+}

applications/solvers/combustion/dieselFoam/dieselFoam.C

@@ -100,7 +100,7 @@ int main(int argc, char *argv[])
         for (label ocorr=1; ocorr <= nOuterCorr; ocorr++)
         {
             #include "YEqn.H"
-            #include "hEqn.H"
+            #include "hsEqn.H"
 
             // --- PISO loop
             for (int corr=1; corr<=nCorr; corr++)

applications/solvers/combustion/fireFoam/Allwclean

@@ -3,7 +3,6 @@ cd ${0%/*} || exit 1    # run from this directory
 set -x
 
 wclean libso combustionModels
-wclean all sensibleEnthalpyCombustionThermophysicalModels
 wclean
 
 # ----------------------------------------------------------------- end-of-file

applications/solvers/combustion/fireFoam/Allwmake

@@ -2,7 +2,6 @@
 cd ${0%/*} || exit 1    # run from this directory
 set -x
 
-sensibleEnthalpyCombustionThermophysicalModels/Allwmake
 wmake libso combustionModels
 wmake
 

applications/solvers/combustion/fireFoam/Make/options

@@ -1,7 +1,5 @@
 EXE_INC = \
     -I./combustionModels/lnInclude \
-    -I./sensibleEnthalpyCombustionThermophysicalModels/basic/lnInclude \
-    -I./sensibleEnthalpyCombustionThermophysicalModels/reactionThermo/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
@@ -11,8 +9,6 @@ EXE_INC = \
 
 
 EXE_LIBS = \
-    -lbasicSensibleEnthalpyThermophysicalModels \
-    -lreactionSensibleEnthalpyThermophysicalModels \
     -lcombustionModels \
     -lspecie \
     -lreactionThermophysicalModels \

applications/solvers/combustion/fireFoam/sensibleEnthalpyCombustionThermophysicalModels/Allwmake

@@ -1,8 +0,0 @@
-#!/bin/sh
-cd ${0%/*} || exit 1    # run from this directory
-set -x
-
-wmake libso basic
-wmake libso reactionThermo
-
-# ----------------------------------------------------------------- end-of-file

applications/solvers/combustion/fireFoam/sensibleEnthalpyCombustionThermophysicalModels/basic/Make/files

@@ -1,4 +0,0 @@
-psiThermo/basicSensiblePsiThermo/basicSensiblePsiThermo.C
-psiThermo/basicSensiblePsiThermo/newBasicSensiblePsiThermo.C
-
-LIB = $(FOAM_LIBBIN)/libbasicSensibleEnthalpyThermophysicalModels

applications/solvers/combustion/fireFoam/sensibleEnthalpyCombustionThermophysicalModels/basic/Make/options

@@ -1,7 +0,0 @@
-EXE_INC = \
-    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-    -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude
-
-LIB_LIBS = \
-    -lfiniteVolume

applications/solvers/combustion/fireFoam/sensibleEnthalpyCombustionThermophysicalModels/basic/psiThermo/basicSensiblePsiThermo/makeBasicSensiblePsiThermo.H

@@ -1,66 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 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
-
-InClass
-    Foam::basicSensiblePsiThermo
-
-Description
-    Macros for creating 'basic' compresibility-based thermo packages
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef makeBasicSensiblePsiThermo_H
-#define makeBasicSensiblePsiThermo_H
-
-#include "basicSensiblePsiThermo.H"
-#include "addToRunTimeSelectionTable.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#define makeBasicSensiblePsiThermo(Cthermo,Mixture,Transport,Thermo,EqnOfState)       \
-                                                                              \
-typedef Cthermo<Mixture<Transport<specieThermo<Thermo<EqnOfState> > > > >     \
-    Cthermo##Mixture##Transport##Thermo##EqnOfState;                          \
-                                                                              \
-defineTemplateTypeNameAndDebugWithName                                        \
-(                                                                             \
-    Cthermo##Mixture##Transport##Thermo##EqnOfState,                          \
-    #Cthermo                                                                  \
-        "<"#Mixture"<"#Transport"<specieThermo<"#Thermo"<"#EqnOfState">>>>>", \
-    0                                                                         \
-);                                                                            \
-                                                                              \
-addToRunTimeSelectionTable                                                    \
-(                                                                             \
-    basicSensiblePsiThermo,                                                           \
-    Cthermo##Mixture##Transport##Thermo##EqnOfState,                          \
-    fvMesh                                                                    \
-)
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

applications/solvers/combustion/fireFoam/sensibleEnthalpyCombustionThermophysicalModels/basic/psiThermo/basicSensiblePsiThermo/newBasicSensiblePsiThermo.C

@@ -1,75 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 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 "basicSensiblePsiThermo.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-Foam::autoPtr<Foam::basicSensiblePsiThermo> Foam::basicSensiblePsiThermo::New
-(
-    const fvMesh& mesh
-)
-{
-    word thermoTypeName;
-
-    // Enclose the creation of the thermophysicalProperties to ensure it is
-    // deleted before the turbulenceModel is created otherwise the dictionary
-    // is entered in the database twice
-    {
-        IOdictionary thermoDict
-        (
-            IOobject
-            (
-                "thermophysicalProperties",
-                mesh.time().constant(),
-                mesh,
-                IOobject::MUST_READ,
-                IOobject::NO_WRITE
-            )
-        );
-
-        thermoDict.lookup("thermoType") >> thermoTypeName;
-    }
-
-    Info<< "Selecting thermodynamics package " << thermoTypeName << endl;
-
-    fvMeshConstructorTable::iterator cstrIter =
-        fvMeshConstructorTablePtr_->find(thermoTypeName);
-
-    if (cstrIter == fvMeshConstructorTablePtr_->end())
-    {
-        FatalErrorIn("basicSensiblePsiThermo::New(const fvMesh&)")
-            << "Unknown basicSensiblePsiThermo type " << thermoTypeName << nl << nl
-            << "Valid basicSensiblePsiThermo types are:" << nl
-            << fvMeshConstructorTablePtr_->toc() << nl
-            << exit(FatalError);
-    }
-
-    return autoPtr<basicSensiblePsiThermo>(cstrIter()(mesh));
-}
-
-
-// ************************************************************************* //

applications/solvers/combustion/fireFoam/sensibleEnthalpyCombustionThermophysicalModels/reactionThermo/Make/files

@@ -1,5 +0,0 @@
-combustionThermo/hsCombustionThermo/hsCombustionThermo.C
-combustionThermo/hsCombustionThermo/newhsCombustionThermo.C
-combustionThermo/hsCombustionThermo/hsCombustionThermos.C
-
-LIB = $(FOAM_LIBBIN)/libreactionSensibleEnthalpyThermophysicalModels

applications/solvers/combustion/fireFoam/sensibleEnthalpyCombustionThermophysicalModels/reactionThermo/Make/options

@@ -1,9 +0,0 @@
-EXE_INC = \
-    -I../basic/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
-    -I$(LIB_SRC)/finiteVolume/lnInclude
-
-LIB_LIBS = \
-    -lfiniteVolume

applications/solvers/combustion/reactingFoam/Make/options

@@ -1,5 +1,4 @@
 EXE_INC = \
-    -I../XiFoam \
     -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \

applications/solvers/combustion/reactingFoam/UEqn.H

@@ -0,0 +1,15 @@
+    fvVectorMatrix UEqn
+    (
+        fvm::ddt(rho, U)
+      + fvm::div(phi, U)
+      + turbulence->divDevRhoReff(U)
+     ==
+        rho*g
+    );
+
+    UEqn.relax();
+
+    if (momentumPredictor)
+    {
+        solve(UEqn == -fvc::grad(p));
+    }

applications/solvers/combustion/reactingFoam/chemistry.H

@@ -21,4 +21,6 @@
     {
         kappa = 1.0;
     }
+
+    chemistrySh = kappa*chemistry.Sh()();
 }

applications/solvers/combustion/reactingFoam/createFields.H

@@ -5,7 +5,7 @@ autoPtr<psiChemistryModel> pChemistry
 );
 psiChemistryModel& chemistry = pChemistry();
 
-hCombustionThermo& thermo = chemistry.thermo();
+hsCombustionThermo& thermo = chemistry.thermo();
 
 basicMultiComponentMixture& composition = thermo.composition();
 PtrList<volScalarField>& Y = composition.Y();
@@ -40,8 +40,8 @@ volVectorField U
 
 volScalarField& p = thermo.p();
 const volScalarField& psi = thermo.psi();
-volScalarField& h = thermo.h();
+volScalarField& hs = thermo.hs();
-
+const volScalarField& T = thermo.T();
 
 #include "compressibleCreatePhi.H"
 
@@ -81,5 +81,18 @@ forAll(Y, i)
 {
     fields.add(Y[i]);
 }
-fields.add(h);
+fields.add(hs);
 
+DimensionedField<scalar, volMesh> chemistrySh
+(
+    IOobject
+    (
+        "chemistry::Sh",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("chemistrySh", dimEnergy/dimTime/dimVolume, 0.0)
+);

applications/solvers/combustion/reactingFoam/hsEqn.H

@@ -0,0 +1,20 @@
+{
+    fvScalarMatrix hsEqn
+    (
+        fvm::ddt(rho, hs)
+      + mvConvection->fvmDiv(phi, hs)
+      - fvm::laplacian(turbulence->alphaEff(), hs)
+//      - fvm::laplacian(turbulence->muEff(), hs)  // unit lewis no.
+     ==
+        DpDt
+      + chemistrySh
+    );
+
+    hsEqn.relax();
+    hsEqn.solve();
+
+    thermo.correct();
+
+    Info<< "T gas min/max   = " << min(T).value() << ", "
+        << max(T).value() << endl;
+}

applications/solvers/combustion/reactingFoam/pEqn.H

@@ -0,0 +1,68 @@
+rho = thermo.rho();
+
+volScalarField rUA = 1.0/UEqn.A();
+U = rUA*UEqn.H();
+
+if (transonic)
+{
+    surfaceScalarField phid
+    (
+        "phid",
+        fvc::interpolate(psi)
+       *(
+            (fvc::interpolate(U) & mesh.Sf())
+          + fvc::ddtPhiCorr(rUA, rho, U, phi)
+        )
+    );
+
+    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    {
+        fvScalarMatrix pEqn
+        (
+            fvm::ddt(psi, p)
+          + fvm::div(phid, p)
+          - fvm::laplacian(rho*rUA, p)
+        );
+
+        pEqn.solve();
+
+        if (nonOrth == nNonOrthCorr)
+        {
+            phi == pEqn.flux();
+        }
+    }
+}
+else
+{
+    phi =
+        fvc::interpolate(rho)
+       *(
+            (fvc::interpolate(U) & mesh.Sf())
+          + fvc::ddtPhiCorr(rUA, rho, U, phi)
+        );
+
+    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    {
+        fvScalarMatrix pEqn
+        (
+            fvm::ddt(psi, p)
+          + fvc::div(phi)
+          - fvm::laplacian(rho*rUA, p)
+        );
+
+        pEqn.solve();
+
+        if (nonOrth == nNonOrthCorr)
+        {
+            phi += pEqn.flux();
+        }
+    }
+}
+
+#include "rhoEqn.H"
+#include "compressibleContinuityErrs.H"
+
+U -= rUA*fvc::grad(p);
+U.correctBoundaryConditions();
+
+DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);

applications/solvers/combustion/reactingFoam/reactingFoam.C

@@ -73,9 +73,7 @@ int main(int argc, char *argv[])
         {
             #include "UEqn.H"
             #include "YEqn.H"
-
+            #include "hsEqn.H"
-            #define Db turbulence->alphaEff()
-            #include "hEqn.H"
 
             // --- PISO loop
             for (int corr=1; corr<=nCorr; corr++)

applications/solvers/combustion/rhoReactingFoam/Make/options

@@ -1,5 +1,4 @@
 EXE_INC = \
-    -I../XiFoam \
     -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \

applications/solvers/combustion/rhoReactingFoam/UEqn.H

@@ -0,0 +1,15 @@
+    fvVectorMatrix UEqn
+    (
+        fvm::ddt(rho, U)
+      + fvm::div(phi, U)
+      + turbulence->divDevRhoReff(U)
+     ==
+        rho*g
+    );
+
+    UEqn.relax();
+
+    if (momentumPredictor)
+    {
+        solve(UEqn == -fvc::grad(p));
+    }

applications/solvers/combustion/rhoReactingFoam/chemistry.H

@@ -21,4 +21,6 @@
     {
         kappa = 1.0;
     }
+
+    chemistrySh = kappa*chemistry.Sh()();
 }

applications/solvers/combustion/rhoReactingFoam/createFields.H

@@ -5,7 +5,7 @@ autoPtr<rhoChemistryModel> pChemistry
 );
 rhoChemistryModel& chemistry = pChemistry();
 
-hReactionThermo& thermo = chemistry.thermo();
+hsReactionThermo& thermo = chemistry.thermo();
 
 basicMultiComponentMixture& composition = thermo.composition();
 PtrList<volScalarField>& Y = composition.Y();
@@ -40,7 +40,8 @@ volVectorField U
 
 volScalarField& p = thermo.p();
 const volScalarField& psi = thermo.psi();
-volScalarField& h = thermo.h();
+volScalarField& hs = thermo.hs();
+const volScalarField& T = thermo.T();
 
 
 #include "compressibleCreatePhi.H"
@@ -81,5 +82,18 @@ forAll(Y, i)
 {
     fields.add(Y[i]);
 }
-fields.add(h);
+fields.add(hs);
 
+DimensionedField<scalar, volMesh> chemistrySh
+(
+    IOobject
+    (
+        "chemistry::Sh",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("chemistrySh", dimEnergy/dimTime/dimVolume, 0.0)
+);

applications/solvers/combustion/rhoReactingFoam/hsEqn.H

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

applications/solvers/combustion/rhoReactingFoam/rhoReactingFoam.C

@@ -74,7 +74,7 @@ int main(int argc, char *argv[])
         {
             #include "UEqn.H"
             #include "YEqn.H"
-            #include "hEqn.H"
+            #include "hsEqn.H"
 
             // --- PISO loop
             for (int corr=1; corr<=nCorr; corr++)

applications/solvers/lagrangian/coalChemistryFoam/chemistry.H

@@ -22,4 +22,6 @@
     {
         kappa = 1.0;
     }
+
+    chemistrySh = kappa*chemistry.Sh()();
 }

applications/solvers/lagrangian/coalChemistryFoam/coalChemistryFoam.C

@@ -84,16 +84,8 @@ int main(int argc, char *argv[])
 
         coalParcels.evolve();
 
-        coalParcels.info();
-
-        Info<< endl;
-
         limestoneParcels.evolve();
 
-        limestoneParcels.info();
-
-        Info<< endl;
-
         #include "chemistry.H"
         #include "rhoEqn.H"
 
@@ -102,16 +94,13 @@ int main(int argc, char *argv[])
         {
             #include "UEqn.H"
             #include "YEqn.H"
-            #include "hEqn.H"
+            #include "hsEqn.H"
 
             // --- PISO loop
             for (int corr=1; corr<=nCorr; corr++)
             {
                 #include "pEqn.H"
             }
-
-            Info<< "T gas min/max   = " << min(T).value() << ", "
-                << max(T).value() << endl;
         }
 
         turbulence->correct();

applications/solvers/lagrangian/coalChemistryFoam/createFields.H

@@ -6,7 +6,7 @@
     );
     psiChemistryModel& chemistry = pChemistry();
 
-    hCombustionThermo& thermo = chemistry.thermo();
+    hsCombustionThermo& thermo = chemistry.thermo();
 
     basicMultiComponentMixture& composition = thermo.composition();
     PtrList<volScalarField>& Y = composition.Y();
@@ -22,7 +22,7 @@
     }
 
     volScalarField& p = thermo.p();
-    volScalarField& h = thermo.h();
+    volScalarField& hs = thermo.hs();
     const volScalarField& T = thermo.T();
     const volScalarField& psi = thermo.psi();
 
@@ -32,7 +32,7 @@
     {
         fields.add(Y[i]);
     }
-    fields.add(h);
+    fields.add(hs);
 
     volScalarField rho
     (
@@ -133,5 +133,19 @@
         "energy",
         mesh,
         dimEnergy/dimTime/dimVolume,
-        "h"
+        "hs"
+    );
+
+    DimensionedField<scalar, volMesh> chemistrySh
+    (
+        IOobject
+        (
+            "chemistry::Sh",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("chemistrySh", dimEnergy/dimTime/dimVolume, 0.0)
     );

applications/solvers/lagrangian/coalChemistryFoam/hEqn.H

@@ -1,22 +0,0 @@
-{
-    fvScalarMatrix hEqn
-    (
-        fvm::ddt(rho, h)
-      + mvConvection->fvmDiv(phi, h)
-      - fvm::laplacian(turbulence->alphaEff(), h)
-     ==
-        DpDt
-     +  coalParcels.Sh()
-     +  limestoneParcels.Sh()
-     +  enthalpySource.Su()
-     +  radiation->Sh(thermo)
-    );
-
-    hEqn.relax();
-
-    hEqn.solve();
-
-    thermo.correct();
-
-    radiation->correct();
-}

applications/solvers/lagrangian/coalChemistryFoam/hsEqn.H

@@ -0,0 +1,26 @@
+{
+    fvScalarMatrix hsEqn
+    (
+        fvm::ddt(rho, hs)
+      + mvConvection->fvmDiv(phi, hs)
+      - fvm::laplacian(turbulence->alphaEff(), hs)
+     ==
+        DpDt
+     +  coalParcels.Sh()
+     +  limestoneParcels.Sh()
+     +  enthalpySource.Su()
+     +  radiation->Shs(thermo)
+     +  chemistrySh
+    );
+
+    hsEqn.relax();
+
+    hsEqn.solve();
+
+    thermo.correct();
+
+    radiation->correct();
+
+    Info<< "T gas min/max   = " << min(T).value() << ", "
+        << max(T).value() << endl;
+}

applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/chemistry.H

@@ -22,4 +22,6 @@
     {
         kappa = 1.0;
     }
+
+    chemistrySh = kappa*chemistry.Sh()();
 }

applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/createFields.H

@@ -6,7 +6,7 @@
     );
     rhoChemistryModel& chemistry = pChemistry();
 
-    hReactionThermo& thermo = chemistry.thermo();
+    hsReactionThermo& thermo = chemistry.thermo();
 
     basicMultiComponentMixture& composition = thermo.composition();
     PtrList<volScalarField>& Y = composition.Y();
@@ -22,7 +22,7 @@
     }
 
     volScalarField& p = thermo.p();
-    volScalarField& h = thermo.h();
+    volScalarField& hs = thermo.hs();
     const volScalarField& T = thermo.T();
     const volScalarField& psi = thermo.psi();
 
@@ -88,4 +88,18 @@
     {
         fields.add(Y[i]);
     }
-    fields.add(h);
+    fields.add(hs);
+
+    DimensionedField<scalar, volMesh> chemistrySh
+    (
+        IOobject
+        (
+            "chemistry::Sh",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("chemistrySh", dimEnergy/dimTime/dimVolume, 0.0)
+    );

applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/hsEqn.H

@@ -32,14 +32,15 @@
     {
         solve
         (
-            fvm::ddt(rho, h)
+            fvm::ddt(rho, hs)
-          + mvConvection->fvmDiv(phi, h)
+          + mvConvection->fvmDiv(phi, hs)
-          - fvm::laplacian(turbulence->alphaEff(), h)
+          - fvm::laplacian(turbulence->alphaEff(), hs)
          ==
             pWork()
           + parcels.Sh()
-          + radiation->Sh(thermo)
+          + radiation->Shs(thermo)
           + energySource.Su()
+          + chemistrySh
         );
 
         thermo.correct();

applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/porousExplicitSourceReactingParcelFoam.C

@@ -33,7 +33,7 @@ Description
     The solver includes:
     - reacting multiphase parcel cloud
     - porous media
-    - point mass sources
+    - mass, momentum and energy sources
     - polynomial based, incompressible thermodynamics (f(T))
 
     Note: ddtPhiCorr not used here when porous zones are active
@@ -89,13 +89,11 @@ int main(int argc, char *argv[])
 
         parcels.evolve();
 
-        parcels.info();
-
         #include "chemistry.H"
         #include "rhoEqn.H"
         #include "UEqn.H"
         #include "YEqn.H"
-        #include "hEqn.H"
+        #include "hsEqn.H"
 
         // --- PISO loop
         for (int corr=0; corr<nCorr; corr++)

applications/solvers/lagrangian/reactingParcelFoam/chemistry.H

@@ -22,4 +22,6 @@
     {
         kappa = 1.0;
     }
+
+    chemistrySh = kappa*chemistry.Sh()();
 }

applications/solvers/lagrangian/reactingParcelFoam/createFields.H

@@ -6,7 +6,7 @@
     );
     psiChemistryModel& chemistry = pChemistry();
 
-    hCombustionThermo& thermo = chemistry.thermo();
+    hsCombustionThermo& thermo = chemistry.thermo();
 
     basicMultiComponentMixture& composition = thermo.composition();
     PtrList<volScalarField>& Y = composition.Y();
@@ -22,7 +22,7 @@
     }
 
     volScalarField& p = thermo.p();
-    volScalarField& h = thermo.h();
+    volScalarField& hs = thermo.hs();
     const volScalarField& T = thermo.T();
     const volScalarField& psi = thermo.psi();
 
@@ -94,4 +94,18 @@
     {
         fields.add(Y[i]);
     }
-    fields.add(h);
+    fields.add(hs);
+
+    DimensionedField<scalar, volMesh> chemistrySh
+    (
+        IOobject
+        (
+            "chemistry::Sh",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("chemistry::Sh", dimEnergy/dimTime/dimVolume, 0.0)
+    );

applications/solvers/lagrangian/reactingParcelFoam/hEqn.H

@@ -1,20 +0,0 @@
-{
-    fvScalarMatrix hEqn
-    (
-        fvm::ddt(rho, h)
-      + mvConvection->fvmDiv(phi, h)
-      - fvm::laplacian(turbulence->alphaEff(), h)
-     ==
-        DpDt
-     +  parcels.Sh()
-     +  radiation->Sh(thermo)
-    );
-
-    hEqn.relax();
-
-    hEqn.solve();
-
-    thermo.correct();
-
-    radiation->correct();
-}

applications/solvers/lagrangian/reactingParcelFoam/hsEqn.H

@@ -0,0 +1,24 @@
+{
+    fvScalarMatrix hEqn
+    (
+        fvm::ddt(rho, hs)
+      + mvConvection->fvmDiv(phi, hs)
+      - fvm::laplacian(turbulence->alphaEff(), hs)
+     ==
+        DpDt
+     +  parcels.Sh()
+     +  radiation->Shs(thermo)
+     +  chemistrySh
+    );
+
+    hEqn.relax();
+
+    hEqn.solve();
+
+    thermo.correct();
+
+    radiation->correct();
+
+    Info<< "T gas min/max   = " << min(T).value() << ", "
+        << max(T).value() << endl;
+}

applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C

@@ -74,8 +74,6 @@ int main(int argc, char *argv[])
 
         parcels.evolve();
 
-        parcels.info();
-
         #include "chemistry.H"
         #include "rhoEqn.H"
 
@@ -88,12 +86,9 @@ int main(int argc, char *argv[])
             // --- PISO loop
             for (int corr=1; corr<=nCorr; corr++)
             {
-                #include "hEqn.H"
+                #include "hsEqn.H"
                 #include "pEqn.H"
             }
-
-            Info<< "T gas min/max   = " << min(T).value() << ", "
-                << max(T).value() << endl;
         }
 
         turbulence->correct();

applications/solvers/lagrangian/uncoupledKinematicParcelFoam/uncoupledKinematicParcelFoam.C

@@ -61,7 +61,6 @@ int main(int argc, char *argv[])
 
         Info<< "Evolving " << kinematicCloud.name() << endl;
         kinematicCloud.evolve();
-        kinematicCloud.info();
 
         runTime.write();
 

applications/utilities/parallelProcessing/reconstructPar/reconstructPar.C

@@ -270,11 +270,31 @@ int main(int argc, char *argv[])
                 procMeshes.boundaryProcAddressing()
             );
 
-            pointReconstructor.reconstructFields<scalar>(objects);
+            pointReconstructor.reconstructFields<scalar>
-            pointReconstructor.reconstructFields<vector>(objects);
+            (
-            pointReconstructor.reconstructFields<sphericalTensor>(objects);
+                objects,
-            pointReconstructor.reconstructFields<symmTensor>(objects);
+                selectedFields
-            pointReconstructor.reconstructFields<tensor>(objects);
+            );
+            pointReconstructor.reconstructFields<vector>
+            (
+                objects,
+                selectedFields
+            );
+            pointReconstructor.reconstructFields<sphericalTensor>
+            (
+                objects,
+                selectedFields
+            );
+            pointReconstructor.reconstructFields<symmTensor>
+            (
+                objects,
+                selectedFields
+            );
+            pointReconstructor.reconstructFields<tensor>
+            (
+                objects,
+                selectedFields
+            );
         }
         else
         {

applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightMesh.C

@@ -30,13 +30,10 @@ License
 #include "fvMesh.H"
 #include "globalMeshData.H"
 #include "PstreamCombineReduceOps.H"
-#include "processorPolyPatch.H"
 #include "cellModeller.H"
 #include "IOmanip.H"
 #include "itoa.H"
 #include "ensightWriteBinary.H"
-#include "globalIndex.H"
-#include "PackedBoolList.H"
 #include "mapDistribute.H"
 
 #include <fstream>
@@ -1111,14 +1108,11 @@ void Foam::ensightMesh::writeAllFacePrimsBinary
             writeEnsDataBinary(key,ensightGeometryFile);
             writeEnsDataBinary(nPrims,ensightGeometryFile);
 
-            if (&prims != NULL)
-            {
             writeFacePrimsBinary
             (
                 UIndirectList<face>(patchFaces, prims)(),
                 ensightGeometryFile
             );
-            }
 
             for (int slave=1; slave<Pstream::nProcs(); slave++)
             {
@@ -1348,17 +1342,18 @@ void Foam::ensightMesh::writeAscii
 
                 // Renumber the patch points/faces into unique points
                 labelList pointToGlobal;
-                labelList uniquePointLabels;
+                labelList uniqueMeshPointLabels;
                 autoPtr<globalIndex> globalPointsPtr =
                 mesh_.globalData().mergePoints
                 (
                     p.meshPoints(),
                     p.meshPointMap(),
                     pointToGlobal,
-                    uniquePointLabels
+                    uniqueMeshPointLabels
                 );
 
-                pointField uniquePoints(p.localPoints(), uniquePointLabels);
+
+                pointField uniquePoints(mesh_.points(), uniqueMeshPointLabels);
                 // Renumber the patch faces
                 faceList patchFaces(p.localFaces());
                 forAll(patchFaces, i)
@@ -1601,16 +1596,17 @@ void Foam::ensightMesh::writeBinary
 
                 // Renumber the patch points/faces into unique points
                 labelList pointToGlobal;
-                labelList uniquePointLabels;
+                labelList uniqueMeshPointLabels;
                 autoPtr<globalIndex> globalPointsPtr =
                 mesh_.globalData().mergePoints
                 (
                     p.meshPoints(),
                     p.meshPointMap(),
                     pointToGlobal,
-                    uniquePointLabels
+                    uniqueMeshPointLabels
                 );
-                pointField uniquePoints(p.localPoints(), uniquePointLabels);
+
+                pointField uniquePoints(mesh_.points(), uniqueMeshPointLabels);
                 // Renumber the patch faces
                 faceList patchFaces(p.localFaces());
                 forAll(patchFaces, i)

src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.C

@@ -1705,7 +1705,7 @@ Foam::autoPtr<Foam::globalIndex> Foam::globalMeshData::mergePoints
     const labelList& meshPoints,
     const Map<label>& meshPointMap,
     labelList& pointToGlobal,
-    labelList& uniquePoints
+    labelList& uniqueMeshPoints
 ) const
 {
     const indirectPrimitivePatch& cpp = coupledPatch();
@@ -1713,125 +1713,175 @@ Foam::autoPtr<Foam::globalIndex> Foam::globalMeshData::mergePoints
     const mapDistribute& pointSlavesMap = globalPointSlavesMap();
 
 
-    // 1. Count number of masters on my processor.
+    // The patch points come in two variants:
-    label nCoupledMaster = 0;
+    // - not on a coupled patch so guaranteed unique
-    label nCoupledSlave = 0;
+    // - on a coupled patch
-    PackedBoolList isMaster(meshPoints.size(), 1);
+    // If the point is on a coupled patch the problem is that the
+    // master-slave structure (globalPointSlaves etc.) assigns one of the
+    // coupled points to be the master but this master point is not
+    // necessarily on the patch itself! (it might just be connected to the
+    // patch point via coupled patches).
+    // So this means that all master point loops should be over the
+    // master-slave structure, not over the patch points and that the unique
+    // point returned is a mesh point.
+    // (unless we want to do the whole master-slave analysis again for the
+    //  current patch only).
 
-    forAll(meshPoints, localPointI)
+
+    // Determine mapping from coupled point to patch point and vice versa
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    labelList patchToCoupled(meshPoints.size(), -1);
+    label nCoupled = 0;
+    labelList coupledToPatch(pointSlavesMap.constructSize(), -1);
+
+    // Note: loop over patch since usually smaller
+    forAll(meshPoints, patchPointI)
     {
-        label meshPointI = meshPoints[localPointI];
+        label meshPointI = meshPoints[patchPointI];
 
         Map<label>::const_iterator iter = cpp.meshPointMap().find(meshPointI);
 
         if (iter != cpp.meshPointMap().end())
         {
-            // My localPointI is a coupled point.
+            patchToCoupled[patchPointI] = iter();
-
+            coupledToPatch[iter()] = patchPointI;
-            label coupledPointI = iter();
+            nCoupled++;
-
-            if (pointSlaves[coupledPointI].size() > 0)
-            {
-                nCoupledMaster++;
-            }
-            else
-            {
-                isMaster[localPointI] = 0;
-                nCoupledSlave++;
-            }
-        }
-    }
-
-    label myUniquePoints = meshPoints.size() - nCoupledSlave;
-
-    Pout<< "Points :" << nl
-        << "    patch            : " << meshPoints.size() << nl
-        << "    of which coupled : " << nCoupledMaster+nCoupledSlave << nl
-        << "    of which master  : " << nCoupledMaster << nl
-        << "    of which slave   : " << nCoupledSlave << nl
-        << endl;
-
-
-    // 2. Create global indexing for unique points.
-    autoPtr<globalIndex> globalPointsPtr(new globalIndex(myUniquePoints));
-
-
-    // 3. Assign global point numbers. Keep slaves unset.
-    pointToGlobal.setSize(meshPoints.size());
-    pointToGlobal = -1;
-    uniquePoints.setSize(myUniquePoints);
-    label nMaster = 0;
-
-    forAll(isMaster, localPointI)
-    {
-        if (isMaster[localPointI])
-        {
-            pointToGlobal[localPointI] = globalPointsPtr().toGlobal(nMaster);
-            uniquePoints[nMaster] = localPointI;
-            nMaster++;
         }
     }
 
+    //Pout<< "Patch:" << nl
+    //    << "    points:" << meshPoints.size() << nl
+    //    << "    of which on coupled patch:" << nCoupled << endl;
 
-    // 4. Push global index for coupled points to slaves.
+
+    // Pull coupled-to-patch information to master
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    pointSlavesMap.distribute(coupledToPatch);
+
+
+    // Check on master whether point is anywhere on patch
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // List of master points that are on the patch
+    DynamicList<label> masterPoints(pointSlaves.size());
+
+    forAll(pointSlaves, coupledPointI)
     {
-        labelList masterToGlobal(pointSlavesMap.constructSize(), -1);
-
-        forAll(meshPoints, localPointI)
-        {
-            label meshPointI = meshPoints[localPointI];
-
-            Map<label>::const_iterator iter = cpp.meshPointMap().find
-            (
-                meshPointI
-            );
-
-            if (iter != cpp.meshPointMap().end())
-            {
-                // My localPointI is a coupled point.
-                label coupledPointI = iter();
-
         const labelList& slaves = pointSlaves[coupledPointI];
 
         if (slaves.size() > 0)
         {
-                    // Duplicate master globalpoint into slave slots
+            // I am master. Is this point on the patch on myself or on any
-                    masterToGlobal[coupledPointI] = pointToGlobal[localPointI];
+            // any slave?
+            if (coupledToPatch[coupledPointI] != -1)
+            {
+                masterPoints.append(coupledPointI);
+            }
+            else
+            {
                 forAll(slaves, i)
                 {
-                        masterToGlobal[slaves[i]] = pointToGlobal[localPointI];
+                    if (coupledToPatch[slaves[i]] != -1)
+                    {
+                        masterPoints.append(coupledPointI);
+                        break;
+                    }
                 }
             }
         }
     }
 
-        // Send back
+
+    // Create global indexing
+    // ~~~~~~~~~~~~~~~~~~~~~~
+    // 1. patch points that are not on coupled patch:
+    //      meshPoints.size()-nCoupled
+    // 2. master points that are on patch:
+    //      masterPoints.size()
+    label myUniquePoints = meshPoints.size()-nCoupled+masterPoints.size();
+    autoPtr<globalIndex> globalPointsPtr(new globalIndex(myUniquePoints));
+
+    //Pout<< "CoupledPatch:" << nl
+    //    << "    points:" << cpp.nPoints() << nl
+    //    << "    of which on patch:" << masterPoints.size() << endl;
+
+
+    // Allocate unique points
+    // ~~~~~~~~~~~~~~~~~~~~~~
+
+    pointToGlobal.setSize(meshPoints.size());
+    pointToGlobal = -1;
+    uniqueMeshPoints.setSize(myUniquePoints);
+
+    // Allocate globals for uncoupled patch points
+    label nMaster = 0;
+    forAll(patchToCoupled, patchPointI)
+    {
+        if (patchToCoupled[patchPointI] == -1)
+        {
+            // Allocate global point
+            label globalI = globalPointsPtr().toGlobal(nMaster);
+            pointToGlobal[patchPointI] = globalI;
+            uniqueMeshPoints[nMaster] = meshPoints[patchPointI];
+            nMaster++;
+        }
+    }
+
+    // Allocate globals for master 
+    labelList masterToGlobal(pointSlavesMap.constructSize(), -456);
+
+    forAll(masterPoints, i)
+    {
+        label coupledPointI = masterPoints[i];
+
+        // Allocate global point
+        label globalI = globalPointsPtr().toGlobal(nMaster);
+        if (coupledToPatch[coupledPointI] != -1)
+        {
+            pointToGlobal[coupledToPatch[coupledPointI]] = globalI;
+        }
+        uniqueMeshPoints[nMaster] = cpp.meshPoints()[coupledPointI];
+        nMaster++;
+
+        // Put global into slave slots
+        const labelList& slaves = pointSlaves[coupledPointI];
+        masterToGlobal[coupledPointI] = globalI;    // not really necessary
+        forAll(slaves, i)
+        {
+            masterToGlobal[slaves[i]] = globalI;
+        }
+    }
+
+
+    if (nMaster != myUniquePoints)
+    {
+        FatalErrorIn("globalMeshData::mergePoints(..)")
+            << "problem." << abort(FatalError);
+    }
+
+
+    // Send back (from slave slots) to originating processor
     pointSlavesMap.reverseDistribute(cpp.nPoints(), masterToGlobal);
 
-        // On slave copy master index into overal map.
+    // On slaves take over global number
-        forAll(meshPoints, localPointI)
+    forAll(patchToCoupled, patchPointI)
     {
-            label meshPointI = meshPoints[localPointI];
+        label coupledPointI = patchToCoupled[patchPointI];
 
-            Map<label>::const_iterator iter = cpp.meshPointMap().find
+        if (coupledPointI != -1)
-            (
-                meshPointI
-            );
-
-            if (iter != cpp.meshPointMap().end())
         {
-                // My localPointI is a coupled point.
-                label coupledPointI = iter();
             const labelList& slaves = pointSlaves[coupledPointI];
 
             if (slaves.size() == 0)
             {
-                    pointToGlobal[localPointI] = masterToGlobal[coupledPointI];
+                pointToGlobal[patchPointI] = masterToGlobal[coupledPointI];
-                }
             }
         }
     }
 
+
     return globalPointsPtr;
 }
 

src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.H

@@ -560,13 +560,17 @@ public:
                 ) const;
 
                 //- Helper for merging patch point data. Takes maps from
-                //  local points to/from mesh
+                //  local points to/from mesh. Determines
+                //  - my unique points. These are mesh points, not patch points
+                //    since the master might not be on the patch.
+                //  - global numbering over all unique indices.
+                //  - the global number for all local points.
                 autoPtr<globalIndex> mergePoints
                 (
                     const labelList& meshPoints,
                     const Map<label>& meshPointMap,
                     labelList& pointToGlobal,
-                    labelList& uniquePoints
+                    labelList& uniqueMeshPoints
                 ) const;
 
 

src/lagrangian/dieselSpray/parcel/parcel.C

@@ -247,7 +247,7 @@ bool Foam::parcel::move(spray& sDB)
         {
             oMass[i] = m()*oYf[i];
             label j = sDB.liquidToGasIndex()[i];
-            oHg += oYf[i]*sDB.gasProperties()[j].H(T());
+            oHg += oYf[i]*sDB.gasProperties()[j].Hs(T());
         }
 
         vector oMom = m()*U();
@@ -273,7 +273,7 @@ bool Foam::parcel::move(spray& sDB)
         {
             nMass[i] = m()*nYf[i];
             label j = sDB.liquidToGasIndex()[i];
-            nHg += nYf[i]*sDB.gasProperties()[j].H(T());
+            nHg += nYf[i]*sDB.gasProperties()[j].Hs(T());
         }
 
         vector nMom = m()*U();
@@ -288,9 +288,7 @@ bool Foam::parcel::move(spray& sDB)
         }
         sDB.sms()[celli] += oMom - nMom;
 
-        sDB.shs()[celli]   +=
+        sDB.shs()[celli] += oTotMass*(oH + oPE) - m()*(nH + nPE);
-            oTotMass*(oH + oPE)
-          - m()*(nH + nPE);
 
         // Remove evaporated mass from stripped mass
         ms() -= ms()*(oTotMass-m())/oTotMass;
@@ -446,7 +444,7 @@ void Foam::parcel::updateParcelProperties
     for (label i=0; i<Nf; i++)
     {
         label j = sDB.liquidToGasIndex()[i];
-        oldhg += Yf0[i]*sDB.gasProperties()[j].H(T());
+        oldhg += Yf0[i]*sDB.gasProperties()[j].Hs(T());
     }
 
     scalar oldhv = fuels.hl(pg, T(), X());
@@ -478,7 +476,7 @@ void Foam::parcel::updateParcelProperties
             for (label i=0; i<Nf; i++)
             {
                 label j = sDB.liquidToGasIndex()[i];
-                newhg += Ynew[i]*sDB.gasProperties()[j].H(Tnew);
+                newhg += Ynew[i]*sDB.gasProperties()[j].Hs(Tnew);
             }
 
             newhv = fuels.hl(pg, Tnew, X());
@@ -616,7 +614,8 @@ void Foam::parcel::updateParcelProperties
                     {
                         if (n>100)
                         {
-                            Info<< "n = " << n << ", T = " << Td << ", pv = " << pAtSurface << endl;
+                            Info<< "n = " << n << ", T = " << Td << ", pv = "
+                                << pAtSurface << endl;
                         }
                     }
                 }

src/lagrangian/dieselSpray/parcel/parcel.H

@@ -185,7 +185,8 @@ public:
             //- Return the names of the liquid components
             inline const List<word>& liquidNames() const;
 
-            //- Return the names of the liquid fuel components - identical with liquidNames
+            //- Return the names of the liquid fuel components
+            //  - identical with liquidNames
             inline const List<word>& fuelNames() const;
 
             //- Return diameter of droplets in parcel

src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C

@@ -260,6 +260,9 @@ void Foam::KinematicCloud<ParcelType>::evolve()
         evolveCloud();
 
         postEvolve();
+
+        info();
+        Info<< endl;
     }
 }
 

src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.C

@@ -272,6 +272,9 @@ void Foam::ReactingCloud<ParcelType>::evolve()
         evolveCloud();
 
         postEvolve();
+
+        info();
+        Info<< endl;
     }
 }
 

src/lagrangian/intermediate/clouds/Templates/ReactingMultiphaseCloud/ReactingMultiphaseCloud.C

@@ -228,6 +228,9 @@ void Foam::ReactingMultiphaseCloud<ParcelType>::evolve()
         evolveCloud();
 
         postEvolve();
+
+        info();
+        Info<< endl;
     }
 }
 

src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C

@@ -160,20 +160,6 @@ Foam::ThermoCloud<ParcelType>::ThermoCloud
         ),
         this->mesh(),
         dimensionedScalar("zero", dimEnergy, 0.0)
-    ),
-    hcTrans_
-    (
-        IOobject
-        (
-            this->name() + "hcTrans",
-            this->db().time().timeName(),
-            this->db(),
-            IOobject::NO_READ,
-            IOobject::NO_WRITE,
-            false
-        ),
-        this->mesh(),
-        dimensionedScalar("zero", dimEnergy, 0.0)
     )
 {
     if (readFields)
@@ -220,7 +206,6 @@ void Foam::ThermoCloud<ParcelType>::resetSourceTerms()
 {
     KinematicCloud<ParcelType>::resetSourceTerms();
     hsTrans_.field() = 0.0;
-    hcTrans_.field() = 0.0;
 }
 
 
@@ -234,6 +219,9 @@ void Foam::ThermoCloud<ParcelType>::evolve()
         evolveCloud();
 
         postEvolve();
+
+        info();
+        Info<< endl;
     }
 }
 

src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.H

@@ -108,15 +108,9 @@ protected:
 
         // Sources
 
-            //- Sensible enthalpy transfer
+            //- Sensible enthalpy transfer [J/kg]
             DimensionedField<scalar, volMesh> hsTrans_;
 
-            //- Chemical enthalpy transfer
-            // - If solving for total enthalpy, the carrier phase enthalpy will
-            //   receive the full enthalpy of reaction via creation of reaction
-            //   products
-            DimensionedField<scalar, volMesh> hcTrans_;
-
 
     // Protected member functions
 
@@ -194,19 +188,10 @@ public:
 
                 // Enthalpy
 
-                    //- Return reference to sensible enthalpy source
+                    //- Sensible enthalpy transfer [J/kg]
                     inline DimensionedField<scalar, volMesh>& hsTrans();
 
-                    //- Return tmp total sensible enthalpy source term
+                    //- Return enthalpy source [J/kg/m3/s]
-                    inline tmp<DimensionedField<scalar, volMesh> > Shs() const;
-
-                    //- Return reference to chemical enthalpy source
-                    inline DimensionedField<scalar, volMesh>& hcTrans();
-
-                    //- Return tmp chemical enthalpy source term
-                    inline tmp<DimensionedField<scalar, volMesh> > Shc() const;
-
-                    //- Return tmp total enthalpy source term
                     inline tmp<DimensionedField<scalar, volMesh> > Sh() const;
 
 

src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloudI.H

@@ -85,82 +85,6 @@ Foam::ThermoCloud<ParcelType>::hsTrans()
 }
 
 
-template<class ParcelType>
-inline Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >
-Foam::ThermoCloud<ParcelType>::Shs() const
-{
-    tmp<DimensionedField<scalar, volMesh> > tShs
-    (
-        new DimensionedField<scalar, volMesh>
-        (
-            IOobject
-            (
-                this->name() + "Shs",
-                this->db().time().timeName(),
-                this->mesh(),
-                IOobject::NO_READ,
-                IOobject::AUTO_WRITE,
-                false
-            ),
-            this->mesh(),
-            dimensionedScalar
-            (
-                "zero",
-                dimMass/dimLength/pow3(dimTime),
-                0.0
-            )
-        )
-    );
-
-    scalarField& Shs = tShs().field();
-    Shs = hsTrans_/(this->mesh().V()*this->db().time().deltaT());
-
-    return tShs;
-}
-
-
-template<class ParcelType>
-inline Foam::DimensionedField<Foam::scalar, Foam::volMesh>&
-Foam::ThermoCloud<ParcelType>::hcTrans()
-{
-    return hcTrans_;
-}
-
-
-template<class ParcelType>
-inline Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >
-Foam::ThermoCloud<ParcelType>::Shc() const
-{
-    tmp<DimensionedField<scalar, volMesh> > tShc
-    (
-        new DimensionedField<scalar, volMesh>
-        (
-            IOobject
-            (
-                this->name() + "Shc",
-                this->db().time().timeName(),
-                this->mesh(),
-                IOobject::NO_READ,
-                IOobject::AUTO_WRITE,
-                false
-            ),
-            this->mesh(),
-            dimensionedScalar
-            (
-                "zero",
-                dimMass/dimLength/pow3(dimTime),
-                0.0
-            )
-        )
-    );
-
-    scalarField& Shc = tShc().field();
-    Shc = hcTrans_/(this->mesh().V()*this->db().time().deltaT());
-
-    return tShc;
-}
-
-
 template<class ParcelType>
 inline Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >
 Foam::ThermoCloud<ParcelType>::Sh() const
@@ -178,19 +102,10 @@ Foam::ThermoCloud<ParcelType>::Sh() const
                 IOobject::AUTO_WRITE,
                 false
             ),
-            this->mesh(),
+            hsTrans_/(this->mesh().V()*this->db().time().deltaT())
-            dimensionedScalar
-            (
-                "zero",
-                dimMass/dimLength/pow3(dimTime),
-                0.0
-            )
         )
     );
 
-    scalarField& Sh = tSh().field();
-    Sh = (hsTrans_ + hcTrans_)/(this->mesh().V()*this->db().time().deltaT());
-
     return tSh;
 }
 
@@ -205,7 +120,7 @@ Foam::ThermoCloud<ParcelType>::Ep() const
         (
             IOobject
             (
-                this->name() + "radiationEp",
+                this->name() + "radiation::Ep",
                 this->db().time().timeName(),
                 this->db(),
                 IOobject::NO_READ,
@@ -248,7 +163,7 @@ Foam::ThermoCloud<ParcelType>::ap() const
         (
             IOobject
             (
-                this->name() + "radiationAp",
+                this->name() + "radiation::ap",
                 this->db().time().timeName(),
                 this->db(),
                 IOobject::NO_READ,
@@ -291,7 +206,7 @@ Foam::ThermoCloud<ParcelType>::sigmap() const
         (
             IOobject
             (
-                this->name() + "radiationSigmap",
+                this->name() + "radiation::sigmap",
                 this->db().time().timeName(),
                 this->db(),
                 IOobject::NO_READ,

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

@@ -267,7 +267,6 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
         YLiquid_,
         dMassPC,
         Sh,
-        dhsTrans,
         Ne,
         NCpW,
         Cs
@@ -296,7 +295,6 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
         canCombust_,
         dMassDV,
         Sh,
-        dhsTrans,
         Ne,
         NCpW,
         Cs
@@ -398,19 +396,11 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
         {
             label gid = td.cloud().composition().localToGlobalCarrierId(GAS, i);
             td.cloud().rhoTrans(gid)[cellI] += np0*dMassGas[i];
-            td.cloud().hcTrans()[cellI] +=
-                np0
-               *dMassGas[i]
-               *td.cloud().mcCarrierThermo().speciesData()[gid].Hc();
         }
         forAll(YLiquid_, i)
         {
             label gid = td.cloud().composition().localToGlobalCarrierId(LIQ, i);
             td.cloud().rhoTrans(gid)[cellI] += np0*dMassLiquid[i];
-            td.cloud().hcTrans()[cellI] +=
-                np0
-               *dMassLiquid[i]
-               *td.cloud().mcCarrierThermo().speciesData()[gid].Hc();
         }
 /*
         // No mapping between solid components and carrier phase
@@ -418,19 +408,11 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
         {
             label gid = td.cloud().composition().localToGlobalCarrierId(SLD, i);
             td.cloud().rhoTrans(gid)[cellI] += np0*dMassSolid[i];
-            td.cloud().hcTrans()[cellI] +=
-                np0
-               *dMassSolid[i]
-               *td.cloud().mcCarrierThermo().speciesData()[gid].Hc();
         }
 */
         forAll(dMassSRCarrier, i)
         {
             td.cloud().rhoTrans(i)[cellI] += np0*dMassSRCarrier[i];
-            td.cloud().hcTrans()[cellI] +=
-                np0
-               *dMassSRCarrier[i]
-               *td.cloud().mcCarrierThermo().speciesData()[i].Hc();
         }
 
         // Update momentum transfer
@@ -476,7 +458,7 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calc
 */
             td.cloud().UTrans()[cellI] += np0*mass1*U1;
             td.cloud().hsTrans()[cellI] +=
-                np0*mass1*HEff(td, pc, T1, idG, idL, idS);
+                np0*mass1*HEff(td, pc, T1, idG, idL, idS); // using total h
         }
     }
 
@@ -520,7 +502,6 @@ void Foam::ReactingMultiphaseParcel<ParcelType>::calcDevolatilisation
     bool& canCombust,
     scalarField& dMassDV,
     scalar& Sh,
-    scalar& dhsTrans,
     scalar& N,
     scalar& NCpW,
     scalarField& Cs

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

@@ -238,7 +238,6 @@ protected:
             bool& canCombust,          // 'can combust' flag
             scalarField& dMassDV,      // mass transfer - local to particle
             scalar& Sh,                // explicit particle enthalpy source
-            scalar& dhsTrans,          // sensible enthalpy transfer to carrier
             scalar& N,                 // flux of species emitted from particle
             scalar& NCpW,              // sum of N*Cp*W of emission species
             scalarField& Cs            // carrier conc. of emission species

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

@@ -287,7 +287,6 @@ void Foam::ReactingParcel<ParcelType>::calc
         Y_,
         dMassPC,
         Sh,
-        dhsTrans,
         Ne,
         NCpW,
         Cs
@@ -341,10 +340,6 @@ void Foam::ReactingParcel<ParcelType>::calc
         {
             label gid = td.cloud().composition().localToGlobalCarrierId(0, i);
             td.cloud().rhoTrans(gid)[cellI] += np0*dMassPC[i];
-            td.cloud().hcTrans()[cellI] +=
-                np0
-               *dMassPC[i]
-               *td.cloud().mcCarrierThermo().speciesData()[gid].Hc();
         }
 
         // Update momentum transfer
@@ -371,7 +366,7 @@ void Foam::ReactingParcel<ParcelType>::calc
                 td.cloud().rhoTrans(gid)[cellI] += np0*mass1*Y_[i];
             }
             td.cloud().UTrans()[cellI] += np0*mass1*U1;
-            td.cloud().hcTrans()[cellI] +=
+            td.cloud().hsTrans()[cellI] +=
                 np0*mass1*td.cloud().composition().H(0, Y_, pc_, T1);
         }
     }
@@ -417,7 +412,6 @@ void Foam::ReactingParcel<ParcelType>::calcPhaseChange
     const scalarField& YComponents,
     scalarField& dMassPC,
     scalar& Sh,
-    scalar& dhsTrans,               // TODO: not used
     scalar& N,
     scalar& NCpW,
     scalarField& Cs
@@ -469,6 +463,7 @@ void Foam::ReactingParcel<ParcelType>::calcPhaseChange
             td.cloud().composition().localToGlobalCarrierId(idPhase, i);
         const label idl = td.cloud().composition().globalIds(idPhase)[i];
 
+        // Calculate enthalpy transfer
         if
         (
             td.cloud().phaseChange().enthalpyTransfer()

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

@@ -209,7 +209,6 @@ protected:
             const scalarField& YComponents, // component mass fractions
             scalarField& dMassPC,      // mass transfer - local to particle
             scalar& Sh,                // explicit particle enthalpy source
-            scalar& dhsTrans,          // sensible enthalpy transfer to carrier
             scalar& N,                 // flux of species emitted from particle
             scalar& NCpW,              // sum of N*Cp*W of emission species
             scalarField& Cs            // carrier conc. of emission species

src/lagrangian/intermediate/submodels/Reacting/CompositionModel/CompositionModel/CompositionModel.C

@@ -410,6 +410,130 @@ Foam::scalar Foam::CompositionModel<CloudType>::H
 }
 
 
+template<class CloudType>
+Foam::scalar Foam::CompositionModel<CloudType>::Hs
+(
+    const label phaseI,
+    const scalarField& Y,
+    const scalar p,
+    const scalar T
+) const
+{
+    const phaseProperties& props = phaseProps_[phaseI];
+    scalar HsMixture = 0.0;
+    switch (props.phase())
+    {
+        case phaseProperties::GAS:
+        {
+            forAll(Y, i)
+            {
+                label gid = props.globalIds()[i];
+                HsMixture += Y[i]*mcCarrierThermo_.speciesData()[gid].Hs(T);
+            }
+            break;
+        }
+        case phaseProperties::LIQUID:
+        {
+            forAll(Y, i)
+            {
+                label gid = props.globalIds()[i];
+                HsMixture +=
+                    Y[i]
+                   *(
+                       this->liquids().properties()[gid].h(p, T)
+                     - this->liquids().properties()[gid].h(p, 298.25)
+                    );
+            }
+            break;
+        }
+        case phaseProperties::SOLID:
+        {
+            forAll(Y, i)
+            {
+                label gid = props.globalIds()[i];
+                HsMixture += Y[i]*this->solids().properties()[gid].cp()*T;
+            }
+            break;
+        }
+        default:
+        {
+            FatalErrorIn
+            (
+                "Foam::scalar Foam::CompositionModel<CloudType>::Hs"
+                "("
+                "    const label, "
+                "    const scalarField&, "
+                "    const scalar, "
+                "    const scalar"
+                ") const"
+            )   << "Unknown phase enumeration" << nl << abort(FatalError);
+        }
+    }
+
+    return HsMixture;
+}
+
+
+template<class CloudType>
+Foam::scalar Foam::CompositionModel<CloudType>::Hc
+(
+    const label phaseI,
+    const scalarField& Y,
+    const scalar p,
+    const scalar T
+) const
+{
+    const phaseProperties& props = phaseProps_[phaseI];
+    scalar HcMixture = 0.0;
+    switch (props.phase())
+    {
+        case phaseProperties::GAS:
+        {
+            forAll(Y, i)
+            {
+                label gid = props.globalIds()[i];
+                HcMixture += Y[i]*mcCarrierThermo_.speciesData()[gid].Hc();
+            }
+            break;
+        }
+        case phaseProperties::LIQUID:
+        {
+            forAll(Y, i)
+            {
+                label gid = props.globalIds()[i];
+                HcMixture +=
+                    Y[i]*this->liquids().properties()[gid].h(p, 298.15);
+            }
+            break;
+        }
+        case phaseProperties::SOLID:
+        {
+            forAll(Y, i)
+            {
+                label gid = props.globalIds()[i];
+                HcMixture += Y[i]*this->solids().properties()[gid].Hf();
+            }
+            break;
+        }
+        default:
+        {
+            FatalErrorIn
+            (
+                "Foam::scalar Foam::CompositionModel<CloudType>::Hc"
+                "("
+                "    const label, "
+                "    const scalarField&, "
+                "    const scalar, "
+                "    const scalar"
+                ") const"
+            )   << "Unknown phase enumeration" << nl << abort(FatalError);
+        }
+    }
+
+    return HcMixture;
+}
+
+
 template<class CloudType>
 Foam::scalar Foam::CompositionModel<CloudType>::cp
 (

src/lagrangian/intermediate/submodels/Reacting/CompositionModel/CompositionModel/CompositionModel.H

@@ -226,7 +226,7 @@ public:
 
         // Evaluation
 
-            //- Return enthalpy for the phase phaseI
+            //- Return total enthalpy for the phase phaseI
             virtual scalar H
             (
                 const label phaseI,
@@ -235,6 +235,24 @@ public:
                 const scalar T
             ) const;
 
+            //- Return sensible enthalpy for the phase phaseI
+            virtual scalar Hs
+            (
+                const label phaseI,
+                const scalarField& Y,
+                const scalar p,
+                const scalar T
+            ) const;
+
+            //- Return chemical enthalpy for the phase phaseI
+            virtual scalar Hc
+            (
+                const label phaseI,
+                const scalarField& Y,
+                const scalar p,
+                const scalar T
+            ) const;
+
             //- Return specific heat caoacity for the phase phaseI
             virtual scalar cp
             (

src/parallel/reconstruct/pointFieldReconstructor.H

@@ -138,7 +138,11 @@ public:
 
         //- Reconstruct and write all fields
         template<class Type>
-        void reconstructFields(const IOobjectList& objects);
+        void reconstructFields
+        (
+            const IOobjectList& objects,
+            const HashSet<word>& selectedFields
+        );
 };
 
 

src/parallel/reconstruct/pointFieldReconstructorReconstructFields.C

@@ -143,7 +143,8 @@ Foam::pointFieldReconstructor::reconstructField(const IOobject& fieldIoObject)
 template<class Type>
 void Foam::pointFieldReconstructor::reconstructFields
 (
-    const IOobjectList& objects
+    const IOobjectList& objects,
+    const HashSet<word>& selectedFields
 )
 {
     word fieldClassName
@@ -157,17 +158,19 @@ void Foam::pointFieldReconstructor::reconstructFields
     {
         Info<< "    Reconstructing " << fieldClassName << "s\n" << endl;
 
-        for
+        forAllConstIter(IOobjectList, fields, fieldIter)
+        {
+            if
             (
-            IOobjectList::iterator fieldIter = fields.begin();
+                !selectedFields.size()
-            fieldIter != fields.end();
+             || selectedFields.found(fieldIter()->name())
-            ++fieldIter
             )
             {
                 Info<< "        " << fieldIter()->name() << endl;
 
                 reconstructField<Type>(*fieldIter())().write();
             }
+        }
 
         Info<< endl;
     }

src/thermophysicalModels/basic/Make/files

@@ -6,11 +6,13 @@ basicThermo/basicThermo.C
 psiThermo/basicPsiThermo/basicPsiThermo.C
 psiThermo/basicPsiThermo/newBasicPsiThermo.C
 psiThermo/hPsiThermo/hPsiThermos.C
+psiThermo/hsPsiThermo/hsPsiThermos.C
 psiThermo/ePsiThermo/ePsiThermos.C
 
 rhoThermo/basicRhoThermo/basicRhoThermo.C
 rhoThermo/basicRhoThermo/newBasicRhoThermo.C
 rhoThermo/hRhoThermo/hRhoThermos.C
+rhoThermo/hsRhoThermo/hsRhoThermos.C
 
 derivedFvPatchFields/fixedEnthalpy/fixedEnthalpyFvPatchScalarField.C
 derivedFvPatchFields/gradientEnthalpy/gradientEnthalpyFvPatchScalarField.C

src/thermophysicalModels/basic/basicThermo/basicThermo.C

@@ -343,6 +343,13 @@ Foam::tmp<Foam::scalarField> Foam::basicThermo::hs
 }
 
 
+Foam::tmp<Foam::volScalarField> Foam::basicThermo::hc() const
+{
+    notImplemented("basicThermo::hc()");
+    return volScalarField::null();
+}
+
+
 Foam::volScalarField& Foam::basicThermo::e()
 {
     notImplemented("basicThermo::e()");

src/thermophysicalModels/basic/basicThermo/basicThermo.H

@@ -145,17 +145,17 @@ public:
             //  Non-const access allowed for transport equations
             virtual volScalarField& h();
 
-            //- Enthalpy [J/kg]
+            //- Total enthalpy [J/kg]
             virtual const volScalarField& h() const;
 
-            //- Enthalpy for cell-set [J/kg]
+            //- Total enthalpy for cell-set [J/kg]
             virtual tmp<scalarField> h
             (
                 const scalarField& T,
                 const labelList& cells
             ) const;
 
-            //- Enthalpy for patch [J/kg]
+            //- Total enthalpy for patch [J/kg]
             virtual tmp<scalarField> h
             (
                 const scalarField& T,
@@ -166,23 +166,26 @@ public:
             //  Non-const access allowed for transport equations
             virtual volScalarField& hs();
 
-            //- Enthalpy [J/kg]
+            //- Sensible enthalpy [J/kg]
             virtual const volScalarField& hs() const;
 
-            //- Enthalpy for cell-set [J/kg]
+            //- Sensible enthalpy for cell-set [J/kg]
             virtual tmp<scalarField> hs
             (
                 const scalarField& T,
                 const labelList& cells
             ) const;
 
-            //- Enthalpy for patch [J/kg]
+            //- Sensible enthalpy for patch [J/kg]
             virtual tmp<scalarField> hs
             (
                 const scalarField& T,
                 const label patchi
             ) const;
 
+            //- Chemical enthalpy [J/kg]
+            virtual tmp<volScalarField> hc() const;
+
             //- Internal energy [J/kg]
             //  Non-const access allowed for transport equations
             virtual volScalarField& e();

src/thermophysicalModels/basic/psiThermo/hsPsiThermo/hsPsiThermo.C

@@ -0,0 +1,347 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "hsPsiThermo.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class MixtureType>
+void Foam::hsPsiThermo<MixtureType>::calculate()
+{
+    const scalarField& hsCells = hs_.internalField();
+    const scalarField& pCells = this->p_.internalField();
+
+    scalarField& TCells = this->T_.internalField();
+    scalarField& psiCells = this->psi_.internalField();
+    scalarField& muCells = this->mu_.internalField();
+    scalarField& alphaCells = this->alpha_.internalField();
+
+    forAll(TCells, celli)
+    {
+        const typename MixtureType::thermoType& mixture_ =
+            this->cellMixture(celli);
+
+        TCells[celli] = mixture_.THs(hsCells[celli], TCells[celli]);
+        psiCells[celli] = mixture_.psi(pCells[celli], TCells[celli]);
+
+        muCells[celli] = mixture_.mu(TCells[celli]);
+        alphaCells[celli] = mixture_.alpha(TCells[celli]);
+    }
+
+    forAll(T_.boundaryField(), patchi)
+    {
+        fvPatchScalarField& pp = this->p_.boundaryField()[patchi];
+        fvPatchScalarField& pT = this->T_.boundaryField()[patchi];
+        fvPatchScalarField& ppsi = this->psi_.boundaryField()[patchi];
+
+        fvPatchScalarField& phs = hs_.boundaryField()[patchi];
+
+        fvPatchScalarField& pmu = this->mu_.boundaryField()[patchi];
+        fvPatchScalarField& palpha = this->alpha_.boundaryField()[patchi];
+
+        if (pT.fixesValue())
+        {
+            forAll(pT, facei)
+            {
+                const typename MixtureType::thermoType& mixture_ =
+                    this->patchFaceMixture(patchi, facei);
+
+                phs[facei] = mixture_.Hs(pT[facei]);
+
+                ppsi[facei] = mixture_.psi(pp[facei], pT[facei]);
+                pmu[facei] = mixture_.mu(pT[facei]);
+                palpha[facei] = mixture_.alpha(pT[facei]);
+            }
+        }
+        else
+        {
+            forAll(pT, facei)
+            {
+                const typename MixtureType::thermoType& mixture_ =
+                    this->patchFaceMixture(patchi, facei);
+
+                pT[facei] = mixture_.THs(phs[facei], pT[facei]);
+
+                ppsi[facei] = mixture_.psi(pp[facei], pT[facei]);
+                pmu[facei] = mixture_.mu(pT[facei]);
+                palpha[facei] = mixture_.alpha(pT[facei]);
+            }
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class MixtureType>
+Foam::hsPsiThermo<MixtureType>::hsPsiThermo(const fvMesh& mesh)
+:
+    basicPsiThermo(mesh),
+    MixtureType(*this, mesh),
+
+    hs_
+    (
+        IOobject
+        (
+            "hs",
+            mesh.time().timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimEnergy/dimMass,
+        this->hBoundaryTypes()
+    )
+{
+    scalarField& hsCells = hs_.internalField();
+    const scalarField& TCells = this->T_.internalField();
+
+    forAll(hsCells, celli)
+    {
+        hsCells[celli] = this->cellMixture(celli).Hs(TCells[celli]);
+    }
+
+    forAll(hs_.boundaryField(), patchi)
+    {
+        hs_.boundaryField()[patchi] ==
+            hs(this->T_.boundaryField()[patchi], patchi);
+    }
+
+    hBoundaryCorrection(hs_);
+
+    calculate();
+
+    // Switch on saving old time
+    this->psi_.oldTime();
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class MixtureType>
+Foam::hsPsiThermo<MixtureType>::~hsPsiThermo()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class MixtureType>
+void Foam::hsPsiThermo<MixtureType>::correct()
+{
+    if (debug)
+    {
+        Info<< "entering hsPsiThermo<MixtureType>::correct()" << endl;
+    }
+
+    // force the saving of the old-time values
+    this->psi_.oldTime();
+
+    calculate();
+
+    if (debug)
+    {
+        Info<< "exiting hsPsiThermo<MixtureType>::correct()" << endl;
+    }
+}
+
+
+template<class MixtureType>
+Foam::tmp<Foam::scalarField> Foam::hsPsiThermo<MixtureType>::hs
+(
+    const scalarField& T,
+    const labelList& cells
+) const
+{
+    tmp<scalarField> ths(new scalarField(T.size()));
+    scalarField& hs = ths();
+
+    forAll(T, celli)
+    {
+        hs[celli] = this->cellMixture(cells[celli]).Hs(T[celli]);
+    }
+
+    return ths;
+}
+
+
+template<class MixtureType>
+Foam::tmp<Foam::scalarField> Foam::hsPsiThermo<MixtureType>::hs
+(
+    const scalarField& T,
+    const label patchi
+) const
+{
+    tmp<scalarField> ths(new scalarField(T.size()));
+    scalarField& hs = ths();
+
+    forAll(T, facei)
+    {
+        hs[facei] = this->patchFaceMixture(patchi, facei).Hs(T[facei]);
+    }
+
+    return ths;
+}
+
+
+template<class MixtureType>
+Foam::tmp<Foam::scalarField> Foam::hsPsiThermo<MixtureType>::Cp
+(
+    const scalarField& T,
+    const label patchi
+) const
+{
+    tmp<scalarField> tCp(new scalarField(T.size()));
+    scalarField& cp = tCp();
+
+    forAll(T, facei)
+    {
+        cp[facei] = this->patchFaceMixture(patchi, facei).Cp(T[facei]);
+    }
+
+    return tCp;
+}
+
+
+template<class MixtureType>
+Foam::tmp<Foam::volScalarField> Foam::hsPsiThermo<MixtureType>::Cp() const
+{
+    const fvMesh& mesh = this->T_.mesh();
+
+    tmp<volScalarField> tCp
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "Cp",
+                mesh.time().timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh,
+            dimensionSet(0, 2, -2, -1, 0),
+            this->T_.boundaryField().types()
+        )
+    );
+
+    volScalarField& cp = tCp();
+
+    forAll(this->T_, celli)
+    {
+        cp[celli] = this->cellMixture(celli).Cp(this->T_[celli]);
+    }
+
+    forAll(this->T_.boundaryField(), patchi)
+    {
+        const fvPatchScalarField& pT = this->T_.boundaryField()[patchi];
+        fvPatchScalarField& pCp = cp.boundaryField()[patchi];
+
+        forAll(pT, facei)
+        {
+            pCp[facei] = this->patchFaceMixture(patchi, facei).Cp(pT[facei]);
+        }
+    }
+
+    return tCp;
+}
+
+
+template<class MixtureType>
+Foam::tmp<Foam::scalarField> Foam::hsPsiThermo<MixtureType>::Cv
+(
+    const scalarField& T,
+    const label patchi
+) const
+{
+    tmp<scalarField> tCv(new scalarField(T.size()));
+    scalarField& cv = tCv();
+
+    forAll(T, facei)
+    {
+        cv[facei] = this->patchFaceMixture(patchi, facei).Cv(T[facei]);
+    }
+
+    return tCv;
+}
+
+
+template<class MixtureType>
+Foam::tmp<Foam::volScalarField> Foam::hsPsiThermo<MixtureType>::Cv() const
+{
+    const fvMesh& mesh = this->T_.mesh();
+
+    tmp<volScalarField> tCv
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "Cv",
+                mesh.time().timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh,
+            dimEnergy/dimMass/dimTemperature
+        )
+    );
+
+    volScalarField& cv = tCv();
+
+    forAll(this->T_, celli)
+    {
+        cv[celli] = this->cellMixture(celli).Cv(this->T_[celli]);
+    }
+
+    forAll(this->T_.boundaryField(), patchi)
+    {
+        cv.boundaryField()[patchi] =
+            Cv(this->T_.boundaryField()[patchi], patchi);
+    }
+
+    return tCv;
+}
+
+
+template<class MixtureType>
+bool Foam::hsPsiThermo<MixtureType>::read()
+{
+    if (basicPsiThermo::read())
+    {
+        MixtureType::read(*this);
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+// ************************************************************************* //

src/thermophysicalModels/basic/psiThermo/hsPsiThermo/hsPsiThermo.H

@@ -0,0 +1,178 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::hsPsiThermo
+
+Description
+    Sensible enthalpy for a mixture based on compressibility
+
+SourceFiles
+    hsPsiThermo.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef hsPsiThermo_H
+#define hsPsiThermo_H
+
+#include "basicPsiThermo.H"
+#include "basicMixture.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class hsPsiThermo Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class MixtureType>
+class hsPsiThermo
+:
+    public basicPsiThermo,
+    public MixtureType
+{
+    // Private data
+
+        //- Sensible enthalpy field [J/kg]
+        volScalarField hs_;
+
+
+    // Private member functions
+
+        //- Calculate the thermo variables
+        void calculate();
+
+        //- Construct as copy (not implemented)
+        hsPsiThermo(const hsPsiThermo<MixtureType>&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("hsPsiThermo");
+
+
+    // Constructors
+
+        //- Construct from mesh
+        hsPsiThermo(const fvMesh&);
+
+
+    //- Destructor
+    virtual ~hsPsiThermo();
+
+
+    // Member functions
+
+        //- Return the compostion of the mixture
+        virtual basicMixture& composition()
+        {
+            return *this;
+        }
+
+        //- Return the compostion of the mixture
+        virtual const basicMixture& composition() const
+        {
+            return *this;
+        }
+
+        //- Update properties
+        virtual void correct();
+
+
+        // Access to thermodynamic state variables
+
+            //- Sensible enthalpy [J/kg]
+            //  Non-const access allowed for transport equations
+            virtual volScalarField& hs()
+            {
+                return hs_;
+            }
+
+            //- Sensible enthalpy [J/kg]
+            virtual const volScalarField& hs() const
+            {
+                return hs_;
+            }
+
+
+        // Fields derived from thermodynamic state variables
+
+            //- Enthalpy for cell-set [J/kg]
+            virtual tmp<scalarField> hs
+            (
+                const scalarField& T,
+                const labelList& cells
+            ) const;
+
+            //- Enthalpy for patch [J/kg]
+            virtual tmp<scalarField> hs
+            (
+                const scalarField& T,
+                const label patchi
+            ) const;
+
+            //- Heat capacity at constant pressure for patch [J/kg/K]
+            virtual tmp<scalarField> Cp
+            (
+                const scalarField& T,
+                const label patchi
+            ) const;
+
+            //- Heat capacity at constant pressure [J/kg/K]
+            virtual tmp<volScalarField> Cp() const;
+
+            //- Heat capacity at constant volume for patch [J/kg/K]
+            virtual tmp<scalarField> Cv
+            (
+                const scalarField& T,
+                const label patchi
+            ) const;
+
+            //- Heat capacity at constant volume [J/kg/K]
+            virtual tmp<volScalarField> Cv() const;
+
+
+        //- Read thermophysicalProperties dictionary
+        virtual bool read();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
+
+#ifdef NoRepository
+#   include "hsPsiThermo.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/thermophysicalModels/basic/psiThermo/hsPsiThermo/hsPsiThermos.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -24,10 +24,7 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "makeHsCombustionThermo.H"
+#include "makeBasicPsiThermo.H"
-
-#include "hsCombustionThermo.H"
-#include "hsPsiMixtureThermo.H"
 
 #include "perfectGas.H"
 
@@ -35,32 +32,47 @@ License
 #include "janafThermo.H"
 #include "specieThermo.H"
 
+#include "constTransport.H"
 #include "sutherlandTransport.H"
 
-#include "veryInhomogeneousMixture.H"
+#include "hsPsiThermo.H"
-
+#include "pureMixture.H"
-#include "multiComponentMixture.H"
-
-#include "thermoPhysicsTypes.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+/* * * * * * * * * * * * * * * private static data * * * * * * * * * * * * * */
 
-
+makeBasicPsiThermo
-makeHsCombustionThermo
 (
-    hsCombustionThermo,
+    hsPsiThermo,
-    hsPsiMixtureThermo,
+    pureMixture,
-    veryInhomogeneousMixture,
+    constTransport,
+    hConstThermo,
+    perfectGas
+);
+
+makeBasicPsiThermo
+(
+    hsPsiThermo,
+    pureMixture,
+    sutherlandTransport,
+    hConstThermo,
+    perfectGas
+);
+
+makeBasicPsiThermo
+(
+    hsPsiThermo,
+    pureMixture,
     sutherlandTransport,
     janafThermo,
     perfectGas
 );
 
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/thermophysicalModels/basic/rhoThermo/hsRhoThermo/hsRhoThermo.C

@@ -0,0 +1,346 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "hsRhoThermo.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class MixtureType>
+void Foam::hsRhoThermo<MixtureType>::calculate()
+{
+    const scalarField& hsCells = this->hs_.internalField();
+    const scalarField& pCells = this->p_.internalField();
+
+    scalarField& TCells = this->T_.internalField();
+    scalarField& psiCells = this->psi_.internalField();
+    scalarField& rhoCells = this->rho_.internalField();
+    scalarField& muCells = this->mu_.internalField();
+    scalarField& alphaCells = this->alpha_.internalField();
+
+    forAll(TCells, celli)
+    {
+        const typename MixtureType::thermoType& mixture_ =
+            this->cellMixture(celli);
+
+        TCells[celli] = mixture_.THs(hsCells[celli], TCells[celli]);
+        psiCells[celli] = mixture_.psi(pCells[celli], TCells[celli]);
+        rhoCells[celli] = mixture_.rho(pCells[celli], TCells[celli]);
+
+        muCells[celli] = mixture_.mu(TCells[celli]);
+        alphaCells[celli] = mixture_.alpha(TCells[celli]);
+    }
+
+    forAll(this->T_.boundaryField(), patchi)
+    {
+        fvPatchScalarField& pp = this->p_.boundaryField()[patchi];
+        fvPatchScalarField& pT = this->T_.boundaryField()[patchi];
+        fvPatchScalarField& ppsi = this->psi_.boundaryField()[patchi];
+        fvPatchScalarField& prho = this->rho_.boundaryField()[patchi];
+
+        fvPatchScalarField& phs = this->hs_.boundaryField()[patchi];
+
+        fvPatchScalarField& pmu = this->mu_.boundaryField()[patchi];
+        fvPatchScalarField& palpha = this->alpha_.boundaryField()[patchi];
+
+        if (pT.fixesValue())
+        {
+            forAll(pT, facei)
+            {
+                const typename MixtureType::thermoType& mixture_ =
+                    this->patchFaceMixture(patchi, facei);
+
+                phs[facei] = mixture_.Hs(pT[facei]);
+
+                ppsi[facei] = mixture_.psi(pp[facei], pT[facei]);
+                prho[facei] = mixture_.rho(pp[facei], pT[facei]);
+                pmu[facei] = mixture_.mu(pT[facei]);
+                palpha[facei] = mixture_.alpha(pT[facei]);
+            }
+        }
+        else
+        {
+            forAll(pT, facei)
+            {
+                const typename MixtureType::thermoType& mixture_ =
+                    this->patchFaceMixture(patchi, facei);
+
+                pT[facei] = mixture_.THs(phs[facei], pT[facei]);
+
+                ppsi[facei] = mixture_.psi(pp[facei], pT[facei]);
+                prho[facei] = mixture_.rho(pp[facei], pT[facei]);
+                pmu[facei] = mixture_.mu(pT[facei]);
+                palpha[facei] = mixture_.alpha(pT[facei]);
+            }
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class MixtureType>
+Foam::hsRhoThermo<MixtureType>::hsRhoThermo(const fvMesh& mesh)
+:
+    basicRhoThermo(mesh),
+    MixtureType(*this, mesh),
+
+    hs_
+    (
+        IOobject
+        (
+            "hs",
+            mesh.time().timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimEnergy/dimMass,
+        this->hBoundaryTypes()
+    )
+{
+    scalarField& hsCells = hs_.internalField();
+    const scalarField& TCells = this->T_.internalField();
+
+    forAll(hsCells, celli)
+    {
+        hsCells[celli] = this->cellMixture(celli).Hs(TCells[celli]);
+    }
+
+    forAll(hs_.boundaryField(), patchi)
+    {
+        hs_.boundaryField()[patchi] ==
+            hs(this->T_.boundaryField()[patchi], patchi);
+    }
+
+    hBoundaryCorrection(hs_);
+
+    calculate();
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class MixtureType>
+Foam::hsRhoThermo<MixtureType>::~hsRhoThermo()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class MixtureType>
+void Foam::hsRhoThermo<MixtureType>::correct()
+{
+    if (debug)
+    {
+        Info<< "entering hsRhoThermo<MixtureType>::correct()" << endl;
+    }
+
+    calculate();
+
+    if (debug)
+    {
+        Info<< "exiting hsRhoThermo<MixtureType>::correct()" << endl;
+    }
+}
+
+
+template<class MixtureType>
+Foam::tmp<Foam::scalarField> Foam::hsRhoThermo<MixtureType>::hs
+(
+    const scalarField& T,
+    const labelList& cells
+) const
+{
+    tmp<scalarField> ths(new scalarField(T.size()));
+    scalarField& hs = ths();
+
+    forAll(T, celli)
+    {
+        hs[celli] = this->cellMixture(cells[celli]).Hs(T[celli]);
+    }
+
+    return ths;
+}
+
+
+template<class MixtureType>
+Foam::tmp<Foam::scalarField> Foam::hsRhoThermo<MixtureType>::hs
+(
+    const scalarField& T,
+    const label patchi
+) const
+{
+    tmp<scalarField> ths(new scalarField(T.size()));
+    scalarField& hs = ths();
+
+    forAll(T, facei)
+    {
+        hs[facei] = this->patchFaceMixture(patchi, facei).Hs(T[facei]);
+    }
+
+    return ths;
+}
+
+
+template<class MixtureType>
+Foam::tmp<Foam::scalarField> Foam::hsRhoThermo<MixtureType>::Cp
+(
+    const scalarField& T,
+    const label patchi
+) const
+{
+    tmp<scalarField> tCp(new scalarField(T.size()));
+    scalarField& cp = tCp();
+
+    forAll(T, facei)
+    {
+        cp[facei] = this->patchFaceMixture(patchi, facei).Cp(T[facei]);
+    }
+
+    return tCp;
+}
+
+
+template<class MixtureType>
+Foam::tmp<Foam::volScalarField> Foam::hsRhoThermo<MixtureType>::Cp() const
+{
+    const fvMesh& mesh = this->T_.mesh();
+
+    tmp<volScalarField> tCp
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "Cp",
+                mesh.time().timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh,
+            dimEnergy/dimMass/dimTemperature,
+            this->T_.boundaryField().types()
+        )
+    );
+
+    volScalarField& cp = tCp();
+
+    forAll(this->T_, celli)
+    {
+        cp[celli] = this->cellMixture(celli).Cp(this->T_[celli]);
+    }
+
+    forAll(this->T_.boundaryField(), patchi)
+    {
+        const fvPatchScalarField& pT = this->T_.boundaryField()[patchi];
+        fvPatchScalarField& pCp = cp.boundaryField()[patchi];
+
+        forAll(pT, facei)
+        {
+            pCp[facei] = this->patchFaceMixture(patchi, facei).Cp(pT[facei]);
+        }
+    }
+
+    return tCp;
+}
+
+
+template<class MixtureType>
+Foam::tmp<Foam::scalarField> Foam::hsRhoThermo<MixtureType>::Cv
+(
+    const scalarField& T,
+    const label patchi
+) const
+{
+    tmp<scalarField> tCv(new scalarField(T.size()));
+    scalarField& cv = tCv();
+
+    forAll(T, facei)
+    {
+        cv[facei] = this->patchFaceMixture(patchi, facei).Cv(T[facei]);
+    }
+
+    return tCv;
+}
+
+
+template<class MixtureType>
+Foam::tmp<Foam::volScalarField> Foam::hsRhoThermo<MixtureType>::Cv() const
+{
+    const fvMesh& mesh = this->T_.mesh();
+
+    tmp<volScalarField> tCv
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "Cv",
+                mesh.time().timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh,
+            dimEnergy/dimMass/dimTemperature
+        )
+    );
+
+    volScalarField& cv = tCv();
+
+    forAll(this->T_, celli)
+    {
+        cv[celli] = this->cellMixture(celli).Cv(this->T_[celli]);
+    }
+
+    forAll(this->T_.boundaryField(), patchi)
+    {
+        cv.boundaryField()[patchi] =
+            Cv(this->T_.boundaryField()[patchi], patchi);
+    }
+
+    return tCv;
+}
+
+
+template<class MixtureType>
+bool Foam::hsRhoThermo<MixtureType>::read()
+{
+    if (basicRhoThermo::read())
+    {
+        MixtureType::read(*this);
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+// ************************************************************************* //

src/thermophysicalModels/basic/rhoThermo/hsRhoThermo/hsRhoThermo.H

@@ -0,0 +1,178 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::hsRhoThermo
+
+Description
+    Sensible enthalpy for a mixture based on density
+
+SourceFiles
+    hsRhoThermo.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef hsRhoThermo_H
+#define hsRhoThermo_H
+
+#include "basicRhoThermo.H"
+#include "basicMixture.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class hsRhoThermo Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class MixtureType>
+class hsRhoThermo
+:
+    public basicRhoThermo,
+    public MixtureType
+{
+    // Private data
+
+        //- Sensible enthalpy field [J/kg]
+        volScalarField hs_;
+
+
+    // Private member functions
+
+        //- Calculate the thermo variables
+        void calculate();
+
+        //- Construct as copy (not implemented)
+        hsRhoThermo(const hsRhoThermo<MixtureType>&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("hsRhoThermo");
+
+
+    // Constructors
+
+        //- Construct from mesh
+        hsRhoThermo(const fvMesh&);
+
+
+    //- Destructor
+    virtual ~hsRhoThermo();
+
+
+    // Member functions
+
+        //- Return the compostion of the combustion mixture
+        virtual basicMixture& composition()
+        {
+            return *this;
+        }
+
+        //- Return the compostion of the combustion mixture
+        virtual const basicMixture& composition() const
+        {
+            return *this;
+        }
+
+        //- Update properties
+        virtual void correct();
+
+
+        // Access to thermodynamic state variables
+
+            //- Sensible enthalpy [J/kg]
+            //  Non-const access allowed for transport equations
+            virtual volScalarField& hs()
+            {
+                return hs_;
+            }
+
+            //- Sensible enthalpy [J/kg]
+            virtual const volScalarField& hs() const
+            {
+                return hs_;
+            }
+
+
+        // Fields derived from thermodynamic state variables
+
+            //- Sensible enthalpy for cell-set [J/kg]
+            virtual tmp<scalarField> hs
+            (
+                const scalarField& T,
+                const labelList& cells
+            ) const;
+
+            //- Sensible enthalpy for patch [J/kg]
+            virtual tmp<scalarField> hs
+            (
+                const scalarField& T,
+                const label patchi
+            ) const;
+
+            //- Heat capacity at constant pressure for patch [J/kg/K]
+            virtual tmp<scalarField> Cp
+            (
+                const scalarField& T,
+                const label patchi
+            ) const;
+
+            //- Heat capacity at constant pressure [J/kg/K]
+            virtual tmp<volScalarField> Cp() const;
+
+            //- Heat capacity at constant volume for patch [J/kg/K]
+            virtual tmp<scalarField> Cv
+            (
+                const scalarField& T,
+                const label patchi
+            ) const;
+
+            //- Heat capacity at constant volume [J/kg/K]
+            virtual tmp<volScalarField> Cv() const;
+
+
+        //- Read thermophysicalProperties dictionary
+        virtual bool read();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
+
+#ifdef NoRepository
+#   include "hsRhoThermo.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/thermophysicalModels/basic/rhoThermo/hsRhoThermo/hsRhoThermos.C

@@ -0,0 +1,80 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "makeBasicRhoThermo.H"
+
+#include "perfectGas.H"
+
+#include "hConstThermo.H"
+#include "janafThermo.H"
+#include "specieThermo.H"
+
+#include "constTransport.H"
+#include "sutherlandTransport.H"
+
+#include "hsRhoThermo.H"
+#include "pureMixture.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/* * * * * * * * * * * * * * * private static data * * * * * * * * * * * * * */
+
+makeBasicRhoThermo
+(
+    hsRhoThermo,
+    pureMixture,
+    constTransport,
+    hConstThermo,
+    perfectGas
+);
+
+makeBasicRhoThermo
+(
+    hsRhoThermo,
+    pureMixture,
+    sutherlandTransport,
+    hConstThermo,
+    perfectGas
+);
+
+makeBasicRhoThermo
+(
+    hsRhoThermo,
+    pureMixture,
+    sutherlandTransport,
+    janafThermo,
+    perfectGas
+);
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/thermophysicalModels/chemistryModel/chemistryModel/ODEChemistryModel/ODEChemistryModel.C

@@ -466,7 +466,7 @@ Foam::ODEChemistryModel<CompType, ThermoType>::tc() const
         this->thermo().rho()
     );
 
-    tmp<volScalarField> tsource
+    tmp<volScalarField> ttc
     (
         new volScalarField
         (
@@ -484,7 +484,7 @@ Foam::ODEChemistryModel<CompType, ThermoType>::tc() const
         )
     );
 
-    scalarField& t = tsource();
+    scalarField& tc = ttc();
 
     label nReaction = reactions_.size();
 
@@ -517,17 +517,58 @@ Foam::ODEChemistryModel<CompType, ThermoType>::tc() const
                 forAll(R.rhs(), s)
                 {
                     scalar sr = R.rhs()[s].stoichCoeff;
-                    t[celli] += sr*pf*cf;
+                    tc[celli] += sr*pf*cf;
                 }
             }
-            t[celli] = nReaction*cSum/t[celli];
+            tc[celli] = nReaction*cSum/tc[celli];
         }
     }
 
 
-    tsource().correctBoundaryConditions();
+    ttc().correctBoundaryConditions();
 
-    return tsource;
+    return ttc;
+}
+
+
+template<class CompType, class ThermoType>
+Foam::tmp<Foam::volScalarField>
+Foam::ODEChemistryModel<CompType, ThermoType>::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->chemistry_)
+    {
+        scalarField& Sh = tSh();
+
+        forAll(Y_, i)
+        {
+            forAll(Sh, cellI)
+            {
+                scalar hi = specieThermo_[i].Hc();
+                Sh[cellI] -= hi*RR_[i][cellI];
+            }
+        }
+    }
+
+    return tSh;
 }
 
 
@@ -545,37 +586,19 @@ Foam::ODEChemistryModel<CompType, ThermoType>::dQ() const
                 this->mesh_.time().timeName(),
                 this->mesh_,
                 IOobject::NO_READ,
-                IOobject::NO_WRITE
+                IOobject::NO_WRITE,
+                false
             ),
             this->mesh_,
-            dimensionedScalar
+            dimensionedScalar("dQ", dimEnergy/dimTime, 0.0),
-            (
+            zeroGradientFvPatchScalarField::typeName
-                "zero",
-                dimensionSet(0, 2, -3 , 0, 0, 0, 0),
-                0.0
-            )
         )
     );
 
     if (this->chemistry_)
     {
-        scalarField& dQ = tdQ();
+        volScalarField& dQ = tdQ();
-
+        dQ.dimensionedInternalField() = this->mesh_.V()*Sh()();
-        scalarField rhoEff(dQ.size(), 0.0);
-
-        forAll(Y_, i)
-        {
-            forAll(dQ, cellI)
-            {
-                scalar Ti = this->thermo().T()[cellI];
-                scalar pi = this->thermo().p()[cellI];
-                rhoEff[cellI] += Y_[i][cellI]*specieThermo_[i].rho(pi, Ti);
-                scalar hi = specieThermo_[i].H(Ti);
-                dQ[cellI] -= hi*RR_[i][cellI];
-            }
-        }
-
-        dQ /= rhoEff;
     }
 
     return tdQ;
@@ -678,6 +701,9 @@ Foam::scalar Foam::ODEChemistryModel<CompType, ThermoType>::solve
 
     scalar deltaTMin = GREAT;
 
+    tmp<volScalarField> thc = this->thermo().hc();
+    const scalarField& hc = thc();
+
     forAll(rho, celli)
     {
         for (label i=0; i<nSpecie_; i++)
@@ -687,7 +713,7 @@ Foam::scalar Foam::ODEChemistryModel<CompType, ThermoType>::solve
 
         scalar rhoi = rho[celli];
         scalar Ti = this->thermo().T()[celli];
-        scalar hi = this->thermo().h()[celli];
+        scalar hi = this->thermo().hs()[celli] + hc[celli];
         scalar pi = this->thermo().p()[celli];
 
         scalarField c(nSpecie_);

src/thermophysicalModels/chemistryModel/chemistryModel/ODEChemistryModel/ODEChemistryModel.H

@@ -39,9 +39,9 @@ SourceFiles
 #ifndef ODEChemistryModel_H
 #define ODEChemistryModel_H
 
-#include "hCombustionThermo.H"
 #include "Reaction.H"
 #include "ODE.H"
+#include "volFieldsFwd.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -49,6 +49,8 @@ namespace Foam
 {
 
 // Forward declaration of classes
+class fvMesh;
+
 template<class CompType, class ThermoType>
 class chemistrySolver;
 
@@ -180,6 +182,9 @@ public:
             //- Return the chemical time scale
             virtual tmp<volScalarField> tc() const;
 
+            //- Return source for enthalpy equation [kg/m/s3]
+            virtual tmp<volScalarField> Sh() const;
+
             //- Return the heat release, i.e. enthalpy/sec [m2/s3]
             virtual tmp<volScalarField> dQ() const;
 

src/thermophysicalModels/chemistryModel/chemistryModel/ODEChemistryModel/ODEChemistryModelI.H

@@ -24,6 +24,7 @@ License
 
 \*---------------------------------------------------------------------------*/
 
+#include "volFields.H"
 #include "zeroGradientFvPatchFields.H"
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
@@ -96,7 +97,7 @@ Foam::ODEChemistryModel<CompType, ThermoType>::RR
                 IOobject::NO_WRITE
             ),
             this->mesh(),
-            dimensionedScalar("zero", dimensionSet(1, -3, -1, 0, 0), 0.0),
+            dimensionedScalar("zero", dimMass/dimVolume/dimTime, 0.0),
             zeroGradientFvPatchScalarField::typeName
         )
     );

src/thermophysicalModels/chemistryModel/chemistryModel/basicChemistryModel/basicChemistryModel.H

@@ -120,7 +120,7 @@ public:
 
             // Fields
 
-                //- Return const access to chemical source terms
+                //- Return const access to chemical source terms [kg/m3/s]
                 virtual tmp<volScalarField> RR(const label i) const = 0;
 
 
@@ -133,7 +133,10 @@ public:
                 //- Return the chemical time scale
                 virtual tmp<volScalarField> tc() const = 0;
 
-                //- Return the heat release
+                //- Return source for enthalpy equation [kg/m/s3]
+                virtual tmp<volScalarField> Sh() const = 0;
+
+                //- Return the heat release, i.e. enthalpy/sec [m2/s3]
                 virtual tmp<volScalarField> dQ() const = 0;
 };
 

src/thermophysicalModels/chemistryModel/chemistryModel/psiChemistryModel/psiChemistryModel.C

@@ -45,7 +45,7 @@ Foam::psiChemistryModel::psiChemistryModel
 )
 :
     basicChemistryModel(mesh),
-    thermo_(hCombustionThermo::NewType(mesh, thermoTypeName))
+    thermo_(hsCombustionThermo::NewType(mesh, thermoTypeName))
 {}
 
 

src/thermophysicalModels/chemistryModel/chemistryModel/psiChemistryModel/psiChemistryModel.H

@@ -41,7 +41,7 @@ SourceFiles
 #include "basicChemistryModel.H"
 #include "autoPtr.H"
 #include "runTimeSelectionTables.H"
-#include "hCombustionThermo.H"
+#include "hsCombustionThermo.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -73,7 +73,7 @@ protected:
     // Protected data
 
         //- Thermo package
-        autoPtr<hCombustionThermo> thermo_;
+        autoPtr<hsCombustionThermo> thermo_;
 
 
 public:
@@ -114,10 +114,10 @@ public:
     // Member Functions
 
         //- Return access to the thermo package
-        inline hCombustionThermo& thermo();
+        inline hsCombustionThermo& thermo();
 
         //- Return const access to the thermo package
-        inline const hCombustionThermo& thermo() const;
+        inline const hsCombustionThermo& thermo() const;
 };
 
 

src/thermophysicalModels/chemistryModel/chemistryModel/psiChemistryModel/psiChemistryModelI.H

@@ -26,13 +26,13 @@ License
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-inline Foam::hCombustionThermo& Foam::psiChemistryModel::thermo()
+inline Foam::hsCombustionThermo& Foam::psiChemistryModel::thermo()
 {
     return thermo_();
 }
 
 
-inline const Foam::hCombustionThermo& Foam::psiChemistryModel::thermo() const
+inline const Foam::hsCombustionThermo& Foam::psiChemistryModel::thermo() const
 {
     return thermo_();
 }

src/thermophysicalModels/chemistryModel/chemistryModel/rhoChemistryModel/rhoChemistryModel.C

@@ -45,7 +45,7 @@ Foam::rhoChemistryModel::rhoChemistryModel
 )
 :
     basicChemistryModel(mesh),
-    thermo_(hReactionThermo::NewType(mesh, thermoTypeName))
+    thermo_(hsReactionThermo::NewType(mesh, thermoTypeName))
 {}
 
 

src/thermophysicalModels/chemistryModel/chemistryModel/rhoChemistryModel/rhoChemistryModel.H

@@ -41,7 +41,7 @@ SourceFiles
 #include "basicChemistryModel.H"
 #include "autoPtr.H"
 #include "runTimeSelectionTables.H"
-#include "hReactionThermo.H"
+#include "hsReactionThermo.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -73,7 +73,7 @@ protected:
     // Protected data
 
         //- Thermo package
-        autoPtr<hReactionThermo> thermo_;
+        autoPtr<hsReactionThermo> thermo_;
 
 
 public:
@@ -114,10 +114,10 @@ public:
     // Member Functions
 
         //- Return access to the thermo package
-        inline hReactionThermo& thermo();
+        inline hsReactionThermo& thermo();
 
         //- Return const access to the thermo package
-        inline const hReactionThermo& thermo() const;
+        inline const hsReactionThermo& thermo() const;
 };
 
 

src/thermophysicalModels/chemistryModel/chemistryModel/rhoChemistryModel/rhoChemistryModelI.H

@@ -26,13 +26,13 @@ License
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-inline Foam::hReactionThermo& Foam::rhoChemistryModel::thermo()
+inline Foam::hsReactionThermo& Foam::rhoChemistryModel::thermo()
 {
     return thermo_();
 }
 
 
-inline const Foam::hReactionThermo& Foam::rhoChemistryModel::thermo() const
+inline const Foam::hsReactionThermo& Foam::rhoChemistryModel::thermo() const
 {
     return thermo_();
 }

src/thermophysicalModels/reactionThermo/Make/files

@@ -8,6 +8,10 @@ combustionThermo/hCombustionThermo/hCombustionThermo.C
 combustionThermo/hCombustionThermo/newhCombustionThermo.C
 combustionThermo/hCombustionThermo/hCombustionThermos.C
 
+combustionThermo/hsCombustionThermo/hsCombustionThermo.C
+combustionThermo/hsCombustionThermo/newhsCombustionThermo.C
+combustionThermo/hsCombustionThermo/hsCombustionThermos.C
+
 combustionThermo/hhuCombustionThermo/hhuCombustionThermo.C
 combustionThermo/hhuCombustionThermo/newhhuCombustionThermo.C
 combustionThermo/hhuCombustionThermo/hhuCombustionThermos.C
@@ -16,6 +20,9 @@ reactionThermo/hReactionThermo/hReactionThermo.C
 reactionThermo/hReactionThermo/newhReactionThermo.C
 reactionThermo/hReactionThermo/hReactionThermos.C
 
+reactionThermo/hsReactionThermo/hsReactionThermo.C
+reactionThermo/hsReactionThermo/newhsReactionThermo.C
+reactionThermo/hsReactionThermo/hsReactionThermos.C
 
 derivedFvPatchFields/fixedUnburntEnthalpy/fixedUnburntEnthalpyFvPatchScalarField.C
 derivedFvPatchFields/gradientUnburntEnthalpy/gradientUnburntEnthalpyFvPatchScalarField.C

src/thermophysicalModels/reactionThermo/combustionThermo/hCombustionThermo/hCombustionThermo.H

@@ -128,9 +128,6 @@ public:
             }
 
 
-        //- Chemical enthalpy [J/kg]
-        virtual tmp<volScalarField> hc() const = 0;
-
         //- Update properties
         virtual void correct() = 0;
 };

src/thermophysicalModels/reactionThermo/combustionThermo/hCombustionThermo/makeCombustionThermo.H

@@ -25,8 +25,6 @@ License
 InClass
     Foam::hCombustionThermo
 
-Description
-
 \*---------------------------------------------------------------------------*/
 
 #ifndef makeCombustionThermo_H

src/thermophysicalModels/reactionThermo/combustionThermo/hsCombustionThermo/hsCombustionThermo.C

@@ -39,7 +39,7 @@ namespace Foam
 
 Foam::hsCombustionThermo::hsCombustionThermo(const fvMesh& mesh)
 :
-    basicSensiblePsiThermo(mesh),
+    basicPsiThermo(mesh),
 
     hs_
     (
@@ -52,7 +52,7 @@ Foam::hsCombustionThermo::hsCombustionThermo(const fvMesh& mesh)
             IOobject::NO_WRITE
         ),
         mesh,
-        dimensionSet(0, 2, -2, 0, 0),
+        dimEnergy/dimMass,
         this->hBoundaryTypes()
     )
 {}

src/thermophysicalModels/reactionThermo/combustionThermo/hsCombustionThermo/hsCombustionThermo.H

@@ -26,7 +26,7 @@ Class
     Foam::hsCombustionThermo
 
 Description
-    Foam::hsCombustionThermo
+    Sensible enthalpy variant of combustionThermo
 
 SourceFiles
     hsCombustionThermo.C
@@ -36,7 +36,7 @@ SourceFiles
 #ifndef hsCombustionThermo_H
 #define hsCombustionThermo_H
 
-#include "basicSensiblePsiThermo.H"
+#include "basicPsiThermo.H"
 #include "basicMultiComponentMixture.H"
 #include "autoPtr.H"
 #include "runTimeSelectionTables.H"
@@ -52,14 +52,14 @@ namespace Foam
 
 class hsCombustionThermo
 :
-    public basicSensiblePsiThermo
+    public basicPsiThermo
 {
 
 protected:
 
     // Protected data
 
-        //- sensible enthalpy field
+        //- Sensible enthalpy field
         volScalarField hs_;
 
 
@@ -127,8 +127,6 @@ public:
                 return hs_;
             }
 
-        //- Chemical enthalpy [J/kg]
-        virtual tmp<volScalarField> hc() const = 0;
 
         //- Update properties
         virtual void correct() = 0;

src/thermophysicalModels/reactionThermo/combustionThermo/hsCombustionThermo/hsCombustionThermos.C

@@ -0,0 +1,153 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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 "makeHsCombustionThermo.H"
+
+#include "hsCombustionThermo.H"
+#include "hsPsiMixtureThermo.H"
+
+#include "perfectGas.H"
+
+#include "hConstThermo.H"
+#include "janafThermo.H"
+#include "specieThermo.H"
+
+#include "constTransport.H"
+#include "sutherlandTransport.H"
+
+#include "dieselMixture.H"
+#include "homogeneousMixture.H"
+#include "inhomogeneousMixture.H"
+#include "veryInhomogeneousMixture.H"
+
+#include "reactingMixture.H"
+#include "multiComponentMixture.H"
+
+#include "thermoPhysicsTypes.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+makeHsCombustionThermo
+(
+    hsCombustionThermo,
+    hsPsiMixtureThermo,
+    homogeneousMixture,
+    constTransport,
+    hConstThermo,
+    perfectGas
+);
+
+makeHsCombustionThermo
+(
+    hsCombustionThermo,
+    hsPsiMixtureThermo,
+    inhomogeneousMixture,
+    constTransport,
+    hConstThermo,
+    perfectGas
+);
+
+makeHsCombustionThermo
+(
+    hsCombustionThermo,
+    hsPsiMixtureThermo,
+    veryInhomogeneousMixture,
+    constTransport,
+    hConstThermo,
+    perfectGas
+);
+
+makeHsCombustionThermo
+(
+    hsCombustionThermo,
+    hsPsiMixtureThermo,
+    homogeneousMixture,
+    sutherlandTransport,
+    janafThermo,
+    perfectGas
+);
+
+makeHsCombustionThermo
+(
+    hsCombustionThermo,
+    hsPsiMixtureThermo,
+    inhomogeneousMixture,
+    sutherlandTransport,
+    janafThermo,
+    perfectGas
+);
+
+makeHsCombustionThermo
+(
+    hsCombustionThermo,
+    hsPsiMixtureThermo,
+    veryInhomogeneousMixture,
+    sutherlandTransport,
+    janafThermo,
+    perfectGas
+);
+
+makeHsCombustionThermo
+(
+    hsCombustionThermo,
+    hsPsiMixtureThermo,
+    dieselMixture,
+    sutherlandTransport,
+    janafThermo,
+    perfectGas
+);
+
+// Multi-component thermo
+
+makeHsCombustionMixtureThermo
+(
+    hsCombustionThermo,
+    hsPsiMixtureThermo,
+    multiComponentMixture,
+    gasThermoPhysics
+);
+
+
+// Multi-component reaction thermo
+
+makeHsCombustionMixtureThermo
+(
+    hsCombustionThermo,
+    hsPsiMixtureThermo,
+    reactingMixture,
+    gasThermoPhysics
+);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/thermophysicalModels/reactionThermo/combustionThermo/hsCombustionThermo/makeHsCombustionThermo.H

@@ -25,15 +25,13 @@ License
 InClass
     Foam::hsCombustionThermo
 
-Description
-
 \*---------------------------------------------------------------------------*/
 
 #ifndef makeHsCombustionThermo_H
 #define makeHsCombustionThermo_H
 
 #include "addToRunTimeSelectionTable.H"
-#include "basicSensiblePsiThermo.H"
+#include "basicPsiThermo.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -53,7 +51,7 @@ defineTemplateTypeNameAndDebugWithName                                        \
                                                                               \
 addToRunTimeSelectionTable                                                    \
 (                                                                             \
-    basicSensiblePsiThermo,                                                   \
+    basicPsiThermo,                                                           \
     MixtureThermo##Mixture##Transport##Thermo##EqnOfState,                    \
     fvMesh                                                                    \
 );                                                                            \
@@ -80,7 +78,7 @@ defineTemplateTypeNameAndDebugWithName                                        \
                                                                               \
 addToRunTimeSelectionTable                                                    \
 (                                                                             \
-    basicSensiblePsiThermo,                                                   \
+    basicPsiThermo,                                                           \
     MixtureThermo##Mixture##ThermoPhys,                                       \
     fvMesh                                                                    \
 );                                                                            \

src/thermophysicalModels/reactionThermo/combustionThermo/mixtureThermos/hPsiMixtureThermo/hPsiMixtureThermo.C

@@ -280,7 +280,7 @@ Foam::hPsiMixtureThermo<MixtureType>::Cp() const
                 IOobject::NO_WRITE
             ),
             mesh,
-            dimensionSet(0, 2, -2, -1, 0)
+            dimEnergy/dimMass/dimTemperature
         )
     );
 

src/thermophysicalModels/reactionThermo/combustionThermo/mixtureThermos/hPsiMixtureThermo/hPsiMixtureThermo.H

@@ -94,12 +94,12 @@ public:
         //- Update properties
         virtual void correct();
 
-        //- Chemical enthalpy [J/kg]
-        virtual tmp<volScalarField> hc() const;
-
 
         // Fields derived from thermodynamic state variables
 
+            //- Chemical enthalpy [J/kg]
+            virtual tmp<volScalarField> hc() const;
+
             //- Enthalpy for cell-set [J/kg]
             virtual tmp<scalarField> h
             (

src/thermophysicalModels/reactionThermo/combustionThermo/mixtureThermos/hhuMixtureThermo/hhuMixtureThermo.C

@@ -296,7 +296,7 @@ Foam::hhuMixtureThermo<MixtureType>::Cp() const
                 IOobject::NO_WRITE
             ),
             mesh,
-            dimensionSet(0, 2, -2, -1, 0)
+            dimEnergy/dimMass/dimTemperature
         )
     );
 

src/thermophysicalModels/reactionThermo/combustionThermo/mixtureThermos/hsPsiMixtureThermo/hsPsiMixtureThermo.C

@@ -205,15 +205,15 @@ Foam::hsPsiMixtureThermo<MixtureType>::hs
     const labelList& cells
 ) const
 {
-    tmp<scalarField> th(new scalarField(T.size()));
+    tmp<scalarField> ths(new scalarField(T.size()));
-    scalarField& hs = th();
+    scalarField& hs = ths();
 
     forAll(T, celli)
     {
         hs[celli] = this->cellMixture(cells[celli]).Hs(T[celli]);
     }
 
-    return th;
+    return ths;
 }
 
 
@@ -225,15 +225,15 @@ Foam::hsPsiMixtureThermo<MixtureType>::hs
     const label patchi
 ) const
 {
-    tmp<scalarField> th(new scalarField(T.size()));
+    tmp<scalarField> ths(new scalarField(T.size()));
-    scalarField& hs = th();
+    scalarField& hs = ths();
 
     forAll(T, facei)
     {
         hs[facei] = this->patchFaceMixture(patchi, facei).Hs(T[facei]);
     }
 
-    return th;
+    return ths;
 }
 
 
@@ -277,7 +277,7 @@ Foam::hsPsiMixtureThermo<MixtureType>::Cp() const
                 IOobject::NO_WRITE
             ),
             mesh,
-            dimensionSet(0, 2, -2, -1, 0)
+            dimEnergy/dimMass/dimTemperature
         )
     );
 

src/thermophysicalModels/reactionThermo/combustionThermo/mixtureThermos/hsPsiMixtureThermo/hsPsiMixtureThermo.H

@@ -95,11 +95,11 @@ public:
         virtual void correct();
 
 
+        // Fields derived from thermodynamic state variables
+
             //- Chemical enthalpy [J/kg]
             virtual tmp<volScalarField> hc() const;
 
-        // Fields derived from thermodynamic state variables
-
             //- Sensible enthalpy for cell-set [J/kg]
             virtual tmp<scalarField> hs
             (

src/thermophysicalModels/reactionThermo/reactionThermo/hReactionThermo/hReactionThermo.H

@@ -128,9 +128,6 @@ public:
             }
 
 
-        //- Chemical enthalpy [J/kg]
-        virtual tmp<volScalarField> hc() const = 0;
-
         //- Update properties
         virtual void correct() = 0;
 };

src/thermophysicalModels/reactionThermo/reactionThermo/hReactionThermo/makeReactionThermo.H

@@ -25,8 +25,6 @@ License
 InClass
     Foam::hReactionThermo
 
-Description
-
 \*---------------------------------------------------------------------------*/
 
 #ifndef makeReactionThermo_H

src/thermophysicalModels/reactionThermo/reactionThermo/hsReactionThermo/hsReactionThermo.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -24,28 +24,43 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "basicSensiblePsiThermo.H"
+#include "hsReactionThermo.H"
+#include "fvMesh.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
-    defineTypeNameAndDebug(basicSensiblePsiThermo, 0);
+    defineTypeNameAndDebug(hsReactionThermo, 0);
-    defineRunTimeSelectionTable(basicSensiblePsiThermo, fvMesh);
+    defineRunTimeSelectionTable(hsReactionThermo, fvMesh);
 }
 
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::basicSensiblePsiThermo::basicSensiblePsiThermo(const fvMesh& mesh)
+Foam::hsReactionThermo::hsReactionThermo(const fvMesh& mesh)
 :
-    basicThermo(mesh)
+    basicRhoThermo(mesh),
+
+    hs_
+    (
+        IOobject
+        (
+            "hs",
+            mesh.time().timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimEnergy/dimMass,
+        this->hBoundaryTypes()
+    )
 {}
 
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
-Foam::basicSensiblePsiThermo::~basicSensiblePsiThermo()
+Foam::hsReactionThermo::~hsReactionThermo()
 {}
 
 

src/thermophysicalModels/reactionThermo/reactionThermo/hsReactionThermo/hsReactionThermo.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -23,21 +23,22 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Class
-    Foam::basicSensiblePsiThermo
+    Foam::hsReactionThermo
 
 Description
-    Basic thermodynamic properties based on compressibility
+    Sensible enthalpy variant for reactionThermo
 
 SourceFiles
-    basicSensiblePsiThermo.C
+    hReactionThermo.C
-    newBasicSensiblePsiThermo.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef basicSensiblePsiThermo_H
+#ifndef hsReactionThermo_H
-#define basicSensiblePsiThermo_H
+#define hsReactionThermo_H
 
-#include "basicThermo.H"
+#include "basicRhoThermo.H"
+#include "basicMultiComponentMixture.H"
+#include "autoPtr.H"
 #include "runTimeSelectionTables.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -46,60 +47,89 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                        Class basicSensiblePsiThermo Declaration
+                     Class hsReactionThermo Declaration
 \*---------------------------------------------------------------------------*/
 
-class basicSensiblePsiThermo
+class hsReactionThermo
 :
-    public basicThermo
+    public basicRhoThermo
 {
 
 protected:
 
-    // Protected member functions
+    // Protected data
 
-        //- Construct as copy (not implemented)
+        //- Sensible enthalpy [J/kg]
-        basicSensiblePsiThermo(const basicSensiblePsiThermo&);
+        volScalarField hs_;
 
 
 public:
 
     //- Runtime type information
-    TypeName("basicSensiblePsiThermo");
+    TypeName("hsReactionThermo");
 
 
-    //- Declare run-time constructor selection table
+    //- Declare run-time constructor selection tables
     declareRunTimeSelectionTable
     (
         autoPtr,
-        basicSensiblePsiThermo,
+        hsReactionThermo,
         fvMesh,
         (const fvMesh& mesh),
         (mesh)
     );
 
+
     // Constructors
 
-        //- Construct from mesh
+        //- Construct from dictionary and mesh
-        basicSensiblePsiThermo(const fvMesh&);
+        hsReactionThermo(const fvMesh&);
 
 
-    //- Selector
+    // Selectors
-    static autoPtr<basicSensiblePsiThermo> New(const fvMesh&);
+
+        //- Standard selection based on fvMesh
+        static autoPtr<hsReactionThermo> New(const fvMesh&);
+
+        //- Select and check that package contains 'thermoType'
+        static autoPtr<hsReactionThermo> NewType
+        (
+            const fvMesh&,
+            const word& thermoType
+        );
+
 
     //- Destructor
-    virtual ~basicSensiblePsiThermo();
+    virtual ~hsReactionThermo();
 
 
     // Member functions
 
-        // Fields derived from thermodynamic state variables
+        //- Return the composition of the multi-component mixture
+        virtual basicMultiComponentMixture& composition() = 0;
 
-            //- Density [kg/m^3] - uses current value of pressure
+        //- Return the composition of the multi-component mixture
-            virtual tmp<volScalarField> rho() const
+        virtual const basicMultiComponentMixture& composition() const = 0;
+
+
+        // Access to thermodynamic state variables
+
+            //- Sensible enthalpy [J/kg]
+            //  Non-const access allowed for transport equations
+            virtual volScalarField& hs()
             {
-                return p_*psi();
+                return hs_;
             }
+
+            //- Sensible enthalpy [J/kg]
+            virtual const volScalarField& hs() const
+            {
+                return hs_;
+            }
+
+
+        //- Update properties
+        virtual void correct() = 0;
 };
 
 

src/thermophysicalModels/reactionThermo/reactionThermo/hsReactionThermo/hsReactionThermos.C

@@ -0,0 +1,171 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "makeHsReactionThermo.H"
+
+#include "hsReactionThermo.H"
+#include "hsRhoMixtureThermo.H"
+
+#include "perfectGas.H"
+
+#include "hConstThermo.H"
+#include "janafThermo.H"
+#include "specieThermo.H"
+
+#include "constTransport.H"
+#include "sutherlandTransport.H"
+
+#include "homogeneousMixture.H"
+#include "inhomogeneousMixture.H"
+#include "veryInhomogeneousMixture.H"
+#include "dieselMixture.H"
+#include "multiComponentMixture.H"
+#include "reactingMixture.H"
+
+#include "thermoPhysicsTypes.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+makeHsReactionThermo
+(
+    hsReactionThermo,
+    hsRhoMixtureThermo,
+    homogeneousMixture,
+    constTransport,
+    hConstThermo,
+    perfectGas
+);
+
+makeHsReactionThermo
+(
+    hsReactionThermo,
+    hsRhoMixtureThermo,
+    inhomogeneousMixture,
+    constTransport,
+    hConstThermo,
+    perfectGas
+);
+
+makeHsReactionThermo
+(
+    hsReactionThermo,
+    hsRhoMixtureThermo,
+    veryInhomogeneousMixture,
+    constTransport,
+    hConstThermo,
+    perfectGas
+);
+
+makeHsReactionThermo
+(
+    hsReactionThermo,
+    hsRhoMixtureThermo,
+    homogeneousMixture,
+    sutherlandTransport,
+    janafThermo,
+    perfectGas
+);
+
+makeHsReactionThermo
+(
+    hsReactionThermo,
+    hsRhoMixtureThermo,
+    inhomogeneousMixture,
+    sutherlandTransport,
+    janafThermo,
+    perfectGas
+);
+
+makeHsReactionThermo
+(
+    hsReactionThermo,
+    hsRhoMixtureThermo,
+    veryInhomogeneousMixture,
+    sutherlandTransport,
+    janafThermo,
+    perfectGas
+);
+
+
+makeHsReactionThermo
+(
+    hsReactionThermo,
+    hsRhoMixtureThermo,
+    dieselMixture,
+    sutherlandTransport,
+    janafThermo,
+    perfectGas
+);
+
+
+// Multi-component thermo
+
+makeHsReactionMixtureThermo
+(
+    hsReactionThermo,
+    hsRhoMixtureThermo,
+    multiComponentMixture,
+    icoPoly8ThermoPhysics
+);
+
+makeHsReactionMixtureThermo
+(
+    hsReactionThermo,
+    hsRhoMixtureThermo,
+    multiComponentMixture,
+    gasThermoPhysics
+);
+
+
+// Multi-component reaction thermo
+
+makeHsReactionMixtureThermo
+(
+    hsReactionThermo,
+    hsRhoMixtureThermo,
+    reactingMixture,
+    icoPoly8ThermoPhysics
+);
+
+makeHsReactionMixtureThermo
+(
+    hsReactionThermo,
+    hsRhoMixtureThermo,
+    reactingMixture,
+    gasThermoPhysics
+);
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //