MRG: Integrated Foundation code to commit 47bd8e1

applications/solvers/combustion/chemFoam/output.H

@@ -1,6 +1,6 @@
     runTime.write();
 
-    Info<< "Sh = " << Sh
+    Info<< "Qdot = " << Qdot
         << ", T = " << thermo.T()[0]
         << ", p = " << thermo.p()[0]
         << ", " << Y[0].name() << " = " << Y[0][0]
@@ -8,4 +8,3 @@
 
     post<< runTime.value() << token::TAB << thermo.T()[0] << token::TAB
         << thermo.p()[0] << endl;
-

applications/solvers/combustion/chemFoam/solveChemistry.H

@@ -1,3 +1,3 @@
 dtChem = chemistry.solve(runTime.deltaT().value());
-    scalar Sh = chemistry.Sh()()[0]/rho[0];
-    integratedHeat += Sh*runTime.deltaT().value();
+scalar Qdot = chemistry.Qdot()()[0]/rho[0];
+integratedHeat += Qdot*runTime.deltaT().value();

applications/solvers/combustion/fireFoam/YEEqn.H

@@ -10,7 +10,7 @@ tmp<fv::convectionScheme<scalar>> mvConvection
 );
 {
     combustion->correct();
-    dQ = combustion->dQ();
+    Qdot = combustion->Qdot();
     volScalarField Yt(0.0*Y[0]);
 
     forAll(Y, i)
@@ -67,7 +67,7 @@ tmp<fv::convectionScheme<scalar>> mvConvection
         )
       - fvm::laplacian(turbulence->alphaEff(), he)
      ==
-        combustion->Sh()
+        Qdot
       + radiation->Sh(thermo)
       + parcels.Sh(he)
       + surfaceFilm.Sh()

applications/solvers/combustion/fireFoam/createFields.H

@@ -131,18 +131,18 @@ Switch solvePyrolysisRegion
     additionalControlsDict.lookupOrDefault<bool>("solvePyrolysisRegion", true)
 );
 
-volScalarField dQ
+volScalarField Qdot
 (
     IOobject
     (
-        "dQ",
+        "Qdot",
         runTime.timeName(),
         mesh,
-        IOobject::NO_READ,
+        IOobject::READ_IF_PRESENT,
         IOobject::AUTO_WRITE
     ),
     mesh,
-    dimensionedScalar("dQ", dimEnergy/dimTime, 0.0)
+    dimensionedScalar("Qdot", dimEnergy/dimVolume/dimTime, 0.0)
 );
 
 

applications/solvers/combustion/reactingFoam/EEqn.H

@@ -17,7 +17,7 @@
         )
       - fvm::laplacian(turbulence->alphaEff(), he)
      ==
-        reaction->Sh()
+        Qdot
       + fvOptions(rho, he)
     );
 

applications/solvers/combustion/reactingFoam/YEqn.H

@@ -11,7 +11,7 @@ tmp<fv::convectionScheme<scalar>> mvConvection
 
 {
     reaction->correct();
-    dQ = reaction->dQ();
+    Qdot = reaction->Qdot();
     volScalarField Yt(0.0*Y[0]);
 
     forAll(Y, i)

applications/solvers/combustion/reactingFoam/createFields.H

@@ -117,18 +117,18 @@ forAll(Y, i)
 }
 fields.add(thermo.he());
 
-volScalarField dQ
+volScalarField Qdot
 (
     IOobject
     (
-        "dQ",
+        "Qdot",
         runTime.timeName(),
         mesh,
-        IOobject::NO_READ,
+        IOobject::READ_IF_PRESENT,
         IOobject::AUTO_WRITE
     ),
     mesh,
-    dimensionedScalar("dQ", dimEnergy/dimTime, 0.0)
+    dimensionedScalar("Qdot", dimEnergy/dimVolume/dimTime, 0.0)
 );
 
 #include "createMRF.H"

applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/createFields.H

@@ -117,18 +117,18 @@ forAll(Y, i)
 }
 fields.add(thermo.he());
 
-volScalarField dQ
+volScalarField Qdot
 (
     IOobject
     (
-        "dQ",
+        "Qdot",
         runTime.timeName(),
         mesh,
-        IOobject::NO_READ,
+        IOobject::READ_IF_PRESENT,
         IOobject::AUTO_WRITE
     ),
     mesh,
-    dimensionedScalar("dQ", dimEnergy/dimTime, 0.0)
+    dimensionedScalar("Qdot", dimEnergy/dimVolume/dimTime, 0.0)
 );
 
 #include "createMRF.H"

applications/solvers/combustion/reactingFoam/rhoReactingFoam/createFields.H

@@ -96,18 +96,18 @@ forAll(Y, i)
 }
 fields.add(thermo.he());
 
-volScalarField dQ
+volScalarField Qdot
 (
     IOobject
     (
-        "dQ",
+        "Qdot",
         runTime.timeName(),
         mesh,
-        IOobject::NO_READ,
+        IOobject::READ_IF_PRESENT,
         IOobject::AUTO_WRITE
     ),
     mesh,
-    dimensionedScalar("dQ", dimEnergy/dimTime, 0.0)
+    dimensionedScalar("Qdot", dimEnergy/dimVolume/dimTime, 0.0)
 );
 
 #include "createMRF.H"

applications/solvers/combustion/reactingFoam/setRDeltaT.H

@@ -43,13 +43,16 @@ License
     // Damping coefficient (1-0)
     scalar rDeltaTDampingCoeff
     (
-        pimpleDict.lookupOrDefault<scalar>("rDeltaTDampingCoeff", 1)
+        pimpleDict.lookupOrDefault<scalar>("rDeltaTDampingCoeff", 1.0)
     );
 
     // Maximum change in cell temperature per iteration
     // (relative to previous value)
     scalar alphaTemp(pimpleDict.lookupOrDefault("alphaTemp", 0.05));
 
+    // Maximum change in cell concentration per iteration
+    // (relative to reference value)
+    scalar alphaY(pimpleDict.lookupOrDefault("alphaY", 1.0));
 
     Info<< "Time scales min/max:" << endl;
 
@@ -68,34 +71,89 @@ License
         rDeltaT.max(1/maxDeltaT);
 
         Info<< "    Flow        = "
-            << gMin(1/rDeltaT.primitiveField()) << ", "
-            << gMax(1/rDeltaT.primitiveField()) << endl;
+            << 1/gMax(rDeltaT.primitiveField()) << ", "
+            << 1/gMin(rDeltaT.primitiveField()) << endl;
     }
 
-    // Reaction source time scale
-    if (alphaTemp < 1.0)
+    // Heat release rate time scale
+    if (alphaTemp < 1)
     {
         volScalarField::Internal rDeltaTT
         (
-            mag(reaction->Sh())/(alphaTemp*rho*thermo.Cp()*T)
+            mag(Qdot)/(alphaTemp*rho*thermo.Cp()*T)
         );
 
         Info<< "    Temperature = "
-            << gMin(1/(rDeltaTT.field() + VSMALL)) << ", "
-            << gMax(1/(rDeltaTT.field() + VSMALL)) << endl;
+            << 1/(gMax(rDeltaTT.field()) + VSMALL) << ", "
+            << 1/(gMin(rDeltaTT.field()) + VSMALL) << endl;
 
-        rDeltaT.ref() = max
+        rDeltaT.ref() = max(rDeltaT(), rDeltaTT);
+    }
+
+    // Reaction rate time scale
+    if (alphaY < 1)
+    {
+        dictionary Yref(pimpleDict.subDict("Yref"));
+
+        volScalarField::Internal rDeltaTY
         (
-            rDeltaT(),
-            rDeltaTT
+            IOobject
+            (
+                "rDeltaTY",
+                runTime.timeName(),
+                mesh
+            ),
+            mesh,
+            dimensionedScalar("rDeltaTY", rDeltaT.dimensions(), 0)
         );
+
+        bool foundY = false;
+
+        forAll(Y, i)
+        {
+            if (i != inertIndex && composition.active(i))
+            {
+                volScalarField& Yi = Y[i];
+
+                if (Yref.found(Yi.name()))
+                {
+                    foundY = true;
+                    scalar Yrefi = readScalar(Yref.lookup(Yi.name()));
+
+                    rDeltaTY.field() = max
+                    (
+                        mag
+                        (
+                            reaction->R(Yi)().source()
+                           /((Yrefi*alphaY)*(rho*mesh.V()))
+                        ),
+                        rDeltaTY
+                    );
+                }
+            }
+        }
+
+        if (foundY)
+        {
+            Info<< "    Composition = "
+                << 1/(gMax(rDeltaTY.field()) + VSMALL) << ", "
+                << 1/(gMin(rDeltaTY.field()) + VSMALL) << endl;
+
+            rDeltaT.ref() = max(rDeltaT(), rDeltaTY);
+        }
+        else
+        {
+            IOWarningIn(args.executable().c_str(), Yref)
+                << "Cannot find any active species in Yref " << Yref
+                << endl;
+        }
     }
 
     // Update tho boundary values of the reciprocal time-step
     rDeltaT.correctBoundaryConditions();
 
     // Spatially smooth the time scale field
-    if (rDeltaTSmoothingCoeff < 1.0)
+    if (rDeltaTSmoothingCoeff < 1)
     {
         fvc::smooth(rDeltaT, rDeltaTSmoothingCoeff);
     }
@@ -105,7 +163,7 @@ License
     // - only increase at a fraction of old time scale
     if
     (
-        rDeltaTDampingCoeff < 1.0
+        rDeltaTDampingCoeff < 1
      && runTime.timeIndex() > runTime.startTimeIndex() + 1
     )
     {
@@ -120,8 +178,8 @@ License
     rDeltaT.correctBoundaryConditions();
 
     Info<< "    Overall     = "
-        << gMin(1/rDeltaT.primitiveField())
-        << ", " << gMax(1/rDeltaT.primitiveField()) << endl;
+        << 1/gMax(rDeltaT.primitiveField())
+        << ", " << 1/gMin(rDeltaT.primitiveField()) << endl;
 }
 
 

applications/solvers/lagrangian/coalChemistryFoam/EEqn.H

@@ -18,7 +18,7 @@
       - fvm::laplacian(turbulence->alphaEff(), he)
      ==
         rho*(U&g)
-      + combustion->Sh()
+      + Qdot
       + coalParcels.Sh(he)
       + limestoneParcels.Sh(he)
       + radiation->Sh(thermo)

applications/solvers/lagrangian/coalChemistryFoam/YEqn.H

@@ -12,7 +12,7 @@ tmp<fv::convectionScheme<scalar>> mvConvection
 
 {
     combustion->correct();
-    dQ = combustion->dQ();
+    Qdot = combustion->Qdot();
     volScalarField Yt(0.0*Y[0]);
 
     forAll(Y, i)

applications/solvers/lagrangian/coalChemistryFoam/createFields.H

@@ -131,18 +131,18 @@ volScalarField dpdt
 Info<< "Creating field kinetic energy K\n" << endl;
 volScalarField K("K", 0.5*magSqr(U));
 
-volScalarField dQ
+volScalarField Qdot
 (
     IOobject
     (
-        "dQ",
+        "Qdot",
         runTime.timeName(),
         mesh,
-        IOobject::NO_READ,
+        IOobject::READ_IF_PRESENT,
         IOobject::AUTO_WRITE
     ),
     mesh,
-    dimensionedScalar("dQ", dimEnergy/dimTime, 0.0)
+    dimensionedScalar("Qdot", dimEnergy/dimVolume/dimTime, 0.0)
 );
 
 #include "createMRF.H"

applications/solvers/lagrangian/coalChemistryFoam/setRDeltaT.H

@@ -80,7 +80,7 @@ License
             (
                (coalParcels.hsTrans() + limestoneParcels.hsTrans())
               /(mesh.V()*runTime.deltaT())
-             + combustion->Sh()()
+             + Qdot
             )
            /(
                alphaTemp

applications/solvers/lagrangian/reactingParcelFilmFoam/EEqn.H

@@ -21,7 +21,7 @@
       + parcels.Sh(he)
       + surfaceFilm.Sh()
       + radiation->Sh(thermo)
-      + combustion->Sh()
+      + Qdot
       + fvOptions(rho, he)
     );
 

applications/solvers/lagrangian/reactingParcelFilmFoam/YEqn.H

@@ -12,7 +12,7 @@ tmp<fv::convectionScheme<scalar>> mvConvection
 
 {
     combustion->correct();
-    dQ = combustion->dQ();
+    Qdot = combustion->Qdot();
     volScalarField Yt(0.0*Y[0]);
 
     forAll(Y, i)

applications/solvers/lagrangian/reactingParcelFilmFoam/createFields.H

@@ -134,18 +134,18 @@ Switch solvePrimaryRegion
     additionalControlsDict.lookup("solvePrimaryRegion")
 );
 
-volScalarField dQ
+volScalarField Qdot
 (
     IOobject
     (
-        "dQ",
+        "Qdot",
         runTime.timeName(),
         mesh,
-        IOobject::NO_READ,
+        IOobject::READ_IF_PRESENT,
         IOobject::AUTO_WRITE
     ),
     mesh,
-    dimensionedScalar("dQ", dimEnergy/dimTime, 0.0)
+    dimensionedScalar("Qdot", dimEnergy/dimVolume/dimTime, 0.0)
 );
 
 #include "createMRF.H"

applications/solvers/lagrangian/reactingParcelFoam/EEqn.H

@@ -20,7 +20,7 @@
         rho*(U&g)
       + parcels.Sh(he)
       + radiation->Sh(thermo)
-      + combustion->Sh()
+      + Qdot
       + fvOptions(rho, he)
     );
 

