merge of thermo branch into master

applications/solvers/combustion/PDRFoam/Make/options

@@ -8,7 +8,7 @@ EXE_INC = \
     -I$(LIB_SRC)/sampling/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/combustion/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
     -I$(LIB_SRC)/turbulenceModels \
     -I$(LIB_SRC)/turbulenceModels/compressible/RAS/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/laminarFlameSpeed/lnInclude \
@@ -23,7 +23,7 @@ EXE_LIBS = \
     -lmeshTools \
     -lcompressibleRASModels \
     -lbasicThermophysicalModels \
-    -lcombustionThermophysicalModels \
+    -lreactionThermophysicalModels \
     -lspecie \
     -llaminarFlameSpeedModels \
     -lfiniteVolume \

applications/solvers/combustion/PDRFoam/PDRFoam.C

@@ -36,7 +36,7 @@ Description
     to be appropriate by comparison with the results from the
     spectral model.
 
-    Strain effects are encorporated directly into the Xi equation
+    Strain effects are incorporated directly into the Xi equation
     but not in the algebraic approximation.  Further work need to be
     done on this issue, particularly regarding the enhanced removal rate
     caused by flame compression.  Analysis using results of the spectral
@@ -70,53 +70,53 @@ Description
 
 int main(int argc, char *argv[])
 {
-#   include "setRootCase.H"
+    #include "setRootCase.H"
 
-#   include "createTime.H"
+    #include "createTime.H"
-#   include "createMesh.H"
+    #include "createMesh.H"
-#   include "readCombustionProperties.H"
+    #include "readCombustionProperties.H"
-#   include "readEnvironmentalProperties.H"
+    #include "readEnvironmentalProperties.H"
-#   include "createFields.H"
+    #include "createFields.H"
-#   include "readPISOControls.H"
+    #include "readPISOControls.H"
-#   include "initContinuityErrs.H"
+    #include "initContinuityErrs.H"
-#   include "readTimeControls.H"
+    #include "readTimeControls.H"
-#   include "CourantNo.H"
+    #include "CourantNo.H"
-#   include "setInitialDeltaT.H"
+    #include "setInitialDeltaT.H"
 
-scalar StCoNum = 0.0;
+    scalar StCoNum = 0.0;
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
 
     while (runTime.run())
     {
-#       include "readTimeControls.H"
+        #include "readTimeControls.H"
-#       include "readPISOControls.H"
+        #include "readPISOControls.H"
-#       include "CourantNo.H"
+        #include "CourantNo.H"
-#       include "setDeltaT.H"
+        #include "setDeltaT.H"
 
         runTime++;
 
         Info<< "\n\nTime = " << runTime.timeName() << endl;
 
-#       include "rhoEqn.H"
+        #include "rhoEqn.H"
-#       include "UEqn.H"
+        #include "UEqn.H"
 
         // --- PISO loop
         for (int corr=1; corr<=nCorr; corr++)
         {
-#           include "bEqn.H"
+            #include "bEqn.H"
-#           include "ftEqn.H"
+            #include "ftEqn.H"
-#           include "huEqn.H"
+            #include "huEqn.H"
-#           include "hEqn.H"
+            #include "hEqn.H"
 
             if (!ign.ignited())
             {
                 hu == h;
             }
 
-#           include "pEqn.H"
+            #include "pEqn.H"
         }
 
         turbulence->correct();

applications/solvers/combustion/PDRFoam/StCourantNo.H

@@ -31,22 +31,24 @@ Description
 \*---------------------------------------------------------------------------*/
 
 {
-scalar meanStCoNum = 0.0;
+    scalar meanStCoNum = 0.0;
 
-if (mesh.nInternalFaces())
+    if (mesh.nInternalFaces())
-{
+    {
         surfaceScalarField SfUfbyDelta =
             mesh.surfaceInterpolation::deltaCoeffs()
         *mag(phiSt/fvc::interpolate(rho));
 
-    StCoNum = max(SfUfbyDelta/mesh.magSf())
+        StCoNum =
-        .value()*runTime.deltaT().value();
+            max(SfUfbyDelta/mesh.magSf()).value()
+           *runTime.deltaT().value();
 
-    meanStCoNum = (sum(SfUfbyDelta)/sum(mesh.magSf()))
+        meanStCoNum =
-        .value()*runTime.deltaT().value();
+            (sum(SfUfbyDelta)/sum(mesh.magSf())).value()
-}
+        *runTime.deltaT().value();
+    }
 
-Info<< "St courant Number mean: " << meanStCoNum
+    Info<< "St courant Number mean: " << meanStCoNum
         << " max: " << StCoNum << endl;
 }
 

applications/solvers/combustion/PDRFoam/bEqn.H

@@ -25,7 +25,7 @@ if (ign.ignited())
 
     // Unburnt gas density
     // ~~~~~~~~~~~~~~~~~~~
-    volScalarField rhou = thermo->rhou();
+    volScalarField rhou = thermo.rhou();
 
     // Calculate flame normal etc.
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~

applications/solvers/combustion/PDRFoam/createFields.H

@@ -1,10 +1,11 @@
     Info<< "Reading thermophysical properties\n" << endl;
 
-    autoPtr<hhuCombustionThermo> thermo
+    autoPtr<hhuCombustionThermo> pThermo
     (
         hhuCombustionThermo::New(mesh)
     );
-    combustionMixture& composition = thermo->composition();
+    hhuCombustionThermo& thermo = pThermo();
+    basicMultiComponentMixture& composition = thermo.composition();
 
     volScalarField rho
     (
@@ -16,13 +17,13 @@
             IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
-        thermo->rho()
+        thermo.rho()
     );
 
-    volScalarField& p = thermo->p();
+    volScalarField& p = thermo.p();
-    const volScalarField& psi = thermo->psi();
+    const volScalarField& psi = thermo.psi();
-    volScalarField& h = thermo->h();
+    volScalarField& h = thermo.h();
-    volScalarField& hu = thermo->hu();
+    volScalarField& hu = thermo.hu();
 
     volScalarField& b = composition.Y("b");
     Info<< "min(b) = " << min(b).value() << endl;
@@ -54,7 +55,7 @@
             rho,
             U,
             phi,
-            thermo()
+            thermo
         )
     );
 

applications/solvers/combustion/PDRFoam/hEqn.H

@@ -8,5 +8,5 @@
         betav*DpDt
     );
 
-    thermo->correct();
+    thermo.correct();
 }

applications/solvers/combustion/PDRFoam/huEqn.H

@@ -13,6 +13,6 @@ if (ign.ignited())
     //+ fvm::Sp(fvc::div(muEff*fvc::grad(b)/(b + 0.001)), hu)
 
      ==
-        betav*DpDt*rho/thermo->rhou()
+        betav*DpDt*rho/thermo.rhou()
     );
 }

applications/solvers/combustion/PDRFoam/laminarFlameSpeed/SCOPE/SCOPELaminarFlameSpeed.H

@@ -196,7 +196,6 @@ public:
 
 
     // Destructor
-
     ~SCOPE();
 
 

applications/solvers/combustion/PDRFoam/pEqn.H

@@ -1,4 +1,4 @@
-rho = thermo->rho();
+rho = thermo.rho();
 
 volScalarField rUA = 1.0/UEqn.A();
 U = invA & UEqn.H();
@@ -8,7 +8,7 @@ if (transonic)
     surfaceScalarField phid
     (
         "phid",
-        fvc::interpolate(thermo->psi())
+        fvc::interpolate(psi)
        *(
             (fvc::interpolate(U) & mesh.Sf())
           + fvc::ddtPhiCorr(rUA, rho, U, phi)

applications/solvers/combustion/XiFoam/Make/options

@@ -2,7 +2,7 @@ EXE_INC = \
     -I$(LIB_SRC)/engine/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/combustion/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
     -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
     -I$(LIB_SRC)/thermophysicalModels/laminarFlameSpeed/lnInclude \
     -I$(LIB_SRC)/sampling/lnInclude \
@@ -13,7 +13,7 @@ EXE_LIBS = \
     -lcompressibleRASModels \
     -lcompressibleLESModels \
     -lbasicThermophysicalModels \
-    -lcombustionThermophysicalModels \
+    -lreactionThermophysicalModels \
     -lspecie \
     -llaminarFlameSpeedModels \
     -lfiniteVolume \

applications/solvers/combustion/XiFoam/XiFoam.C

@@ -109,7 +109,7 @@ int main(int argc, char *argv[])
 
         turbulence->correct();
 
-        rho = thermo->rho();
+        rho = thermo.rho();
 
         runTime.write();
 

applications/solvers/combustion/XiFoam/bEqn.H

@@ -6,7 +6,7 @@ if (ign.ignited())
 
     // Unburnt gas density
     // ~~~~~~~~~~~~~~~~~~~
-    volScalarField rhou = thermo->rhou();
+    volScalarField rhou = thermo.rhou();
 
     // Calculate flame normal etc.
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -76,7 +76,7 @@ if (ign.ignited())
 
     volScalarField epsilon = pow(uPrimeCoef, 3)*turbulence->epsilon();
 
-    volScalarField tauEta = sqrt(thermo->muu()/(rhou*epsilon));
+    volScalarField tauEta = sqrt(thermo.muu()/(rhou*epsilon));
 
     volScalarField Reta = up/
     (

applications/solvers/combustion/XiFoam/createFields.H

@@ -1,10 +1,11 @@
     Info<< "Reading thermophysical properties\n" << endl;
 
-    autoPtr<hhuCombustionThermo> thermo
+    autoPtr<hhuCombustionThermo> pThermo
     (
         hhuCombustionThermo::New(mesh)
     );
-    combustionMixture& composition = thermo->composition();
+    hhuCombustionThermo& thermo = pThermo();
+    basicMultiComponentMixture& composition = thermo.composition();
 
     volScalarField rho
     (
@@ -16,18 +17,18 @@
             IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
-        thermo->rho()
+        thermo.rho()
     );
 
-    volScalarField& p = thermo->p();
+    volScalarField& p = thermo.p();
-    const volScalarField& psi = thermo->psi();
+    const volScalarField& psi = thermo.psi();
-    volScalarField& h = thermo->h();
+    volScalarField& h = thermo.h();
-    volScalarField& hu = thermo->hu();
+    volScalarField& hu = thermo.hu();
 
     volScalarField& b = composition.Y("b");
     Info<< "min(b) = " << min(b).value() << endl;
 
-    const volScalarField& T = thermo->T();
+    const volScalarField& T = thermo.T();
 
 
     Info<< "\nReading field U\n" << endl;
@@ -55,7 +56,7 @@
             rho,
             U,
             phi,
-            thermo()
+            thermo
         )
     );
 

applications/solvers/combustion/XiFoam/hEqn.H

@@ -8,5 +8,5 @@
         DpDt
     );
 
