merge with master + fix conflicts in tuts/mesh

.gitignore

@@ -7,6 +7,7 @@
 *.orig
 *.bak
 \#*\#
+.directory
 
 # CVS recovered versions - anywhere
 .#*
@@ -48,7 +49,9 @@ doc/[Dd]oxygen/man
 # source packages - anywhere
 *.tar.bz2
 *.tar.gz
+*.tar
 *.tgz
+*.gtgz
 
 # ignore the persistent .build tag in the main directory
 /.build

applications/solvers/DNS/dnsFoam/dnsFoam.C

@@ -41,17 +41,16 @@ Description
 
 int main(int argc, char *argv[])
 {
+    #include "setRootCase.H"
 
-#   include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMeshNoClear.H"
+    #include "readTransportProperties.H"
+    #include "createFields.H"
+    #include "readTurbulenceProperties.H"
+    #include "initContinuityErrs.H"
 
-#   include "createTime.H"
-#   include "createMeshNoClear.H"
-#   include "readTransportProperties.H"
-#   include "createFields.H"
-#   include "readTurbulenceProperties.H"
-#   include "initContinuityErrs.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< nl << "Starting time loop" << endl;
 
@@ -59,7 +58,7 @@ int main(int argc, char *argv[])
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-#       include "readPISOControls.H"
+        #include "readPISOControls.H"
 
         force.internalField() = ReImSum
         (
@@ -69,12 +68,12 @@ int main(int argc, char *argv[])
             )
         );
 
