From 593befd959664ff0696aa36109fe76c00a998a84 Mon Sep 17 00:00:00 2001
From: Will Bainbridge <http://cfd.direct>
Date: Thu, 16 Dec 2021 11:08:19 +0000
Subject: [PATCH] chemistryReductionMethod: Reduce code duplication

---
 .../chemistryModel/chemistryModel.C           |  18 +-
 .../chemistryModel/reduction/DAC/DAC.C        |  77 +--------
 .../chemistryModel/reduction/DAC/DAC.H        |   5 +-
 .../chemistryModel/reduction/DRG/DRG.C        |  71 +-------
 .../chemistryModel/reduction/DRG/DRG.H        |   5 +-
 .../chemistryModel/reduction/DRGEP/DRGEP.C    |  74 +-------
 .../chemistryModel/reduction/DRGEP/DRGEP.H    |   5 +-
 .../chemistryModel/reduction/EFA/EFA.C        |  69 +-------
 .../chemistryModel/reduction/EFA/EFA.H        |   5 +-
 .../chemistryModel/reduction/PFA/PFA.C        |  66 +-------
 .../chemistryModel/reduction/PFA/PFA.H        |   5 +-
 .../chemistryReduction/chemistryReduction.C   | 116 -------------
 .../chemistryReduction/chemistryReduction.H   | 159 ------------------
 .../chemistryReduction/chemistryReductionI.H  |  53 ------
 .../chemistryReductionMethod.C                | 128 +++++++++++++-
 .../chemistryReductionMethod.H                |  87 ++++++++--
 .../chemistryReductionMethodI.H               |  43 +++--
 .../noChemistryReduction.C                    |   1 -
 .../noChemistryReduction.H                    |   1 -
 19 files changed, 263 insertions(+), 725 deletions(-)
 delete mode 100644 src/thermophysicalModels/chemistryModel/chemistryModel/reduction/chemistryReduction/chemistryReduction.C
 delete mode 100644 src/thermophysicalModels/chemistryModel/chemistryModel/reduction/chemistryReduction/chemistryReduction.H
 delete mode 100644 src/thermophysicalModels/chemistryModel/chemistryModel/reduction/chemistryReduction/chemistryReductionI.H

diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel.C b/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel.C
index 92e0f95329..bc25f8ce49 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel.C
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/chemistryModel.C
@@ -791,28 +791,20 @@ Foam::scalar Foam::chemistryModel<ThermoType>::solve
         {
             if (mechRedActive_)
             {
+                // Compute concentrations
                 for (label i=0; i<nSpecie_; i++)
                 {
-                    const scalar Yi = Yvf_[i][celli];
-                    c_[i] = rho0*Yi/specieThermos_[i].W();
+                    c_[i] = rho0*Y_[i]/specieThermos_[i].W();
                 }
 
                 // Reduce mechanism change the number of species (only active)
-                mechRed_.reduceMechanism
-                (
-                    p,
-                    T,
-                    c_,
-                    sc_,
-                    cTos_,
-                    sToc_,
-                    celli
-                );
+                mechRed_.reduceMechanism(p, T, c_, cTos_, sToc_, celli);
 
+                // Set the simplified mass fraction field
                 sY_.setSize(nSpecie_);
                 for (label i=0; i<nSpecie_; i++)
                 {
-                    sY_[i] = specieThermos_[sToc(i)].W()*sc_[i]/rho0;
+                    sY_[i] = Y_[sToc(i)];
                 }
             }
 
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/reduction/DAC/DAC.C b/src/thermophysicalModels/chemistryModel/chemistryModel/reduction/DAC/DAC.C
index b740b219ec..b35858ad72 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/reduction/DAC/DAC.C
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/reduction/DAC/DAC.C
@@ -34,7 +34,7 @@ Foam::chemistryReductionMethods::DAC<ThermoType>::DAC
     chemistryModel<ThermoType>& chemistry
 )
 :
-    chemistryReduction<ThermoType>(dict, chemistry),
+    chemistryReductionMethod<ThermoType>(dict, chemistry),
     searchInitSet_(),
     zprime_(0),
     nbCLarge_(3),
@@ -237,13 +237,12 @@ void Foam::chemistryReductionMethods::DAC<ThermoType>::reduceMechanism
     const scalar p,
     const scalar T,
     const scalarField& c,
-    DynamicField<scalar>& sc,
     List<label>& ctos,
     DynamicList<label>& stoc,
     const label li
 )
 {
-    chemistryReduction<ThermoType>::initReduceMechanism();
+    chemistryReductionMethod<ThermoType>::initReduceMechanism();
 
     scalarField c1(this->chemistry_.nEqns(), 0.0);
     for(label i=0; i<this->nSpecie(); i++)
@@ -494,7 +493,6 @@ 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);
-    label speciesNumber = 0;
 
     // Set all species to inactive and activate them according
     // to rAB and initial set
@@ -517,17 +515,14 @@ void Foam::chemistryReductionMethods::DAC<ThermoType>::reduceMechanism
             // When phiLarge and phiProgress >= phiTol then
             // CO, HO2 and fuel are in the SIS
             Q.push(COId_);
-            speciesNumber++;
             this->activeSpecies_[COId_] = true;
             Rvalue[COId_] = 1.0;
             Q.push(HO2Id_);
-            speciesNumber++;
             this->activeSpecies_[HO2Id_] = true;
             Rvalue[HO2Id_] = 1.0;
             forAll(fuelSpeciesID_,i)
             {
                 Q.push(fuelSpeciesID_[i]);
-                speciesNumber++;
                 this->activeSpecies_[fuelSpeciesID_[i]] = true;
                 Rvalue[fuelSpeciesID_[i]] = 1.0;
             }
@@ -538,11 +533,9 @@ void Foam::chemistryReductionMethods::DAC<ThermoType>::reduceMechanism
             // When phiLarge < phiTol and phiProgress >= phiTol then
             // CO, HO2 are in the SIS
             Q.push(COId_);