-    thermo->correct();
+    thermo.correct();
 }

applications/solvers/combustion/XiFoam/huEqn.H

@@ -13,6 +13,6 @@ if (ign.ignited())
     //+ fvm::Sp(fvc::div(muEff*fvc::grad(b)/(b + 0.001)), hu)
 
      ==
-        DpDt*rho/thermo->rhou()
+        DpDt*rho/thermo.rhou()
     );
 }

applications/solvers/combustion/XiFoam/pEqn.H

@@ -1,4 +1,4 @@
-rho = thermo->rho();
+rho = thermo.rho();
 
 volScalarField rUA = 1.0/UEqn.A();
 U = rUA*UEqn.H();
@@ -8,7 +8,7 @@ if (transonic)
     surfaceScalarField phid
     (
         "phid",
-        fvc::interpolate(thermo->psi())
+        fvc::interpolate(psi)
        *(
             (fvc::interpolate(U) & mesh.Sf())
           + fvc::ddtPhiCorr(rUA, rho, U, phi)
@@ -35,8 +35,8 @@ if (transonic)
 else
 {
     phi =
-        fvc::interpolate(rho)*
+        fvc::interpolate(rho)
-        (
+       *(
             (fvc::interpolate(U) & mesh.Sf())
           + fvc::ddtPhiCorr(rUA, rho, U, phi)
         );

applications/solvers/combustion/coalChemistryFoam/additionalOutput.H

@@ -1,48 +0,0 @@
-{
-    tmp<volScalarField> tdQ
-    (
-        new volScalarField
-        (
-            IOobject
-            (
-                "dQ",
-                runTime.timeName(),
-                mesh,
-                IOobject::NO_READ,
-                IOobject::AUTO_WRITE
-            ),
-            mesh,
-            dimensionedScalar
-            (
-                "zero",
-                dimensionSet(1, -3, -1, 0, 0, 0, 0),
-                0.0
-            )
-        )
-    );
-
-    scalarField& dQ = tdQ();
-
-    scalarField cp(dQ.size(), 0.0);
-
-    forAll(Y, i)
-    {
-        volScalarField RRi = chemistry.RR(i);
-
-        forAll(h, celli)
-        {
-            scalar Ti = T[celli];
-            cp[celli] += Y[i][celli]*chemistry.specieThermo()[i].Cp(Ti);
-            scalar hi = chemistry.specieThermo()[i].h(Ti);
-            scalar RR = RRi[celli];
-            dQ[celli] -= hi*RR;
-        }
-    }
-
-    forAll(dQ, celli)
-    {
-        dQ[celli] /= cp[celli];
-    }
-
-    tdQ().write();
-}

applications/solvers/combustion/coldEngineFoam/Make/options

@@ -4,7 +4,7 @@ EXE_INC = \
     -I$(LIB_SRC)/engine/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/combustion/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
     -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
     -I$(LIB_SRC)/finiteVolume/lnInclude
 

applications/solvers/combustion/coldEngineFoam/coldEngineFoam.C

@@ -33,7 +33,7 @@ Description
 #include "fvCFD.H"
 #include "engineTime.H"
 #include "engineMesh.H"
-#include "basicThermo.H"
+#include "basicPsiThermo.H"
 #include "turbulenceModel.H"
 #include "OFstream.H"
 
@@ -41,27 +41,27 @@ Description
 
 int main(int argc, char *argv[])
 {
-#   include "setRootCase.H"
+    #include "setRootCase.H"
 
-#   include "createEngineTime.H"
+    #include "createEngineTime.H"
-#   include "createEngineMesh.H"
+    #include "createEngineMesh.H"
-#   include "createFields.H"
+    #include "createFields.H"
-#   include "initContinuityErrs.H"
+    #include "initContinuityErrs.H"
-#   include "readEngineTimeControls.H"
+    #include "readEngineTimeControls.H"
-#   include "compressibleCourantNo.H"
+    #include "compressibleCourantNo.H"
-#   include "setInitialDeltaT.H"
+    #include "setInitialDeltaT.H"
-#   include "startSummary.H"
+    #include "startSummary.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
 
     while (runTime.run())
     {
-#       include "readPISOControls.H"
+        #include "readPISOControls.H"
-#       include "readEngineTimeControls.H"
+        #include "readEngineTimeControls.H"
-#       include "compressibleCourantNo.H"
+        #include "compressibleCourantNo.H"
-#       include "setDeltaT.H"
+        #include "setDeltaT.H"
 
         runTime++;
 
@@ -70,22 +70,22 @@ int main(int argc, char *argv[])
 
         mesh.move();
 
-#       include "rhoEqn.H"
+        #include "rhoEqn.H"
 
-#       include "UEqn.H"
+        #include "UEqn.H"
 
         // --- PISO loop
         for (int corr=1; corr<=nCorr; corr++)
         {
-#           include "hEqn.H"
+            #include "hEqn.H"
-#           include "pEqn.H"
+            #include "pEqn.H"
         }
 
         turbulence->correct();
 
         runTime.write();
 
-#       include "logSummary.H"
+        #include "logSummary.H"
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
             << "  ClockTime = " << runTime.elapsedClockTime() << " s"

applications/solvers/combustion/coldEngineFoam/createFields.H

@@ -1,9 +1,10 @@
     Info<< "Reading thermophysical properties\n" << endl;
 
-    autoPtr<basicThermo> thermo
+    autoPtr<basicPsiThermo> pThermo
     (
-        basicThermo::New(mesh)
+        basicPsiThermo::New(mesh)
     );
+    basicPsiThermo& thermo = pThermo();
 
     volScalarField rho
     (
@@ -15,13 +16,13 @@
             IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
-        thermo->rho()
+        thermo.rho()
     );
 
-    volScalarField& p = thermo->p();
+    volScalarField& p = thermo.p();
-    const volScalarField& psi = thermo->psi();
+    const volScalarField& psi = thermo.psi();
-    volScalarField& h = thermo->h();
+    volScalarField& h = thermo.h();
-    const volScalarField& T = thermo->T();
+    const volScalarField& T = thermo.T();
 
 
     Info<< "\nReading field U\n" << endl;
@@ -38,7 +39,7 @@
         mesh
     );
 
-#   include "compressibleCreatePhi.H"
+    #include "compressibleCreatePhi.H"
 
 
     Info<< "Creating turbulence model\n" << endl;
@@ -49,7 +50,7 @@
             rho,
             U,
             phi,
-            thermo()
+            thermo
         )
     );
 

applications/solvers/combustion/coldEngineFoam/hEqn.H

@@ -8,5 +8,5 @@
         DpDt
     );
 
-    thermo->correct();
+    thermo.correct();
 }

applications/solvers/combustion/dieselEngineFoam/Make/options

@@ -7,10 +7,10 @@ EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/liquidMixture/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/thermophysicalFunctions/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/combustion/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/laminarFlameSpeed/lnInclude \
-    -I$(LIB_SRC)/../applications/solvers/combustion/XiFoam \
+    -I$(LIB_SRC)/../applications/solvers/reactionThermo/XiFoam \
     -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
     -I$(LIB_SRC)/ODE/lnInclude \
     -I$(LIB_SRC)/engine/lnInclude \
@@ -20,7 +20,7 @@ EXE_LIBS = \
     -lengine \
     -lcompressibleRASModels \
     -lcompressibleLESModels \
-    -lcombustionThermophysicalModels \
+    -lreactionThermophysicalModels \
     -lfiniteVolume \
     -llagrangian \
     -ldieselSpray \

applications/solvers/combustion/dieselEngineFoam/YEqn.H

@@ -14,7 +14,7 @@ tmp<fv::convectionScheme<scalar> > mvConvection
     label inertIndex = -1;
     volScalarField Yt = 0.0*Y[0];
 
-    for(label i=0; i<Y.size(); i++)
+    for (label i=0; i<Y.size(); i++)
     {
         if (Y[i].name() != inertSpecie)
         {
@@ -42,5 +42,4 @@ tmp<fv::convectionScheme<scalar> > mvConvection
 
     Y[inertIndex] = scalar(1) - Yt;
     Y[inertIndex].max(0.0);
-
 }

applications/solvers/combustion/dieselEngineFoam/createFields.H

@@ -1,13 +1,17 @@
 Info<< nl << "Reading thermophysicalProperties" << endl;
-autoPtr<hCombustionThermo> thermo
-(
-    hCombustionThermo::New(mesh)
-);
 
-combustionMixture& composition = thermo->composition();
+autoPtr<psiChemistryModel> pChemistry
+(
+    psiChemistryModel::New(mesh)
+);
+psiChemistryModel& chemistry = pChemistry();
+
+hCombustionThermo& thermo = chemistry.thermo();
+
+basicMultiComponentMixture& composition = thermo.composition();
 PtrList<volScalarField>& Y = composition.Y();
 
-word inertSpecie(thermo->lookup("inertSpecie"));
+word inertSpecie(thermo.lookup("inertSpecie"));
 
 volScalarField rho
 (
@@ -17,7 +21,7 @@ volScalarField rho
         runTime.timeName(),
         mesh
     ),
-    thermo->rho()
+    thermo.rho()
 );
 
 Info<< "Reading field U\n" << endl;
@@ -35,10 +39,10 @@ volVectorField U
 );
 
 
-volScalarField& p = thermo->p();
+volScalarField& p = thermo.p();
-const volScalarField& psi = thermo->psi();
+const volScalarField& psi = thermo.psi();
-const volScalarField& T = thermo->T();
+const volScalarField& T = thermo.T();
-volScalarField& h = thermo->h();
+volScalarField& h = thermo.h();
 
 
 #include "compressibleCreatePhi.H"
@@ -65,7 +69,7 @@ autoPtr<compressible::turbulenceModel> turbulence
         rho,
         U,
         phi,
-        thermo()
+        thermo
     )
 );
 
@@ -73,31 +77,11 @@ Info<< "Creating field DpDt\n" << endl;
 volScalarField DpDt =
     fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
 
-Info << "Constructing chemical mechanism" << endl;
-chemistryModel chemistry
-(
-    thermo(),
-    rho
-);
 
 multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
 
-for(label i=0; i<Y.size(); i++)
+forAll (Y, i)
 {
     fields.add(Y[i]);
 }
 fields.add(h);
-
-volScalarField dQ
-(
-    IOobject
-    (
-        "dQ",
-        runTime.timeName(),
-        mesh,
-        IOobject::NO_READ,
-        IOobject::AUTO_WRITE
-    ),
-    mesh,
-    dimensionedScalar("zero", dimensionSet(1,-3,-1,0,0,0,0), 0.0)
-);

applications/solvers/combustion/dieselEngineFoam/createSpray.H

