TDACChemistryModel: simplified, rationalized and automated the handling of variableTimeStep

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

@@ -25,6 +25,7 @@ License
 
 #include "TDACChemistryModel.H"
 #include "UniformField.H"
+#include "localEulerDdtScheme.H"
 #include "clockTime.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
@@ -37,6 +38,11 @@ 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),

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

@@ -85,6 +85,8 @@ class TDACChemistryModel
 {
     // Private member data
 
+        bool variableTimeStep_;
+
         label timeSteps_;
 
         // Mechanism reduction
@@ -159,11 +161,11 @@ public:
 
     // Member Functions
 
-        inline label timeSteps() const
+        //- Return true if the time-step is variable and/or non-uniform
-        {
+        inline bool variableTimeStep() const;
-            return timeSteps_;
-        }
 
+        //- 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;

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

@@ -103,25 +103,8 @@ Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::ISAT
         {
             scaleFactor_[Ysize + 2] = readScalar(scaleDict.lookup("deltaT"));
         }
-        else
-        {
-            // When the variableTimeStep option is false, if the application
-            // has variable time step activated, the maximum lifetime of a
-            // chemPoint should be 1 time step.
-            bool adjustTimeStep =
-                runTime_.controlDict().lookupOrDefault("adjustTimeStep", false);
-            if (chPMaxLifeTime_ > 1 && adjustTimeStep)
-            {
-                WarningInFunction
-                    << " variableTimeStep is not activate for ISAT while"
-                    << " the time step might be adjusted by the application."
-                    << nl
-                    << " This might lead to errors in the chemistry." << nl
-                    << " To avoid this warning either set chPMaxLifeTime to 1"
-                    << " or activate variableTimeStep." << endl;
-            }
-        }
     }
+
     if (this->variableTimeStep())
     {
         nAdditionalEqns_ = 3;

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

@@ -28,19 +28,15 @@ 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),
-    variableTimeStep_(false),
     nAdditionalEqns_(0)
-{
+{}
-}
 
 
 template<class CompType, class ThermoType>
@@ -56,10 +52,9 @@ Foam::binaryNode<CompType, ThermoType>::binaryNode
     nodeLeft_(nullptr),
     nodeRight_(nullptr),
     parent_(parent),
-    variableTimeStep_(elementLeft->variableTimeStep()),
     v_(elementLeft->completeSpaceSize(), 0)
 {
-    if (this->variableTimeStep_)
+    if (elementLeft->variableTimeStep())
     {
         nAdditionalEqns_ = 3;
     }
@@ -72,41 +67,11 @@ Foam::binaryNode<CompType, ThermoType>::binaryNode
     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()),
-    variableTimeStep_
-    (
-        this->coeffsDict_.lookupOrDefault("variableTimeStep", false)
-    ),
-    v_(bn->v()),
-    a_(bn->a())
-{
-    if (this->variableTimeStep_)
-    {
-        nAdditionalEqns_ = 3;
-    }
-    else
-    {
-        nAdditionalEqns_ = 2;
-    }
-}
-
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-
 template<class CompType, class ThermoType>
-void
+void Foam::binaryNode<CompType, ThermoType>::calcV
-Foam::binaryNode<CompType, ThermoType>::calcV
 (
     chemPointISAT<CompType, ThermoType>*& elementLeft,
     chemPointISAT<CompType, ThermoType>*& elementRight,
@@ -116,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();
@@ -143,7 +109,7 @@ Foam::binaryNode<CompType, ThermoType>::calcV
         }
         if (!mechReductionActive || (mechReductionActive && !(outOfIndexI)))
         {
-            v[i]=0.0;
+            v[i] = 0;
             for (label j=0; j<elementLeft->completeSpaceSize(); j++)
             {
                 label sj = j;
@@ -173,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];
@@ -183,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);
         }
@@ -198,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();
+    scalar a = 0;
-    for (label i=0; i<completeSpaceSize; i++)
+    forAll(phih, i)
     {
         a += v_[i]*phih[i];
     }
+
     return a;
 }
 

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

