diff --git a/etc/caseDicts/solvers/chemistry/TDAC/chemistryProperties.cfg b/etc/caseDicts/solvers/chemistry/TDAC/chemistryProperties.cfg
index 4ba574f81a..befb1a76dc 100644
--- a/etc/caseDicts/solvers/chemistry/TDAC/chemistryProperties.cfg
+++ b/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;
 
diff --git a/etc/caseDicts/solvers/chemistry/TDAC/chemistryPropertiesFlame.cfg b/etc/caseDicts/solvers/chemistry/TDAC/chemistryPropertiesFlame.cfg
index 03f4018327..a3d5e2d47a 100644
--- a/etc/caseDicts/solvers/chemistry/TDAC/chemistryPropertiesFlame.cfg
+++ b/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;
 
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModel.C b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModel.C
index 8749e7aeee..6d76516c4f 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModel.C
+++ b/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_),
-    tabulationResults_
-    (
-        IOobject
-        (
-            thermo.phasePropertyName("TabulationResults"),
-            this->time().timeName(),
-            this->mesh(),
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        this->mesh(),
-        scalar(0)
-    )
+    simplifiedToCompleteIndex_(this->nSpecie_)
 {
     const basicSpecieMixture& composition = this->thermo().composition();
 
@@ -94,20 +81,7 @@ Foam::TDACChemistryModel<ThermoType>::TDACChemistryModel
         *this
     );
 
-    if (mechRed_->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())
+    if (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]);
-                if (growOrAdd)
-                {
-                    this->setTabulationResultsAdd(celli);
-                    addNewLeafCpuTime_ += clockTime_.timeIncrement();
-                }
-                else
-                {
-                    this->setTabulationResultsGrow(celli);
-                    growCpuTime_ += clockTime_.timeIncrement();
-                }
+                tabulation_->add(phiq, Rphiq, celli, rhoi, deltaT[celli]);
             }
 
             // 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())
-    {
-        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;
-    }
+    mechRed_->update();
+    tabulation_->update();
 
     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;
-}
-
-
 // ************************************************************************* //
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModel.H b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModel.H
index a423ddcc17..917dfce7ab 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModel.H
+++ b/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
 
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModelI.H b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModelI.H
index d045ec6490..44e51422a9 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModelI.H
+++ b/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;
-    }
-}
-
-
 // ************************************************************************* //
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DAC/DAC.C b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DAC/DAC.C
index f4ae85d21c..12329ad3d7 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DAC/DAC.C
+++ b/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),
-    sH_(this->nSpecie_,0),
-    sO_(this->nSpecie_,0),
-    sN_(this->nSpecie_,0),
+    sC_(this->nSpecie(),0),
+    sH_(this->nSpecie(),0),
+    sO_(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_+1] = p;
+    c1[this->nSpecie()] = T;
+    c1[this->nSpecie()+1] = p;
 
     // Compute the rAB matrix
-    RectangularMatrix<scalar> rABNum(this->nSpecie_,this->nSpecie_,0.0);
-    scalarField PA(this->nSpecie_,0.0);
-    scalarField CA(this->nSpecie_,0.0);
+    RectangularMatrix<scalar> rABNum(this->nSpecie(),this->nSpecie(),0.0);
+    scalarField PA(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->NsSimp_+1] = p;
-    this->chemistry_.setNsDAC(this->NsSimp_);
+    simplifiedC[this->nActiveSpecies_] = T;
+    simplifiedC[this->nActiveSpecies_+1] = p;
+    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();
 }
 
 
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DAC/DAC.H b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DAC/DAC.H
index afc0e9c659..0a726a3890 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DAC/DAC.H
+++ b/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
 
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DRG/DRG.C b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DRG/DRG.C
index b37f5fd358..3918d8c828 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DRG/DRG.C
+++ b/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_+1] = p;
+    c1[this->nSpecie()] = T;
+    c1[this->nSpecie()+1] = p;
 
     // Compute the rAB matrix
