updates for new thermo and general code tidying

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,23 +31,25 @@ Description
 \*---------------------------------------------------------------------------*/
 
 {
-scalar meanStCoNum = 0.0;
+    scalar meanStCoNum = 0.0;
 
-if (mesh.nInternalFaces())
+    if (mesh.nInternalFaces())
-{
+    {
-    surfaceScalarField SfUfbyDelta = 
+        surfaceScalarField SfUfbyDelta =
-        mesh.surfaceInterpolation::deltaCoeffs()
+            mesh.surfaceInterpolation::deltaCoeffs()
-       *mag(phiSt/fvc::interpolate(rho));
+        *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;
+        << " 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,8 +196,7 @@ public:
 
 
     // Destructor
-
+    ~SCOPE();
-        ~SCOPE();
 
 
     // Member functions

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/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/createFields.H

@@ -84,7 +84,7 @@
         mesh
     );
 
-#   include "compressibleCreatePhi.H"
+    #include "compressibleCreatePhi.H"
 
     DimensionedField<scalar, volMesh> kappa
     (

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/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/dieselEngineFoam/spraySummary.H

@@ -1,30 +1,30 @@
-        label Nparcels = dieselSpray.size();
+    label Nparcels = dieselSpray.size();
-        reduce(Nparcels, sumOp<label>());
+    reduce(Nparcels, sumOp<label>());
 
-        Info<< "\nNumber of parcels in system.... | "
+    Info<< "\nNumber of parcels in system.... | "
-            << Nparcels << endl
+        << Nparcels << endl
-            << "Injected liquid mass........... | "
+        << "Injected liquid mass........... | "
-            << 1e6*dieselSpray.injectedMass(runTime.value()) << " mg" << endl
+        << 1e6*dieselSpray.injectedMass(runTime.value()) << " mg" << endl
-            << "Liquid Mass in system.......... | "
+        << "Liquid Mass in system.......... | "
-            << 1e6*dieselSpray.liquidMass() << " mg" << endl
+        << 1e6*dieselSpray.liquidMass() << " mg" << endl
-            << "SMD, Dmax...................... | "
+        << "SMD, Dmax...................... | "
-            << dieselSpray.smd()*1e6 << " mu, "
+        << dieselSpray.smd()*1e6 << " mu, "
-            << dieselSpray.maxD()*1e6 << " mu"
+        << dieselSpray.maxD()*1e6 << " mu"
-            << endl;
+        << endl;
 
-        scalar evapMass =
+    scalar evapMass =
-            dieselSpray.injectedMass(runTime.value())
+        dieselSpray.injectedMass(runTime.value())
-          - dieselSpray.liquidMass();
+    - dieselSpray.liquidMass();
 
-        scalar gasMass = fvc::domainIntegrate(rho).value();
+    scalar gasMass = fvc::domainIntegrate(rho).value();
 
-        if (dieselSpray.twoD())
+    if (dieselSpray.twoD())
-        {
+    {
-            gasMass *= 2.0*mathematicalConstant::pi/dieselSpray.angleOfWedge();
+        gasMass *= 2.0*mathematicalConstant::pi/dieselSpray.angleOfWedge();
-        }
+    }
 
-        scalar addedMass = gasMass - gasMass0;
+    scalar addedMass = gasMass - gasMass0;
 
-        Info<< "Added gas mass................. | " << 1e6*addedMass << " mg"
+    Info<< "Added gas mass................. | " << 1e6*addedMass << " mg"
-            << nl << "Evaporation Continuity Error... | "
+        << nl << "Evaporation Continuity Error... | "
-            << 1e6*(addedMass - evapMass) << " mg" << endl;
+        << 1e6*(addedMass - evapMass) << " mg" << endl;

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/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,8 +8,8 @@ if (transonic)
     surfaceScalarField phid
     (
         "phid",
-        fvc::interpolate(thermo->psi())
+        fvc::interpolate(psi)
-        *((fvc::interpolate(U) & mesh.Sf()) - fvc::meshPhi(rho, U))
+       *((fvc::interpolate(U) & mesh.Sf()) - fvc::meshPhi(rho, U))
     );
 
     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)

applications/solvers/combustion/engineFoam/readEngineTimeControls.H

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

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

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

applications/solvers/heatTransfer/buoyantPisoFoam/createFields.H

@@ -1,9 +1,10 @@
     Info<< "Reading thermophysical properties\n" << endl;
 
-    autoPtr<basicRhoThermo> thermo
+    autoPtr<basicRhoThermo> pThermo
     (
         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,7 +60,7 @@
         fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p)
     );
 
-    thermo->correct();
+    thermo.correct();
 
     dimensionedScalar initialMass = fvc::domainIntegrate(rho);
 

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,11 +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;
+    thermo.rho() -= psi*p;
 
     volScalarField rUA = 1.0/UEqn.A();
     surfaceScalarField rhorUAf("(rho*(1|A(U)))", fvc::interpolate(rho*rUA));
@@ -48,7 +48,7 @@
     }
 
     // Second part of thermodynamic density update
