thermophysicalModels: Implementation of the full algebraic Jacobian

including third-body and pressure dependent derivatives, and derivative of the
temperature term.  The complete Jacobian is more robust than the incomplete and
partially approximate form used previously and improves the efficiency of the
stiff ODE solvers which rely on the Jacobian.

Reaction rate evaluation moved from the chemistryModel to specie library to
simplfy support for alternative reaction rate expressions and associated
Jacobian terms.

Temperature clipping included in the Reaction class. This is inactive by default
but for most cases it is advised to provide temperature limits (high and
low). These are provided in the foamChemistryFile with the keywords Thigh and
Tlow. When using chemkinToFoam these values are set to the limits of the Janaf
thermodynamic data.  With the new Jacobian this temperature clipping has proved
very beneficial for stability and for some cases essential.

Improvement of the TDAC MRU list better integrated in add and grow functions.

To get the most out of this significant development it is important to re-tune
the ODE integration tolerances, in particular the absTol in the odeCoeffs
sub-dictionary of the chemistryProperties dictionary:

odeCoeffs
{
    solver          seulex;
    absTol          1e-12;
    relTol          0.01;
}

Typically absTol can now be set to 1e-8 and relTol to 0.1 except for ignition
time problems, and with theses settings the integration is still robust but for
many cases a lot faster than previously.

Code development and integration undertaken by
Francesco Contino
Henry G. Weller, CFD Direct

applications/utilities/thermophysical/chemkinToFoam/chemkinToFoam.C

@@ -67,15 +67,13 @@ int main(int argc, char *argv[])
 
     chemkinReader cr(species, args[1], args[3], args[2], newFormat);
 
-
     OFstream reactionsFile(args[4]);
-    reactionsFile
+    reactionsFile.writeKeyword("elements")
-        << "elements" << cr.elementNames() << token::END_STATEMENT << nl << nl;
+        << cr.elementNames() << token::END_STATEMENT << nl << nl;
-    reactionsFile
+    reactionsFile.writeKeyword("species")
-        << "species" << cr.species() << token::END_STATEMENT << nl << nl;
+        << cr.species() << token::END_STATEMENT << nl << nl;
     cr.reactions().write(reactionsFile);
 
-
     // Temporary hack to splice the specie composition data into the thermo file
     // pending complete integration into the thermodynamics structure
 
@@ -103,6 +101,14 @@ int main(int argc, char *argv[])
 
     thermoDict.write(OFstream(args[5])(), false);
 
+    reactionsFile << nl;
+
+    reactionsFile.writeKeyword("Tlow")
+        << Reaction<gasHThermoPhysics>::TlowDefault
+        << token::END_STATEMENT << nl;
+    reactionsFile.writeKeyword("Thigh")
+        << Reaction<gasHThermoPhysics>::ThighDefault
+        << token::END_STATEMENT << nl << nl;
 
     Info<< "End\n" << endl;
 

src/thermophysicalModels/chemistryModel/chemistryModel/StandardChemistryModel/StandardChemistryModel.C

@@ -56,7 +56,7 @@ Foam::StandardChemistryModel<ReactionThermo, ThermoType>::StandardChemistryModel
         BasicChemistryModel<ReactionThermo>::template lookupOrDefault<scalar>
         (
             "Treact",
-            0.0
+            0
         )
     ),
     RR_(nSpecie_),
@@ -80,7 +80,7 @@ Foam::StandardChemistryModel<ReactionThermo, ThermoType>::StandardChemistryModel
                     IOobject::NO_WRITE
                 ),
                 thermo.p().mesh(),
-                dimensionedScalar("zero", dimMass/dimVolume/dimTime, 0.0)
+                dimensionedScalar("zero", dimMass/dimVolume/dimTime, 0)
             )
         );
     }
@@ -109,8 +109,6 @@ void Foam::StandardChemistryModel<ReactionThermo, ThermoType>::omega
     scalarField& dcdt
 ) const
 {
-    scalar pf, cf, pr, cr;
-    label lRef, rRef;
 
     dcdt = Zero;
 
@@ -118,24 +116,7 @@ void Foam::StandardChemistryModel<ReactionThermo, ThermoType>::omega
     {
         const Reaction<ThermoType>& R = reactions_[i];
 
-        const scalar omegai = omega
+        R.omega(p, T, c, dcdt);
-        (
-            R, c, T, p, pf, cf, lRef, pr, cr, rRef
-        );
-
-        forAll(R.lhs(), s)
-        {
-            const label si = R.lhs()[s].index;
-            const scalar sl = R.lhs()[s].stoichCoeff;
-            dcdt[si] -= sl*omegai;
-        }
-
-        forAll(R.rhs(), s)
-        {
-            const label si = R.rhs()[s].index;
-            const scalar sr = R.rhs()[s].stoichCoeff;
-            dcdt[si] += sr*omegai;
-        }
     }
 }
 
@@ -156,123 +137,11 @@ Foam::scalar Foam::StandardChemistryModel<ReactionThermo, ThermoType>::omegaI
 ) const
 {
     const Reaction<ThermoType>& R = reactions_[index];
-    scalar w = omega(R, c, T, p, pf, cf, lRef, pr, cr, rRef);
+    scalar w = R.omega(p, T, c, pf, cf, lRef, pr, cr, rRef);
     return(w);
 }
 
 
-template<class ReactionThermo, class ThermoType>
-Foam::scalar Foam::StandardChemistryModel<ReactionThermo, ThermoType>::omega
-(
-    const Reaction<ThermoType>& R,
-    const scalarField& c,
-    const scalar T,
-    const scalar p,
-    scalar& pf,
-    scalar& cf,
-    label& lRef,
-    scalar& pr,
-    scalar& cr,
-    label& rRef
-) const
-{
-    const scalar kf = R.kf(p, T, c);
-    const scalar kr = R.kr(kf, p, T, c);
-
-    pf = 1.0;
-    pr = 1.0;
-
-    const label Nl = R.lhs().size();
-    const label Nr = R.rhs().size();
-
-    label slRef = 0;
-    lRef = R.lhs()[slRef].index;
-
-    pf = kf;
-    for (label s = 1; s < Nl; s++)
-    {
-        const label si = R.lhs()[s].index;
-
-        if (c[si] < c[lRef])
-        {
-            const scalar exp = R.lhs()[slRef].exponent;
-            pf *= pow(max(c[lRef], 0.0), exp);
-            lRef = si;
-            slRef = s;
-        }
-        else
-        {
-            const scalar exp = R.lhs()[s].exponent;
-            pf *= pow(max(c[si], 0.0), exp);
-        }
-    }
-    cf = max(c[lRef], 0.0);
-
-    {
-        const scalar exp = R.lhs()[slRef].exponent;
-        if (exp < 1.0)
-        {
-            if (cf > small)
-            {
-                pf *= pow(cf, exp - 1.0);
-            }
-            else
-            {
-                pf = 0.0;
-            }
-        }
-        else
-        {
-            pf *= pow(cf, exp - 1.0);
-        }
-    }
-
-    label srRef = 0;
-    rRef = R.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 = R.rhs()[s].index;
-        if (c[si] < c[rRef])
-        {
-            const scalar exp = R.rhs()[srRef].exponent;
-            pr *= pow(max(c[rRef], 0.0), exp);
-            rRef = si;
-            srRef = s;
-        }
-        else
-        {
-            const scalar exp = R.rhs()[s].exponent;
-            pr *= pow(max(c[si], 0.0), exp);
-        }
-    }
-    cr = max(c[rRef], 0.0);
-
-    {
-        const scalar exp = R.rhs()[srRef].exponent;
-        if (exp < 1.0)
-        {
-            if (cr > small)
-            {
-                pr *= pow(cr, exp - 1.0);
-            }
-            else
-            {
-                pr = 0.0;
-            }
-        }
-        else
-        {
-            pr *= pow(cr, exp - 1.0);
-        }
-    }
-
-    return pf*cf - pr*cr;
-}
-
-
 template<class ReactionThermo, class ThermoType>
 void Foam::StandardChemistryModel<ReactionThermo, ThermoType>::derivatives
 (
@@ -286,29 +155,29 @@ void Foam::StandardChemistryModel<ReactionThermo, ThermoType>::derivatives
 
     forAll(c_, i)
     {
-        c_[i] = max(c[i], 0.0);
+        c_[i] = max(c[i], 0);
     }
 
     omega(c_, T, p, dcdt);
 
     // Constant pressure
     // dT/dt = ...
-    scalar rho = 0.0;
+    scalar rho = 0;
-    scalar cSum = 0.0;
+    scalar cSum = 0;
     for (label i = 0; i < nSpecie_; i++)
     {
         const scalar W = specieThermo_[i].W();
         cSum += c_[i];
         rho += W*c_[i];
     }
-    scalar cp = 0.0;
+    scalar cp = 0;
     for (label i=0; i<nSpecie_; i++)
     {
         cp += c_[i]*specieThermo_[i].cp(p, T);
     }
     cp /= rho;
 
-    scalar dT = 0.0;
+    scalar dT = 0;
     for (label i = 0; i < nSpecie_; i++)
     {
         const scalar hi = specieThermo_[i].ha(p, T);
@@ -319,7 +188,7 @@ void Foam::StandardChemistryModel<ReactionThermo, ThermoType>::derivatives
     dcdt[nSpecie_] = -dT;
 
     // dp/dt = ...
-    dcdt[nSpecie_ + 1] = 0.0;
+    dcdt[nSpecie_ + 1] = 0;
 }
 
 
@@ -329,7 +198,7 @@ void Foam::StandardChemistryModel<ReactionThermo, ThermoType>::jacobian
     const scalar t,
     const scalarField& c,
     scalarField& dcdt,
-    scalarSquareMatrix& dfdc
+    scalarSquareMatrix& J
 ) const
 {
     const scalar T = c[nSpecie_];
@@ -337,131 +206,67 @@ void Foam::StandardChemistryModel<ReactionThermo, ThermoType>::jacobian
 
     forAll(c_, i)
     {
-        c_[i] = max(c[i], 0.0);
+        c_[i] = max(c[i], 0);
     }
 
-    dfdc = Zero;
+    J = Zero;
-
+    dcdt = Zero;
-    // Length of the first argument must be nSpecie_
-    omega(c_, T, p, dcdt);
 
+    // To compute the species derivatives of the temperature term,
+    // the enthalpies of the individual species is needed
+    scalarField hi(nSpecie_);
+    scalarField cpi(nSpecie_);
+    for (label i = 0; i < nSpecie_; i++)
+    {
+        hi[i] = specieThermo_[i].ha(p, T);
+        cpi[i] = specieThermo_[i].cp(p, T);
+    }
+    scalar omegaI = 0;
+    List<label> dummy;
     forAll(reactions_, ri)
     {
         const Reaction<ThermoType>& R = reactions_[ri];
-
+        scalar kfwd, kbwd;
-        const scalar kf0 = R.kf(p, T, c_);
+        R.dwdc(p, T, c_, J, dcdt, omegaI, kfwd, kbwd, false, dummy);
-        const scalar kr0 = R.kr(kf0, p, T, c_);
+        R.dwdT(p, T, c_, omegaI, kfwd, kbwd, J, false, dummy, nSpecie_);
-
-        forAll(R.lhs(), j)
-        {
-            const label sj = R.lhs()[j].index;
-            scalar kf = kf0;
-            forAll(R.lhs(), i)
-            {
-                const label si = R.lhs()[i].index;
-                const scalar el = R.lhs()[i].exponent;
-                if (i == j)
-                {
-                    if (el < 1.0)
-                    {
-                        if (c_[si] > small)
-                        {
-                            kf *= el*pow(c_[si], el - 1.0);
-                        }
-                        else
-                        {
-                            kf = 0.0;
-                        }
-                    }
-                    else
-                    {
-                        kf *= el*pow(c_[si], el - 1.0);
-                    }
-                }
-                else
-                {
-                    kf *= pow(c_[si], el);
-                }
     }
 
-            forAll(R.lhs(), i)
+    // The species derivatives of the temperature term are partially computed
-            {
+    // while computing dwdc, they are completed hereunder:
-                const label si = R.lhs()[i].index;
+    scalar cpMean = 0;
-                const scalar sl = R.lhs()[i].stoichCoeff;
+    scalar dcpdTMean = 0;
-                dfdc(si, sj) -= sl*kf;
-            }
-            forAll(R.rhs(), i)
-            {
-                const label si = R.rhs()[i].index;
-                const scalar sr = R.rhs()[i].stoichCoeff;
-                dfdc(si, sj) += sr*kf;
-            }
-        }
-
-        forAll(R.rhs(), j)
-        {
-            const label sj = R.rhs()[j].index;
-            scalar kr = kr0;
-            forAll(R.rhs(), i)
-            {
-                const label si = R.rhs()[i].index;
-                const scalar er = R.rhs()[i].exponent;
-                if (i == j)
-                {
-                    if (er < 1.0)
-                    {
-                        if (c_[si] > small)
-                        {
-                            kr *= er*pow(c_[si], er - 1.0);
-                        }
-                        else
-                        {
-                            kr = 0.0;
-                        }
-                    }
-                    else
-                    {
-                        kr *= er*pow(c_[si], er - 1.0);
-                    }
-                }
-                else
-                {
-                    kr *= pow(c_[si], er);
-                }
-            }
-
-            forAll(R.lhs(), i)
-            {
-                const label si = R.lhs()[i].index;
-                const scalar sl = R.lhs()[i].stoichCoeff;
-                dfdc(si, sj) += sl*kr;
-            }
-            forAll(R.rhs(), i)
-            {
-                const label si = R.rhs()[i].index;
-                const scalar sr = R.rhs()[i].stoichCoeff;
-                dfdc(si, sj) -= sr*kr;
-            }
-        }
-    }
-
-    // Calculate the dcdT elements numerically
-    const scalar delta = 1.0e-3;
-
-    omega(c_, T + delta, p, dcdt_);
     for (label i=0; i<nSpecie_; i++)
     {
-        dfdc(i, nSpecie_) = dcdt_[i];
+        cpMean += c_[i]*cpi[i]; // J/(m3.K)
+        dcpdTMean += c_[i]*specieThermo_[i].dcpdT(p, T);
     }
-
+    scalar dTdt = 0.0;
-    omega(c_, T - delta, p, dcdt_);
     for (label i=0; i<nSpecie_; i++)
     {
-        dfdc(i, nSpecie_) = 0.5*(dfdc(i, nSpecie_) - dcdt_[i])/delta;
+        dTdt += hi[i]*dcdt[i]; // J/(m3.s)
+    }
+    dTdt /= -cpMean; // K/s
+
+    for (label i = 0; i < nSpecie_; i++)
+    {
+        J(nSpecie_, i) = 0;
+        for (label j = 0; j < nSpecie_; j++)
+        {
+            J(nSpecie_, i) += hi[j]*J(j, i);
+        }
+        J(nSpecie_, i) += cpi[i]*dTdt; // J/(mol.s)
+        J(nSpecie_, i) /= -cpMean;    // K/s / (mol/m3)
     }
 
-    dfdc(nSpecie_, nSpecie_) = 0;
+    // ddT of dTdt
-    dfdc(nSpecie_ + 1, nSpecie_) = 0;
+    J(nSpecie_, nSpecie_) = 0;
+    for (label i = 0; i < nSpecie_; i++)
+    {
+        J(nSpecie_, nSpecie_) += cpi[i]*dcdt[i] + hi[i]*J(i, nSpecie_);
+    }
+    J(nSpecie_, nSpecie_) += dTdt*dcpdTMean;
+    J(nSpecie_, nSpecie_) /= -cpMean;
+    J(nSpecie_, nSpecie_) += dTdt/T;
 }
 
 
@@ -509,7 +314,7 @@ Foam::StandardChemistryModel<ReactionThermo, ThermoType>::tc() const
             const scalar Ti = T[celli];
             const scalar pi = p[celli];
 
-            scalar cSum = 0.0;
+            scalar cSum = 0;
 
             for (label i=0; i<nSpecie_; i++)
             {
@@ -521,7 +326,7 @@ Foam::StandardChemistryModel<ReactionThermo, ThermoType>::tc() const
             {
                 const Reaction<ThermoType>& R = reactions_[i];
 
-                omega(R, c_, Ti, pi, pf, cf, lRef, pr, cr, rRef);
+                R.omega(pi, Ti, c_, pf, cf, lRef, pr, cr, rRef);
 
                 forAll(R.rhs(), s)
                 {
@@ -557,7 +362,7 @@ Foam::StandardChemistryModel<ReactionThermo, ThermoType>::Qdot() const
                 false
             ),
             this->mesh_,
-            dimensionedScalar("zero", dimEnergy/dimVolume/dimTime, 0.0)
+            dimensionedScalar("zero", dimEnergy/dimVolume/dimTime, 0)
         )
     );
 
@@ -600,7 +405,7 @@ Foam::StandardChemistryModel<ReactionThermo, ThermoType>::calculateRR
                 IOobject::NO_WRITE
             ),
             this->mesh(),
-            dimensionedScalar("zero", dimMass/dimVolume/dimTime, 0.0)
+            dimensionedScalar("zero", dimMass/dimVolume/dimTime, 0)
         )
     );
 
@@ -628,7 +433,7 @@ Foam::StandardChemistryModel<ReactionThermo, ThermoType>::calculateRR
         }
 
         const Reaction<ThermoType>& R = reactions_[ri];
-        const scalar omegai = omega(R, c_, Ti, pi, pf, cf, lRef, pr, cr, rRef);
+        const scalar omegai = R.omega(pi, Ti, c_, pf, cf, lRef, pr, cr, rRef);
 
         forAll(R.lhs(), s)
         {

src/thermophysicalModels/chemistryModel/chemistryModel/StandardChemistryModel/StandardChemistryModel.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -163,22 +163,6 @@ public:
             scalarField& dcdt
         ) const;
 
-        //- Return the reaction rate for reaction r and the reference
-        //  species and charateristic times
-        virtual scalar omega
-        (
-            const Reaction<ThermoType>& r,
-            const scalarField& c,
-            const scalar T,
-            const scalar p,
-            scalar& pf,
-            scalar& cf,
-            label& lRef,
-            scalar& pr,
-            scalar& cr,
-            label& rRef
-        ) const;
-
 
         //- Return the reaction rate for iReaction and the reference
         //  species and charateristic times