-    RectangularMatrix<scalar> rABNum(this->nSpecie_,this->nSpecie_,0.0);
-    scalarField rABDen(this->nSpecie_,0.0);
+    RectangularMatrix<scalar> rABNum(this->nSpecie(),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->chemistry_.simplifiedC().setSize(this->NsSimp_+2);
-    this->chemistry_.simplifiedToCompleteIndex().setSize(this->NsSimp_);
+    this->nActiveSpecies_ = speciesNumber;
+    this->chemistry_.simplifiedC().setSize(this->nActiveSpecies_+2);
+    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->NsSimp_+1] = p;
-    this->chemistry_.setNsDAC(this->NsSimp_);
+    this->chemistry_.simplifiedC()[this->nActiveSpecies_] = T;
+    this->chemistry_.simplifiedC()[this->nActiveSpecies_+1] = p;
+    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();
 }
 
 
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DRG/DRG.H b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DRG/DRG.H
index 30ca47c028..5cc92a4d69 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DRG/DRG.H
+++ b/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
 
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DRGEP/DRGEP.C b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DRGEP/DRGEP.C
index 95571643f9..ed3979399e 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DRGEP/DRGEP.C
+++ b/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),
-    sH_(this->nSpecie_,0),
-    sO_(this->nSpecie_,0),
-    sN_(this->nSpecie_,0),
+    sC_(this->nSpecie(), 0),
+    sH_(this->nSpecie(), 0),
+    sO_(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_+1] = p;
+    c1[this->nSpecie()] = T;
+    c1[this->nSpecie()+1] = p;
 
     // Compute the rAB matrix
-    RectangularMatrix<scalar> rABNum(this->nSpecie_,this->nSpecie_,0.0);
-    scalarField PA(this->nSpecie_,0.0);
-    scalarField CA(this->nSpecie_,0.0);
+    RectangularMatrix<scalar> rABNum(this->nSpecie(),this->nSpecie(),0.0);
+    scalarField PA(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->NsSimp_+1] = p;
-    this->chemistry_.setNsDAC(this->NsSimp_);
+    simplifiedC[this->nActiveSpecies_] = T;
+    simplifiedC[this->nActiveSpecies_+1] = p;
+    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();
 }
 
 
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DRGEP/DRGEP.H b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DRGEP/DRGEP.H
index ebc23def91..66d2935161 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DRGEP/DRGEP.H
+++ b/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
 
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/EFA/EFA.C b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/EFA/EFA.C
index 40c5d750db..8bc89747c1 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/EFA/EFA.C
+++ b/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),
-    sC_(this->nSpecie_,0),
-    sH_(this->nSpecie_,0),
-    sO_(this->nSpecie_,0),
-    sN_(this->nSpecie_,0),
+    chemistryReduction<ThermoType>(dict, chemistry),
+    sC_(this->nSpecie(),0),
+    sH_(this->nSpecie(),0),
+    sO_(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_+1] = p;
+    c1[this->nSpecie()] = T;
+    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<scalar> CFluxAB(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> NFluxAB(this->nSpecie_, this->nSpecie_, 0.0);
+    RectangularMatrix<label> rABPos(this->nSpecie(), this->nSpecie(), -1);
+    RectangularMatrix<scalar> CFluxAB(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> 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->NsSimp_+1] = p;
-    this->chemistry_.setNsDAC(this->NsSimp_);
+    simplifiedC[this->nActiveSpecies_] = T;
+    simplifiedC[this->nActiveSpecies_+1] = p;
+    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();
 }
 
 
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/EFA/EFA.H b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/EFA/EFA.H
index adad83f5a7..8505ed567b 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/EFA/EFA.H
+++ b/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
 
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/PFA/PFA.C b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/PFA/PFA.C
index 84dd1fa12b..f9969bb1e5 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/PFA/PFA.C
+++ b/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_+1] = p;
+    c1[this->nSpecie()] = T;
+    c1[this->nSpecie()+1] = p;
 
     // Compute the rAB matrix