-    thermo->rho() += psi*p;
+    thermo.rho() += psi*p;
 
     U += rUA*fvc::reconstruct((phi - phiU)/rhorUAf);
     U.correctBoundaryConditions();
@@ -65,7 +65,7 @@
         p +=
             (initialMass - fvc::domainIntegrate(psi*p))
             /fvc::domainIntegrate(psi);
-        thermo->rho() = psi*p;
+        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/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/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,7 +69,7 @@
             rho,
             U,
             phi,
-            thermo()
+            thermo
         )
     );
 
@@ -74,13 +77,6 @@
     volScalarField DpDt =
         fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
 
-    Info << "Constructing chemical mechanism" << endl;
-    chemistryModel chemistry
-    (
-        thermo(),
-        rho
-    );
-
     multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
 
     forAll (Y, i)

applications/solvers/lagrangian/reactingParcelFoam/hEqn.H

@@ -7,14 +7,14 @@
      ==
         DpDt
      +  parcels.Sh()
-     +  radiation->Sh(thermo())
+     +  radiation->Sh(thermo)
     );
 
     hEqn.relax();
 
     hEqn.solve();
 
-    thermo->correct();
+    thermo.correct();
 
     radiation->correct();
 }

applications/solvers/lagrangian/reactingParcelFoam/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/lagrangian/reactingParcelFoam/reactingParcelFoam.C

@@ -35,9 +35,9 @@ Description
 #include "hCombustionThermo.H"
 #include "turbulenceModel.H"
 #include "BasicReactingCloud.H"
-#include "chemistryModel.H"
+#include "psiChemistryModel.H"
 #include "chemistrySolver.H"
-#include "thermoPhsyicsTypes.H"
+#include "thermoPhysicsTypes.H"
 #include "radiationModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -100,12 +100,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/reactingParcelFoam/rhoEqn.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2007 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

applications/solvers/lagrangian/trackedReactingParcelFoam/Make/options

@@ -1,4 +1,7 @@
+REACTINGPARCELFOAM = ../reactingParcelFoam
+
 EXE_INC = \
+    -I$(REACTINGPARCELFOAM) \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I${LIB_SRC}/meshTools/lnInclude \
     -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \

applications/solvers/lagrangian/trackedReactingParcelFoam/UEqn.H

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

applications/solvers/lagrangian/trackedReactingParcelFoam/YEqn.H

@@ -26,7 +26,7 @@ tmp<fv::convectionScheme<scalar> > mvConvection
               + mvConvection->fvmDiv(phi, Yi)
               - fvm::laplacian(turbulence->muEff(), Yi)
               ==
-                reactingParcels.Srho(i)
+                parcels.Srho(i)
               + kappa*chemistry.RR(i)().dimensionedInternalField()
               + pointMassSources.Su(i)
             );

applications/solvers/lagrangian/trackedReactingParcelFoam/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/trackedReactingParcelFoam/chemistry.H

@@ -1,25 +0,0 @@
-{
-    Info << "Solving chemistry" << endl;
-
-    chemistry.solve
-    (
-        runTime.value() - runTime.deltaT().value(),
-        runTime.deltaT().value()
-    );
-
-    // turbulent time scale
-    if (turbulentReaction)
-    {
-        DimensionedField<scalar, volMesh> tk =
-            Cmix*sqrt(turbulence->muEff()/rho/turbulence->epsilon());
-        DimensionedField<scalar, volMesh> tc =
-            chemistry.tc()().dimensionedInternalField();
-
-        // Chalmers PaSR model
-        kappa = (runTime.deltaT() + tc)/(runTime.deltaT() + tc + tk);
-    }
-    else
-    {
-        kappa = 1.0;
-    }
-}

applications/solvers/lagrangian/trackedReactingParcelFoam/createClouds.H

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

applications/solvers/lagrangian/trackedReactingParcelFoam/createFields.H

