From ad82628b14856addbffe2dd018c3c6d998067ddc Mon Sep 17 00:00:00 2001
From: Will Bainbridge <http://cfd.direct>
Date: Thu, 22 Apr 2021 11:36:54 +0100
Subject: [PATCH] Reaction: Simplified rate methods

The linearisation of the reaction rate relative to the concentration of
a reference specie is not required anywhere. This data has been removed
from the rate method's output and the internals of the rate method
simplified accordingly. The clipping in the rate methods has also been
simplified and made consistent across the different rate and rate
derivative methods.
---
 .../StandardChemistryModel.C                  |  20 +-
 .../TDACChemistryModel/TDACChemistryModel.C   |  34 ++--
 .../TDACChemistryModel/reduction/DAC/DAC.C    |   5 +-
 .../TDACChemistryModel/reduction/DRG/DRG.C    |   5 +-
 .../reduction/DRGEP/DRGEP.C                   |   5 +-
 .../TDACChemistryModel/reduction/EFA/EFA.C    |   7 +-
 .../TDACChemistryModel/reduction/PFA/PFA.C    |   5 +-
 .../reaction/Reactions/Reaction/Reaction.C    | 173 ++++--------------
 .../reaction/Reactions/Reaction/Reaction.H    |   8 +-
 9 files changed, 68 insertions(+), 194 deletions(-)

diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/StandardChemistryModel/StandardChemistryModel.C b/src/thermophysicalModels/chemistryModel/chemistryModel/StandardChemistryModel/StandardChemistryModel.C
index a8864b5256..9308097051 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/StandardChemistryModel/StandardChemistryModel.C
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/StandardChemistryModel/StandardChemistryModel.C
@@ -255,8 +255,7 @@ Foam::StandardChemistryModel<ThermoType>::tc() const
 
     const label nReaction = reactions_.size();
 