@@ -254,7 +238,7 @@ public:
                 const scalar t,
                 const scalarField& c,
                 scalarField& dcdt,
-                scalarSquareMatrix& dfdc
+                scalarSquareMatrix& J
             ) const;
 
             virtual void solve

src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModel.C

@@ -168,9 +168,9 @@ void Foam::TDACChemistryModel<ReactionThermo, ThermoType>::omega
         {
             const Reaction<ThermoType>& R = this->reactions_[i];
 
-            scalar omegai = omega
+            scalar omegai = R.omega
             (
-                R, c, T, p, pf, cf, lRef, pr, cr, rRef
+                p, T, c, pf, cf, lRef, pr, cr, rRef
             );
 
             forAll(R.lhs(), s)
@@ -232,17 +232,17 @@ Foam::scalar Foam::TDACChemistryModel<ReactionThermo, ThermoType>::omega
         if (c[si] < c[lRef])
         {
             const scalar exp = R.lhs()[slRef].exponent;
-            pf *= pow(max(c[lRef], 0.0), exp);
+            pf *= pow(max(c[lRef], 0), exp);
             lRef = si;
             slRef = s;
         }
         else
         {
             const scalar exp = R.lhs()[s].exponent;
-            pf *= pow(max(c[si], 0.0), exp);
+            pf *= pow(max(c[si], 0), exp);
         }
     }
-    cf = max(c[lRef], 0.0);
+    cf = max(c[lRef], 0);
 
     {
         const scalar exp = R.lhs()[slRef].exponent;
@@ -274,17 +274,17 @@ Foam::scalar Foam::TDACChemistryModel<ReactionThermo, ThermoType>::omega
         if (c[si] < c[rRef])
         {
             const scalar exp = R.rhs()[srRef].exponent;
-            pr *= pow(max(c[rRef], 0.0), exp);
+            pr *= pow(max(c[rRef], 0), exp);
             rRef = si;
             srRef = s;
         }
         else
         {
             const scalar exp = R.rhs()[s].exponent;
-            pr *= pow(max(c[si], 0.0), exp);
+            pr *= pow(max(c[si], 0), exp);
         }
     }
-    cr = max(c[rRef], 0.0);
+    cr = max(c[rRef], 0);
 
     {
         const scalar exp = R.rhs()[srRef].exponent;
@@ -334,14 +334,14 @@ void Foam::TDACChemistryModel<ReactionThermo, ThermoType>::derivatives
         // efficiencies
         for (label i=0; i<NsDAC_; i++)
         {
-            this->c_[simplifiedToCompleteIndex_[i]] = max(c[i], 0.0);
+            this->c_[simplifiedToCompleteIndex_[i]] = max(c[i], 0);
         }
     }
     else
     {
         for (label i=0; i<this->nSpecie(); i++)
         {
-            this->c_[i] = max(c[i], 0.0);
+            this->c_[i] = max(c[i], 0);
         }
     }
 
@@ -398,7 +398,8 @@ void Foam::TDACChemistryModel<ReactionThermo, ThermoType>::jacobian
 (
     const scalar t,
     const scalarField& c,
-    scalarSquareMatrix& dfdc
+    scalarField& dcdt,
+    scalarSquareMatrix& J
 ) const
 {
     const bool reduced = mechRed_->active();
@@ -416,183 +417,119 @@ void Foam::TDACChemistryModel<ReactionThermo, ThermoType>::jacobian
         this->c_ = completeC_;
         for (label i=0; i<NsDAC_; i++)
         {
-            this->c_[simplifiedToCompleteIndex_[i]] = max(c[i], 0.0);
+            this->c_[simplifiedToCompleteIndex_[i]] = max(c[i], 0);
         }
     }
     else
     {
         forAll(this->c_, i)
         {
-            this->c_[i] = max(c[i], 0.0);
+            this->c_[i] = max(c[i], 0);
         }
     }
 
-    dfdc = Zero;
+    J = Zero;
+    dcdt = Zero;
+    scalarField hi(this->c_.size());
+    scalarField cpi(this->c_.size());
+    forAll(hi, i)
+    {
+        hi[i] = this->specieThermo_[i].ha(p, T);
+        cpi[i] = this->specieThermo_[i].cp(p, T);
+    }
+
+    scalar omegaI = 0;
 
     forAll(this->reactions_, ri)
     {
         if (!reactionsDisabled_[ri])
         {
             const Reaction<ThermoType>& R = this->reactions_[ri];
+            scalar kfwd, kbwd;
+            R.dwdc
+            (
+                p,
+                T,
+                this->c_,
+                J,
+                dcdt,
+                omegaI,
+                kfwd,
+                kbwd,
+                reduced,
+                completeToSimplifiedIndex_
+            );
+            R.dwdT
+            (
+                p,
+                T,
+                this->c_,
+                omegaI,
+                kfwd,
+                kbwd,
+                J,
+                reduced,
+                completeToSimplifiedIndex_,
+                this->nSpecie_
+            );
+        }
+    }
 
-            const scalar kf0 = R.kf(p, T, this->c_);
+    // The species derivatives of the temperature term are partially computed
-            const scalar kr0 = R.kr(kf0, p, T, this->c_);
+    // while computing dwdc, they are completed hereunder:
-
+    scalar cpMean = 0;
-            forAll(R.lhs(), j)
+    scalar dcpdTMean = 0;
+    forAll(this->c_, i)
+    {
+        cpMean += this->c_[i]*cpi[i]; // J/(m3.K)
+        // Already multiplied by rho
+        dcpdTMean += this->c_[i]*this->specieThermo_[i].dcpdT(p, T);
+    }
+
+    scalar dTdt = 0;
+    forAll(hi, i)
     {
-                label sj = R.lhs()[j].index;
         if (reduced)
         {
-                    sj = completeToSimplifiedIndex_[sj];
+            const label si = completeToSimplifiedIndex_[i];
-                }
+            if (si != -1)
-                scalar kf = kf0;
-                forAll(R.lhs(), i)
             {
-                    const label si = R.lhs()[i].index;
+                dTdt += hi[i]*dcdt[si]; // J/(m3.s)
-                    const scalar el = R.lhs()[i].exponent;
-                    if (i == j)
-                    {
-                        if (el < 1)
-                        {
-                            if (this->c_[si] > small)
-                            {
-                                kf *= el*pow(this->c_[si], el - 1);
-                            }
-                            else
-                            {
-                                kf = 0;
             }
         }
         else
         {
-                            kf *= el*pow(this->c_[si], el - 1);
+            dTdt += hi[i]*dcdt[i]; // J/(m3.s)
         }
     }
-                    else
+    dTdt /= -cpMean; // K/s
+    dcdt[this->nSpecie_] = dTdt;
+
+    for (label i = 0; i < this->nSpecie_; i++)
     {
-                        kf *= pow(this->c_[si], el);
+        J(this->nSpecie_, i) = 0;
+        for (label j = 0; j < this->nSpecie_; j++)
+        {
+            const label sj = reduced ? simplifiedToCompleteIndex_[j] : j;
+            J(this->nSpecie_, i) += hi[sj]*J(j, i);
         }
+        const label si = reduced ? simplifiedToCompleteIndex_[i] : i;
+        J(this->nSpecie_, i) += cpi[si]*dTdt; // J/(mol.s)
+        J(this->nSpecie_, i) /= -cpMean;    // K/s / (mol/m3)
     }
 
-                forAll(R.lhs(), i)
+    // ddT of dTdt
+    J(this->nSpecie_, this->nSpecie_) = 0;
+    for (label i = 0; i < this->nSpecie_; i++)
     {
-                    label si = R.lhs()[i].index;
+        const label si = reduced ? simplifiedToCompleteIndex_[i] : i;
-                    if (reduced)
+        J(this->nSpecie_, this->nSpecie_) +=
-                    {
+            cpi[si]*dcdt[i]
-                        si = completeToSimplifiedIndex_[si];
+          + hi[si]*J(i, this->nSpecie_);
     }
-                    const scalar sl = R.lhs()[i].stoichCoeff;
+    J(this->nSpecie_, this->nSpecie_) += dTdt*dcpdTMean;
-                    dfdc(si, sj) -= sl*kf;
+    J(this->nSpecie_, this->nSpecie_) /= -cpMean;
-                }
+    J(this->nSpecie_, this->nSpecie_) += dTdt/T;
-                forAll(R.rhs(), i)
-                {
-                    label si = R.rhs()[i].index;
-                    if (reduced)
-                    {
-                        si = completeToSimplifiedIndex_[si];
-                    }
-                    const scalar sr = R.rhs()[i].stoichCoeff;
-                    dfdc(si, sj) += sr*kf;
-                }
-            }
-
-            forAll(R.rhs(), j)
-            {
-                label sj = R.rhs()[j].index;
-                if (reduced)
-                {
-                    sj = completeToSimplifiedIndex_[sj];
-                }
-                scalar kr = kr0;
-                forAll(R.rhs(), i)
-                {
-                    const label si = R.rhs()[i].index;
-                    const scalar er = R.rhs()[i].exponent;
-                    if (i == j)
-                    {
-                        if (er < 1)
-                        {
-                            if (this->c_[si] > small)
-                            {
-                                kr *= er*pow(this->c_[si], er - 1);
-                            }
-                            else
-                            {
-                                kr = 0;
-                            }
-                        }
-                        else
-                        {
-                            kr *= er*pow(this->c_[si], er - 1);
-                        }
-                    }
-                    else
-                    {
-                        kr *= pow(this->c_[si], er);
-                    }
-                }
-
-                forAll(R.lhs(), i)
-                {
-                    label si = R.lhs()[i].index;
-                    if (reduced)
-                    {
-                        si = completeToSimplifiedIndex_[si];
-                    }
-                    const scalar sl = R.lhs()[i].stoichCoeff;
-                    dfdc(si, sj) += sl*kr;
-                }
-                forAll(R.rhs(), i)
-                {
-                    label si = R.rhs()[i].index;
-                    if (reduced)
-                    {
-                        si = completeToSimplifiedIndex_[si];
-                    }
-                    const scalar sr = R.rhs()[i].stoichCoeff;
-                    dfdc(si, sj) -= sr*kr;
-                }
-            }
-        }
-    }
-
-    // Calculate the dcdT elements numerically
-    const scalar delta = 1e-3;
-
-    omega(this->c_, T + delta, p, this->dcdt_);
-    for (label i=0; i<this->nSpecie_; i++)
-    {
-        dfdc(i, this->nSpecie_) = this->dcdt_[i];
-    }
-
-    omega(this->c_, T - delta, p, this->dcdt_);
-    for (label i=0; i<this->nSpecie_; i++)
-    {
-        dfdc(i, this->nSpecie_) =
-            0.5*(dfdc(i, this->nSpecie_) - this->dcdt_[i])/delta;
-    }
-
-    dfdc(this->nSpecie_, this->nSpecie_) = 0;
-    dfdc(this->nSpecie_ + 1, this->nSpecie_) = 0;
-}
-
-
-template<class ReactionThermo, class ThermoType>
-void Foam::TDACChemistryModel<ReactionThermo, ThermoType>::jacobian
-(
-    const scalar t,
-    const scalarField& c,
-    scalarField& dcdt,
-    scalarSquareMatrix& dfdc
-) const
-{
-    jacobian(t, c, dfdc);
-
-    const scalar T = c[this->nSpecie_];
-    const scalar p = c[this->nSpecie_ + 1];
-
-    // Note: Uses the c_ field initialized by the call to jacobian above
-    omega(this->c_, T, p, dcdt);
 }
 
 
@@ -915,7 +852,7 @@ setTabulationResultsAdd
     const label celli
 )
 {
-    tabulationResults_[celli] = 0.0;
+    tabulationResults_[celli] = 0;
 }
 
 
@@ -923,7 +860,7 @@ template<class ReactionThermo, class ThermoType>
 void Foam::TDACChemistryModel<ReactionThermo, ThermoType>::
 setTabulationResultsGrow(const label celli)
 {
-    tabulationResults_[celli] = 1.0;
+    tabulationResults_[celli] = 1;
 }
 
 
@@ -934,7 +871,7 @@ setTabulationResultsRetrieve
     const label celli
 )
 {
-    tabulationResults_[celli] = 2.0;
+    tabulationResults_[celli] = 2;
 }
 
 

src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/TDACChemistryModel.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2016-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2016-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -221,20 +221,12 @@ public:
                 scalarField& dcdt
             ) const;
 
-            //- Pure jacobian function for tabulation
-            void jacobian
-            (
-                const scalar t,
-                const scalarField& c,
-                scalarSquareMatrix& dfdc
-            ) const;
-
             virtual void jacobian
             (
                 const scalar t,
                 const scalarField& c,
                 scalarField& dcdt,
-                scalarSquareMatrix& dfdc
+                scalarSquareMatrix& J
             ) const;
 
             virtual void solve

src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/ISAT/ISAT.C

@@ -65,7 +65,7 @@ Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::ISAT
     (
         this->coeffsDict_.lookupOrDefault
         (
-            "minBalanceThreshold",0.1*chemisTree_.maxNLeafs()
+            "minBalanceThreshold", 0.1*chemisTree_.maxNLeafs()
         )
     ),
     MRURetrieve_(this->coeffsDict_.lookupOrDefault("MRURetrieve", false)),
@@ -168,20 +168,11 @@ void Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::addToMRU
         else // chemPoint not yet in the list, iter is last
         {
             if (MRUList_.size() == maxMRUSize_)
-            {
-                if (iter() == MRUList_.last())
             {
                 MRUList_.remove(iter);
                 MRUList_.insert(phi0);
             }
             else
-                {
-                    FatalErrorInFunction
-                        << "wrong MRUList construction"
-                        << exit(FatalError);
-                }
-            }
-            else
             {
                 MRUList_.insert(phi0);
             }
@@ -238,13 +229,13 @@ void Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::calcNewC
 
                 // As we use an approximation of A, Rphiq should be checked for
                 // negative values
-                Rphiq[i] = max(0.0,Rphiq[i]);
+                Rphiq[i] = max(0, Rphiq[i]);
             }
             // The species is not active A(i, j) = I(i, j)
             else
             {
                 Rphiq[i] += dphi[i];
-                Rphiq[i] = max(0.0,Rphiq[i]);
+                Rphiq[i] = max(0, Rphiq[i]);
             }
         }
         else // Mechanism reduction is not active
@@ -255,7 +246,7 @@ void Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::calcNewC
             }
             // As we use a first order gradient matrix, Rphiq should be checked
             // for negative values
-            Rphiq[i] = max(0.0,Rphiq[i]);
+            Rphiq[i] = max(0, Rphiq[i]);
         }
     }
 }
@@ -322,6 +313,9 @@ Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::cleanAndBalance()
         }
         x = xtmp;
     }
+
+    MRUList_.clear();
+
     // Check if the tree should be balanced according to criterion:
     //  -the depth of the tree bigger than a*log2(size), log2(size) being the
     //      ideal depth (e.g. 4 leafs can be stored in a tree of depth 2)
@@ -333,7 +327,6 @@ Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::cleanAndBalance()
     )
     {
         chemisTree_.balance();
-        MRUList_.clear();
         treeModified = true;
     }
 
@@ -365,7 +358,7 @@ void Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::computeA
         Rcq[i] = rhoi*Rphiq[s2c]/this->chemistry_.specieThermo()[s2c].W();
     }
     Rcq[speciesNumber] = Rphiq[Rphiq.size() - nAdditionalEqns_];
-    Rcq[speciesNumber+1] = Rphiq[Rphiq.size() - nAdditionalEqns_ + 1];
+    Rcq[speciesNumber + 1] = Rphiq[Rphiq.size() - nAdditionalEqns_ + 1];
     if (this->variableTimeStep())
     {
         Rcq[speciesNumber + 2] = Rphiq[Rphiq.size() - nAdditionalEqns_ + 2];
@@ -373,15 +366,15 @@ void Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::computeA
 
     // Aaa is computed implicitly,
     // A is given by A = C(psi0, t0+dt), where C is obtained through solving
-    // d/dt C(psi0,t) = J(psi(t))C(psi0,t)
+    // d/dt C(psi0, t) = J(psi(t))C(psi0, t)
     // If we solve it implicitly:
-    // (C(psi0, t0+dt) - C(psi0,t0))/dt = J(psi(t0+dt))C(psi0,t0+dt)
+    // (C(psi0, t0+dt) - C(psi0, t0))/dt = J(psi(t0+dt))C(psi0, t0+dt)
     // The Jacobian is thus computed according to the mapping
     // C(psi0,t0+dt)*(I-dt*J(psi(t0+dt))) = C(psi0, t0)
     // A = C(psi0,t0)/(I-dt*J(psi(t0+dt)))
     // where C(psi0,t0) = I
-
+    scalarField dcdt(speciesNumber + 2, Zero);
-    this->chemistry_.jacobian(runTime_.value(), Rcq, A);
+    this->chemistry_.jacobian(runTime_.value(), Rcq, dcdt, A);
 
     // The jacobian is computed according to the molar concentration
     // the following conversion allows the code to use A with mass fraction
@@ -410,21 +403,48 @@ void Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::computeA
         // Columns for pressure and temperature
         A(i, speciesNumber) *=
             -dt*this->chemistry_.specieThermo()[si].W()/rhoi;
-        A(i, speciesNumber+1) *=
+        A(i, speciesNumber + 1) *=
             -dt*this->chemistry_.specieThermo()[si].W()/rhoi;
     }
 