@@ -1,14 +1,15 @@
 Info << "Constructing Spray" << endl;
 
-PtrList<specieProperties> gasProperties(Y.size());
+PtrList<gasThermoPhysics> gasProperties(Y.size());
 forAll(gasProperties, i)
 {
     gasProperties.set
     (
         i,
-        new specieProperties
+        new gasThermoPhysics
         (
-            dynamic_cast<const reactingMixture&>(thermo()).speciesData()[i]
+            dynamic_cast<const reactingMixture<gasThermoPhysics>&>
+                (thermo).speciesData()[i]
         )
     );
 }

applications/solvers/combustion/dieselEngineFoam/dieselEngineFoam.C

@@ -36,40 +36,41 @@ Description
 #include "hCombustionThermo.H"
 #include "turbulenceModel.H"
 #include "spray.H"
-#include "chemistryModel.H"
+#include "psiChemistryModel.H"
 #include "chemistrySolver.H"
 #include "multivariateScheme.H"
 #include "Switch.H"
 #include "OFstream.H"
 #include "volPointInterpolation.H"
+#include "thermoPhysicsTypes.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 int main(int argc, char *argv[])
 {
-#   include "setRootCase.H"
+    #include "setRootCase.H"
-#   include "createEngineTime.H"
+    #include "createEngineTime.H"
-#   include "createEngineMesh.H"
+    #include "createEngineMesh.H"
-#   include "createFields.H"
+    #include "createFields.H"
-#   include "readEnvironmentalProperties.H"
+    #include "readEnvironmentalProperties.H"
-#   include "readCombustionProperties.H"
+    #include "readCombustionProperties.H"
-#   include "createSpray.H"
+    #include "createSpray.H"
-#   include "initContinuityErrs.H"
+    #include "initContinuityErrs.H"
-#   include "readEngineTimeControls.H"
+    #include "readEngineTimeControls.H"
-#   include "compressibleCourantNo.H"
+    #include "compressibleCourantNo.H"
-#   include "setInitialDeltaT.H"
+    #include "setInitialDeltaT.H"
-#   include "startSummary.H"
+    #include "startSummary.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info << "\nStarting time loop\n" << endl;
 
     while (runTime.run())
     {
-#       include "readPISOControls.H"
+        #include "readPISOControls.H"
-#       include "readEngineTimeControls.H"
+        #include "readEngineTimeControls.H"
-#       include "compressibleCourantNo.H"
+        #include "compressibleCourantNo.H"
-#       include "setDeltaT.H"
+        #include "setDeltaT.H"
 
         runTime++;
 
@@ -101,27 +102,27 @@ int main(int argc, char *argv[])
             kappa = (runTime.deltaT() + tc)/(runTime.deltaT() + tc + tk);
         }
 
-#       include "rhoEqn.H"
+        #include "rhoEqn.H"
-#       include "UEqn.H"
+        #include "UEqn.H"
 
         for (label ocorr=1; ocorr <= nOuterCorr; ocorr++)
         {
-#           include "YEqn.H"
+            #include "YEqn.H"
-#           include "hEqn.H"
+            #include "hEqn.H"
 
             // --- PISO loop
             for (int corr=1; corr<=nCorr; corr++)
             {
-#               include "pEqn.H"
+                #include "pEqn.H"
             }
         }
 
         turbulence->correct();
 
-#       include "logSummary.H"
+        #include "logSummary.H"
-#       include "spraySummary.H"
+        #include "spraySummary.H"
 
-        rho = thermo->rho();
+        rho = thermo.rho();
 
         runTime.write();
 

applications/solvers/combustion/dieselEngineFoam/hEqn.H

@@ -9,32 +9,5 @@
      + dieselSpray.heatTransferSource()
     );
 
-    thermo->correct();
+    thermo.correct();
-
-    forAll(dQ, i)
-    {
-        dQ[i] = 0.0;
-    }
-
-    scalarField cp(dQ.size(), 0.0);
-
-    forAll(Y, i)
-    {
-        volScalarField RRi = chemistry.RR(i);
-
-        forAll(h, celli)
-        {
-            scalar Ti = T[celli];
-            cp[celli] += Y[i][celli]*chemistry.specieThermo()[i].Cp(Ti);
-            scalar hi = chemistry.specieThermo()[i].h(Ti);
-            scalar RR = RRi[celli];
-            dQ[celli] -= hi*RR;
-        }
-
-    }
-
-    forAll(dQ, celli)
-    {
-        dQ[celli] /= cp[celli];
-    }
 }

applications/solvers/combustion/dieselEngineFoam/pEqn.H

@@ -1,4 +1,4 @@
-rho = thermo->rho();
+rho = thermo.rho();
 
 volScalarField A = UEqn.A();
 U = UEqn.H()/A;
@@ -8,7 +8,7 @@ if (transonic)
     surfaceScalarField phid
     (
         "phid",
-        fvc::interpolate(thermo->psi())
+        fvc::interpolate(psi)
        *((fvc::interpolate(U) & mesh.Sf()) - fvc::meshPhi(rho, U))
     );
 

applications/solvers/combustion/dieselEngineFoam/rhoEqn.H

@@ -44,7 +44,7 @@ volScalarField Sevap
     dimensionedScalar("zero", dimensionSet(1, -3, -1, 0, 0), 0.0)
 );
 