-    scalar pf, cf, pr, cr;
-    label lRef, rRef;
+    scalar omegaf, omegar;
 
     if (this->chemistry_)
     {
@@ -279,12 +278,11 @@ Foam::StandardChemistryModel<ThermoType>::tc() const
             forAll(reactions_, i)
             {
                 const Reaction<ThermoType>& R = reactions_[i];
-
-                R.omega(pi, Ti, c_, celli, pf, cf, lRef, pr, cr, rRef);
+                R.omega(pi, Ti, c_, celli, omegaf, omegar);
 
                 forAll(R.rhs(), s)
                 {
-                    tc[celli] += R.rhs()[s].stoichCoeff*pf*cf;
+                    tc[celli] += R.rhs()[s].stoichCoeff*omegaf;
                 }
             }
 
@@ -360,8 +358,7 @@ Foam::StandardChemistryModel<ThermoType>::calculateRR
 
     reactionEvaluationScope scope(*this);
 
-    scalar pf, cf, pr, cr;
-    label lRef, rRef;
+    scalar omegaf, omegar;
 
     forAll(rho, celli)
     {
@@ -376,16 +373,13 @@ Foam::StandardChemistryModel<ThermoType>::calculateRR
         }
 
         const Reaction<ThermoType>& R = reactions_[ri];
-        const scalar omegai = R.omega
-        (
-            pi, Ti, c_, celli, pf, cf, lRef, pr, cr, rRef
-        );
+        const scalar omegaI = R.omega(pi, Ti, c_, celli, omegaf, omegar);
 
         forAll(R.lhs(), s)
         {
             if (si == R.lhs()[s].index)
             {
-                RR[celli] -= R.lhs()[s].stoichCoeff*omegai;
+                RR[celli] -= R.lhs()[s].stoichCoeff*omegaI;
             }
         }
 
@@ -393,7 +387,7 @@ Foam::StandardChemistryModel<ThermoType>::calculateRR
         {
             if (si == R.rhs()[s].index)
             {
-                RR[celli] += R.rhs()[s].stoichCoeff*omegai;
+                RR[celli] += R.rhs()[s].stoichCoeff*omegaI;
             }
         }
 
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModel.C b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModel.C
index dd14ae8e29..4d51cf081b 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModel.C
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModel.C
@@ -135,8 +135,7 @@ void Foam::TDACChemistryModel<ThermoType>::omega
 {
     const bool reduced = mechRed_->active();
 
-    scalar pf, cf, pr, cr;
-    label lRef, rRef;
+    scalar omegaf, omegar;
 
     dcdt = Zero;
 
@@ -145,33 +144,26 @@ void Foam::TDACChemistryModel<ThermoType>::omega
         if (!reactionsDisabled_[i])
         {
             const Reaction<ThermoType>& R = this->reactions_[i];
-
-            scalar omegai = R.omega
-            (
-                p, T, c, li, pf, cf, lRef, pr, cr, rRef
-            );
+            const scalar omegaI = R.omega(p, T, c, li, omegaf, omegar);
 
             forAll(R.lhs(), s)
             {
-                label si = R.lhs()[s].index;
-                if (reduced)
-                {
-                    si = completeToSimplifiedIndex_[si];
-                }
-
+                const label si =
+                    reduced
+                  ? completeToSimplifiedIndex_[R.lhs()[s].index]
+                  : R.lhs()[s].index;
                 const scalar sl = R.lhs()[s].stoichCoeff;
-                dcdt[si] -= sl*omegai;
+                dcdt[si] -= sl*omegaI;
             }
+
             forAll(R.rhs(), s)
             {
-                label si = R.rhs()[s].index;
-                if (reduced)
-                {
-                    si = completeToSimplifiedIndex_[si];
-                }
-
+                const label si =
+                    reduced
+                  ? completeToSimplifiedIndex_[R.rhs()[s].index]
+                  : R.rhs()[s].index;
                 const scalar sr = R.rhs()[s].stoichCoeff;
-                dcdt[si] += sr*omegai;
+                dcdt[si] += sr*omegaI;
             }
         }
     }
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DAC/DAC.C b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DAC/DAC.C
index 68f179e2c3..a4dc6fc309 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DAC/DAC.C
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DAC/DAC.C
@@ -263,14 +263,13 @@ void Foam::chemistryReductionMethods::DAC<ThermoType>::reduceMechanism
     // Index of the other species involved in the rABNum
     RectangularMatrix<label> rABOtherSpec(this->nSpecie_, this->nSpecie_, -1);
 
-    scalar pf, cf, pr, cr;
-    label lRef, rRef;
     forAll(this->chemistry_.reactions(), i)
     {
         const Reaction<ThermoType>& R = this->chemistry_.reactions()[i];
 
         // for each reaction compute omegai
-        scalar omegai = R.omega(p, T, c1, li, pf, cf, lRef, pr, cr, rRef);
+        scalar omegaf, omegar;
+        const scalar omegai = R.omega(p, T, c1, li, omegaf, omegar);
 
         // Then for each pair of species composing this reaction,
         // compute the rAB matrix (separate the numerator and
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DRG/DRG.C b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DRG/DRG.C
index a807df2ca2..35011e41ef 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DRG/DRG.C
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DRG/DRG.C
@@ -86,14 +86,13 @@ void Foam::chemistryReductionMethods::DRG<ThermoType>::reduceMechanism
     // Index of the other species involved in the rABNum
     RectangularMatrix<label> rABOtherSpec(this->nSpecie_, this->nSpecie_, -1);
 
-    scalar pf, cf, pr, cr;
-    label lRef, rRef;
     forAll(this->chemistry_.reactions(), i)
     {
         const Reaction<ThermoType>& R = this->chemistry_.reactions()[i];
 
         // For each reaction compute omegai
-        scalar omegai = R.omega(p, T, c1, li, pf, cf, lRef, pr, cr, rRef);
+        scalar omegaf, omegar;
+        const scalar omegai = R.omega(p, T, c1, li, omegaf, omegar);
 
         // Then for each pair of species composing this reaction,
         // compute the rAB matrix (separate the numerator and
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DRGEP/DRGEP.C b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DRGEP/DRGEP.C
index ce765a8e48..c4c45efc78 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DRGEP/DRGEP.C
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DRGEP/DRGEP.C
@@ -131,15 +131,14 @@ void Foam::chemistryReductionMethods::DRGEP<ThermoType>::reduceMechanism
     // Index of the other species involved in the rABNum
     RectangularMatrix<label> rABOtherSpec(this->nSpecie_, this->nSpecie_, -1);
 
-    scalar pf, cf, pr, cr;
-    label lRef, rRef;
     scalarField omegaV(this->chemistry_.reactions().size());
     forAll(this->chemistry_.reactions(), i)
     {
         const Reaction<ThermoType>& R = this->chemistry_.reactions()[i];
 
         // for each reaction compute omegai
-        scalar omegai = R.omega(p, T, c1, li, pf, cf, lRef, pr, cr, rRef);
+        scalar omegaf, omegar;
+        const scalar omegai = R.omega(p, T, c1, li, omegaf, omegar);
         omegaV[i] = omegai;
 
         // then for each pair of species composing this reaction,
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/EFA/EFA.C b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/EFA/EFA.C
index a42532b61b..709364751a 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/EFA/EFA.C
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/EFA/EFA.C
@@ -127,16 +127,15 @@ void Foam::chemistryReductionMethods::EFA<ThermoType>::reduceMechanism
     // Index of the other species involved in the rABNum
     RectangularMatrix<label> rABOtherSpec(this->nSpecie_, this->nSpecie_, -1);
 
-    scalar pf, cf, pr, cr;
-    label lRef, rRef;
     forAll(this->chemistry_.reactions(), i)
     {
         const Reaction<ThermoType>& R = this->chemistry_.reactions()[i];
 
         // for each reaction compute omegai
-        R.omega(p, T, c1, li, pf, cf, lRef, pr, cr, rRef);
+        scalar omegaf, omegar;
+        R.omega(p, T, c1, li, omegaf, omegar);
 
-        scalar fr = mag(pf*cf)+mag(pr*cr);
+        scalar fr = mag(omegaf) + mag(omegar);
         scalar NCi(0.0),NHi(0.0),NOi(0.0),NNi(0.0);
         forAll(R.lhs(),s)
         {
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/PFA/PFA.C b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/PFA/PFA.C
index 7375d0edbb..5c0a4be2dc 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/PFA/PFA.C
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/PFA/PFA.C
@@ -88,14 +88,13 @@ void Foam::chemistryReductionMethods::PFA<ThermoType>::reduceMechanism
     // Index of the other species involved in the rABNum
     RectangularMatrix<label> rABOtherSpec(this->nSpecie_, this->nSpecie_, -1);
 
-    scalar pf, cf, pr, cr;
-    label lRef, rRef;
     forAll(this->chemistry_.reactions(), i)
     {
         const Reaction<ThermoType>& R = this->chemistry_.reactions()[i];
 
         // for each reaction compute omegai
-        scalar omegai = R.omega(p, T, c1, li, pf, cf, lRef, pr, cr, rRef);
+        scalar omegaf, omegar;
+        const scalar omegai = R.omega(p, T, c1, li, omegaf, omegar);
 
         // then for each pair of species composing this reaction,
         // compute the rAB matrix (separate the numerator and
diff --git a/src/thermophysicalModels/specie/reaction/Reactions/Reaction/Reaction.C b/src/thermophysicalModels/specie/reaction/Reactions/Reaction/Reaction.C
index ac6bdf0a19..bcf0b373cd 100644
--- a/src/thermophysicalModels/specie/reaction/Reactions/Reaction/Reaction.C
+++ b/src/thermophysicalModels/specie/reaction/Reactions/Reaction/Reaction.C
@@ -295,13 +295,9 @@ void Foam::Reaction<ReactionThermo>::omega
     scalarField& dcdt
 ) const
 {
-    scalar pf, cf, pr, cr;
-    label lRef, rRef;
+    scalar omegaf, omegar;
 
-    scalar omegaI = omega
-    (
-        p, T, c, li, pf, cf, lRef, pr, cr, rRef
-    );
+    const scalar omegaI = omega(p, T, c, li, omegaf, omegar);
 
     forAll(lhs(), i)
     {
@@ -325,111 +321,29 @@ Foam::scalar Foam::Reaction<ReactionThermo>::omega
     const scalar T,
     const scalarField& c,
     const label li,
-    scalar& pf,
-    scalar& cf,
-    label& lRef,
-    scalar& pr,
-    scalar& cr,
-    label& rRef
+    scalar& omegaf,
+    scalar& omegar
 ) const
 {
-
     scalar clippedT = min(max(T, this->Tlow()), this->Thigh());
 
-    const scalar kf = this->kf(p, clippedT, c, li);
-    const scalar kr = this->kr(kf, p, clippedT, c, li);
+    omegaf = this->kf(p, clippedT, c, li);
+    omegar = this->kr(omegaf, p, clippedT, c, li);
 
-    pf = 1;
-    pr = 1;
-
-    const label Nl = lhs().size();
-    const label Nr = rhs().size();
-
-    label slRef = 0;
-    lRef = lhs()[slRef].index;
-
-    pf = kf;
-    for (label s = 1; s < Nl; s++)
+    forAll(lhs(), i)
     {
-        const label si = lhs()[s].index;
-
-        if (c[si] < c[lRef])
-        {
-            const scalar exp = lhs()[slRef].exponent;
-            pf *= pow(max(c[lRef], 0), exp);
-            lRef = si;
-            slRef = s;
-        }
-        else
-        {
-            const scalar exp = lhs()[s].exponent;
-            pf *= pow(max(c[si], 0), exp);
-        }
+        const label si = lhs()[i].index;
+        const scalar el = lhs()[i].exponent;
+        omegaf *= c[si] >= small || el >= 1 ? pow(max(c[si], 0), el) : 0;
     }
-    cf = max(c[lRef], 0);
-
+    forAll(rhs(), i)
     {
-        const scalar exp = lhs()[slRef].exponent;
-        if (exp < 1)
-        {
-            if (cf > small)
-            {
-                pf *= pow(cf, exp - 1);
-            }
-            else
-            {
-                pf = 0;
-            }
-        }
-        else
-        {
-            pf *= pow(cf, exp - 1);
-        }
+        const label si = rhs()[i].index;
+        const scalar er = rhs()[i].exponent;
+        omegar *= c[si] >= small || er >= 1 ? pow(max(c[si], 0), er) : 0;
     }
 
-    label srRef = 0;
-    rRef = rhs()[srRef].index;
-
-    // Find the matrix element and element position for the rhs
-    pr = kr;
-    for (label s = 1; s < Nr; s++)
-    {
-        const label si = rhs()[s].index;
-        if (c[si] < c[rRef])
-        {
-            const scalar exp = rhs()[srRef].exponent;
-            pr *= pow(max(c[rRef], 0), exp);
-            rRef = si;
-            srRef = s;
-        }
-        else
-        {
-            const scalar exp = rhs()[s].exponent;
-            pr *= pow(max(c[si], 0), exp);
-        }
-    }
-    cr = max(c[rRef], 0);
-
-    {
-        const scalar exp = rhs()[srRef].exponent;
-        if (exp < 1)
-        {
-            if (cr > small)
-            {
-                pr *= pow(cr, exp - 1);
-            }
-            else
-            {
-                pr = 0;
-            }
-        }
-        else
-        {
-            pr *= pow(cr, exp - 1);
-        }
-    }
-
-    return pf*cf - pr*cr;
+    return omegaf - omegar;
 }
 
 
@@ -449,10 +363,9 @@ void Foam::Reaction<ReactionThermo>::dwdc
     const List<label>& c2s
 ) const
 {
-    scalar pf, cf, pr, cr;
-    label lRef, rRef;
+    scalar omegaf, omegar;
 
-    omegaI = omega(p, T, c, li, pf, cf, lRef, pr, cr, rRef);
+    omegaI = omega(p, T, c, li, omegaf, omegar);
 
     forAll(lhs(), i)
     {
@@ -480,25 +393,17 @@ void Foam::Reaction<ReactionThermo>::dwdc
             const scalar el = lhs()[i].exponent;
             if (i == j)
             {
-                if (el < 1)
-                {
-                    if (c[si] > SMALL)
-                    {
-                        kf *= el*pow(c[si] + VSMALL, el - 1);
-                    }
-                    else
-                    {
-                        kf = 0;
-                    }
-                }
-                else
-                {
-                    kf *= el*pow(c[si], el - 1);
-                }
+                kf *=
+                    c[si] >= small || el >= 1
+                  ? el*pow(max(c[si], 0), el - 1)
+                  : 0;
             }
             else
             {
-                kf *= pow(c[si], el);
+                kf *=
+                    c[si] >= small || el >= 1
+                  ? pow(max(c[si], 0), el)
+                  : 0;
             }
         }
 
@@ -526,25 +431,17 @@ void Foam::Reaction<ReactionThermo>::dwdc
             const scalar er = rhs()[i].exponent;
             if (i == j)
             {
-                if (er < 1)
-                {
-                    if (c[si] > SMALL)
-                    {
-                        kr *= er*pow(c[si] + VSMALL, er - 1);
-                    }
-                    else
-                    {
-                        kr = 0;
-                    }
-                }
-                else
-                {
-                    kr *= er*pow(c[si], er - 1);
-                }
+                kr *=
+                    c[si] >= small || er >= 1
+                  ? er*pow(max(c[si], 0), er - 1)
+                  : 0;
             }
             else
             {
-                kr *= pow(c[si], er);
+                kr *=
+                    c[si] >= small || er >= 1
+                  ? pow(max(c[si], 0), er)
+                  : 0;
             }
         }
 
@@ -620,13 +517,13 @@ void Foam::Reaction<ReactionThermo>::dwdT
     {
         const label si = lhs()[i].index;
         const scalar el = lhs()[i].exponent;
-        dkfdT *= pow(c[si], el);
+        dkfdT *= c[si] >= small || el >= 1 ? pow(max(c[si], 0), el) : 0;
     }
     forAll(rhs(), i)
     {
         const label si = rhs()[i].index;
         const scalar er = rhs()[i].exponent;
-        dkrdT *= pow(c[si], er);
+        dkrdT *= c[si] >= small || er >= 1 ? pow(max(c[si], 0), er) : 0;
     }
 
     const scalar dqidT = dkfdT - dkrdT;
diff --git a/src/thermophysicalModels/specie/reaction/Reactions/Reaction/Reaction.H b/src/thermophysicalModels/specie/reaction/Reactions/Reaction/Reaction.H
index a778592cae..734b55b18e 100644
--- a/src/thermophysicalModels/specie/reaction/Reactions/Reaction/Reaction.H
+++ b/src/thermophysicalModels/specie/reaction/Reactions/Reaction/Reaction.H
@@ -241,12 +241,8 @@ public:
                 const scalar T,
                 const scalarField& c,
                 const label li,
-                scalar& pf,
-                scalar& cf,
-                label& lRef,
-                scalar& pr,
-                scalar& cr,
-                label& rRef
+                scalar& omegaf,
+                scalar& omegar
             ) const;
 
         // Reaction rate coefficients