-            speciesNumber++;
             this->activeSpecies_[COId_] = true;
             Rvalue[COId_] = 1.0;
             Q.push(HO2Id_);
-            speciesNumber++;
             this->activeSpecies_[HO2Id_] = true;
             Rvalue[HO2Id_] = 1.0;
 
@@ -551,7 +544,6 @@ void Foam::chemistryReductionMethods::DAC<ThermoType>::reduceMechanism
                 forAll(fuelSpeciesID_,i)
                 {
                     Q.push(fuelSpeciesID_[i]);
-                    speciesNumber++;
                     this->activeSpecies_[fuelSpeciesID_[i]] = true;
                     Rvalue[fuelSpeciesID_[i]] = 1.0;
                 }
@@ -562,12 +554,10 @@ void Foam::chemistryReductionMethods::DAC<ThermoType>::reduceMechanism
             // When phiLarge and phiProgress< phiTol then
             // CO2, H2O are in the SIS
             Q.push(CO2Id_);
-            speciesNumber++;
             this->activeSpecies_[CO2Id_] = true;
             Rvalue[CO2Id_] = 1.0;
 
             Q.push(H2OId_);
-            speciesNumber++;
             this->activeSpecies_[H2OId_] = true;
             Rvalue[H2OId_] = 1.0;
             if (forceFuelInclusion_)
@@ -575,7 +565,6 @@ void Foam::chemistryReductionMethods::DAC<ThermoType>::reduceMechanism
                 forAll(fuelSpeciesID_,i)
                 {
                     Q.push(fuelSpeciesID_[i]);
-                    speciesNumber++;
                     this->activeSpecies_[fuelSpeciesID_[i]] = true;
                     Rvalue[fuelSpeciesID_[i]] = 1.0;
                 }
@@ -585,7 +574,6 @@ void Foam::chemistryReductionMethods::DAC<ThermoType>::reduceMechanism
         if (T>NOxThreshold_ && NOId_!=-1)
         {
             Q.push(NOId_);
-            speciesNumber++;
             this->activeSpecies_[NOId_] = true;
             Rvalue[NOId_] = 1.0;
         }
@@ -596,7 +584,6 @@ void Foam::chemistryReductionMethods::DAC<ThermoType>::reduceMechanism
         {
             label q = SIS[i];
             this->activeSpecies_[q] = true;
-            speciesNumber++;
             Q.push(q);
             Rvalue[q] = 1.0;
         }
@@ -631,7 +618,6 @@ void Foam::chemistryReductionMethods::DAC<ThermoType>::reduceMechanism
                         if (!this->activeSpecies_[otherSpec])
                         {
                             this->activeSpecies_[otherSpec] = true;
-                            speciesNumber++;
                         }
                     }
                 }
@@ -639,64 +625,7 @@ void Foam::chemistryReductionMethods::DAC<ThermoType>::reduceMechanism
         }
     }
 
-    // Put a flag on the reactions containing at least one removed species
-    forAll(this->chemistry_.reactions(), i)
-    {
-        const Reaction<ThermoType>& R = this->chemistry_.reactions()[i];
-        this->reactionsDisabled_[i] = false;
-
-        forAll(R.lhs(), s)
-        {
-            label ss = R.lhs()[s].index;
-            if (!this->activeSpecies_[ss])
-            {
-                // Flag the reaction to disable it
-                this->reactionsDisabled_[i] = true;
-                break;
-            }
-        }
-        if (!this->reactionsDisabled_[i])
-        {
-            forAll(R.rhs(), s)
-            {
-                label ss = R.rhs()[s].index;
-                if (!this->activeSpecies_[ss])
-                {
-                    // Flag the reaction to disable it
-                    this->reactionsDisabled_[i] = true;
-                    break;
-                }
-            }
-        }
-    }
-
-    this->nActiveSpecies_ = speciesNumber;
-    sc.setSize(this->nActiveSpecies_ + 2);
-    stoc.setSize(this->nActiveSpecies_);
-
-    label j = 0;
-    for (label i=0; i<this->nSpecie(); i++)
-    {
-        if (this->activeSpecies_[i])
-        {
-            stoc[j] = i;
-            sc[j] = c[i];
-            ctos[i] = j++;
-            if (!this->chemistry_.active(i))
-            {
-                this->chemistry_.setActive(i);
-            }
-        }
-        else
-        {
-            ctos[i] = -1;
-        }
-     }
-
-    sc[this->nActiveSpecies_] = T;
-    sc[this->nActiveSpecies_ + 1] = p;
-
-    chemistryReduction<ThermoType>::endReduceMechanism();
+    chemistryReductionMethod<ThermoType>::endReduceMechanism(ctos, stoc);
 }
 
 
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/reduction/DAC/DAC.H b/src/thermophysicalModels/chemistryModel/chemistryModel/reduction/DAC/DAC.H
index f42f5ecf73..2dbd166652 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/reduction/DAC/DAC.H
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/reduction/DAC/DAC.H
@@ -89,7 +89,7 @@ SourceFiles
 #ifndef DAC_H
 #define DAC_H
 
-#include "chemistryReduction.H"
+#include "chemistryReductionMethod.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -104,7 +104,7 @@ namespace chemistryReductionMethods
 template<class ThermoType>
 class DAC
 :
-    public chemistryReduction<ThermoType>
+    public chemistryReductionMethod<ThermoType>
 {
     // Private Data
 
@@ -152,7 +152,6 @@ public:
             const scalar p,
             const scalar T,
             const scalarField& c,
-            DynamicField<scalar>& sc,
             List<label>& ctos,
             DynamicList<label>& stoc,
             const label li
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/reduction/DRG/DRG.C b/src/thermophysicalModels/chemistryModel/chemistryModel/reduction/DRG/DRG.C
index 8a5be3d1e3..d6bec8b44d 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/reduction/DRG/DRG.C
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/reduction/DRG/DRG.C
@@ -34,7 +34,7 @@ Foam::chemistryReductionMethods::DRG<ThermoType>::DRG
     chemistryModel<ThermoType>& chemistry
 )
 :