-for(label i=0; i<Y.size(); i++)
+for (label i=0; i<Y.size(); i++)
 {
     if (dieselSpray.isLiquidFuel()[i])
     {

applications/solvers/combustion/dieselFoam/Make/options

@@ -8,17 +8,17 @@ EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/liquidMixture/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/thermophysicalFunctions/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/combustion/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/laminarFlameSpeed/lnInclude \
-    -I$(LIB_SRC)/../applications/solvers/combustion/XiFoam \
+    -I$(LIB_SRC)/../applications/solvers/reactionThermo/XiFoam \
     -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
     -I$(LIB_SRC)/ODE/lnInclude
 
 EXE_LIBS = \
     -lcompressibleRASModels \
     -lcompressibleLESModels \
-    -lcombustionThermophysicalModels \
+    -lreactionThermophysicalModels \
     -llagrangian \
     -ldieselSpray \
     -lliquids \

applications/solvers/combustion/dieselFoam/dieselFoam.C

@@ -34,7 +34,7 @@ Description
 #include "hCombustionThermo.H"
 #include "turbulenceModel.H"
 #include "spray.H"
-#include "chemistryModel.H"
+#include "psiChemistryModel.H"
 #include "chemistrySolver.H"
 
 #include "multivariateScheme.H"
@@ -46,28 +46,27 @@ Description
 
 int main(int argc, char *argv[])
 {
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "createFields.H"
+    #include "readEnvironmentalProperties.H"
+    #include "readCombustionProperties.H"
+    #include "createSpray.H"
+    #include "initContinuityErrs.H"
+    #include "readTimeControls.H"
+    #include "compressibleCourantNo.H"
+    #include "setInitialDeltaT.H"
 
-#   include "setRootCase.H"
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-#   include "createTime.H"
-#   include "createMesh.H"
-#   include "createFields.H"
-#   include "readEnvironmentalProperties.H"
-#   include "readCombustionProperties.H"
-#   include "createSpray.H"
-#   include "initContinuityErrs.H"
-#   include "readTimeControls.H"
-#   include "compressibleCourantNo.H"
-#   include "setInitialDeltaT.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info << "\nStarting time loop\n" << endl;
 
     while (runTime.run())
     {
-#       include "readPISOControls.H"
+        #include "readPISOControls.H"
-#       include "compressibleCourantNo.H"
+        #include "compressibleCourantNo.H"
-#       include "setDeltaT.H"
+        #include "setDeltaT.H"
 
         runTime++;
         Info<< "Time = " << runTime.timeName() << nl << endl;
@@ -94,26 +93,26 @@ int main(int argc, char *argv[])
             kappa = (runTime.deltaT() + tc)/(runTime.deltaT()+tc+tk);
         }
 
-#       include "rhoEqn.H"
+        #include "rhoEqn.H"
-#       include "UEqn.H"
+        #include "UEqn.H"
 
         for (label ocorr=1; ocorr <= nOuterCorr; ocorr++)
         {
-#           include "YEqn.H"
+            #include "YEqn.H"
-#           include "hEqn.H"
+            #include "hEqn.H"
 
             // --- PISO loop
             for (int corr=1; corr<=nCorr; corr++)
             {
-#               include "pEqn.H"
+                #include "pEqn.H"
             }
         }
 
         turbulence->correct();
 
-#       include "spraySummary.H"
+        #include "spraySummary.H"
 
-        rho = thermo->rho();
+        rho = thermo.rho();
 
         runTime.write();
 

applications/solvers/combustion/dieselFoam/pEqn.H

@@ -1,4 +1,4 @@
-rho = thermo->rho();
+rho = thermo.rho();
 
 volScalarField rUA = 1.0/UEqn.A();
 U = rUA*UEqn.H();
@@ -8,7 +8,7 @@ if (transonic)
     surfaceScalarField phid
     (
         "phid",
-        fvc::interpolate(thermo->psi())
+        fvc::interpolate(psi)
        *(
             (fvc::interpolate(U) & mesh.Sf())
           + fvc::ddtPhiCorr(rUA, rho, U, phi)
@@ -37,8 +37,8 @@ if (transonic)
 else
 {
     phi =
-        fvc::interpolate(rho)*
+        fvc::interpolate(rho)
-        (
+       *(
             (fvc::interpolate(U) & mesh.Sf())
           + fvc::ddtPhiCorr(rUA, rho, U, phi)
         );

applications/solvers/combustion/engineFoam/Make/options

@@ -3,7 +3,7 @@ EXE_INC = \
     -I$(LIB_SRC)/engine/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/combustion/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
     -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
     -I$(LIB_SRC)/thermophysicalModels/laminarFlameSpeed/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude
@@ -13,7 +13,7 @@ EXE_LIBS = \
     -lcompressibleRASModels \
     -lcompressibleLESModels \
     -lbasicThermophysicalModels \
-    -lcombustionThermophysicalModels \
+    -lreactionThermophysicalModels \
     -lspecie \
     -llaminarFlameSpeedModels \
     -lfiniteVolume

applications/solvers/combustion/engineFoam/engineFoam.C

@@ -63,29 +63,29 @@ Description
 
 int main(int argc, char *argv[])
 {
-#   include "setRootCase.H"
+    #include "setRootCase.H"
 
-#   include "createEngineTime.H"
+    #include "createEngineTime.H"
-#   include "createEngineMesh.H"
+    #include "createEngineMesh.H"
-#   include "readPISOControls.H"
+    #include "readPISOControls.H"
-#   include "readCombustionProperties.H"
+    #include "readCombustionProperties.H"
-#   include "createFields.H"
+    #include "createFields.H"
-#   include "initContinuityErrs.H"
+    #include "initContinuityErrs.H"
-#   include "readEngineTimeControls.H"
+    #include "readEngineTimeControls.H"
-#   include "compressibleCourantNo.H"
+    #include "compressibleCourantNo.H"
-#   include "setInitialDeltaT.H"
+    #include "setInitialDeltaT.H"
-#   include "startSummary.H"
+    #include "startSummary.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info << "\nStarting time loop\n" << endl;
 
     while (runTime.run())
     {
-#       include "readPISOControls.H"
+        #include "readPISOControls.H"
-#       include "readEngineTimeControls.H"
+        #include "readEngineTimeControls.H"
-#       include "compressibleCourantNo.H"
+        #include "compressibleCourantNo.H"
-#       include "setDeltaT.H"
+        #include "setDeltaT.H"
 
         runTime++;
 
@@ -93,31 +93,31 @@ int main(int argc, char *argv[])
 
         mesh.move();
 
-#       include "rhoEqn.H"
+        #include "rhoEqn.H"
 
-#       include "UEqn.H"
+        #include "UEqn.H"
 
         // --- PISO loop
         for (int corr=1; corr<=nCorr; corr++)
         {
-#           include "ftEqn.H"
+            #include "ftEqn.H"
-#           include "bEqn.H"
+            #include "bEqn.H"
-#           include "huEqn.H"
+            #include "huEqn.H"
-#           include "hEqn.H"
+            #include "hEqn.H"
 
             if (!ign.ignited())
             {
                 hu == h;
             }
 
-#           include "pEqn.H"
+            #include "pEqn.H"
         }
 
         turbulence->correct();
 
-#       include "logSummary.H"
+        #include "logSummary.H"
 
-        rho = thermo->rho();
+        rho = thermo.rho();
 
         runTime.write();
 

applications/solvers/combustion/engineFoam/pEqn.H

@@ -1,4 +1,4 @@
-rho = thermo->rho();
+rho = thermo.rho();
 
 volScalarField rUA = 1.0/UEqn.A();
 U = rUA*UEqn.H();
@@ -8,7 +8,7 @@ if (transonic)
     surfaceScalarField phid
     (
         "phid",
-        fvc::interpolate(thermo->psi())
+        fvc::interpolate(psi)
        *((fvc::interpolate(U) & mesh.Sf()) - fvc::meshPhi(rho, U))
     );
 

applications/solvers/combustion/engineFoam/readEngineTimeControls.H

@@ -1,3 +1,3 @@
-#include "readTimeControls.H"
+    #include "readTimeControls.H"
 
     maxDeltaT = runTime.userTimeToTime(maxDeltaT);

applications/solvers/combustion/reactingFoam/Make/options

@@ -2,7 +2,7 @@ EXE_INC = \
     -I../XiFoam \
     -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/combustion/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
     -I$(LIB_SRC)/ODE/lnInclude \
@@ -11,7 +11,7 @@ EXE_INC = \
 EXE_LIBS = \
     -lcompressibleRASModels \
     -lcompressibleLESModels \
-    -lcombustionThermophysicalModels \
+    -lreactionThermophysicalModels \
     -lspecie \
     -lbasicThermophysicalModels \
     -lchemistryModel \

applications/solvers/combustion/reactingFoam/YEqn.H

@@ -13,7 +13,7 @@ tmp<fv::convectionScheme<scalar> > mvConvection
     label inertIndex = -1;
     volScalarField Yt = 0.0*Y[0];
 
-    for(label i=0; i<Y.size(); i++)
+    for (label i=0; i<Y.size(); i++)
     {
         if (Y[i].name() != inertSpecie)
         {

applications/solvers/combustion/reactingFoam/createFields.H

@@ -1,13 +1,16 @@
 Info<< nl << "Reading thermophysicalProperties" << endl;
-autoPtr<hCombustionThermo> thermo
+autoPtr<psiChemistryModel> pChemistry
 (
-    hCombustionThermo::New(mesh)
+    psiChemistryModel::New(mesh)
 );
+psiChemistryModel& chemistry = pChemistry();
 
-combustionMixture& composition = thermo->composition();
+hCombustionThermo& thermo = chemistry.thermo();
+
+basicMultiComponentMixture& composition = thermo.composition();
 PtrList<volScalarField>& Y = composition.Y();
 
-word inertSpecie(thermo->lookup("inertSpecie"));
+word inertSpecie(thermo.lookup("inertSpecie"));
 
 volScalarField rho
 (
@@ -17,7 +20,7 @@ volScalarField rho
         runTime.timeName(),
         mesh
     ),
-    thermo->rho()
+    thermo.rho()
 );
 
 Info<< "Reading field U\n" << endl;
@@ -35,9 +38,9 @@ volVectorField U
 );
 
 
-volScalarField& p = thermo->p();
+volScalarField& p = thermo.p();
-const volScalarField& psi = thermo->psi();
+const volScalarField& psi = thermo.psi();
-volScalarField& h = thermo->h();
+volScalarField& h = thermo.h();
 
 
 #include "compressibleCreatePhi.H"
@@ -64,7 +67,7 @@ autoPtr<compressible::turbulenceModel> turbulence
         rho,
         U,
         phi,
-        thermo()
+        thermo
     )
 );
 
@@ -72,31 +75,11 @@ Info<< "Creating field DpDt\n" << endl;
 volScalarField DpDt =
     fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
 
-Info << "Constructing chemical mechanism" << endl;
-chemistryModel chemistry
-(
-    thermo(),
-    rho
-);
-
 multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
 
-for(label i=0; i<Y.size(); i++)
+forAll (Y, i)
 {
     fields.add(Y[i]);
 }
 fields.add(h);
 
-volScalarField dQ
-(
-    IOobject
-    (
-        "dQ",
-        runTime.timeName(),
-        mesh,
-        IOobject::NO_READ,
-        IOobject::AUTO_WRITE
-    ),
-    mesh,
-    dimensionedScalar("zero", dimensionSet(1, -3, -1, 0, 0, 0, 0), 0.0)
-);

applications/solvers/combustion/reactingFoam/reactingFoam.C

@@ -33,7 +33,7 @@ Description
 #include "fvCFD.H"
 #include "hCombustionThermo.H"
 #include "turbulenceModel.H"
-#include "chemistryModel.H"
+#include "psiChemistryModel.H"
 #include "chemistrySolver.H"
 #include "multivariateScheme.H"
 
@@ -41,52 +41,52 @@ Description
 
 int main(int argc, char *argv[])
 {
-#   include "setRootCase.H"
+    #include "setRootCase.H"
-#   include "createTime.H"
+    #include "createTime.H"
-#   include "createMesh.H"
+    #include "createMesh.H"
-#   include "readChemistryProperties.H"
+    #include "readChemistryProperties.H"
-#   include "readEnvironmentalProperties.H"
+    #include "readEnvironmentalProperties.H"
-#   include "createFields.H"
+    #include "createFields.H"
-#   include "initContinuityErrs.H"
+    #include "initContinuityErrs.H"
-#   include "readTimeControls.H"
+    #include "readTimeControls.H"
-#   include "compressibleCourantNo.H"
+    #include "compressibleCourantNo.H"
-#   include "setInitialDeltaT.H"
+    #include "setInitialDeltaT.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info << "\nStarting time loop\n" << endl;
 
     while (runTime.run())
     {
-#       include "readTimeControls.H"
+        #include "readTimeControls.H"
-#       include "readPISOControls.H"
+        #include "readPISOControls.H"
-#       include "compressibleCourantNo.H"
+        #include "compressibleCourantNo.H"
-#       include "setDeltaT.H"
+        #include "setDeltaT.H"
 
         runTime++;
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-#       include "chemistry.H"
+        #include "chemistry.H"
-#       include "rhoEqn.H"
+        #include "rhoEqn.H"
-#       include "UEqn.H"
+        #include "UEqn.H"
 
         for (label ocorr=1; ocorr <= nOuterCorr; ocorr++)
         {
-#           include "YEqn.H"
+            #include "YEqn.H"
 
-#           define Db turbulence->alphaEff()
+            #define Db turbulence->alphaEff()
-#           include "hEqn.H"
+            #include "hEqn.H"
 
             // --- PISO loop
             for (int corr=1; corr<=nCorr; corr++)
             {
-#               include "pEqn.H"
+                #include "pEqn.H"
             }
         }
 
         turbulence->correct();
 
-        rho = thermo->rho();
+        rho = thermo.rho();
 
         runTime.write();
 

applications/solvers/combustion/reactingFoam/readChemistryProperties.H

@@ -8,7 +8,8 @@ IOdictionary chemistryProperties
         runTime.constant(),
         mesh,
         IOobject::MUST_READ,
-        IOobject::NO_WRITE
+        IOobject::NO_WRITE,
+        false
     )
 );
 

applications/solvers/compressible/rhoCentralFoam/compressibleCourantNo.H

@@ -38,11 +38,11 @@ if (mesh.nInternalFaces())
     surfaceScalarField amaxSfbyDelta =
         mesh.surfaceInterpolation::deltaCoeffs()*amaxSf;
 
-    CoNum = max(amaxSfbyDelta/mesh.magSf())
+    CoNum = max(amaxSfbyDelta/mesh.magSf()).value()*runTime.deltaT().value();
-        .value()*runTime.deltaT().value();
 
-    meanCoNum = (sum(amaxSfbyDelta)/sum(mesh.magSf()))
+    meanCoNum =
-        .value()*runTime.deltaT().value();
+        (sum(amaxSfbyDelta)/sum(mesh.magSf())).value()
+       *runTime.deltaT().value();
 }
 
 Info<< "Mean and max Courant Numbers = "

applications/solvers/compressible/rhoCentralFoam/createFields.H

@@ -1,15 +1,16 @@
 Info<< "Reading thermophysical properties\n" << endl;
 
-autoPtr<basicThermo> thermo
+autoPtr<basicPsiThermo> pThermo
 (
-    basicThermo::New(mesh)
+    basicPsiThermo::New(mesh)
 );
+basicPsiThermo& thermo = pThermo();
 
-volScalarField& p = thermo->p();
+volScalarField& p = thermo.p();
-volScalarField& h = thermo->h();
+volScalarField& h = thermo.h();
-const volScalarField& T = thermo->T();
+const volScalarField& T = thermo.T();
-const volScalarField& psi = thermo->psi();
+const volScalarField& psi = thermo.psi();
-const volScalarField& mu = thermo->mu();
+const volScalarField& mu = thermo.mu();
 
 bool inviscid(true);
 if (max(mu.internalField()) > 0.0)
@@ -42,7 +43,7 @@ volScalarField rho
         IOobject::NO_READ,
         IOobject::AUTO_WRITE
     ),
-    thermo->rho(),
+    thermo.rho(),
     rhoBoundaryTypes
 );
 

applications/solvers/compressible/rhoCentralFoam/rhoBoundaryTypes.H

@@ -3,10 +3,7 @@ wordList rhoBoundaryTypes = pbf.types();
 
 forAll(rhoBoundaryTypes, patchi)
 {
-    if
+    if (rhoBoundaryTypes[patchi] == "waveTransmissive")
-    (
-        rhoBoundaryTypes[patchi] == "waveTransmissive"
-    )
     {
         rhoBoundaryTypes[patchi] = zeroGradientFvPatchScalarField::typeName;
     }

applications/solvers/compressible/rhoCentralFoam/rhoCentralFoam.C

@@ -32,7 +32,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "basicThermo.H"
+#include "basicPsiThermo.H"
 #include "zeroGradientFvPatchFields.H"
 #include "fixedRhoFvPatchScalarField.H"
 
@@ -40,18 +40,17 @@ Description
 
 int main(int argc, char *argv[])
 {
+    #include "setRootCase.H"
 
-#   include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "createFields.H"
+    #include "readThermophysicalProperties.H"
+    #include "readTimeControls.H"
 
-#   include "createTime.H"
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-#   include "createMesh.H"
-#   include "createFields.H"
-#   include "readThermophysicalProperties.H"
-#   include "readTimeControls.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+    #include "readFluxScheme.H"
-
-#   include "readFluxScheme.H"
 
     dimensionedScalar v_zero("v_zero",dimVolume/dimTime, 0.0);
 
@@ -91,7 +90,7 @@ int main(int argc, char *argv[])
         surfaceScalarField phiv_pos = U_pos & mesh.Sf();
         surfaceScalarField phiv_neg = U_neg & mesh.Sf();
 
-        volScalarField c = sqrt(thermo->Cp()/thermo->Cv()*rPsi);
+        volScalarField c = sqrt(thermo.Cp()/thermo.Cv()*rPsi);
         surfaceScalarField cSf_pos = fvc::interpolate(c, pos, "reconstruct(T)")*mesh.magSf();
         surfaceScalarField cSf_neg = fvc::interpolate(c, neg, "reconstruct(T)")*mesh.magSf();
 
@@ -102,9 +101,9 @@ int main(int argc, char *argv[])
 
         surfaceScalarField amaxSf("amaxSf", max(mag(am), mag(ap)));
 
-#       include "compressibleCourantNo.H"
+        #include "compressibleCourantNo.H"
-#       include "readTimeControls.H"
+        #include "readTimeControls.H"
-#       include "setDeltaT.H"
+        #include "setDeltaT.H"
 
         runTime++;
 
@@ -183,7 +182,7 @@ int main(int argc, char *argv[])
 
         h = (rhoE + p)/rho - 0.5*magSqr(U);
         h.correctBoundaryConditions();
-        thermo->correct();
+        thermo.correct();
         rhoE.boundaryField() =
             rho.boundaryField()*
             (
@@ -193,15 +192,15 @@ int main(int argc, char *argv[])
 
         if (!inviscid)
         {
-            volScalarField k("k", thermo->Cp()*mu/Pr);
+            volScalarField k("k", thermo.Cp()*mu/Pr);
             solve
             (
                 fvm::ddt(rho, h) - fvc::ddt(rho, h)
-              - fvm::laplacian(thermo->alpha(), h)
+              - fvm::laplacian(thermo.alpha(), h)
-              + fvc::laplacian(thermo->alpha(), h)
+              + fvc::laplacian(thermo.alpha(), h)
               - fvc::laplacian(k, T)
             );
-            thermo->correct();
+            thermo.correct();
             rhoE = rho*(h + 0.5*magSqr(U)) - p;
         }
 

applications/solvers/compressible/rhoPimpleFoam/createFields.H

@@ -1,13 +1,14 @@
     Info<< "Reading thermophysical properties\n" << endl;
 
-    autoPtr<basicThermo> thermo
+    autoPtr<basicPsiThermo> pThermo
     (
-        basicThermo::New(mesh)
+        basicPsiThermo::New(mesh)
     );
+    basicPsiThermo& thermo = pThermo();
 
-    volScalarField& p = thermo->p();
+    volScalarField& p = thermo.p();
-    volScalarField& h = thermo->h();
+    volScalarField& h = thermo.h();
-    const volScalarField& psi = thermo->psi();
+    const volScalarField& psi = thermo.psi();
 
     volScalarField rho
     (
@@ -19,7 +20,7 @@
             IOobject::READ_IF_PRESENT,
             IOobject::AUTO_WRITE
         ),
-        thermo->rho()
+        thermo.rho()
     );
 
     Info<< "Reading field U\n" << endl;
@@ -51,7 +52,7 @@
             rho,
             U,
             phi,
-            thermo()
+            thermo
         )
     );
 

applications/solvers/compressible/rhoPimpleFoam/hEqn.H

@@ -19,5 +19,5 @@
         hEqn.solve();
     }
 
-    thermo->correct();
+    thermo.correct();
 }

applications/solvers/compressible/rhoPimpleFoam/pEqn.H

@@ -1,4 +1,4 @@
-rho = thermo->rho();
+rho = thermo.rho();
 
 volScalarField rUA = 1.0/UEqn().A();
 U = rUA*UEqn().H();
@@ -13,7 +13,7 @@ if (transonic)
     surfaceScalarField phid
     (
         "phid",
-        fvc::interpolate(thermo->psi())
+        fvc::interpolate(psi)
        *(
             (fvc::interpolate(U) & mesh.Sf())
           + fvc::ddtPhiCorr(rUA, rho, U, phi)
@@ -99,7 +99,7 @@ else
     // Explicitly relax pressure for momentum corrector
     p.relax();
 
-    rho = thermo->rho();
+    rho = thermo.rho();
     rho.relax();
     Info<< "rho max/min : " << max(rho).value()
         << " " << min(rho).value() << endl;
@@ -117,7 +117,7 @@ bound(p, pMin);
 /*
 if (closedVolume)
 {
-    p += (initialMass - fvc::domainIntegrate(thermo->psi()*p))
+    p += (initialMass - fvc::domainIntegrate(psi*p))
-        /fvc::domainIntegrate(thermo->psi());
+        /fvc::domainIntegrate(psi);
 }
 */

applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C

@@ -35,7 +35,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "basicThermo.H"
+#include "basicPsiThermo.H"
 #include "turbulenceModel.H"
 #include "bound.H"
 

applications/solvers/compressible/rhoPisoFoam/createFields.H

@@ -1,13 +1,14 @@
     Info<< "Reading thermophysical properties\n" << endl;
 
-    autoPtr<basicThermo> thermo
+    autoPtr<basicPsiThermo> pThermo
     (
-        basicThermo::New(mesh)
+        basicPsiThermo::New(mesh)
     );
+    basicPsiThermo& thermo = pThermo();
 
-    volScalarField& p = thermo->p();
+    volScalarField& p = thermo.p();
-    volScalarField& h = thermo->h();
+    volScalarField& h = thermo.h();
-    const volScalarField& psi = thermo->psi();
+    const volScalarField& psi = thermo.psi();
 
     volScalarField rho
     (
@@ -19,7 +20,7 @@
             IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
-        thermo->rho()
+        thermo.rho()
     );
 
     Info<< "\nReading field U\n" << endl;
@@ -47,7 +48,7 @@
             rho,
             U,
             phi,
-            thermo()
+            thermo
         )
     );
 

applications/solvers/compressible/rhoPisoFoam/hEqn.H

@@ -8,5 +8,5 @@
         DpDt
     );
 
-    thermo->correct();
+    thermo.correct();
 }