-    // For temperature and pressure, only unity on the diagonal
+    // For the temperature and pressure lines, ddc(dTdt)
-    A(speciesNumber, speciesNumber) = 1;
+    // should be converted in ddY(dTdt)
-    A(speciesNumber + 1, speciesNumber + 1) = 1;
+    for (label i=0; i<speciesNumber; i++)
+    {
+        label si = i;
+        if (mechRedActive)
+        {
+            si = this->chemistry_.simplifiedToCompleteIndex()[i];
+        }
+
+        A(speciesNumber, i) *=
+            -dt*rhoi/this->chemistry_.specieThermo()[si].W();
+        A(speciesNumber + 1, i) *=
+            -dt*rhoi/this->chemistry_.specieThermo()[si].W();
+    }
+
+    A(speciesNumber, speciesNumber) = -dt*A(speciesNumber, speciesNumber) + 1;
+
+    A(speciesNumber + 1, speciesNumber + 1) =
+        -dt*A(speciesNumber + 1, speciesNumber + 1) + 1;
+
     if (this->variableTimeStep())
     {
-        A[speciesNumber + 2][speciesNumber + 2] = 1;
+        A(speciesNumber + 2, speciesNumber + 2) = 1;
     }
 
     // Inverse of (I-dt*J(psi(t0+dt)))
     LUscalarMatrix LUA(A);
     LUA.inv(A);
+
+    // After inversion, lines of p and T are set to 0 except diagonal.  This
+    // avoid skewness of the ellipsoid of accuracy and potential issues in the
+    // binary tree.
+    for (label i=0; i<speciesNumber; i++)
+    {
+        A(speciesNumber, i) = 0;
+        A(speciesNumber + 1, i) = 0;
+    }
 }
 
 
@@ -498,7 +518,7 @@ bool Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::retrieve
         }
         lastSearch_->lastTimeUsed() = this->chemistry_.timeSteps();
         addToMRU(phi0);
-        calcNewC(phi0,phiq, Rphiq);
+        calcNewC(phi0, phiq, Rphiq);
         nRetrieved_++;
         return true;
     }
@@ -525,11 +545,12 @@ Foam::label Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::add
     // option is on, the code first tries to grow the point hold by lastSearch_
     if (lastSearch_ && growPoints_)
     {
-        if (grow(lastSearch_,phiq, Rphiq))
+        if (grow(lastSearch_, phiq, Rphiq))
         {
             nGrowth_++;
             growthOrAddFlag = 0;
-            // the structure of the tree is not modified, return false
+            addToMRU(lastSearch_);
+            //the structure of the tree is not modified, return false
             return growthOrAddFlag;
         }
     }
@@ -605,7 +626,10 @@ Foam::label Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::add
         scaleFactor_.size(),
         lastSearch_ // lastSearch_ may be nullptr (handled by binaryTree)
     );
-
+    if (lastSearch_ != nullptr)
+    {
+        addToMRU(lastSearch_);
+    }
     nAdd_++;
 
     return growthOrAddFlag;

src/thermophysicalModels/radiation/derivedFvPatchFields/radiationCoupledBase/radiationCoupledBase.C

@@ -28,7 +28,7 @@ License
 #include "mappedPatchBase.H"
 #include "fvPatchFieldMapper.H"
 #include "radiationModel.H"
-#include "noRadiation.H"
+#include "opaqueSolid.H"
 #include "absorptionEmissionModel.H"
 
 // * * * * * * * * * * * * * Static Member Data  * * * * * * * * * * * * * * //
@@ -151,29 +151,20 @@ Foam::scalarField Foam::radiationCoupledBase::emissivity() const
 
             const polyMesh& nbrMesh = mpp.sampleMesh();
 
-            const radiation::radiationModel& radiation =
+            const radiation::opaqueSolid& radiation =
-                nbrMesh.lookupObject<radiation::radiationModel>
+                nbrMesh.lookupObject<radiation::opaqueSolid>
                 (
                     "radiationProperties"
                 );
 
-            if (isType<radiation::noRadiation>(radiation))
-            {
-                FatalErrorInFunction
-                    << "No radiation model defined for region "
-                    << nbrMesh.name() << nl
-                    << "    required for option 'solidRadiation' "
-                       "of patchField type " << type()
-                    << " on patch " << patch_.name()
-                    << abort(FatalError);
-            }
-
             const fvMesh& nbrFvMesh = refCast<const fvMesh>(nbrMesh);
 
             const fvPatch& nbrPatch =
                 nbrFvMesh.boundary()[mpp.samplePolyPatch().index()];
 
-
+            // NOTE: for an opaqueSolid the absorptionEmission model returns the
+            // emissivity of the surface rather than the emission coefficient
+            // and the input specification MUST correspond to this.
             scalarField emissivity
             (
                 radiation.absorptionEmission().e()().boundaryField()

src/thermophysicalModels/reactionThermo/chemistryReaders/chemkinReader/chemkinLexer.L

@@ -651,6 +651,18 @@ bool finishReaction = false;
     }
 
 <readThermoSpecieName>{end} {
+        Reaction<gasHThermoPhysics>::TlowDefault = max
+        (
+            Reaction<gasHThermoPhysics>::TlowDefault,
+            currentLowT
+        );
+
+        Reaction<gasHThermoPhysics>::ThighDefault = min
+        (
+            Reaction<gasHThermoPhysics>::ThighDefault,
+            currentHighT
+        );
+
         BEGIN(INITIAL);
     }
 

src/thermophysicalModels/reactionThermo/chemistryReaders/chemkinReader/chemkinReader.C

@@ -26,6 +26,7 @@ License
 #include "chemkinReader.H"
 #include <fstream>
 #include "atomicWeights.H"
+#include "ReactionProxy.H"
 #include "IrreversibleReaction.H"
 #include "ReversibleReaction.H"
 #include "NonEquilibriumReversibleReaction.H"
@@ -182,7 +183,7 @@ void Foam::chemkinReader::addReactionType
                 new IrreversibleReaction
                 <Reaction, gasHThermoPhysics, ReactionRateType>
                 (
-                    Reaction<gasHThermoPhysics>
+                    ReactionProxy<gasHThermoPhysics>
                     (
                         speciesTable_,
                         lhs.shrink(),
@@ -202,7 +203,7 @@ void Foam::chemkinReader::addReactionType
                 new ReversibleReaction
                 <Reaction, gasHThermoPhysics, ReactionRateType>
                 (
-                    Reaction<gasHThermoPhysics>
+                    ReactionProxy<gasHThermoPhysics>
                     (
                         speciesTable_,
                         lhs.shrink(),
@@ -490,7 +491,7 @@ void Foam::chemkinReader::addReaction
                     new NonEquilibriumReversibleReaction
                         <Reaction, gasHThermoPhysics, ArrheniusReactionRate>
                     (
-                        Reaction<gasHThermoPhysics>
+                        ReactionProxy<gasHThermoPhysics>
                         (
                             speciesTable_,
                             lhs.shrink(),
@@ -546,7 +547,7 @@ void Foam::chemkinReader::addReaction
                         thirdBodyArrheniusReactionRate
                     >
                     (
-                        Reaction<gasHThermoPhysics>
+                        ReactionProxy<gasHThermoPhysics>
                         (
                             speciesTable_,
                             lhs.shrink(),
@@ -646,9 +647,13 @@ void Foam::chemkinReader::addReaction
                 reactions_.append
                 (
                     new NonEquilibriumReversibleReaction
-                        <Reaction, gasHThermoPhysics, LandauTellerReactionRate>
+                    <
+                        Reaction,
+                        gasHThermoPhysics,
+                        LandauTellerReactionRate
+                    >
                     (
-                        Reaction<gasHThermoPhysics>
+                        ReactionProxy<gasHThermoPhysics>
                         (
                             speciesTable_,
                             lhs.shrink(),
@@ -780,6 +785,9 @@ void Foam::chemkinReader::read
     const fileName& transportFileName
 )
 {
+    Reaction<gasHThermoPhysics>::TlowDefault = 0;
+    Reaction<gasHThermoPhysics>::ThighDefault = great;
+
     transportDict_.read(IFstream(transportFileName)());
 
     if (thermoFileName != fileName::null)

src/thermophysicalModels/solidSpecie/reaction/Reactions/solidReaction/solidReaction.H

@@ -129,28 +129,6 @@ public:
         );
 
 
-        //- Construct and return a clone
-        virtual autoPtr<Reaction<ReactionThermo>> clone() const
-        {
-            return autoPtr<Reaction<ReactionThermo>>
-            (
-                new solidReaction<ReactionThermo>(*this)
-            );
-        }
-
-        //- Construct and return a clone with new speciesTable
-        virtual autoPtr<Reaction<ReactionThermo>> clone
-        (
-            const speciesTable& species
-        ) const
-        {
-            return autoPtr<Reaction<ReactionThermo>>
-            (
-                new solidReaction<ReactionThermo>(*this, species)
-            );
-        }
-
-
     //- Destructor
     virtual ~solidReaction()
     {}
@@ -160,12 +138,13 @@ public:
 
         // Access
 
-
+            //- Access to the gas components of the left hand side
-            // - Access to gas components of the reaction
             virtual const List<specieCoeffs>& grhs() const;
+
+            //- Access to the gas components of the right hand side
             virtual const List<specieCoeffs>& glhs() const;
 
-            // - Access to gas specie list
+            //- Access to gas specie list
             virtual const speciesTable& gasSpecies() const;
 
 

src/thermophysicalModels/solidSpecie/reaction/reactionRate/solidArrheniusReactionRate/solidArrheniusReactionRate.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -104,6 +104,37 @@ public:
             const scalarField& c
         ) const;
 
+        inline scalar ddT
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c
+        ) const;
+
+        //- Third-body efficiencies (beta = 1-alpha)
+        //  non-empty only for third-body reactions
+        //  with enhanced molecularity (alpha != 1)
+        inline const List<Tuple2<label, scalar>>& beta() const;
+
+        //- Species concentration derivative of the pressure dependent term
+        //  By default this value is 1 as it multiplies the third-body term
+        inline void dcidc
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c,
+            scalarField& dcidc
+        ) const;
+
+        //- Temperature derivative of the pressure dependent term
+        //  By default this value is 0 since ddT of molecularity is approx.0
+        inline scalar dcidT
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c
+        ) const;
+
 
          //- Write to stream
         inline void write(Ostream& os) const;

src/thermophysicalModels/solidSpecie/reaction/reactionRate/solidArrheniusReactionRate/solidArrheniusReactionRateI.H

@@ -60,18 +60,60 @@ inline Foam::scalar Foam::solidArrheniusReactionRate::operator()
     const scalarField&
 ) const
 {
-    scalar ak = A_;
-
     if (T < Tcrit_)
     {
-        ak *= 0.0;
+        return 0;
     }
     else
     {
-        ak *= exp(-Ta_/T);
+        return A_*exp(-Ta_/T);
     }
+}
 
-    return ak;
+
+inline Foam::scalar Foam::solidArrheniusReactionRate::ddT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField&
+) const
+{
+    if (T < Tcrit_)
+    {
+        return 0;
+    }
+    else
+    {
+        return A_*exp(-Ta_/T)*Ta_/sqr(T);
+    }
+}
+
+
+inline const Foam::List<Foam::Tuple2<Foam::label, Foam::scalar>>&
+Foam::solidArrheniusReactionRate::beta() const
+{
+    return NullObjectRef<List<Tuple2<label, scalar>>>();
+}
+
+
+inline void Foam::solidArrheniusReactionRate::dcidc
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c,
+    scalarField& dcidc
+) const
+{}
+
+
+inline Foam::scalar Foam::solidArrheniusReactionRate::dcidT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    return 0;
 }
 
 

src/thermophysicalModels/specie/reaction/Reactions/IrreversibleReaction/IrreversibleReaction.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -106,6 +106,160 @@ Foam::scalar Foam::IrreversibleReaction
 }
 
 
+template
+<
+    template<class> class ReactionType,
+    class ReactionThermo,
+    class ReactionRate
+>
+Foam::scalar Foam::IrreversibleReaction
+<
+    ReactionType,
+    ReactionThermo,
+    ReactionRate
+>::kr
+(
+    const scalar kfwd,
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    return 0;
+}
+
+
+template
+<
+    template<class> class ReactionType,
+    class ReactionThermo,
+    class ReactionRate
+>
+Foam::scalar Foam::IrreversibleReaction
+<
+    ReactionType,
+    ReactionThermo,
+    ReactionRate
+>::kr
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    return 0;
+}
+
+
+template
+<
+    template<class> class ReactionType,
+    class ReactionThermo,
+    class ReactionRate
+>
+Foam::scalar Foam::IrreversibleReaction
+<
+    ReactionType,
+    ReactionThermo,
+    ReactionRate
+>::dkfdT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    return k_.ddT(p, T, c);
+}
+
+
+template
+<
+    template<class> class ReactionType,
+    class ReactionThermo,
+    class ReactionRate
+>
+Foam::scalar Foam::IrreversibleReaction
+<
+    ReactionType,
+    ReactionThermo,
+    ReactionRate
+>::dkrdT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c,
+    const scalar dkfdT,
+    const scalar kr
+) const
+{
+    return 0;
+}
+
+
+template
+<
+    template<class> class ReactionType,
+    class ReactionThermo,
+    class ReactionRate
+>
+const Foam::List<Foam::Tuple2<Foam::label, Foam::scalar>>&
+Foam::IrreversibleReaction
+<
+    ReactionType,
+    ReactionThermo,
+    ReactionRate
+>::beta() const
+{
+    return k_.beta();
+}
+
+
+template
+<
+    template<class> class ReactionType,
+    class ReactionThermo,
+    class ReactionRate
+>
+void Foam::IrreversibleReaction
+<
+    ReactionType,
+    ReactionThermo,
+    ReactionRate
+>::dcidc
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c,
+    scalarField& dcidc
+) const
+{
+    k_.dcidc(p, T, c, dcidc);
+}
+
+
+template
+<
+    template<class> class ReactionType,
+    class ReactionThermo,
+    class ReactionRate
+>
+Foam::scalar Foam::IrreversibleReaction
+<
+    ReactionType,
+    ReactionThermo,
+    ReactionRate
+>::dcidT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    return k_.dcidT(p, T, c);
+}
+
+
 template
 <
     template<class> class ReactionType,

src/thermophysicalModels/specie/reaction/Reactions/IrreversibleReaction/IrreversibleReaction.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,8 +25,7 @@ Class
     Foam::IrreversibleReaction
 
 Description
-    Simple extension of Reaction to handle reversible reactions using
+    Simple extension of Reaction to handle irreversible reactions
-    equilibrium thermodynamics.
 
 SourceFiles
     IrreversibleReaction.C
@@ -164,6 +163,71 @@ public:
                 const scalarField& c
             ) const;
 
+            //- Reverse rate constant from the given forward rate constant
+            //  Returns 0
+            virtual scalar kr
+            (
+                const scalar kfwd,
+                const scalar p,
+                const scalar T,
+                const scalarField& c
+            ) const;
+
+            //- Reverse rate constant
+            //  Returns 0
+            virtual scalar kr
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c
+            ) const;
+
+
+        // IrreversibleReaction Jacobian functions
+
+            //- Temperature derivative of forward rate
+            virtual scalar dkfdT
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c
+            ) const;
+
+            //- Temperature derivative of reverse rate
+            //  Returns 0
+            virtual scalar dkrdT
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c,
+                const scalar dkfdT,
+                const scalar kr
+            ) const;
+
+            //- Third-body efficiencies (beta = 1-alpha)
+            //  non-empty only for third-body reactions
+            //  with enhanced molecularity (alpha != 1)
+            virtual const List<Tuple2<label, scalar>>& beta() const;
+
+            //- Species concentration derivative of the pressure dependent term
+            //  By default this value is 1 as it multiplies the third-body term
+            virtual void dcidc
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c,
+                scalarField& dcidc
+            ) const;
+
+            //- Temperature derivative of the pressure dependent term
+            //  By default this value is 0 since ddT of molecularity is approx.0
+            virtual scalar dcidT
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c
+            ) const;
+
 
         //- Write
         virtual void write(Ostream&) const;

src/thermophysicalModels/specie/reaction/Reactions/NonEquilibriumReversibleReaction/NonEquilibriumReversibleReaction.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -178,6 +178,117 @@ Foam::NonEquilibriumReversibleReaction
 }
 
 
+template
+<
+    template<class> class ReactionType,
+    class ReactionThermo,
+    class ReactionRate
+>
+Foam::scalar
+Foam::NonEquilibriumReversibleReaction
+<
+    ReactionType,
+    ReactionThermo,
+    ReactionRate
+>::dkfdT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    return fk_.ddT(p, T, c);
+}
+
+
+template
+<
+    template<class> class ReactionType,
+    class ReactionThermo,
+    class ReactionRate
+>
+Foam::scalar
+Foam::NonEquilibriumReversibleReaction
+<
+    ReactionType,
+    ReactionThermo,
+    ReactionRate
+>::dkrdT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c,
+    const scalar dkfdT,
+    const scalar kr
+) const
+{
+    return rk_.ddT(p, T, c);
+}
+
+
+template
+<
+    template<class> class ReactionType,
+    class ReactionThermo,
+    class ReactionRate
+>
+const Foam::List<Foam::Tuple2<Foam::label, Foam::scalar>>&
+Foam::NonEquilibriumReversibleReaction
+<
+    ReactionType,
+    ReactionThermo,
+    ReactionRate
+>::beta() const
+{
+    return rk_.beta();
+}
+
+
+template
+<
+    template<class> class ReactionType,
+    class ReactionThermo,
+    class ReactionRate
+>
+void Foam::NonEquilibriumReversibleReaction
+<
+    ReactionType,
+    ReactionThermo,
+    ReactionRate
+>::dcidc
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c,
+    scalarField& dcidc
+) const
+{
+    fk_.dcidc(p, T, c, dcidc);
+}
+
+
+template
+<
+    template<class> class ReactionType,
+    class ReactionThermo,
+    class ReactionRate
+>
+Foam::scalar Foam::NonEquilibriumReversibleReaction
+<
+    ReactionType,
+    ReactionThermo,
+    ReactionRate
+>::dcidT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    return fk_.dcidT(p, T, c);
+}
+
+
 template
 <
     template<class> class ReactionType,

src/thermophysicalModels/specie/reaction/Reactions/NonEquilibriumReversibleReaction/NonEquilibriumReversibleReaction.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -167,6 +167,51 @@ public:
             ) const;
 
 
