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.

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;
+    scalar omegaf, omegar;
-    label lRef, rRef;
 
     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, omegaf, omegar);
-                R.omega(pi, Ti, c_, celli, pf, cf, lRef, pr, cr, rRef);
 
                 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;
+    scalar omegaf, omegar;
-    label lRef, rRef;
 
     forAll(rho, celli)
     {
@@ -376,16 +373,13 @@ Foam::StandardChemistryModel<ThermoType>::calculateRR
         }
 
         const Reaction<ThermoType>& R = reactions_[ri];
-        const scalar omegai = R.omega
+        const scalar omegaI = R.omega(pi, Ti, c_, celli, omegaf, omegar);
-        (
-            pi, Ti, c_, celli, pf, cf, lRef, pr, cr, rRef
-        );
 
         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;
             }
         }
 

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;
+    scalar omegaf, omegar;
-    label lRef, rRef;
 
     dcdt = Zero;
 
@@ -145,33 +144,26 @@ void Foam::TDACChemistryModel<ThermoType>::omega
         if (!reactionsDisabled_[i])
         {
             const Reaction<ThermoType>& R = this->reactions_[i];
-
+            const scalar omegaI = R.omega(p, T, c, li, omegaf, omegar);
-            scalar omegai = R.omega
-            (
-                p, T, c, li, pf, cf, lRef, pr, cr, rRef
-            );
 
             forAll(R.lhs(), s)
             {
-                label si = R.lhs()[s].index;
+                const label si =
-                if (reduced)
+                    reduced
-                {
+                  ? completeToSimplifiedIndex_[R.lhs()[s].index]
-                    si = completeToSimplifiedIndex_[si];
+                  : 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;
+                const label si =
-                if (reduced)
+                    reduced
-                {
+                  ? completeToSimplifiedIndex_[R.rhs()[s].index]
-                    si = completeToSimplifiedIndex_[si];
+                  : R.rhs()[s].index;
-                }
-
                 const scalar sr = R.rhs()[s].stoichCoeff;
-                dcdt[si] += sr*omegai;
+                dcdt[si] += sr*omegaI;
             }
         }
     }

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

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

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,

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)
         {

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

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;
+    scalar omegaf, omegar;
-    label lRef, rRef;
 
-    scalar omegaI = omega
+    const scalar omegaI = omega(p, T, c, li, omegaf, omegar);
-    (
-        p, T, c, li, pf, cf, lRef, pr, cr, rRef
-    );
 
     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& omegaf,
-    scalar& cf,
+    scalar& omegar
-    label& lRef,
-    scalar& pr,
-    scalar& cr,
-    label& rRef
 ) const
 {
-
     scalar clippedT = min(max(T, this->Tlow()), this->Thigh());
 
-    const scalar kf = this->kf(p, clippedT, c, li);
+    omegaf = this->kf(p, clippedT, c, li);
-    const scalar kr = this->kr(kf, p, clippedT, c, li);
+    omegar = this->kr(omegaf, p, clippedT, c, li);
 
-    pf = 1;
+    forAll(lhs(), i)
-    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++)
     {
-        const label si = lhs()[s].index;
+        const label si = lhs()[i].index;
-
+        const scalar el = lhs()[i].exponent;
-        if (c[si] < c[lRef])
+        omegaf *= c[si] >= small || el >= 1 ? pow(max(c[si], 0), el) : 0;
-        {
-            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);
-        }
     }
-    cf = max(c[lRef], 0);
+    forAll(rhs(), i)
-
     {
-        const scalar exp = lhs()[slRef].exponent;
+        const label si = rhs()[i].index;
-        if (exp < 1)
+        const scalar er = rhs()[i].exponent;
-        {
+        omegar *= c[si] >= small || er >= 1 ? pow(max(c[si], 0), er) : 0;
-            if (cf > small)
-            {
-                pf *= pow(cf, exp - 1);
-            }
-            else
-            {
-                pf = 0;
-            }
-        }
-        else
-        {
-            pf *= pow(cf, exp - 1);
-        }
     }
 
-    label srRef = 0;
+    return omegaf - omegar;
-    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;
 }
 
 
@@ -449,10 +363,9 @@ void Foam::Reaction<ReactionThermo>::dwdc
     const List<label>& c2s
 ) const
 {
-    scalar pf, cf, pr, cr;
+    scalar omegaf, omegar;
-    label lRef, rRef;
 
-    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)
+                kf *=
-                {
+                    c[si] >= small || el >= 1
-                    if (c[si] > SMALL)
+                  ? el*pow(max(c[si], 0), el - 1)
-                    {
+                  : 0;
-                        kf *= el*pow(c[si] + VSMALL, el - 1);
-                    }
-                    else
-                    {
-                        kf = 0;
-                    }
-                }
-                else
-                {
-                    kf *= el*pow(c[si], el - 1);
-                }
             }
             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)
+                kr *=
-                {
+                    c[si] >= small || er >= 1
-                    if (c[si] > SMALL)
+                  ? er*pow(max(c[si], 0), er - 1)
-                    {
+                  : 0;
-                        kr *= er*pow(c[si] + VSMALL, er - 1);
-                    }
-                    else
-                    {
-                        kr = 0;
-                    }
-                }
-                else
-                {
-                    kr *= er*pow(c[si], er - 1);
-                }
             }
             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;

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& omegaf,
-                scalar& cf,
+                scalar& omegar
-                label& lRef,
-                scalar& pr,
-                scalar& cr,
-                label& rRef
             ) const;
 
         // Reaction rate coefficients