applications/solvers/compressible/rhoPisoFoam/pEqn.H

@@ -1,4 +1,4 @@
-rho = thermo->rho();
+rho = thermo.rho();
 
 volScalarField rUA = 1.0/UEqn.A();
 U = rUA*UEqn.H();
@@ -8,7 +8,7 @@ if (transonic)
     surfaceScalarField phid
     (
         "phid",
-        fvc::interpolate(thermo->psi())
+        fvc::interpolate(psi)
        *(
             (fvc::interpolate(U) & mesh.Sf())
           + fvc::ddtPhiCorr(rUA, rho, U, phi)

applications/solvers/compressible/rhoPisoFoam/rhoPisoFoam.C

@@ -31,7 +31,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "basicThermo.H"
+#include "basicPsiThermo.H"
 #include "turbulenceModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -77,7 +77,7 @@ int main(int argc, char *argv[])
 
         turbulence->correct();
 
-        rho = thermo->rho();
+        rho = thermo.rho();
 
         runTime.write();
 

applications/solvers/compressible/rhoPorousSimpleFoam/createFields.H

@@ -1,9 +1,10 @@
     Info<< "Reading thermophysical properties\n" << endl;
 
-    autoPtr<basicThermo> thermo
+    autoPtr<basicPsiThermo> pThermo
     (
-        basicThermo::New(mesh)
+        basicPsiThermo::New(mesh)
     );
+    basicPsiThermo& thermo = pThermo();
 
     volScalarField rho
     (
@@ -15,11 +16,12 @@
             IOobject::READ_IF_PRESENT,
             IOobject::AUTO_WRITE
         ),
-        thermo->rho()
+        thermo.rho()
     );
 
-    volScalarField& p = thermo->p();
+    volScalarField& p = thermo.p();
-    volScalarField& h = thermo->h();
+    volScalarField& h = thermo.h();
+    const volScalarField& psi = thermo.psi();
 
 
     Info<< "Reading field U\n" << endl;
@@ -36,7 +38,7 @@
         mesh
     );
 
-#   include "compressibleCreatePhi.H"
+    #include "compressibleCreatePhi.H"
 
 
     label pRefCell = 0;
@@ -56,7 +58,7 @@
             rho,
             U,
             phi,
-            thermo()
+            thermo
         )
     );
 

applications/solvers/compressible/rhoPorousSimpleFoam/hEqn.H

@@ -14,5 +14,5 @@
     eqnResidual = hEqn.solve().initialResidual();
     maxResidual = max(eqnResidual, maxResidual);
 
-    thermo->correct();
+    thermo.correct();
 }

applications/solvers/compressible/rhoPorousSimpleFoam/pEqn.H

@@ -65,10 +65,10 @@ bound(p, pMin);
 // to obey overall mass continuity
 if (closedVolume)
 {
-    p += (initialMass - fvc::domainIntegrate(thermo->psi()*p))
+    p += (initialMass - fvc::domainIntegrate(psi*p))
-        /fvc::domainIntegrate(thermo->psi());
+        /fvc::domainIntegrate(psi);
 }
 