+        // ReversibleReaction Jacobian functions
+
+            //- Temperature derivative of forward rate
+            virtual scalar dkfdT
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c
+            ) const;
+
+            //- Temperature derivative of backward rate
+            virtual scalar dkrdT
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c,
+                const scalar dkfdT,
+                const scalar kr
+            ) const;
+
+            //- Third-body efficiencies (beta = 1-alpha)
+            //  non-empty only for third-body reactions
+            //  with enhanced molecularity (alpha != 1)
+            virtual const List<Tuple2<label, scalar>>& beta() const;
+
+            //- Species concentration derivative of the pressure dependent term
+            //  By default this value is 1 as it multiplies the third-body term
+            virtual void dcidc
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c,
+                scalarField& dcidc
+            ) const;
+
+            //- Temperature derivative of the pressure dependent term
+            //  By default this value is 0 since ddT of molecularity is approx.0
+            virtual scalar dcidT
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c
+            ) const;
+
+
         //- Write
         virtual void write(Ostream&) const;
 };

src/thermophysicalModels/specie/reaction/Reactions/Reaction/Reaction.C

@@ -29,7 +29,13 @@ License
 // * * * * * * * * * * * * * * * * Static Data * * * * * * * * * * * * * * * //
 
 template<class ReactionThermo>
-Foam::label Foam::Reaction<ReactionThermo>::nUnNamedReactions = 0;
+Foam::label Foam::Reaction<ReactionThermo>::nUnNamedReactions(0);
+
+template<class ReactionThermo>
+Foam::scalar Foam::Reaction<ReactionThermo>::TlowDefault(0);
+
+template<class ReactionThermo>
+Foam::scalar Foam::Reaction<ReactionThermo>::ThighDefault(great);
 
 // * * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * //
 
@@ -160,6 +166,8 @@ Foam::Reaction<ReactionThermo>::Reaction
     ReactionThermo::thermoType(*thermoDatabase[species[0]]),
     name_("un-named-reaction-" + Foam::name(getNewReactionID())),
     species_(species),
+    Tlow_(TlowDefault),
+    Thigh_(ThighDefault),
     lhs_(lhs),
     rhs_(rhs)
 {
@@ -177,6 +185,8 @@ Foam::Reaction<ReactionThermo>::Reaction
     ReactionThermo::thermoType(r),
     name_(r.name() + "Copy"),
     species_(species),
+    Tlow_(r.Tlow()),
+    Thigh_(r.Thigh()),
     lhs_(r.lhs_),
     rhs_(r.rhs_)
 {}
@@ -197,7 +207,7 @@ Foam::Reaction<ReactionThermo>::specieCoeffs::specieCoeffs
     }
     else
     {
-        stoichCoeff = 1.0;
+        stoichCoeff = 1;
     }
 
     exponent = stoichCoeff;
@@ -330,7 +340,9 @@ Foam::Reaction<ReactionThermo>::Reaction
 :
     ReactionThermo::thermoType(*thermoDatabase[species[0]]),
     name_(dict.dictName()),
-    species_(species)
+    species_(species),
+    Tlow_(dict.lookupOrDefault<scalar>("Tlow", TlowDefault)),
+    Thigh_(dict.lookupOrDefault<scalar>("Thigh", ThighDefault))
 {
     setLRhs
     (
@@ -388,39 +400,402 @@ void Foam::Reaction<ReactionThermo>::write(Ostream& os) const
 
 
 template<class ReactionThermo>
-Foam::scalar Foam::Reaction<ReactionThermo>::kf
+void Foam::Reaction<ReactionThermo>::ddot
 (
     const scalar p,
     const scalar T,
-    const scalarField& c
+    const scalarField& c,
+    scalarField& d
 ) const
 {
-    return 0.0;
 }
 
 
 template<class ReactionThermo>
-Foam::scalar Foam::Reaction<ReactionThermo>::kr
+void Foam::Reaction<ReactionThermo>::fdot
 (
+    const scalar p,
+    const scalar T,
+    const scalarField& c,
+    scalarField& f
+) const
+{
+}
+
+
+template<class ReactionThermo>
+void Foam::Reaction<ReactionThermo>::omega
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c,
+    scalarField& dcdt
+) const
+{
+    scalar pf, cf, pr, cr;
+    label lRef, rRef;
+
+    scalar omegaI = omega
+    (
+        p, T, c, pf, cf, lRef, pr, cr, rRef
+    );
+
+    forAll(lhs_, i)
+    {
+        const label si = lhs_[i].index;
+        const scalar sl = lhs_[i].stoichCoeff;
+        dcdt[si] -= sl*omegaI;
+    }
+    forAll(rhs_, i)
+    {
+        const label si = rhs_[i].index;
+        const scalar sr = rhs_[i].stoichCoeff;
+        dcdt[si] += sr*omegaI;
+    }
+}
+
+
+template<class ReactionThermo>
+Foam::scalar Foam::Reaction<ReactionThermo>::omega
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c,
+    scalar& pf,
+    scalar& cf,
+    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);
+    const scalar kr = this->kr(kf, p, clippedT, c);
+
+    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++)
+    {
+        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);
+        }
+    }
+    cf = max(c[lRef], 0);
+
+    {
+        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);
+        }
+    }
+
+    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;
+}
+
+
+template<class ReactionThermo>
+void Foam::Reaction<ReactionThermo>::dwdc
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c,
+    scalarSquareMatrix& J,
+    scalarField& dcdt,
+    scalar& omegaI,
+    scalar& kfwd,
+    scalar& kbwd,
+    const bool reduced,
+    const List<label>& c2s
+) const
+{
+
+    scalar pf, cf, pr, cr;
+    label lRef, rRef;
+
+    omegaI = omega(p, T, c, pf, cf, lRef, pr, cr, rRef);
+
+    forAll(lhs_, i)
+    {
+        const label si = reduced ? c2s[lhs_[i].index] : lhs_[i].index;
+        const scalar sl = lhs_[i].stoichCoeff;
+        dcdt[si] -= sl*omegaI;
+    }
+    forAll(rhs_, i)
+    {
+        const label si = reduced ? c2s[rhs_[i].index] : rhs_[i].index;
+        const scalar sr = rhs_[i].stoichCoeff;
+        dcdt[si] += sr*omegaI;
+    }
+
+    kfwd = this->kf(p, T, c);
+    kbwd = this->kr(kfwd, p, T, c);
+
+    forAll(lhs_, j)
+    {
+        const label sj = reduced ? c2s[lhs_[j].index] : lhs_[j].index;
+        scalar kf = kfwd;
+        forAll(lhs_, i)
+        {
+            const label si = lhs_[i].index;
+            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);
+                }
+            }
+            else
+            {
+                kf *= pow(c[si], el);
+            }
+        }
+
+        forAll(lhs_, i)
+        {
+            const label si = reduced ? c2s[lhs_[i].index] : lhs_[i].index;
+            const scalar sl = lhs_[i].stoichCoeff;
+            J(si, sj) -= sl*kf;
+        }
+        forAll(rhs_, i)
+        {
+            const label si = reduced ? c2s[rhs_[i].index] : rhs_[i].index;
+            const scalar sr = rhs_[i].stoichCoeff;
+            J(si, sj) += sr*kf;
+        }
+    }
+
+    forAll(rhs_, j)
+    {
+        const label sj = reduced ? c2s[rhs_[j].index] : rhs_[j].index;
+        scalar kr = kbwd;
+        forAll(rhs_, i)
+        {
+            const label si = rhs_[i].index;
+            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);
+                }
+            }
+            else
+            {
+                kr *= pow(c[si], er);
+            }
+        }
+
+        forAll(lhs_, i)
+        {
+            const label si = reduced ? c2s[lhs_[i].index] : lhs_[i].index;
+            const scalar sl = lhs_[i].stoichCoeff;
+            J(si, sj) += sl*kr;
+        }
+        forAll(rhs_, i)
+        {
+            const label si = reduced ? c2s[rhs_[i].index] : rhs_[i].index;
+            const scalar sr = rhs_[i].stoichCoeff;
+            J(si, sj) -= sr*kr;
+        }
+    }
+
+    // When third-body species are involved, additional terms are added
+    // beta function returns an empty list when third-body are not involved
+    if (this->beta().size() > 0)
+    {
+        scalarField dcidc(this->beta().size());
+        this->dcidc(p, T, c, dcidc);
+        forAll(this->beta(), j)
+        {
+            label sj = this-> beta()[j].first();
+            sj = reduced ? c2s[sj] : sj;
+            if (sj != -1)
+            {
+                forAll(lhs_, i)
+                {
+                    const label si =
+                        reduced ? c2s[lhs_[i].index] : lhs_[i].index;
+                    const scalar sl = lhs_[i].stoichCoeff;
+                    J(si, sj) -= sl*dcidc[j]*omegaI;
+                }
+                forAll(rhs_, i)
+                {
+                    const label si =
+                        reduced ? c2s[rhs_[i].index] : rhs_[i].index;
+                    const scalar sr = rhs_[i].stoichCoeff;
+                    J(si, sj) += sr*dcidc[j]*omegaI;
+                }
+            }
+        }
+    }
+}
+
+
+template<class ReactionThermo>
+void Foam::Reaction<ReactionThermo>::dwdT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c,
+    const scalar omegaI,
     const scalar kfwd,
-    const scalar p,
+    const scalar kbwd,
-    const scalar T,
+    scalarSquareMatrix& J,
-    const scalarField& c
+    const bool reduced,
+    const List<label>& c2s,
+    const label indexT
 ) const
 {
-    return 0.0;
+    scalar kf = kfwd;
-}
+    scalar kr = kbwd;
 
+    scalar dkfdT = this->dkfdT(p, T, c);
+    scalar dkrdT = this->dkrdT(p, T, c, dkfdT, kr);
 
-template<class ReactionThermo>
+    scalar sumExp = 0.0;
-Foam::scalar Foam::Reaction<ReactionThermo>::kr
+    forAll(lhs_, i)
-(
+    {
-    const scalar p,
+        const label si = lhs_[i].index;
-    const scalar T,
+        const scalar el = lhs_[i].exponent;
-    const scalarField& c
+        const scalar cExp = pow(c[si], el);
-) const
+        dkfdT *= cExp;
-{
+        kf *= cExp;
-    return 0.0;
+        sumExp += el;
+    }
+    kf *= -sumExp/T;
+
+    sumExp = 0.0;
+    forAll(rhs_, i)
+    {
+        const label si = rhs_[i].index;
+        const scalar er = rhs_[i].exponent;
+        const scalar cExp = pow(c[si], er);
+        dkrdT *= cExp;
+        kr *= cExp;
+        sumExp += er;
+    }
+    kr *= -sumExp/T;
+
+    // dqidT includes the third-body (or pressure dependent) effect
+    scalar dqidT = dkfdT - dkrdT + kf - kr;
+
+    // For reactions including third-body efficiencies or pressure dependent
+    // reaction, an additional term is needed
+    scalar dcidT = this->dcidT(p, T, c);
+    dcidT *= omegaI;
+
+    // J(i, indexT) = sum_reactions nu_i dqdT
+    forAll(lhs_, i)
+    {
+        const label si = reduced ? c2s[lhs_[i].index] : lhs_[i].index;
+        const scalar sl = lhs_[i].stoichCoeff;
+        J(si, indexT) -= sl*(dqidT + dcidT);
+    }
+    forAll(rhs_, i)
+    {
+        const label si = reduced ? c2s[rhs_[i].index] : rhs_[i].index;
+        const scalar sr = rhs_[i].stoichCoeff;
+        J(si, indexT) += sr*(dqidT + dcidT);
+    }
 }
 
 

src/thermophysicalModels/specie/reaction/Reactions/Reaction/Reaction.H

@@ -40,6 +40,8 @@ SourceFiles
 #include "speciesTable.H"
 #include "HashPtrTable.H"
 #include "scalarField.H"
+#include "simpleMatrix.H"
+#include "Tuple2.H"
 #include "typeInfo.H"
 #include "runTimeSelectionTables.H"
 
@@ -84,6 +86,8 @@ public:
         //- Number of un-named reactions
         static label nUnNamedReactions;
 
+        //- Default temperature limits of applicability of reaction rates
+        static scalar TlowDefault, ThighDefault;
 
     // Public data types
 
@@ -134,6 +138,9 @@ private:
         //- List of specie names present in reaction system
         const speciesTable& species_;
 
+        //- Temperature limits of applicability of reaction rates
+        scalar Tlow_, Thigh_;
+
         //- Specie info for the left-hand-side of the reaction
         List<specieCoeffs> lhs_;
 
@@ -201,25 +208,13 @@ public:
         );
 
         //- Construct and return a clone
-        virtual autoPtr<Reaction<ReactionThermo>> clone() const //
+        virtual autoPtr<Reaction<ReactionThermo>> clone() const = 0;
-        {
-            return autoPtr<Reaction<ReactionThermo>>
-            (
-                new Reaction<ReactionThermo>(*this)
-            );
-        }
 
         //- Construct and return a clone with new speciesTable
         virtual autoPtr<Reaction<ReactionThermo>> clone
         (
             const speciesTable& species
-        ) const
+        ) const = 0;
-        {
-            return autoPtr<Reaction<ReactionThermo>>
-            (
-                new Reaction<ReactionThermo>(*this, species)
-            );
-        }
 
 
     // Selectors
@@ -242,14 +237,22 @@ public:
 
         // Access
 
+            inline word& name();
             inline const word& name() const;
 
-            // - Access to basic components of the reaction
+            inline scalar Tlow() const;
+            inline scalar Thigh() const;
+
+            //- Access to the components of the left hand side
             inline const List<specieCoeffs>& lhs() const;
+
+            //- Access to the components of the right hand side
             inline const List<specieCoeffs>& rhs() const;
 
-            // - Access to gas components of the reaction
+            //- Access to the gas components of the left hand side
             virtual const List<specieCoeffs>& grhs() const;
+
+            //- Access to the gas components of the right hand side
             virtual const List<specieCoeffs>& glhs() const;
 
             //- Access to specie list
@@ -268,6 +271,49 @@ public:
             );
 
 
+        // Reaction rate coefficients
+
+            //- Forward reaction rate
+            void ddot
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c,
+                scalarField& d
+            ) const;
+
+            //- Backward reaction rate
+            void fdot
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c,
+                scalarField& f
+            ) const;
+
+            //- Net reaction rate for individual species
+            void omega
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c,
+                scalarField& dcdt
+            ) const;
+
+            //- Net reaction rate
+            scalar omega
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c,
+                scalar& pf,
+                scalar& cf,
+                label& lRef,
+                scalar& pr,
+                scalar& cr,
+                label& rRef
+            ) const;
+
         // Reaction rate coefficients
 
             //- Forward rate constant
@@ -276,7 +322,7 @@ public:
                 const scalar p,
                 const scalar T,
                 const scalarField& c
-            ) const;
+            ) const = 0;
 
             //- Reverse rate constant from the given forward rate constant
             virtual scalar kr
@@ -285,18 +331,91 @@ public:
                 const scalar p,
                 const scalar T,
                 const scalarField& c
-            ) const;
+            ) const = 0;
 
-            //- Reverse rate constant.
+            //- Reverse rate constant
-            //  Note this evaluates the forward rate constant and divides by the
-            //  equilibrium constant
             virtual scalar kr
             (
                 const scalar p,
                 const scalar T,
                 const scalarField& c
+            ) const = 0;
+
+
+        // Jacobian coefficients
+
+            //- Derivative of the net reaction rate for each species involved
+            //  w.r.t. the species concentration
+            void dwdc
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c,
+                scalarSquareMatrix& J,
+                scalarField& dcdt,
+                scalar& omegaI,
+                scalar& kfwd,
+                scalar& kbwd,
+                const bool reduced,
+                const List<label>& c2s
             ) const;
 
+            //- Derivative of the net reaction rate for each species involved
+            //  w.r.t. the temperature
+            void dwdT
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c,
+                const scalar omegaI,
+                const scalar kfwd,
+                const scalar kbwd,
+                scalarSquareMatrix& J,
+                const bool reduced,
+                const List<label>& c2s,
+                const label indexT
+            ) const;
+
+            //- Temperature derivative of forward rate
+            virtual scalar dkfdT
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c
+            ) const = 0;
+
+            //- Temperature derivative of reverse rate
+            virtual scalar dkrdT
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c,
+                const scalar dkfdT,
+                const scalar kr
+            ) const = 0;
+
+            //- Third-body efficiencies (beta = 1-alpha)
+            //  non-empty only for third-body reactions
+            //  with enhanced molecularity (alpha != 1)
+            virtual const List<Tuple2<label, scalar>>& beta() const = 0;
+
+            //- Species concentration derivative of the pressure dependent term
+            virtual void dcidc
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c,
+                scalarField& dcidc
+            ) const = 0;
+
+            //- Temperature derivative of the pressure dependent term
+            virtual scalar dcidT
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c
+            ) const = 0;
+
 
         //- Write
         virtual void write(Ostream&) const;

src/thermophysicalModels/specie/reaction/Reactions/Reaction/ReactionI.H

@@ -32,6 +32,13 @@ namespace Foam
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+template<class ReactionThermo>
+inline word& Reaction<ReactionThermo>::name()
+{
+    return name_;
+}
+
+
 template<class ReactionThermo>
 inline const word& Reaction<ReactionThermo>::name() const
 {
@@ -39,6 +46,20 @@ inline const word& Reaction<ReactionThermo>::name() const
 }
 
 
+template<class ReactionThermo>
+inline scalar Reaction<ReactionThermo>::Tlow() const
+{
+    return Tlow_;
+}
+
+
+template<class ReactionThermo>
+inline scalar Reaction<ReactionThermo>::Thigh() const
+{
+    return Thigh_;
+}
+
+
 template<class ReactionThermo>
 inline const List<typename Reaction<ReactionThermo>::specieCoeffs>&
 Reaction<ReactionThermo>::lhs() const