-    RectangularMatrix<scalar> PAB(this->nSpecie_,this->nSpecie_,0.0);
-    RectangularMatrix<scalar> CAB(this->nSpecie_,this->nSpecie_,0.0);
-    scalarField PA(this->nSpecie_,0.0);
-    scalarField CA(this->nSpecie_,0.0);
+    RectangularMatrix<scalar> PAB(this->nSpecie(),this->nSpecie(),0.0);
+    RectangularMatrix<scalar> CAB(this->nSpecie(),this->nSpecie(),0.0);
+    scalarField PA(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> CAB2nd(this->nSpecie_,this->nSpecie_,0.0);
+    RectangularMatrix<scalar> PAB2nd(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->NsSimp_+1] = p;
-    this->chemistry_.setNsDAC(this->NsSimp_);
+    simplifiedC[this->nActiveSpecies_] = T;
+    simplifiedC[this->nActiveSpecies_+1] = p;
+    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();
 }
 
 
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/PFA/PFA.H b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/PFA/PFA.H
index 47d3b99c6c..1db17e4411 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/PFA/PFA.H
+++ b/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
 
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/chemistryReduction/chemistryReduction.C b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/chemistryReduction/chemistryReduction.C
new file mode 100644
index 0000000000..066c12ad15
--- /dev/null
+++ b/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;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/chemistryReduction/chemistryReduction.H b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/chemistryReduction/chemistryReduction.H
new file mode 100644
index 0000000000..628eb671c4
--- /dev/null
+++ b/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
+
+// ************************************************************************* //
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/chemistryReduction/chemistryReductionI.H b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/chemistryReduction/chemistryReductionI.H
new file mode 100644
index 0000000000..f06c9eca12
--- /dev/null
+++ b/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_;
+}
+
+
+// ************************************************************************* //
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/chemistryReductionMethod/chemistryReductionMethod.C b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/chemistryReductionMethod/chemistryReductionMethod.C
index 3e13f860e6..19b437662b 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/chemistryReductionMethod/chemistryReductionMethod.C
+++ b/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())
 {}
 
 
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/chemistryReductionMethod/chemistryReductionMethod.H b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/chemistryReductionMethod/chemistryReductionMethod.H
index 25cac1d1fe..d6ae5f852b 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/chemistryReductionMethod/chemistryReductionMethod.H
+++ b/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
-    const dictionary coeffsDict_;
-
-    //- 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_;
+        //- Number of active species
+        label nActiveSpecies_;
 
 
 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?
-        inline bool log() const;
-
-        //- Return the active species
-        inline const List<bool>& activeSpecies() const;
+        //- Return the number of species
+        inline label nSpecie();
 
         //- Return the number of active species
-        inline label NsSimp();
-
-        //- Return the initial number of species
-        inline label nSpecie();
+        inline label nActiveSpecies();
 
         //- 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
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/chemistryReductionMethod/chemistryReductionMethodI.H b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/chemistryReductionMethod/chemistryReductionMethodI.H
index 7be38b7184..ccfc826124 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/chemistryReductionMethod/chemistryReductionMethodI.H
+++ b/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_;
-}
-
-
-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_;
+    return nActiveSpecies_;
 }
 
 
@@ -69,12 +48,4 @@ inline Foam::label Foam::chemistryReductionMethod<ThermoType>::nSpecie()
 }
 
 
-template<class ThermoType>
-inline Foam::scalar
-Foam::chemistryReductionMethod<ThermoType>::tolerance() const
-{
-    return tolerance_;
-}
-
-
 // ************************************************************************* //
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/chemistryReductionMethod/chemistryReductionMethodNew.C b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/chemistryReductionMethod/chemistryReductionMethodNew.C
index dd22297b3d..7aa52f5de7 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/chemistryReductionMethod/chemistryReductionMethodNew.C
+++ b/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"));
-
-    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())
+    if (dict.found("reduction"))
     {
-        FatalErrorInFunction
-            << "Unknown " << typeName_() << " type " << methodName << endl
-            << endl;
+        const dictionary& reductionDict(dict.subDict("reduction"));
 
-        const wordList names(dictionaryConstructorTablePtr_->sortedToc());
+        const word methodName(reductionDict.lookup("method"));
 
-        wordList thisCmpts;
-        thisCmpts.append(word::null);
-        thisCmpts.append
-        (
-            basicThermo::splitThermoName(ThermoType::typeName(), 5)
-        );
+        Info<< "Selecting chemistry reduction method " << methodName << endl;
 
-        wordList validNames;
-        forAll(names, i)
+        const word methodTypeName =
+            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
-            << "Valid " << typeName_() << " types are:" << validNames << endl
-            << exit(FatalError);
+        return autoPtr<chemistryReductionMethod<ThermoType>>
+        (
+            cstrIter()(dict, chemistry)
+        );
+    }
+    else
+    {
+        return autoPtr<chemistryReductionMethod<ThermoType>>
+        (
+            new chemistryReductionMethods::none<ThermoType>(dict, chemistry)
+        );
     }