-    chemistryReduction<ThermoType>(dict, chemistry),
+    chemistryReductionMethod<ThermoType>(dict, chemistry),
     searchInitSet_()
 {
     const wordHashSet initSet(this->coeffsDict_.lookup("initialSet"));
@@ -60,13 +60,12 @@ void Foam::chemistryReductionMethods::DRG<ThermoType>::reduceMechanism
     const scalar p,
     const scalar T,
     const scalarField& c,
-    DynamicField<scalar>& sc,
     List<label>& ctos,
     DynamicList<label>& stoc,
     const label li
 )
 {
-    chemistryReduction<ThermoType>::initReduceMechanism();
+    chemistryReductionMethod<ThermoType>::initReduceMechanism();
 
     scalarField c1(this->nSpecie()+2, 0.0);
 
@@ -209,8 +208,6 @@ void Foam::chemistryReductionMethods::DRG<ThermoType>::reduceMechanism
     }
     // rii = 0.0 by definition
 
-    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++)
@@ -226,7 +223,6 @@ void Foam::chemistryReductionMethods::DRG<ThermoType>::reduceMechanism
     {
         label q = searchInitSet_[i];
         this->activeSpecies_[q] = true;
-        speciesNumber++;
         Q.push(q);
     }
 
@@ -258,73 +254,12 @@ void Foam::chemistryReductionMethods::DRG<ThermoType>::reduceMechanism
                 {
                     Q.push(otherSpec);
                     this->activeSpecies_[otherSpec] = true;
-                    speciesNumber++;
                 }
             }
         }
     }
 
-    // Put a flag on the reactions containing at least one removed species
-    forAll(this->chemistry_.reactions(), i)
-    {
-        const Reaction<ThermoType>& R = this->chemistry_.reactions()[i];
-        this->reactionsDisabled_[i] = false;
-
-        forAll(R.lhs(), s)
-        {
-            label ss = R.lhs()[s].index;
-
-            // The species is inactive then the reaction is removed
-            if (!this->activeSpecies_[ss])
-            {
-                // Flag the reaction to disable it
-                this->reactionsDisabled_[i] = true;
-                break;
-            }
-        }
-
-        // If the reaction has not been disabled yet
-        if (!this->reactionsDisabled_[i])
-        {
-            forAll(R.rhs(), s)
-            {
-                label ss = R.rhs()[s].index;
-                if (!this->activeSpecies_[ss])
-                {
-                    this->reactionsDisabled_[i] = true;
-                    break;
-                }
-            }
-        }
-    }
-
-    this->nActiveSpecies_ = speciesNumber;
-    sc.setSize(this->nActiveSpecies_ + 2);
-    stoc.setSize(this->nActiveSpecies_);
-
-    label j = 0;
-    for (label i=0; i<this->nSpecie(); i++)
-    {
-        if (this->activeSpecies_[i])
-        {
-            stoc[j] = i;
-            sc[j] = c[i];
-            ctos[i] = j++;
-            if (!this->chemistry_.active(i))
-            {
-                this->chemistry_.setActive(i);
-            }
-        }
-        else
-        {
-            ctos[i] = -1;
-        }
-    }
-
-    sc[this->nActiveSpecies_] = T;
-    sc[this->nActiveSpecies_+1] = p;
-
-    chemistryReduction<ThermoType>::endReduceMechanism();
+    chemistryReductionMethod<ThermoType>::endReduceMechanism(ctos, stoc);
 }
 
 
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/reduction/DRG/DRG.H b/src/thermophysicalModels/chemistryModel/chemistryModel/reduction/DRG/DRG.H
index 5a63f82d02..77b6356009 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/reduction/DRG/DRG.H
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/reduction/DRG/DRG.H
@@ -35,7 +35,7 @@ SourceFiles
 #ifndef DRG_H
 #define DRG_H
 
-#include "chemistryReduction.H"
+#include "chemistryReductionMethod.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -51,7 +51,7 @@ namespace chemistryReductionMethods
 template<class ThermoType>
 class DRG
 :
-    public chemistryReduction<ThermoType>
+    public chemistryReductionMethod<ThermoType>
 {
     // Private Data
 
@@ -86,7 +86,6 @@ public:
             const scalar p,
             const scalar T,
             const scalarField& c,
-            DynamicField<scalar>& sc,
             List<label>& ctos,
             DynamicList<label>& stoc,
             const label li
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/reduction/DRGEP/DRGEP.C b/src/thermophysicalModels/chemistryModel/chemistryModel/reduction/DRGEP/DRGEP.C
index 8b89f04647..dbda10a20a 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/reduction/DRGEP/DRGEP.C
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/reduction/DRGEP/DRGEP.C
@@ -35,7 +35,7 @@ Foam::chemistryReductionMethods::DRGEP<ThermoType>::DRGEP
     chemistryModel<ThermoType>& chemistry
 )
 :
-    chemistryReduction<ThermoType>(dict, chemistry),
+    chemistryReductionMethod<ThermoType>(dict, chemistry),
     searchInitSet_(),
     sC_(this->nSpecie(), 0),
     sH_(this->nSpecie(), 0),
@@ -43,7 +43,7 @@ Foam::chemistryReductionMethods::DRGEP<ThermoType>::DRGEP
     sN_(this->nSpecie(), 0),
     NGroupBased_(50)
 {
-    chemistryReduction<ThermoType>::initReduceMechanism();
+    chemistryReductionMethod<ThermoType>::initReduceMechanism();
 
     const wordHashSet initSet(this->coeffsDict_.lookup("initialSet"));
     forAllConstIter(wordHashSet, initSet, iter)
@@ -106,7 +106,6 @@ void Foam::chemistryReductionMethods::DRGEP<ThermoType>::reduceMechanism
     const scalar p,
     const scalar T,
     const scalarField& c,
-    DynamicField<scalar>& sc,
     List<label>& ctos,
     DynamicList<label>& stoc,
     const label li
@@ -323,15 +322,16 @@ 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);
-    label speciesNumber = 0;
     List<bool> disabledSpecies(this->nSpecie(),false);
 
     // set all species to inactive and activate them according
     // to rAB and initial set