-#       include "globalProperties.H"
+        #include "globalProperties.H"
 
         fvVectorMatrix UEqn
         (
-            fvm::ddt(U) 
-          + fvm::div(phi, U) 
+            fvm::ddt(U)
+          + fvm::div(phi, U)
           - fvm::laplacian(nu, U)
          ==
             force
@@ -90,7 +89,7 @@ int main(int argc, char *argv[])
             volScalarField rUA = 1.0/UEqn.A();
 
             U = rUA*UEqn.H();
-            phi = (fvc::interpolate(U) & mesh.Sf()) 
+            phi = (fvc::interpolate(U) & mesh.Sf())
                 + fvc::ddtPhiCorr(rUA, U, phi);
 
             fvScalarMatrix pEqn
@@ -102,7 +101,7 @@ int main(int argc, char *argv[])
 
             phi -= pEqn.flux();
 
-#           include "continuityErrs.H"
+            #include "continuityErrs.H"
 
             U -= rUA*fvc::grad(p);
             U.correctBoundaryConditions();

applications/solvers/Lagrangian/kinematicParcelFoam/Make/files

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

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,52 @@ Description
 
 int main(int argc, char *argv[])
 {
-#   include "setRootCase.H"
+    #include "setRootCase.H"
 
-#   include "createTime.H"
-#   include "createMesh.H"
-#   include "readCombustionProperties.H"
-#   include "readEnvironmentalProperties.H"
-#   include "createFields.H"
-#   include "readPISOControls.H"
-#   include "initContinuityErrs.H"
-#   include "readTimeControls.H"
-#   include "CourantNo.H"
-#   include "setInitialDeltaT.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "readCombustionProperties.H"
+    #include "readEnvironmentalProperties.H"
+    #include "createFields.H"
+    #include "initContinuityErrs.H"
+    #include "readTimeControls.H"
+    #include "CourantNo.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 "readPISOControls.H"
-#       include "CourantNo.H"
-#       include "setDeltaT.H"
+        #include "readTimeControls.H"
+        #include "readPISOControls.H"
+        #include "CourantNo.H"
+        #include "setDeltaT.H"
 
         runTime++;
 
         Info<< "\n\nTime = " << runTime.timeName() << endl;
 
-#       include "rhoEqn.H"
-#       include "UEqn.H"
+        #include "rhoEqn.H"
+        #include "UEqn.H"
 
         // --- PISO loop
         for (int corr=1; corr<=nCorr; corr++)
         {
-#           include "bEqn.H"
-#           include "ftEqn.H"
-#           include "huEqn.H"
-#           include "hEqn.H"
+            #include "bEqn.H"
+            #include "ftEqn.H"
+            #include "huEqn.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())
-{
-    surfaceScalarField SfUfbyDelta = 
-        mesh.surfaceInterpolation::deltaCoeffs()
-       *mag(phiSt/fvc::interpolate(rho));
+    if (mesh.nInternalFaces())
+    {
+        surfaceScalarField SfUfbyDelta =
+            mesh.surfaceInterpolation::deltaCoeffs()
+        *mag(phiSt/fvc::interpolate(rho));
 
-    StCoNum = max(SfUfbyDelta/mesh.magSf())
-        .value()*runTime.deltaT().value();
+        StCoNum =
+            max(SfUfbyDelta/mesh.magSf()).value()
+           *runTime.deltaT().value();
 
-    meanStCoNum = (sum(SfUfbyDelta)/sum(mesh.magSf()))
-        .value()*runTime.deltaT().value();
-}
+        meanStCoNum =
+            (sum(SfUfbyDelta)/sum(mesh.magSf())).value()
+        *runTime.deltaT().value();
+    }
 
-Info<< "St courant Number mean: " << meanStCoNum
-    << " max: " << StCoNum << endl;
+    Info<< "St courant Number mean: " << meanStCoNum
+        << " max: " << StCoNum << endl;
 }
 
 // ************************************************************************* //

applications/solvers/combustion/PDRFoam/bEqn.H

@@ -1,7 +1,7 @@
 tmp<fv::convectionScheme<scalar> > mvConvection
 (
     fv::convectionScheme<scalar>::New
-    (    
+    (
         mesh,
         fields,
         phi,
@@ -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();
-    const volScalarField& psi = thermo->psi();
-    volScalarField& h = thermo->h();
-    volScalarField& hu = thermo->hu();
+    volScalarField& p = thermo.p();
+    const volScalarField& psi = thermo.psi();
+    volScalarField& h = thermo.h();
+    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)
@@ -34,7 +34,7 @@ if (transonic)
 }
 else
 {
-    phi = 
+    phi =
         fvc::interpolate(rho)*
         (
             (fvc::interpolate(U) & mesh.Sf())

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

@@ -61,55 +61,54 @@ Description
 
 int main(int argc, char *argv[])
 {
-#   include "setRootCase.H"
+    #include "setRootCase.H"
 
-#   include "createTime.H"
-#   include "createMesh.H"
-#   include "readCombustionProperties.H"
-#   include "readEnvironmentalProperties.H"
-#   include "createFields.H"
-#   include "readPISOControls.H"
-#   include "initContinuityErrs.H"
-#   include "readTimeControls.H"
-#   include "compressibleCourantNo.H"
-#   include "setInitialDeltaT.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "readCombustionProperties.H"
+    #include "readEnvironmentalProperties.H"
+    #include "createFields.H"
+    #include "initContinuityErrs.H"
+    #include "readTimeControls.H"
+    #include "compressibleCourantNo.H"
+    #include "setInitialDeltaT.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
 
     while (runTime.run())
     {
-#       include "readTimeControls.H"
-#       include "readPISOControls.H"
-#       include "compressibleCourantNo.H"
-#       include "setDeltaT.H"
+        #include "readTimeControls.H"
+        #include "readPISOControls.H"
+        #include "compressibleCourantNo.H"
+        #include "setDeltaT.H"
 
         runTime++;
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-#       include "rhoEqn.H"
-#       include "UEqn.H"
+        #include "rhoEqn.H"
+        #include "UEqn.H"
 
         // --- PISO loop
         for (int corr=1; corr<=nCorr; corr++)
         {
-#           include "ftEqn.H"
-#           include "bEqn.H"
-#           include "huEqn.H"
-#           include "hEqn.H"
+            #include "ftEqn.H"
+            #include "bEqn.H"
+            #include "huEqn.H"
+            #include "hEqn.H"
 
             if (!ign.ignited())
             {
                 hu == h;
             }
 
-#           include "pEqn.H"
+            #include "pEqn.H"
         }
 
         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/
     (
@@ -180,7 +180,7 @@ if (ign.ignited())
         // with a linear correction function to give a plausible profile for Xi
         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-        volScalarField XiEqStar = 
+        volScalarField XiEqStar =
             scalar(1.001) + XiCoef*sqrt(up/(Su + SuMin))*Reta;
 
         volScalarField XiEq =

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();
-    const volScalarField& psi = thermo->psi();
-    volScalarField& h = thermo->h();
-    volScalarField& hu = thermo->hu();
+    volScalarField& p = thermo.p();
+    const volScalarField& psi = thermo.psi();
+    volScalarField& h = thermo.h();
+    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)
@@ -34,9 +34,9 @@ if (transonic)
 }
 else
 {
-    phi = 
-        fvc::interpolate(rho)*
-        (
+    phi =
+        fvc::interpolate(rho)
+       *(
             (fvc::interpolate(U) & mesh.Sf())
           + fvc::ddtPhiCorr(rUA, rho, U, phi)
         );

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 "createEngineMesh.H"
-#   include "createFields.H"
-#   include "initContinuityErrs.H"
-#   include "readEngineTimeControls.H"
-#   include "compressibleCourantNo.H"
-#   include "setInitialDeltaT.H"
-#   include "startSummary.H"
+    #include "createEngineTime.H"
+    #include "createEngineMesh.H"
+    #include "createFields.H"
+    #include "initContinuityErrs.H"
+    #include "readEngineTimeControls.H"
+    #include "compressibleCourantNo.H"
+    #include "setInitialDeltaT.H"
+    #include "startSummary.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
 
     while (runTime.run())
     {
-#       include "readPISOControls.H"
-#       include "readEngineTimeControls.H"
-#       include "compressibleCourantNo.H"
-#       include "setDeltaT.H"
+        #include "readPISOControls.H"
+        #include "readEngineTimeControls.H"
+        #include "compressibleCourantNo.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 "pEqn.H"
+            #include "hEqn.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();
-    const volScalarField& psi = thermo->psi();
-    volScalarField& h = thermo->h();
-    const volScalarField& T = thermo->T();
+    volScalarField& p = thermo.p();
+    const volScalarField& psi = thermo.psi();
+    volScalarField& h = thermo.h();
+    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)
         {
@@ -39,8 +39,7 @@ tmp<fv::convectionScheme<scalar> > mvConvection
             inertIndex = i;
         }
     }
-    
+
     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();
-const volScalarField& psi = thermo->psi();
-const volScalarField& T = thermo->T();
-volScalarField& h = thermo->h();
+volScalarField& p = thermo.p();
+const volScalarField& psi = thermo.psi();
+const volScalarField& T = thermo.T();
+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 "createEngineTime.H"
-#   include "createEngineMesh.H"
-#   include "createFields.H"
-#   include "readEnvironmentalProperties.H"
-#   include "readCombustionProperties.H"
-#   include "createSpray.H"
-#   include "initContinuityErrs.H"
-#   include "readEngineTimeControls.H"
-#   include "compressibleCourantNo.H"
-#   include "setInitialDeltaT.H"
-#   include "startSummary.H"
+    #include "setRootCase.H"
+    #include "createEngineTime.H"
+    #include "createEngineMesh.H"
+    #include "createFields.H"
+    #include "readEnvironmentalProperties.H"
+    #include "readCombustionProperties.H"
+    #include "createSpray.H"
+    #include "initContinuityErrs.H"
+    #include "readEngineTimeControls.H"
+    #include "compressibleCourantNo.H"
+    #include "setInitialDeltaT.H"
+    #include "startSummary.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info << "\nStarting time loop\n" << endl;
 
     while (runTime.run())
     {
-#       include "readPISOControls.H"
-#       include "readEngineTimeControls.H"
-#       include "compressibleCourantNo.H"
-#       include "setDeltaT.H"
+        #include "readPISOControls.H"
+        #include "readEngineTimeControls.H"
+        #include "compressibleCourantNo.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 "UEqn.H"
+        #include "rhoEqn.H"
+        #include "UEqn.H"
 
         for (label ocorr=1; ocorr <= nOuterCorr; ocorr++)
         {
-#           include "YEqn.H"
-#           include "hEqn.H"
+            #include "YEqn.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 "spraySummary.H"
+        #include "logSummary.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();
-
-    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];
-    }
+    thermo.correct();
 }

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

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 "compressibleCourantNo.H"
-#       include "setDeltaT.H"
+        #include "readPISOControls.H"
+        #include "compressibleCourantNo.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 "UEqn.H"
+        #include "rhoEqn.H"
+        #include "UEqn.H"
 
         for (label ocorr=1; ocorr <= nOuterCorr; ocorr++)
         {
-#           include "YEqn.H"
-#           include "hEqn.H"
+            #include "YEqn.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)
@@ -36,9 +36,9 @@ if (transonic)
 }
 else
 {
-    phi = 
-        fvc::interpolate(rho)*
-        (
+    phi =
+        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,28 @@ Description
 
 int main(int argc, char *argv[])
 {
-#   include "setRootCase.H"
+    #include "setRootCase.H"
 
-#   include "createEngineTime.H"
-#   include "createEngineMesh.H"
-#   include "readPISOControls.H"
-#   include "readCombustionProperties.H"
-#   include "createFields.H"
-#   include "initContinuityErrs.H"
-#   include "readEngineTimeControls.H"
-#   include "compressibleCourantNo.H"
-#   include "setInitialDeltaT.H"
-#   include "startSummary.H"
+    #include "createEngineTime.H"
+    #include "createEngineMesh.H"
+    #include "readCombustionProperties.H"
+    #include "createFields.H"
+    #include "initContinuityErrs.H"
+    #include "readEngineTimeControls.H"
+    #include "compressibleCourantNo.H"
+    #include "setInitialDeltaT.H"
+    #include "startSummary.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info << "\nStarting time loop\n" << endl;
 
     while (runTime.run())
     {
-#       include "readPISOControls.H"
-#       include "readEngineTimeControls.H"
-#       include "compressibleCourantNo.H"
-#       include "setDeltaT.H"
+        #include "readPISOControls.H"
+        #include "readEngineTimeControls.H"
+        #include "compressibleCourantNo.H"
+        #include "setDeltaT.H"
 
         runTime++;
 
@@ -93,31 +92,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 "bEqn.H"
-#           include "huEqn.H"
-#           include "hEqn.H"
+            #include "ftEqn.H"
+            #include "bEqn.H"
+            #include "huEqn.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(U) & mesh.Sf()) - fvc::meshPhi(rho, U))
+        fvc::interpolate(psi)
+       *((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/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)
         {
@@ -37,7 +37,7 @@ tmp<fv::convectionScheme<scalar> > mvConvection
             inertIndex = i;
         }
     }
-    
+
     Y[inertIndex] = scalar(1) - Yt;
     Y[inertIndex].max(0.0);
 }

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,10 +38,9 @@ volVectorField U
 );
 
 
-volScalarField& p = thermo->p();
-const volScalarField& psi = thermo->psi();
-const volScalarField& T = thermo->T();
-volScalarField& h = thermo->h();
+volScalarField& p = thermo.p();
+const volScalarField& psi = thermo.psi();
+volScalarField& h = thermo.h();
 
 
 #include "compressibleCreatePhi.H"
@@ -65,7 +67,7 @@ autoPtr<compressible::turbulenceModel> turbulence
         rho,
         U,
         phi,
-        thermo()
+        thermo
     )
 );
 
@@ -73,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 "createTime.H"
-#   include "createMesh.H"
-#   include "readChemistryProperties.H"
-#   include "readEnvironmentalProperties.H"
-#   include "createFields.H"
-#   include "initContinuityErrs.H"
-#   include "readTimeControls.H"
-#   include "compressibleCourantNo.H"
-#   include "setInitialDeltaT.H"
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "readChemistryProperties.H"
+    #include "readEnvironmentalProperties.H"
+    #include "createFields.H"
+    #include "initContinuityErrs.H"
+    #include "readTimeControls.H"
+    #include "compressibleCourantNo.H"
+    #include "setInitialDeltaT.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info << "\nStarting time loop\n" << endl;
 
     while (runTime.run())
     {
-#       include "readTimeControls.H"
-#       include "readPISOControls.H"
-#       include "compressibleCourantNo.H"
-#       include "setDeltaT.H"
+        #include "readTimeControls.H"
+        #include "readPISOControls.H"
+        #include "compressibleCourantNo.H"
+        #include "setDeltaT.H"
 
         runTime++;
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-#       include "chemistry.H"
-#       include "rhoEqn.H"
-#       include "UEqn.H"
+        #include "chemistry.H"
+        #include "rhoEqn.H"
+        #include "UEqn.H"
 
         for (label ocorr=1; ocorr <= nOuterCorr; ocorr++)
         {
-#           include "YEqn.H"
+            #include "YEqn.H"
 
-#           define Db turbulence->alphaEff()
-#           include "hEqn.H"
+            #define Db turbulence->alphaEff()
+            #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/combustion/rhoReactingFoam/Make/files

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

applications/solvers/combustion/rhoReactingFoam/Make/options

@@ -0,0 +1,19 @@
+EXE_INC = \
+    -I../XiFoam \
+    -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
+    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
+    -I$(LIB_SRC)/ODE/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude
+
+EXE_LIBS = \
+    -lcompressibleRASModels \
+    -lcompressibleLESModels \
+    -lreactionThermophysicalModels \
+    -lspecie \
+    -lbasicThermophysicalModels \
+    -lchemistryModel \
+    -lODE \
+    -lfiniteVolume

applications/solvers/combustion/rhoReactingFoam/YEqn.H

@@ -0,0 +1,43 @@
+tmp<fv::convectionScheme<scalar> > mvConvection
+(
+    fv::convectionScheme<scalar>::New
+    (
+        mesh,
+        fields,
+        phi,
+        mesh.divScheme("div(phi,Yi_h)")
+    )
+);
+
+{
+    label inertIndex = -1;
+    volScalarField Yt = 0.0*Y[0];
+
+    for (label i=0; i<Y.size(); i++)
+    {
+        if (Y[i].name() != inertSpecie)
+        {
+            volScalarField& Yi = Y[i];
+
+            solve
+            (
+                fvm::ddt(rho, Yi)
+              + mvConvection->fvmDiv(phi, Yi)
+              - fvm::laplacian(turbulence->muEff(), Yi)
+             ==
+                kappa*chemistry.RR(i),
+                mesh.solver("Yi")
+            );
+
+            Yi.max(0.0);
+            Yt += Yi;
+        }
+        else
+        {
+            inertIndex = i;
+        }
+    }
+
+    Y[inertIndex] = scalar(1) - Yt;
+    Y[inertIndex].max(0.0);
+}

applications/solvers/combustion/rhoReactingFoam/chemistry.H

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

applications/solvers/combustion/rhoReactingFoam/createFields.H

@@ -0,0 +1,85 @@
+Info<< nl << "Reading thermophysicalProperties" << endl;
+autoPtr<rhoChemistryModel> pChemistry
+(
+    rhoChemistryModel::New(mesh)
+);
+rhoChemistryModel& chemistry = pChemistry();
+
+hReactionThermo& thermo = chemistry.thermo();
+
+basicMultiComponentMixture& composition = thermo.composition();
+PtrList<volScalarField>& Y = composition.Y();
+
+word inertSpecie(thermo.lookup("inertSpecie"));
+
+volScalarField rho
+(
+    IOobject
+    (
+        "rho",
+        runTime.timeName(),
+        mesh
+    ),
+    thermo.rho()
+);
+
+Info<< "Reading field U\n" << endl;
+volVectorField U
+(
+    IOobject
+    (
+        "U",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);
+
+
+volScalarField& p = thermo.p();
+const volScalarField& psi = thermo.psi();
+volScalarField& h = thermo.h();
+
+
+#include "compressibleCreatePhi.H"
+
+volScalarField kappa
+(
+    IOobject
+    (
+        "kappa",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("zero", dimless, 0.0)
+);
+
+Info << "Creating turbulence model.\n" << nl;
+autoPtr<compressible::turbulenceModel> turbulence
+(
+    compressible::turbulenceModel::New
+    (
+        rho,
+        U,
+        phi,
+        thermo
+    )
+);
+
+Info<< "Creating field DpDt\n" << endl;
+volScalarField DpDt =
+    fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
+
+multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
+
+forAll (Y, i)
+{
+    fields.add(Y[i]);
+}
+fields.add(h);
+

applications/solvers/combustion/rhoReactingFoam/pEqn.H

@@ -0,0 +1,93 @@
+{
+    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();
+    U = rUA*UEqn.H();
+
+    if (transonic)
+    {
+        surfaceScalarField phiv =
+            (fvc::interpolate(U) & mesh.Sf())
+          + fvc::ddtPhiCorr(rUA, rho, U, phi);
+
+        phi = fvc::interpolate(rho)*phiv;
+
+        surfaceScalarField phid
+        (
+            "phid",
+            fvc::interpolate(thermo.psi())*phiv
+        );
+
+        for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+        {
+            fvScalarMatrix pEqn
+            (
+                fvc::ddt(rho) + fvc::div(phi)
+              + correction(fvm::ddt(psi, p) + fvm::div(phid, p))
+              - fvm::laplacian(rho*rUA, p)
+            );
+
+            if (ocorr == nOuterCorr && corr == nCorr && nonOrth == nNonOrthCorr)
+            {
+                pEqn.solve(mesh.solver(p.name() + "Final"));
+            }
+            else
+            {
+                pEqn.solve();
+            }
+
+            if (nonOrth == nNonOrthCorr)
+            {
+                phi += pEqn.flux();
+            }
+        }
+    }
+    else
+    {
+        phi =
+            fvc::interpolate(rho)
+           *(
+                (fvc::interpolate(U) & mesh.Sf())
+              + fvc::ddtPhiCorr(rUA, rho, U, phi)
+            );
+
+        for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+        {
+            fvScalarMatrix pEqn
+            (
+                fvc::ddt(rho) + psi*correction(fvm::ddt(p))
+              + fvc::div(phi)
+              - fvm::laplacian(rho*rUA, p)
+            );
+
+            if (ocorr == nOuterCorr && corr == nCorr && nonOrth == nNonOrthCorr)
+            {
+                pEqn.solve(mesh.solver(p.name() + "Final"));
+            }
+            else
+            {
+                pEqn.solve();
+            }
+
+            if (nonOrth == nNonOrthCorr)
+            {
+                phi += pEqn.flux();
+            }
+        }
+    }
+
+    // Second part of thermodynamic density update
+    thermo.rho() += psi*p;
+
+    #include "rhoEqn.H"
+    #include "compressibleContinuityErrs.H"
+
+    U -= rUA*fvc::grad(p);
+    U.correctBoundaryConditions();
+
+    DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
+}

applications/solvers/combustion/rhoReactingFoam/readChemistryProperties.H

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

applications/solvers/combustion/rhoReactingFoam/rhoReactingFoam.C

@@ -0,0 +1,102 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Application
+    rhoReactingFoam
+
+Description
+    Chemical reaction code using density based thermodynamics package.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "hReactionThermo.H"
+#include "turbulenceModel.H"
+#include "rhoChemistryModel.H"
+#include "chemistrySolver.H"
+#include "multivariateScheme.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "readChemistryProperties.H"
+    #include "readEnvironmentalProperties.H"
+    #include "createFields.H"
+    #include "initContinuityErrs.H"
+    #include "readTimeControls.H"
+    #include "compressibleCourantNo.H"
+    #include "setInitialDeltaT.H"
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info << "\nStarting time loop\n" << endl;
+
+    while (runTime.run())
+    {
+        #include "readTimeControls.H"
+        #include "readPISOControls.H"
+        #include "compressibleCourantNo.H"
+        #include "setDeltaT.H"
+
+        runTime++;
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        #include "chemistry.H"
+        #include "rhoEqn.H"
+        #include "UEqn.H"
+
+        for (label ocorr=1; ocorr <= nOuterCorr; ocorr++)
+        {
+            #include "YEqn.H"
+            #include "hEqn.H"
+
+            // --- PISO loop
+            for (int corr=1; corr<=nCorr; corr++)
+            {
+                #include "pEqn.H"
+            }
+        }
+
+        turbulence->correct();
+
+        rho = thermo.rho();
+
+        runTime.write();
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/solvers/compressible/rhoCentralFoam/compressibleCourantNo.H

@@ -38,11 +38,11 @@ if (mesh.nInternalFaces())
     surfaceScalarField amaxSfbyDelta =
         mesh.surfaceInterpolation::deltaCoeffs()*amaxSf;
 
-    CoNum = max(amaxSfbyDelta/mesh.magSf())
-        .value()*runTime.deltaT().value();
+    CoNum = max(amaxSfbyDelta/mesh.magSf()).value()*runTime.deltaT().value();
 
-    meanCoNum = (sum(amaxSfbyDelta)/sum(mesh.magSf()))
-        .value()*runTime.deltaT().value();
+    meanCoNum =
+        (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& h = thermo->h();
-const volScalarField& T = thermo->T();
-const volScalarField& psi = thermo->psi();
-const volScalarField& mu = thermo->mu();
+volScalarField& p = thermo.p();
+volScalarField& h = thermo.h();
+const volScalarField& T = thermo.T();
+const volScalarField& psi = thermo.psi();
+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
-    (
-        rhoBoundaryTypes[patchi] == "waveTransmissive"
-    )
+    if (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 "readTimeControls.H"
-#       include "setDeltaT.H"
+        #include "compressibleCourantNo.H"
+        #include "readTimeControls.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)
-              + fvc::laplacian(thermo->alpha(), h)
+              - fvm::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& h = thermo->h();
-    const volScalarField& psi = thermo->psi();
+    volScalarField& p = thermo.p();
+    volScalarField& h = thermo.h();
+    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)
@@ -58,8 +58,6 @@ else
         //+ fvc::ddtPhiCorr(rUA, rho, U, phi)
         );
 
-    //bool closedVolume = adjustPhi(phi, U, p);
-
     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
         // Pressure corrector
@@ -99,7 +97,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 +115,7 @@ bound(p, pMin);
 /*
 if (closedVolume)
 {
-    p += (initialMass - fvc::domainIntegrate(thermo->psi()*p))
-        /fvc::domainIntegrate(thermo->psi());
+    p += (initialMass - fvc::domainIntegrate(psi*p))
+        /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& h = thermo->h();
-    const volScalarField& psi = thermo->psi();
+    volScalarField& p = thermo.p();
+    volScalarField& h = thermo.h();
+    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)
@@ -34,7 +34,7 @@ if (transonic)
 }
 else
 {
-    phi = 
+    phi =
         fvc::interpolate(rho)*
         (
             (fvc::interpolate(U) & mesh.Sf())

applications/solvers/compressible/rhoPisoFoam/rhoPisoFoam.C

@@ -31,7 +31,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "basicThermo.H"
+#include "basicPsiThermo.H"
 #include "turbulenceModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -43,14 +43,13 @@ int main(int argc, char *argv[])
     #include "createTime.H"
     #include "createMesh.H"
     #include "createFields.H"
-    #include "readPISOControls.H"
     #include "initContinuityErrs.H"
     #include "readTimeControls.H"
     #include "compressibleCourantNo.H"
     #include "setInitialDeltaT.H"
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
 
@@ -77,7 +76,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& h = thermo->h();
+    volScalarField& p = thermo.p();
+    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))
-        /fvc::domainIntegrate(thermo->psi());
+    p += (initialMass - fvc::domainIntegrate(psi*p))
+        /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"
 
@@ -40,14 +40,13 @@ Description
 
 int main(int argc, char *argv[])
 {
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "createFields.H"
+    #include "initContinuityErrs.H"
 
-#   include "setRootCase.H"
-#   include "createTime.H"
-#   include "createMesh.H"
-#   include "createFields.H"
-#   include "initContinuityErrs.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
 
@@ -55,17 +54,17 @@ int main(int argc, char *argv[])
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-#       include "readSIMPLEControls.H"
-#       include "initConvergenceCheck.H"
+        #include "readSIMPLEControls.H"
+        #include "initConvergenceCheck.H"
 
         p.storePrevIter();
         rho.storePrevIter();
 
         // Pressure-velocity SIMPLE corrector
         {
-#           include "UEqn.H"
-#           include "hEqn.H"
-#           include "pEqn.H"
+            #include "UEqn.H"
+            #include "hEqn.H"
+            #include "pEqn.H"
         }
 
         turbulence->correct();
@@ -75,7 +74,7 @@ int main(int argc, char *argv[])
             << "  ClockTime = " << runTime.elapsedClockTime() << " s"
             << nl << endl;
 
-#       include "convergenceCheck.H"
+        #include "convergenceCheck.H"
     }
 
     Info<< "End\n" << endl;

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& h = thermo->h();
-
+    volScalarField& p = thermo.p();
+    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))
-        /fvc::domainIntegrate(thermo->psi());
+    p += (initialMass - fvc::domainIntegrate(psi*p))
+        /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/rhopSonicFoam/rhopSonicFoam.C

@@ -42,14 +42,13 @@ Description
 
 int main(int argc, char *argv[])
 {
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "readThermodynamicProperties.H"
+    #include "createFields.H"
 
-#   include "setRootCase.H"
-#   include "createTime.H"
-#   include "createMesh.H"
-#   include "readThermodynamicProperties.H"
-#   include "createFields.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
 
@@ -57,10 +56,10 @@ int main(int argc, char *argv[])
     {
         Info<< "Time = " << runTime.value() << nl << endl;
 
-#       include "readPISOControls.H"
+        #include "readPISOControls.H"
         scalar HbyAblend = readScalar(piso.lookup("HbyAblend"));
 
-#       include "readTimeControls.H"
+        #include "readTimeControls.H"
 
         scalar CoNum = max
         (
@@ -70,7 +69,7 @@ int main(int argc, char *argv[])
 
         Info<< "Max Courant Number = " << CoNum << endl;
 
-#       include "setDeltaT.H"
+        #include "setDeltaT.H"
 
         for (int outerCorr=0; outerCorr<nOuterCorr; outerCorr++)
         {
@@ -145,9 +144,9 @@ int main(int argc, char *argv[])
 
                 phi -= phiGradp;
 
-#               include "resetPhiPatches.H"
+                #include "resetPhiPatches.H"
 
-                surfaceScalarField rhof = 
+                surfaceScalarField rhof =
                     mvConvection.interpolationScheme()()(rho)()
                    .interpolate(rho);
 
@@ -165,7 +164,7 @@ int main(int argc, char *argv[])
 
                 phi += phiGradp + pEqn.flux();
                 rho = psi*p;
-                rhof = 
+                rhof =
                     mvConvection.interpolationScheme()()(rho)()
                    .interpolate(rho);
                 phiv = phi/rhof;

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& h = thermo->h();
-    const volScalarField& psi = thermo->psi();
+    volScalarField& p = thermo.p();
+    volScalarField& e = thermo.e();
+    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"
 
@@ -72,7 +72,7 @@ int main(int argc, char *argv[])
 
         solve(UEqn == -fvc::grad(p));
 
-        #include "hEqn.H"
+        #include "eEqn.H"
 
 
         // --- PISO loop
@@ -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/UEqn.H

@@ -0,0 +1,8 @@
+fvVectorMatrix UEqn
+(
+    fvm::ddt(rho, U)
+    + fvm::div(phi, U)
+    + turbulence->divDevRhoReff(U)
+);
+
+solve(UEqn == -fvc::grad(p));

applications/solvers/compressible/sonicFoam/compressibleContinuityErrs.H

@@ -1,12 +0,0 @@
-{
-#   include "rhoEqn.H"
-}
-{
-    scalar sumLocalContErr = (sum(mag(rho - psi*p))/sum(rho)).value();
-    scalar globalContErr = (sum(rho - psi*p)/sum(rho)).value();
-    cumulativeContErr += globalContErr;
-
-    Info<< "time step continuity errors : sum local = " << sumLocalContErr
-         << ", global = " << globalContErr
-         << ", cumulative = " << cumulativeContErr << endl;
-}

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& h = thermo->h();
-    const volScalarField& psi = thermo->psi();
+    volScalarField& p = thermo.p();
+    volScalarField& e = thermo.e();
+    const volScalarField& psi = thermo.psi();
 
     volScalarField rho
     (
@@ -17,7 +18,7 @@
             runTime.timeName(),
             mesh
         ),
-        thermo->rho()
+        thermo.rho()
     );
 
     Info<< "Reading field U\n" << endl;
@@ -34,7 +35,7 @@
         mesh
     );
 
-#   include "compressibleCreatePhi.H"
+    #include "compressibleCreatePhi.H"
 
 
     Info<< "Creating turbulence model\n" << endl;
@@ -45,10 +46,6 @@
             rho,
             U,
             phi,
-            thermo()
+            thermo
         )
     );
-
-    Info<< "Creating field DpDt\n" << endl;
-    volScalarField DpDt =
-        fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);

applications/solvers/compressible/sonicFoam/eEqn.H

@@ -0,0 +1,12 @@
+{
+    solve
+    (
+        fvm::ddt(rho, e)
+      + fvm::div(phi, e)
+      - fvm::laplacian(turbulence->alphaEff(), e)
+      ==
+      - p*fvc::div(phi/fvc::interpolate(rho))
+    );
+
+    thermo.correct();
+}

applications/solvers/compressible/sonicFoam/hEqn.H

@@ -1,12 +0,0 @@
-{
-        solve
-        (
-            fvm::ddt(rho, h)
-          + fvm::div(phi, h)
-          - fvm::laplacian(turbulence->alphaEff(), h)
-         ==
-            DpDt
-        );
-
-        thermo->correct();
-}

applications/solvers/compressible/sonicFoam/pEqn.H

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

applications/solvers/compressible/sonicFoam/readThermodynamicProperties.H

@@ -1,23 +0,0 @@
-    Info<< "Reading thermodynamicProperties\n" << endl;
-
-    IOdictionary thermodynamicProperties
-    (
-        IOobject
-        (
-            "thermodynamicProperties",
-            runTime.constant(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE
-        )
-    );
-
-    dimensionedScalar R
-    (
-        thermodynamicProperties.lookup("R")
-    );
-
-    dimensionedScalar Cv
-    (
-        thermodynamicProperties.lookup("Cv")
-    );

applications/solvers/compressible/sonicFoam/readTransportProperties.H

@@ -1,18 +0,0 @@
-    Info<< "Reading transportProperties\n" << endl;
-
-    IOdictionary transportProperties
-    (
-        IOobject
-        (
-            "transportProperties",
-            runTime.constant(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE
-        )
-    );
-
-    dimensionedScalar mu
-    (
-        transportProperties.lookup("mu")
-    );

applications/solvers/compressible/sonicFoam/sonicFoam.C

@@ -32,7 +32,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "basicThermo.H"
+#include "basicPsiThermo.H"
 #include "turbulenceModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -58,64 +58,21 @@ int main(int argc, char *argv[])
 
         #include "rhoEqn.H"
 
-        fvVectorMatrix UEqn
-        (
-            fvm::ddt(rho, U)
-          + fvm::div(phi, U)
-          + turbulence->divDevRhoReff(U)
-        );
+        #include "UEqn.H"
 
-        solve(UEqn == -fvc::grad(p));
-
-        #include "hEqn.H"
+        #include "eEqn.H"
 
 
         // --- PISO loop
 
         for (int corr=0; corr<nCorr; corr++)
         {
-            volScalarField rUA = 1.0/UEqn.A();
-            U = rUA*UEqn.H();
-
-            surfaceScalarField phid
-            (
-                "phid",
-                fvc::interpolate(thermo->psi())
-               *(
-                   (fvc::interpolate(U) & mesh.Sf())
-                 + fvc::ddtPhiCorr(rUA, rho, U, phi)
-               )
-            );
-
-            for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
-            {
-                fvScalarMatrix pEqn
-                (
-                    fvm::ddt(psi, p)
-                  + fvm::div(phid, p)
-                  - fvm::laplacian(rho*rUA, p)
-                );
-
-                pEqn.solve();
-
-                if (nonOrth == nNonOrthCorr)
-                {
-                    phi = pEqn.flux();
-                }
-            }
-
-             #include "compressibleContinuityErrs.H"
-
-            U -= rUA*fvc::grad(p);
-            U.correctBoundaryConditions();
+            #include "pEqn.H"
         }
 
-        DpDt =
-            fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
-
         turbulence->correct();
 
-        rho = psi*p;
+        rho = thermo.rho();
 
         runTime.write();
 

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 "compressibleCourantNo.H"
+        #include "readPISOControls.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/discreteMethods/dsmc/dsmcFoam/dsmcFoam.C

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

applications/solvers/discreteMethods/molecularDynamics/mdEquilibrationFoam/mdEquilibrationFoam.C

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

applications/solvers/discreteMethods/molecularDynamics/mdEquilibrationFoam/readmdEquilibrationDict.H

@@ -14,5 +14,5 @@ IOdictionary mdEquilibrationDict
 
 scalar targetTemperature = readScalar
 (
-    mdEquilibrationDict.lookup("equilibrationTargetTemperature")
+    mdEquilibrationDict.lookup("targetTemperature")
 );

applications/solvers/discreteMethods/molecularDynamics/mdFoam/mdFoam.C

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

applications/solvers/electromagnetics/electrostaticFoam/electrostaticFoam.C

@@ -36,14 +36,13 @@ Description
 
 int main(int argc, char *argv[])
 {
+    #include "setRootCase.H"
 
-#   include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "createFields.H"
 
-#   include "createTime.H"
-#   include "createMesh.H"
-#   include "createFields.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting iteration loop\n" << endl;
 

applications/solvers/electromagnetics/mhdFoam/mhdFoam.C

@@ -26,7 +26,7 @@ Application
     mhdFoam
 
 Description
-    Solver for magnetohydrodynamics (MHD): incompressible, laminar flow of a 
+    Solver for magnetohydrodynamics (MHD): incompressible, laminar flow of a
     conducting fluid under the influence of a magnetic field.
 
     An applied magnetic field H acts as a driving force,
@@ -58,27 +58,25 @@ Description
 
 int main(int argc, char *argv[])
 {
+    #include "setRootCase.H"
 
-#   include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "createFields.H"
+    #include "initContinuityErrs.H"
 
-#   include "createTime.H"
-#   include "createMesh.H"
-#   include "createFields.H"
-#   include "initContinuityErrs.H"
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< nl << "Starting time loop" << endl;
 
     while (runTime.loop())
     {
-#       include "readPISOControls.H"
-#       include "readBPISOControls.H"
+        #include "readPISOControls.H"
+        #include "readBPISOControls.H"
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-#       include "CourantNo.H"
+        #include "CourantNo.H"
 
         {
             fvVectorMatrix UEqn
@@ -101,7 +99,7 @@ int main(int argc, char *argv[])
 
                 U = rUA*UEqn.H();
 
-                phi = (fvc::interpolate(U) & mesh.Sf()) 
+                phi = (fvc::interpolate(U) & mesh.Sf())
                     + fvc::ddtPhiCorr(rUA, U, phi);
 
                 for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
@@ -120,7 +118,7 @@ int main(int argc, char *argv[])
                     }
                 }
 
-#               include "continuityErrs.H"
+                #include "continuityErrs.H"
 
                 U -= rUA*fvc::grad(p);
                 U.correctBoundaryConditions();
@@ -154,7 +152,7 @@ int main(int argc, char *argv[])
 
             phiB -= pBEqn.flux();
 
-#           include "magneticFieldErr.H"
+            #include "magneticFieldErr.H"
         }
 
         runTime.write();

applications/solvers/financial/financialFoam/financialFoam.C

@@ -38,14 +38,13 @@ Description
 
 int main(int argc, char *argv[])
 {
+    #include "setRootCase.H"
 
-#include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "createFields.H"
 
-#include "createTime.H"
-#include "createMesh.H"
-#include "createFields.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< nl << "Calculating value(price of comodities)" << endl;