-
-    return autoPtr<chemistryReductionMethod<ThermoType>>
-    (
-        cstrIter()(dict, chemistry)
-    );
 }
 
 
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/noChemistryReduction/noChemistryReduction.C b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/noChemistryReduction/noChemistryReduction.C
index c6800f4eb1..636c4b80f3 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/noChemistryReduction/noChemistryReduction.C
+++ b/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  * * * * * * * * * * * * * * * //
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/noChemistryReduction/noChemistryReduction.H b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/noChemistryReduction/noChemistryReduction.H
index 821c5fbf53..b20bd84894 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/noChemistryReduction/noChemistryReduction.H
+++ b/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()
+        {}
 };
 
 
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/ISAT/ISAT.C b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/ISAT/ISAT.C
index 05d9b36978..5d7a744c6a 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/ISAT/ISAT.C
+++ b/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)),
-    maxMRUSize_(this->coeffsDict_.lookupOrDefault("maxMRUSize", 0)),
+    MRURetrieve_(coeffsDict_.lookupOrDefault("MRURetrieve", false)),
+    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"));
-        label Ysize = this->chemistry_.Y().size();
-        scalar otherScaleFactor = scaleDict.lookup<scalar>("otherSpecies");
-        for (label i=0; i<Ysize; i++)
+        if (!scaleDict.found(chemistry_.Y()[i].member()))
         {
-            if (!scaleDict.found(this->chemistry_.Y()[i].member()))
-            {
-                scaleFactor_[i] = otherScaleFactor;
-            }
-            else
-            {
-                scaleFactor_[i] =
-                    scaleDict.lookup<scalar>(this->chemistry_.Y()[i].member());
-            }
+            scaleFactor_[i] = otherScaleFactor;
+        }
+        else
+        {
+            scaleFactor_[i] =
+                scaleDict.lookup<scalar>(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
-    bool mechRedActive = this->chemistry_.mechRed()->active();
+    label nEqns = chemistry_.nEqns(); // Species, T, p
+    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();
-    label speciesNumber = this->chemistry_.nSpecie();
-    scalarField Rcq(this->chemistry_.nEqns() + 1);
+    bool mechRedActive = chemistry_.mechRed()->active();
+    label speciesNumber = chemistry_.nSpecie();
+    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()
-               /this->chemistry_.specieThermos()[sj].W();
+              -dt*chemistry_.specieThermos()[si].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 =
-            this->chemistry_.timeSteps() - phi0->timeTag();
+        scalar elapsedTimeSteps = 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;
+}
+
+
 // ************************************************************************* //
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/ISAT/ISAT.H b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/ISAT/ISAT.H
index 5cf5a62b26..a7a6d5d15b 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/ISAT/ISAT.H
+++ b/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
-        virtual inline label size()
+        //- Return the number of chemistry evaluations
+        inline label timeSteps() const
         {
-            return chemisTree_.size();
+            return timeSteps_;
         }
 
         virtual void writePerformance();
@@ -235,10 +277,9 @@ public:
             const scalar deltaT
         );
 
