From 1932addc2c213b5fd95b88277661b6fd4ea85669 Mon Sep 17 00:00:00 2001
From: Henry Weller <http://cfd.direct>
Date: Tue, 7 Sep 2021 21:03:55 +0100
Subject: [PATCH] TDACChemistryModel: Rationalised interface to reduction and
 tabulation methods

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

This is the first step towards merging TDACChemistryModel with standardChemistryModel.
---
 .../chemistry/TDAC/chemistryProperties.cfg    |   6 -
 .../TDAC/chemistryPropertiesFlame.cfg         |   6 -
 .../TDACChemistryModel/TDACChemistryModel.C   | 144 ++-----------
 .../TDACChemistryModel/TDACChemistryModel.H   |  41 +---
 .../TDACChemistryModel/TDACChemistryModelI.H  |  25 ---
 .../TDACChemistryModel/reduction/DAC/DAC.C    |  60 +++---
 .../TDACChemistryModel/reduction/DAC/DAC.H    |   4 +-
 .../TDACChemistryModel/reduction/DRG/DRG.C    |  44 ++--
 .../TDACChemistryModel/reduction/DRG/DRG.H    |   6 +-
 .../reduction/DRGEP/DRGEP.C                   |  68 +++---
 .../reduction/DRGEP/DRGEP.H                   |   6 +-
 .../TDACChemistryModel/reduction/EFA/EFA.C    |  62 +++---
 .../TDACChemistryModel/reduction/EFA/EFA.H    |   8 +-
 .../TDACChemistryModel/reduction/PFA/PFA.C    |  56 ++---
 .../TDACChemistryModel/reduction/PFA/PFA.H    |   8 +-
 .../chemistryReduction/chemistryReduction.C   | 112 ++++++++++
 .../chemistryReduction/chemistryReduction.H   | 156 ++++++++++++++
 .../chemistryReduction/chemistryReductionI.H  |  53 +++++
 .../chemistryReductionMethod.C                |  11 +-
 .../chemistryReductionMethod.H                |  51 ++---
 .../chemistryReductionMethodI.H               |  35 +--
 .../chemistryReductionMethodNew.C             |  89 ++++----
 .../noChemistryReduction.C                    |   6 +-
 .../noChemistryReduction.H                    |  13 +-
 .../TDACChemistryModel/tabulation/ISAT/ISAT.C | 204 +++++++++++++-----
 .../TDACChemistryModel/tabulation/ISAT/ISAT.H |  55 ++++-
 .../tabulation/ISAT/binaryTree/binaryTree.C   |  14 +-
 .../tabulation/ISAT/binaryTree/binaryTree.H   |  11 +-
 .../ISAT/chemPointISAT/chemPointISAT.C        |  32 +--
 .../ISAT/chemPointISAT/chemPointISAT.H        |  14 +-
 .../chemistryTabulationMethod.C               |   9 +-
 .../chemistryTabulationMethod.H               |  41 +---
 .../chemistryTabulationMethodNew.C            |  89 ++++----
 .../noChemistryTabulation.C                   |   6 +-
 .../noChemistryTabulation.H                   |  25 ++-
 35 files changed, 903 insertions(+), 667 deletions(-)
 create mode 100644 src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/chemistryReduction/chemistryReduction.C
 create mode 100644 src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/chemistryReduction/chemistryReduction.H
 create mode 100644 src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/chemistryReduction/chemistryReductionI.H

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()
+        {}
 };