TDACChemistryModel: Rationalised interface to reduction and tabulation methods

Most tabulation and reduction specific code in TDACChemistryModel has been moved
into the appropriate tabulation and reduction classes and the "none" methods of
each simplified so that they can be instantiated without any user input.  This
allows the "none" tabulation and reduction methods to be selected automatically
when the "tabulation" or "reduction" entries are not present in the
chemistryProperties dictionary so that the TDACChemistryModel can be run in the
same manner as the standardChemistryModel when tabulation and reduction are not
required.

This is the first step towards merging TDACChemistryModel with standardChemistryModel.

etc/caseDicts/solvers/chemistry/TDAC/chemistryProperties.cfg

@@ -8,9 +8,6 @@
 
 reduction
 {
-    // Activate reduction
-    active  on;
-
     // Switch logging of the reduction statistics and performance
     log         on;
 
@@ -41,9 +38,6 @@ reduction
 
 tabulation
 {
-    // Activate tabulation
-    active      on;
-
     // Switch logging of the tabulation statistics and performance
     log         on;
 

etc/caseDicts/solvers/chemistry/TDAC/chemistryPropertiesFlame.cfg

@@ -8,9 +8,6 @@
 
 reduction
 {
-    // Activate reduction
-    active  on;
-
     // Switch logging of the reduction statistics and performance
     log         on;
 
@@ -41,9 +38,6 @@ reduction
 
 tabulation
 {
-    // Activate tabulation
-    active      on;
-
     // Switch logging of the tabulation statistics and performance
     log         on;
 

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

@@ -37,25 +37,12 @@ Foam::TDACChemistryModel<ThermoType>::TDACChemistryModel
 )
 :
     standardChemistryModel<ThermoType>(thermo),
-    timeSteps_(0),
+    log_(this->lookupOrDefault("log", false)),
     NsDAC_(this->nSpecie_),
     completeC_(this->nSpecie_, 0),
     reactionsDisabled_(this->reactions_.size(), false),
     completeToSimplifiedIndex_(this->nSpecie_, -1),
-    simplifiedToCompleteIndex_(this->nSpecie_),
+    simplifiedToCompleteIndex_(this->nSpecie_)
-    tabulationResults_
-    (
-        IOobject
-        (
-            thermo.phasePropertyName("TabulationResults"),
-            this->time().timeName(),
-            this->mesh(),
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        this->mesh(),
-        scalar(0)
-    )
 {
     const basicSpecieMixture& composition = this->thermo().composition();
 
@@ -94,20 +81,7 @@ Foam::TDACChemistryModel<ThermoType>::TDACChemistryModel
         *this
     );
 
-    if (mechRed_->log())
+    if (log_)
-    {
-        cpuReduceFile_ = logFile("cpu_reduce.out");
-        nActiveSpeciesFile_ = logFile("nActiveSpecies.out");
-    }
-
-    if (tabulation_->log())
-    {
-        cpuAddFile_ = logFile("cpu_add.out");
-        cpuGrowFile_ = logFile("cpu_grow.out");
-        cpuRetrieveFile_ = logFile("cpu_retrieve.out");
-    }
-
-    if (mechRed_->log() || tabulation_->log())
     {
         cpuSolveFile_ = logFile("cpu_solve.out");
     }
@@ -362,7 +336,6 @@ void Foam::TDACChemistryModel<ThermoType>::jacobian
     d2TdtdT /= ccp;
 
     // d(dpdt)/dc = 0 (pressure is assumed constant)
-
     // d(dpdt)/dT = 0 (pressure is assumed constant)
 }
 
@@ -374,27 +347,15 @@ Foam::scalar Foam::TDACChemistryModel<ThermoType>::solve
     const DeltaTType& deltaT
 )
 {
-    // Increment counter of time-step
-    timeSteps_++;
-
     const bool reduced = mechRed_->active();
+    tabulation_->reset();
 
     const basicSpecieMixture& composition = this->thermo().composition();
 
     // CPU time analysis
-    const clockTime clockTime_= clockTime();
+    const clockTime clockTime_ = clockTime();
     clockTime_.timeIncrement();
-    scalar reduceMechCpuTime_ = 0;
-    scalar addNewLeafCpuTime_ = 0;
-    scalar growCpuTime_ = 0;
     scalar solveChemistryCpuTime_ = 0;
-    scalar searchISATCpuTime_ = 0;
-
-    this->resetTabulationResults();
-
-    // Average number of active species
-    scalar nActiveSpecies = 0;
-    scalar nAvg = 0;
 
     basicChemistryModel::correct();
 
@@ -445,15 +406,13 @@ Foam::scalar Foam::TDACChemistryModel<ThermoType>::solve
         // When tabulation is active (short-circuit evaluation for retrieve)
         // It first tries to retrieve the solution of the system with the
         // information stored through the tabulation method
-        if (tabulation_->active() && tabulation_->retrieve(phiq, Rphiq))
+        if (tabulation_->retrieve(phiq, Rphiq))
         {
             // Retrieved solution stored in Rphiq
             for (label i=0; i<this->nSpecie(); i++)
             {
                 c[i] = rhoi*Rphiq[i]/this->specieThermos_[i].W();
             }
-
-            searchISATCpuTime_ += clockTime_.timeIncrement();
         }
         // This position is reached when tabulation is not used OR
         // if the solution is not retrieved.
@@ -468,9 +427,6 @@ Foam::scalar Foam::TDACChemistryModel<ThermoType>::solve
             {
                 // Reduce mechanism change the number of species (only active)
                 mechRed_->reduceMechanism(pi, Ti, c, celli);
-                nActiveSpecies += mechRed_->NsSimp();
-                nAvg++;
-                reduceMechCpuTime_ += clockTime_.timeIncrement();
             }
 
             // Calculate the chemical source terms
@@ -506,13 +462,14 @@ Foam::scalar Foam::TDACChemistryModel<ThermoType>::solve
                 timeLeft -= dt;
             }
 
+            if (log_)
             {
                 solveChemistryCpuTime_ += clockTime_.timeIncrement();
             }
 
             // If tabulation is used, we add the information computed here to
             // the stored points (either expand or add)
-            if (tabulation_->active())
+            if (tabulation_->tabulates())
             {
                 forAll(c, i)
                 {
@@ -522,18 +479,7 @@ Foam::scalar Foam::TDACChemistryModel<ThermoType>::solve
                 Rphiq[Rphiq.size()-2] = pi;
                 Rphiq[Rphiq.size()-1] = deltaT[celli];
 
-                label growOrAdd =
+                tabulation_->add(phiq, Rphiq, celli, rhoi, deltaT[celli]);
-                    tabulation_->add(phiq, Rphiq, celli, rhoi, deltaT[celli]);
-                if (growOrAdd)
-                {
-                    this->setTabulationResultsAdd(celli);
-                    addNewLeafCpuTime_ += clockTime_.timeIncrement();
-                }
-                else
-                {
-                    this->setTabulationResultsGrow(celli);
-                    growCpuTime_ += clockTime_.timeIncrement();
-                }
             }
 
             // When operations are done and if mechanism reduction is active,
@@ -557,51 +503,15 @@ Foam::scalar Foam::TDACChemistryModel<ThermoType>::solve
         }
     }
 
-    if (mechRed_->log() || tabulation_->log())
+    if (log_)
     {
         cpuSolveFile_()
             << this->time().timeOutputValue()
             << "    " << solveChemistryCpuTime_ << endl;
     }
 
-    if (mechRed_->log())
+    mechRed_->update();
-    {
+    tabulation_->update();
-        cpuReduceFile_()
-            << this->time().timeOutputValue()
-            << "    " << reduceMechCpuTime_ << endl;
-    }
-
-    if (tabulation_->active())
-    {
-        // Every time-step, look if the tabulation should be updated
-        tabulation_->update();
-
-        // Write the performance of the tabulation
-        tabulation_->writePerformance();
-
-        if (tabulation_->log())
-        {
-            cpuRetrieveFile_()
-                << this->time().timeOutputValue()
-                << "    " << searchISATCpuTime_ << endl;
-
-            cpuGrowFile_()
-                << this->time().timeOutputValue()
-                << "    " << growCpuTime_ << endl;
-
-            cpuAddFile_()
-                << this->time().timeOutputValue()
-                << "    " << addNewLeafCpuTime_ << endl;
-        }
-    }
-
-    if (reduced && nAvg && mechRed_->log())
-    {
-        // Write average number of species
-        nActiveSpeciesFile_()
-            << this->time().timeOutputValue()
-            << "    " << nActiveSpecies/nAvg << endl;
-    }
 
     if (Pstream::parRun())
     {
@@ -647,34 +557,4 @@ Foam::scalar Foam::TDACChemistryModel<ThermoType>::solve
 }
 
 
-template<class ThermoType>
-void Foam::TDACChemistryModel<ThermoType>::
-setTabulationResultsAdd
-(
-    const label celli
-)
-{
-    tabulationResults_[celli] = 0;
-}
-
-
-template<class ThermoType>
-void Foam::TDACChemistryModel<ThermoType>::
-setTabulationResultsGrow(const label celli)
-{
-    tabulationResults_[celli] = 1;
-}
-
-
-template<class ThermoType>
-void Foam::TDACChemistryModel<ThermoType>::
-setTabulationResultsRetrieve
-(
-    const label celli
-)
-{
-    tabulationResults_[celli] = 2;
-}
-
-
 // ************************************************************************* //

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

@@ -85,7 +85,8 @@ class TDACChemistryModel
 {
     // Private member data
 
-        label timeSteps_;
+        //- Switch to select performance logging
+        Switch log_;
 
         // Mechanism reduction
         label NsDAC_;
@@ -94,35 +95,16 @@ class TDACChemistryModel
         Field<bool> reactionsDisabled_;
         Field<label> completeToSimplifiedIndex_;
         DynamicList<label> simplifiedToCompleteIndex_;
+
+        //- Reduction
         autoPtr<chemistryReductionMethod<ThermoType>> mechRed_;
 
-        // Tabulation
+        //- Tabulation
         autoPtr<chemistryTabulationMethod<ThermoType>> tabulation_;
 
-        // Log file for the average time spent reducing the chemistry
-        autoPtr<OFstream> cpuReduceFile_;
-
-        // Write average number of species
-        autoPtr<OFstream> nActiveSpeciesFile_;
-
-        //- Log file for the average time spent adding tabulated data
-        autoPtr<OFstream> cpuAddFile_;
-
-        //- Log file for the average time spent growing tabulated data
-        autoPtr<OFstream> cpuGrowFile_;
-
-        //- Log file for the average time spent retrieving tabulated data
-        autoPtr<OFstream> cpuRetrieveFile_;
-
         //- 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
 
@@ -154,9 +136,6 @@ public:
 
     // Member Functions
 
-        //- 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;
 
@@ -241,16 +220,6 @@ public:
 
             inline autoPtr<chemistryReductionMethod<ThermoType>>& mechRed();
 
-            inline tmp<volScalarField> tabulationResults() const;
-
-            void setTabulationResultsAdd(const label celli);
-
-            void setTabulationResultsGrow(const label celli);
-
-            void setTabulationResultsRetrieve(const label celli);
-
-            inline void resetTabulationResults();
-
 
     // Member Operators
 

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

@@ -25,13 +25,6 @@ License
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-template<class ThermoType>
-inline Foam::label Foam::TDACChemistryModel<ThermoType>::timeSteps() const
-{
-    return timeSteps_;
-}
-
-
 template<class ThermoType>
 inline Foam::autoPtr<Foam::OFstream>
 Foam::TDACChemistryModel<ThermoType>::logFile(const word& name) const
@@ -63,14 +56,6 @@ Foam::TDACChemistryModel<ThermoType>::mechRed()
 }
 
 
-template<class ThermoType>
-inline Foam::tmp<Foam::volScalarField>
-Foam::TDACChemistryModel<ThermoType>::tabulationResults() const
-{
-    return tabulationResults_;
-}
-
-
 template<class ThermoType>
 inline void Foam::TDACChemistryModel<ThermoType>::setActive(const label i)
 {
@@ -145,14 +130,4 @@ inline Foam::scalarField& Foam::TDACChemistryModel<ThermoType>::simplifiedC()
 }
 
 
-template<class ThermoType>
-void Foam::TDACChemistryModel<ThermoType>::resetTabulationResults()
-{
-    forAll(tabulationResults_, tabi)
-    {
-        tabulationResults_[tabi] = 2;
-    }
-}
-
-
 // ************************************************************************* //

src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DAC/DAC.C

@@ -34,14 +34,14 @@ Foam::chemistryReductionMethods::DAC<ThermoType>::DAC
     TDACChemistryModel<ThermoType>& chemistry
 )
 :