-rho = thermo->rho();
+rho = thermo.rho();
 rho.relax();
 Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl;

applications/solvers/compressible/rhoPorousSimpleFoam/rhoPorousSimpleFoam.C

@@ -27,12 +27,12 @@ Application
 
 Description
     Steady-state solver for turbulent flow of compressible fluids with
-    implicit or explicit porosity treatment
+    RANS turbulence modelling, and implicit or explicit porosity treatment
 
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "basicThermo.H"
+#include "basicPsiThermo.H"
 #include "RASModel.H"
 #include "porousZones.H"
 

applications/solvers/compressible/rhoSimpleFoam/createFields.H

@@ -1,9 +1,10 @@
     Info<< "Reading thermophysical properties\n" << endl;
 
-    autoPtr<basicThermo> thermo
+    autoPtr<basicPsiThermo> pThermo
     (
-        basicThermo::New(mesh)
+        basicPsiThermo::New(mesh)
     );
+    basicPsiThermo& thermo = pThermo();
 
     volScalarField rho
     (
@@ -15,12 +16,12 @@
             IOobject::READ_IF_PRESENT,
             IOobject::AUTO_WRITE
         ),
-        thermo->rho()
+        thermo.rho()
     );
 
-    volScalarField& p = thermo->p();
+    volScalarField& p = thermo.p();
-    volScalarField& h = thermo->h();
+    volScalarField& h = thermo.h();
-
+    const volScalarField& psi = thermo.psi();
 
     Info<< "Reading field U\n" << endl;
     volVectorField U
@@ -56,7 +57,7 @@
             rho,
             U,
             phi,
-            thermo()
+            thermo
         )
     );
 

applications/solvers/compressible/rhoSimpleFoam/hEqn.H

@@ -14,5 +14,5 @@
     eqnResidual = hEqn.solve().initialResidual();
     maxResidual = max(eqnResidual, maxResidual);
 
-    thermo->correct();
+    thermo.correct();
 }

applications/solvers/compressible/rhoSimpleFoam/pEqn.H

@@ -1,4 +1,4 @@
-rho = thermo->rho();
+rho = thermo.rho();
 
 volScalarField rUA = 1.0/UEqn().A();
 U = rUA*UEqn().H();
@@ -11,7 +11,7 @@ if (transonic)
     surfaceScalarField phid
     (
         "phid",
-        fvc::interpolate(thermo->psi())*(fvc::interpolate(U) & mesh.Sf())
+        fvc::interpolate(psi)*(fvc::interpolate(U) & mesh.Sf())
     );
 
     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
@@ -82,7 +82,7 @@ else
 // Explicitly relax pressure for momentum corrector
 p.relax();
 
-rho = thermo->rho();
+rho = thermo.rho();
 rho.relax();
 Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl;
 
@@ -95,6 +95,6 @@ bound(p, pMin);
 // to obey overall mass continuity
 if (closedVolume)
 {
-    p += (initialMass - fvc::domainIntegrate(thermo->psi()*p))
+    p += (initialMass - fvc::domainIntegrate(psi*p))
-        /fvc::domainIntegrate(thermo->psi());
+        /fvc::domainIntegrate(psi);
 }

applications/solvers/compressible/rhoSimpleFoam/rhoSimpleFoam.C

@@ -26,13 +26,13 @@ Application
     rhoSimpleFoam
 
 Description
-    Steady-state SIMPLE solver for laminar or turbulent flow of
+    Steady-state SIMPLE solver for laminar or turbulent RANS flow of
     compressible fluids.
 
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "basicThermo.H"
+#include "basicPsiThermo.H"
 #include "RASModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

applications/solvers/compressible/sonicDyMFoam/compressibleContinuityErrs.H

@@ -1,5 +1,5 @@
 {
-#   include "rhoEqn.H"
+    #include "rhoEqn.H"
 }
 {
     scalar sumLocalContErr =

applications/solvers/compressible/sonicDyMFoam/createFields.H

@@ -1,13 +1,14 @@
     Info<< "Reading thermophysical properties\n" << endl;
 
-    autoPtr<basicThermo> thermo
+    autoPtr<basicPsiThermo> pThermo
     (
-        basicThermo::New(mesh)
+        basicPsiThermo::New(mesh)
     );
+    basicPsiThermo& thermo = pThermo();
 
-    volScalarField& p = thermo->p();
+    volScalarField& p = thermo.p();
-    volScalarField& h = thermo->h();
+    volScalarField& h = thermo.h();
-    const volScalarField& psi = thermo->psi();
+    const volScalarField& psi = thermo.psi();
 
     volScalarField rho
     (
@@ -17,7 +18,7 @@
             runTime.timeName(),
             mesh
         ),
-        thermo->rho()
+        thermo.rho()
     );
 
     Info<< "Reading field U\n" << endl;
@@ -45,7 +46,7 @@
             rho,
             U,
             phi,
-            thermo()
+            thermo
         )
     );
 

applications/solvers/compressible/sonicDyMFoam/sonicDyMFoam.C

@@ -32,7 +32,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "basicThermo.H"
+#include "basicPsiThermo.H"
 #include "turbulenceModel.H"
 #include "motionSolver.H"
 
@@ -84,8 +84,8 @@ int main(int argc, char *argv[])
             surfaceScalarField phid
             (
                 "phid",
-                fvc::interpolate(psi)*
+                fvc::interpolate(psi)
-                (
+               *(
                     (fvc::interpolate(U) & mesh.Sf()) - fvc::meshPhi(rho, U)
                 )
             );

applications/solvers/compressible/sonicFoam/createFields.H

@@ -1,13 +1,14 @@
     Info<< "Reading thermophysical properties\n" << endl;
 
-    autoPtr<basicThermo> thermo
+    autoPtr<basicPsiThermo> pThermo
     (
-        basicThermo::New(mesh)
+        basicPsiThermo::New(mesh)
     );
+    basicPsiThermo& thermo = pThermo();
 
-    volScalarField& p = thermo->p();
+    volScalarField& p = thermo.p();
-    volScalarField& h = thermo->h();
+    volScalarField& h = thermo.h();
-    const volScalarField& psi = thermo->psi();
+    const volScalarField& psi = thermo.psi();
 
     volScalarField rho
     (
@@ -17,7 +18,7 @@
             runTime.timeName(),
             mesh
         ),
-        thermo->rho()
+        thermo.rho()
     );
 
     Info<< "Reading field U\n" << endl;
@@ -45,7 +46,7 @@
             rho,
             U,
             phi,
-            thermo()
+            thermo
         )
     );
 

applications/solvers/compressible/sonicFoam/hEqn.H

@@ -8,5 +8,5 @@
             DpDt
         );
 
-        thermo->correct();
+        thermo.correct();
 }

applications/solvers/compressible/sonicFoam/sonicFoam.C

@@ -32,7 +32,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "basicThermo.H"
+#include "basicPsiThermo.H"
 #include "turbulenceModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -80,7 +80,7 @@ int main(int argc, char *argv[])
             surfaceScalarField phid
             (
                 "phid",
-                fvc::interpolate(thermo->psi())
+                fvc::interpolate(psi)
                *(
                    (fvc::interpolate(U) & mesh.Sf())
                  + fvc::ddtPhiCorr(rUA, rho, U, phi)

applications/solvers/compressible/sonicLiquidFoam/createFields.H

@@ -41,4 +41,4 @@
     );
 
 
-#   include "compressibleCreatePhi.H"
+    #include "compressibleCreatePhi.H"

applications/solvers/compressible/sonicLiquidFoam/sonicLiquidFoam.C

@@ -37,16 +37,15 @@ Description
 
 int main(int argc, char *argv[])
 {
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "readThermodynamicProperties.H"
+    #include "readTransportProperties.H"
+    #include "createFields.H"
+    #include "initContinuityErrs.H"
 
-#   include "setRootCase.H"
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-#   include "createTime.H"
-#   include "createMesh.H"
-#   include "readThermodynamicProperties.H"
-#   include "readTransportProperties.H"
-#   include "createFields.H"
-#   include "initContinuityErrs.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
 
@@ -54,10 +53,10 @@ int main(int argc, char *argv[])
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-#       include "readPISOControls.H"
+        #include "readPISOControls.H"
-#       include "compressibleCourantNo.H"
+        #include "compressibleCourantNo.H"
 
-#       include "rhoEqn.H"
+        #include "rhoEqn.H"
 
         fvVectorMatrix UEqn
         (
@@ -79,8 +78,8 @@ int main(int argc, char *argv[])
             surfaceScalarField phid
             (
                 "phid",
-                psi*
+                psi
-                (
+               *(
                     (fvc::interpolate(U) & mesh.Sf())
                   + fvc::ddtPhiCorr(rUA, rho, U, phi)
                 )
@@ -100,7 +99,7 @@ int main(int argc, char *argv[])
 
             phi += pEqn.flux();
 
-#           include "compressibleContinuityErrs.H"
+            #include "compressibleContinuityErrs.H"
 
             U -= rUA*fvc::grad(p);
             U.correctBoundaryConditions();

applications/solvers/heatTransfer/buoyantPisoFoam/buoyantPisoFoam.C

@@ -28,12 +28,12 @@ Application
 Description
     Transient Solver for buoyant, turbulent flow of compressible fluids for
     ventilation and heat-transfer. Turbulence is modelled using a run-time
-    selectable compressible RAS model.
+    selectable compressible RAS or LES model.
 
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "basicThermo.H"
+#include "basicRhoThermo.H"
 #include "turbulenceModel.H"
 #include "fixedGradientFvPatchFields.H"
 
@@ -81,6 +81,8 @@ int main(int argc, char *argv[])
 
         turbulence->correct();
 
+        rho = thermo.rho();
+
         runTime.write();
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"

applications/solvers/heatTransfer/buoyantPisoFoam/createFields.H

@@ -1,9 +1,10 @@
     Info<< "Reading thermophysical properties\n" << endl;
 
-    autoPtr<basicThermo> thermo
+    autoPtr<basicRhoThermo> pThermo
     (
-        basicThermo::New(mesh)
+        basicRhoThermo::New(mesh)
     );
+    basicRhoThermo& thermo = pThermo();
 
     volScalarField rho
     (
@@ -15,12 +16,12 @@
             IOobject::NO_READ,
             IOobject::NO_WRITE
         ),
-        thermo->rho()
+        thermo.rho()
     );
 
-    volScalarField& p = thermo->p();
+    volScalarField& p = thermo.p();
-    volScalarField& h = thermo->h();
+    volScalarField& h = thermo.h();
-    const volScalarField& psi = thermo->psi();
+    const volScalarField& psi = thermo.psi();
 
 
     Info<< "Reading field U\n" << endl;
@@ -37,7 +38,7 @@
         mesh
     );
 
-#   include "compressibleCreatePhi.H"
+    #include "compressibleCreatePhi.H"
 
 
     Info<< "Creating turbulence model\n" << endl;
@@ -48,7 +49,7 @@
             rho,
             U,
             phi,
-            thermo()
+            thermo
         )
     );
 
@@ -59,6 +60,9 @@
         fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p)
     );
 