@@ -1,111 +0,0 @@
-    Info<< "Reading thermophysical properties" << nl << endl;
-
-    autoPtr<hCombustionThermo> thermo
-    (
-        hCombustionThermo::New(mesh)
-    );
-
-    combustionMixture& composition = thermo->composition();
-    PtrList<volScalarField>& Y = composition.Y();
-
-    word inertSpecie(thermo->lookup("inertSpecie"));
-
-    volScalarField& p = thermo->p();
-    volScalarField& h = thermo->h();
-    const volScalarField& T = thermo->T();
-    const volScalarField& psi = thermo->psi();
-
-    volScalarField rho
-    (
-        IOobject
-        (
-            "rho",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        thermo->rho()
-    );
-
-    Info<< "Reading field U" << nl << endl;
-    volVectorField U
-    (
-        IOobject
-        (
-            "U",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-#   include "compressibleCreatePhi.H"
-
-    DimensionedField<scalar, volMesh> kappa
-    (
-        IOobject
-        (
-            "kappa",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh,
-        dimensionedScalar("zero", dimless, 0.0)
-    );
-
-    Info<< "Creating turbulence model" << nl << endl;
-    autoPtr<compressible::turbulenceModel> turbulence
-    (
-        compressible::turbulenceModel::New
-        (
-            rho,
-            U,
-            phi,
-            thermo()
-        )
-    );
-
-    Info<< "Creating field DpDt" << nl << endl;
-    volScalarField DpDt =
-        fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
-
-    Info << "Constructing chemical mechanism" << nl << endl;
-    chemistryModel chemistry
-    (
-        thermo(),
-        rho
-    );
-
-    multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
-
-    forAll (Y, i)
-    {
-        fields.add(Y[i]);
-    }
-    fields.add(h);
-
-    Info<< "Creating porous zones" << nl << endl;
-    porousZones pZones(mesh);
-    Switch pressureImplicitPorosity(false);
-
-    label nUCorr = 0;
-    if (pZones.size())
-    {
-        // nUCorrectors for pressureImplicitPorosity
-        if (mesh.solutionDict().subDict("PISO").found("nUCorrectors"))
-        {
-            mesh.solutionDict().subDict("PISO").lookup("nUCorrectors")
-                >> nUCorr;
-        }
-
-        if (nUCorr > 0)
-        {
-            pressureImplicitPorosity = true;
-        }
-    }
-

applications/solvers/lagrangian/trackedReactingParcelFoam/createPorousZones.H

@@ -0,0 +1,21 @@
+    Info<< "Creating porous zones" << nl << endl;
+
+    porousZones pZones(mesh);
+    Switch pressureImplicitPorosity(false);
+
+    label nUCorr = 0;
+    if (pZones.size())
+    {
+        // nUCorrectors for pressureImplicitPorosity
+        if (mesh.solutionDict().subDict("PISO").found("nUCorrectors"))
+        {
+            mesh.solutionDict().subDict("PISO").lookup("nUCorrectors")
+                >> nUCorr;
+        }
+
+        if (nUCorr > 0)
+        {
+            pressureImplicitPorosity = true;
+        }
+    }
+

applications/solvers/lagrangian/trackedReactingParcelFoam/hEqn.H

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

applications/solvers/lagrangian/trackedReactingParcelFoam/pEqn.H

@@ -1,4 +1,4 @@
-rho = thermo->rho();
+rho = thermo.rho();
 
 if (pressureImplicitPorosity)
 {
@@ -14,7 +14,7 @@ if (transonic)
     surfaceScalarField phid
     (
         "phid",
-        fvc::interpolate(thermo->psi())
+        fvc::interpolate(psi)
        *(
             (fvc::interpolate(U) & mesh.Sf())
 //          + fvc::ddtPhiCorr(rUA, rho, U, phi)
@@ -40,7 +40,7 @@ if (transonic)
           + fvm::div(phid, p)
           - lapTerm()
          ==
-            reactingParcels.Srho()
+            parcels.Srho()
           + pointMassSources.Su()
         );
 
@@ -80,7 +80,7 @@ else
           + fvc::div(phi)
           - lapTerm()
          ==
-            reactingParcels.Srho()
+            parcels.Srho()
           + pointMassSources.Su()
         );
 

applications/solvers/lagrangian/trackedReactingParcelFoam/readChemistryProperties.H

@@ -1,22 +0,0 @@
-Info<< "Reading chemistry properties\n" << endl;
-
-IOdictionary chemistryProperties
-(
-    IOobject
-    (
-        "chemistryProperties",
-        runTime.constant(),
-        mesh,
-        IOobject::MUST_READ,
-        IOobject::NO_WRITE
-    )
-);
-
-Switch turbulentReaction(chemistryProperties.lookup("turbulentReaction"));
-
-dimensionedScalar Cmix("Cmix", dimless, 1.0);
-
-if (turbulentReaction)
-{
-    chemistryProperties.lookup("Cmix") >> Cmix;
-}

applications/solvers/lagrangian/trackedReactingParcelFoam/rhoEqn.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2007 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -36,7 +36,7 @@ Description
         fvm::ddt(rho)
       + fvc::div(phi)
       ==
-        reactingParcels.Srho()
+        parcels.Srho()
       + pointMassSources.Su()
     );
 }

applications/solvers/lagrangian/trackedReactingParcelFoam/trackedReactingParcelFoam.C

@@ -32,9 +32,9 @@ Description
 #include "hCombustionThermo.H"
 #include "turbulenceModel.H"
 #include "BasicTrackedReactingCloud.H"
-#include "chemistryModel.H"
+#include "psiChemistryModel.H"
 #include "chemistrySolver.H"
-#include "reactingThermoTypes.H"
+#include "thermoPhysicsTypes.H"
 #include "radiationModel.H"
 #include "porousZones.H"
 #include "timeActivatedExplicitMulticomponentPointSource.H"
@@ -53,6 +53,7 @@ int main(int argc, char *argv[])
     #include "createRadiationModel.H"
     #include "createClouds.H"
     #include "createMulticomponentPointSources.H"
+    #include "createPorousZones.H"
     #include "readPISOControls.H"
     #include "initContinuityErrs.H"
     #include "readTimeControls.H"
@@ -74,9 +75,9 @@ int main(int argc, char *argv[])
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        reactingParcels.evolve();
+        parcels.evolve();
 
-        reactingParcels.info();
+        parcels.info();
 
         #include "chemistry.H"
         #include "rhoEqn.H"
@@ -100,12 +101,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/uncoupledKinematicParcelFoam/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,7 +16,7 @@
             IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
-        thermo->rho()
+        thermo.rho()
     );
 
     Info<< "\nReading field U\n" << endl;
@@ -42,7 +43,7 @@
             rho,
             U,
             phi,
-            thermo()
+            thermo
         )
     );
 