-    chemistryReductionMethod<ThermoType>(dict, chemistry),
+    chemistryReduction<ThermoType>(dict, chemistry),
     searchInitSet_(),
     zprime_(0),
     nbCLarge_(3),
-    sC_(this->nSpecie_,0),
+    sC_(this->nSpecie(),0),
-    sH_(this->nSpecie_,0),
+    sH_(this->nSpecie(),0),
-    sO_(this->nSpecie_,0),
+    sO_(this->nSpecie(),0),
-    sN_(this->nSpecie_,0),
+    sN_(this->nSpecie(),0),
     CO2Id_(-1),
     COId_(-1),
     HO2Id_(-1),
@@ -117,7 +117,7 @@ Foam::chemistryReductionMethods::DAC<ThermoType>::DAC
             this->coeffsDict_.template lookup<scalar>("NOxThreshold");
     }
 
-    for (label i=0; i<this->nSpecie_; i++)
+    for (label i=0; i<this->nSpecie(); i++)
     {
         const List<specieElement>& curSpecieComposition =
             chemistry.mixture().specieComposition(i);
@@ -240,28 +240,30 @@ void Foam::chemistryReductionMethods::DAC<ThermoType>::reduceMechanism
     const label li
 )
 {
+    chemistryReduction<ThermoType>::initReduceMechanism();
+
     scalarField& completeC(this->chemistry_.completeC());
     scalarField c1(this->chemistry_.nEqns(), 0.0);
-    for(label i=0; i<this->nSpecie_; i++)
+    for(label i=0; i<this->nSpecie(); i++)
     {
         c1[i] = c[i];
         completeC[i] = c[i];
     }
 
-    c1[this->nSpecie_] = T;
+    c1[this->nSpecie()] = T;
-    c1[this->nSpecie_+1] = p;
+    c1[this->nSpecie()+1] = p;
 
     // Compute the rAB matrix
-    RectangularMatrix<scalar> rABNum(this->nSpecie_,this->nSpecie_,0.0);
+    RectangularMatrix<scalar> rABNum(this->nSpecie(),this->nSpecie(),0.0);
-    scalarField PA(this->nSpecie_,0.0);
+    scalarField PA(this->nSpecie(),0.0);
-    scalarField CA(this->nSpecie_,0.0);
+    scalarField CA(this->nSpecie(),0.0);
 
     // Number of initialised rAB for each lines
-    Field<label> NbrABInit(this->nSpecie_,0);
+    Field<label> NbrABInit(this->nSpecie(),0);
     // Position of the initialised rAB, -1 when not initialised
-    RectangularMatrix<label> rABPos(this->nSpecie_, this->nSpecie_, -1);
+    RectangularMatrix<label> rABPos(this->nSpecie(), this->nSpecie(), -1);
     // Index of the other species involved in the rABNum
-    RectangularMatrix<label> rABOtherSpec(this->nSpecie_, this->nSpecie_, -1);
+    RectangularMatrix<label> rABOtherSpec(this->nSpecie(), this->nSpecie(), -1);
 
     forAll(this->chemistry_.reactions(), i)
     {
@@ -288,7 +290,7 @@ void Foam::chemistryReductionMethods::DAC<ThermoType>::reduceMechanism
         {
             label ss = R.lhs()[s].index;
             scalar sl = -R.lhs()[s].stoichCoeff; // vAi = v''-v' => here -v'
-            List<bool> deltaBi(this->nSpecie_, false);
+            List<bool> deltaBi(this->nSpecie(), false);
             FIFOStack<label> usedIndex;
             forAll(R.lhs(), j)
             {
@@ -351,7 +353,7 @@ void Foam::chemistryReductionMethods::DAC<ThermoType>::reduceMechanism
         {
             label ss = R.rhs()[s].index;
             scalar sl = R.rhs()[s].stoichCoeff; // vAi = v''-v' => here v''
-            List<bool> deltaBi(this->nSpecie_, false);
+            List<bool> deltaBi(this->nSpecie(), false);
             FIFOStack<label> usedIndex;
             forAll(R.lhs(), j)
             {
@@ -453,7 +455,7 @@ void Foam::chemistryReductionMethods::DAC<ThermoType>::reduceMechanism
         scalarList Na(nElements,0.0);
         scalarList Nal(nElements,0.0); // for large hydrocarbons
 
-        for (label i=0; i<this->nSpecie_; i++)
+        for (label i=0; i<this->nSpecie(); i++)
         {
             // Complete combustion products are not considered
             if
@@ -490,12 +492,12 @@ void Foam::chemistryReductionMethods::DAC<ThermoType>::reduceMechanism
 
     // Using the rAB matrix (numerator and denominator separated)
     // compute the R value according to the search initiating set
-    scalarField Rvalue(this->nSpecie_,0.0);
+    scalarField Rvalue(this->nSpecie(),0.0);
     label speciesNumber = 0;
 
     // Set all species to inactive and activate them according
     // to rAB and initial set
-    for (label i=0; i<this->nSpecie_; i++)
+    for (label i=0; i<this->nSpecie(); i++)
     {
         this->activeSpecies_[i] = false;
     }
@@ -667,15 +669,15 @@ void Foam::chemistryReductionMethods::DAC<ThermoType>::reduceMechanism
         }
     }
 
-    this->NsSimp_ = speciesNumber;
+    this->nActiveSpecies_ = speciesNumber;
     scalarField& simplifiedC(this->chemistry_.simplifiedC());
-    simplifiedC.setSize(this->NsSimp_+2);
+    simplifiedC.setSize(this->nActiveSpecies_+2);
     DynamicList<label>& s2c(this->chemistry_.simplifiedToCompleteIndex());
-    s2c.setSize(this->NsSimp_);
+    s2c.setSize(this->nActiveSpecies_);
     Field<label>& c2s(this->chemistry_.completeToSimplifiedIndex());
 
     label j = 0;
-    for (label i=0; i<this->nSpecie_; i++)
+    for (label i=0; i<this->nSpecie(); i++)
     {
         if (this->activeSpecies_[i])
         {
@@ -693,13 +695,15 @@ void Foam::chemistryReductionMethods::DAC<ThermoType>::reduceMechanism
         }
      }
 
-    simplifiedC[this->NsSimp_] = T;
+    simplifiedC[this->nActiveSpecies_] = T;
-    simplifiedC[this->NsSimp_+1] = p;
+    simplifiedC[this->nActiveSpecies_+1] = p;
-    this->chemistry_.setNsDAC(this->NsSimp_);
+    this->chemistry_.setNsDAC(this->nActiveSpecies_);
 
     // Change temporary Ns in chemistryModel
     // to make the function nEqns working
-    this->chemistry_.setNSpecie(this->NsSimp_);
+    this->chemistry_.setNSpecie(this->nActiveSpecies_);
+
+    chemistryReduction<ThermoType>::endReduceMechanism();
 }
 
 

src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DAC/DAC.H

@@ -89,7 +89,7 @@ SourceFiles
 #ifndef DAC_H
 #define DAC_H
 
-#include "chemistryReductionMethod.H"
+#include "chemistryReduction.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -104,7 +104,7 @@ namespace chemistryReductionMethods
 template<class ThermoType>
 class DAC
 :
-    public chemistryReductionMethod<ThermoType>
+    public chemistryReduction<ThermoType>
 {
     // Private Data
 

src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DRG/DRG.C

@@ -34,7 +34,7 @@ Foam::chemistryReductionMethods::DRG<ThermoType>::DRG
     TDACChemistryModel<ThermoType>& chemistry
 )
 :
-    chemistryReductionMethod<ThermoType>(dict, chemistry),
+    chemistryReduction<ThermoType>(dict, chemistry),
     searchInitSet_()
 {
     const wordHashSet initSet(this->coeffsDict_.lookup("initialSet"));
@@ -63,28 +63,30 @@ void Foam::chemistryReductionMethods::DRG<ThermoType>::reduceMechanism
     const label li
 )
 {
-    scalarField c1(this->nSpecie_+2, 0.0);
+    chemistryReduction<ThermoType>::initReduceMechanism();
 
-    for(label i=0; i<this->nSpecie_; i++)
+    scalarField c1(this->nSpecie()+2, 0.0);
+
+    for(label i=0; i<this->nSpecie(); i++)
     {
         c1[i] = c[i];
     }
 
-    c1[this->nSpecie_] = T;
+    c1[this->nSpecie()] = T;
-    c1[this->nSpecie_+1] = p;
+    c1[this->nSpecie()+1] = p;
 
     // Compute the rAB matrix
-    RectangularMatrix<scalar> rABNum(this->nSpecie_,this->nSpecie_,0.0);
+    RectangularMatrix<scalar> rABNum(this->nSpecie(),this->nSpecie(),0.0);
-    scalarField rABDen(this->nSpecie_,0.0);
+    scalarField rABDen(this->nSpecie(),0.0);
 
     // Number of initialised rAB for each lines
-    Field<label> NbrABInit(this->nSpecie_,0);
+    Field<label> NbrABInit(this->nSpecie(),0);
 
     // Position of the initialised rAB, -1 when not initialised
-    RectangularMatrix<label> rABPos(this->nSpecie_, this->nSpecie_, -1);
+    RectangularMatrix<label> rABPos(this->nSpecie(), this->nSpecie(), -1);
 
     // Index of the other species involved in the rABNum
-    RectangularMatrix<label> rABOtherSpec(this->nSpecie_, this->nSpecie_, -1);
+    RectangularMatrix<label> rABOtherSpec(this->nSpecie(), this->nSpecie(), -1);
 
     forAll(this->chemistry_.reactions(), i)
     {
@@ -152,7 +154,7 @@ void Foam::chemistryReductionMethods::DRG<ThermoType>::reduceMechanism
             // Absolute value of aggregated value
             scalar curwA = ((wA[id]>=0) ? wA[id] : -wA[id]);
 
-            List<bool> deltaBi(this->nSpecie_, false);
+            List<bool> deltaBi(this->nSpecie(), false);
             FIFOStack<label> usedIndex;
             forAll(R.lhs(), j)
             {
@@ -208,7 +210,7 @@ void Foam::chemistryReductionMethods::DRG<ThermoType>::reduceMechanism
 
     // Set all species to inactive and activate them according
     // to rAB and initial set
-    for (label i=0; i<this->nSpecie_; i++)
+    for (label i=0; i<this->nSpecie(); i++)
     {
         this->activeSpecies_[i] = false;
     }
@@ -293,12 +295,12 @@ void Foam::chemistryReductionMethods::DRG<ThermoType>::reduceMechanism
         }
     }
 
-    this->NsSimp_ = speciesNumber;
+    this->nActiveSpecies_ = speciesNumber;
-    this->chemistry_.simplifiedC().setSize(this->NsSimp_+2);
+    this->chemistry_.simplifiedC().setSize(this->nActiveSpecies_+2);
-    this->chemistry_.simplifiedToCompleteIndex().setSize(this->NsSimp_);
+    this->chemistry_.simplifiedToCompleteIndex().setSize(this->nActiveSpecies_);
 
     label j = 0;
-    for (label i=0; i<this->nSpecie_; i++)
+    for (label i=0; i<this->nSpecie(); i++)
     {
         if (this->activeSpecies_[i])
         {
@@ -316,13 +318,15 @@ void Foam::chemistryReductionMethods::DRG<ThermoType>::reduceMechanism
         }
     }
 
-    this->chemistry_.simplifiedC()[this->NsSimp_] = T;
+    this->chemistry_.simplifiedC()[this->nActiveSpecies_] = T;
-    this->chemistry_.simplifiedC()[this->NsSimp_+1] = p;
+    this->chemistry_.simplifiedC()[this->nActiveSpecies_+1] = p;
-    this->chemistry_.setNsDAC(this->NsSimp_);
+    this->chemistry_.setNsDAC(this->nActiveSpecies_);
 
     // Change temporary Ns in chemistryModel
     // to make the function nEqns working
-    this->chemistry_.setNSpecie(this->NsSimp_);
+    this->chemistry_.setNSpecie(this->nActiveSpecies_);
+
+    chemistryReduction<ThermoType>::endReduceMechanism();
 }
 
 

src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DRG/DRG.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2016-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2016-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -35,7 +35,7 @@ SourceFiles
 #ifndef DRG_H
 #define DRG_H
 
-#include "chemistryReductionMethod.H"
+#include "chemistryReduction.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -51,7 +51,7 @@ namespace chemistryReductionMethods
 template<class ThermoType>
 class DRG
 :
-    public chemistryReductionMethod<ThermoType>
+    public chemistryReduction<ThermoType>
 {
     // Private Data
 

src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DRGEP/DRGEP.C

@@ -35,14 +35,16 @@ Foam::chemistryReductionMethods::DRGEP<ThermoType>::DRGEP
     TDACChemistryModel<ThermoType>& chemistry
 )
 :
-    chemistryReductionMethod<ThermoType>(dict, chemistry),
+    chemistryReduction<ThermoType>(dict, chemistry),
     searchInitSet_(),
-    sC_(this->nSpecie_,0),
+    sC_(this->nSpecie(), 0),
-    sH_(this->nSpecie_,0),
+    sH_(this->nSpecie(), 0),
-    sO_(this->nSpecie_,0),
+    sO_(this->nSpecie(), 0),
-    sN_(this->nSpecie_,0),
+    sN_(this->nSpecie(), 0),
     NGroupBased_(50)
 {
+    chemistryReduction<ThermoType>::initReduceMechanism();
+
     const wordHashSet initSet(this->coeffsDict_.lookup("initialSet"));
     forAllConstIter(wordHashSet, initSet, iter)
     {
@@ -54,7 +56,7 @@ Foam::chemistryReductionMethods::DRGEP<ThermoType>::DRGEP
         NGroupBased_ = this->coeffsDict_.template lookup<label>("NGroupBased");
     }
 
-    for (label i=0; i<this->nSpecie_; i++)
+    for (label i=0; i<this->nSpecie(); i++)
     {
         const List<specieElement>& curSpecieComposition =
             chemistry.mixture().specieComposition(i);
@@ -110,26 +112,26 @@ void Foam::chemistryReductionMethods::DRGEP<ThermoType>::reduceMechanism
     scalarField& completeC(this->chemistry_.completeC());
     scalarField c1(this->chemistry_.nEqns(), 0.0);
 
-    for (label i=0; i<this->nSpecie_; i++)
+    for (label i=0; i<this->nSpecie(); i++)
     {
         c1[i] = c[i];
         completeC[i] = c[i];
     }
 
-    c1[this->nSpecie_] = T;
+    c1[this->nSpecie()] = T;
-    c1[this->nSpecie_+1] = p;
+    c1[this->nSpecie()+1] = p;
 
     // Compute the rAB matrix
-    RectangularMatrix<scalar> rABNum(this->nSpecie_,this->nSpecie_,0.0);
+    RectangularMatrix<scalar> rABNum(this->nSpecie(),this->nSpecie(),0.0);
-    scalarField PA(this->nSpecie_,0.0);
+    scalarField PA(this->nSpecie(),0.0);
-    scalarField CA(this->nSpecie_,0.0);
+    scalarField CA(this->nSpecie(),0.0);
 
     // Number of initialised rAB for each lines
-    Field<label> NbrABInit(this->nSpecie_,0);
+    Field<label> NbrABInit(this->nSpecie(),0);
     // Position of the initialised rAB, -1 when not initialised
-    RectangularMatrix<label> rABPos(this->nSpecie_, this->nSpecie_, -1);
+    RectangularMatrix<label> rABPos(this->nSpecie(), this->nSpecie(), -1);
     // Index of the other species involved in the rABNum
-    RectangularMatrix<label> rABOtherSpec(this->nSpecie_, this->nSpecie_, -1);
+    RectangularMatrix<label> rABOtherSpec(this->nSpecie(), this->nSpecie(), -1);
 
     scalarField omegaV(this->chemistry_.reactions().size());
     forAll(this->chemistry_.reactions(), i)
@@ -151,7 +153,7 @@ void Foam::chemistryReductionMethods::DRGEP<ThermoType>::reduceMechanism
         {
             label ss = R.lhs()[s].index;
             scalar sl = -R.lhs()[s].stoichCoeff; // vAi = v''-v' => here -v'
-            List<bool> deltaBi(this->nSpecie_, false);
+            List<bool> deltaBi(this->nSpecie(), false);
             FIFOStack<label> usedIndex;
             forAll(R.lhs(), j)
             {
@@ -211,7 +213,7 @@ void Foam::chemistryReductionMethods::DRGEP<ThermoType>::reduceMechanism
         {
             label ss = R.rhs()[s].index;
             scalar sl = R.rhs()[s].stoichCoeff; // vAi = v''-v' => here v''
-            List<bool> deltaBi(this->nSpecie_, false);
+            List<bool> deltaBi(this->nSpecie(), false);
             FIFOStack<label> usedIndex;
             forAll(R.lhs(), j)
             {
@@ -305,7 +307,7 @@ void Foam::chemistryReductionMethods::DRGEP<ThermoType>::reduceMechanism
     scalarList Ca(nElements,0.0);
 
     // for (label q=0; q<SIS.size(); q++)
-    for (label i=0; i<this->nSpecie_; i++)
+    for (label i=0; i<this->nSpecie(); i++)
     {
         Pa[0] += sC_[i]*max(0.0,(PA[i]-CA[i]));
         Ca[0] += sC_[i]*max(0.0,-(PA[i]-CA[i]));
@@ -319,13 +321,13 @@ void Foam::chemistryReductionMethods::DRGEP<ThermoType>::reduceMechanism
 
     // Using the rAB matrix (numerator and denominator separated)
     // compute the R value according to the search initiating set
-    scalarField Rvalue(this->nSpecie_,0.0);
+    scalarField Rvalue(this->nSpecie(),0.0);
     label speciesNumber = 0;
-    List<bool> disabledSpecies(this->nSpecie_,false);
+    List<bool> disabledSpecies(this->nSpecie(),false);
 
     // set all species to inactive and activate them according
     // to rAB and initial set
-    for (label i=0; i<this->nSpecie_; i++)
+    for (label i=0; i<this->nSpecie(); i++)
     {
         this->activeSpecies_[i] = false;
     }
@@ -453,7 +455,7 @@ void Foam::chemistryReductionMethods::DRGEP<ThermoType>::reduceMechanism
 
     // Group-based reduction
     // number of species disabled in the first step
-    label NDisabledSpecies(this->nSpecie_-speciesNumber);
+    label NDisabledSpecies(this->nSpecie()-speciesNumber);
 
     // while the number of removed species is greater than NGroupBased, the rAB
     // are reevaluated according to the group based definition for each loop the
@@ -507,7 +509,7 @@ void Foam::chemistryReductionMethods::DRGEP<ThermoType>::reduceMechanism
             {
                 label ss = R.lhs()[s].index;
                 scalar sl = -R.lhs()[s].stoichCoeff; // vAi = v''-v' => here -v'
-                List<bool> deltaBi(this->nSpecie_, false);
+                List<bool> deltaBi(this->nSpecie(), false);
                 bool alreadyDisabled(false);
                 FIFOStack<label> usedIndex;
                 forAll(R.lhs(), j)
@@ -561,7 +563,7 @@ void Foam::chemistryReductionMethods::DRGEP<ThermoType>::reduceMechanism
             {
                 label ss = R.rhs()[s].index;
                 scalar sl = R.rhs()[s].stoichCoeff; // vAi = v''-v' => here v''
-                List<bool> deltaBi(this->nSpecie_, false);
+                List<bool> deltaBi(this->nSpecie(), false);
                 bool alreadyDisabled(false);
                 FIFOStack<label> usedIndex;
                 forAll(R.lhs(), j)
@@ -686,15 +688,15 @@ void Foam::chemistryReductionMethods::DRGEP<ThermoType>::reduceMechanism
         }
     }
 
-    this->NsSimp_ = speciesNumber;
+    this->nActiveSpecies_ = speciesNumber;
     scalarField& simplifiedC(this->chemistry_.simplifiedC());
-    simplifiedC.setSize(this->NsSimp_+2);
+    simplifiedC.setSize(this->nActiveSpecies_+2);
     DynamicList<label>& s2c(this->chemistry_.simplifiedToCompleteIndex());
-    s2c.setSize(this->NsSimp_);
+    s2c.setSize(this->nActiveSpecies_);
     Field<label>& c2s(this->chemistry_.completeToSimplifiedIndex());
 
     label j = 0;
-    for (label i=0; i<this->nSpecie_; i++)
+    for (label i=0; i<this->nSpecie(); i++)
     {
         if (this->activeSpecies_[i])
         {
@@ -711,12 +713,14 @@ void Foam::chemistryReductionMethods::DRGEP<ThermoType>::reduceMechanism
             c2s[i] = -1;
         }
     }
-    simplifiedC[this->NsSimp_] = T;
+    simplifiedC[this->nActiveSpecies_] = T;
-    simplifiedC[this->NsSimp_+1] = p;
+    simplifiedC[this->nActiveSpecies_+1] = p;
-    this->chemistry_.setNsDAC(this->NsSimp_);
+    this->chemistry_.setNsDAC(this->nActiveSpecies_);
     // change temporary Ns in chemistryModel
     // to make the function nEqns working
-    this->chemistry_.setNSpecie(this->NsSimp_);
+    this->chemistry_.setNSpecie(this->nActiveSpecies_);
+
+    chemistryReduction<ThermoType>::endReduceMechanism();
 }
 
 

src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DRGEP/DRGEP.H

@@ -103,7 +103,7 @@ SourceFiles
 #ifndef DRGEP_H
 #define DRGEP_H
 
-#include "chemistryReductionMethod.H"
+#include "chemistryReduction.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -113,13 +113,13 @@ namespace chemistryReductionMethods
 {
 
 /*---------------------------------------------------------------------------*\
-                             Class ode Declaration
+                             Class DRGEP Declaration
 \*---------------------------------------------------------------------------*/
 
 template<class ThermoType>
 class DRGEP
 :
-    public chemistryReductionMethod<ThermoType>
+    public chemistryReduction<ThermoType>
 {
     // Private Data
 

src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/EFA/EFA.C

@@ -35,14 +35,14 @@ Foam::chemistryReductionMethods::EFA<ThermoType>::EFA
     TDACChemistryModel<ThermoType>& chemistry
 )
 :
-    chemistryReductionMethod<ThermoType>(dict, chemistry),
+    chemistryReduction<ThermoType>(dict, chemistry),
-    sC_(this->nSpecie_,0),
+    sC_(this->nSpecie(),0),
-    sH_(this->nSpecie_,0),
+    sH_(this->nSpecie(),0),
-    sO_(this->nSpecie_,0),
+    sO_(this->nSpecie(),0),
-    sN_(this->nSpecie_,0),
+    sN_(this->nSpecie(),0),
     sortPart_(0.05)
 {
-    for (label i=0; i<this->nSpecie_; i++)
+    for (label i=0; i<this->nSpecie(); i++)
     {
         const List<specieElement>& curSpecieComposition =
             chemistry.mixture().specieComposition(i);
@@ -99,33 +99,35 @@ void Foam::chemistryReductionMethods::EFA<ThermoType>::reduceMechanism
     const label li
 )
 {
+    chemistryReduction<ThermoType>::initReduceMechanism();
+
     scalarField& completeC(this->chemistry_.completeC());
     scalarField c1(this->chemistry_.nEqns(), 0.0);
 
-    for (label i=0; i<this->nSpecie_; i++)
+    for (label i=0; i<this->nSpecie(); i++)
     {
         c1[i] = c[i];
         completeC[i] = c[i];
     }
 
-    c1[this->nSpecie_] = T;
+    c1[this->nSpecie()] = T;
-    c1[this->nSpecie_+1] = p;
+    c1[this->nSpecie()+1] = p;
 
 
     // Number of initialised rAB for each lines
-    Field<label> NbrABInit(this->nSpecie_,0);
+    Field<label> NbrABInit(this->nSpecie(),0);
 
     // Position of the initialised rAB, -1 when not initialised
-    RectangularMatrix<label> rABPos(this->nSpecie_, this->nSpecie_, -1);
+    RectangularMatrix<label> rABPos(this->nSpecie(), this->nSpecie(), -1);
-    RectangularMatrix<scalar> CFluxAB(this->nSpecie_, this->nSpecie_, 0.0);
+    RectangularMatrix<scalar> CFluxAB(this->nSpecie(), this->nSpecie(), 0.0);
-    RectangularMatrix<scalar> HFluxAB(this->nSpecie_, this->nSpecie_, 0.0);
+    RectangularMatrix<scalar> HFluxAB(this->nSpecie(), this->nSpecie(), 0.0);
-    RectangularMatrix<scalar> OFluxAB(this->nSpecie_, this->nSpecie_, 0.0);
+    RectangularMatrix<scalar> OFluxAB(this->nSpecie(), this->nSpecie(), 0.0);
-    RectangularMatrix<scalar> NFluxAB(this->nSpecie_, this->nSpecie_, 0.0);
+    RectangularMatrix<scalar> NFluxAB(this->nSpecie(), this->nSpecie(), 0.0);
     scalar CFlux(0.0), HFlux(0.0), OFlux(0.0), NFlux(0.0);
     label nbPairs(0);
 
     // Index of the other species involved in the rABNum
-    RectangularMatrix<label> rABOtherSpec(this->nSpecie_, this->nSpecie_, -1);
+    RectangularMatrix<label> rABOtherSpec(this->nSpecie(), this->nSpecie(), -1);
 
     forAll(this->chemistry_.reactions(), i)
     {
@@ -243,7 +245,7 @@ void Foam::chemistryReductionMethods::EFA<ThermoType>::reduceMechanism
     // Select species according to the total flux cutoff (1-tolerance)
     // of the flux is included
     label speciesNumber = 0;
-    for (label i=0; i<this->nSpecie_; i++)
+    for (label i=0; i<this->nSpecie(); i++)
     {
         this->activeSpecies_[i] = false;
     }
@@ -254,7 +256,7 @@ void Foam::chemistryReductionMethods::EFA<ThermoType>::reduceMechanism
         labelList source(nbPairs);
         labelList sink(nbPairs);
         label nP(0);
-        for (int A=0; A<this->nSpecie_ ; A++)
+        for (int A=0; A<this->nSpecie() ; A++)
         {
             for (int j=0; j<NbrABInit[A]; j++)
             {
@@ -320,7 +322,7 @@ void Foam::chemistryReductionMethods::EFA<ThermoType>::reduceMechanism
         labelList source(nbPairs);
         labelList sink(nbPairs);
         label nP(0);
-        for (int A=0; A<this->nSpecie_ ; A++)
+        for (int A=0; A<this->nSpecie() ; A++)
         {
             for (int j=0; j<NbrABInit[A]; j++)
             {
@@ -383,7 +385,7 @@ void Foam::chemistryReductionMethods::EFA<ThermoType>::reduceMechanism
         labelList source(nbPairs);
         labelList sink(nbPairs);
         label nP(0);
-        for (int A=0; A<this->nSpecie_ ; A++)
+        for (int A=0; A<this->nSpecie() ; A++)
         {
             for (int j=0; j<NbrABInit[A]; j++)
             {
@@ -445,7 +447,7 @@ void Foam::chemistryReductionMethods::EFA<ThermoType>::reduceMechanism
         labelList source(nbPairs);
         labelList sink(nbPairs);
         label nP(0);
-        for (int A=0; A<this->nSpecie_ ; A++)
+        for (int A=0; A<this->nSpecie() ; A++)
         {
             for (int j=0; j<NbrABInit[A]; j++)
             {
@@ -529,15 +531,15 @@ void Foam::chemistryReductionMethods::EFA<ThermoType>::reduceMechanism
         }
     }
 
-    this->NsSimp_ = speciesNumber;
+    this->nActiveSpecies_ = speciesNumber;
     scalarField& simplifiedC(this->chemistry_.simplifiedC());
-    simplifiedC.setSize(this->NsSimp_+2);
+    simplifiedC.setSize(this->nActiveSpecies_+2);
     DynamicList<label>& s2c(this->chemistry_.simplifiedToCompleteIndex());
-    s2c.setSize(this->NsSimp_);
+    s2c.setSize(this->nActiveSpecies_);
     Field<label>& c2s(this->chemistry_.completeToSimplifiedIndex());
     label j = 0;
 
-    for (label i=0; i<this->nSpecie_; i++)
+    for (label i=0; i<this->nSpecie(); i++)
     {
         if (this->activeSpecies_[i])
         {
@@ -554,12 +556,14 @@ void Foam::chemistryReductionMethods::EFA<ThermoType>::reduceMechanism
             c2s[i] = -1;
         }
     }
-    simplifiedC[this->NsSimp_] = T;
+    simplifiedC[this->nActiveSpecies_] = T;
-    simplifiedC[this->NsSimp_+1] = p;
+    simplifiedC[this->nActiveSpecies_+1] = p;
-    this->chemistry_.setNsDAC(this->NsSimp_);
+    this->chemistry_.setNsDAC(this->nActiveSpecies_);
     // change temporary Ns in chemistryModel
     // to make the function nEqns working
-    this->chemistry_.setNSpecie(this->NsSimp_);
+    this->chemistry_.setNSpecie(this->nActiveSpecies_);
+
+    chemistryReduction<ThermoType>::endReduceMechanism();
 }
 
 

src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/EFA/EFA.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2016-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2016-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -34,7 +34,7 @@ SourceFiles
 #ifndef EFA_H
 #define EFA_H
 
-#include "chemistryReductionMethod.H"
+#include "chemistryReduction.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -44,13 +44,13 @@ namespace chemistryReductionMethods
 {
 
 /*---------------------------------------------------------------------------*\
-                             Class ode Declaration
+                             Class EFA Declaration
 \*---------------------------------------------------------------------------*/
 
 template<class ThermoType>
 class EFA
 :
-    public chemistryReductionMethod<ThermoType>
+    public chemistryReduction<ThermoType>
 {
     // Private Data
 

src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/PFA/PFA.C

@@ -34,7 +34,7 @@ Foam::chemistryReductionMethods::PFA<ThermoType>::PFA
     TDACChemistryModel<ThermoType>& chemistry
 )
 :
-    chemistryReductionMethod<ThermoType>(dict, chemistry),
+    chemistryReduction<ThermoType>(dict, chemistry),
     searchInitSet_()
 {
     const wordHashSet initSet(this->coeffsDict_.lookup("initialSet"));
@@ -63,30 +63,32 @@ void Foam::chemistryReductionMethods::PFA<ThermoType>::reduceMechanism
     const label li
 )
 {
+    chemistryReduction<ThermoType>::initReduceMechanism();
+
     scalarField& completeC(this->chemistry_.completeC());
     scalarField c1(this->chemistry_.nEqns(), 0.0);
 
-    for (label i=0; i<this->nSpecie_; i++)
+    for (label i=0; i<this->nSpecie(); i++)
     {
         c1[i] = c[i];
         completeC[i] = c[i];
     }
 
-    c1[this->nSpecie_] = T;
+    c1[this->nSpecie()] = T;
-    c1[this->nSpecie_+1] = p;
+    c1[this->nSpecie()+1] = p;
 
     // Compute the rAB matrix
-    RectangularMatrix<scalar> PAB(this->nSpecie_,this->nSpecie_,0.0);
+    RectangularMatrix<scalar> PAB(this->nSpecie(),this->nSpecie(),0.0);
-    RectangularMatrix<scalar> CAB(this->nSpecie_,this->nSpecie_,0.0);
+    RectangularMatrix<scalar> CAB(this->nSpecie(),this->nSpecie(),0.0);
-    scalarField PA(this->nSpecie_,0.0);
+    scalarField PA(this->nSpecie(),0.0);
-    scalarField CA(this->nSpecie_,0.0);
+    scalarField CA(this->nSpecie(),0.0);
 
     // Number of initialised rAB for each lines
-    Field<label> NbrABInit(this->nSpecie_,0);
+    Field<label> NbrABInit(this->nSpecie(),0);
     // Position of the initialised rAB, -1 when not initialised
-    RectangularMatrix<label> rABPos(this->nSpecie_, this->nSpecie_, -1);
+    RectangularMatrix<label> rABPos(this->nSpecie(), this->nSpecie(), -1);
     // Index of the other species involved in the rABNum
-    RectangularMatrix<label> rABOtherSpec(this->nSpecie_, this->nSpecie_, -1);
+    RectangularMatrix<label> rABOtherSpec(this->nSpecie(), this->nSpecie(), -1);
 
     forAll(this->chemistry_.reactions(), i)
     {
@@ -149,7 +151,7 @@ void Foam::chemistryReductionMethods::PFA<ThermoType>::reduceMechanism
         {
             label curID = wAID[id];
             scalar curwA = wA[id];
-            List<bool> deltaBi(this->nSpecie_, false);
+            List<bool> deltaBi(this->nSpecie(), false);
             FIFOStack<label> usedIndex;
             forAll(R.lhs(),j)
             {
@@ -236,19 +238,19 @@ void Foam::chemistryReductionMethods::PFA<ThermoType>::reduceMechanism
     // is a connection of second generation and it will be aggregated in the
     // final step to evaluate the total connection strength (or path flux).
     // Compute rsecond=rAri*rriB with A!=ri!=B
-    RectangularMatrix<scalar> PAB2nd(this->nSpecie_,this->nSpecie_,0.0);
+    RectangularMatrix<scalar> PAB2nd(this->nSpecie(),this->nSpecie(),0.0);
-    RectangularMatrix<scalar> CAB2nd(this->nSpecie_,this->nSpecie_,0.0);
+    RectangularMatrix<scalar> CAB2nd(this->nSpecie(),this->nSpecie(),0.0);
 
     // Number of initialised rAB for each lines
-    Field<label> NbrABInit2nd(this->nSpecie_, 0);
+    Field<label> NbrABInit2nd(this->nSpecie(), 0);
 
     // Position of the initialised rAB, -1 when not initialised
-    RectangularMatrix<label> rABPos2nd(this->nSpecie_, this->nSpecie_, -1);
+    RectangularMatrix<label> rABPos2nd(this->nSpecie(), this->nSpecie(), -1);
 
     // Index of the other species involved in the rABNum
     RectangularMatrix<label> rABOtherSpec2nd
     (
-        this->nSpecie_, this->nSpecie_, -1
+        this->nSpecie(), this->nSpecie(), -1
     );
 
     forAll(NbrABInit, A)
@@ -291,7 +293,7 @@ void Foam::chemistryReductionMethods::PFA<ThermoType>::reduceMechanism
 
     // set all species to inactive and activate them according
     // to rAB and initial set
-    for (label i=0; i<this->nSpecie_; i++)
+    for (label i=0; i<this->nSpecie(); i++)
     {
         this->activeSpecies_[i] = false;
     }
@@ -388,15 +390,15 @@ void Foam::chemistryReductionMethods::PFA<ThermoType>::reduceMechanism
         }
     }
 
-    this->NsSimp_ = speciesNumber;
+    this->nActiveSpecies_ = speciesNumber;
     scalarField& simplifiedC(this->chemistry_.simplifiedC());
-    simplifiedC.setSize(this->NsSimp_+2);
+    simplifiedC.setSize(this->nActiveSpecies_+2);
     DynamicList<label>& s2c(this->chemistry_.simplifiedToCompleteIndex());
-    s2c.setSize(this->NsSimp_);
+    s2c.setSize(this->nActiveSpecies_);
     Field<label>& c2s(this->chemistry_.completeToSimplifiedIndex());
 
     label j = 0;
-    for (label i=0; i<this->nSpecie_; i++)
+    for (label i=0; i<this->nSpecie(); i++)
     {
         if (this->activeSpecies_[i])
         {
@@ -413,12 +415,14 @@ void Foam::chemistryReductionMethods::PFA<ThermoType>::reduceMechanism
             c2s[i] = -1;
         }
     }
-    simplifiedC[this->NsSimp_] = T;
+    simplifiedC[this->nActiveSpecies_] = T;
-    simplifiedC[this->NsSimp_+1] = p;
+    simplifiedC[this->nActiveSpecies_+1] = p;
-    this->chemistry_.setNsDAC(this->NsSimp_);
+    this->chemistry_.setNsDAC(this->nActiveSpecies_);
     // change temporary Ns in chemistryModel
     // to make the function nEqns working
-    this->chemistry_.setNSpecie(this->NsSimp_);
+    this->chemistry_.setNSpecie(this->nActiveSpecies_);
+
+    chemistryReduction<ThermoType>::endReduceMechanism();
 }
 
 

src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/PFA/PFA.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2016-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2016-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -35,7 +35,7 @@ SourceFiles
 #ifndef PFA_H
 #define PFA_H
 
-#include "chemistryReductionMethod.H"
+#include "chemistryReduction.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -45,13 +45,13 @@ namespace chemistryReductionMethods
 {
 
 /*---------------------------------------------------------------------------*\
-                             Class ode Declaration
+                             Class PFA Declaration
 \*---------------------------------------------------------------------------*/
 
 template<class ThermoType>
 class PFA
 :
-    public chemistryReductionMethod<ThermoType>
+    public chemistryReduction<ThermoType>
 {
     // Private Data
 

src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/chemistryReduction/chemistryReduction.C

@@ -0,0 +1,112 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2021 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "chemistryReduction.H"
+#include "TDACChemistryModel.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class ThermoType>
+Foam::chemistryReduction<ThermoType>::chemistryReduction
+(
+    const Foam::IOdictionary& dict,
+    Foam::TDACChemistryModel<ThermoType>& chemistry
+)
+:
+    chemistryReductionMethod<ThermoType>(dict, chemistry),
+    coeffsDict_(dict.subDict("reduction")),
+    chemistry_(chemistry),
+    activeSpecies_(chemistry.nSpecie(), false),
+    log_(coeffsDict_.lookupOrDefault<Switch>("log", false)),
+    tolerance_(coeffsDict_.lookupOrDefault<scalar>("tolerance", 1e-4)),
+    clockTime_(clockTime()),
+    sumnActiveSpecies_(0),
+    sumn_(0),
+    reduceMechCpuTime_(0)
+{
+    if (log_)
+    {
+        cpuReduceFile_ = chemistry.logFile("cpu_reduce.out");
+        nActiveSpeciesFile_ = chemistry.logFile("nActiveSpecies.out");
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class ThermoType>
+Foam::chemistryReduction<ThermoType>::~chemistryReduction()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class ThermoType>
+void Foam::chemistryReduction<ThermoType>::initReduceMechanism()
+{
+    if (log_)
+    {
+        clockTime_.timeIncrement();
+    }
+}
+
+
+template<class ThermoType>
+void Foam::chemistryReduction<ThermoType>::endReduceMechanism()
+{
+    if (log_)
+    {
+        sumnActiveSpecies_ += this->nActiveSpecies_;
+        sumn_++;
+        reduceMechCpuTime_ += clockTime_.timeIncrement();
+    }
+}
+
+
+template<class ThermoType>
+void Foam::chemistryReduction<ThermoType>::update()
+{
+    if (log_)
+    {
+        cpuReduceFile_()
+            << this->chemistry_.time().timeOutputValue()
+            << "    " << reduceMechCpuTime_ << endl;
+
+        if (sumn_)
+        {
+            // Write average number of species
+            nActiveSpeciesFile_()
+                << this->chemistry_.time().timeOutputValue()
+                << "    " << sumnActiveSpecies_/sumn_ << endl;
+        }
+
+        sumnActiveSpecies_ = 0;
+        sumn_ = 0;
+        reduceMechCpuTime_ = 0;
+    }
+}
+
+
+// ************************************************************************* //

src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/chemistryReduction/chemistryReduction.H

@@ -0,0 +1,156 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2021 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/>.
+
+Class
+    Foam::chemistryReduction
+
+Description
+    An abstract class for methods of chemical mechanism reduction
+
+SourceFiles
+    chemistryReduction.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef chemistryReduction_H
+#define chemistryReduction_H
+
+#include "chemistryReductionMethod.H"
+#include "Switch.H"
+#include "int64.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                       Class chemistrySolver Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class ThermoType>
+class chemistryReduction
+:
+    public chemistryReductionMethod<ThermoType>
+{
+
+protected:
+
+    // Protected Data
+
+        //- Dictionary that store the algorithm data
+        const dictionary coeffsDict_;
+
+        TDACChemistryModel<ThermoType>& chemistry_;
+
+        //- List of active species (active = true)
+        List<bool> activeSpecies_;
+
+        void initReduceMechanism();
+
+        void endReduceMechanism();
+
+
+private:
+
+        //- Switch to select performance logging
+        Switch log_;
+
+        //- Tolerance for the mechanism reduction algorithm
+        scalar tolerance_;
+
+        const clockTime clockTime_;
+
+        int64_t sumnActiveSpecies_;
+        int64_t sumn_;
+        scalar reduceMechCpuTime_;
+
+        // Log file for the average time spent reducing the chemistry
+        autoPtr<OFstream> cpuReduceFile_;
+
+        // Write average number of species
+        autoPtr<OFstream> nActiveSpeciesFile_;
+
+
+public:
+
+    // Constructors
+
+        //- Construct from components
+        chemistryReduction
+        (
+            const IOdictionary& dict,
+            TDACChemistryModel<ThermoType>& chemistry
+        );
+
+
+    //- Destructor
+    virtual ~chemistryReduction();
+
+
+    // Member Functions
+
+        //- Return mechanism reduction active
+        virtual bool active() const
+        {
+            return true;
+        }
+
+        //- Return the tolerance
+        inline scalar tolerance() const;
+
+        //- Return the active species
+        inline const List<bool>& activeSpecies() const;
+
+        //- Reduce the mechanism
+        virtual void reduceMechanism
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c,
+            const label li
+        ) = 0;
+
+        virtual void update();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "chemistryReductionI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "chemistryReduction.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/chemistryReduction/chemistryReductionI.H

@@ -0,0 +1,53 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2021 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/>.
+
+Class
+    Foam::chemistryReduction
+
+Description
+    An abstract class for reducing chemical mechanism
+
+SourceFiles
+    chemistryReduction.C
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class ThermoType>
+inline Foam::scalar
+Foam::chemistryReduction<ThermoType>::tolerance() const
+{
+    return tolerance_;
+}
+
+
+template<class ThermoType>
+inline const Foam::List<bool>&
+Foam::chemistryReduction<ThermoType>::activeSpecies() const
+{
+    return activeSpecies_;
+}
+
+
+// ************************************************************************* //

src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/chemistryReductionMethod/chemistryReductionMethod.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2016-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2016-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -35,15 +35,8 @@ Foam::chemistryReductionMethod<ThermoType>::chemistryReductionMethod
     Foam::TDACChemistryModel<ThermoType>& chemistry
 )
 :
-    dict_(dict),
-    coeffsDict_(dict.subDict("reduction")),
-    active_(coeffsDict_.lookupOrDefault<Switch>("active", false)),
-    log_(coeffsDict_.lookupOrDefault<Switch>("log", false)),
-    chemistry_(chemistry),
-    activeSpecies_(chemistry.nSpecie(), false),
-    NsSimp_(chemistry.nSpecie()),
     nSpecie_(chemistry.nSpecie()),
-    tolerance_(coeffsDict_.lookupOrDefault<scalar>("tolerance", 1e-4))
+    nActiveSpecies_(chemistry.nSpecie())
 {}
 
 

src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/chemistryReductionMethod/chemistryReductionMethod.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2016-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2016-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -36,7 +36,6 @@ SourceFiles
 #define chemistryReductionMethod_H
 
 #include "IOdictionary.H"
-#include "Switch.H"
 #include "scalarField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -54,33 +53,18 @@ class TDACChemistryModel;
 template<class ThermoType>
 class chemistryReductionMethod
 {
+    // Private Data
+
+        //- Cache the number of species
+        const label nSpecie_;
+
 
 protected:
 
-    const IOdictionary& dict_;
+    // Protected Data
 
-    //- Dictionary that store the algorithm data
+        //- Number of active species
-    const dictionary coeffsDict_;
+        label nActiveSpecies_;
-
-    //- Is mechanism reduction active?
-    Switch active_;
-
-    //- Switch to select performance logging
-    Switch log_;
-
-    TDACChemistryModel<ThermoType>& chemistry_;
-
-    //- List of active species (active = true)
-    List<bool> activeSpecies_;
-
-    //- Number of active species
-    label NsSimp_;
-
-    //- Number of species
-    const label nSpecie_;
-
-    //- Tolerance for the mechanism reduction algorithm
-    scalar tolerance_;
 
 
 public:
@@ -129,19 +113,13 @@ public:
     // Member Functions
 
         //- Is mechanism reduction active?
-        inline bool active() const;
+        virtual bool active() const = 0;
 
-        //- Is performance data logging enabled?
+        //- Return the number of species
-        inline bool log() const;
+        inline label nSpecie();
-
-        //- Return the active species
-        inline const List<bool>& activeSpecies() const;
 
         //- Return the number of active species
-        inline label NsSimp();
+        inline label nActiveSpecies();
-
-        //- Return the initial number of species
-        inline label nSpecie();
 
         //- Return the tolerance
         inline scalar tolerance() const;
@@ -154,8 +132,11 @@ public:
             const scalarField& c,
             const label li
         ) = 0;
+
+        virtual void update() = 0;
 };
 
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/chemistryReductionMethod/chemistryReductionMethodI.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2016-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2016-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -35,30 +35,9 @@ SourceFiles
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class ThermoType>
-inline bool Foam::chemistryReductionMethod<ThermoType>::active() const
+inline Foam::label Foam::chemistryReductionMethod<ThermoType>::nActiveSpecies()
 {
-    return active_;
+    return nActiveSpecies_;
-}
-
-
-template<class ThermoType>
-inline bool Foam::chemistryReductionMethod<ThermoType>::log() const
-{
-    return active_ && log_;
-}
-
-
-template<class ThermoType>
-inline const Foam::List<bool>&
-Foam::chemistryReductionMethod<ThermoType>::activeSpecies() const
-{
-    return activeSpecies_;
-}
-
-template<class ThermoType>
-inline Foam::label Foam::chemistryReductionMethod<ThermoType>::NsSimp()
-{
-    return NsSimp_;
 }
 
 
@@ -69,12 +48,4 @@ inline Foam::label Foam::chemistryReductionMethod<ThermoType>::nSpecie()
 }
 
 
-template<class ThermoType>
-inline Foam::scalar
-Foam::chemistryReductionMethod<ThermoType>::tolerance() const
-{
-    return tolerance_;
-}
-
-
 // ************************************************************************* //

src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/chemistryReductionMethod/chemistryReductionMethodNew.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2016-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2016-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -24,6 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "chemistryReductionMethod.H"
+#include "noChemistryReduction.H"
 #include "Time.H"
 
 // * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * //
@@ -36,53 +37,67 @@ Foam::chemistryReductionMethod<ThermoType>::New
     TDACChemistryModel<ThermoType>& chemistry
 )
 {
-    const dictionary& reductionDict(dict.subDict("reduction"));
+    if (dict.found("reduction"))
-
-    const word methodName(reductionDict.lookup("method"));
-
-    Info<< "Selecting chemistry reduction method " << methodName << endl;
-
-    const word methodTypeName =
-        methodName + '<' + ThermoType::typeName() + '>';
-
-    typename dictionaryConstructorTable::iterator cstrIter =
-        dictionaryConstructorTablePtr_->find(methodTypeName);
-
-    if (cstrIter == dictionaryConstructorTablePtr_->end())
     {
-        FatalErrorInFunction
+        const dictionary& reductionDict(dict.subDict("reduction"));
-            << "Unknown " << typeName_() << " type " << methodName << endl
-            << endl;
 
-        const wordList names(dictionaryConstructorTablePtr_->sortedToc());
+        const word methodName(reductionDict.lookup("method"));
 
-        wordList thisCmpts;
+        Info<< "Selecting chemistry reduction method " << methodName << endl;
-        thisCmpts.append(word::null);
-        thisCmpts.append
-        (
-            basicThermo::splitThermoName(ThermoType::typeName(), 5)
-        );
 
-        wordList validNames;
+        const word methodTypeName =
-        forAll(names, i)
+            methodName + '<' + ThermoType::typeName() + '>';
+
+        typename dictionaryConstructorTable::iterator cstrIter =
+            dictionaryConstructorTablePtr_->find(methodTypeName);
+
+        if (cstrIter == dictionaryConstructorTablePtr_->end())
         {
-            const wordList cmpts(basicThermo::splitThermoName(names[i], 6));
+            FatalErrorInFunction
+                << "Unknown " << typeName_() << " type " << methodName << endl
+                << endl;
 
-            if (SubList<word>(cmpts, 5, 1) == SubList<word>(thisCmpts, 5, 1))
+            const wordList names(dictionaryConstructorTablePtr_->sortedToc());
+
+            wordList thisCmpts;
+            thisCmpts.append(word::null);
+            thisCmpts.append
+            (
+                basicThermo::splitThermoName(ThermoType::typeName(), 5)
+            );
+
+            wordList validNames;
+            forAll(names, i)
             {
-                validNames.append(cmpts[0]);
+                const wordList cmpts(basicThermo::splitThermoName(names[i], 6));
+
+                if
+                (
+                    SubList<word>(cmpts, 5, 1) == SubList<word>(thisCmpts, 5, 1)
+                )
+                {
+                    validNames.append(cmpts[0]);
+                }
             }
+
+            FatalErrorInFunction
+                << "Valid " << typeName_()
+                << " types are:" << validNames << endl
+                << exit(FatalError);
         }
 
-        FatalErrorInFunction
+        return autoPtr<chemistryReductionMethod<ThermoType>>
-            << "Valid " << typeName_() << " types are:" << validNames << endl
+        (
-            << exit(FatalError);
+            cstrIter()(dict, chemistry)
+        );
+    }
+    else
+    {
+        return autoPtr<chemistryReductionMethod<ThermoType>>
+        (
+            new chemistryReductionMethods::none<ThermoType>(dict, chemistry)
+        );
     }
-
-    return autoPtr<chemistryReductionMethod<ThermoType>>
-    (
-        cstrIter()(dict, chemistry)
-    );
 }
 
 

src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/noChemistryReduction/noChemistryReduction.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2016-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2016-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -35,9 +35,7 @@ Foam::chemistryReductionMethods::none<ThermoType>::none
 )
 :
     chemistryReductionMethod<ThermoType>(dict, chemistry)
-{
+{}
-    this->active_ = false;
-}
 
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //

src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/noChemistryReduction/noChemistryReduction.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2016-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2016-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,6 +25,8 @@ Class
     Foam::chemistryReductionMethods::none
 
 Description
+    A chemistryReductionMethod which does nothing to allow reduction to be
+    switched-off.
 
 SourceFiles
     noChemistryReduction.C
@@ -75,6 +77,12 @@ public:
 
     // Member Functions
 
+        //- Return mechanism reduction inactive
+        virtual bool active() const
+        {
+            return false;
+        }
+
         //- Reduce the mechanism
         virtual void reduceMechanism
         (
@@ -83,6 +91,9 @@ public:
             const scalarField& c,
             const label li
         );
+
+        virtual void update()
+        {}
 };
 
 

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

@@ -40,71 +40,95 @@ Foam::chemistryTabulationMethods::ISAT<ThermoType>::ISAT
         chemistryProperties,
         chemistry
     ),
-    chemisTree_(chemistry, this->coeffsDict_),
+    coeffsDict_(chemistryProperties.subDict("tabulation")),
+    chemistry_(chemistry),
+    log_(coeffsDict_.lookupOrDefault<Switch>("log", false)),
+    chemisTree_(*this, coeffsDict_),
     scaleFactor_(chemistry.nEqns() + 1, 1),
     runTime_(chemistry.time()),
+    timeSteps_(0),
     chPMaxLifeTime_
     (
-        this->coeffsDict_.lookupOrDefault("chPMaxLifeTime", INT_MAX)
+        coeffsDict_.lookupOrDefault("chPMaxLifeTime", INT_MAX)
     ),
-    maxGrowth_(this->coeffsDict_.lookupOrDefault("maxGrowth", INT_MAX)),
+    maxGrowth_(coeffsDict_.lookupOrDefault("maxGrowth", INT_MAX)),
     checkEntireTreeInterval_
     (
-        this->coeffsDict_.lookupOrDefault("checkEntireTreeInterval", INT_MAX)
+        coeffsDict_.lookupOrDefault("checkEntireTreeInterval", INT_MAX)
     ),
     maxDepthFactor_
     (
-        this->coeffsDict_.lookupOrDefault
+        coeffsDict_.lookupOrDefault
         (
             "maxDepthFactor",
             (chemisTree_.maxNLeafs() - 1)
-           /(log(scalar(chemisTree_.maxNLeafs()))/log(2.0))
+           /(Foam::log(scalar(chemisTree_.maxNLeafs()))/Foam::log(2.0))
         )
     ),
     minBalanceThreshold_
     (
-        this->coeffsDict_.lookupOrDefault
+        coeffsDict_.lookupOrDefault
         (
             "minBalanceThreshold", 0.1*chemisTree_.maxNLeafs()
         )
     ),
-    MRURetrieve_(this->coeffsDict_.lookupOrDefault("MRURetrieve", false)),
+    MRURetrieve_(coeffsDict_.lookupOrDefault("MRURetrieve", false)),
-    maxMRUSize_(this->coeffsDict_.lookupOrDefault("maxMRUSize", 0)),
+    maxMRUSize_(coeffsDict_.lookupOrDefault("maxMRUSize", 0)),
     lastSearch_(nullptr),
-    growPoints_(this->coeffsDict_.lookupOrDefault("growPoints", true)),
+    growPoints_(coeffsDict_.lookupOrDefault("growPoints", true)),
+    tolerance_(coeffsDict_.lookupOrDefault("tolerance", 1e-4)),
     nRetrieved_(0),
     nGrowth_(0),
     nAdd_(0),
+    addNewLeafCpuTime_(0),
+    growCpuTime_(0),
+    searchISATCpuTime_(0),
+    clockTime_(clockTime()),
+    tabulationResults_
+    (
+        IOobject
+        (
+            chemistry.thermo().phasePropertyName("TabulationResults"),
+            chemistry.time().timeName(),
+            chemistry.mesh(),
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        chemistry.mesh(),
+        scalar(0)
+    ),
+
     cleaningRequired_(false)
 {
-    if (this->active_)
+    dictionary scaleDict(coeffsDict_.subDict("scaleFactor"));
+    label Ysize = chemistry_.Y().size();
+    scalar otherScaleFactor = scaleDict.lookup<scalar>("otherSpecies");
+    for (label i=0; i<Ysize; i++)
     {
-        dictionary scaleDict(this->coeffsDict_.subDict("scaleFactor"));
+        if (!scaleDict.found(chemistry_.Y()[i].member()))
-        label Ysize = this->chemistry_.Y().size();
-        scalar otherScaleFactor = scaleDict.lookup<scalar>("otherSpecies");
-        for (label i=0; i<Ysize; i++)
         {
-            if (!scaleDict.found(this->chemistry_.Y()[i].member()))
+            scaleFactor_[i] = otherScaleFactor;
-            {
+        }
-                scaleFactor_[i] = otherScaleFactor;
+        else
-            }
+        {
-            else
+            scaleFactor_[i] =
-            {
+                scaleDict.lookup<scalar>(chemistry_.Y()[i].member());
-                scaleFactor_[i] =
-                    scaleDict.lookup<scalar>(this->chemistry_.Y()[i].member());
-            }
         }
-        scaleFactor_[Ysize] = scaleDict.lookup<scalar>("Temperature");
-        scaleFactor_[Ysize + 1] = scaleDict.lookup<scalar>("Pressure");
-        scaleFactor_[Ysize + 2] = scaleDict.lookup<scalar>("deltaT");
     }
+    scaleFactor_[Ysize] = scaleDict.lookup<scalar>("Temperature");
+    scaleFactor_[Ysize + 1] = scaleDict.lookup<scalar>("Pressure");
+    scaleFactor_[Ysize + 2] = scaleDict.lookup<scalar>("deltaT");
 
-    if (this->log())
+    if (log_)
     {
         nRetrievedFile_ = chemistry.logFile("found_isat.out");
         nGrowthFile_ = chemistry.logFile("growth_isat.out");
         nAddFile_ = chemistry.logFile("add_isat.out");
         sizeFile_ = chemistry.logFile("size_isat.out");
+
+        cpuAddFile_ = chemistry.logFile("cpu_add.out");
+        cpuGrowFile_ = chemistry.logFile("cpu_grow.out");
+        cpuRetrieveFile_ = chemistry.logFile("cpu_retrieve.out");
     }
 }
 
@@ -174,8 +198,8 @@ void Foam::chemistryTabulationMethods::ISAT<ThermoType>::calcNewC
     scalarField& Rphiq
 )
 {
-    label nEqns = this->chemistry_.nEqns(); // Species, T, p
+    label nEqns = chemistry_.nEqns(); // Species, T, p
-    bool mechRedActive = this->chemistry_.mechRed()->active();
+    bool mechRedActive = chemistry_.mechRed()->active();
     Rphiq = phi0->Rphi();
     scalarField dphi(phiq-phi0->phi());
     const scalarSquareMatrix& gradientsMatrix = phi0->A();
@@ -285,7 +309,7 @@ Foam::chemistryTabulationMethods::ISAT<ThermoType>::cleanAndBalance()
         chemPointISAT<ThermoType>* xtmp =
             chemisTree_.treeSuccessor(x);
 
-        scalar elapsedTimeSteps = this->chemistry_.timeSteps() - x->timeTag();
+        scalar elapsedTimeSteps = timeSteps() - x->timeTag();
 
         if ((elapsedTimeSteps > chPMaxLifeTime_) || (x->nGrowth() > maxGrowth_))
         {
@@ -304,7 +328,7 @@ Foam::chemistryTabulationMethods::ISAT<ThermoType>::cleanAndBalance()
     (
         chemisTree_.size() > minBalanceThreshold_
      && chemisTree_.depth() >
-        maxDepthFactor_*log(scalar(chemisTree_.size()))/log(2.0)
+        maxDepthFactor_*Foam::log(scalar(chemisTree_.size()))/Foam::log(2.0)
     )
     {
         chemisTree_.balance();
@@ -327,17 +351,17 @@ void Foam::chemistryTabulationMethods::ISAT<ThermoType>::computeA
     const scalar dt
 )
 {
-    bool mechRedActive = this->chemistry_.mechRed()->active();
+    bool mechRedActive = chemistry_.mechRed()->active();
-    label speciesNumber = this->chemistry_.nSpecie();
+    label speciesNumber = chemistry_.nSpecie();
-    scalarField Rcq(this->chemistry_.nEqns() + 1);
+    scalarField Rcq(chemistry_.nEqns() + 1);
     for (label i=0; i<speciesNumber; i++)
     {
         label s2c = i;
         if (mechRedActive)
         {
-            s2c = this->chemistry_.simplifiedToCompleteIndex()[i];
+            s2c = chemistry_.simplifiedToCompleteIndex()[i];
         }
-        Rcq[i] = rhoi*Rphiq[s2c]/this->chemistry_.specieThermos()[s2c].W();
+        Rcq[i] = rhoi*Rphiq[s2c]/chemistry_.specieThermos()[s2c].W();
     }
     Rcq[speciesNumber] = Rphiq[Rphiq.size() - 3];
     Rcq[speciesNumber + 1] = Rphiq[Rphiq.size() - 2];
@@ -353,7 +377,7 @@ void Foam::chemistryTabulationMethods::ISAT<ThermoType>::computeA
     // A = C(psi0,t0)/(I-dt*J(psi(t0+dt)))
     // where C(psi0,t0) = I
     scalarField dcdt(speciesNumber + 2, Zero);
-    this->chemistry_.jacobian(runTime_.value(), Rcq, li, dcdt, A);
+    chemistry_.jacobian(runTime_.value(), Rcq, li, dcdt, A);
 
     // The jacobian is computed according to the molar concentration
     // the following conversion allows the code to use A with mass fraction
@@ -363,7 +387,7 @@ void Foam::chemistryTabulationMethods::ISAT<ThermoType>::computeA
 
         if (mechRedActive)
         {
-            si = this->chemistry_.simplifiedToCompleteIndex()[i];
+            si = chemistry_.simplifiedToCompleteIndex()[i];
         }
 
         for (label j=0; j<speciesNumber; j++)
@@ -371,19 +395,19 @@ void Foam::chemistryTabulationMethods::ISAT<ThermoType>::computeA
             label sj = j;
             if (mechRedActive)
             {
-                sj = this->chemistry_.simplifiedToCompleteIndex()[j];
+                sj = chemistry_.simplifiedToCompleteIndex()[j];
             }
             A(i, j) *=
-              -dt*this->chemistry_.specieThermos()[si].W()
+              -dt*chemistry_.specieThermos()[si].W()
-               /this->chemistry_.specieThermos()[sj].W();
+               /chemistry_.specieThermos()[sj].W();
         }
 
         A(i, i) += 1;
         // Columns for pressure and temperature
         A(i, speciesNumber) *=
-            -dt*this->chemistry_.specieThermos()[si].W()/rhoi;
+            -dt*chemistry_.specieThermos()[si].W()/rhoi;
         A(i, speciesNumber + 1) *=
-            -dt*this->chemistry_.specieThermos()[si].W()/rhoi;
+            -dt*chemistry_.specieThermos()[si].W()/rhoi;
     }
 
     // For the temperature and pressure lines, ddc(dTdt)
@@ -393,13 +417,13 @@ void Foam::chemistryTabulationMethods::ISAT<ThermoType>::computeA
         label si = i;
         if (mechRedActive)
         {
-            si = this->chemistry_.simplifiedToCompleteIndex()[i];
+            si = chemistry_.simplifiedToCompleteIndex()[i];
         }
 
         A(speciesNumber, i) *=
-            -dt*rhoi/this->chemistry_.specieThermos()[si].W();
+            -dt*rhoi/chemistry_.specieThermos()[si].W();
         A(speciesNumber + 1, i) *=
-            -dt*rhoi/this->chemistry_.specieThermos()[si].W();
+            -dt*rhoi/chemistry_.specieThermos()[si].W();
     }
 
     A(speciesNumber, speciesNumber) = -dt*A(speciesNumber, speciesNumber) + 1;
@@ -433,6 +457,11 @@ bool Foam::chemistryTabulationMethods::ISAT<ThermoType>::retrieve
     scalarField& Rphiq
 )
 {
+    if (log_)
+    {
+        clockTime_.timeIncrement();
+    }
+
     bool retrieved(false);
     chemPointISAT<ThermoType>* phi0;
 
@@ -482,8 +511,7 @@ bool Foam::chemistryTabulationMethods::ISAT<ThermoType>::retrieve
     if (retrieved)
     {
         phi0->increaseNumRetrieve();
-        scalar elapsedTimeSteps =
+        scalar elapsedTimeSteps = timeSteps() - phi0->timeTag();
-            this->chemistry_.timeSteps() - phi0->timeTag();
 
         // Raise a flag when the chemPoint has been used more than the allowed
         // number of time steps
@@ -492,7 +520,7 @@ bool Foam::chemistryTabulationMethods::ISAT<ThermoType>::retrieve
             cleaningRequired_ = true;
             phi0->toRemove() = true;
         }
-        lastSearch_->lastTimeUsed() = this->chemistry_.timeSteps();
+        lastSearch_->lastTimeUsed() = timeSteps();
         addToMRU(phi0);
         calcNewC(phi0, phiq, Rphiq);
         nRetrieved_++;
@@ -504,6 +532,11 @@ bool Foam::chemistryTabulationMethods::ISAT<ThermoType>::retrieve
         // or if the tree is empty
         return false;
     }
+
+    if (log_)
+    {
+        searchISATCpuTime_ += clockTime_.timeIncrement();
+    }
 }
 
 
@@ -517,7 +550,13 @@ Foam::label Foam::chemistryTabulationMethods::ISAT<ThermoType>::add
     const scalar deltaT
 )
 {
+    if (log_)
+    {
+        clockTime_.timeIncrement();
+    }
+
     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_)
@@ -527,7 +566,15 @@ Foam::label Foam::chemistryTabulationMethods::ISAT<ThermoType>::add
             nGrowth_++;
             growthOrAddFlag = 0;
             addToMRU(lastSearch_);
-            //the structure of the tree is not modified, return false
+
+            tabulationResults_[li] = 1;
+
+            if (log_)
+            {
+                growCpuTime_ += clockTime_.timeIncrement();
+            }
+
+            // the structure of the tree is not modified, return false
             return growthOrAddFlag;
         }
     }
@@ -575,7 +622,7 @@ Foam::label Foam::chemistryTabulationMethods::ISAT<ThermoType>::add
                      tempList[i]->Rphi(),
                      tempList[i]->A(),
                      scaleFactor(),
-                     this->tolerance(),
+                     tolerance_,
                      scaleFactor_.size(),
                      nulPhi
                 );
@@ -589,7 +636,7 @@ Foam::label Foam::chemistryTabulationMethods::ISAT<ThermoType>::add
     }
 
     // Compute the A matrix needed to store the chemPoint.
-    label ASize = this->chemistry_.nEqns() + 1;
+    label ASize = chemistry_.nEqns() + 1;
     scalarSquareMatrix A(ASize, Zero);
     computeA(A, Rphiq, li, rho, deltaT);
 
@@ -599,7 +646,7 @@ Foam::label Foam::chemistryTabulationMethods::ISAT<ThermoType>::add
         Rphiq,
         A,
         scaleFactor(),
-        this->tolerance(),
+        tolerance_,
         scaleFactor_.size(),
         lastSearch_ // lastSearch_ may be nullptr (handled by binaryTree)
     );
@@ -609,6 +656,13 @@ Foam::label Foam::chemistryTabulationMethods::ISAT<ThermoType>::add
     }
     nAdd_++;
 
+    tabulationResults_[li] = 0;
+
+    if (log_)
+    {
+        addNewLeafCpuTime_ += clockTime_.timeIncrement();
+    }
+
     return growthOrAddFlag;
 }
 
@@ -616,7 +670,7 @@ Foam::label Foam::chemistryTabulationMethods::ISAT<ThermoType>::add
 template<class ThermoType>
 void Foam::chemistryTabulationMethods::ISAT<ThermoType>::writePerformance()
 {
-    if (this->log())
+    if (log_)
     {
         nRetrievedFile_()
             << runTime_.timeOutputValue() << "    " << nRetrieved_ << endl;
@@ -631,9 +685,47 @@ void Foam::chemistryTabulationMethods::ISAT<ThermoType>::writePerformance()
         nAdd_ = 0;
 
         sizeFile_()
-            << runTime_.timeOutputValue() << "    " << this->size() << endl;
+            << runTime_.timeOutputValue() << "    "
+            << chemisTree_.size() << endl;
+
+        cpuRetrieveFile_()
+            << runTime_.timeOutputValue()
+            << "    " << searchISATCpuTime_ << endl;
+        searchISATCpuTime_ = 0;
+
+        cpuGrowFile_()
+            << runTime_.timeOutputValue()
+            << "    " << growCpuTime_ << endl;
+        growCpuTime_ = 0;
+
+        cpuAddFile_()
+            << runTime_.timeOutputValue()
+            << "    " << addNewLeafCpuTime_ << endl;
+        addNewLeafCpuTime_ = 0;
     }
 }
 
 
+template<class ThermoType>
+void Foam::chemistryTabulationMethods::ISAT<ThermoType>::reset()
+{
+    // Increment counter of time-step
+    timeSteps_++;
+
+    forAll(tabulationResults_, i)
+    {
+        tabulationResults_[i] = 2;
+    }
+}
+
+
+template<class ThermoType>
+bool Foam::chemistryTabulationMethods::ISAT<ThermoType>::update()
+{
+    bool updated = cleanAndBalance();
+    writePerformance();
+    return updated;
+}
+
+
 // ************************************************************************* //

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

@@ -61,6 +61,13 @@ class ISAT
 {
     // Private Data
 
+        const dictionary coeffsDict_;
+
+        TDACChemistryModel<ThermoType>& chemistry_;
+
+        //- Switch to select performance logging
+        Switch log_;
+
         //- List of the stored 'points' organised in a binary tree
         binaryTree<ThermoType> chemisTree_;
 
@@ -69,6 +76,8 @@ class ISAT
 
         const Time& runTime_;
 
+        label timeSteps_;
+
         //- Lifetime (number of time steps) of a stored point
         label chPMaxLifeTime_;
 
@@ -100,16 +109,38 @@ class ISAT
         //- Switch to allow growth (on by default)
         Switch growPoints_;
 
+        scalar tolerance_;
+
         // Statistics on ISAT usage
         label nRetrieved_;
         label nGrowth_;
         label nAdd_;
+        scalar addNewLeafCpuTime_;
+        scalar growCpuTime_;
+        scalar searchISATCpuTime_;
+
+        const clockTime clockTime_;
 
         autoPtr<OFstream> nRetrievedFile_;
         autoPtr<OFstream> nGrowthFile_;
         autoPtr<OFstream> nAddFile_;
         autoPtr<OFstream> sizeFile_;
 
+        //- Log file for the average time spent adding tabulated data
+        autoPtr<OFstream> cpuAddFile_;
+
+        //- Log file for the average time spent growing tabulated data
+        autoPtr<OFstream> cpuGrowFile_;
+
+        //- Log file for the average time spent retrieving tabulated data
+        autoPtr<OFstream> cpuRetrieveFile_;
+
+        // Field containing information about tabulation:
+        // 0 -> add (direct integration)
+        // 1 -> grow
+        // 2 -> retrieve
+        volScalarField::Internal tabulationResults_;
+
         bool cleaningRequired_;
 
 
@@ -196,6 +227,17 @@ public:
 
     // Member Functions
 
+        //- Return true as this tabulation method tabulates
+        virtual bool tabulates()
+        {
+            return true;
+        }
+
+        TDACChemistryModel<ThermoType>& chemistry()
+        {
+            return chemistry_;
+        }
+
         inline binaryTree<ThermoType>& chemisTree()
         {
             return chemisTree_;
@@ -206,10 +248,10 @@ public:
             return scaleFactor_;
         }
 
-        //- Return the size of the binary tree
+        //- Return the number of chemistry evaluations
-        virtual inline label size()
+        inline label timeSteps() const
         {
-            return chemisTree_.size();
+            return timeSteps_;
         }
 
         virtual void writePerformance();
@@ -235,10 +277,9 @@ public:
             const scalar deltaT
         );
 
-        virtual bool update()
+        virtual void reset();
-        {
+
-            return cleanAndBalance();
+        virtual bool update();
-        }
 };
 
 

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

@@ -337,11 +337,11 @@ void Foam::binaryTree<ThermoType>::deleteAllNode(bn* subTreeRoot)
 template<class ThermoType>
 Foam::binaryTree<ThermoType>::binaryTree
 (
-    TDACChemistryModel<ThermoType>& chemistry,
+    chemistryTabulationMethods::ISAT<ThermoType>& table,
     dictionary coeffsDict
 )
 :
-    chemistry_(chemistry),
+    table_(table),
     root_(nullptr),
     maxNLeafs_(coeffsDict.lookup<label>("maxNLeafs")),
     size_(0),
@@ -394,7 +394,7 @@ void Foam::binaryTree<ThermoType>::insertNewLeaf
         chP* newChemPoint =
             new chP
             (
-                chemistry_,
+                table_,
                 phiq,
                 Rphiq,
                 A,
@@ -404,7 +404,7 @@ void Foam::binaryTree<ThermoType>::insertNewLeaf
                 coeffsDict_,
                 root_
             );
-        root_->leafLeft()=newChemPoint;
+        root_->leafLeft() = newChemPoint;
     }
     else // at least one point stored
     {
@@ -421,7 +421,7 @@ void Foam::binaryTree<ThermoType>::insertNewLeaf
         chP* newChemPoint =
             new chP
             (
-                chemistry_,
+                table_,
                 phiq,
                 Rphiq,
                 A,
@@ -634,7 +634,7 @@ void Foam::binaryTree<ThermoType>::deleteLeaf(chP*& phi0)
 template<class ThermoType>
 void Foam::binaryTree<ThermoType>::balance()
 {
-    scalarField mean(chemistry_.nEqns(),0.0);
+    scalarField mean(table_.chemistry().nEqns(), 0.0);
 
     //1) walk through the entire tree by starting with the tree's most left
     // chemPoint
@@ -652,7 +652,7 @@ void Foam::binaryTree<ThermoType>::balance()
     mean /= size_;
 
     //3) compute the variance for each space direction
-    List<scalar> variance(chemistry_.nEqns(),0.0);
+    List<scalar> variance(table_.chemistry().nEqns(),0.0);
     forAll(chemPoints, j)
     {
         const scalarField& phij = chemPoints[j]->phi();

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

@@ -51,8 +51,9 @@ Description
 namespace Foam
 {
 
-template<class ThermoType>
+/*---------------------------------------------------------------------------*\
-class TDACChemistryModel;
+                           Class binaryTree Declaration
+\*---------------------------------------------------------------------------*/
 
 template<class ThermoType>
 class binaryTree
@@ -69,8 +70,8 @@ private:
 
     // Private Data
 
-        //- Reference to the chemistryModel
+        //- Reference to the ISAT table
-        TDACChemistryModel<ThermoType>& chemistry_;
+        chemistryTabulationMethods::ISAT<ThermoType>& table_;
 
         //- Root node of the binary tree
         bn *root_;
@@ -134,7 +135,7 @@ public:
         //- Construct from dictionary and chemistryOnLineLibrary
         binaryTree
         (
-            TDACChemistryModel<ThermoType>& chemistry,
+            chemistryTabulationMethods::ISAT<ThermoType>& table,
             dictionary coeffsDict
         );
 

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

@@ -199,7 +199,7 @@ void Foam::chemPointISAT<ThermoType>::rotate
 template<class ThermoType>
 Foam::chemPointISAT<ThermoType>::chemPointISAT
 (
-    TDACChemistryModel<ThermoType>& chemistry,
+    chemistryTabulationMethods::ISAT<ThermoType>& table,
     const scalarField& phi,
     const scalarField& Rphi,
     const scalarSquareMatrix& A,
@@ -210,7 +210,7 @@ Foam::chemPointISAT<ThermoType>::chemPointISAT
     binaryNode<ThermoType>* node
 )
 :
-    chemistry_(chemistry),
+    table_(table),
     phi_(phi),
     Rphi_(Rphi),
     A_(A),
@@ -218,10 +218,10 @@ Foam::chemPointISAT<ThermoType>::chemPointISAT
     node_(node),
     completeSpaceSize_(completeSpaceSize),
     nGrowth_(0),
-    nActiveSpecies_(chemistry.mechRed()->NsSimp()),
+    nActiveSpecies_(table_.chemistry().mechRed()->nActiveSpecies()),
     simplifiedToCompleteIndex_(nActiveSpecies_),
-    timeTag_(chemistry_.timeSteps()),
+    timeTag_(table.timeSteps()),
-    lastTimeUsed_(chemistry_.timeSteps()),
+    lastTimeUsed_(table.timeSteps()),
     toRemove_(false),
     maxNumNewDim_(coeffsDict.lookupOrDefault("maxNumNewDim",0)),
     printProportion_(coeffsDict.lookupOrDefault("printProportion",false)),
@@ -237,18 +237,18 @@ Foam::chemPointISAT<ThermoType>::chemPointISAT
     idT_ = completeSpaceSize - 3;
     idp_ = completeSpaceSize - 2;
 
-    bool isMechRedActive = chemistry_.mechRed()->active();
+    bool isMechRedActive = table_.chemistry().mechRed()->active();
     if (isMechRedActive)
     {
         for (label i=0; i<completeSpaceSize-3; i++)
         {
             completeToSimplifiedIndex_[i] =
-                chemistry.completeToSimplifiedIndex()[i];
+                table_.chemistry().completeToSimplifiedIndex()[i];
         }
         for (label i=0; i<nActiveSpecies_; i++)
         {
             simplifiedToCompleteIndex_[i] =
-                chemistry.simplifiedToCompleteIndex()[i];
+                table_.chemistry().simplifiedToCompleteIndex()[i];
         }
     }
 
@@ -307,7 +307,7 @@ Foam::chemPointISAT<ThermoType>::chemPointISAT
     Foam::chemPointISAT<ThermoType>& p
 )
 :
-    chemistry_(p.chemistry()),
+    table_(p.table_),
     phi_(p.phi()),
     Rphi_(p.Rphi()),
     LT_(p.LT()),
@@ -340,7 +340,7 @@ template<class ThermoType>
 bool Foam::chemPointISAT<ThermoType>::inEOA(const scalarField& phiq)
 {
     scalarField dphi(phiq-phi());
-    bool isMechRedActive = chemistry_.mechRed()->active();
+    bool isMechRedActive = table_.chemistry().mechRed()->active();
 
     const label dim =
         isMechRedActive ? nActiveSpecies_ : completeSpaceSize() - 3;
@@ -452,7 +452,7 @@ bool Foam::chemPointISAT<ThermoType>::inEOA(const scalarField& phiq)
             }
             else
             {
-                propName = chemistry_.Y()[maxIndex].member();
+                propName = table_.chemistry().Y()[maxIndex].member();
             }
             Info<< "Direction maximum impact to error in ellipsoid: "
                 << propName << endl;
@@ -480,7 +480,7 @@ bool Foam::chemPointISAT<ThermoType>::checkSolution
     scalarField dphi(phiq - phi());
     const scalarField& scaleFactorV(scaleFactor());
     const scalarSquareMatrix& Avar(A());
-    bool isMechRedActive = chemistry_.mechRed()->active();
+    bool isMechRedActive = table_.chemistry().mechRed()->active();
     scalar dRl = 0;
 
     const label dim =
@@ -540,7 +540,7 @@ bool Foam::chemPointISAT<ThermoType>::grow(const scalarField& phiq)
 {
     scalarField dphi(phiq - phi());
     label initNActiveSpecies(nActiveSpecies_);
-    bool isMechRedActive = chemistry_.mechRed()->active();
+    bool isMechRedActive = table_.chemistry().mechRed()->active();
 
     if (isMechRedActive)
     {
@@ -556,7 +556,7 @@ bool Foam::chemPointISAT<ThermoType>::grow(const scalarField& phiq)
             if
             (
                 completeToSimplifiedIndex_[i] == -1
-             && chemistry_.completeToSimplifiedIndex()[i]!=-1
+             && table_.chemistry().completeToSimplifiedIndex()[i]!=-1
             )
             {
                 activeAdded++;
@@ -567,7 +567,7 @@ bool Foam::chemPointISAT<ThermoType>::grow(const scalarField& phiq)
             if
             (
                 completeToSimplifiedIndex_[i] != -1
-             && chemistry_.completeToSimplifiedIndex()[i] == -1
+             && table_.chemistry().completeToSimplifiedIndex()[i] == -1
             )
             {
                 activeAdded++;
@@ -579,7 +579,7 @@ bool Foam::chemPointISAT<ThermoType>::grow(const scalarField& phiq)
             if
             (
                 completeToSimplifiedIndex_[i] == -1
-             && chemistry_.completeToSimplifiedIndex()[i] == -1
+             && table_.chemistry().completeToSimplifiedIndex()[i] == -1
              && dphi[i] != 0
             )
             {

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

@@ -135,6 +135,12 @@ class binaryNode;
 template<class ThermoType>
 class TDACChemistryModel;
 
+namespace chemistryTabulationMethods
+{
+    template<class ThermoType>
+    class ISAT;
+}
+
 
 /*---------------------------------------------------------------------------*\
                         Class chemPointISAT Declaration
@@ -145,8 +151,8 @@ class chemPointISAT
 {
     // Private Data
 
-        //- Pointer to the chemistryModel object
+        //- Reference to the ISAT table
-        TDACChemistryModel<ThermoType>& chemistry_;
+        chemistryTabulationMethods::ISAT<ThermoType>& table_;
 
         //- Vector storing the composition, temperature and pressure
         //  and deltaT if a variable time step is set on
@@ -245,7 +251,7 @@ public:
         //- Construct from components
         chemPointISAT
         (
-            TDACChemistryModel<ThermoType>& chemistry,
+            chemistryTabulationMethods::ISAT<ThermoType>& table,
             const scalarField& phi,
             const scalarField& Rphi,
             const scalarSquareMatrix& A,
@@ -272,7 +278,7 @@ public:
         //- Access to the TDACChemistryModel
         inline TDACChemistryModel<ThermoType>& chemistry()
         {
-            return chemistry_;
+            return table_.chemistry();
         }
 
         inline label nGrowth()

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2016-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2016-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -34,13 +34,6 @@ Foam::chemistryTabulationMethod<ThermoType>::chemistryTabulationMethod
     const dictionary& dict,
     TDACChemistryModel<ThermoType>& chemistry
 )
-:
-    dict_(dict),
-    coeffsDict_(dict.subDict("tabulation")),
-    active_(coeffsDict_.lookupOrDefault<Switch>("active", false)),
-    log_(coeffsDict_.lookupOrDefault<Switch>("log", false)),
-    chemistry_(chemistry),
-    tolerance_(coeffsDict_.lookupOrDefault<scalar>("tolerance", 1e-4))
 {}
 
 

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

@@ -37,7 +37,6 @@ SourceFiles
 
 #include "IOdictionary.H"
 #include "scalarField.H"
-#include "Switch.H"
 #include "runTimeSelectionTables.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -58,23 +57,6 @@ template<class ThermoType>
 class chemistryTabulationMethod
 {
 
-protected:
-
-    const dictionary& dict_;
-
-    const dictionary coeffsDict_;
-
-    //- Is tabulation active?
-    Switch active_;
-
-    //- Switch to select performance logging
-    Switch log_;
-
-    TDACChemistryModel<ThermoType>& chemistry_;
-
-    scalar tolerance_;
-
-
 public:
 
     //- Runtime type information
@@ -120,24 +102,9 @@ public:
 
     // Member Functions
 
-        inline bool active()
+        //- Return true if the tabulation method tabulates
-        {
+        //  otherwise return false
-            return active_;
+        virtual bool tabulates() = 0;
-        }
-
-        inline bool log()
-        {
-            return active_ && log_;
-        }
-
-        inline scalar tolerance() const
-        {
-            return tolerance_;
-        }
-
-        virtual label size() = 0;
-
-        virtual void writePerformance() = 0;
 
         // Retrieve function: (only virtual here)
         // Try to retrieve a stored point close enough (according to tolerance)
@@ -165,6 +132,8 @@ public:
         // The underlying structure of the tabulation is updated/cleaned
         // to increase the performance of the retrieve
         virtual bool update() = 0;
+
+        virtual void reset() = 0;
 };
 
 

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2016-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2016-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -24,6 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "chemistryTabulationMethod.H"
+#include "noChemistryTabulation.H"
 #include "Time.H"
 
 // * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * //
@@ -36,53 +37,67 @@ Foam::chemistryTabulationMethod<ThermoType>::New
     TDACChemistryModel<ThermoType>& chemistry
 )
 {
-    const dictionary& tabulationDict(dict.subDict("tabulation"));
+    if (dict.found("tabulation"))
-
-    const word methodName(tabulationDict.lookup("method"));
-
-    Info<< "Selecting chemistry tabulation method " << methodName << endl;
-
-    const word methodTypeName =
-        methodName + '<' + ThermoType::typeName() + '>';
-
-    typename dictionaryConstructorTable::iterator cstrIter =
-        dictionaryConstructorTablePtr_->find(methodTypeName);
-
-    if (cstrIter == dictionaryConstructorTablePtr_->end())
     {
-        FatalErrorInFunction
+        const dictionary& tabulationDict(dict.subDict("tabulation"));
-            << "Unknown " << typeName_() << " type " << methodName << endl
-            << endl;
 
-        const wordList names(dictionaryConstructorTablePtr_->sortedToc());
+        const word methodName(tabulationDict.lookup("method"));
 
-        wordList thisCmpts;
+        Info<< "Selecting chemistry tabulation method " << methodName << endl;
-        thisCmpts.append(word::null);
-        thisCmpts.append
-        (
-            basicThermo::splitThermoName(ThermoType::typeName(), 5)
-        );
 
-        wordList validNames;
+        const word methodTypeName =
-        forAll(names, i)
+            methodName + '<' + ThermoType::typeName() + '>';
+
+        typename dictionaryConstructorTable::iterator cstrIter =
+            dictionaryConstructorTablePtr_->find(methodTypeName);
+
+        if (cstrIter == dictionaryConstructorTablePtr_->end())
         {
-            const wordList cmpts(basicThermo::splitThermoName(names[i], 6));
+            FatalErrorInFunction
+                << "Unknown " << typeName_() << " type " << methodName << endl
+                << endl;
 
-            if (SubList<word>(cmpts, 5, 1) == SubList<word>(thisCmpts, 5, 1))
+            const wordList names(dictionaryConstructorTablePtr_->sortedToc());
+
+            wordList thisCmpts;
+            thisCmpts.append(word::null);
+            thisCmpts.append
+            (
+                basicThermo::splitThermoName(ThermoType::typeName(), 5)
+            );
+
+            wordList validNames;
+            forAll(names, i)
             {
-                validNames.append(cmpts[0]);
+                const wordList cmpts(basicThermo::splitThermoName(names[i], 6));
+
+                if
+                (
+                    SubList<word>(cmpts, 5, 1) == SubList<word>(thisCmpts, 5, 1)
+                )
+                {
+                    validNames.append(cmpts[0]);
+                }
             }
+
+            FatalErrorInFunction
+                << "Valid " << typeName_()
+                << " types are:" << validNames << endl
+                << exit(FatalError);
         }
 
-        FatalErrorInFunction
+        return autoPtr<chemistryTabulationMethod<ThermoType>>
-            << "Valid " << typeName_() << " types are:" << validNames << endl
+        (
-            << exit(FatalError);
+            cstrIter()(dict, chemistry)
+        );
+    }
+    else
+    {
+        return autoPtr<chemistryTabulationMethod<ThermoType>>
+        (
+            new chemistryTabulationMethods::none<ThermoType>(dict, chemistry)
+        );
     }
-
-    return autoPtr<chemistryTabulationMethod<ThermoType>>
-    (
-        cstrIter()(dict, chemistry)
-    );
 }
 
 

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2016-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2016-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -39,9 +39,7 @@ Foam::chemistryTabulationMethods::none<ThermoType>::none
         chemistryProperties,
         chemistry
     )
-{
+{}
-    this->active_ = false;
-}
 
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2016-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2016-2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,6 +25,11 @@ Class
     Foam::chemistryTabulationMethods::none
 
 Description
+    A chemistryTabulationMethod which does nothing to allow tabulation to be
+    switched-off.
+
+SourceFiles
+    noChemistryTabulation.C
 
 \*---------------------------------------------------------------------------*/
 
@@ -73,16 +78,10 @@ public:
 
     // Member Functions
 
-        //- Return the size of the binary tree
+        //- Return false as this tabulation method does not tabulate
-        virtual label size()
+        virtual bool tabulates()
         {
-            NotImplemented;
+            return false;
-            return 0;
-        }
-
-        virtual void writePerformance()
-        {
-            NotImplemented;
         }
 
         //- Find the closest stored leaf of phiQ and store the result in
@@ -93,7 +92,6 @@ public:
             scalarField& Rphiq
         )
         {
-            NotImplemented;
             return false;
         }
 
@@ -110,15 +108,16 @@ public:
             const scalar deltaT
         )
         {
-            NotImplemented;
             return false;
         }
 
         virtual bool update()
         {
-            NotImplemented;
             return false;
         }
+
+        virtual void reset()
+        {}
 };