+    label NActiveSpecies = 0;
     for (label i=0; i<this->nSpecie(); i++)
     {
         this->activeSpecies_[i] = false;
     }
+
     // Initialise the FIFOStack for search set
     FIFOStack<label> Q;
     const labelList& SIS(this->searchInitSet_);
@@ -384,7 +384,7 @@ void Foam::chemistryReductionMethods::DRGEP<ThermoType>::reduceMechanism
         if (alphaA > this->tolerance())
         {
             this->activeSpecies_[q] = true;
-            speciesNumber++;
+            NActiveSpecies ++;
             Q.push(q);
             QStart.append(q);
             alphaQ.append(1.0);
@@ -413,7 +413,7 @@ void Foam::chemistryReductionMethods::DRGEP<ThermoType>::reduceMechanism
         Q.push(specID);
         QStart.append(specID);
         alphaQ.append(1.0);
-        speciesNumber++;
+        NActiveSpecies ++;
         Rvalue[specID] = 1.0;
         this->activeSpecies_[specID] = true;
     }
@@ -446,7 +446,7 @@ void Foam::chemistryReductionMethods::DRGEP<ThermoType>::reduceMechanism
                         if (!this->activeSpecies_[otherSpec])
                         {
                             this->activeSpecies_[otherSpec] = true;
-                            speciesNumber++;
+                            NActiveSpecies ++;
                         }
                     }
                 }
@@ -456,7 +456,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() - NActiveSpecies);
 
     // while the number of removed species is greater than NGroupBased, the rAB
     // are reevaluated according to the group based definition for each loop the
@@ -648,7 +648,6 @@ void Foam::chemistryReductionMethods::DRGEP<ThermoType>::reduceMechanism
                                 if (!this->activeSpecies_[otherSpec])
                                 {
                                     this->activeSpecies_[otherSpec] = true;
-                                    speciesNumber++;
                                     NDisabledSpecies--;
                                 }
                             }
@@ -659,62 +658,7 @@ void Foam::chemistryReductionMethods::DRGEP<ThermoType>::reduceMechanism
         }
     }
 
-    // End of group-based reduction
-
-    // Put a flag on the reactions containing at least one removed species
-    forAll(this->chemistry_.reactions(), i)
-    {
-        const Reaction<ThermoType>& R = this->chemistry_.reactions()[i];
-        this->reactionsDisabled_[i] = false;
-        forAll(R.lhs(), s)
-        {
-            label ss = R.lhs()[s].index;
-            if (!this->activeSpecies_[ss])
-            {
-                this->reactionsDisabled_[i] = true;
-                break;
-            }
-        }
-        if (!this->reactionsDisabled_[i])
-        {
-            forAll(R.rhs(), s)
-            {
-                label ss = R.rhs()[s].index;
-                if (!this->activeSpecies_[ss])
-                {
-                    this->reactionsDisabled_[i] = true;
-                    break;
-                }
-            }
-        }
-    }
-
-    this->nActiveSpecies_ = speciesNumber;
-    sc.setSize(this->nActiveSpecies_+2);
-    stoc.setSize(this->nActiveSpecies_);
-
-    label j = 0;
-    for (label i=0; i<this->nSpecie(); i++)
-    {
-        if (this->activeSpecies_[i])
-        {
-            stoc[j] = i;
-            sc[j] = c[i];
-            ctos[i] = j++;
-            if (!this->chemistry_.active(i))
-            {
-                this->chemistry_.setActive(i);
-            }
-        }
-        else
-        {
-            ctos[i] = -1;
-        }
-    }
-    sc[this->nActiveSpecies_] = T;
-    sc[this->nActiveSpecies_+1] = p;
-
-    chemistryReduction<ThermoType>::endReduceMechanism();
+    chemistryReductionMethod<ThermoType>::endReduceMechanism(ctos, stoc);
 }
 
 
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/reduction/DRGEP/DRGEP.H b/src/thermophysicalModels/chemistryModel/chemistryModel/reduction/DRGEP/DRGEP.H
index 03c838c761..1bafddd53f 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/reduction/DRGEP/DRGEP.H
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/reduction/DRGEP/DRGEP.H
@@ -103,7 +103,7 @@ SourceFiles
 #ifndef DRGEP_H
 #define DRGEP_H
 
-#include "chemistryReduction.H"
+#include "chemistryReductionMethod.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -119,7 +119,7 @@ namespace chemistryReductionMethods
 template<class ThermoType>
 class DRGEP
 :
-    public chemistryReduction<ThermoType>
+    public chemistryReductionMethod<ThermoType>
 {
     // Private Data
 
@@ -157,7 +157,6 @@ public:
             const scalar p,
             const scalar T,
             const scalarField& c,
-            DynamicField<scalar>& sc,
             List<label>& ctos,
             DynamicList<label>& stoc,
             const label li
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/reduction/EFA/EFA.C b/src/thermophysicalModels/chemistryModel/chemistryModel/reduction/EFA/EFA.C
index 1ae37fc1b1..563d7b15ff 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/reduction/EFA/EFA.C
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/reduction/EFA/EFA.C
@@ -35,7 +35,7 @@ Foam::chemistryReductionMethods::EFA<ThermoType>::EFA
     chemistryModel<ThermoType>& chemistry
 )
 :
-    chemistryReduction<ThermoType>(dict, chemistry),
+    chemistryReductionMethod<ThermoType>(dict, chemistry),
     sC_(this->nSpecie(),0),
     sH_(this->nSpecie(),0),
     sO_(this->nSpecie(),0),