applications/solvers/lagrangian/reactingParcelFoam/YEqn.H

@@ -11,7 +11,7 @@ tmp<fv::convectionScheme<scalar>> mvConvection
 
 {
     combustion->correct();
-    dQ = combustion->dQ();
+    Qdot = combustion->Qdot();
     volScalarField Yt(0.0*Y[0]);
 
     forAll(Y, i)

applications/solvers/lagrangian/reactingParcelFoam/createFields.H

@@ -121,18 +121,18 @@ forAll(Y, i)
 }
 fields.add(thermo.he());
 
-volScalarField dQ
+volScalarField Qdot
 (
     IOobject
     (
-        "dQ",
+        "Qdot",
         runTime.timeName(),
         mesh,
-        IOobject::NO_READ,
+        IOobject::READ_IF_PRESENT,
         IOobject::AUTO_WRITE
     ),
     mesh,
-    dimensionedScalar("dQ", dimEnergy/dimTime, 0.0)
+    dimensionedScalar("Qdot", dimEnergy/dimVolume/dimTime, 0.0)
 );
 
 #include "createMRF.H"

applications/solvers/lagrangian/reactingParcelFoam/setRDeltaT.H

@@ -79,7 +79,7 @@ License
             mag
             (
                parcels.hsTrans()/(mesh.V()*runTime.deltaT())
-             + combustion->Sh()()
+             + Qdot
             )
            /(
                alphaTemp

applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/EEqn.H

@@ -14,7 +14,7 @@
         rho*(U&g)
       + parcels.Sh(he)
       + radiation->Sh(thermo)
-      + combustion->Sh()
+      + Qdot
       + fvOptions(rho, he)
     );
 

applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/YEqn.H

@@ -11,7 +11,7 @@ tmp<fv::convectionScheme<scalar>> mvConvection
 
 {
     combustion->correct();
-    dQ = combustion->dQ();
+    Qdot = combustion->Qdot();
     volScalarField Yt(0.0*Y[0]);
 
     forAll(Y, i)

applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/createFields.H

@@ -103,18 +103,18 @@ forAll(Y, i)
 }
 fields.add(thermo.he());
 
-volScalarField dQ
+volScalarField Qdot
 (
     IOobject
     (
-        "dQ",
+        "Qdot",
         runTime.timeName(),
         mesh,
-        IOobject::NO_READ,
+        IOobject::READ_IF_PRESENT,
         IOobject::AUTO_WRITE
     ),
     mesh,
-    dimensionedScalar("dQ", dimEnergy/dimTime, 0.0)
+    dimensionedScalar("Qdot", dimEnergy/dimVolume/dimTime, 0.0)
 );
 
 #include "createMRF.H"

applications/solvers/lagrangian/sprayFoam/createFields.H

@@ -117,18 +117,18 @@ forAll(Y, i)
 }
 fields.add(thermo.he());
 
-volScalarField dQ
+volScalarField Qdot
 (
     IOobject
     (
-        "dQ",
+        "Qdot",
         runTime.timeName(),
         mesh,
         IOobject::NO_READ,
         IOobject::AUTO_WRITE
     ),
     mesh,
-    dimensionedScalar("dQ", dimEnergy/dimTime, 0.0)
+    dimensionedScalar("Qdot", dimEnergy/dimVolume/dimTime, 0.0)
 );
 
 #include "createMRF.H"

applications/solvers/lagrangian/sprayFoam/sprayDyMFoam/sprayDyMFoam.C

@@ -138,7 +138,7 @@ int main(int argc, char *argv[])
 
         if (runTime.write())
         {
-            combustion->dQ()().write();
+            combustion->Qdot()().write();
         }
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"

applications/solvers/lagrangian/sprayFoam/sprayEngineFoam/sprayEngineFoam.C

@@ -113,7 +113,7 @@ int main(int argc, char *argv[])
 
         if (runTime.write())
         {
-            combustion->dQ()().write();
+            combustion->Qdot()().write();
         }
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"

applications/solvers/lagrangian/sprayFoam/sprayFoam.C

@@ -78,6 +78,8 @@ int main(int argc, char *argv[])
 
         parcels.evolve();
 
+        if (pimple.solveFlow())
+        {
             #include "rhoEqn.H"
 
             // --- Pressure-velocity PIMPLE corrector loop
@@ -103,7 +105,15 @@ int main(int argc, char *argv[])
 
             if (runTime.write())
             {
-            combustion->dQ()().write();
+                combustion->Qdot()().write();
+            }
+        }
+        else
+        {
+            if (runTime.writeTime())
+            {
+                parcels.write();
+            }
         }
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/AnisothermalPhaseModel/AnisothermalPhaseModel.C

@@ -135,7 +135,7 @@ Foam::AnisothermalPhaseModel<BasePhaseModel>::heEqn()
             he
         )
      ==
-        this->Sh()
+        this->Qdot()
     );
 
     // Add the appropriate pressure-work term

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/InertPhaseModel/InertPhaseModel.C

@@ -65,7 +65,7 @@ Foam::InertPhaseModel<BasePhaseModel>::R
 
 template<class BasePhaseModel>
 Foam::tmp<Foam::volScalarField>
-Foam::InertPhaseModel<BasePhaseModel>::Sh() const
+Foam::InertPhaseModel<BasePhaseModel>::Qdot() const
 {
     return tmp<volScalarField>
     (
@@ -73,7 +73,7 @@ Foam::InertPhaseModel<BasePhaseModel>::Sh() const
         (
             IOobject
             (
-                IOobject::groupName("Sh", this->name()),
+                IOobject::groupName("Qdot", this->name()),
                 this->mesh().time().timeName(),
                 this->mesh()
             ),

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/InertPhaseModel/InertPhaseModel.H

@@ -78,8 +78,8 @@ public:
                 volScalarField& Yi
             ) const;
 
-            //- Return the reaction heat source
-            virtual tmp<volScalarField> Sh() const;
+            //- Return the heat release rate
+            virtual tmp<volScalarField> Qdot() const;
 };
 
 

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/ReactingPhaseModel/ReactingPhaseModel.C

@@ -105,9 +105,9 @@ Foam::ReactingPhaseModel<BasePhaseModel, ReactionType>::R
 
 template<class BasePhaseModel, class ReactionType>
 Foam::tmp<Foam::volScalarField>
-Foam::ReactingPhaseModel<BasePhaseModel, ReactionType>::Sh() const
+Foam::ReactingPhaseModel<BasePhaseModel, ReactionType>::Qdot() const
 {
-    return reaction_->Sh();
+    return reaction_->Qdot();
 }
 
 

applications/solvers/multiphase/reactingEulerFoam/phaseSystems/phaseModel/ReactingPhaseModel/ReactingPhaseModel.H

@@ -87,8 +87,8 @@ public:
             volScalarField& Yi
         ) const;
 
-        //- Return the reacton heat source
-        virtual tmp<volScalarField> Sh() const;
+        //- Return heat release rate
+        virtual tmp<volScalarField> Qdot() const;
 };
 
 

applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pU/pEqn.H

@@ -7,8 +7,7 @@ volScalarField rAU1
     1.0
    /(
         U1Eqn.A()
-      + max(phase1.residualAlpha() - alpha1, scalar(0))
-       *rho1/runTime.deltaT()
+      + byDt(max(phase1.residualAlpha() - alpha1, scalar(0))*rho1)
     )
 );
 volScalarField rAU2
@@ -17,8 +16,7 @@ volScalarField rAU2
     1.0
    /(
         U2Eqn.A()
-      + max(phase2.residualAlpha() - alpha2, scalar(0))
-       *rho2/runTime.deltaT()
+      + byDt(max(phase2.residualAlpha() - alpha2, scalar(0))*rho2)
     )
 );
 
@@ -101,8 +99,8 @@ while (pimple.correct())
         rAU1
        *(
             U1Eqn.H()
-          + max(phase1.residualAlpha() - alpha1, scalar(0))
-           *rho1*U1.oldTime()/runTime.deltaT()
+          + byDt(max(phase1.residualAlpha() - alpha1, scalar(0))*rho1)
+           *U1.oldTime()
         );
 
     volVectorField HbyA2
@@ -114,8 +112,8 @@ while (pimple.correct())
         rAU2
        *(
             U2Eqn.H()
-         +  max(phase2.residualAlpha() - alpha2, scalar(0))
-           *rho2*U2.oldTime()/runTime.deltaT()
+         +  byDt(max(phase2.residualAlpha() - alpha2, scalar(0))*rho2)
+           *U2.oldTime()
         );
 
     surfaceScalarField ghSnGradRho
@@ -175,11 +173,9 @@ while (pimple.correct())
     (
         IOobject::groupName("phiHbyA", phase1.name()),
         fvc::flux(HbyA1)
-      + phiCorrCoeff1*fvc::interpolate(alpha1.oldTime()*rho1.oldTime()*rAU1)
-       *(
-            MRF.absolute(phi1.oldTime())
-          - fvc::flux(U1.oldTime())
-        )/runTime.deltaT()
+      + phiCorrCoeff1
+       *fvc::interpolate(byDt(alpha1.oldTime()*rho1.oldTime()*rAU1))
+       *(MRF.absolute(phi1.oldTime()) - fvc::flux(U1.oldTime()))
       - phiF1()
       - phig1
     );
@@ -189,11 +185,9 @@ while (pimple.correct())
     (
         IOobject::groupName("phiHbyA", phase2.name()),
         fvc::flux(HbyA2)
-      + phiCorrCoeff2*fvc::interpolate(alpha2.oldTime()*rho2.oldTime()*rAU2)
-       *(
-            MRF.absolute(phi2.oldTime())
-          - fvc::flux(U2.oldTime())
-        )/runTime.deltaT()
+      + phiCorrCoeff2
+       *fvc::interpolate(byDt(alpha2.oldTime()*rho2.oldTime()*rAU2))
+       *(MRF.absolute(phi2.oldTime()) - fvc::flux(U2.oldTime()))
       - phiF2()
       - phig2
     );

applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/DDtU.H

@@ -1,2 +1,2 @@
-ddtPhi1 = fvc::ddt(phi1);
-ddtPhi2 = fvc::ddt(phi2);
+tddtPhi1 = fvc::ddt(phi1);
+tddtPhi2 = fvc::ddt(phi2);

applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/createDDtU.H

@@ -1,2 +1,7 @@
-surfaceScalarField ddtPhi1(fvc::ddt(phi1));
-surfaceScalarField ddtPhi2(fvc::ddt(phi2));
+tmp<surfaceScalarField> tddtPhi1;
+tmp<surfaceScalarField> tddtPhi2;
+
+if (faceMomentum)
+{
+    #include "DDtU.H"
+}

applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/pUf/pEqn.H

@@ -159,7 +159,7 @@ while (pimple.correct())
             (alphaRhof10 + Vmf)
            *byDt(MRF.absolute(phi1.oldTime()))
           + fvc::flux(U1Eqn.H())
-          + Vmf*ddtPhi2
+          + Vmf*tddtPhi2()
           + Kdf*MRF.absolute(phi2)
           - Ff1()
         );
@@ -177,7 +177,7 @@ while (pimple.correct())
             (alphaRhof20 + Vmf)
            *byDt(MRF.absolute(phi2.oldTime()))
           + fvc::flux(U2Eqn.H())