src/thermophysicalModels/specie/reaction/Reactions/ReactionList/ReactionList.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -84,6 +84,13 @@ bool Foam::ReactionList<ThermoType>::readReactionDict()
 {
     const dictionary& reactions(dict_.subDict("reactions"));
 
+    // Set general temperature limits from the dictionary
+    Reaction<ThermoType>::TlowDefault =
+        dict_.lookupOrDefault<scalar>("Tlow", 0);
+
+    Reaction<ThermoType>::ThighDefault =
+        dict_.lookupOrDefault<scalar>("Thigh", great);
+
     forAllConstIter(dictionary, reactions, iter)
     {
         const word reactionName = iter().keyword();

src/thermophysicalModels/specie/reaction/Reactions/ReactionProxy/ReactionProxy.C

@@ -0,0 +1,207 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018 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 "ReactionProxy.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class ReactionThermo>
+Foam::ReactionProxy<ReactionThermo>::ReactionProxy
+(
+    const speciesTable& species,
+    const List<typename Reaction<ReactionThermo>::specieCoeffs>& lhs,
+    const List<typename Reaction<ReactionThermo>::specieCoeffs>& rhs,
+    const HashPtrTable<ReactionThermo>& thermoDatabase
+)
+:
+    Reaction<ReactionThermo>
+    (
+        species,
+        lhs,
+        rhs,
+        thermoDatabase
+    )
+{}
+
+
+template<class ReactionThermo>
+Foam::ReactionProxy<ReactionThermo>::ReactionProxy
+(
+    const Reaction<ReactionThermo>& r,
+    const speciesTable& species
+)
+:
+    Reaction<ReactionThermo>
+    (
+        r,
+        species
+    )
+{}
+
+
+template<class ReactionThermo>
+Foam::ReactionProxy<ReactionThermo>::ReactionProxy
+(
+    const speciesTable& species,
+    const HashPtrTable<ReactionThermo>& thermoDatabase,
+    const dictionary& dict
+)
+:
+    Reaction<ReactionThermo>
+    (
+        species,
+        thermoDatabase,
+        dict
+    )
+{}
+
+
+template<class ReactionThermo>
+Foam::autoPtr<Foam::Reaction<ReactionThermo>>
+Foam::ReactionProxy<ReactionThermo>::clone() const
+{
+    NotImplemented;
+    return autoPtr<Foam::Reaction<ReactionThermo>>();
+}
+
+
+template<class ReactionThermo>
+Foam::autoPtr<Foam::Reaction<ReactionThermo>>
+Foam::ReactionProxy<ReactionThermo>::clone
+(
+    const speciesTable& species
+) const
+{
+    NotImplemented;
+    return autoPtr<Reaction<ReactionThermo>>();
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class ReactionThermo>
+Foam::scalar Foam::ReactionProxy<ReactionThermo>::kf
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    NotImplemented;
+    return 0;
+}
+
+
+template<class ReactionThermo>
+Foam::scalar Foam::ReactionProxy<ReactionThermo>::kr
+(
+    const scalar kfwd,
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    NotImplemented;
+    return 0;
+}
+
+
+template<class ReactionThermo>
+Foam::scalar Foam::ReactionProxy<ReactionThermo>::kr
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    NotImplemented;
+    return 0;
+}
+
+
+template<class ReactionThermo>
+Foam::scalar Foam::ReactionProxy<ReactionThermo>::dkfdT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    NotImplemented;
+    return 0;
+}
+
+
+template<class ReactionThermo>
+Foam::scalar Foam::ReactionProxy<ReactionThermo>::dkrdT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c,
+    const scalar dkfdT,
+    const scalar kr
+) const
+{
+    NotImplemented;
+    return 0;
+}
+
+
+template<class ReactionThermo>
+const Foam::List<Foam::Tuple2<Foam::label, Foam::scalar>>&
+Foam::ReactionProxy<ReactionThermo>::beta() const
+{
+    NotImplemented;
+    return NullObjectRef<List<Tuple2<label, scalar>>>();
+}
+
+
+template<class ReactionThermo>
+void Foam::ReactionProxy<ReactionThermo>::dcidc
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c,
+    scalarField& dcidc
+) const
+{
+    NotImplemented;
+}
+
+
+template<class ReactionThermo>
+Foam::scalar Foam::ReactionProxy<ReactionThermo>::dcidT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    NotImplemented;
+    return 0;
+}
+
+
+// ************************************************************************* //

src/thermophysicalModels/specie/reaction/Reactions/ReactionProxy/ReactionProxy.H

@@ -0,0 +1,189 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018 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::ReactionProxy
+
+Description
+    Proxy version of Reaction which provides dummy implementations of the
+    abstract virtual functions.
+
+    Used for read-construction and format conversion.
+
+SourceFiles
+    ReactionProxy.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef ReactionProxy_H
+#define ReactionProxy_H
+
+#include "Reaction.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                        Class ReactionProxy Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class ReactionThermo>
+class ReactionProxy
+:
+    public Reaction<ReactionThermo>
+{
+
+public:
+
+    // Constructors
+
+        //- Construct from components
+        ReactionProxy
+        (
+            const speciesTable& species,
+            const List<typename Reaction<ReactionThermo>::specieCoeffs>& lhs,
+            const List<typename Reaction<ReactionThermo>::specieCoeffs>& rhs,
+            const HashPtrTable<ReactionThermo>& thermoDatabase
+        );
+
+        //- Construct as copy given new speciesTable
+        ReactionProxy
+        (
+            const Reaction<ReactionThermo>&,
+            const speciesTable& species
+        );
+
+        //- Construct from dictionary
+        ReactionProxy
+        (
+            const speciesTable& species,
+            const HashPtrTable<ReactionThermo>& thermoDatabase,
+            const dictionary& dict
+        );
+
+        //- Construct and return a clone
+        virtual autoPtr<Reaction<ReactionThermo>> clone() const;
+
+        //- Construct and return a clone with new speciesTable
+        virtual autoPtr<Reaction<ReactionThermo>> clone
+        (
+            const speciesTable& species
+        ) const;
+
+
+    //- Destructor
+    virtual ~ReactionProxy()
+    {}
+
+
+    // Member Functions
+
+        // Reaction rate coefficients
+
+            //- Forward rate constant
+            virtual scalar kf
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c
+            ) const;
+
+            //- Reverse rate constant from the given forward rate constant
+            virtual scalar kr
+            (
+                const scalar kfwd,
+                const scalar p,
+                const scalar T,
+                const scalarField& c
+            ) const;
+
+            //- Reverse rate constant
+            virtual scalar kr
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c
+            ) const;
+
+
+        // Jacobian coefficients
+
+            //- Temperature derivative of forward rate
+            virtual scalar dkfdT
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c
+            ) const;
+
+            //- Temperature derivative of reverse rate
+            virtual scalar dkrdT
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c,
+                const scalar dkfdT,
+                const scalar kr
+            ) const;
+
+            //- Third-body efficiencies (beta = 1-alpha)
+            //  non-empty only for third-body reactions
+            //  with enhanced molecularity (alpha != 1)
+            virtual const List<Tuple2<label, scalar>>& beta() const;
+
+            //- Species concentration derivative of the pressure dependent term
+            virtual void dcidc
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c,
+                scalarField& dcidc
+            ) const;
+
+            //- Temperature derivative of the pressure dependent term
+            virtual scalar dcidT
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c
+            ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "ReactionProxy.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/thermophysicalModels/specie/reaction/Reactions/ReversibleReaction/ReversibleReaction.C

@@ -151,6 +151,117 @@ Foam::scalar Foam::ReversibleReaction
 }
 
 
+template
+<
+    template<class> class ReactionType,
+    class ReactionThermo,
+    class ReactionRate
+>
+Foam::scalar Foam::ReversibleReaction
+<
+    ReactionType,
+    ReactionThermo,
+    ReactionRate
+>::dkfdT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    return k_.ddT(p, T, c);
+}
+
+
+template
+<
+    template<class> class ReactionType,
+    class ReactionThermo,
+    class ReactionRate
+>
+Foam::scalar Foam::ReversibleReaction
+<
+    ReactionType,
+    ReactionThermo,
+    ReactionRate
+>::dkrdT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c,
+    const scalar dkfdT,
+    const scalar kr
+) const
+{
+    scalar Kc = max(this->Kc(p, T), rootSmall);
+
+    return dkfdT/Kc - kr*this->dKcdTbyKc(p, T);
+}
+
+
+template
+<
+    template<class> class ReactionType,
+    class ReactionThermo,
+    class ReactionRate
+>
+const Foam::List<Foam::Tuple2<Foam::label, Foam::scalar>>&
+Foam::ReversibleReaction
+<
+    ReactionType,
+    ReactionThermo,
+    ReactionRate
+>::beta() const
+{
+    return k_.beta();
+}
+
+
+template
+<
+    template<class> class ReactionType,
+    class ReactionThermo,
+    class ReactionRate
+>
+void Foam::ReversibleReaction
+<
+    ReactionType,
+    ReactionThermo,
+    ReactionRate
+>::dcidc
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c,
+    scalarField& dcidc
+) const
+{
+    k_.dcidc(p, T, c, dcidc);
+}
+
+
+template
+<
+    template<class> class ReactionType,
+    class ReactionThermo,
+    class ReactionRate
+>
+Foam::scalar Foam::ReversibleReaction
+<
+    ReactionType,
+    ReactionThermo,
+    ReactionRate
+>::dcidT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    return k_.dcidT(p, T, c);
+}
+
+
 template
 <
     template<class> class ReactionType,

src/thermophysicalModels/specie/reaction/Reactions/ReversibleReaction/ReversibleReaction.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -180,6 +180,51 @@ public:
             ) const;
 
 
+        // ReversibleReaction Jacobian functions
+
+            //- Temperature derivative of forward rate
+            virtual scalar dkfdT
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c
+            ) const;
+
+            //- Temperature derivative of backward rate
+            virtual scalar dkrdT
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c,
+                const scalar dkfdT,
+                const scalar kr
+            ) const;
+
+            //- Third-body efficiencies (beta = 1-alpha)
+            //  non-empty only for third-body reactions
+            //  with enhanced molecularity (alpha != 1)
+            virtual const List<Tuple2<label, scalar>>& beta() const;
+
+            //- Species concentration derivative of the pressure dependent term
+            //  By default this value is 1 as it multiplies the third-body term
+            virtual void dcidc
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c,
+                scalarField& dcidc
+            ) const;
+
+            //- Temperature derivative of the pressure dependent term
+            //  By default this value is 0 since ddT of molecularity is approx.0
+            virtual scalar dcidT
+            (
+                const scalar p,
+                const scalar T,
+                const scalarField& c
+            ) const;
+
+
         //- Write
         virtual void write(Ostream&) const;
 };

src/thermophysicalModels/specie/reaction/reactionRate/ArrheniusReactionRate/ArrheniusReactionRate.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -100,6 +100,37 @@ public:
             const scalarField& c
         ) const;
 
+        inline scalar ddT
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c
+        ) const;
+
+        //- Third-body efficiencies (beta = 1-alpha)
+        //  non-empty only for third-body reactions
+        //  with enhanced molecularity (alpha != 1)
+        inline const List<Tuple2<label, scalar>>& beta() const;
+
+        //- Species concentration derivative of the pressure dependent term
+        //  By default this value is 1 as it multiplies the third-body term
+        inline void dcidc
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c,
+            scalarField& dcidc
+        ) const;
+
+        //- Temperature derivative of the pressure dependent term
+        //  By default this value is 0 since ddT of molecularity is approx.0
+        inline scalar dcidT
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c
+        ) const;
+
         //- Write to stream
         inline void write(Ostream& os) const;
 

src/thermophysicalModels/specie/reaction/reactionRate/ArrheniusReactionRate/ArrheniusReactionRateI.H

@@ -75,6 +75,57 @@ inline Foam::scalar Foam::ArrheniusReactionRate::operator()
 }
 
 
+inline Foam::scalar Foam::ArrheniusReactionRate::ddT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField&
+) const
+{
+    scalar ak = A_;
+
+    if (mag(beta_) > vSmall)
+    {
+        ak *= pow(T, beta_);
+    }
+
+    if (mag(Ta_) > vSmall)
+    {
+        ak *= exp(-Ta_/T);
+    }
+
+    return ak*(beta_+Ta_/T)/T;
+}
+
+
+inline const Foam::List<Foam::Tuple2<Foam::label, Foam::scalar>>&
+Foam::ArrheniusReactionRate::beta() const
+{
+    return NullObjectRef<List<Tuple2<label, scalar>>>();
+}
+
+
+inline void Foam::ArrheniusReactionRate::dcidc
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c,
+    scalarField& dcidc
+) const
+{}
+
+
+inline Foam::scalar Foam::ArrheniusReactionRate::dcidT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    return 0;
+}
+
+
 inline void Foam::ArrheniusReactionRate::write(Ostream& os) const
 {
     os.writeKeyword("A") << A_ << token::END_STATEMENT << nl;

src/thermophysicalModels/specie/reaction/reactionRate/ChemicallyActivatedReactionRate/ChemicallyActivatedReactionRate.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -69,6 +69,7 @@ class ChemicallyActivatedReactionRate
         ReactionRate kInf_;
         ChemicallyActivationFunction F_;
         thirdBodyEfficiencies thirdBodyEfficiencies_;
+        List<Tuple2<label, scalar>> beta_;
 
 
 public:
@@ -109,6 +110,40 @@ public:
             const scalarField& c
         ) const;
 
+        inline scalar ddT
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c
+        ) const;
+
+        //- Third-body efficiencies (beta = 1-alpha)
+        //  non-empty only for third-body reactions
+        //  with enhanced molecularity (alpha != 1)
+        inline const List<Tuple2<label, scalar>>& beta() const
+        {
+            return beta_;
+        }
+
+        //- Species concentration derivative of the pressure dependent term
+        //  By default this value is 1 as it multiplies the third-body term
+        inline void dcidc
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c,
+            scalarField& dcidc
+        ) const;
+
+        //- Temperature derivative of the pressure dependent term
+        //  By default this value is 0 since ddT of molecularity is approx.0
+        inline scalar dcidT
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c
+        ) const;
+
         //- Write to stream
         inline void write(Ostream& os) const;
 

src/thermophysicalModels/specie/reaction/reactionRate/ChemicallyActivatedReactionRate/ChemicallyActivatedReactionRateI.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -42,7 +42,12 @@ inline Foam::ChemicallyActivatedReactionRate
     kInf_(kInf),
     F_(F),
     thirdBodyEfficiencies_(tbes)
-{}
+{
+    forAll(tbes, i)
+    {
+        beta_.append(Tuple2<label, scalar>(i, tbes[i]));
+    }
+}
 
 
 template<class ReactionRate, class ChemicallyActivationFunction>
@@ -60,7 +65,19 @@ inline Foam::ChemicallyActivatedReactionRate
     kInf_(species, dict),
     F_(dict),
     thirdBodyEfficiencies_(species, dict)
-{}
+{
+    forAll(thirdBodyEfficiencies_, i)
+    {
+        beta_.append
+        (
+            Tuple2<label, scalar>
+            (
+                i,
+                thirdBodyEfficiencies_[i]
+            )
+        );
+    }
+}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
@@ -77,15 +94,108 @@ inline Foam::scalar Foam::ChemicallyActivatedReactionRate
     const scalarField& c
 ) const
 {
-    scalar k0 = k0_(p, T, c);
+    const scalar k0 = k0_(p, T, c);
-    scalar kInf = kInf_(p, T, c);
+    const scalar kInf = kInf_(p, T, c);
-
+    const scalar Pr = k0*thirdBodyEfficiencies_.M(c)/kInf;
-    scalar Pr = k0*thirdBodyEfficiencies_.M(c)/kInf;
 
     return k0*(1/(1 + Pr))*F_(T, Pr);
 }
 
 
+template<class ReactionRate, class ChemicallyActivationFunction>
+inline Foam::scalar Foam::ChemicallyActivatedReactionRate
+<
+    ReactionRate,
+    ChemicallyActivationFunction
+>::ddT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    const scalar k0 = k0_(p, T, c);
+    const scalar kInf = kInf_(p, T, c);
+    const scalar Pr = k0*thirdBodyEfficiencies_.M(c)/kInf;
+
+    return (1/(1 + Pr))*F_(T, Pr)*k0_.ddT(p, T, c);
+}
+
+
+template<class ReactionRate, class ChemicallyActivationFunction>
+inline void Foam::ChemicallyActivatedReactionRate
+<
+    ReactionRate,
+    ChemicallyActivationFunction
+>::dcidc
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c,
+    scalarField& dcidc
+) const
+{
+    const scalar M = thirdBodyEfficiencies_.M(c);
+
+    if (M > small)
+    {
+        const scalar k0 = k0_(p, T, c);
+        const scalar kInf = kInf_(p, T, c);
+
+        const scalar Pr = k0*M/kInf;
+        const scalar F = F_(T, Pr);
+
+        forAll(beta_, i)
+        {
+            const scalar dPrdci = -beta_[i].second()*k0/kInf;
+            const scalar dFdci = F_.ddc(Pr, F, dPrdci, T);
+            dcidc[i] = (-dPrdci/(1 + Pr) + dFdci/F);
+        }
+    }
+    else
+    {
+        forAll(beta_, i)
+        {
+            dcidc[i] = 0;
+        }
+    }
+}
+
+
+template<class ReactionRate, class ChemicallyActivationFunction>
+inline Foam::scalar Foam::ChemicallyActivatedReactionRate
+<
+    ReactionRate,
+    ChemicallyActivationFunction
+>::dcidT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    const scalar M = thirdBodyEfficiencies_.M(c);
+
+    if (M > small)
+    {
+        const scalar k0 = k0_(p, T, c);
+        const scalar kInf = kInf_(p, T, c);
+
+        const scalar Pr = k0*thirdBodyEfficiencies_.M(c)/kInf;
+        const scalar F = F_(T, Pr);
+        const scalar dPrdT =
+            Pr*(k0_.ddT(p, T, c)/k0 - kInf_.ddT(p, T, c)/kInf - 1/T);
+        const scalar dFdT = F_.ddT(Pr, F, dPrdT, T);
+
+        return (-dPrdT/(1 + Pr) + dFdT/F);
+    }
+    else
+    {
+        return 0;
+    }
+}
+
+
 template<class ReactionRate, class ChemicallyActivationFunction>
 inline void Foam::ChemicallyActivatedReactionRate
 <

src/thermophysicalModels/specie/reaction/reactionRate/FallOffReactionRate/FallOffReactionRate.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -68,6 +68,7 @@ class FallOffReactionRate
         ReactionRate kInf_;
         FallOffFunction F_;
         thirdBodyEfficiencies thirdBodyEfficiencies_;