@@ -96,13 +96,12 @@ void Foam::chemistryReductionMethods::EFA<ThermoType>::reduceMechanism
     const scalar p,
     const scalar T,
     const scalarField& c,
-    DynamicField<scalar>& sc,
     List<label>& ctos,
     DynamicList<label>& stoc,
     const label li
 )
 {
-    chemistryReduction<ThermoType>::initReduceMechanism();
+    chemistryReductionMethod<ThermoType>::initReduceMechanism();
 
     scalarField c1(this->chemistry_.nEqns(), 0.0);
 
@@ -245,7 +244,6 @@ 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++)
     {
         this->activeSpecies_[i] = false;
@@ -300,12 +298,10 @@ void Foam::chemistryReductionMethods::EFA<ThermoType>::reduceMechanism
                 if (!this->activeSpecies_[source[idx[startPoint+i]]])
                 {
                     this->activeSpecies_[source[idx[startPoint+i]]] = true;
-                    speciesNumber++;
                 }
                 if (!this->activeSpecies_[sink[idx[startPoint+i]]])
                 {
                     this->activeSpecies_[sink[idx[startPoint+i]]] = true;
-                    speciesNumber++;
                 }
                 if (cumFlux >= threshold)
                 {
@@ -363,12 +359,10 @@ void Foam::chemistryReductionMethods::EFA<ThermoType>::reduceMechanism
                 if (!this->activeSpecies_[source[idx[startPoint+i]]])
                 {
                     this->activeSpecies_[source[idx[startPoint+i]]] = true;
-                    speciesNumber++;
                 }
                 if (!this->activeSpecies_[sink[idx[startPoint+i]]])
                 {
                     this->activeSpecies_[sink[idx[startPoint+i]]] = true;
-                    speciesNumber++;
                 }
                 if (cumFlux >= threshold)
                 {
@@ -425,12 +419,10 @@ void Foam::chemistryReductionMethods::EFA<ThermoType>::reduceMechanism
                 if (!this->activeSpecies_[source[idx[startPoint+i]]])
                 {
                     this->activeSpecies_[source[idx[startPoint+i]]] = true;
-                    speciesNumber++;
                 }
                 if (!this->activeSpecies_[sink[idx[startPoint+i]]])
                 {
                     this->activeSpecies_[sink[idx[startPoint+i]]] = true;
-                    speciesNumber++;
                 }
                 if (cumFlux >= threshold)
                 {
@@ -487,12 +479,10 @@ void Foam::chemistryReductionMethods::EFA<ThermoType>::reduceMechanism
                 if (!this->activeSpecies_[source[idx[startPoint+i]]])
                 {
                     this->activeSpecies_[source[idx[startPoint+i]]] = true;
-                    speciesNumber++;
                 }
                 if (!this->activeSpecies_[sink[idx[startPoint+i]]])
                 {
                     this->activeSpecies_[sink[idx[startPoint+i]]] = true;
-                    speciesNumber++;
                 }
                 if (cumFlux >= threshold)
                 {
@@ -504,60 +494,7 @@ void Foam::chemistryReductionMethods::EFA<ThermoType>::reduceMechanism
         }
     }
 
-    // Put a flag on the reactions containing at least one removed species
-    forAll(this->chemistry_.reactions(), i)
-    {
-        const Reaction<ThermoType>& R = this->chemistry_.reactions()[i];
-        this->reactionsDisabled_[i] = false;
-        forAll(R.lhs(), s)
-        {
-            label ss = R.lhs()[s].index;
-            if (!this->activeSpecies_[ss])
-            {
-                this->reactionsDisabled_[i] = true;
-                break;
-            }
-        }
-        if (!this->reactionsDisabled_[i])
-        {
-            forAll(R.rhs(), s)
-            {
-                label ss = R.rhs()[s].index;
-                if (!this->activeSpecies_[ss])
-                {
-                    this->reactionsDisabled_[i] = true;
-                    break;
-                }
-            }
-        }
-    }
-
-    this->nActiveSpecies_ = speciesNumber;
-    sc.setSize(this->nActiveSpecies_+2);
-    stoc.setSize(this->nActiveSpecies_);
-    label j = 0;
-
-    for (label i=0; i<this->nSpecie(); i++)
-    {
-        if (this->activeSpecies_[i])
-        {
-            stoc[j] = i;
-            sc[j] = c[i];
-            ctos[i] = j++;
-            if (!this->chemistry_.active(i))
-            {
-                this->chemistry_.setActive(i);
-            }
-        }
-        else
-        {
-            ctos[i] = -1;
-        }
-    }
-    sc[this->nActiveSpecies_] = T;
-    sc[this->nActiveSpecies_+1] = p;
-
-    chemistryReduction<ThermoType>::endReduceMechanism();
+    chemistryReductionMethod<ThermoType>::endReduceMechanism(ctos, stoc);
 }
 
 
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/reduction/EFA/EFA.H b/src/thermophysicalModels/chemistryModel/chemistryModel/reduction/EFA/EFA.H
index c8379e80cd..9fee8f4e64 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/reduction/EFA/EFA.H
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/reduction/EFA/EFA.H
@@ -34,7 +34,7 @@ SourceFiles
 #ifndef EFA_H
 #define EFA_H
 
-#include "chemistryReduction.H"
+#include "chemistryReductionMethod.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -50,7 +50,7 @@ namespace chemistryReductionMethods
 template<class ThermoType>
 class EFA
 :
-    public chemistryReduction<ThermoType>
+    public chemistryReductionMethod<ThermoType>
 {
     // Private Data
 
@@ -85,7 +85,6 @@ public:
             const scalar p,
             const scalar T,
             const scalarField& c,
-            DynamicField<scalar>& sc,
             List<label>& ctos,
             DynamicList<label>& stoc,
             const label li
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/reduction/PFA/PFA.C b/src/thermophysicalModels/chemistryModel/chemistryModel/reduction/PFA/PFA.C
index b48d130ab2..11e60097ba 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/reduction/PFA/PFA.C
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/reduction/PFA/PFA.C
@@ -34,7 +34,7 @@ Foam::chemistryReductionMethods::PFA<ThermoType>::PFA
     chemistryModel<ThermoType>& chemistry
 )
 :
-    chemistryReduction<ThermoType>(dict, chemistry),
+    chemistryReductionMethod<ThermoType>(dict, chemistry),
     searchInitSet_()
 {
     const wordHashSet initSet(this->coeffsDict_.lookup("initialSet"));
@@ -60,13 +60,12 @@ void Foam::chemistryReductionMethods::PFA<ThermoType>::reduceMechanism
     const scalar p,
     const scalar T,
     const scalarField& c,
-    DynamicField<scalar>& sc,
     List<label>& ctos,
     DynamicList<label>& stoc,
     const label li
 )
 {
-    chemistryReduction<ThermoType>::initReduceMechanism();
+    chemistryReductionMethod<ThermoType>::initReduceMechanism();
 
     scalarField c1(this->chemistry_.nEqns(), 0.0);
 
@@ -290,8 +289,6 @@ void Foam::chemistryReductionMethods::PFA<ThermoType>::reduceMechanism
     }
 
     // Using the rAB matrix (numerator and denominator separated)
-    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++)
@@ -307,7 +304,6 @@ void Foam::chemistryReductionMethods::PFA<ThermoType>::reduceMechanism
     {
         label q = SIS[i];
         this->activeSpecies_[q] = true;
-        speciesNumber++;
         Q.push(q);
     }
 
@@ -338,7 +334,6 @@ void Foam::chemistryReductionMethods::PFA<ThermoType>::reduceMechanism
                 {
                     Q.push(otherSpec);
                     this->activeSpecies_[otherSpec] = true;
-                    speciesNumber++;
                 }
 
             }