@@ -55,6 +56,6 @@
         kinematicCloudName,
         rho,
         U,
-        thermo().mu(),
+        thermo.mu(),
         g
     );

applications/solvers/lagrangian/uncoupledKinematicParcelFoam/uncoupledKinematicParcelFoam.C

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

applications/solvers/multiphase/bubbleFoam/bubbleFoam.C

@@ -40,16 +40,15 @@ Description
 
 int main(int argc, char *argv[])
 {
+    #include "setRootCase.H"
 
-#   include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "readEnvironmentalProperties.H"
+    #include "createFields.H"
+    #include "initContinuityErrs.H"
 
-#   include "createTime.H"
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-#   include "createMesh.H"
-#   include "readEnvironmentalProperties.H"
-#   include "createFields.H"
-#   include "initContinuityErrs.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
 
@@ -57,30 +56,30 @@ int main(int argc, char *argv[])
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-#       include "readBubbleFoamControls.H"
+        #include "readBubbleFoamControls.H"
-#       include "CourantNo.H"
+        #include "CourantNo.H"
 
-#       include "alphaEqn.H"
+        #include "alphaEqn.H"
-#       include "liftDragCoeffs.H"
+        #include "liftDragCoeffs.H"
-#       include "UEqns.H"
+        #include "UEqns.H"
 
         // --- PISO loop
 
         for (int corr=0; corr<nCorr; corr++)
         {
-#           include "pEqn.H"
+            #include "pEqn.H"
 
             if (correctAlpha)
             {
-#               include "alphaEqn.H"
+                #include "alphaEqn.H"
             }
         }
 
-#       include "DDtU.H"
+        #include "DDtU.H"
 
-#       include "kEpsilon.H"
+        #include "kEpsilon.H"
 
-#       include "write.H"
+        #include "write.H"
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
             << "  ClockTime = " << runTime.elapsedClockTime() << " s"

applications/solvers/multiphase/bubbleFoam/createFields.H

@@ -141,8 +141,8 @@
         transportProperties.lookup("Ct")
     );
 
-#   include "createPhia.H"
+    #include "createPhia.H"
-#   include "createPhib.H"
+    #include "createPhib.H"
 
     surfaceScalarField phi
     (
@@ -157,7 +157,7 @@
     );
 
 
-#   include "createRASTurbulence.H"
+    #include "createRASTurbulence.H"
 
     Info<< "Calculating field DDtUa and DDtUb\n" << endl;