+        List<Tuple2<label, scalar>> beta_;
 
 
 public:
@@ -106,6 +107,38 @@ public:
             const scalarField& c
         ) const;
 
+        inline scalar ddT
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c
+        ) const;
+
+        //- Third-body efficiencies (beta = 1-alpha)
+        //  non-empty only for third-body reactions
+        //  with enhanced molecularity (alpha != 1)
+        inline const List<Tuple2<label, scalar>>& beta() const
+        {
+            return beta_;
+        }
+
+        //- Species concentration derivative of the pressure dependent term
+        inline void dcidc
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c,
+            scalarField& dcidc
+        ) const;
+
+        //- Temperature derivative of the pressure dependent term
+        inline scalar dcidT
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c
+        ) const;
+
         //- Write to stream
         inline void write(Ostream& os) const;
 

src/thermophysicalModels/specie/reaction/reactionRate/FallOffReactionRate/FallOffReactionRateI.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -39,7 +39,12 @@ FallOffReactionRate
     kInf_(kInf),
     F_(F),
     thirdBodyEfficiencies_(tbes)
-{}
+{
+    forAll(tbes, i)
+    {
+        beta_.append(Tuple2<label, scalar>(i, tbes[i]));
+    }
+}
 
 
 template<class ReactionRate, class FallOffFunction>
@@ -54,7 +59,19 @@ FallOffReactionRate
     kInf_(species, dict.subDict("kInf")),
     F_(dict.subDict("F")),
     thirdBodyEfficiencies_(species, dict.subDict("thirdBodyEfficiencies"))
-{}
+{
+    forAll(thirdBodyEfficiencies_, i)
+    {
+        beta_.append
+        (
+            Tuple2<label, scalar>
+            (
+                i,
+                thirdBodyEfficiencies_[i]
+            )
+        );
+    }
+}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
@@ -68,15 +85,97 @@ Foam::FallOffReactionRate<ReactionRate, FallOffFunction>::operator()
     const scalarField& c
 ) const
 {
-    scalar k0 = k0_(p, T, c);
+    const scalar k0 = k0_(p, T, c);
-    scalar kInf = kInf_(p, T, c);
+    const scalar kInf = kInf_(p, T, c);
-
+    const scalar Pr = k0*thirdBodyEfficiencies_.M(c)/kInf;
-    scalar Pr = k0*thirdBodyEfficiencies_.M(c)/kInf;
 
     return kInf*(Pr/(1 + Pr))*F_(T, Pr);
 }
 
 
+template<class ReactionRate, class FallOffFunction>
+inline Foam::scalar
+Foam::FallOffReactionRate<ReactionRate, FallOffFunction>::ddT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    const scalar k0 = k0_(p, T, c);
+    const scalar kInf = kInf_(p, T, c);
+    const scalar Pr = k0*thirdBodyEfficiencies_.M(c)/kInf;
+
+    return (Pr/(1 + Pr))*F_(T, Pr)*kInf_.ddT(p, T, c);
+}
+
+
+template<class ReactionRate, class FallOffFunction>
+inline void Foam::FallOffReactionRate<ReactionRate, FallOffFunction>::dcidc
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c,
+    scalarField& dcidc
+) const
+{
+    const scalar M = thirdBodyEfficiencies_.M(c);
+
+    if (M > small)
+    {
+        const scalar k0 = k0_(p, T, c);
+        const scalar kInf = kInf_(p, T, c);
+        const scalar Pr = k0*M/kInf;
+        const scalar F = F_(T, Pr);
+
+        forAll(beta_, i)
+        {
+            const scalar dPrdci = -beta_[i].second()*k0/kInf;
+            const scalar dFdci = F_.ddc(Pr, F, dPrdci, T);
+            dcidc[i] = (dPrdci/(Pr*(1 + Pr)) + dFdci/F);
+        }
+    }
+    else
+    {
+        forAll(beta_, i)
+        {
+            dcidc[i] = 0;
+        }
+    }
+}
+
+
+template<class ReactionRate, class FallOffFunction>
+inline Foam::scalar
+Foam::FallOffReactionRate<ReactionRate, FallOffFunction>::dcidT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    const scalar M = thirdBodyEfficiencies_.M(c);
+
+    if (M > small)
+    {
+        const scalar k0 = k0_(p, T, c);
+        const scalar kInf = kInf_(p, T, c);
+
+        const scalar Pr = k0*thirdBodyEfficiencies_.M(c)/kInf;
+        const scalar F = F_(T, Pr);
+        const scalar dPrdT =
+            Pr*(k0_.ddT(p, T, c)/k0 - kInf_.ddT(p, T, c)/kInf - 1/T);
+        const scalar dFdT = F_.ddT(Pr, F, dPrdT, T);
+
+        return (dPrdT/(Pr*(1 + Pr)) + dFdT/F);
+    }
+    else
+    {
+        return 0;
+    }
+}
+
+
 template<class ReactionRate, class FallOffFunction>
 inline void Foam::FallOffReactionRate<ReactionRate, FallOffFunction>::write
 (

src/thermophysicalModels/specie/reaction/reactionRate/JanevReactionRate/JanevReactionRate.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -103,6 +103,35 @@ public:
             const scalarField& c
         ) const;
 
+        inline scalar ddT
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c
+        ) const;
+
+        //- Third-body efficiencies (beta = 1-alpha)
+        //  non-empty only for third-body reactions
+        //  with enhanced molecularity (alpha != 1)
+        inline const List<Tuple2<label, scalar>>& beta() const;
+
+        //- Species concentration derivative of the pressure dependent term
+        inline void dcidc
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c,
+            scalarField& dcidc
+        ) const;
+
+        //- Temperature derivative of the pressure dependent term
+        inline scalar dcidT
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c
+        ) const;
+
         //- Write to stream
         inline void write(Ostream& os) const;
 

src/thermophysicalModels/specie/reaction/reactionRate/JanevReactionRate/JanevReactionRateI.H

@@ -69,7 +69,7 @@ inline Foam::scalar Foam::JanevReactionRate::operator()
         lta *= pow(T, beta_);
     }
 
-    scalar expArg = 0.0;
+    scalar expArg = 0;
 
     if (mag(Ta_) > vSmall)
     {
@@ -89,6 +89,75 @@ inline Foam::scalar Foam::JanevReactionRate::operator()
 }
 
 
+inline Foam::scalar Foam::JanevReactionRate::ddT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField&
+) const
+{
+    scalar lta = A_;
+
+    if (mag(beta_) > vSmall)
+    {
+        lta *= pow(T, beta_);
+    }
+
+    scalar expArg = 0;
+
+    if (mag(Ta_) > vSmall)
+    {
+        expArg -= Ta_/T;
+    }
+
+    scalar lnT = log(T);
+
+    for (int n=0; n<nb_; n++)
+    {
+        expArg += b_[n]*pow(lnT, n);
+    }
+
+    scalar deriv = b_[1];
+
+    for (int n=2; n<nb_; n++)
+    {
+        deriv += n*b_[n]*pow(lnT, n-1);
+    }
+
+    lta *= exp(expArg);
+
+    return lta*(beta_+Ta_/T+deriv)/T;
+}
+
+
+inline const Foam::List<Foam::Tuple2<Foam::label, Foam::scalar>>&
+Foam::JanevReactionRate::beta() const
+{
+    return NullObjectRef<List<Tuple2<label, scalar>>>();
+}
+
+
+inline void Foam::JanevReactionRate::dcidc
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c,
+    scalarField& dcidc
+) const
+{}
+
+
+inline Foam::scalar Foam::JanevReactionRate::dcidT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    return 0;
+}
+
+
 inline void Foam::JanevReactionRate::write(Ostream& os) const
 {
     os.writeKeyword("A") << A_ << nl;

src/thermophysicalModels/specie/reaction/reactionRate/LandauTellerReactionRate/LandauTellerReactionRate.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -102,6 +102,35 @@ public:
             const scalarField& c
         ) const;
 
+        inline scalar ddT
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c
+        ) const;
+
+        //- Third-body efficiencies (beta = 1-alpha)
+        //  non-empty only for third-body reactions
+        //  with enhanced molecularity (alpha != 1)
+        inline const List<Tuple2<label, scalar>>& beta() const;
+
+        //- Species concentration derivative of the pressure dependent term
+        inline void dcidc
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c,
+            scalarField& dcidc
+        ) const;
+
+        //- Temperature derivative of the pressure dependent term
+        inline scalar dcidT
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c
+        ) const;
+
         //- Write to stream
         inline void write(Ostream& os) const;
 

src/thermophysicalModels/specie/reaction/reactionRate/LandauTellerReactionRate/LandauTellerReactionRateI.H

@@ -72,7 +72,7 @@ inline Foam::scalar Foam::LandauTellerReactionRate::operator()
         lta *= pow(T, beta_);
     }
 
-    scalar expArg = 0.0;
+    scalar expArg = 0;
 
     if (mag(Ta_) > vSmall)
     {
@@ -98,6 +98,81 @@ inline Foam::scalar Foam::LandauTellerReactionRate::operator()
 }
 
 
+inline Foam::scalar Foam::LandauTellerReactionRate::ddT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField&
+) const
+{
+    scalar lta = A_;
+
+    if (mag(beta_) > vSmall)
+    {
+        lta *= pow(T, beta_);
+    }
+
+    scalar expArg = 0;
+    scalar deriv = 0;
+
+    if (mag(Ta_) > vSmall)
+    {
+        scalar TaT = Ta_/T;
+        expArg -= TaT;
+        deriv += TaT;
+    }
+
+    if (mag(B_) > vSmall)
+    {
+        scalar BT = B_/cbrt(T);
+        expArg += BT;
+        deriv -= BT/3;
+    }
+
+    if (mag(C_) > vSmall)
+    {
+        scalar CT = C_/pow(T, 2.0/3.0);
+        expArg += CT;
+        deriv -= 2*CT/3;
+    }
+
+    if (mag(expArg) > vSmall)
+    {
+        lta *= exp(expArg);
+    }
+
+    return lta*(beta_+deriv)/T;
+}
+
+
+inline const Foam::List<Foam::Tuple2<Foam::label, Foam::scalar>>&
+Foam::LandauTellerReactionRate::beta() const
+{
+    return NullObjectRef<List<Tuple2<label, scalar>>>();
+}
+
+
+inline void Foam::LandauTellerReactionRate::dcidc
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c,
+    scalarField& dcidc
+) const
+{}
+
+
+inline Foam::scalar Foam::LandauTellerReactionRate::dcidT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    return 0;
+}
+
+
 inline void Foam::LandauTellerReactionRate::write(Ostream& os) const
 {
     os.writeKeyword("A") << A_ << token::END_STATEMENT << nl;

src/thermophysicalModels/specie/reaction/reactionRate/LangmuirHinshelwood/LangmuirHinshelwoodReactionRate.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -104,6 +104,35 @@ public:
             const scalarField& c
         ) const;
 
+        inline scalar ddT
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c
+        ) const;
+
+        //- Third-body efficiencies (beta = 1-alpha)
+        //  non-empty only for third-body reactions
+        //  with enhanced molecularity (alpha != 1)
+        inline const List<Tuple2<label, scalar>>& beta() const;
+
+        //- Species concentration derivative of the pressure dependent term
+        inline void dcidc
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c,
+            scalarField& dcidc
+        ) const;
+
+        //- Temperature derivative of the pressure dependent term
+        inline scalar dcidT
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c
+        ) const;
+
         //- Write to stream
         inline void write(Ostream& os) const;
 

src/thermophysicalModels/specie/reaction/reactionRate/LangmuirHinshelwood/LangmuirHinshelwoodReactionRateI.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -89,6 +89,77 @@ inline Foam::scalar Foam::LangmuirHinshelwoodReactionRate::operator()
 }
 
 
+inline Foam::scalar Foam::LangmuirHinshelwoodReactionRate::ddT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    scalar den =
+    (
+        T
+       *sqr(1 + A_[1]*exp(-Ta_[1]/T)*c[co_] + A_[2]*exp(-Ta_[2]/T)*c[c3h6_])
+       *(1 + A_[3]*exp(-Ta_[3]/T)*sqr(c[co_])*sqr(c[c3h6_]))
+       *(1 + A_[4]*exp(-Ta_[4]/T)*pow(c[no_], 0.7))
+    );
+    scalar rate = A_[0]*exp(-Ta_[0]/T)/ den;
+
+    scalar derivDen =
+    (
+        sqr(1 + A_[1]*exp(-Ta_[1]/T)*c[co_] + A_[2]*exp(-Ta_[2]/T)*c[c3h6_])
+       *(1 + A_[3]*exp(-Ta_[3]/T)*sqr(c[co_])*sqr(c[c3h6_]))
+       *(1 + A_[4]*exp(-Ta_[4]/T)*pow(c[no_], 0.7))
+      + 2*T*(1 + A_[1]*exp(-Ta_[1]/T)*c[co_] + A_[2]*exp(-Ta_[2]/T)*c[c3h6_])
+       *
+        (
+            A_[1]*exp(-Ta_[1]/T)*c[co_]*Ta_[1]/sqr(T)
+          + A_[2]*exp(-Ta_[2]/T)*c[c3h6_]*Ta_[2]/sqr(T)
+        )
+       *(1 + A_[3]*exp(-Ta_[3]/T)*sqr(c[co_])*sqr(c[c3h6_]))
+       *(1 + A_[4]*exp(-Ta_[4]/T)*pow(c[no_], 0.7))
+      + T
+       *sqr(1 + A_[1]*exp(-Ta_[1]/T)*c[co_] + A_[2]*exp(-Ta_[2]/T)*c[c3h6_])
+       *(A_[3]*exp(-Ta_[3]/T)*Ta_[3]*sqr(c[co_])*sqr(c[c3h6_])/sqr(T))
+       *(1 + A_[4]*exp(-Ta_[4]/T)*pow(c[no_], 0.7))
+      + T
+       *sqr(1 + A_[1]*exp(-Ta_[1]/T)*c[co_] + A_[2]*exp(-Ta_[2]/T)*c[c3h6_])
+       *(1 + A_[3]*exp(-Ta_[3]/T)*sqr(c[co_])*sqr(c[c3h6_]))
+       *(A_[4]*exp(-Ta_[4]/T)*Ta_[4]*pow(c[no_], 0.7))/sqr(T)
+    );
+
+    return rate*(Ta_[0]/sqr(T) - derivDen/den);
+}
+
+
+inline const Foam::List<Foam::Tuple2<Foam::label, Foam::scalar>>&
+Foam::LangmuirHinshelwoodReactionRate::beta() const
+{
+    return NullObjectRef<List<Tuple2<label, scalar>>>();
+}
+
+
+inline void Foam::LangmuirHinshelwoodReactionRate::dcidc
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c,
+    scalarField& dcidc
+) const
+{}
+
+
+inline Foam::scalar Foam::LangmuirHinshelwoodReactionRate::dcidT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    return 0;
+}
+
+
 inline void Foam::LangmuirHinshelwoodReactionRate::write(Ostream& os) const
 {
     FixedList<Tuple2<scalar, scalar>, n_> coeffs;

src/thermophysicalModels/specie/reaction/reactionRate/fallOffFunctions/LindemannFallOffFunction/LindemannFallOffFunction.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -80,6 +80,22 @@ public:
             const scalar Pr
         ) const;
 
+        inline scalar ddT
+        (
+            const scalar Pr,
+            const scalar F,
+            const scalar dPrdT,
+            const scalar T
+        ) const;
+
+        inline scalar ddc
+        (
+            const scalar Pr,
+            const scalar F,
+            const scalar dPrdc,
+            const scalar T
+        ) const;
+
         //- Write to stream
         inline void write(Ostream& os) const;
 

src/thermophysicalModels/specie/reaction/reactionRate/fallOffFunctions/LindemannFallOffFunction/LindemannFallOffFunctionI.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -48,6 +48,30 @@ inline Foam::scalar Foam::LindemannFallOffFunction::operator()
 }
 
 
+inline Foam::scalar Foam::LindemannFallOffFunction::ddT
+(
+    const scalar,
+    const scalar,
+    const scalar,
+    const scalar
+) const
+{
+    return 0;
+}
+
+
+inline Foam::scalar Foam::LindemannFallOffFunction::ddc
+(
+    const scalar,
+    const scalar,
+    const scalar,
+    const scalar
+) const
+{
+    return 0;
+}
+
+
 inline void Foam::LindemannFallOffFunction::write(Ostream& os) const
 {}
 

src/thermophysicalModels/specie/reaction/reactionRate/fallOffFunctions/SRIFallOffFunction/SRIFallOffFunction.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -91,6 +91,22 @@ public:
             const scalar Pr
         ) const;
 
+        inline scalar ddT
+        (
+            const scalar Pr,
+            const scalar F,
+            const scalar dPrdT,
+            const scalar T
+        ) const;
+
+        inline scalar ddc
+        (
+            const scalar Pr,
+            const scalar F,
+            const scalar dPrdc,
+            const scalar T
+        ) const;
+
         //- Write to stream
         inline void write(Ostream& os) const;
 

src/thermophysicalModels/specie/reaction/reactionRate/fallOffFunctions/SRIFallOffFunction/SRIFallOffFunctionI.H

@@ -60,11 +60,48 @@ inline Foam::scalar Foam::SRIFallOffFunction::operator()
     const scalar Pr
 ) const
 {
-    scalar X = 1.0/(1.0 + sqr(log10(max(Pr, small))));
+    scalar X = 1.0/(1 + sqr(log10(max(Pr, small))));
     return d_*pow(a_*exp(-b_/T) + exp(-T/c_), X)*pow(T, e_);
 }
 
 