-          + Vmf*ddtPhi1
+          + Vmf*tddtPhi1()
           + Kdf*MRF.absolute(phi1)
           - Ff2()
        );

applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/reactingTwoPhaseEulerFoam.C

@@ -45,6 +45,18 @@ Description
 
 namespace Foam
 {
+    tmp<volScalarField> byDt(const volScalarField& vf)
+    {
+        if (fv::localEulerDdt::enabled(vf.mesh()))
+        {
+            return fv::localEulerDdt::localRDeltaT(vf.mesh())*vf;
+        }
+        else
+        {
+            return vf/vf.mesh().time().deltaT();
+        }
+    }
+
     tmp<surfaceScalarField> byDt(const surfaceScalarField& sf)
     {
         if (fv::localEulerDdt::enabled(sf.mesh()))
@@ -92,7 +104,7 @@ int main(int argc, char *argv[])
         )
     );
 
-    #include "createRDeltaTf.H"
+    #include "pUf/createRDeltaTf.H"
     #include "pUf/createDDtU.H"
 
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/derivedFvPatchFields/wallBoilingSubModels/partitioningModels/cosine/cosine.C

@@ -83,7 +83,7 @@ cosine::fLiquid
                     1 - cos
                     (
                         constant::mathematical::pi
-                       *(alphaLiquid1_ - alphaLiquid)
+                       *(alphaLiquid - alphaLiquid0_)
                        /(alphaLiquid1_ - alphaLiquid0_)
                     )
                 )

applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/derivedFvPatchFields/wallBoilingSubModels/partitioningModels/cosine/cosine.H

@@ -32,13 +32,11 @@ Description
       - alphaLiquid0 0.05
 
     \verbatim
-        Ioilev, A., Samigulin, M., Ustinenko, V., Kucherova, P., Tentner,
-        A., Lo, S., & Splawski, A. (2007).
-        Advances in the modeling of cladding heat transfer
-        and critical heat flux in boiling water reactor fuel assemblies.
-        In Proc. 12th International Topical Meeting on Nuclear Reactor
-        Thermal Hydraulics (NURETH-12),
-        Pittsburgh, Pennsylvania, USA.
+        Tentner, A., Lo, S., & Kozlov, V. (2006).
+        Advances in computational fluid dynamics modeling
+        of two-phase flow in boiling water reactor fuel assemblies.
+        In International Conference of Nuclear Engineering,
+        Miami, Florida, USA.
     \endverbatim
 
 SourceFiles

applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/twoPhaseCompressibleTurbulenceModels/derivedFvPatchFields/wallBoilingSubModels/partitioningModels/linear/linear.H

@@ -31,10 +31,8 @@ Description
       - alphaLiquid1 0.1
       - alphaLiquid0 0.05
 
-    Reference:
     \verbatim
-        Ioilev, A., Samigulin, M., Ustinenko, V., Kucherova, P., Tentner, A.,
-        Lo, S., & Splawski, A. (2007).
+        Ioilev, A., Samigulin, M., Ustinenko (2007).
         Advances in the modeling of cladding heat transfer
         and critical heat flux in boiling water reactor fuel assemblies.
         In Proc. 12th International Topical Meeting on

applications/test/fileNameClean/Test-fileNameClean.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -65,7 +65,7 @@ int main(int argc, char *argv[])
     argList::noBanner();
     argList::noParallel();
     argList::validArgs.insert("fileName .. fileNameN");
-    argList::addOption("istream", "fileName", "test Istream values");
+    argList::addOption("istream", "file", "test Istream values");
 
     argList args(argc, argv, false, true);
 

applications/utilities/mesh/conversion/fluentMeshToFoam/fluentMeshToFoam.L

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -874,6 +874,12 @@ int main(int argc, char *argv[])
         "factor",
         "geometry scaling factor - default is 1"
     );
+    argList::addOption
+    (
+        "2D",
+        "thickness",
+        "use when converting a 2-D mesh (applied before scale)"
+    );
     argList::addBoolOption
     (
         "writeSets",
@@ -960,22 +966,26 @@ int main(int argc, char *argv[])
     // Construct shapes from face lists
     cellShapeList cellShapes(nCells);
 
-    // Extrude 2-D mesh into 3-D
-
     Info<< "dimension of grid: " << dimensionOfGrid << endl;
     faceList frontAndBackFaces;
 
     if (dimensionOfGrid == 2)
     {
-        const scalar extrusionFactor = 0.01;
+        // Extrude 2-D mesh into 3-D, in z-direction
 
-        boundBox box(max(points), min(points));
+        scalar twoDThickness = 1.0;
 
-        const scalar zOffset = extrusionFactor*box.mag();
+        if (!args.optionReadIfPresent("2D", twoDThickness))
+        {
+            const scalar extrusionFactor = 0.02; //0.01 in each direction
+            boundBox box(points);
+            twoDThickness = extrusionFactor*box.mag();
+        }
 
-        // two-dimensional grid. Extrude in z-direction
-        Info<< "Grid is 2-D. Extruding in z-direction by: "
-            << 2*zOffset << endl;
+        Info<< "Grid is 2-D. Extruding in z-direction by: " << twoDThickness
+            << endl;
+
+        const scalar zOffset = twoDThickness / 2;
 
         pointField oldPoints = points;
 

applications/utilities/mesh/conversion/plot3dToFoam/plot3dToFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -81,7 +81,7 @@ int main(int argc, char *argv[])
     (
         "2D",
         "thickness",
-        "use when converting a 2-D geometry"
+        "use when converting a 2-D mesh (applied before scale)"
     );
 
     argList args(argc, argv);

applications/utilities/mesh/generation/foamyMesh/cellSizeAndAlignmentGrid/Make/options

@@ -29,7 +29,7 @@ EXE_LIBS =  \
     -lfiniteVolume    \
     -lmeshTools \
     -ldecompositionMethods \
-    -L$(FOAM_LIBBIN)/dummy -lptscotchDecomp \
+    -L$(FOAM_LIBBIN)/dummy -lscotchDecomp -lptscotchDecomp \
     -ledgeMesh \
     -ltriSurface \
     -ldynamicMesh \

applications/utilities/mesh/generation/foamyMesh/foamyQuadMesh/Make/options

@@ -35,6 +35,6 @@ EXE_LIBS = \
     -ltriSurface \
     -ldynamicMesh \
     -ldecompositionMethods \
-    -L$(FOAM_LIBBIN)/dummy -lptscotchDecomp \
+    -L$(FOAM_LIBBIN)/dummy -lscotchDecomp -lptscotchDecomp \
     -lsampling \
     -lfileFormats

applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  | Copyright (C) 2015-2016 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -710,7 +710,7 @@ int main(int argc, char *argv[])
     Foam::argList::addOption
     (
         "outFile",
-        "fileName",
+        "file",
         "name of the file to save the simplified surface to"
     );
     #include "addProfilingOption.H"

applications/utilities/mesh/manipulation/renumberMesh/Make/options

@@ -17,4 +17,5 @@ EXE_LIBS = \
     -lrenumberMethods \
     -lreconstruct \
     $(LINK_FLAGS) \
-    -ldecompositionMethods -L$(FOAM_LIBBIN)/dummy -lmetisDecomp -lscotchDecomp
+    -ldecompositionMethods \
+    -L$(FOAM_LIBBIN)/dummy -lmetisDecomp -lscotchDecomp

applications/utilities/miscellaneous/foamDictionary/foamDictionary.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2016-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -396,7 +396,7 @@ int main(int argc, char *argv[])
             );
             if (entPtr)
             {
-                Info<< *entPtr << endl;
+                Info<< *entPtr;
             }
         }
         else if (args.optionFound("remove"))
@@ -467,7 +467,11 @@ int main(int argc, char *argv[])
                             const tokenList& tokens = entPtr->stream();
                             forAll(tokens, i)
                             {
-                                Info<< tokens[i] << token::SPACE;
+                                Info<< tokens[i];
+                                if (i < tokens.size() - 1)
+                                {
+                                    Info<< token::SPACE;
+                                }
                             }
                             Info<< endl;
                         }
@@ -478,7 +482,7 @@ int main(int argc, char *argv[])
                     }
                     else
                     {
-                        Info<< *entPtr << endl;
+                        Info<< *entPtr;
                     }
                 }
             }

bin/stressComponents