@@ -67,9 +67,6 @@ public:
     //- Parent node
     binaryNode<CompType, ThermoType>* parent_;
 
-    //- Switch to allow variable time step (off by default)
-    bool variableTimeStep_;
-
     //- Number of equations in addition to the species eqs.
     label nAdditionalEqns_;
 
@@ -137,11 +134,6 @@ public:
             chemPointISAT<CompType, ThermoType>* elementRight,
             binaryNode<CompType, ThermoType>* parent
         );
-        //- Construct from another binary node
-        binaryNode
-        (
-            binaryNode<CompType, ThermoType> *bn
-        );
 
 
     // Member functions

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

@@ -227,18 +227,14 @@ Foam::chemPointISAT<CompType, ThermoType>::chemPointISAT
     printProportion_(coeffsDict.lookupOrDefault("printProportion",false)),
     numRetrieve_(0),
     nLifeTime_(0),
-    variableTimeStep_
-    (
-        coeffsDict.lookupOrDefault("variableTimeStep", false)
-    ),
     completeToSimplifiedIndex_
     (
-        completeSpaceSize - (2 + (variableTimeStep_ == 1 ? 1 : 0))
+        completeSpaceSize - (2 + (variableTimeStep() == 1 ? 1 : 0))
     )
 {
     tolerance_ = tolerance;
 
-    if (this->variableTimeStep_)
+    if (variableTimeStep())
     {
         nAdditionalEqns_ = 3;
         iddeltaT_ = completeSpaceSize - 1;
@@ -342,12 +338,11 @@ Foam::chemPointISAT<CompType, ThermoType>::chemPointISAT
     maxNumNewDim_(p.maxNumNewDim()),
     numRetrieve_(0),
     nLifeTime_(0),
-    variableTimeStep_(p.variableTimeStep()),
     completeToSimplifiedIndex_(p.completeToSimplifiedIndex())
 {
     tolerance_ = p.tolerance();
 
-    if (this->variableTimeStep_)
+    if (variableTimeStep())
     {
         nAdditionalEqns_ = 3;
         idT_ = completeSpaceSize() - 3;
@@ -407,7 +402,7 @@ bool Foam::chemPointISAT<CompType, ThermoType>::inEOA(const scalarField& phiq)
 
             temp += LT_(si, dim)*dphi[idT_];
             temp += LT_(si, dim+1)*dphi[idp_];
-            if (variableTimeStep_)
+            if (variableTimeStep())
             {
                 temp += LT_(si, dim+2)*dphi[iddeltaT_];
             }
@@ -426,7 +421,7 @@ bool Foam::chemPointISAT<CompType, ThermoType>::inEOA(const scalarField& phiq)
     }
 
     // Temperature
-    if (variableTimeStep_)
+    if (variableTimeStep())
     {
         epsTemp +=
             sqr
@@ -447,7 +442,7 @@ bool Foam::chemPointISAT<CompType, ThermoType>::inEOA(const scalarField& phiq)
     }
 
     // Pressure
-    if (variableTimeStep_)
+    if (variableTimeStep())
     {
         epsTemp +=
             sqr
@@ -461,7 +456,7 @@ bool Foam::chemPointISAT<CompType, ThermoType>::inEOA(const scalarField& phiq)
         epsTemp += sqr(LT_(dim+1, dim+1)*dphi[idp_]);
     }
 
-    if (variableTimeStep_)
+    if (variableTimeStep())
     {
         epsTemp += sqr(LT_[dim+2][dim+2]*dphi[iddeltaT_]);
     }
@@ -477,7 +472,7 @@ bool Foam::chemPointISAT<CompType, ThermoType>::inEOA(const scalarField& phiq)
         propEps[idp_] =
             sqr(LT_(dim+1, dim+1)*dphi[idp_]);
 
-        if (variableTimeStep_)
+        if (variableTimeStep())
         {
             propEps[iddeltaT_] =
                 sqr(LT_[dim+2][dim+2]*dphi[iddeltaT_]);
@@ -572,7 +567,7 @@ bool Foam::chemPointISAT<CompType, ThermoType>::checkSolution
                 }
                 dRl += Avar(si, nActiveSpecies_)*dphi[idT_];
                 dRl += Avar(si, nActiveSpecies_+1)*dphi[idp_];
-                if (variableTimeStep_)
+                if (variableTimeStep())
                 {
                     dRl += Avar(si, nActiveSpecies_+2)*dphi[iddeltaT_];
                 }
@@ -719,7 +714,8 @@ bool Foam::chemPointISAT<CompType, ThermoType>::grow(const scalarField& phiq)
                 LTvar(initNActiveSpecies+1, initNActiveSpecies+1);
             A_(nActiveSpecies_+1, nActiveSpecies_+1)=
                 Avar(initNActiveSpecies+1, initNActiveSpecies+1);
-            if (variableTimeStep_)
+
+            if (variableTimeStep())
             {
                 LT_(nActiveSpecies_+2, nActiveSpecies_+2)=
                     LTvar(initNActiveSpecies+2, initNActiveSpecies+2);
@@ -755,9 +751,11 @@ bool Foam::chemPointISAT<CompType, ThermoType>::grow(const scalarField& phiq)
             }
             phiTilde[i] += LT_(i, j)*dphi[sj];
         }
+
         phiTilde[i] += LT_(i, dim-nAdditionalEqns_)*dphi[idT_];
         phiTilde[i] += LT_(i, dim-nAdditionalEqns_+1)*dphi[idp_];
-        if (variableTimeStep_)
+
+        if (variableTimeStep())
         {
             phiTilde[i] += LT_(i, dim-nAdditionalEqns_ + 2)*dphi[iddeltaT_];
         }
@@ -827,7 +825,7 @@ simplifiedToCompleteIndex
     {
         return completeSpaceSize_-nAdditionalEqns_ + 1;
     }
-    else if (variableTimeStep_ && (i == nActiveSpecies_ + 2))
+    else if (variableTimeStep() && (i == nActiveSpecies_ + 2))
     {
         return completeSpaceSize_-nAdditionalEqns_ + 2;
     }

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

@@ -198,9 +198,6 @@ class chemPointISAT
         //   to still live according to the maxChPLifeTime_ parameter
         label nLifeTime_;
 
-        //- Switch to allow variable time step (off by default)
-        Switch variableTimeStep_;
-
         List<label> completeToSimplifiedIndex_;
 
         //- Number of equations in addition to the species eqs.
@@ -397,9 +394,9 @@ public:
             return nLifeTime_;
         }
 
-        inline Switch variableTimeStep()
+        inline bool variableTimeStep() const
         {
-            return variableTimeStep_;
+            return chemistry_.variableTimeStep();
         }
 
         // ISAT functions

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

@@ -40,8 +40,7 @@ Foam::chemistryTabulationMethod<CompType, ThermoType>::chemistryTabulationMethod
     active_(coeffsDict_.lookupOrDefault<Switch>("active", false)),
     log_(coeffsDict_.lookupOrDefault<Switch>("log", false)),
     chemistry_(chemistry),
-    tolerance_(coeffsDict_.lookupOrDefault<scalar>("tolerance", 1e-4)),
+    tolerance_(coeffsDict_.lookupOrDefault<scalar>("tolerance", 1e-4))
-    variableTimeStep_(coeffsDict_.lookupOrDefault("variableTimeStep", false))
 {}
 
 

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

@@ -74,8 +74,6 @@ protected:
 
     scalar tolerance_;
 
-    //- Switch to allow variable time step (off by default)
-    const Switch variableTimeStep_;
 
 public:
 
@@ -134,7 +132,7 @@ public:
 
         inline bool variableTimeStep()
         {
-            return variableTimeStep_;
+            return chemistry_.variableTimeStep();
         }
 
         inline scalar tolerance() const

src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/makeChemistryTabulationMethods.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

tutorials/combustion/reactingFoam/laminar/counterFlowFlame2DLTS_GRI_TDAC/constant/chemistryProperties

@@ -78,8 +78,6 @@ tabulation
 
     printNumRetrieve   off;
 
-    variableTimeStep   on;
-
     // Tolerance used for retrieve and grow
     tolerance   1e-3;