+inline Foam::scalar Foam::SRIFallOffFunction::ddT
+(
+    const scalar Pr,
+    const scalar F,
+    const scalar dPrdT,
+    const scalar T
+) const
+{
+    scalar X = 1.0/(1 + sqr(log10(max(Pr, small))));
+    scalar dXdPr = -X*X*2*log10(Pr)/Pr/log(10.0);
+    return
+        F
+       *(
+            e_/T
+          + X
+           *(a_*b_*exp(-b_/T)/sqr(T) - exp(-T/c_)/c_)
+           /(a_*exp(-b_/T) + exp(-T/c_))
+          + dXdPr*dPrdT*log(a_*exp(-b_/T) + exp(-T/c_))
+        );
+}
+
+
+inline Foam::scalar Foam::SRIFallOffFunction::ddc
+(
+    const scalar Pr,
+    const scalar F,
+    const scalar dPrdc,
+    const scalar T
+) const
+{
+    scalar X = 1.0/(1 + sqr(log10(max(Pr, small))));
+    scalar dXdPr = -X*X*2*log10(Pr)/Pr/log(10.0);
+    scalar dXdc = dXdPr*dPrdc;
+    return F*dXdc*log(a_*exp(-b_/T) + exp(-T/c_));
+}
+
+
 inline void Foam::SRIFallOffFunction::write(Ostream& os) const
 {
     os.writeKeyword("a") << a_ << token::END_STATEMENT << nl;

src/thermophysicalModels/specie/reaction/reactionRate/fallOffFunctions/TroeFallOffFunction/TroeFallOffFunction.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -91,6 +91,22 @@ public:
             const scalar Pr
         ) const;
 
+        inline scalar ddT
+        (
+            const scalar Pr,
+            const scalar F,
+            const scalar dPrdT,
+            const scalar T
+        ) const;
+
+        inline scalar ddc
+        (
+            const scalar Pr,
+            const scalar F,
+            const scalar dPrdc,
+            const scalar T
+        ) const;
+
         //- Write to stream
         inline void write(Ostream& os) const;
 

src/thermophysicalModels/specie/reaction/reactionRate/fallOffFunctions/TroeFallOffFunction/TroeFallOffFunctionI.H

@@ -71,7 +71,102 @@ inline Foam::scalar Foam::TroeFallOffFunction::operator()
     scalar n = 0.75 - 1.27*logFcent;
 
     scalar logPr = log10(max(Pr, small));
-    return pow(10.0, logFcent/(1.0 + sqr((logPr + c)/(n - d*(logPr + c)))));
+    return pow(10, logFcent/(1 + sqr((logPr + c)/(n - d*(logPr + c)))));
+}
+
+
+inline Foam::scalar Foam::TroeFallOffFunction::ddT
+(
+    const scalar Pr,
+    const scalar F,
+    const scalar dPrdT,
+    const scalar T
+) const
+{
+    scalar logPr = log10(max(Pr, small));
+    scalar logTen = log(10.0);
+    scalar Fcent =
+    (
+        max
+        (
+            (1 - alpha_)*exp(-T/Tsss_) + alpha_*exp(-T/Ts_) + exp(-Tss_/T),
+            small
+        )
+    );
+    scalar logFcent = log10(Fcent);
+
+    scalar dFcentdT =
+    (
+        (alpha_ - 1)*exp(-T/Tsss_)/Tsss_
+      - alpha_*exp(-T/Ts_)/Ts_
+      + Tss_*exp(-Tss_/T)/sqr(T)
+    );
+
+    scalar d = 0.14;
+    scalar dlogFcentdT = dFcentdT/Fcent/logTen;
+    scalar c = -0.4 - 0.67*logFcent;
+    scalar dcdT = -0.67*dlogFcentdT;
+    scalar n = 0.75 - 1.27*logFcent;
+    scalar dndT = -1.27*dlogFcentdT;
+
+    scalar dlogPrdT = dPrdT/Pr/logTen;
+
+    scalar dParentdT =
+        2.0*(logPr + c)/sqr(n - d*(logPr + c))
+       *(
+            (dlogPrdT + dcdT)
+          - (logPr + c)*(dndT - d*(dlogPrdT + dcdT))/(n - d*(logPr + c))
+        );
+
+    return
+    (
+        F*logTen
+       *(
+            dlogFcentdT/(1.0 + sqr((logPr + c)/(n - d*(logPr + c))))
+          - logFcent*dParentdT/sqr(1.0 + sqr((logPr + c)/(n - d*(logPr + c))))
+        )
+    );
+}
+
+
+inline Foam::scalar Foam::TroeFallOffFunction::ddc
+(
+    const scalar Pr,
+    const scalar F,
+    const scalar dPrdc,
+    const scalar T
+) const
+{
+    scalar logPr = log10(max(Pr, small));
+    scalar logTen = log(10.0);
+    scalar logFcent = log10
+    (
+        max
+        (
+            (1 - alpha_)*exp(-T/Tsss_) + alpha_*exp(-T/Ts_) + exp(-Tss_/T),
+            small
+        )
+    );
+
+    scalar dlogPrdc = dPrdc/Pr/logTen;
+    scalar d = 0.14;
+    scalar c = -0.4 - 0.67*logFcent;
+    scalar n = 0.75 - 1.27*logFcent;
+
+    scalar dParentdc =
+        2.0*(logPr + c)/sqr(n - d*(logPr + c))
+       *(
+            (dlogPrdc)
+          - (logPr + c)*(-d*(dlogPrdc))/(n - d*(logPr + c))
+        );
+
+    return
+    (
+        F*logTen
+       *(
+          - logFcent*dParentdc/sqr(1.0 + sqr((logPr + c)/(n - d*(logPr + c))))
+        )
+    );
 }
 
 

src/thermophysicalModels/specie/reaction/reactionRate/infiniteReactionRate/infiniteReactionRate.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -87,6 +87,35 @@ public:
             const scalarField& c
         ) const;
 
+        inline scalar ddT
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c
+        ) const;
+
+        //- Third-body efficiencies (beta = 1-alpha)
+        //  non-empty only for third-body reactions
+        //  with enhanced molecularity (alpha != 1)
+        inline const List<Tuple2<label, scalar>>& beta() const;
+
+        //- Species concentration derivative of the pressure dependent term
+        inline void dcidc
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c,
+            scalarField& dcidc
+        ) const;
+
+        //- Temperature derivative of the pressure dependent term
+        inline scalar dcidT
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c
+        ) const;
+
         //- Write to stream
         inline void write(Ostream& os) const;
 

src/thermophysicalModels/specie/reaction/reactionRate/infiniteReactionRate/infiniteReactionRateI.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -53,6 +53,44 @@ inline Foam::scalar Foam::infiniteReactionRate::operator()
     return (1);
 }
 
+inline Foam::scalar Foam::infiniteReactionRate::ddT
+(
+    const scalar p,
+    const scalar,
+    const scalarField&
+) const
+{
+    return (0);
+}
+
+
+inline const Foam::List<Foam::Tuple2<Foam::label, Foam::scalar>>&
+Foam::infiniteReactionRate::beta() const
+{
+    return NullObjectRef<List<Tuple2<label, scalar>>>();
+}
+
+
+inline void Foam::infiniteReactionRate::dcidc
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c,
+    scalarField& dcidc
+) const
+{}
+
+
+inline Foam::scalar Foam::infiniteReactionRate::dcidT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    return 0;
+}
+
 
 inline Foam::Ostream& Foam::operator<<
 (

src/thermophysicalModels/specie/reaction/reactionRate/powerSeries/powerSeriesReactionRate.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -103,6 +103,35 @@ public:
             const scalarField& c
         ) const;
 
+        inline scalar ddT
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c
+        ) const;
+
+        //- Third-body efficiencies (beta = 1-alpha)
+        //  non-empty only for third-body reactions
+        //  with enhanced molecularity (alpha != 1)
+        inline const List<Tuple2<label, scalar>>& beta() const;
+
+        //- Species concentration derivative of the pressure dependent term
+        inline void dcidc
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c,
+            scalarField& dcidc
+        ) const;
+
+        //- Temperature derivative of the pressure dependent term
+        inline scalar dcidT
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c
+        ) const;
+
         //- Write to stream
         inline void write(Ostream& os) const;
 

src/thermophysicalModels/specie/reaction/reactionRate/powerSeries/powerSeriesReactionRateI.H

@@ -69,7 +69,7 @@ inline Foam::scalar Foam::powerSeriesReactionRate::operator()
         lta *= pow(T, beta_);
     }
 
-    scalar expArg = 0.0;
+    scalar expArg = 0;
 
     forAll(coeffs_, n)
     {
@@ -82,6 +82,64 @@ inline Foam::scalar Foam::powerSeriesReactionRate::operator()
 }
 
 
+inline Foam::scalar Foam::powerSeriesReactionRate::ddT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField&
+) const
+{
+    scalar lta = A_;
+
+    if (mag(beta_) > vSmall)
+    {
+        lta *= pow(T, beta_);
+    }
+
+    scalar expArg = 0;
+    scalar deriv = 0;
+
+    forAll(coeffs_, n)
+    {
+        scalar cT = coeffs_[n]/pow(T, n + 1);
+        expArg += cT;
+        deriv -= (n + 1)*cT;
+    }
+
+    lta *= exp(expArg);
+
+    return lta*(beta_+deriv)/T;
+}
+
+
+inline const Foam::List<Foam::Tuple2<Foam::label, Foam::scalar>>&
+Foam::powerSeriesReactionRate::beta() const
+{
+    return NullObjectRef<List<Tuple2<label, scalar>>>();
+}
+
+
+inline void Foam::powerSeriesReactionRate::dcidc
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c,
+    scalarField& dcidc
+) const
+{}
+
+
+inline Foam::scalar Foam::powerSeriesReactionRate::dcidT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    return 0;
+}
+
+
 inline void Foam::powerSeriesReactionRate::write(Ostream& os) const
 {
     os.writeKeyword("A") << A_ << token::END_STATEMENT << nl;

src/thermophysicalModels/specie/reaction/reactionRate/thirdBodyArrheniusReactionRate/thirdBodyArrheniusReactionRate.H

@@ -61,6 +61,7 @@ class thirdBodyArrheniusReactionRate
     // Private data
 
         thirdBodyEfficiencies thirdBodyEfficiencies_;
+        List<Tuple2<label, scalar>> beta_;
 
 
 public:
@@ -99,6 +100,38 @@ public:
             const scalarField& c
         ) const;
 
+        inline scalar ddT
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c
+        ) const;
+
+        //- Third-body efficiencies (beta = 1-alpha)
+        //  non-empty only for third-body reactions
+        //  with enhanced molecularity (alpha != 1)
+        inline const List<Tuple2<label, scalar>>& beta() const
+        {
+            return beta_;
+        }
+
+        //- Species concentration derivative of the pressure dependent term
+        inline void dcidc
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c,
+            scalarField& dcidc
+        ) const;
+
+        //- Temperature derivative of the pressure dependent term
+        inline scalar dcidT
+        (
+            const scalar p,
+            const scalar T,
+            const scalarField& c
+        ) const;
+
         //- Write to stream
         inline void write(Ostream& os) const;
 

src/thermophysicalModels/specie/reaction/reactionRate/thirdBodyArrheniusReactionRate/thirdBodyArrheniusReactionRateI.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -35,7 +35,12 @@ inline Foam::thirdBodyArrheniusReactionRate::thirdBodyArrheniusReactionRate
 :
     ArrheniusReactionRate(A, beta, Ta),
     thirdBodyEfficiencies_(tbes)
-{}
+{
+    forAll(tbes, i)
+    {
+        beta_.append(Tuple2<label, scalar>(i, tbes[i]));
+    }
+}
 
 
 inline Foam::thirdBodyArrheniusReactionRate::thirdBodyArrheniusReactionRate
@@ -50,7 +55,19 @@ inline Foam::thirdBodyArrheniusReactionRate::thirdBodyArrheniusReactionRate
         dict
     ),
     thirdBodyEfficiencies_(species, dict)
-{}
+{
+    forAll(thirdBodyEfficiencies_, i)
+    {
+        beta_.append
+        (
+            Tuple2<label, scalar>
+            (
+                i,
+                thirdBodyEfficiencies_[i]
+            )
+        );
+    }
+}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
@@ -68,6 +85,46 @@ inline Foam::scalar Foam::thirdBodyArrheniusReactionRate::operator()
 }
 
 
+inline Foam::scalar Foam::thirdBodyArrheniusReactionRate::ddT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    return
+        thirdBodyEfficiencies_.M(c)
+       *ArrheniusReactionRate::ddT(p, T, c);
+}
+
+
+inline void Foam::thirdBodyArrheniusReactionRate::dcidc
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c,
+    scalarField& dcidc
+) const
+{
+    scalar M = thirdBodyEfficiencies_.M(c);
+    forAll(beta_, i)
+    {
+        dcidc[i] = beta_[i].second()/max(M, small);
+    }
+}
+
+
+inline Foam::scalar Foam::thirdBodyArrheniusReactionRate::dcidT
+(
+    const scalar p,
+    const scalar T,
+    const scalarField& c
+) const
+{
+    return -1.0/T;
+}
+
+
 inline void Foam::thirdBodyArrheniusReactionRate::write(Ostream& os) const
 {
     ArrheniusReactionRate::write(os);

src/thermophysicalModels/specie/reaction/reactionRate/thirdBodyEfficiencies/thirdBodyEfficienciesI.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -83,7 +83,7 @@ inline Foam::thirdBodyEfficiencies::thirdBodyEfficiencies
 
 inline Foam::scalar Foam::thirdBodyEfficiencies::M(const scalarList& c) const
 {
-    scalar M = 0.0;
+    scalar M = 0;
     forAll(*this, i)
     {
         M += operator[](i)*c[i];

src/thermophysicalModels/specie/thermo/eConst/eConstThermo.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -155,6 +155,16 @@ public:
             inline scalar S(const scalar p, const scalar T) const;
 
 
+        // Derivative term used for Jacobian
+
+            //- Derivative of Gibbs free energy w.r.t. temperature
+            inline scalar dGdT(const scalar p, const scalar T) const;
+
+            //- Temperature derivative of heat capacity at constant pressure
+            inline scalar dCpdT(const scalar p, const scalar T) const;
+
+
+
         // I-O
 
             //- Write to Ostream

src/thermophysicalModels/specie/thermo/eConst/eConstThermoI.H

@@ -139,6 +139,28 @@ inline Foam::scalar Foam::eConstThermo<EquationOfState>::S
 }
 
 
+template<class EquationOfState>
+inline Foam::scalar Foam::eConstThermo<EquationOfState>::dGdT
+(
+    const scalar p,
+    const scalar T
+) const
+{
+    return 0;
+}
+
+
+template<class EquationOfState>
+inline Foam::scalar Foam::eConstThermo<EquationOfState>::dCpdT
+(
+    const scalar p,
+    const scalar T
+) const
+{
+    return 0;
+}
+
+
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
 template<class EquationOfState>

src/thermophysicalModels/specie/thermo/hConst/hConstThermo.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -148,6 +148,15 @@ public:
             inline scalar S(const scalar p, const scalar T) const;
 
 
+        // Derivative term used for Jacobian
+
+            //- Derivative of Gibbs free energy w.r.t. temperature
+            inline scalar dGdT(const scalar p, const scalar T) const;
+
+            //- Temperature derivative of heat capacity at constant pressure
+            inline scalar dCpdT(const scalar p, const scalar T) const;
+
+
         // I-O
 
             //- Write to Ostream

src/thermophysicalModels/specie/thermo/hConst/hConstThermoI.H

@@ -136,6 +136,26 @@ inline Foam::scalar Foam::hConstThermo<EquationOfState>::S
 }
 
 
+template<class EquationOfState>
+inline Foam::scalar Foam::hConstThermo<EquationOfState>::dGdT
+(
+    const scalar p, const scalar T
+) const
+{
+    return 0;
+}
+
+
+template<class EquationOfState>
+inline Foam::scalar Foam::hConstThermo<EquationOfState>::dCpdT
+(
+    const scalar p, const scalar T
+) const
+{
+    return 0;
+}
+
+
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
 template<class EquationOfState>

src/thermophysicalModels/specie/thermo/hPolynomial/hPolynomialThermo.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -198,6 +198,15 @@ public:
             inline scalar S(const scalar p, const scalar T) const;
 
 
+        // Derivative term used for Jacobian
+
+            //- Derivative of Gibbs free energy w.r.t. temperature
+            inline scalar dGdT(const scalar p, const scalar T) const;
+
+            //- Temperature derivative of heat capacity at constant pressure
+            inline scalar dCpdT(const scalar p, const scalar T) const;
+
+
         // I-O
 
             //- Write to Ostream

src/thermophysicalModels/specie/thermo/hPolynomial/hPolynomialThermoI.H

@@ -126,6 +126,36 @@ inline Foam::scalar Foam::hPolynomialThermo<EquationOfState, PolySize>::S
 }
 
 
+template<class EquationOfState, int PolySize>
+inline Foam::scalar Foam::hPolynomialThermo<EquationOfState, PolySize>::dGdT
+(
+    const scalar p,
+    const scalar T
+) const
+{
+    return
+    (
+        hCoeffs_.derivative(T)
+      - T*sCoeffs_.derivative(T)
+      - sCoeffs_.value(T)
+    );
+}
+
+
+template<class EquationOfState, int PolySize>
+inline Foam::scalar Foam::hPolynomialThermo<EquationOfState, PolySize>::dCpdT
+(
+    const scalar p,
+    const scalar T
+) const
+{
+    return
+    (
+        CpCoeffs_.derivative(T)
+    );
+}
+
+
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
 template<class EquationOfState, int PolySize>

src/thermophysicalModels/specie/thermo/hPower/hPowerThermo.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -169,6 +169,15 @@ public:
             inline scalar S(const scalar p, const scalar T) const;
 
 
+        // Derivative term used for Jacobian
+
+            //- Derivative of Gibbs free energy w.r.t. temperature
+            inline scalar dGdT(const scalar p, const scalar T) const;
+
+            //- Temperature derivative of heat capacity at constant pressure
+            inline scalar dCpdT(const scalar p, const scalar T) const;
+
+
     // Member operators
 
         inline void operator+=(const hPowerThermo&);

src/thermophysicalModels/specie/thermo/hPower/hPowerThermoI.H

@@ -164,6 +164,30 @@ inline Foam::scalar Foam::hPowerThermo<EquationOfState>::S
 }
 
 