@@ -0,0 +1,39 @@
+#!/bin/sh
+#------------------------------------------------------------------------------
+# =========                 |
+# \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+#  \\    /   O peration     |
+#   \\  /    A nd           | Copyright (C) 2016-2017 OpenFOAM Foundation
+#    \\/     M anipulation  |
+#-------------------------------------------------------------------------------
+# License
+#     This file is part of OpenFOAM.
+#
+#     OpenFOAM is free software: you can redistribute it and/or modify it
+#     under the terms of the GNU General Public License as published by
+#     the Free Software Foundation, either version 3 of the License, or
+#     (at your option) any later version.
+#
+#     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+#     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+#     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+#     for more details.
+#
+#     You should have received a copy of the GNU General Public License
+#     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+#
+# Script
+#     stressComponents
+#
+# Description
+#     Script to suggest using the new "-postProcess" solver option.
+#
+#------------------------------------------------------------------------------
+Script=${0##*/}
+
+echo $Script "has been superceded by the -postProcess solver option:"
+echo "<solverName> -funcs '(R components(turbulenceProperties:R))'"
+echo "e.g."
+echo "simpleFoam -postProcess -funcs '(R components(turbulenceProperties:R))'"
+
+#------------------------------------------------------------------------------

bin/wallGradU

@@ -0,0 +1,37 @@
+#!/bin/sh
+#------------------------------------------------------------------------------
+# =========                 |
+# \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+#  \\    /   O peration     |
+#   \\  /    A nd           | Copyright (C) 2017 OpenFOAM Foundation
+#    \\/     M anipulation  |
+#-------------------------------------------------------------------------------
+# License
+#     This file is part of OpenFOAM.
+#
+#     OpenFOAM is free software: you can redistribute it and/or modify it
+#     under the terms of the GNU General Public License as published by
+#     the Free Software Foundation, either version 3 of the License, or
+#     (at your option) any later version.
+#
+#     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+#     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+#     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+#     for more details.
+#
+#     You should have received a copy of the GNU General Public License
+#     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+#
+# Script
+#     wallGradU
+#
+# Description
+#     Script to suggest using the new "postProcess" utility.
+#
+#------------------------------------------------------------------------------
+Script=${0##*/}
+
+echo $Script "has been superceded by the postProcess utility:"
+echo "    postProcess -func 'grad(U)'"
+
+#------------------------------------------------------------------------------

bin/wdot

@@ -0,0 +1,37 @@
+#!/bin/sh
+#------------------------------------------------------------------------------
+# =========                 |
+# \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+#  \\    /   O peration     |
+#   \\  /    A nd           | Copyright (C) 2017 OpenFOAM Foundation
+#    \\/     M anipulation  |
+#-------------------------------------------------------------------------------
+# License
+#     This file is part of OpenFOAM.
+#
+#     OpenFOAM is free software: you can redistribute it and/or modify it
+#     under the terms of the GNU General Public License as published by
+#     the Free Software Foundation, either version 3 of the License, or
+#     (at your option) any later version.
+#
+#     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+#     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+#     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+#     for more details.
+#
+#     You should have received a copy of the GNU General Public License
+#     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+#
+# Script
+#     wdot
+#
+# Description
+#     Script to suggest using the new "postProcess" utility.
+#
+#------------------------------------------------------------------------------
+Script=${0##*/}
+
+echo $Script "has been superceded by the postProcess utility:"
+echo "postProcess -func XiReactionRate"
+
+#------------------------------------------------------------------------------

etc/caseDicts/postProcessing/forces/forceCoeffs.cfg

@@ -9,6 +9,5 @@
 #includeEtc "caseDicts/postProcessing/forces/forces.cfg"
 
 type        forceCoeffs;
-rhoInf      1;           // Redundant for incompressible
 
 // ************************************************************************* //

etc/caseDicts/postProcessing/forces/forceCoeffsCompressible

@@ -25,6 +25,6 @@ dragDir     (1 0 0);
 CofR        (0 0 0);
 pitchAxis   (0 1 0);
 
-#includeEtc "caseDicts/postProcessing/forces/forceCoeffsCompressible.cfg"
+#includeEtc "caseDicts/postProcessing/forces/forceCoeffs.cfg"
 
 // ************************************************************************* //

etc/caseDicts/postProcessing/forces/forceCoeffsIncompressible

@@ -24,6 +24,6 @@ dragDir     (1 0 0);
 CofR        (0 0 0);
 pitchAxis   (0 1 0);
 
-#includeEtc "caseDicts/postProcessing/forces/forceCoeffs.cfg"
+#includeEtc "caseDicts/postProcessing/forces/forceCoeffsIncompressible.cfg"
 
 // ************************************************************************* //

etc/caseDicts/postProcessing/forces/forceCoeffsIncompressible.cfg

@@ -8,4 +8,7 @@
 
 #includeEtc "caseDicts/postProcessing/forces/forceCoeffs.cfg"
 
+rho         rhoInf;
+rhoInf      1;    // Redundant, but currently read in
+
 // ************************************************************************* //

etc/caseDicts/postProcessing/forces/forces.cfg

@@ -12,7 +12,6 @@ libs            ("libforces.so");
 writeControl    timeStep;
 writeInterval   1;
 
-rho             rhoInf; // Incompressible solver
 log             off;
 
 // ************************************************************************* //

etc/caseDicts/postProcessing/forces/forcesCompressible

@@ -16,6 +16,6 @@ patches     (patch1 patch2);
 CofR        (0 0 0);
 pitchAxis   (0 1 0);
 
-#includeEtc "caseDicts/postProcessing/forces/forcesCompressible.cfg"
+#includeEtc "caseDicts/postProcessing/forces/forces.cfg"
 
 // ************************************************************************* //

etc/caseDicts/postProcessing/forces/forcesIncompressible

@@ -11,7 +11,7 @@ Description
 
 \*---------------------------------------------------------------------------*/
 
-#includeEtc "caseDicts/postProcessing/forces/forces.cfg"
+#includeEtc "caseDicts/postProcessing/forces/forcesIncompressible.cfg"
 
 rhoInf      1.225;    // Fluid density
 patches     (patch1 patch2);

etc/caseDicts/postProcessing/forces/forcesIncompressible.cfg

@@ -8,6 +8,6 @@
 
 #includeEtc "caseDicts/postProcessing/forces/forces.cfg"
 
-rhoInf          1;  // Redundant
+rho             rhoInf;
 
 // ************************************************************************* //

src/OpenFOAM/interpolations/interpolation2DTable/interpolation2DTable.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -93,7 +93,7 @@ Foam::interpolation2DTable<Type>::interpolation2DTable(const dictionary& dict)
 :
     List<Tuple2<scalar, List<Tuple2<scalar, Type>>>>(),
     boundsHandling_(wordToBoundsHandling(dict.lookup("outOfBounds"))),
-    fileName_(dict.lookup("fileName")),
+    fileName_(dict.lookup("file")),
     reader_(tableReader<Type>::New(dict))
 {
     readTable();
@@ -445,7 +445,7 @@ void Foam::interpolation2DTable<Type>::checkOrder() const
 template<class Type>
 void Foam::interpolation2DTable<Type>::write(Ostream& os) const
 {
-    os.writeKeyword("fileName")
+    os.writeKeyword("file")
         << fileName_ << token::END_STATEMENT << nl;
     os.writeKeyword("outOfBounds")
         << boundsHandlingToWord(boundsHandling_) << token::END_STATEMENT << nl;

src/OpenFOAM/interpolations/interpolationLookUpTable/interpolationLookUpTable.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -284,7 +284,7 @@ Foam::interpolationLookUpTable<Type>::interpolationLookUpTable
 )
 :
     List<scalarField>(),
-    fileName_(fileName(dict.lookup("fileName")).expand()),
+    fileName_(fileName(dict.lookup("file")).expand()),
     dim_(0),
     min_(0.0),
     delta_(0.0),

src/OpenFOAM/interpolations/interpolationTable/interpolationTable.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -97,7 +97,7 @@ Foam::interpolationTable<Type>::interpolationTable(const dictionary& dict)
 :
     List<Tuple2<scalar, Type>>(),
     boundsHandling_(wordToBoundsHandling(dict.lookup("outOfBounds"))),
-    fileName_(dict.lookup("fileName")),
+    fileName_(dict.lookup("file")),
     reader_(tableReader<Type>::New(dict))
 {
     readTable();
@@ -229,7 +229,7 @@ void Foam::interpolationTable<Type>::check() const
 template<class Type>
 void Foam::interpolationTable<Type>::write(Ostream& os) const
 {
-    os.writeKeyword("fileName")
+    os.writeKeyword("file")
         << fileName_ << token::END_STATEMENT << nl;
     os.writeKeyword("outOfBounds")
         << boundsHandlingToWord(boundsHandling_) << token::END_STATEMENT << nl;

src/OpenFOAM/primitives/functions/Function1/CSV/CSV.C

@@ -2,8 +2,8 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2016-2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -215,7 +215,7 @@ Foam::Function1Types::CSV<Type>::CSV
     componentColumns_(coeffs_.lookup("componentColumns")),
     separator_(coeffs_.lookupOrDefault<string>("separator", string(","))[0]),
     mergeSeparators_(readBool(coeffs_.lookup("mergeSeparators"))),
-    fName_(fName != fileName::null ? fName : coeffs_.lookup("fileName"))
+    fName_(fName != fileName::null ? fName : coeffs_.lookup("file"))
 {
     if (componentColumns_.size() != pTraits<Type>::nComponents)
     {
@@ -283,7 +283,7 @@ void Foam::Function1Types::CSV<Type>::writeData(Ostream& os) const
 
     os.writeEntry("separator", string(separator_));
     os.writeEntry("mergeSeparators", mergeSeparators_);
-    os.writeEntry("fileName",        fName_);
+    os.writeEntry("file", fName_);
 
     os.endBlock() << flush;
 }

src/OpenFOAM/primitives/functions/Function1/TableFile/TableFile.C

@@ -2,8 +2,8 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+     \\/     M anipulation  | Copyright (C) 2016-2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -38,7 +38,7 @@ Foam::Function1Types::TableFile<Type>::TableFile
     fName_("none")
 {
     const dictionary coeffs(dict.subDict(entryName + "Coeffs"));
-    coeffs.lookup("fileName") >> fName_;
+    coeffs.lookup("file") >> fName_;
 
     fileName expandedFile(fName_);
     IFstream is(expandedFile.expand());
@@ -86,7 +86,7 @@ void Foam::Function1Types::TableFile<Type>::writeData(Ostream& os) const
     // the values themselves
     TableBase<Type>::writeEntries(os);
 
-    os.writeEntry("fileName", fName_);
+    os.writeEntry("file", fName_);
 
     os.endBlock() << flush;
 }

src/OpenFOAM/primitives/functions/Polynomial/Polynomial.H

@@ -27,17 +27,19 @@ Class
 Description
     Polynomial templated on size (order):
 
-        poly = sum(coeff_[i]*x^i) logCoeff*log(x)
+    \verbatim
+        poly = sum(coeffs[i]*x^i) + logCoeff*log(x)
+    \endverbatim
 
-    where 0 \<= i \<= N
+    where <tt> 0 <= i <= N </tt>
 
     - integer powers, starting at zero
-    - value(x) to evaluate the poly for a given value
-    - derivative(x) returns derivative at value
-    - integral(x1, x2) returns integral between two scalar values
-    - integral() to return a new, integral coeff polynomial
+    - \c value(x) to evaluate the poly for a given value
+    - \c derivative(x) returns derivative at value
+    - \c integral(x1, x2) returns integral between two scalar values
+    - \c integral() to return a new, integral coeff polynomial
       - increases the size (order)
-    - integralMinus1() to return a new, integral coeff polynomial where
+    - \c integralMinus1() to return a new, integral coeff polynomial where
       the base poly starts at order -1
 
 SourceFiles

src/OpenFOAM/primitives/functions/Polynomial/polynomialFunction.H

@@ -27,16 +27,18 @@ Class
 Description
     Polynomial function representation
 
-        poly = logCoeff*log(x) + sum(coeff_[i]*x^i)
+    \verbatim
+        poly = logCoeff*log(x) + sum(coeffs[i]*x^i)
+    \endverbatim
 
-    where 0 \<= i \<= N
+    where <tt> 0 <= i <= N </tt>
 
     - integer powers, starting at zero
-    - value(x) to evaluate the poly for a given value
-    - integrate(x1, x2) between two scalar values
-    - integral() to return a new, integral coeff polynomial
+    - \c value(x) to evaluate the poly for a given value
+    - \c integrate(x1, x2) between two scalar values
+    - \c integral() to return a new, integral coeff polynomial
       - increases the size (order)
-    - integralMinus1() to return a new, integral coeff polynomial where
+    - \c integralMinus1() to return a new, integral coeff polynomial where
       the base poly starts at order -1
 
 See also

src/combustionModels/PaSR/PaSR.C

@@ -116,29 +116,14 @@ Foam::combustionModels::PaSR<Type>::R(volScalarField& Y) const
 
 template<class Type>
 Foam::tmp<Foam::volScalarField>
-Foam::combustionModels::PaSR<Type>::dQ() const
+Foam::combustionModels::PaSR<Type>::Qdot() const
 {
     return tmp<volScalarField>
     (
         new volScalarField
         (
             IOobject::groupName("PaSR:dQ", this->phaseName_),
-            kappa_*laminar<Type>::dQ()
-        )
-    );
-}
-
-
-template<class Type>
-Foam::tmp<Foam::volScalarField>
-Foam::combustionModels::PaSR<Type>::Sh() const
-{
-    return tmp<volScalarField>
-    (
-        new volScalarField
-        (
-            IOobject::groupName("PaSR:Sh", this->phaseName_),
-            kappa_*laminar<Type>::Sh()
+            kappa_*laminar<Type>::Qdot()
         )
     );
 }

src/combustionModels/PaSR/PaSR.H

@@ -110,11 +110,8 @@ public:
         //- Fuel consumption rate matrix.
         virtual tmp<fvScalarMatrix> R(volScalarField& Y) const;
 
-        //- Heat release rate calculated from fuel consumption rate matrix
-        virtual tmp<volScalarField> dQ() const;
-
-        //-  Return source for enthalpy equation [kg/m/s3]
-        virtual tmp<volScalarField> Sh() const;
+        //- Heat release rate [kg/m/s3]
+        virtual tmp<volScalarField> Qdot() const;
 
         //- Update properties from given dictionary
         virtual bool read();

src/combustionModels/combustionModel/combustionModel.C

@@ -81,7 +81,6 @@ Foam::combustionModel::~combustionModel()
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
-
 bool Foam::combustionModel::read()
 {
     if (regIOobject::read())
@@ -97,25 +96,4 @@ bool Foam::combustionModel::read()
 }
 
 
-Foam::tmp<Foam::volScalarField> Foam::combustionModel::Sh() const
-{
-    return tmp<Foam::volScalarField>
-    (
-        new volScalarField
-        (
-            IOobject
-            (
-                IOobject::groupName("Sh", phaseName_),
-                mesh_.time().timeName(),
-                mesh_,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            mesh_,
-            dimensionedScalar("zero", dimEnergy/dimVolume/dimTime, 0.0)
-        )
-    );
-}
-
-
 // ************************************************************************* //

src/combustionModels/combustionModel/combustionModel.H

@@ -140,11 +140,8 @@ public:
         //- Fuel consumption rate matrix, i.e. source term for fuel equation
         virtual tmp<fvScalarMatrix> R(volScalarField& Y) const = 0;
 
-        //- Heat release rate calculated from fuel consumption rate matrix
-        virtual tmp<volScalarField> dQ() const = 0;
-
-        //-  Return source for enthalpy equation [kg/m/s3]
-        virtual tmp<volScalarField> Sh() const;
+        //- Heat release rate [kg/m/s3]
+        virtual tmp<volScalarField> Qdot() const = 0;
 
         //- Update properties from given dictionary
         virtual bool read();

src/combustionModels/laminar/laminar.C

@@ -136,15 +136,15 @@ Foam::combustionModels::laminar<Type>::R(volScalarField& Y) const
 
 template<class Type>
 Foam::tmp<Foam::volScalarField>
-Foam::combustionModels::laminar<Type>::dQ() const
+Foam::combustionModels::laminar<Type>::Qdot() const
 {
-    tmp<volScalarField> tdQ
+    tmp<volScalarField> tQdot
     (
         new volScalarField
         (
             IOobject
             (
-                IOobject::groupName(typeName + ":dQ", this->phaseName_),
+                IOobject::groupName(typeName + ":Qdot", this->phaseName_),
                 this->mesh().time().timeName(),
                 this->mesh(),
                 IOobject::NO_READ,
@@ -152,47 +152,16 @@ Foam::combustionModels::laminar<Type>::dQ() const
                 false
             ),
             this->mesh(),
-            dimensionedScalar("dQ", dimEnergy/dimTime, 0.0)
+            dimensionedScalar("Qdot", dimEnergy/dimVolume/dimTime, 0.0)
         )
     );
 
     if (this->active())
     {
-        tdQ.ref() = this->chemistryPtr_->dQ();
+        tQdot.ref() = this->chemistryPtr_->Qdot();
     }
 
-    return tdQ;
-}
-
-
-template<class Type>
-Foam::tmp<Foam::volScalarField>
-Foam::combustionModels::laminar<Type>::Sh() const
-{
-    tmp<volScalarField> tSh
-    (
-        new volScalarField
-        (
-            IOobject
-            (
-                IOobject::groupName(typeName + ":Sh", this->phaseName_),
-                this->mesh().time().timeName(),
-                this->mesh(),
-                IOobject::NO_READ,
-                IOobject::NO_WRITE,
-                false
-            ),
-            this->mesh(),
-            dimensionedScalar("zero", dimEnergy/dimTime/dimVolume, 0.0)
-        )
-    );
-
-    if (this->active())
-    {
-        tSh.ref() = this->chemistryPtr_->Sh();
-    }
-
-    return tSh;
+    return tQdot;
 }
 
 

src/combustionModels/laminar/laminar.H

@@ -108,11 +108,8 @@ public:
         //- Fuel consumption rate matrix.
         virtual tmp<fvScalarMatrix> R(volScalarField& Y) const;
 
-        //- Heat release rate calculated from fuel consumption rate matrix
-        virtual tmp<volScalarField> dQ() const;
-
-        //-  Return source for enthalpy equation [kg/m/s3]
-        virtual tmp<volScalarField> Sh() const;
+        //- Heat release rate [kg/m/s3]
+        virtual tmp<volScalarField> Qdot() const;
 
         //- Update properties from given dictionary
         virtual bool read();

src/combustionModels/noCombustion/noCombustion.C

@@ -75,15 +75,15 @@ Foam::combustionModels::noCombustion<CombThermoType>::R
 
 template<class CombThermoType>
 Foam::tmp<Foam::volScalarField>
-Foam::combustionModels::noCombustion<CombThermoType>::dQ() const
+Foam::combustionModels::noCombustion<CombThermoType>::Qdot() const
 {
-    tmp<volScalarField> tdQ
+    tmp<volScalarField> tQdot
     (
         new volScalarField
         (
             IOobject
             (
-                IOobject::groupName("dQ", this->phaseName_),
+                IOobject::groupName("Qdot", this->phaseName_),
                 this->mesh().time().timeName(),
                 this->mesh(),
                 IOobject::NO_READ,
@@ -91,37 +91,11 @@ Foam::combustionModels::noCombustion<CombThermoType>::dQ() const
                 false
             ),
             this->mesh(),
-            dimensionedScalar("dQ", dimEnergy/dimTime, 0.0)
+            dimensionedScalar("Qdot", dimEnergy/dimVolume/dimTime, 0.0)
         )
     );
 
-    return tdQ;
-}
-
-
-template<class CombThermoType>
-Foam::tmp<Foam::volScalarField>
-Foam::combustionModels::noCombustion<CombThermoType>::Sh() const
-{
-    tmp<volScalarField> tSh
-    (
-        new volScalarField
-        (
-            IOobject
-            (
-                IOobject::groupName("Sh", this->phaseName_),
-                this->mesh().time().timeName(),
-                this->mesh(),
-                IOobject::NO_READ,
-                IOobject::NO_WRITE,
-                false
-            ),
-            this->mesh(),
-            dimensionedScalar("zero", dimEnergy/dimTime/dimVolume, 0.0)
-        )
-    );
-
-    return tSh;
+    return tQdot;
 }
 
 

src/combustionModels/noCombustion/noCombustion.H

@@ -92,11 +92,8 @@ public:
         //- Fuel consumption rate matrix
         virtual tmp<fvScalarMatrix> R(volScalarField& Y) const;
 
-        //- Heat release rate calculated from fuel consumption rate matrix
-        virtual tmp<volScalarField> dQ() const;
-
-        //-  Return source for enthalpy equation [kg/m/s3]
-        virtual tmp<volScalarField> Sh() const;
+        //- Heat release rate [kg/m/s3]
+        virtual tmp<volScalarField> Qdot() const;
 
         //- Update properties from given dictionary
         virtual bool read();

src/combustionModels/singleStepCombustion/singleStepCombustion.C

@@ -127,7 +127,7 @@ tmp<fvScalarMatrix> singleStepCombustion<CombThermoType, ThermoType>::R
 
 template<class CombThermoType, class ThermoType>
 tmp<volScalarField>
-singleStepCombustion<CombThermoType, ThermoType>::Sh() const
+singleStepCombustion<CombThermoType, ThermoType>::Qdot() const
 {
     const label fuelI = singleMixturePtr_->fuelIndex();
     volScalarField& YFuel =
@@ -137,37 +137,6 @@ singleStepCombustion<CombThermoType, ThermoType>::Sh() const
 }
 
 
-template<class CombThermoType, class ThermoType>
-tmp<volScalarField>
-singleStepCombustion<CombThermoType, ThermoType>::dQ() const
-{
-    tmp<volScalarField> tdQ
-    (
-        new volScalarField
-        (
-            IOobject
-            (
-                IOobject::groupName("dQ", this->phaseName_),
-                this->mesh_.time().timeName(),
-                this->mesh_,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE,
-                false
-            ),
-            this->mesh_,
-            dimensionedScalar("dQ", dimEnergy/dimTime, 0.0)
-        )
-    );
-
-    if (this->active())
-    {
-        volScalarField& dQ = tdQ.ref();
-        dQ.ref() = this->mesh().V()*Sh()();
-    }
-    return tdQ;
-}
-
-
 template<class CombThermoType, class ThermoType>
 bool singleStepCombustion<CombThermoType, ThermoType>::read()
 {

src/combustionModels/singleStepCombustion/singleStepCombustion.H

@@ -102,11 +102,8 @@ public:
         //- Fuel consumption rate matrix
         virtual tmp<fvScalarMatrix> R(volScalarField& Y) const;
 
-        //- Heat release rate calculated from fuel consumption rate matrix
-        virtual tmp<volScalarField> dQ() const;
-
-        //- Sensible enthalpy source term
-        virtual tmp<volScalarField> Sh() const;
+        //- Heat release rate [kg/m/s3]
+        virtual tmp<volScalarField> Qdot() const;
 
         //- Update properties from given dictionary
         virtual bool read();

src/combustionModels/zoneCombustion/zoneCombustion.C

@@ -165,17 +165,9 @@ Foam::combustionModels::zoneCombustion<Type>::R(volScalarField& Y) const
 
 template<class Type>
 Foam::tmp<Foam::volScalarField>
-Foam::combustionModels::zoneCombustion<Type>::dQ() const
+Foam::combustionModels::zoneCombustion<Type>::Qdot() const
 {
-    return filter(combustionModelPtr_->dQ());
-}
-
-
-template<class Type>
-Foam::tmp<Foam::volScalarField>
-Foam::combustionModels::zoneCombustion<Type>::Sh() const
-{
-    return filter(combustionModelPtr_->Sh());
+    return filter(combustionModelPtr_->Qdot());
 }
 
 

src/combustionModels/zoneCombustion/zoneCombustion.H

@@ -114,11 +114,8 @@ public:
         //- Fuel consumption rate matrix.
         virtual tmp<fvScalarMatrix> R(volScalarField& Y) const;
 
-        //- Heat release rate calculated from fuel consumption rate matrix
-        virtual tmp<volScalarField> dQ() const;
-
-        //-  Return source for enthalpy equation [kg/m/s3]
-        virtual tmp<volScalarField> Sh() const;
+        //- Heat release rate [kg/m/s3]
+        virtual tmp<volScalarField> Qdot() const;
 
         //- Update properties from given dictionary
         virtual bool read();

src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.C

@@ -40,8 +40,9 @@ void Foam::pimpleControl::read()
 {
     solutionControl::read(false);
 
-    // Read solution controls
     const dictionary& pimpleDict = dict();
+
+    solveFlow_ = pimpleDict.lookupOrDefault<Switch>("solveFlow", true);
     nCorrPIMPLE_ = pimpleDict.lookupOrDefault<label>("nOuterCorrectors", 1);
     nCorrPISO_ = pimpleDict.lookupOrDefault<label>("nCorrectors", 1);
     turbOnFinalIterOnly_ =
@@ -123,6 +124,7 @@ bool Foam::pimpleControl::criteriaSatisfied()
 Foam::pimpleControl::pimpleControl(fvMesh& mesh, const word& dictName)
 :
     solutionControl(mesh, dictName),
+    solveFlow_(true),
     nCorrPIMPLE_(0),
     nCorrPISO_(0),
     corrPISO_(0),

src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.H

@@ -69,6 +69,9 @@ protected:
 
         // Solution controls
 
+            //- Flag to indicate whether to solve for the flow
+            bool solveFlow_;
+
             //- Maximum number of PIMPLE correctors
             label nCorrPIMPLE_;
 
@@ -131,22 +134,25 @@ public:
             //- PIMPLE loop
             virtual bool loop();
 
-            //- Pressure corrector loop
+            //- Pressure corrector loop control
             inline bool correct();
 
-            //- Helper function to identify when to store the intial residuals
+            //- Return true to store the intial residuals
             inline bool storeInitialResiduals() const;
 
-            //- Helper function to identify first PIMPLE (outer) iteration
+            //- Return true for first PIMPLE (outer) iteration
             inline bool firstIter() const;
 
-            //- Helper function to identify final PIMPLE (outer) iteration
+            //- Return true fore final PIMPLE (outer) iteration
             inline bool finalIter() const;
 
-            //- Helper function to identify final inner iteration
+            //- Return true for final inner iteration
             inline bool finalInnerIter() const;
 
-            //- Helper function to identify whether to solve for turbulence
+            //- Return true to solve for flow
+            inline bool solveFlow() const;
+
+            //- Return true to solve for turbulence
             inline bool turbCorr() const;
 };
 

src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControlI.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -92,6 +92,12 @@ inline bool Foam::pimpleControl::finalInnerIter() const
 }
 
 
+inline bool Foam::pimpleControl::solveFlow() const
+{
+    return solveFlow_;
+}
+
+
 inline bool Foam::pimpleControl::turbCorr() const
 {
     return !turbOnFinalIterOnly_ || finalIter();

src/finiteVolume/fields/fvPatchFields/derived/prghTotalPressure/prghTotalPressureFvPatchScalarField.H

@@ -33,6 +33,9 @@ Description
 
         \f[
             p_rgh = p - \rho g.(h - hRef)
+        \f]
+
+        \f[
             p = p0 - 0.5 \rho |U|^2
         \f]
 

src/finiteVolume/fvMesh/fvMesh.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -267,7 +267,7 @@ Foam::fvMesh::fvMesh(const IOobject& io)
         InfoInFunction << "Constructing fvMesh from IOobject" << endl;
     }
 
-    // Check the existance of the cell volumes and read if present
+    // Check the existence of the cell volumes and read if present
     // and set the storage of V00
     if (isFile(time().timePath()/"V0"))
     {
@@ -288,7 +288,7 @@ Foam::fvMesh::fvMesh(const IOobject& io)
         V00();
     }
 
-    // Check the existance of the mesh fluxes, read if present and set the
+    // Check the existence of the mesh fluxes, read if present and set the
     // mesh to be moving
     if (isFile(time().timePath()/"meshPhi"))
     {

src/functionObjects/forces/forces/forces.C

@@ -870,7 +870,7 @@ bool Foam::functionObjects::forces::read(const dictionary& dict)
         // Reference density needed for incompressible calculations
         if (rhoName_ == "rhoInf")
         {
-            rhoRef_ = readScalar(dict.lookup("rhoInf"));
+            dict.lookup("rhoInf") >> rhoRef_;
         }
 
         // Reference pressure, 0 by default

src/functionObjects/utilities/abort/abort.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -125,7 +125,7 @@ bool Foam::functionObjects::abort::read(const dictionary& dict)
         action_ = nextWrite;
     }
 
-    if (dict.readIfPresent("fileName", abortFile_))
+    if (dict.readIfPresent("file", abortFile_))
     {
         abortFile_.expand();
     }

src/fvOptions/sources/derived/effectivenessHeatExchangerSource/effectivenessHeatExchangerSource.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2017 OpenFOAM Foundation
      \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -124,7 +124,7 @@ Usage
     \endverbatim
 
 Note
-- the table with name "fileName" should have the same units as the
+- the table with name "file" should have the same units as the
   secondary mass flow rate and kg/s for phi
 - faceZone is the faces at the inlet of the cellzone, it needs to be
   created with flip map flags. It is used to integrate the net mass flow

src/fvOptions/sources/derived/rotorDiskSource/profileModel/profileModel.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -51,7 +51,7 @@ Foam::profileModel::profileModel(const dictionary& dict, const word& name)
     name_(name),
     fName_(fileName::null)
 {
-    dict.readIfPresent("fileName", fName_);
+    dict.readIfPresent("file", fName_);
 }
 
 // * * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * //

src/meshTools/searchableSurface/triSurfaceMesh.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -35,83 +35,59 @@ License
 
 namespace Foam
 {
-
     defineTypeNameAndDebug(triSurfaceMesh, 0);
     addToRunTimeSelectionTable(searchableSurface, triSurfaceMesh, dict);
-
     word triSurfaceMesh::meshSubDir = "triSurface";
-
 }
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-//// Special version of Time::findInstance that does not check headerOk
-//// to search for instances of raw files
-//Foam::word Foam::triSurfaceMesh::findRawInstance
-//(
-//    const Time& runTime,
-//    const fileName& dir,
-//    const word& name
-//)
-//{
-//    // Check current time first
-//    if (isFile(runTime.path()/runTime.timeName()/dir/name))
-//    {
-//        return runTime.timeName();
-//    }
-//    instantList ts = runTime.times();
-//    label instanceI;
-//
-//    for (instanceI = ts.size()-1; instanceI >= 0; --instanceI)
-//    {
-//        if (ts[instanceI].value() <= runTime.timeOutputValue())
-//        {
-//            break;
-//        }
-//    }
-//
-//    // continue searching from here
-//    for (; instanceI >= 0; --instanceI)
-//    {
-//        if (isFile(runTime.path()/ts[instanceI].name()/dir/name))
-//        {
-//            return ts[instanceI].name();
-//        }
-//    }
-//
-//
-//    // not in any of the time directories, try constant
-//
-//    // Note. This needs to be a hard-coded constant, rather than the
-//    // constant function of the time, because the latter points to
-//    // the case constant directory in parallel cases
-//
-//    if (isFile(runTime.path()/runTime.constant()/dir/name))
-//    {
-//        return runTime.constant();
-//    }
-//
-//    FatalErrorInFunction
-//        << "Cannot find file \"" << name << "\" in directory "
-//        << runTime.constant()/dir
-//        << exit(FatalError);
-//
-//    return runTime.constant();
-//}
-
-
-const Foam::fileName& Foam::triSurfaceMesh::checkFile
-(
-    const fileName& fName,
-    const fileName& objectName
-)
+Foam::fileName Foam::triSurfaceMesh::checkFile(const IOobject& io)
 {
+    const fileName fName(io.filePath());
     if (fName.empty())
     {
         FatalErrorInFunction
             << "Cannot find triSurfaceMesh starting from "
-            << objectName << exit(FatalError);
+            << io.objectPath() << exit(FatalError);
     }
+
+    return fName;
+}
+
+
+Foam::fileName Foam::triSurfaceMesh::checkFile
+(
+    const IOobject& io,
+    const dictionary& dict
+)
+{
+    fileName fName;
+    if (dict.readIfPresent("file", fName, false, false))
+    {
+        fName.expand();
+        if (!fName.isAbsolute())
+        {
+            fName = io.objectPath().path()/fName;
+        }
+        if (!exists(fName))
+        {
+            FatalErrorInFunction
+                << "Cannot find triSurfaceMesh at " << fName
+                << exit(FatalError);
+        }
+    }
+    else
+    {
+        fName = io.filePath();
+        if (!exists(fName))
+        {
+            FatalErrorInFunction
+                << "Cannot find triSurfaceMesh starting from "
+                << io.objectPath() << exit(FatalError);
+        }
+    }
+
     return fName;
 }
 
@@ -227,7 +203,7 @@ Foam::triSurfaceMesh::triSurfaceMesh(const IOobject& io, const triSurface& s)
         (
             io.name(),
             io.instance(),
-            io.local(),  //"triSurfaceFields",
+            io.local(),
             io.db(),
             io.readOpt(),
             io.writeOpt(),
@@ -235,7 +211,7 @@ Foam::triSurfaceMesh::triSurfaceMesh(const IOobject& io, const triSurface& s)
         )
     ),
     triSurface(s),
-    triSurfaceRegionSearch(static_cast<const triSurface&>(*this)),
+    triSurfaceRegionSearch(s),
     minQuality_(-1),
     surfaceClosed_(-1),
     outsideVolType_(volumeType::UNKNOWN)
@@ -250,43 +226,21 @@ Foam::triSurfaceMesh::triSurfaceMesh(const IOobject& io)
 :
     // Find instance for triSurfaceMesh
     searchableSurface(io),
-    //(
-    //    IOobject
-    //    (
-    //        io.name(),
-    //        io.time().findInstance(io.local(), word::null),
-    //        io.local(),
-    //        io.db(),
-    //        io.readOpt(),
-    //        io.writeOpt(),
-    //        io.registerObject()
-    //    )
-    //),
     // Reused found instance in objectRegistry
     objectRegistry
     (
         IOobject
         (
             io.name(),
-            static_cast<const searchableSurface&>(*this).instance(),
-            io.local(), //"triSurfaceFields",
+            searchableSurface::instance(),
+            io.local(),
             io.db(),
             io.readOpt(),
             io.writeOpt(),
             false       // searchableSurface already registered under name
         )
     ),
-    triSurface
-    (
-        checkFile
-        (
-            typeFilePath<searchableSurface>
-            (
-                static_cast<const searchableSurface&>(*this)
-            ),
-            searchableSurface::objectPath()
-        )
-    ),
+    triSurface(checkFile(static_cast<const searchableSurface&>(*this))),
     triSurfaceRegionSearch(static_cast<const triSurface&>(*this)),
     minQuality_(-1),
     surfaceClosed_(-1),
@@ -305,51 +259,32 @@ Foam::triSurfaceMesh::triSurfaceMesh
 )
 :
     searchableSurface(io),
-    //(
-    //    IOobject
-    //    (
-    //        io.name(),
-    //        io.time().findInstance(io.local(), word::null),
-    //        io.local(),
-    //        io.db(),
-    //        io.readOpt(),
-    //        io.writeOpt(),
-    //        io.registerObject()
-    //    )
-    //),
     // Reused found instance in objectRegistry
     objectRegistry
     (
         IOobject
         (
             io.name(),
-            static_cast<const searchableSurface&>(*this).instance(),
-            io.local(), //"triSurfaceFields",
+            searchableSurface::instance(),
+            io.local(),
             io.db(),
             io.readOpt(),
             io.writeOpt(),
             false       // searchableSurface already registered under name
         )
     ),
-    triSurface
-    (
-        checkFile
-        (
-            typeFilePath<searchableSurface>
-            (
-                static_cast<const searchableSurface&>(*this)
-            ),
-            searchableSurface::objectPath()
-        )
-    ),
+    triSurface(checkFile(static_cast<const searchableSurface&>(*this), dict)),
     triSurfaceRegionSearch(static_cast<const triSurface&>(*this), dict),
     minQuality_(-1),
     surfaceClosed_(-1),
     outsideVolType_(volumeType::UNKNOWN)
 {
+    // Reading from supplied file name instead of objectPath/filePath
+    dict.readIfPresent("file", fName_, false, false);
+
     scalar scaleFactor = 0;
 
-    // allow rescaling of the surface points
+    // Allow rescaling of the surface points
     // eg, CAD geometries are often done in millimeters
     if (dict.readIfPresent("scale", scaleFactor) && scaleFactor > 0)
     {
@@ -495,7 +430,7 @@ Foam::triSurfaceMesh::edgeTree() const
             label nPoints;
             PatchTools::calcBounds
             (
-                static_cast<const triSurface&>(*this),
+                *this,
                 bb,
                 nPoints
             );
@@ -552,7 +487,6 @@ const Foam::wordList& Foam::triSurfaceMesh::regions() const
 }
 
 
-// Find out if surface is closed.
 bool Foam::triSurfaceMesh::hasVolumeType() const
 {
     if (surfaceClosed_ == -1)
@@ -660,7 +594,7 @@ void Foam::triSurfaceMesh::getNormal
     vectorField& normal
 ) const
 {
-    const triSurface& s = static_cast<const triSurface&>(*this);
+    const triSurface& s = *this;
     const pointField& pts = s.points();
 
     normal.setSize(info.size());
@@ -865,7 +799,24 @@ bool Foam::triSurfaceMesh::writeObject
             runTime.timeName();
     }
 
-    fileName fullPath(searchableSurface::objectPath());
+    fileName fullPath;
+    if (fName_.size())
+    {
+        // Override file name
+
+        fullPath = fName_;
+
+        fullPath.expand();
+        if (!fullPath.isAbsolute())
+        {
+            // Add directory from regIOobject
+            fullPath = searchableSurface::objectPath().path()/fullPath;
+        }
+    }
+    else
+    {
+        fullPath = searchableSurface::objectPath();
+    }
 
     if (!mkDir(fullPath.path()))
     {

src/meshTools/searchableSurface/triSurfaceMesh.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -68,10 +68,11 @@ class triSurfaceMesh
     public triSurface,
     public triSurfaceRegionSearch
 {
-private:
-
     // Private member data
 
+        //- Supplied fileName override
+        fileName fName_;
+
         //- Optional min triangle quality. Triangles below this get
         //  ignored for normal calculation
         scalar minQuality_;
@@ -91,20 +92,11 @@ private:
 
     // Private Member Functions
 
-        ////- Helper: find instance of files without header
-        //static word findRawInstance
-        //(
-        //    const Time&,
-        //    const fileName&,
-        //    const word&
-        //);
+        //- Return fileName to load IOobject from
+        static fileName checkFile(const IOobject& io);
 
-        //- Check file existence
-        static const fileName& checkFile
-        (
-            const fileName& fName,
-            const fileName& objectName
-        );
+        //- Return fileName to load IOobject from. Optional override of fileName
+        static fileName checkFile(const IOobject&, const dictionary&);
 
         //- Helper function for isSurfaceClosed
         static bool addFaceToEdge

src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.C

@@ -318,7 +318,7 @@ void reactingOneDim::solveEnergy()
       + fvc::laplacian(alpha, h_)
       - fvc::laplacian(kappa(), T())
      ==
-        chemistrySh_
+        chemistryQdot_
       + solidChemistry_->RRsHs()
     );
 
@@ -371,7 +371,7 @@ void reactingOneDim::calculateMassTransfer()
 
     if (infoOutput_)
     {
-        totalHeatRR_ = fvc::domainIntegrate(chemistrySh_);
+        totalHeatRR_ = fvc::domainIntegrate(chemistryQdot_);
 
         addedGasMass_ +=
             fvc::domainIntegrate(solidChemistry_->RRg())*time_.deltaT();
@@ -440,11 +440,11 @@ reactingOneDim::reactingOneDim
         dimensionedScalar("zero", dimEnergy/dimTime, 0.0)
     ),
 
-    chemistrySh_
+    chemistryQdot_
     (
         IOobject
         (
-            "chemistrySh",
+            "chemistryQdot",
             time().timeName(),
             regionMesh(),
             IOobject::NO_READ,
@@ -540,11 +540,11 @@ reactingOneDim::reactingOneDim
         dimensionedScalar("zero", dimEnergy/dimTime, 0.0)
     ),
 
-    chemistrySh_
+    chemistryQdot_
     (
         IOobject
         (
-            "chemistrySh",
+            "chemistryQdot",
             time().timeName(),
             regionMesh(),
             IOobject::NO_READ,
@@ -699,7 +699,7 @@ void reactingOneDim::evolveRegion()
 
     solveContinuity();
 
-    chemistrySh_ = solidChemistry_->Sh()();
+    chemistryQdot_ = solidChemistry_->Qdot()();
 
     updateFields();
 

src/regionModels/pyrolysisModels/reactingOneDim/reactingOneDim.H

@@ -119,8 +119,8 @@ protected:
             //- Sensible enthalpy gas flux [J/m2/s]
             volScalarField phiHsGas_;
 
-            //- Heat release [J/s/m3]
-            volScalarField chemistrySh_;
+            //- Heat release rate [J/s/m3]
+            volScalarField chemistryQdot_;
 
 
         // Source term fields

src/thermophysicalModels/chemistryModel/Make/options

@@ -7,7 +7,8 @@ EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/thermophysicalFunctions/lnInclude \
     -I$(LIB_SRC)/turbulenceModels/compressible/lnInclude \
     -I$(LIB_SRC)/ODE/lnInclude \
-    -I$(LIB_SRC)/finiteVolume/lnInclude
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude
 
 LIB_LIBS = \
     -lcompressibleTransportModels \
@@ -16,4 +17,5 @@ LIB_LIBS = \
     -lspecie \
     -lthermophysicalFunctions \
     -lODE \
-    -lfiniteVolume
+    -lfiniteVolume \
+    -lmeshTools

src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModel.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2016-2017 OpenFOAM Foundation
      \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -25,6 +25,7 @@ License
 
 #include "TDACChemistryModel.H"
 #include "UniformField.H"
+#include "localEulerDdtScheme.H"
 #include "clockTime.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
@@ -37,12 +38,31 @@ Foam::TDACChemistryModel<CompType, ThermoType>::TDACChemistryModel
 )
 :
     chemistryModel<CompType, ThermoType>(mesh, phaseName),
+    variableTimeStep_
+    (
+        mesh.time().controlDict().lookupOrDefault("adjustTimeStep", false)
+     || fv::localEulerDdt::enabled(mesh)
+    ),
+    timeSteps_(0),
     NsDAC_(this->nSpecie_),
-    completeC_(this->nSpecie_,0.0),
+    completeC_(this->nSpecie_, 0),
     reactionsDisabled_(this->reactions_.size(), false),
+    specieComp_(this->nSpecie_),
     completeToSimplifiedIndex_(this->nSpecie_, -1),
     simplifiedToCompleteIndex_(this->nSpecie_),
-    specieComp_(this->nSpecie_)
+    tabulationResults_
+    (
+        IOobject
+        (
+            "TabulationResults",
+            this->time().timeName(),
+            this->mesh(),
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        scalar(0)
+    )
 {
     basicMultiComponentMixture& composition = this->thermo().composition();
 
@@ -579,8 +599,13 @@ Foam::scalar Foam::TDACChemistryModel<CompType, ThermoType>::solve
     const DeltaTType& deltaT
 )
 {
+    // Increment counter of time-step
+    timeSteps_++;
+
     const bool reduced = mechRed_->active();
 
+    label nAdditionalEqn = (tabulation_->variableTimeStep() ? 1 : 0);
+
     basicMultiComponentMixture& composition = this->thermo().composition();
 
     // CPU time analysis
@@ -625,9 +650,9 @@ Foam::scalar Foam::TDACChemistryModel<CompType, ThermoType>::solve
     scalarField c0(this->nSpecie_);
 
     // Composition vector (Yi, T, p)
-    scalarField phiq(this->nEqns());
+    scalarField phiq(this->nEqns() + nAdditionalEqn);
 
-    scalarField Rphiq(this->nEqns());
+    scalarField Rphiq(this->nEqns() + nAdditionalEqn);
 
     forAll(rho, celli)
     {
@@ -643,6 +668,11 @@ Foam::scalar Foam::TDACChemistryModel<CompType, ThermoType>::solve
         }
         phiq[this->nSpecie()]=Ti;
         phiq[this->nSpecie() + 1]=pi;
+        if (tabulation_->variableTimeStep())
+        {
+            phiq[this->nSpecie() + 2] = deltaT[celli];
+        }
+
 
         // Initialise time progress
         scalar timeLeft = deltaT[celli];
@@ -668,7 +698,7 @@ Foam::scalar Foam::TDACChemistryModel<CompType, ThermoType>::solve
         // This position is reached when tabulation is not used OR
         // if the solution is not retrieved.
         // In the latter case, it adds the information to the tabulation
-        // (it will either expand the current data or add a new stored poin).
+        // (it will either expand the current data or add a new stored point).
         else
         {
             clockTime_.timeIncrement();
@@ -720,12 +750,28 @@ Foam::scalar Foam::TDACChemistryModel<CompType, ThermoType>::solve
                 {
                     Rphiq[i] = c[i]/rhoi*this->specieThermo_[i].W();
                 }
-
+                if (tabulation_->variableTimeStep())
+                {
+                    Rphiq[Rphiq.size()-3] = Ti;
+                    Rphiq[Rphiq.size()-2] = pi;
+                    Rphiq[Rphiq.size()-1] = deltaT[celli];
+                }
+                else
+                {
                     Rphiq[Rphiq.size()-2] = Ti;
                     Rphiq[Rphiq.size()-1] = pi;
-                tabulation_->add(phiq, Rphiq, rhoi);
                 }
-
+                label growOrAdd =
+                    tabulation_->add(phiq, Rphiq, rhoi, deltaT[celli]);
+                if (growOrAdd)
+                {
+                    this->setTabulationResultsAdd(celli);
+                }
+                else
+                {
+                    this->setTabulationResultsGrow(celli);
+                }
+            }
             addNewLeafCpuTime_ += clockTime_.timeIncrement();
 
             // When operations are done and if mechanism reduction is active,
@@ -840,4 +886,35 @@ Foam::scalar Foam::TDACChemistryModel<CompType, ThermoType>::solve
 }
 
 
+template<class CompType, class ThermoType>
+void Foam::TDACChemistryModel<CompType, ThermoType>::setTabulationResultsAdd
+(
+    const label celli
+)
+{
+    tabulationResults_[celli] = 0.0;
+}
+
+
+template<class CompType, class ThermoType>
+void Foam::TDACChemistryModel<CompType, ThermoType>::setTabulationResultsGrow
+(
+    const label celli
+)
+{
+    tabulationResults_[celli] = 1.0;
+}
+
+
+template<class CompType, class ThermoType>
+void Foam::TDACChemistryModel<CompType, ThermoType>::
+setTabulationResultsRetrieve
+(
+    const label celli
+)
+{
+    tabulationResults_[celli] = 2.0;
+}
+
+
 // ************************************************************************* //

src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModel.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2016-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -85,14 +85,18 @@ class TDACChemistryModel
 {
     // Private member data
 
+        bool variableTimeStep_;
+
+        label timeSteps_;
+
         // Mechanism reduction
         label NsDAC_;
         scalarField completeC_;
         scalarField simplifiedC_;
         Field<bool> reactionsDisabled_;
+        List<List<specieElement>> specieComp_;
         Field<label> completeToSimplifiedIndex_;
         DynamicList<label> simplifiedToCompleteIndex_;
-        List<List<specieElement>> specieComp_;
         autoPtr<chemistryReductionMethod<CompType, ThermoType>> mechRed_;
 
         // Tabulation
@@ -113,6 +117,12 @@ class TDACChemistryModel
         //- Log file for average time spent solving the chemistry
         autoPtr<OFstream> cpuSolveFile_;
 
+        // Field containing information about tabulation:
+        // 0 -> add (direct integration)
+        // 1 -> grow
+        // 2 -> retrieve
+        volScalarField tabulationResults_;
+
 
     // Private Member Functions
 
@@ -151,6 +161,12 @@ public:
 
     // Member Functions
 
+        //- Return true if the time-step is variable and/or non-uniform
+        inline bool variableTimeStep() const;
+
+        //- Return the number of chemistry evaluations (used by ISAT)
+        inline label timeSteps() const;
+
         //- Create and return a TDAC log file of the given name
         inline autoPtr<OFstream> logFile(const word& name) const;
 
@@ -256,6 +272,17 @@ public:
 
             inline autoPtr<chemistryReductionMethod<CompType, ThermoType>>&
                 mechRed();
+
+            tmp<volScalarField> tabulationResults() const
+            {
+                return tabulationResults_;
+            }
+
+            void setTabulationResultsAdd(const label celli);
+
+            void setTabulationResultsGrow(const label celli);
+
+            void setTabulationResultsRetrieve(const label celli);
 };
 
 

src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModelI.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2016-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,6 +25,22 @@ License
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+template<class CompType, class ThermoType>
+inline bool
+Foam::TDACChemistryModel<CompType, ThermoType>::variableTimeStep() const
+{
+    return variableTimeStep_;
+}
+
+
+template<class CompType, class ThermoType>
+inline Foam::label
+Foam::TDACChemistryModel<CompType, ThermoType>::timeSteps() const
+{
+    return timeSteps_;
+}
+
+
 template<class CompType, class ThermoType>
 inline Foam::autoPtr<Foam::OFstream>
 Foam::TDACChemistryModel<CompType, ThermoType>::logFile(const word& name) const

src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/ISAT/ISAT.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2016-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -41,16 +41,11 @@ Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::ISAT
         chemistry
     ),
     chemisTree_(chemistry, this->coeffsDict_),
-    scaleFactor_(chemistry.nEqns(), 1.0),
+    scaleFactor_(chemistry.nEqns() + ((this->variableTimeStep()) ? 1 : 0), 1),
     runTime_(chemistry.time()),
     chPMaxLifeTime_
     (
-        this->coeffsDict_.lookupOrDefault
-        (
-            "chPMaxLifeTime",
-            (runTime_.endTime().value()-runTime_.startTime().value())
-           /runTime_.deltaT().value()
-        )
+        this->coeffsDict_.lookupOrDefault("chPMaxLifeTime", INT_MAX)
     ),
     maxGrowth_(this->coeffsDict_.lookupOrDefault("maxGrowth", INT_MAX)),
     checkEntireTreeInterval_
@@ -104,6 +99,19 @@ Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::ISAT
         }
         scaleFactor_[Ysize] = readScalar(scaleDict.lookup("Temperature"));
         scaleFactor_[Ysize + 1] = readScalar(scaleDict.lookup("Pressure"));
+        if (this->variableTimeStep())
+        {
+            scaleFactor_[Ysize + 2] = readScalar(scaleDict.lookup("deltaT"));
+        }
+    }
+
+    if (this->variableTimeStep())
+    {
+        nAdditionalEqns_ = 3;
+    }
+    else
+    {
+        nAdditionalEqns_ = 2;
     }
 
     if (this->log())
@@ -190,7 +198,7 @@ void Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::calcNewC
     scalarField& Rphiq
 )
 {
-    label nEqns = this->chemistry_.nEqns(); // Full set of species
+    label nEqns = this->chemistry_.nEqns(); // Species, T, p
     bool mechRedActive = this->chemistry_.mechRed()->active();
     Rphiq = phi0->Rphi();
     scalarField dphi(phiq-phi0->phi());
@@ -199,7 +207,7 @@ void Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::calcNewC
 
     // Rphiq[i]=Rphi0[i]+A(i, j)dphi[j]
     // where Aij is dRi/dphi_j
-    for (label i=0; i<nEqns-2; i++)
+    for (label i=0; i<nEqns-nAdditionalEqns_; i++)
     {
         if (mechRedActive)
         {
@@ -218,7 +226,16 @@ void Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::calcNewC
                 Rphiq[i] +=
                     gradientsMatrix(si, phi0->nActiveSpecies())*dphi[nEqns - 2];
                 Rphiq[i] +=
-                    gradientsMatrix(si, phi0->nActiveSpecies()+1)*dphi[nEqns-1];
+                    gradientsMatrix(si, phi0->nActiveSpecies() + 1)
+                   *dphi[nEqns - 1];
+
+                if (this->variableTimeStep())
+                {
+                    Rphiq[i] +=
+                        gradientsMatrix(si, phi0->nActiveSpecies() + 2)
+                       *dphi[nEqns];
+                }
+
                 // As we use an approximation of A, Rphiq should be checked for
                 // negative values
                 Rphiq[i] = max(0.0,Rphiq[i]);
@@ -260,10 +277,11 @@ bool Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::grow
 
     // Raise a flag when the chemPoint used has been grown more than the
     // allowed number of time
-    if (!phi0->toRemove() && phi0->nGrowth() > maxGrowth_)
+    if (phi0->nGrowth() > maxGrowth_)
     {
         cleaningRequired_ = true;
         phi0->toRemove() = true;
+        return false;
     }
 
     // If the solution RphiQ is still within the tolerance we try to grow it
@@ -294,14 +312,10 @@ Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::cleanAndBalance()
     {
         chemPointISAT<CompType, ThermoType>* xtmp =
             chemisTree_.treeSuccessor(x);
-        // timeOutputValue returns timeToUserTime(value()), therefore, it should
-        // be compare with timeToUserTime(deltaT)
-        scalar elapsedTime = runTime_.timeOutputValue() - x->timeTag();
-        scalar maxElapsedTime =
-            chPMaxLifeTime_
-          * runTime_.timeToUserTime(runTime_.deltaTValue());
 
-        if ((elapsedTime > maxElapsedTime) || (x->nGrowth() > maxGrowth_))
+        scalar elapsedTimeSteps = this->chemistry_.timeSteps() - x->timeTag();
+
+        if ((elapsedTimeSteps > chPMaxLifeTime_) || (x->nGrowth() > maxGrowth_))
         {
             chemisTree_.deleteLeaf(x);
             treeModified = true;
@@ -334,13 +348,13 @@ void Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::computeA
 (
     scalarSquareMatrix& A,
     const scalarField& Rphiq,
-    const scalar rhoi
+    const scalar rhoi,
+    const scalar dt
 )
 {
-    scalar dt = runTime_.deltaTValue();
     bool mechRedActive = this->chemistry_.mechRed()->active();
     label speciesNumber = this->chemistry_.nSpecie();
-    scalarField Rcq(this->chemistry_.nEqns());
+    scalarField Rcq(this->chemistry_.nEqns() + nAdditionalEqns_ - 2);
     for (label i=0; i<speciesNumber; i++)
     {
         label s2c = i;
@@ -350,8 +364,12 @@ void Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::computeA
         }
         Rcq[i] = rhoi*Rphiq[s2c]/this->chemistry_.specieThermo()[s2c].W();
     }
-    Rcq[speciesNumber] = Rphiq[Rphiq.size()-2];
-    Rcq[speciesNumber+1] = Rphiq[Rphiq.size()-1];
+    Rcq[speciesNumber] = Rphiq[Rphiq.size() - nAdditionalEqns_];
+    Rcq[speciesNumber+1] = Rphiq[Rphiq.size() - nAdditionalEqns_ + 1];
+    if (this->variableTimeStep())
+    {
+        Rcq[speciesNumber + 2] = Rphiq[Rphiq.size() - nAdditionalEqns_ + 2];
+    }
 
     // Aaa is computed implicitely,
     // A is given by A = C(psi0, t0+dt), where C is obtained through solving
@@ -399,6 +417,10 @@ void Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::computeA
     // For temperature and pressure, only unity on the diagonal
     A(speciesNumber, speciesNumber) = 1;
     A(speciesNumber + 1, speciesNumber + 1) = 1;
+    if (this->variableTimeStep())
+    {
+        A[speciesNumber + 2][speciesNumber + 2] = 1;
+    }
 
     // Inverse of (I-dt*J(psi(t0+dt)))
     LUscalarMatrix LUA(A);
@@ -463,20 +485,18 @@ bool Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::retrieve
 
     if (retrieved)
     {
-        scalar elapsedTime =
-            runTime_.timeOutputValue() - phi0->timeTag();
-        scalar maxElapsedTime =
-            chPMaxLifeTime_
-          * runTime_.timeToUserTime(runTime_.deltaTValue());
+        phi0->increaseNumRetrieve();
+        scalar elapsedTimeSteps =
+            this->chemistry_.timeSteps() - phi0->timeTag();
 
         // Raise a flag when the chemPoint has been used more than the allowed
         // number of time steps
-        if (elapsedTime > maxElapsedTime && !phi0->toRemove())
+        if (elapsedTimeSteps > chPMaxLifeTime_ && !phi0->toRemove())
         {
             cleaningRequired_ = true;
             phi0->toRemove() = true;
         }
-        lastSearch_->lastTimeUsed() = runTime_.timeOutputValue();
+        lastSearch_->lastTimeUsed() = this->chemistry_.timeSteps();
         addToMRU(phi0);
         calcNewC(phi0,phiq, Rphiq);
         nRetrieved_++;
@@ -492,13 +512,15 @@ bool Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::retrieve
 
 
 template<class CompType, class ThermoType>
-bool Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::add
+Foam::label Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::add
 (
     const scalarField& phiq,
     const scalarField& Rphiq,
-    const scalar rho
+    const scalar rho,
+    const scalar deltaT
 )
 {
+    label growthOrAddFlag = 1;
     // If lastSearch_ holds a valid pointer to a chemPoint AND the growPoints_
     // option is on, the code first tries to grow the point hold by lastSearch_
     if (lastSearch_ && growPoints_)
@@ -506,13 +528,12 @@ bool Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::add
         if (grow(lastSearch_,phiq, Rphiq))
         {
             nGrowth_++;
+            growthOrAddFlag = 0;
             //the structure of the tree is not modified, return false
-            return false;
+            return growthOrAddFlag;
         }
     }
 
-    bool treeCleanedOrCleared(false);
-
     // If the code reach this point, it is either because lastSearch_ is not
     // valid, OR because growPoints_ is not on, OR because the grow operation
     // has failed. In the three cases, a new point is added to the tree.
@@ -567,16 +588,12 @@ bool Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::add
         // The structure has been changed, it will force the binary tree to
         // perform a new search and find the most appropriate point still stored
         lastSearch_ = nullptr;
-
-        // Either cleanAndBalance has changed the tree or it has been cleared
-        // in any case treeCleanedOrCleared should be set to true
-        treeCleanedOrCleared = true;
     }
 
     // Compute the A matrix needed to store the chemPoint.
-    label ASize = this->chemistry_.nEqns(); // Reduced when mechRed is active
+    label ASize = this->chemistry_.nEqns() + nAdditionalEqns_ - 2;
     scalarSquareMatrix A(ASize, Zero);
-    computeA(A, Rphiq, rho);
+    computeA(A, Rphiq, rho, deltaT);
 
     chemisTree().insertNewLeaf
     (
@@ -591,7 +608,7 @@ bool Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::add
 
     nAdd_++;
 
-    return treeCleanedOrCleared;
+    return growthOrAddFlag;
 }
 
 

src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/ISAT/ISAT.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2016-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -112,6 +112,9 @@ class ISAT
 
         bool cleaningRequired_;
 
+        //- Number of equations in addition to the species eqs.
+        label nAdditionalEqns_;
+
 
     // Private Member Functions
 
@@ -168,7 +171,8 @@ class ISAT
         (
             scalarSquareMatrix& A,
             const scalarField& Rphiq,
-            const scalar rho
+            const scalar rho,
+            const scalar dt
         );
 
 
@@ -224,11 +228,12 @@ public:
         //  This function can grow an existing point or add a new leaf to the
         //  binary tree Input : phiq the new composition to store Rphiq the
         //  mapping of the new composition point
-        virtual bool add
+        virtual label add
         (
             const scalarField& phiq,
             const scalarField& Rphiq,
-            const scalar rho
+            const scalar rho,
+            const scalar deltaT
         );
 
         virtual bool update()

src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/ISAT/binaryNode/binaryNode.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2016-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -28,15 +28,14 @@ License
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class CompType, class ThermoType>
-Foam::binaryNode<CompType, ThermoType>::binaryNode
-(
-)
+Foam::binaryNode<CompType, ThermoType>::binaryNode()
 :
     leafLeft_(nullptr),
     leafRight_(nullptr),
     nodeLeft_(nullptr),
     nodeRight_(nullptr),
-    parent_(nullptr)
+    parent_(nullptr),
+    nAdditionalEqns_(0)
 {}
 
 
@@ -53,34 +52,26 @@ Foam::binaryNode<CompType, ThermoType>::binaryNode
     nodeLeft_(nullptr),
     nodeRight_(nullptr),
     parent_(parent),
-    v_(elementLeft->completeSpaceSize(),0.0)
+    v_(elementLeft->completeSpaceSize(), 0)
 {
+    if (elementLeft->variableTimeStep())
+    {
+        nAdditionalEqns_ = 3;
+    }
+    else
+    {
+        nAdditionalEqns_ = 2;
+    }
+
     calcV(elementLeft, elementRight, v_);
     a_ = calcA(elementLeft, elementRight);
 }
 
-template<class CompType, class ThermoType>
-Foam::binaryNode<CompType, ThermoType>::binaryNode
-(
-    binaryNode<CompType, ThermoType> *bn
-)
-:
-    leafLeft_(bn->leafLeft()),
-    leafRight_(bn->leafRight()),
-    nodeLeft_(bn->nodeLeft()),
-    nodeRight_(bn->nodeRight()),
-    parent_(bn->parent()),
-    v_(bn->v()),
-    a_(bn->a())
-{}
-
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-
 template<class CompType, class ThermoType>
-void
-Foam::binaryNode<CompType, ThermoType>::calcV
+void Foam::binaryNode<CompType, ThermoType>::calcV
 (
     chemPointISAT<CompType, ThermoType>*& elementLeft,
     chemPointISAT<CompType, ThermoType>*& elementRight,
@@ -90,7 +81,8 @@ Foam::binaryNode<CompType, ThermoType>::calcV
     // LT is the transpose of the L matrix
     scalarSquareMatrix& LT = elementLeft->LT();
     bool mechReductionActive = elementLeft->chemistry().mechRed()->active();
-    // difference of composition in the full species domain
+
+    // Difference of composition in the full species domain
     scalarField phiDif(elementRight->phi() - elementLeft->phi());
     const scalarField& scaleFactor(elementLeft->scaleFactor());
     scalar epsTol = elementLeft->tolerance();
@@ -102,7 +94,7 @@ Foam::binaryNode<CompType, ThermoType>::calcV
         bool outOfIndexI = true;
         if (mechReductionActive)
         {
-            if (i<elementLeft->completeSpaceSize()-2)
+            if (i<elementLeft->completeSpaceSize() - nAdditionalEqns_)
             {
                 si = elementLeft->completeToSimplifiedIndex()[i];
                 outOfIndexI = (si == -1);
@@ -110,20 +102,21 @@ Foam::binaryNode<CompType, ThermoType>::calcV
             else // temperature and pressure
             {
                 outOfIndexI = false;
-                label dif = i-(elementLeft->completeSpaceSize()-2);
+                const label dif =
+                    i - (elementLeft->completeSpaceSize() - nAdditionalEqns_);
                 si = elementLeft->nActiveSpecies() + dif;
             }
         }
         if (!mechReductionActive || (mechReductionActive && !(outOfIndexI)))
         {
-            v[i]=0.0;
+            v[i] = 0;
             for (label j=0; j<elementLeft->completeSpaceSize(); j++)
             {
                 label sj = j;
                 bool outOfIndexJ = true;
                 if (mechReductionActive)
                 {
-                    if (j<elementLeft->completeSpaceSize()-2)
+                    if (j < elementLeft->completeSpaceSize() - nAdditionalEqns_)
                     {
                         sj = elementLeft->completeToSimplifiedIndex()[j];
                         outOfIndexJ = (sj==-1);
@@ -131,7 +124,12 @@ Foam::binaryNode<CompType, ThermoType>::calcV
                     else
                     {
                         outOfIndexJ = false;
-                        label dif = j-(elementLeft->completeSpaceSize()-2);
+                        const label dif =
+                            j
+                          - (
+                                elementLeft->completeSpaceSize()
+                              - nAdditionalEqns_
+                            );
                         sj = elementLeft->nActiveSpecies() + dif;
                     }
                 }
@@ -141,7 +139,7 @@ Foam::binaryNode<CompType, ThermoType>::calcV
                   ||(mechReductionActive && !(outOfIndexJ))
                 )
                 {
-                    // since L is a lower triangular matrix k=0->min(i, j)
+                    // Since L is a lower triangular matrix k=0->min(i, j)
                     for (label k=0; k<=min(si, sj); k++)
                     {
                         v[i] += LT(k, si)*LT(k, sj)*phiDif[j];
@@ -151,7 +149,7 @@ Foam::binaryNode<CompType, ThermoType>::calcV
         }
         else
         {
-            // when it is an inactive species the diagonal element of LT is
+            // When it is an inactive species the diagonal element of LT is
             // 1/(scaleFactor*epsTol)
             v[i] = phiDif[i]/sqr(scaleFactor[i]*epsTol);
         }
@@ -166,13 +164,13 @@ Foam::scalar Foam::binaryNode<CompType, ThermoType>::calcA
     chemPointISAT<CompType, ThermoType>* elementRight
 )
 {
-    scalar a = 0.0;
     scalarField phih((elementLeft->phi() + elementRight->phi())/2);
-    label completeSpaceSize = elementLeft->completeSpaceSize();
-    for (label i=0; i<completeSpaceSize; i++)
+    scalar a = 0;
+    forAll(phih, i)
     {
         a += v_[i]*phih[i];
     }
+
     return a;
 }