-        virtual bool update()
-        {
-            return cleanAndBalance();
-        }
+        virtual void reset();
+
+        virtual bool update();
 };
 
 
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/ISAT/binaryTree/binaryTree.C b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/ISAT/binaryTree/binaryTree.C
index ba5a3c6bd1..9095c2b12c 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/ISAT/binaryTree/binaryTree.C
+++ b/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();
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/ISAT/binaryTree/binaryTree.H b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/ISAT/binaryTree/binaryTree.H
index e8e55b5a03..d75ec0bcf3 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/ISAT/binaryTree/binaryTree.H
+++ b/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
-        TDACChemistryModel<ThermoType>& chemistry_;
+        //- Reference to the ISAT table
+        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
         );
 
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/ISAT/chemPointISAT/chemPointISAT.C b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/ISAT/chemPointISAT/chemPointISAT.C
index 88b99dc4bb..383ccc4e2f 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/ISAT/chemPointISAT/chemPointISAT.C
+++ b/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()),
-    lastTimeUsed_(chemistry_.timeSteps()),
+    timeTag_(table.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
             )
             {
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/ISAT/chemPointISAT/chemPointISAT.H b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/ISAT/chemPointISAT/chemPointISAT.H
index c7c2ff5661..513a741fbc 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/ISAT/chemPointISAT/chemPointISAT.H
+++ b/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
-        TDACChemistryModel<ThermoType>& chemistry_;
+        //- Reference to the ISAT table
+        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()
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/chemistryTabulationMethod/chemistryTabulationMethod.C b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/chemistryTabulationMethod/chemistryTabulationMethod.C
index 48bf2f7087..8e64853621 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/chemistryTabulationMethod/chemistryTabulationMethod.C
+++ b/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))
 {}
 
 
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/chemistryTabulationMethod/chemistryTabulationMethod.H b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/chemistryTabulationMethod/chemistryTabulationMethod.H
index 0084dd4cad..38606df9ba 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/chemistryTabulationMethod/chemistryTabulationMethod.H
+++ b/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 active_;
-        }
-
-        inline bool log()
-        {
-            return active_ && log_;
-        }
-
-        inline scalar tolerance() const
-        {
-            return tolerance_;
-        }
-
-        virtual label size() = 0;
-
-        virtual void writePerformance() = 0;
+        //- Return true if the tabulation method tabulates
+        //  otherwise return false
+        virtual bool tabulates() = 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;
 };
 
 
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/chemistryTabulationMethod/chemistryTabulationMethodNew.C b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/chemistryTabulationMethod/chemistryTabulationMethodNew.C
index 8e4680f162..4033a63f82 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/chemistryTabulationMethod/chemistryTabulationMethodNew.C
+++ b/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"));
-
-    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())
+    if (dict.found("tabulation"))
     {
-        FatalErrorInFunction
-            << "Unknown " << typeName_() << " type " << methodName << endl
-            << endl;
+        const dictionary& tabulationDict(dict.subDict("tabulation"));
 
-        const wordList names(dictionaryConstructorTablePtr_->sortedToc());
+        const word methodName(tabulationDict.lookup("method"));
 
-        wordList thisCmpts;
-        thisCmpts.append(word::null);
-        thisCmpts.append
-        (
-            basicThermo::splitThermoName(ThermoType::typeName(), 5)
-        );
+        Info<< "Selecting chemistry tabulation method " << methodName << endl;
 
-        wordList validNames;
-        forAll(names, i)
+        const word methodTypeName =
+            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
-            << "Valid " << typeName_() << " types are:" << validNames << endl
-            << exit(FatalError);
+        return autoPtr<chemistryTabulationMethod<ThermoType>>
+        (
+            cstrIter()(dict, chemistry)
+        );
+    }
+    else
+    {
+        return autoPtr<chemistryTabulationMethod<ThermoType>>
+        (
+            new chemistryTabulationMethods::none<ThermoType>(dict, chemistry)
+        );
     }
-
-    return autoPtr<chemistryTabulationMethod<ThermoType>>
-    (
-        cstrIter()(dict, chemistry)
-    );
 }
 
 
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/noChemistryTabulation/noChemistryTabulation.C b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/noChemistryTabulation/noChemistryTabulation.C
index e086839143..9d2cfb4a26 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/noChemistryTabulation/noChemistryTabulation.C
+++ b/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  * * * * * * * * * * * * * * * //
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/noChemistryTabulation/noChemistryTabulation.H b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/noChemistryTabulation/noChemistryTabulation.H
index a962160dbb..3b0ece2cbe 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/noChemistryTabulation/noChemistryTabulation.H
+++ b/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
-        virtual label size()
+        //- Return false as this tabulation method does not tabulate
+        virtual bool tabulates()
         {
-            NotImplemented;
-            return 0;
-        }
-
-        virtual void writePerformance()
-        {
-            NotImplemented;
+            return false;
         }
 
         //- 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()
+        {}
 };