-    thermo->correct();
+    thermo.correct();
 
     dimensionedScalar initialMass = fvc::domainIntegrate(rho);
+
+    dimensionedScalar totalVolume = sum(mesh.V());
+

applications/solvers/heatTransfer/buoyantPisoFoam/hEqn.H

@@ -11,5 +11,5 @@
     hEqn.relax();
     hEqn.solve();
 
-    thermo->correct();
+    thermo.correct();
 }

applications/solvers/heatTransfer/buoyantPisoFoam/pEqn.H

@@ -1,7 +1,11 @@
 {
     bool closedVolume = p.needReference();
 
-    rho = thermo->rho();
+    rho = thermo.rho();
+
+    // Thermodynamic density needs to be updated by psi*d(p) after the
+    // pressure solution - done in 2 parts. Part 1:
+    thermo.rho() -= psi*p;
 
     volScalarField rUA = 1.0/UEqn.A();
     surfaceScalarField rhorUAf("(rho*(1|A(U)))", fvc::interpolate(rho*rUA));
@@ -23,7 +27,7 @@
     {
         fvScalarMatrix pEqn
         (
-            fvm::ddt(psi, p)
+            fvc::ddt(rho) + psi*correction(fvm::ddt(p))
           + fvc::div(phi)
           - fvm::laplacian(rhorUAf, p)
         );
@@ -43,6 +47,9 @@
         }
     }
 
+    // Second part of thermodynamic density update
+    thermo.rho() += psi*p;
+
     U += rUA*fvc::reconstruct((phi - phiU)/rhorUAf);
     U.correctBoundaryConditions();
 
@@ -55,8 +62,10 @@
     // to obey overall mass continuity
     if (closedVolume)
     {
-        p += (initialMass - fvc::domainIntegrate(thermo->psi()*p))
+        p +=
-            /fvc::domainIntegrate(thermo->psi());
+            (initialMass - fvc::domainIntegrate(psi*p))
-        rho = thermo->rho();
+            /fvc::domainIntegrate(psi);
+        thermo.rho() = psi*p;
+        rho += (initialMass - fvc::domainIntegrate(rho))/totalVolume;
     }
 }

applications/solvers/heatTransfer/buoyantSimpleFoam/buoyantSimpleFoam.C

@@ -31,7 +31,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "basicThermo.H"
+#include "basicPsiThermo.H"
 #include "RASModel.H"
 #include "fixedGradientFvPatchFields.H"
 

applications/solvers/heatTransfer/buoyantSimpleFoam/createFields.H

@@ -1,9 +1,10 @@
     Info<< "Reading thermophysical properties\n" << endl;
 
-    autoPtr<basicThermo> thermo
+    autoPtr<basicPsiThermo> pThermo
     (
-        basicThermo::New(mesh)
+        basicPsiThermo::New(mesh)
     );
+    basicPsiThermo& thermo = pThermo();
 
     volScalarField rho
     (
@@ -15,11 +16,12 @@
             IOobject::NO_READ,
             IOobject::NO_WRITE
         ),
-        thermo->rho()
+        thermo.rho()
     );
 
-    volScalarField& p = thermo->p();
+    volScalarField& p = thermo.p();
-    volScalarField& h = thermo->h();
+    volScalarField& h = thermo.h();
+    const volScalarField& psi = thermo.psi();
 
 
     Info<< "Reading field U\n" << endl;
@@ -36,7 +38,7 @@
         mesh
     );
 
-#   include "compressibleCreatePhi.H"
+    #include "compressibleCreatePhi.H"
 
 
     Info<< "Creating turbulence model\n" << endl;
@@ -47,11 +49,11 @@
             rho,
             U,
             phi,
-            thermo()
+            thermo
         )
     );
 
-    thermo->correct();
+    thermo.correct();
 
     label pRefCell = 0;
     scalar pRefValue = 0.0;

applications/solvers/heatTransfer/buoyantSimpleFoam/hEqn.H

@@ -14,5 +14,5 @@
     eqnResidual = hEqn.solve().initialResidual();
     maxResidual = max(eqnResidual, maxResidual);
 
-    thermo->correct();
+    thermo.correct();
 }

applications/solvers/heatTransfer/buoyantSimpleFoam/initConvergenceCheck.H

@@ -1,4 +1,4 @@
-// initialize values for convergence checks
+    // initialize values for convergence checks
 
     scalar eqnResidual = 1, maxResidual = 0;
     scalar convergenceCriterion = 0;

applications/solvers/heatTransfer/buoyantSimpleFoam/pEqn.H

@@ -1,5 +1,5 @@
 {
-    rho = thermo->rho();
+    rho = thermo.rho();
 
     volScalarField rUA = 1.0/UEqn().A();
     surfaceScalarField rhorUAf("(rho*(1|A(U)))", fvc::interpolate(rho*rUA));
@@ -39,8 +39,8 @@
             // to obey overall mass continuity
             if (closedVolume)
             {
-                p += (initialMass - fvc::domainIntegrate(thermo->psi()*p))
+                p += (initialMass - fvc::domainIntegrate(psi*p))
-                    /fvc::domainIntegrate(thermo->psi());
+                    /fvc::domainIntegrate(psi);
             }
 
             // Calculate the conservative fluxes
@@ -58,7 +58,7 @@
 
     #include "continuityErrs.H"
 
-    rho = thermo->rho();
+    rho = thermo.rho();
     rho.relax();
     Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value()
         << endl;

applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/buoyantSimpleRadiationFoam.C

@@ -32,7 +32,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "basicThermo.H"
+#include "basicPsiThermo.H"
 #include "RASModel.H"
 #include "fixedGradientFvPatchFields.H"
 #include "radiationModel.H"
@@ -61,7 +61,7 @@ int main(int argc, char *argv[])
 #       include "readSIMPLEControls.H"
 #       include "initConvergenceCheck.H"
 
-        pd.storePrevIter();
+        p.storePrevIter();
         rho.storePrevIter();
 
         // Pressure-velocity SIMPLE corrector

applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/hEqn.H

@@ -7,7 +7,7 @@
      ==
         fvc::div(phi/fvc::interpolate(rho)*fvc::interpolate(p))
       - p*fvc::div(phi/fvc::interpolate(rho))
-      + radiation->Sh(thermo())
+      + radiation->Sh(thermo)
     );
 
     hEqn.relax();
@@ -15,7 +15,7 @@
     eqnResidual = hEqn.solve().initialResidual();
     maxResidual = max(eqnResidual, maxResidual);
 
-    thermo->correct();
+    thermo.correct();
 
     radiation->correct();
 }

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C

@@ -32,7 +32,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "basicThermo.H"
+#include "basicPsiThermo.H"
 #include "turbulenceModel.H"
 #include "fixedGradientFvPatchFields.H"
 #include "regionProperties.H"
@@ -42,37 +42,36 @@ Description
 
 int main(int argc, char *argv[])
 {
-
+    #include "setRootCase.H"
-#   include "setRootCase.H"
+    #include "createTime.H"
-#   include "createTime.H"
 
     regionProperties rp(runTime);
 
-#   include "createFluidMeshes.H"
+    #include "createFluidMeshes.H"
-#   include "createSolidMeshes.H"
+    #include "createSolidMeshes.H"
 
-#   include "createFluidFields.H"
+    #include "createFluidFields.H"
-#   include "createSolidFields.H"
+    #include "createSolidFields.H"
 
-#   include "initContinuityErrs.H"
+    #include "initContinuityErrs.H"
 
-#   include "readTimeControls.H"
+    #include "readTimeControls.H"
 
     if (fluidRegions.size())
     {
-#       include "compressibleMultiRegionCourantNo.H"
+        #include "compressibleMultiRegionCourantNo.H"
-#       include "setInitialDeltaT.H"
+        #include "setInitialDeltaT.H"
     }
 
     while (runTime.run())
     {
-#       include "readTimeControls.H"
+        #include "readTimeControls.H"
-#       include "readPIMPLEControls.H"
+        #include "readPIMPLEControls.H"
 
         if (fluidRegions.size())
         {
-#           include "compressibleMultiRegionCourantNo.H"
+            #include "compressibleMultiRegionCourantNo.H"
-#           include "setDeltaT.H"
+            #include "setDeltaT.H"
         }
 
         runTime++;
@@ -83,8 +82,8 @@ int main(int argc, char *argv[])
         {
             forAll(fluidRegions, i)
             {
-#               include "setRegionFluidFields.H"
+                #include "setRegionFluidFields.H"
-#               include "storeOldFluidFields.H"
+                #include "storeOldFluidFields.H"
             }
         }
 
@@ -96,18 +95,18 @@ int main(int argc, char *argv[])
             {
                 Info<< "\nSolving for fluid region "
                     << fluidRegions[i].name() << endl;
-#               include "setRegionFluidFields.H"
+                #include "setRegionFluidFields.H"
-#               include "readFluidMultiRegionPIMPLEControls.H"
+                #include "readFluidMultiRegionPIMPLEControls.H"
-#               include "solveFluid.H"
+                #include "solveFluid.H"
             }
 
             forAll(solidRegions, i)
             {
                 Info<< "\nSolving for solid region "
                     << solidRegions[i].name() << endl;
-#               include "setRegionSolidFields.H"
+                #include "setRegionSolidFields.H"
-#               include "readSolidMultiRegionPIMPLEControls.H"
+                #include "readSolidMultiRegionPIMPLEControls.H"
-#               include "solveSolid.H"
+                #include "solveSolid.H"
             }
         }
 

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H

@@ -1,5 +1,5 @@
     // Initialise fluid field pointer lists
-    PtrList<basicThermo> thermoFluid(fluidRegions.size());
+    PtrList<basicPsiThermo> thermoFluid(fluidRegions.size());
     PtrList<volScalarField> rhoFluid(fluidRegions.size());
     PtrList<volScalarField> KFluid(fluidRegions.size());
     PtrList<volVectorField> UFluid(fluidRegions.size());
@@ -19,7 +19,7 @@
         thermoFluid.set
         (
             i,
-            basicThermo::New(fluidRegions[i]).ptr()
+            basicPsiThermo::New(fluidRegions[i]).ptr()
         );
 
         Info<< "    Adding to rhoFluid\n" << endl;

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H

@@ -1,6 +1,6 @@
     const fvMesh& mesh = fluidRegions[i];
 
-    basicThermo& thermo = thermoFluid[i];
+    basicPsiThermo& thermo = thermoFluid[i];
     volScalarField& rho = rhoFluid[i];
     volScalarField& K = KFluid[i];
     volVectorField& U = UFluid[i];

applications/solvers/lagrangian/coalChemistryFoam/Make/options

@@ -13,7 +13,7 @@ EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/solids/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/solidMixture/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/thermophysicalFunctions/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/combustion/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
     -I$(LIB_SRC)/ODE/lnInclude
@@ -33,7 +33,7 @@ EXE_LIBS = \
     -lsolids \
     -lsolidMixture \
     -lthermophysicalFunctions \
-    -lcombustionThermophysicalModels \
+    -lreactionThermophysicalModels \
     -lchemistryModel \
     -lradiation \
     -lODE

applications/solvers/lagrangian/coalChemistryFoam/coalChemistryFoam.C

@@ -36,9 +36,9 @@ Description
 #include "turbulenceModel.H"
 #include "basicThermoCloud.H"
 #include "CoalCloud.H"
-#include "chemistryModel.H"
+#include "psiChemistryModel.H"
 #include "chemistrySolver.H"
-#include "reactingThermoTypes.H"
+#include "thermoPhysicsTypes.H"
 #include "timeActivatedExplicitCellSource.H"
 #include "radiationModel.H"
 
@@ -113,12 +113,9 @@ int main(int argc, char *argv[])
 
         turbulence->correct();
 
-        rho = thermo->rho();
+        rho = thermo.rho();
 
-        if (runTime.write())
+        runTime.write();
-        {
-            #include "additionalOutput.H"
-        }
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
             << "  ClockTime = " << runTime.elapsedClockTime() << " s"

applications/solvers/lagrangian/coalChemistryFoam/createClouds.H

@@ -1,11 +1,11 @@
 Info<< "\nConstructing coal cloud" << endl;
-CoalCloud<specieReactingProperties> coalParcels
+CoalCloud<gasThermoPhysics> coalParcels
 (
     "coalCloud1",
     rho,
     U,
     g,
-    thermo()
+    thermo
 );
 
 Info<< "\nConstructing limestone cloud" << endl;
@@ -15,5 +15,5 @@ basicThermoCloud limestoneParcels
     rho,
     U,
     g,
-    thermo()
+    thermo
 );