+template<class EquationOfState>
+inline Foam::scalar Foam::hPowerThermo<EquationOfState>::dGdT
+(
+    const scalar p, const scalar T
+) const
+{
+    // To be implemented
+    NotImplemented;
+    return 0;
+}
+
+
+template<class EquationOfState>
+inline Foam::scalar Foam::hPowerThermo<EquationOfState>::dCpdT
+(
+    const scalar p, const scalar T
+) const
+{
+    // To be implemented
+    NotImplemented;
+    return 0;
+}
+
+
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
 template<class EquationOfState>

src/thermophysicalModels/specie/thermo/hRefConst/hRefConstThermo.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2015-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2015-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -152,6 +152,15 @@ public:
             inline scalar S(const scalar p, const scalar T) const;
 
 
+        // Derivative term used for Jacobian
+
+            //- Derivative of Gibbs free energy w.r.t. temperature
+            inline scalar dGdT(const scalar p, const scalar T) const;
+
+            //- Temperature derivative of heat capacity at constant pressure
+            inline scalar dCpdT(const scalar p, const scalar T) const;
+
+
         // I-O
 
             //- Write to Ostream

src/thermophysicalModels/specie/thermo/hRefConst/hRefConstThermoI.H

@@ -142,6 +142,26 @@ inline Foam::scalar Foam::hRefConstThermo<EquationOfState>::S
 }
 
 
+template<class EquationOfState>
+inline Foam::scalar Foam::hRefConstThermo<EquationOfState>::dGdT
+(
+    const scalar p, const scalar T
+) const
+{
+    return 0;
+}
+
+
+template<class EquationOfState>
+inline Foam::scalar Foam::hRefConstThermo<EquationOfState>::dCpdT
+(
+    const scalar p, const scalar T
+) const
+{
+    return 0;
+}
+
+
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
 template<class EquationOfState>

src/thermophysicalModels/specie/thermo/janaf/janafThermo.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -187,6 +187,15 @@ public:
             inline scalar S(const scalar p, const scalar T) const;
 
 
+        // Derivative term used for Jacobian
+
+            //- Derivative of Gibbs free energy w.r.t. temperature
+            inline scalar dGdT(const scalar p, const scalar T) const;
+
+            //- Temperature derivative of heat capacity at constant pressure
+            inline scalar dCpdT(const scalar p, const scalar T) const;
+
+
         // I-O
 
             //- Write to Ostream

src/thermophysicalModels/specie/thermo/janaf/janafThermoI.H

@@ -237,6 +237,31 @@ inline Foam::scalar Foam::janafThermo<EquationOfState>::S
 }
 
 
+template<class EquationOfState>
+inline Foam::scalar Foam::janafThermo<EquationOfState>::dGdT
+(
+    const scalar p,
+    const scalar T
+) const
+{
+    const coeffArray& a = coeffs(T);
+    return -((a[0] + a[5]/T)/T + a[1]/2 + T*(a[2]/3 + T*(a[3]/4 + T*a[4]/5)));
+}
+
+
+template<class EquationOfState>
+inline Foam::scalar Foam::janafThermo<EquationOfState>::dCpdT
+(
+    const scalar p,
+    const scalar T
+) const
+{
+    const coeffArray& a = coeffs(T);
+    return
+        (((4*a[4]*T + 3*a[3])*T + 2*a[2])*T + a[1]);
+}
+
+
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
 template<class EquationOfState>

src/thermophysicalModels/specie/thermo/thermo/thermo.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -22,7 +22,7 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::thermo
+    Foam::species::thermo
 
 Description
     Basic thermodynamics type based on the use of fitting functions for
@@ -318,6 +318,15 @@ public:
             ) const;
 
 
+        // Derivative term used for Jacobian
+
+            //- Derivative of B (acooding to Niemeyer et al.) w.r.t. temperature
+            inline scalar dKcdTbyKc(const scalar p, const scalar T) const;
+
+            //- Derivative of cp w.r.t. temperature
+            inline scalar dcpdT(const scalar p, const scalar T) const;
+
+
         // I-O
 
             //- Write to Ostream

src/thermophysicalModels/specie/thermo/thermo/thermoI.H

@@ -444,6 +444,35 @@ inline Foam::scalar Foam::species::thermo<Thermo, Type>::TEa
 }
 
 
+template<class Thermo, template<class> class Type>
+inline Foam::scalar
+Foam::species::thermo<Thermo, Type>::dKcdTbyKc
+(
+    const scalar p,
+    const scalar T
+) const
+{
+    const scalar nm = this->Y()/this->W();
+
+    if (equal(nm, small))
+    {
+        return -this->dGdT(Pstd, T)*this->Y()/RR;
+    }
+    else
+    {
+        return -(nm/T + this->dGdT(Pstd, T)*this->Y()/RR);
+    }
+}
+
+
+template<class Thermo, template<class> class Type>
+inline Foam::scalar
+Foam::species::thermo<Thermo, Type>::dcpdT(const scalar p, const scalar T) const
+{
+    return this->dCpdT(p, T)*this->W();;
+}
+
+
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
 template<class Thermo, template<class> class Type>

tutorials/combustion/chemFoam/nc7h16/constant/chemistryProperties

@@ -33,7 +33,7 @@ EulerImplicitCoeffs
 odeCoeffs
 {
     solver          seulex;
-    absTol          1e-14;
+    absTol          1e-12;
     relTol          1e-1;
 }
 

tutorials/combustion/fireFoam/LES/oppositeBurningPanels/constant/panelRegion/chemistryProperties

@@ -28,7 +28,7 @@ initialChemicalTimeStep 1e-07;
 odeCoeffs
 {
     solver          seulex;
-    absTol          1e-12;
+    absTol          1e-8;
     relTol          0.01;
 }
 

tutorials/combustion/reactingFoam/RAS/DLR_A_LTS/constant/chemistryProperties

@@ -35,9 +35,9 @@ initialChemicalTimeStep 1e-07;
 
 odeCoeffs
 {
-    solver          Rosenbrock34;  // Rosenbrock34, seulex or rodas23
+    solver          seulex;
-    absTol          1e-12;
+    absTol          1e-08;
-    relTol          1e-7;
+    relTol          0.1;
 }
 
 reduction
@@ -86,7 +86,7 @@ tabulation
     printNumRetrieve   off;
 
     // Tolerance used for retrieve and grow
-    tolerance   1e-4;
+    tolerance   0.003;
 
     // ISAT is the only method currently available
     method    ISAT;
@@ -97,24 +97,24 @@ tabulation
         otherSpecies 1;
         Temperature  1000;
         Pressure     1e15;
-        deltaT       0.5;
+        deltaT       1;
     }
 
     // Maximum number of leafs stored in the binary tree
-    maxNLeafs  2000;
+    maxNLeafs  5000;
 
     // Maximum life time of the leafs (in time steps) used in unsteady
     // simulations to force renewal of the stored chemPoints and keep the tree
     // small
-    chPMaxLifeTime 100;
+    chPMaxLifeTime 1000;
 
     // Maximum number of growth allowed on a chemPoint to avoid distorted
     // chemPoints
-    maxGrowth  10;
+    maxGrowth  100;
 
     // Number of time steps between analysis of the tree to remove old
     // chemPoints or try to balance it
-    checkEntireTreeInterval  5;
+    checkEntireTreeInterval  500;
 
     // Parameters used to decide whether to balance or not if the tree's depth
     // is larger than maxDepthFactor*log2(nLeafs) then balance the tree
@@ -123,7 +123,6 @@ tabulation
     // Try to balance the tree only if the size of the tree is greater
     minBalanceThreshold 30;
 
-
     // Activate the use of a MRU (most recently used) list
     MRURetrieve false;
 

tutorials/combustion/reactingFoam/RAS/DLR_A_LTS/constant/chemistryProperties.test

@@ -1,142 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  dev                                   |
-|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       dictionary;
-    location    "constant";
-    object      chemistryProperties;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-chemistryType
-{
-    solver            ode;
-    method            TDAC;
-}
-
-chemistry           on;
-
-importantSpecies
-{
-    CH4;
-    H2O;
-    O2;
-    CO2;
-}
-
-initialChemicalTimeStep 1e-07;
-
-odeCoeffs
-{
-    solver          Rosenbrock34;  // Rosenbrock34, seulex or rodas23
-    absTol          1e-12;
-    relTol          1e-7;
-}
-
-reduction
-{
-    // Activate reduction
-    active  on;
-
-    // Switch logging of the reduction statistics and performance
-    log         on;
-
-    // Tolerance depends on the reduction method, see details for each method
-    tolerance   1e-4;
-
-    // Available methods: DRG, DAC, DRGEP, PFA, EFA
-    method DAC;
-
-    // Search initiating set (SIS) of species, needed for most methods
-    initialSet
-    {
-        CO;
-        CH4;
-        HO2;
-    }
-
-    // For DAC, option to automatically change the SIS switch from HO2 to H2O
-    // and CO to CO2, + disable fuel
-    automaticSIS    off;
-
-    // When automaticSIS, the method needs to know the fuel
-    fuelSpecies
-    {
-        CH4 1;
-    }
-}
-
-tabulation
-{
-    // Activate tabulation
-    active      on;
-
-    // Switch logging of the tabulation statistics and performance
-    log         on;
-
-    printProportion    off;
-
-    printNumRetrieve   off;
-
-    // Tolerance used for retrieve and grow
-    tolerance   0.003;
-
-    // ISAT is the only method currently available
-    method    ISAT;
-
-    // Scale factors used in the definition of the ellipsoid of accuracy
-    scaleFactor
-    {
-        otherSpecies 1;
-        Temperature  10000;
-        Pressure     1e15;
-        deltaT       1;
-    }
-
-    // Maximum number of leafs stored in the binary tree
-    maxNLeafs  5000;
-
-    // Maximum life time of the leafs (in time steps) used in unsteady
-    // simulations to force renewal of the stored chemPoints and keep the tree
-    // small
-    chPMaxLifeTime 1000;
-
-    // Maximum number of growth allowed on a chemPoint to avoid distorted
-    // chemPoints
-    maxGrowth  100;
-
-    // Number of time steps between analysis of the tree to remove old
-    // chemPoints or try to balance it
-    checkEntireTreeInterval  500;
-
-    // Parameters used to decide whether to balance or not if the tree's depth
-    // is larger than maxDepthFactor*log2(nLeafs) then balance the tree
-    maxDepthFactor   2;
-
-    // Try to balance the tree only if the size of the tree is greater
-    minBalanceThreshold 30;
-
-    // Activate the use of a MRU (most recently used) list
-    MRURetrieve false;
-
-    // Maximum size of the MRU list
-    maxMRUSize 0;
-
-    // Allow to grow points
-    growPoints  true;
-
-    // When mechanism reduction is used, new dimensions might be added
-    // maxNumNewDim set the maximum number of new dimensions added during a
-    // growth
-    maxNumNewDim 10;
-}
-
-
-// ************************************************************************* //

tutorials/combustion/reactingFoam/RAS/DLR_A_LTS/constant/reactionsGRI

@@ -3608,3 +3608,6 @@ reactions
         Ta              0;
     }
 }
+
+Tlow            250;
+Thigh           5000;

tutorials/combustion/reactingFoam/RAS/SandiaD_LTS/Allrun

@@ -14,16 +14,16 @@ runApplication chemkinToFoam \
 runApplication blockMesh
 runApplication setFields
 
-
 # Run the application without chemistry until 1500 to let the flow field develop
+foamDictionary -entry "startTime" -set "0" system/controlDict
 foamDictionary -entry "writeInterval" -set "1500" system/controlDict
 foamDictionary -entry "endTime" -set "1500" system/controlDict
 foamDictionary -entry "chemistry" -set "off" constant/chemistryProperties
 
 runApplication $application
 
-
 # Run with chemistry until flame reach its full size
+foamDictionary -entry "startTime" -set "1500" system/controlDict
 foamDictionary -entry "writeInterval" -set "100" system/controlDict
 foamDictionary -entry "endTime" -set "5000" system/controlDict
 foamDictionary -entry "chemistry" -set "on" constant/chemistryProperties

tutorials/combustion/reactingFoam/RAS/SandiaD_LTS/constant/chemistryProperties

@@ -20,7 +20,7 @@ chemistryType
     method          TDAC;
 }
 
-chemistry       off;
+chemistry       on;
 
 importantSpecies
 {
@@ -35,8 +35,8 @@ initialChemicalTimeStep 1e-07;
 odeCoeffs
 {
     solver          seulex;
-    absTol          1e-12;
+    absTol          1e-08;
-    relTol          1e-07;
+    relTol          0.1;
 }
 
 reduction

tutorials/combustion/reactingFoam/RAS/SandiaD_LTS/constant/reactionsGRI

@@ -3608,3 +3608,6 @@ reactions
         Ta              0;
     }
 }
+
+Tlow            250;
+Thigh           5000;

tutorials/combustion/reactingFoam/RAS/SandiaD_LTS/system/controlDict

@@ -16,19 +16,19 @@ FoamFile
 
 application     reactingFoam;
 
-startFrom       latestTime;
+startFrom       startTime;
 
-startTime       0;
+startTime       1500;
 
 stopAt          endTime;
 
-endTime         1500;
+endTime         5000;
 
 deltaT          1;
 
 writeControl    runTime;
 
-writeInterval   1500;
+writeInterval   100;
 
 purgeWrite      0;
 

tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D/constant/chemistryProperties

@@ -33,7 +33,7 @@ EulerImplicitCoeffs
 odeCoeffs
 {
     solver          Rosenbrock43;
-    absTol          1e-12;
+    absTol          1e-8;
     relTol          0.01;
 }
 

tutorials/combustion/reactingFoam/laminar/counterFlowFlame2DLTS/constant/chemistryProperties

@@ -33,7 +33,7 @@ EulerImplicitCoeffs
 odeCoeffs
 {
     solver          Rosenbrock43;
-    absTol          1e-12;
+    absTol          1e-8;
     relTol          0.01;
 }
 

tutorials/combustion/reactingFoam/laminar/counterFlowFlame2DLTS_GRI_TDAC/constant/chemistryProperties

@@ -24,11 +24,12 @@ chemistryType
 chemistry       on;
 
 initialChemicalTimeStep 1e-7;
+//maxChemicalTimeStep 1E-3;
 
 odeCoeffs
 {
     solver          seulex;
-    absTol          1e-12;
+    absTol          1e-8;
     relTol          1e-1;
 }
 
@@ -78,7 +79,7 @@ tabulation
     printNumRetrieve   off;
 
     // Tolerance used for retrieve and grow
-    tolerance   1e-3;
+    tolerance   3e-3;
 
     // ISAT is the only method currently available
     method    ISAT;
@@ -87,7 +88,7 @@ tabulation
     scaleFactor
     {
         otherSpecies 1;
-        Temperature  25000;
+        Temperature  10000;
         Pressure     1e15;
         deltaT       1;
     }

tutorials/combustion/reactingFoam/laminar/counterFlowFlame2DLTS_GRI_TDAC/constant/reactionsGRI

@@ -5588,3 +5588,6 @@ reactions
         }
     }
 }
+
+Thigh   5000;
+Tlow    200;

tutorials/combustion/reactingFoam/laminar/counterFlowFlame2DLTS_GRI_TDAC/system/controlDict

@@ -23,7 +23,7 @@ startTime       0;
 
 stopAt          endTime;
 
-endTime         1000;
+endTime         1500;
 
 deltaT          1;
 

tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D_GRI/constant/chemistryProperties

@@ -7,13 +7,12 @@
 \*---------------------------------------------------------------------------*/
 FoamFile
 {
-    version     2.0;
+    version         2;
     format          ascii;
     class           dictionary;
     location        "constant";
     object          chemistryProperties;
 }
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 chemistryType
 {
@@ -22,14 +21,14 @@ chemistryType
 
 chemistry       on;
 
-initialChemicalTimeStep 1e-7;
+initialChemicalTimeStep 1e-07;
-maxChemicalTimeStep     1e-4;
 
 odeCoeffs
 {
     solver          seulex;
-    absTol          1e-12;
+    absTol          1e-8;
-    relTol          1e-1;
+    relTol          0.1;
 }
 
+
 // ************************************************************************* //

tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D_GRI/constant/reactionsGRI

@@ -5588,3 +5588,6 @@ reactions
         }
     }
 }
+
+Thigh   5000;
+Tlow    200;

tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D_GRI_TDAC/constant/chemistryProperties

@@ -24,12 +24,12 @@ chemistryType
 chemistry       on;
 
 initialChemicalTimeStep 1e-7;
-maxChemicalTimeStep     1e-4;
+//maxChemicalTimeStep 1E-3;
 
 odeCoeffs
 {
     solver          seulex;
-    absTol          1e-12;
+    absTol          1e-8;
     relTol          1e-1;
 }
 
@@ -79,7 +79,7 @@ tabulation
     printNumRetrieve   off;
 
     // Tolerance used for retrieve and grow
-    tolerance   1e-4;
+    tolerance   3e-3;
 
     // ISAT is the only method currently available
     method    ISAT;
@@ -88,7 +88,7 @@ tabulation
     scaleFactor
     {
         otherSpecies 1;
-        Temperature  2500;
+        Temperature  10000;
         Pressure     1e15;
         deltaT       1;
     }

tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D_GRI_TDAC/constant/chemistryProperties.new

@@ -28,7 +28,7 @@ initialChemicalTimeStep 1e-7;
 odeCoeffs
 {
     solver          seulex;
-    absTol          1e-12;
+    absTol          1e-8;
     relTol          1e-1;
 }
 

tutorials/combustion/reactingFoam/laminar/counterFlowFlame2D_GRI_TDAC/constant/reactionsGRI

@@ -5588,3 +5588,6 @@ reactions
         }
     }
 }
+
+Thigh   5000;
+Tlow    200;

tutorials/lagrangian/reactingParcelFoam/counterFlowFlame2DLTS/constant/chemistryProperties

@@ -33,7 +33,7 @@ EulerImplicitCoeffs
 odeCoeffs
 {
     solver          Rosenbrock43;
-    absTol          1e-12;
+    absTol          1e-8;
     relTol          0.01;
 }
 

tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/bubbleColumnEvaporatingReacting/constant/chemistryProperties.gas

@@ -33,7 +33,7 @@ EulerImplicitCoeffs
 odeCoeffs
 {
     solver          Rosenbrock43;
-    absTol          1e-12;
+    absTol          1e-8;
     relTol          0.01;
 }