@@ -356,67 +351,12 @@ void Foam::chemistryReductionMethods::PFA<ThermoType>::reduceMechanism
                 {
                     Q.push(otherSpec);
                     this->activeSpecies_[otherSpec] = true;
-                    speciesNumber++;
                 }
             }
         }
     }
 
-    // Put a flag on the reactions containing at least one removed species
-    forAll(this->chemistry_.reactions(), i)
-    {
-        const Reaction<ThermoType>& R = this->chemistry_.reactions()[i];
-        this->reactionsDisabled_[i] = false;
-
-        forAll(R.lhs(), s)
-        {
-            label ss = R.lhs()[s].index;
-            if (!this->activeSpecies_[ss])
-            {
-                this->reactionsDisabled_[i] = true;
-                break;
-            }
-        }
-        if (!this->reactionsDisabled_[i])
-        {
-            forAll(R.rhs(), s)
-            {
-                label ss = R.rhs()[s].index;
-                if (!this->activeSpecies_[ss])
-                {
-                    this->reactionsDisabled_[i] = true;
-                    break;
-                }
-            }
-        }
-    }
-
-    this->nActiveSpecies_ = speciesNumber;
-    sc.setSize(this->nActiveSpecies_+2);
-    stoc.setSize(this->nActiveSpecies_);
-
-    label j = 0;
-    for (label i=0; i<this->nSpecie(); i++)
-    {
-        if (this->activeSpecies_[i])
-        {
-            stoc[j] = i;
-            sc[j] = c[i];
-            ctos[i] = j++;
-            if (!this->chemistry_.active(i))
-            {
-                this->chemistry_.setActive(i);
-            }
-        }
-        else
-        {
-            ctos[i] = -1;
-        }
-    }
-    sc[this->nActiveSpecies_] = T;
-    sc[this->nActiveSpecies_+1] = p;
-
-    chemistryReduction<ThermoType>::endReduceMechanism();
+    chemistryReductionMethod<ThermoType>::endReduceMechanism(ctos, stoc);
 }
 
 
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/reduction/PFA/PFA.H b/src/thermophysicalModels/chemistryModel/chemistryModel/reduction/PFA/PFA.H
index 56ebbb0229..4b36ccbfa8 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/reduction/PFA/PFA.H
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/reduction/PFA/PFA.H
@@ -35,7 +35,7 @@ SourceFiles
 #ifndef PFA_H
 #define PFA_H
 
-#include "chemistryReduction.H"
+#include "chemistryReductionMethod.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -51,7 +51,7 @@ namespace chemistryReductionMethods
 template<class ThermoType>
 class PFA
 :