applications/solvers/lagrangian/coalChemistryFoam/createFields.H

@@ -1,19 +1,30 @@
     Info<< "Reading thermophysical properties\n" << endl;
 
-    autoPtr<hCombustionThermo> thermo
+    autoPtr<psiChemistryModel> pChemistry
     (
-        hCombustionThermo::New(mesh)
+        psiChemistryModel::New(mesh)
     );
+    psiChemistryModel& chemistry = pChemistry();
 
-    combustionMixture& composition = thermo->composition();
+    hCombustionThermo& thermo = chemistry.thermo();
+
+    basicMultiComponentMixture& composition = thermo.composition();
     PtrList<volScalarField>& Y = composition.Y();
 
-    word inertSpecie(thermo->lookup("inertSpecie"));
+    word inertSpecie(thermo.lookup("inertSpecie"));
 
-    volScalarField& p = thermo->p();
+    volScalarField& p = thermo.p();
-    volScalarField& h = thermo->h();
+    volScalarField& h = thermo.h();
-    const volScalarField& T = thermo->T();
+    const volScalarField& T = thermo.T();
-    const volScalarField& psi = thermo->psi();
+    const volScalarField& psi = thermo.psi();
+
+    multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
+
+    forAll (Y, i)
+    {
+        fields.add(Y[i]);
+    }
+    fields.add(h);
 
     volScalarField rho
     (
@@ -25,7 +36,7 @@
             IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
-        thermo->rho()
+        thermo.rho()
     );
 
     // lagrangian effective density field - used externally (optional)
@@ -73,7 +84,7 @@
         mesh
     );
 
-#   include "compressibleCreatePhi.H"
+    #include "compressibleCreatePhi.H"
 
     DimensionedField<scalar, volMesh> kappa
     (
@@ -97,13 +108,16 @@
             rho,
             U,
             phi,
-            thermo()
+            thermo
         )
     );
 
     Info<< "Creating field DpDt\n" << endl;
-    volScalarField DpDt =
+    volScalarField DpDt
-        fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
+    (
+        "DpDt",
+        fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p)
+    );
 
     Info<< "\nConstructing explicit enthalpy cell source" << endl;
     timeActivatedExplicitCellSource enthalpySource
@@ -112,18 +126,3 @@
         mesh,
         h.dimensions()*phi.dimensions()/mesh.V().dimensions()
     );
-
-    Info << "Constructing chemical mechanism" << endl;
-    chemistryModel chemistry
-    (
-        thermo(),
-        rho
-    );
-
-    multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
-
-    forAll (Y, i)
-    {
-        fields.add(Y[i]);
-    }
-    fields.add(h);

applications/solvers/lagrangian/coalChemistryFoam/hEqn.H

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

applications/solvers/lagrangian/coalChemistryFoam/pEqn.H

@@ -1,4 +1,4 @@
-rho = thermo->rho();
+rho = thermo.rho();
 
 volScalarField rUA = 1.0/UEqn.A();
 U = rUA*UEqn.H();
@@ -8,7 +8,7 @@ if (transonic)
     surfaceScalarField phid
     (
         "phid",
-        fvc::interpolate(thermo->psi())
+        fvc::interpolate(psi)
        *(
             (fvc::interpolate(U) & mesh.Sf())
           + fvc::ddtPhiCorr(rUA, rho, U, phi)

applications/solvers/lagrangian/coalChemistryFoam/readChemistryProperties.H

@@ -8,7 +8,8 @@ IOdictionary chemistryProperties
         runTime.constant(),
         mesh,
         IOobject::MUST_READ,
-        IOobject::NO_WRITE
+        IOobject::NO_WRITE,
+        false
     )
 );
 

applications/solvers/lagrangian/reactingParcelFoam/Make/options

@@ -12,7 +12,7 @@ EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/solids/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/solidMixture/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/thermophysicalFunctions/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/combustion/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
     -I$(LIB_SRC)/ODE/lnInclude
@@ -31,7 +31,7 @@ EXE_LIBS = \
     -lsolids \
     -lsolidMixture \
     -lthermophysicalFunctions \
-    -lcombustionThermophysicalModels \
+    -lreactionThermophysicalModels \
     -lchemistryModel \
     -lradiation \
     -lODE

applications/solvers/lagrangian/reactingParcelFoam/YEqn.H

@@ -26,7 +26,8 @@ tmp<fv::convectionScheme<scalar> > mvConvection
               - fvm::laplacian(turbulence->muEff(), Yi)
               ==
                 parcels.Srho(i)
-              + kappa*chemistry.RR(i)().dimensionedInternalField()
+              + kappa*chemistry.RR(i)().dimensionedInternalField(),
+                mesh.solver("Yi")
             );
 
             Yi.max(0.0);

applications/solvers/lagrangian/reactingParcelFoam/additionalOutput.H

@@ -1,48 +0,0 @@
-{
-    tmp<volScalarField> tdQ
-    (
-        new volScalarField
-        (
-            IOobject
-            (
-                "dQ",
-                runTime.timeName(),
-                mesh,
-                IOobject::NO_READ,
-                IOobject::AUTO_WRITE
-            ),
-            mesh,
-            dimensionedScalar
-            (
-                "zero",
-                dimensionSet(1, -3, -1, 0, 0, 0, 0),
-                0.0
-            )
-        )
-    );
-
-    scalarField& dQ = tdQ();
-
-    scalarField cp(dQ.size(), 0.0);
-
-    forAll(Y, i)
-    {
-        volScalarField RRi = chemistry.RR(i);
-
-        forAll(h, celli)
-        {
-            scalar Ti = T[celli];
-            cp[celli] += Y[i][celli]*chemistry.specieThermo()[i].Cp(Ti);
-            scalar hi = chemistry.specieThermo()[i].h(Ti);
-            scalar RR = RRi[celli];
-            dQ[celli] -= hi*RR;
-        }
-    }
-
-    forAll(dQ, celli)
-    {
-        dQ[celli] /= cp[celli];
-    }
-
-    tdQ().write();
-}

applications/solvers/lagrangian/reactingParcelFoam/createClouds.H

@@ -1,9 +1,9 @@
 Info<< "\nConstructing reacting cloud" << endl;
-BasicReactingCloud<specieReactingProperties> parcels
+BasicReactingCloud<gasThermoPhysics> parcels
 (
     "reactingCloud1",
     rho,
     U,
     g,
-    thermo()
+    thermo
 );

applications/solvers/lagrangian/reactingParcelFoam/createFields.H

@@ -1,19 +1,22 @@
     Info<< "Reading thermophysical properties\n" << endl;
 
-    autoPtr<hCombustionThermo> thermo
+    autoPtr<psiChemistryModel> pChemistry
     (
-        hCombustionThermo::New(mesh)
+        psiChemistryModel::New(mesh)
     );
+    psiChemistryModel& chemistry = pChemistry();
 
-    combustionMixture& composition = thermo->composition();
+    hCombustionThermo& thermo = chemistry.thermo();
+
+    basicMultiComponentMixture& composition = thermo.composition();
     PtrList<volScalarField>& Y = composition.Y();
 
-    word inertSpecie(thermo->lookup("inertSpecie"));
+    word inertSpecie(thermo.lookup("inertSpecie"));
 
-    volScalarField& p = thermo->p();
+    volScalarField& p = thermo.p();
-    volScalarField& h = thermo->h();
+    volScalarField& h = thermo.h();
-    const volScalarField& T = thermo->T();
+    const volScalarField& T = thermo.T();
-    const volScalarField& psi = thermo->psi();
+    const volScalarField& psi = thermo.psi();
 
     volScalarField rho
     (
@@ -25,7 +28,7 @@
             IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
-        thermo->rho()
+        thermo.rho()
     );
 
     Info<< "\nReading field U\n" << endl;
@@ -42,7 +45,7 @@
         mesh
     );
 
-#   include "compressibleCreatePhi.H"
+    #include "compressibleCreatePhi.H"
 
     DimensionedField<scalar, volMesh> kappa
     (
@@ -66,19 +69,15 @@
             rho,
             U,
             phi,
-            thermo()
+            thermo
         )
     );
 
     Info<< "Creating field DpDt\n" << endl;
-    volScalarField DpDt =
+    volScalarField DpDt
-        fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
-
-    Info << "Constructing chemical mechanism" << endl;
-    chemistryModel chemistry
     (
-        thermo(),
+        "DpDt",
-        rho
+        fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p)
     );
 
     multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;