-    public chemistryReduction<ThermoType>
+    public chemistryReductionMethod<ThermoType>
 {
     // Private Data
 
@@ -87,7 +87,6 @@ public:
             const scalar p,
             const scalar T,
             const scalarField& c,
-            DynamicField<scalar>& sc,
             List<label>& ctos,
             DynamicList<label>& stoc,
             const label li
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/reduction/chemistryReduction/chemistryReduction.C b/src/thermophysicalModels/chemistryModel/chemistryModel/reduction/chemistryReduction/chemistryReduction.C
deleted file mode 100644
index 7aec6aabd7..0000000000
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/reduction/chemistryReduction/chemistryReduction.C
+++ /dev/null
@@ -1,116 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  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 "chemistryModel.H"
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-template<class ThermoType>
-Foam::chemistryReduction<ThermoType>::chemistryReduction
-(
-    const Foam::IOdictionary& dict,
-    Foam::chemistryModel<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()
-{
-    // Change temporary Ns in chemistryModel
-    // to make the function nEqns working
-    this->chemistry_.setNSpecie(this->nActiveSpecies_);
-
-    if (log_)
-    {
-        sumnActiveSpecies_ += this->nActiveSpecies_;
-        sumn_++;
-        reduceMechCpuTime_ += clockTime_.timeIncrement();
-    }
-}
-
-
-template<class ThermoType>
-void Foam::chemistryReduction<ThermoType>::update()
-{
-    if (log_)
-    {
-        cpuReduceFile_()
-            << this->chemistry_.time().userTimeValue()
-            << "    " << reduceMechCpuTime_ << endl;
-
-        if (sumn_)
-        {
-            // Write average number of species
-            nActiveSpeciesFile_()
-                << this->chemistry_.time().userTimeValue()
-                << "    " << sumnActiveSpecies_/sumn_ << endl;
-        }
-
-        sumnActiveSpecies_ = 0;
-        sumn_ = 0;
-        reduceMechCpuTime_ = 0;
-    }
-}
-
-
-// ************************************************************************* //
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/reduction/chemistryReduction/chemistryReduction.H b/src/thermophysicalModels/chemistryModel/chemistryModel/reduction/chemistryReduction/chemistryReduction.H
deleted file mode 100644
index cebf1682a2..0000000000
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/reduction/chemistryReduction/chemistryReduction.H
+++ /dev/null
@@ -1,159 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  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_;
-
-        chemistryModel<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,
-            chemistryModel<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,
-            DynamicField<scalar>& sc,
-            List<label>& ctos,
-            DynamicList<label>& stoc,
-            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/reduction/chemistryReduction/chemistryReductionI.H b/src/thermophysicalModels/chemistryModel/chemistryModel/reduction/chemistryReduction/chemistryReductionI.H
deleted file mode 100644
index f06c9eca12..0000000000
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/reduction/chemistryReduction/chemistryReductionI.H
+++ /dev/null
@@ -1,53 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  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/reduction/chemistryReductionMethod/chemistryReductionMethod.C b/src/thermophysicalModels/chemistryModel/chemistryModel/reduction/chemistryReductionMethod/chemistryReductionMethod.C
index 055f1542f9..8b248afb05 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/reduction/chemistryReductionMethod/chemistryReductionMethod.C
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/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-2021 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -35,10 +35,25 @@ Foam::chemistryReductionMethod<ThermoType>::chemistryReductionMethod
     Foam::chemistryModel<ThermoType>& chemistry
 )
 :
+    coeffsDict_(dict.subDict("reduction")),
+    chemistry_(chemistry),
     nSpecie_(chemistry.nSpecie()),
     nActiveSpecies_(chemistry.nSpecie()),
-    reactionsDisabled_(chemistry.reactions().size(), false)
-{}
+    reactionsDisabled_(chemistry.nReaction(), false),
+    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  * * * * * * * * * * * * * * * //
@@ -48,4 +63,111 @@ Foam::chemistryReductionMethod<ThermoType>::~chemistryReductionMethod()
 {}
 
 
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class ThermoType>
+void Foam::chemistryReductionMethod<ThermoType>::initReduceMechanism()
+{
+    if (log_)
+    {
+        clockTime_.timeIncrement();
+    }
+}
+
+
+template<class ThermoType>
+void Foam::chemistryReductionMethod<ThermoType>::endReduceMechanism
+(
+    List<label>& ctos,
+    DynamicList<label>& stoc
+)
+{
+    // Disable reactions containing removed species
+    forAll(chemistry_.reactions(), i)
+    {
+        const Reaction<ThermoType>& R = chemistry_.reactions()[i];
+        reactionsDisabled_[i] = false;
+
+        forAll(R.lhs(), s)
+        {
+            label ss = R.lhs()[s].index;
+            if (!activeSpecies_[ss])
+            {
+                reactionsDisabled_[i] = true;
+                break;
+            }
+        }
+
+        if (!reactionsDisabled_[i])
+        {
+            forAll(R.rhs(), s)
+            {
+                label ss = R.rhs()[s].index;
+                if (!activeSpecies_[ss])
+                {
+                    reactionsDisabled_[i] = true;
+                    break;
+                }
+            }
+        }
+    }
+
+    // Set the total number of active species
+    nActiveSpecies_ = count(activeSpecies_, true);
+
+    // Set the indexing arrays
+    stoc.setSize(nActiveSpecies_);
+    for (label i=0, j=0; i<nSpecie(); i++)
+    {
+        if (activeSpecies_[i])
+        {
+            stoc[j] = i;
+            ctos[i] = j++;
+            if (!chemistry_.active(i))
+            {
+                chemistry_.setActive(i);
+            }
+        }
+        else
+        {
+            ctos[i] = -1;
+        }
+    }
+
+    // Change the number of species in the chemistry model
+    chemistry_.setNSpecie(nActiveSpecies_);
+
+    if (log_)
+    {
+        sumnActiveSpecies_ += nActiveSpecies_;
+        sumn_++;
+        reduceMechCpuTime_ += clockTime_.timeIncrement();
+    }
+}
+
+
+template<class ThermoType>
+void Foam::chemistryReductionMethod<ThermoType>::update()
+{
+    if (log_)
+    {
+        cpuReduceFile_()
+            << chemistry_.time().userTimeValue()
+            << "    " << reduceMechCpuTime_ << endl;
+
+        if (sumn_)
+        {
+            // Write average number of species
+            nActiveSpeciesFile_()
+                << chemistry_.time().userTimeValue()
+                << "    " << sumnActiveSpecies_/sumn_ << endl;
+        }
+
+        sumnActiveSpecies_ = 0;
+        sumn_ = 0;
+        reduceMechCpuTime_ = 0;
+    }
+}
+
+
 // ************************************************************************* //
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/reduction/chemistryReductionMethod/chemistryReductionMethod.H b/src/thermophysicalModels/chemistryModel/chemistryModel/reduction/chemistryReductionMethod/chemistryReductionMethod.H
index cc7f8f2ad5..aba16b5f28 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/reduction/chemistryReductionMethod/chemistryReductionMethod.H
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/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-2021 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -29,6 +29,9 @@ Description
 
 SourceFiles
     chemistryReductionMethod.C
+    chemistryReductionMethods.C
+    chemistryReductionMethodNew.C
+    chemistryReductionMethodI.H
 
 \*---------------------------------------------------------------------------*/
 
@@ -37,6 +40,9 @@ SourceFiles
 
 #include "IOdictionary.H"
 #include "DynamicField.H"
+#include "Switch.H"
+#include "clockTime.H"
+#include "OFstream.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -47,27 +53,72 @@ template<class ThermoType>
 class chemistryModel;
 
 /*---------------------------------------------------------------------------*\
-                       Class chemistrySolver Declaration
+                   Class chemistryReductionMethod Declaration
 \*---------------------------------------------------------------------------*/
 
 template<class ThermoType>
 class chemistryReductionMethod
 {
-    // Private Data
-
-        //- Cache the number of species
-        const label nSpecie_;
-
-
 protected:
 
     // Protected Data
 
+        //- Dictionary that store the algorithm data
+        const dictionary coeffsDict_;
+
+        //- Reference to the chemistry model
+        chemistryModel<ThermoType>& chemistry_;
+
+        //- Total number of species
+        const label nSpecie_;
+
         //- Number of active species
         label nActiveSpecies_;
 
+        //- List of disabled reactions (disabled = true)
         Field<bool> reactionsDisabled_;
 
+        //- List of active species (active = true)
+        List<bool> activeSpecies_;
+
+
+    //- Protected Member Functions
+
+        //- Initialise reduction of the mechanism
+        void initReduceMechanism();
+
+        //- End reduction of the mechanism
+        void endReduceMechanism(List<label>& ctos, DynamicList<label>& stoc);
+
+
+private:
+
+    // Private Data
+
+        //- Switch to select performance logging
+        Switch log_;
+
+        //- Tolerance for the mechanism reduction algorithm
+        scalar tolerance_;
+
+        //- Clock time for logging
+        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:
 
@@ -114,8 +165,11 @@ public:
 
     // Member Functions
 
-        //- Is mechanism reduction active?
-        virtual bool active() const = 0;
+        //- Return mechanism reduction active
+        virtual bool active() const
+        {
+            return true;
+        }
 
         //- Return the number of species
         inline label nSpecie();
@@ -126,10 +180,11 @@ public:
         //- Return the tolerance
         inline scalar tolerance() const;
 
-        inline bool reactionDisabled(const label i) const
-        {
-            return reactionsDisabled_[i];
-        }
+        //- Return the active species
+        inline const List<bool>& activeSpecies() const;
+
+        //- Return whether or not a reaction is disabled
+        inline bool reactionDisabled(const label i) const;
 
         //- Reduce the mechanism
         virtual void reduceMechanism
@@ -137,13 +192,13 @@ public:
             const scalar p,
             const scalar T,
             const scalarField& c,
-            DynamicField<scalar>& sc,
             List<label>& ctos,
             DynamicList<label>& stoc,
             const label li
         ) = 0;
 
-        virtual void update() = 0;
+        //- ...
+        virtual void update();
 };
 
 
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/reduction/chemistryReductionMethod/chemistryReductionMethodI.H b/src/thermophysicalModels/chemistryModel/chemistryModel/reduction/chemistryReductionMethod/chemistryReductionMethodI.H
index 3102238a78..b0fb13e63e 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/reduction/chemistryReductionMethod/chemistryReductionMethodI.H
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/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-2021 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2021 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -21,19 +21,19 @@ License
     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::chemistryReductionMethod
-
-Description
-    An abstract class for reducing chemical mechanism
-
-SourceFiles
-    chemistryReductionMethod.C
-
 \*---------------------------------------------------------------------------*/
 
+#include "chemistryReductionMethod.H"
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+template<class ThermoType>
+inline Foam::label Foam::chemistryReductionMethod<ThermoType>::nSpecie()
+{
+    return nSpecie_;
+}
+
+
 template<class ThermoType>
 inline Foam::label
 Foam::chemistryReductionMethod<ThermoType>::nActiveSpecies() const
@@ -43,9 +43,28 @@ Foam::chemistryReductionMethod<ThermoType>::nActiveSpecies() const
 
 
 template<class ThermoType>
-inline Foam::label Foam::chemistryReductionMethod<ThermoType>::nSpecie()
+inline const Foam::List<bool>&
+Foam::chemistryReductionMethod<ThermoType>::activeSpecies() const
 {
-    return nSpecie_;
+    return activeSpecies_;
+}
+
+
+template<class ThermoType>
+inline Foam::scalar
+Foam::chemistryReductionMethod<ThermoType>::tolerance() const
+{
+    return tolerance_;
+}
+
+
+template<class ThermoType>
+inline bool Foam::chemistryReductionMethod<ThermoType>::reactionDisabled
+(
+    const label i
+) const
+{
+    return reactionsDisabled_[i];
 }
 
 
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/reduction/noChemistryReduction/noChemistryReduction.C b/src/thermophysicalModels/chemistryModel/chemistryModel/reduction/noChemistryReduction/noChemistryReduction.C
index 98d8852d67..1d664e2db3 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/reduction/noChemistryReduction/noChemistryReduction.C
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/reduction/noChemistryReduction/noChemistryReduction.C
@@ -53,7 +53,6 @@ void Foam::chemistryReductionMethods::none<ThermoType>::reduceMechanism
     const scalar p,
     const scalar T,
     const scalarField& c,
-    DynamicField<scalar>& sc,
     List<label>& ctos,
     DynamicList<label>& stoc,
     const label li
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/reduction/noChemistryReduction/noChemistryReduction.H b/src/thermophysicalModels/chemistryModel/chemistryModel/reduction/noChemistryReduction/noChemistryReduction.H
index 05a39ce8d0..9c0bc77d5c 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/reduction/noChemistryReduction/noChemistryReduction.H
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/reduction/noChemistryReduction/noChemistryReduction.H
@@ -89,7 +89,6 @@ public:
             const scalar p,
             const scalar T,
             const scalarField& c,
-            DynamicField<scalar>& sc,
             List<label>& ctos,
             DynamicList<label>& stoc,
             const label li