chemistryModel, reactions, ODESolver: Added "li" argument to all functions which evaluate the reaction rate

In chemistryModel "li" is set to the current cell index but for other reacting
systems it should be set to the current index of the element for which the
reaction system is being evaluated.

In the ODESolver "li" is the current index of the element for which the ODE
system is being solved if there is a list of related systems being solved,
otherwise it can be set to 0.

src/ODE/ODESolvers/Euler/Euler.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2013-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2019 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -68,6 +68,7 @@ Foam::scalar Foam::Euler::solve
 (
     const scalar x0,
     const scalarField& y0,
+    const label li,
     const scalarField& dydx0,
     const scalar dx,
     scalarField& y
@@ -93,10 +94,11 @@ void Foam::Euler::solve
 (
     scalar& x,
     scalarField& y,
+    const label li,
     scalar& dxTry
 ) const
 {
-    adaptiveSolver::solve(odes_, x, y, dxTry);
+    adaptiveSolver::solve(odes_, x, y, li, dxTry);
 }
 
 

src/ODE/ODESolvers/Euler/Euler.H

@@ -97,6 +97,7 @@ public:
         (
             const scalar x0,
             const scalarField& y0,
+            const label li,
             const scalarField& dydx0,
             const scalar dx,
             scalarField& y
@@ -107,6 +108,7 @@ public:
         (
             scalar& x,
             scalarField& y,
+            const label li,
             scalar& dxTry
         ) const;
 };

src/ODE/ODESolvers/EulerSI/EulerSI.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2013-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2019 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -78,12 +78,13 @@ Foam::scalar Foam::EulerSI::solve
 (
     const scalar x0,
     const scalarField& y0,
+    const label li,
     const scalarField& dydx0,
     const scalar dx,
     scalarField& y
 ) const
 {
-    odes_.jacobian(x0, y0, dfdx_, dfdy_);
+    odes_.jacobian(x0, y0, li, dfdx_, dfdy_);
 
     for (label i=0; i<n_; i++)
     {
@@ -118,10 +119,11 @@ void Foam::EulerSI::solve
 (
     scalar& x,
     scalarField& y,
+    const label li,
     scalar& dxTry
 ) const
 {
-    adaptiveSolver::solve(odes_, x, y, dxTry);
+    adaptiveSolver::solve(odes_, x, y, li, dxTry);
 }
 
 

src/ODE/ODESolvers/EulerSI/EulerSI.H

@@ -104,6 +104,7 @@ public:
         (
             const scalar x0,
             const scalarField& y0,
+            const label li,
             const scalarField& dydx0,
             const scalar dx,
             scalarField& y
@@ -114,6 +115,7 @@ public:
         (
             scalar& x,
             scalarField& y,
+            const label li,
             scalar& dxTry
         ) const;
 };

src/ODE/ODESolvers/ODESolver/ODESolver.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2019 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -116,11 +116,12 @@ void Foam::ODESolver::solve
 (
     scalar& x,
     scalarField& y,
+    const label li,
     scalar& dxTry
 ) const
 {
     stepState step(dxTry);
-    solve(x, y, step);
+    solve(x, y, li, step);
     dxTry = step.dxTry;
 }
 
@@ -129,11 +130,12 @@ void Foam::ODESolver::solve
 (
     scalar& x,
     scalarField& y,
+    const label li,
     stepState& step
 ) const
 {
     scalar x0 = x;
-    solve(x, y, step.dxTry);
+    solve(x, y, li, step.dxTry);
     step.dxDid = x - x0;
 }
 
@@ -143,6 +145,7 @@ void Foam::ODESolver::solve
     const scalar xStart,
     const scalar xEnd,
     scalarField& y,
+    const label li,
     scalar& dxTry
 ) const
 {
@@ -164,7 +167,7 @@ void Foam::ODESolver::solve
         }
 
         // Integrate as far as possible up to step.dxTry
-        solve(x, y, step);
+        solve(x, y, li, step);
 
         // Check if reached xEnd
         if ((x - xEnd)*(xEnd - xStart) >= 0)

src/ODE/ODESolvers/ODESolver/ODESolver.H

@@ -183,35 +183,42 @@ public:
 
         inline void resizeMatrix(scalarSquareMatrix& m) const;
 
-        //- Solve the ODE system as far as possible up to dxTry
+        //- Solve the ODE system from the current state xStart, y
-        //  adjusting the step as necessary to provide a solution within
+        //  and the optional index into the list of systems to solve li
-        //  the specified tolerance.
+        //  as far as possible up to dxTry adjusting the step as necessary
+        //  to provide a solution within the specified tolerance.
         //  Update the state and return an estimate for the next step in dxTry
         virtual void solve
         (
             scalar& x,
             scalarField& y,
+            const label li,
             scalar& dxTry
         ) const;
 
-        //- Solve the ODE system as far as possible up to dxTry
+        //- Solve the ODE system from the current state xStart, y
-        //  adjusting the step as necessary to provide a solution within
+        //  and the optional index into the list of systems to solve li
-        //  the specified tolerance.
+        //  as far as possible up to step.dxTry adjusting the step as necessary
-        //  Update the state and return an estimate for the next step in dxTry
+        //  to provide a solution within the specified tolerance.
+        //  Update the state and return the step actually taken in step.dxDid
+        //  and an estimate for the next step in step.dxTry
         virtual void solve
         (
             scalar& x,
             scalarField& y,
+            const label li,
             stepState& step
         ) const;
 
-        //- Solve the ODE system from xStart to xEnd, update the state
+        //- Solve the ODE system from the current state xStart, y
-        //  and return an estimate for the next step in dxTry
+        //  and the optional index into the list of systems to solve li
+        //  to xEnd and return an estimate for the next step in dxTry
         virtual void solve
         (
             const scalar xStart,
             const scalar xEnd,
             scalarField& y,
+            const label li,
             scalar& dxEst
         ) const;
 

src/ODE/ODESolvers/RKCK45/RKCK45.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2013-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2019 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -114,6 +114,7 @@ Foam::scalar Foam::RKCK45::solve
 (
     const scalar x0,
     const scalarField& y0,
+    const label li,
     const scalarField& dydx0,
     const scalar dx,
     scalarField& y
@@ -124,21 +125,21 @@ Foam::scalar Foam::RKCK45::solve
         yTemp_[i] = y0[i] + a21*dx*dydx0[i];
     }
 
-    odes_.derivatives(x0 + c2*dx, yTemp_, k2_);
+    odes_.derivatives(x0 + c2*dx, yTemp_, li, k2_);
 
     forAll(yTemp_, i)
     {
         yTemp_[i] = y0[i] + dx*(a31*dydx0[i] + a32*k2_[i]);
     }
 
-    odes_.derivatives(x0 + c3*dx, yTemp_, k3_);
+    odes_.derivatives(x0 + c3*dx, yTemp_, li, k3_);
 
     forAll(yTemp_, i)
     {
         yTemp_[i] = y0[i] + dx*(a41*dydx0[i] + a42*k2_[i] + a43*k3_[i]);
     }
 
-    odes_.derivatives(x0 + c4*dx, yTemp_, k4_);
+    odes_.derivatives(x0 + c4*dx, yTemp_, li, k4_);
 
     forAll(yTemp_, i)
     {
@@ -146,7 +147,7 @@ Foam::scalar Foam::RKCK45::solve
           + dx*(a51*dydx0[i] + a52*k2_[i] + a53*k3_[i] + a54*k4_[i]);
     }
 
-    odes_.derivatives(x0 + c5*dx, yTemp_, k5_);
+    odes_.derivatives(x0 + c5*dx, yTemp_, li, k5_);
 
     forAll(yTemp_, i)
     {
@@ -155,7 +156,7 @@ Foam::scalar Foam::RKCK45::solve
            *(a61*dydx0[i] + a62*k2_[i] + a63*k3_[i] + a64*k4_[i] + a65*k5_[i]);
     }
 
-    odes_.derivatives(x0 + c6*dx, yTemp_, k6_);
+    odes_.derivatives(x0 + c6*dx, yTemp_, li, k6_);
 
     forAll(y, i)
     {
@@ -178,10 +179,11 @@ void Foam::RKCK45::solve
 (
     scalar& x,
     scalarField& y,
+    const label li,
     scalar& dxTry
 ) const
 {
-    adaptiveSolver::solve(odes_, x, y, dxTry);
+    adaptiveSolver::solve(odes_, x, y, li, dxTry);
 }
 
 

src/ODE/ODESolvers/RKCK45/RKCK45.H

@@ -117,6 +117,7 @@ public:
         (
             const scalar x0,
             const scalarField& y0,
+            const label li,
             const scalarField& dydx0,
             const scalar dx,
             scalarField& y
@@ -127,6 +128,7 @@ public:
         (
             scalar& x,
             scalarField& y,
+            const label li,
             scalar& dxTry
         ) const;
 };

src/ODE/ODESolvers/RKDP45/RKDP45.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2013-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2019 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -119,6 +119,7 @@ Foam::scalar Foam::RKDP45::solve
 (
     const scalar x0,
     const scalarField& y0,
+    const label li,
     const scalarField& dydx0,
     const scalar dx,
     scalarField& y
@@ -129,21 +130,21 @@ Foam::scalar Foam::RKDP45::solve
         yTemp_[i] = y0[i] + a21*dx*dydx0[i];
     }
 
-    odes_.derivatives(x0 + c2*dx, yTemp_, k2_);
+    odes_.derivatives(x0 + c2*dx, yTemp_, li, k2_);
 
     forAll(yTemp_, i)
     {
         yTemp_[i] = y0[i] + dx*(a31*dydx0[i] + a32*k2_[i]);
     }
 
-    odes_.derivatives(x0 + c3*dx, yTemp_, k3_);
+    odes_.derivatives(x0 + c3*dx, yTemp_, li, k3_);
 
     forAll(yTemp_, i)
     {
         yTemp_[i] = y0[i] + dx*(a41*dydx0[i] + a42*k2_[i] + a43*k3_[i]);
     }
 
-    odes_.derivatives(x0 + c4*dx, yTemp_, k4_);
+    odes_.derivatives(x0 + c4*dx, yTemp_, li, k4_);
 
     forAll(yTemp_, i)
     {
@@ -151,7 +152,7 @@ Foam::scalar Foam::RKDP45::solve
           + dx*(a51*dydx0[i] + a52*k2_[i] + a53*k3_[i] + a54*k4_[i]);
     }
 
-    odes_.derivatives(x0 + c5*dx, yTemp_, k5_);
+    odes_.derivatives(x0 + c5*dx, yTemp_, li, k5_);
 
     forAll(yTemp_, i)
     {
@@ -160,7 +161,7 @@ Foam::scalar Foam::RKDP45::solve
            *(a61*dydx0[i] + a62*k2_[i] + a63*k3_[i] + a64*k4_[i] + a65*k5_[i]);
     }
 
-    odes_.derivatives(x0 + dx, yTemp_, k6_);
+    odes_.derivatives(x0 + dx, yTemp_, li, k6_);
 
     forAll(y, i)
     {
@@ -169,7 +170,7 @@ Foam::scalar Foam::RKDP45::solve
     }
 
     // Reuse k2_ for the derivative of the new state
-    odes_.derivatives(x0 + dx, y, k2_);
+    odes_.derivatives(x0 + dx, y, li, k2_);
 
     forAll(err_, i)
     {
@@ -187,10 +188,11 @@ void Foam::RKDP45::solve
 (
     scalar& x,
     scalarField& y,
+    const label li,
     scalar& dxTry
 ) const
 {
-    adaptiveSolver::solve(odes_, x, y, dxTry);
+    adaptiveSolver::solve(odes_, x, y, li, dxTry);
 }
 
 

src/ODE/ODESolvers/RKDP45/RKDP45.H

@@ -120,6 +120,7 @@ public:
         (
             const scalar x0,
             const scalarField& y0,
+            const label li,
             const scalarField& dydx0,
             const scalar dx,
             scalarField& y
@@ -130,6 +131,7 @@ public:
         (
             scalar& x,
             scalarField& y,
+            const label li,
             scalar& dxTry
         ) const;
 };

src/ODE/ODESolvers/RKF45/RKF45.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2013-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2019 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -115,6 +115,7 @@ Foam::scalar Foam::RKF45::solve
 (
     const scalar x0,
     const scalarField& y0,
+    const label li,
     const scalarField& dydx0,
     const scalar dx,
     scalarField& y
@@ -125,21 +126,21 @@ Foam::scalar Foam::RKF45::solve
         yTemp_[i] = y0[i] + a21*dx*dydx0[i];
     }
 
-    odes_.derivatives(x0 + c2*dx, yTemp_, k2_);
+    odes_.derivatives(x0 + c2*dx, yTemp_, li, k2_);
 
     forAll(yTemp_, i)
     {
         yTemp_[i] = y0[i] + dx*(a31*dydx0[i] + a32*k2_[i]);
     }
 
-    odes_.derivatives(x0 + c3*dx, yTemp_, k3_);
+    odes_.derivatives(x0 + c3*dx, yTemp_, li, k3_);
 
     forAll(yTemp_, i)
     {
         yTemp_[i] = y0[i] + dx*(a41*dydx0[i] + a42*k2_[i] + a43*k3_[i]);
     }
 
-    odes_.derivatives(x0 + c4*dx, yTemp_, k4_);
+    odes_.derivatives(x0 + c4*dx, yTemp_, li, k4_);
 
     forAll(yTemp_, i)
     {
@@ -147,7 +148,7 @@ Foam::scalar Foam::RKF45::solve
           + dx*(a51*dydx0[i] + a52*k2_[i] + a53*k3_[i] + a54*k4_[i]);
     }
 
-    odes_.derivatives(x0 + c5*dx, yTemp_, k5_);
+    odes_.derivatives(x0 + c5*dx, yTemp_, li, k5_);
 
     forAll(yTemp_, i)
     {
@@ -156,7 +157,7 @@ Foam::scalar Foam::RKF45::solve
            *(a61*dydx0[i] + a62*k2_[i] + a63*k3_[i] + a64*k4_[i] + a65*k5_[i]);
     }
 
-    odes_.derivatives(x0 + c6*dx, yTemp_, k6_);
+    odes_.derivatives(x0 + c6*dx, yTemp_, li, k6_);
 
     // Calculate the 5th-order solution
     forAll(y, i)
@@ -183,10 +184,11 @@ void Foam::RKF45::solve
 (
     scalar& x,
     scalarField& y,
+    const label li,
     scalar& dxTry
 ) const
 {
-    adaptiveSolver::solve(odes_, x, y, dxTry);
+    adaptiveSolver::solve(odes_, x, y, li, dxTry);
 }
 
 

src/ODE/ODESolvers/RKF45/RKF45.H

@@ -1,4 +1,4 @@
-/*---------------------------------------------------------------------------*\
+ /*---------------------------------------------------------------------------*\
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
@@ -122,6 +122,7 @@ public:
         (
             const scalar x0,
             const scalarField& y0,
+            const label li,
             const scalarField& dydx0,
             const scalar dx,
             scalarField& y
@@ -132,6 +133,7 @@ public:
         (
             scalar& x,
             scalarField& y,
+            const label li,
             scalar& dxTry
         ) const;
 };

src/ODE/ODESolvers/Rosenbrock12/Rosenbrock12.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2013-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2019 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -94,12 +94,13 @@ Foam::scalar Foam::Rosenbrock12::solve
 (
     const scalar x0,
     const scalarField& y0,
+    const label li,
     const scalarField& dydx0,
     const scalar dx,
     scalarField& y
 ) const
 {
-    odes_.jacobian(x0, y0, dfdx_, dfdy_);
+    odes_.jacobian(x0, y0, li, dfdx_, dfdy_);
 
     for (label i=0; i<n_; i++)
     {
@@ -127,7 +128,7 @@ Foam::scalar Foam::Rosenbrock12::solve
         y[i] = y0[i] + a21*k1_[i];
     }
 
-    odes_.derivatives(x0 + c2*dx, y, dydx_);
+    odes_.derivatives(x0 + c2*dx, y, li, dydx_);
 
     forAll(k2_, i)
     {
@@ -151,10 +152,11 @@ void Foam::Rosenbrock12::solve
 (
     scalar& x,
     scalarField& y,
+    const label li,
     scalar& dxTry
 ) const
 {
-    adaptiveSolver::solve(odes_, x, y, dxTry);
+    adaptiveSolver::solve(odes_, x, y, li, dxTry);
 }
 
 

src/ODE/ODESolvers/Rosenbrock12/Rosenbrock12.H

@@ -111,6 +111,7 @@ public:
         (
             const scalar x0,
             const scalarField& y0,
+            const label li,
             const scalarField& dydx0,
             const scalar dx,
             scalarField& y
@@ -121,6 +122,7 @@ public:
         (
             scalar& x,
             scalarField& y,
+            const label li,
             scalar& dxTry
         ) const;
 };

src/ODE/ODESolvers/Rosenbrock23/Rosenbrock23.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2013-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2019 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -108,12 +108,13 @@ Foam::scalar Foam::Rosenbrock23::solve
 (
     const scalar x0,
     const scalarField& y0,
+    const label li,
     const scalarField& dydx0,
     const scalar dx,
     scalarField& y
 ) const
 {
-    odes_.jacobian(x0, y0, dfdx_, dfdy_);
+    odes_.jacobian(x0, y0, li, dfdx_, dfdy_);
 
     for (label i=0; i<n_; i++)
     {
@@ -141,7 +142,7 @@ Foam::scalar Foam::Rosenbrock23::solve
         y[i] = y0[i] + a21*k1_[i];
     }
 
-    odes_.derivatives(x0 + c2*dx, y, dydx_);
+    odes_.derivatives(x0 + c2*dx, y, li, dydx_);
 
     forAll(k2_, i)
     {
@@ -174,10 +175,11 @@ void Foam::Rosenbrock23::solve
 (
     scalar& x,
     scalarField& y,
+    const label li,
     scalar& dxTry
 ) const
 {
-    adaptiveSolver::solve(odes_, x, y, dxTry);
+    adaptiveSolver::solve(odes_, x, y, li, dxTry);
 }
 
 

src/ODE/ODESolvers/Rosenbrock23/Rosenbrock23.H

@@ -113,6 +113,7 @@ public:
         (
             const scalar x0,
             const scalarField& y0,
+            const label li,
             const scalarField& dydx0,
             const scalar dx,
             scalarField& y
@@ -123,6 +124,7 @@ public:
         (
             scalar& x,
             scalarField& y,
+            const label li,
             scalar& dxTry
         ) const;
 };

src/ODE/ODESolvers/Rosenbrock34/Rosenbrock34.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2013-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2019 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -154,12 +154,13 @@ Foam::scalar Foam::Rosenbrock34::solve
 (
     const scalar x0,
     const scalarField& y0,
+    const label li,
     const scalarField& dydx0,
     const scalar dx,
     scalarField& y
 ) const
 {
-    odes_.jacobian(x0, y0, dfdx_, dfdy_);
+    odes_.jacobian(x0, y0, li, dfdx_, dfdy_);
 
     for (label i=0; i<n_; i++)
     {
@@ -187,7 +188,7 @@ Foam::scalar Foam::Rosenbrock34::solve
         y[i] = y0[i] + a21*k1_[i];
     }
 
-    odes_.derivatives(x0 + c2*dx, y, dydx_);
+    odes_.derivatives(x0 + c2*dx, y, li, dydx_);
 
     forAll(k2_, i)
     {
@@ -202,7 +203,7 @@ Foam::scalar Foam::Rosenbrock34::solve
         y[i] = y0[i] + a31*k1_[i] + a32*k2_[i];
     }
 
-    odes_.derivatives(x0 + c3*dx, y, dydx_);
+    odes_.derivatives(x0 + c3*dx, y, li, dydx_);
 
     forAll(k3_, i)
     {
@@ -235,10 +236,11 @@ void Foam::Rosenbrock34::solve
 (
     scalar& x,
     scalarField& y,
+    const label li,
     scalar& dxTry
 ) const
 {
-    adaptiveSolver::solve(odes_, x, y, dxTry);
+    adaptiveSolver::solve(odes_, x, y, li, dxTry);
 }
 
 

src/ODE/ODESolvers/Rosenbrock34/Rosenbrock34.H

@@ -124,6 +124,7 @@ public:
         (
             const scalar x0,
             const scalarField& y0,
+            const label li,
             const scalarField& dydx0,
             const scalar dx,
             scalarField& y
@@ -134,6 +135,7 @@ public:
         (
             scalar& x,
             scalarField& y,
+            const label li,
             scalar& dxTry
         ) const;
 };

src/ODE/ODESolvers/SIBS/SIBS.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2019 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -92,10 +92,11 @@ void Foam::SIBS::solve
 (
     scalar& x,
     scalarField& y,
+    const label li,
     scalar& dxTry
 ) const
 {
-    odes_.derivatives(x, y, dydx0_);
+    odes_.derivatives(x, y, li, dydx0_);
 
     scalar h = dxTry;
     bool exitflag = false;
@@ -143,7 +144,7 @@ void Foam::SIBS::solve
     label k = 0;
     yTemp_ = y;
 
-    odes_.jacobian(x, y, dfdx_, dfdy_);
+    odes_.jacobian(x, y, li, dfdx_, dfdy_);
 
     if (x != xNew_ || h != dxTry)
     {
@@ -170,7 +171,7 @@ void Foam::SIBS::solve
                     << exit(FatalError);
             }
 
-            SIMPR(x, yTemp_, dydx0_, dfdx_, dfdy_, h, nSeq_[k], ySeq_);
+            SIMPR(x, yTemp_, li, dydx0_, dfdx_, dfdy_, h, nSeq_[k], ySeq_);
             scalar xest = sqr(h/nSeq_[k]);
 
             polyExtrapolate(k, xest, ySeq_, y, yErr_, x_p_, d_p_);

src/ODE/ODESolvers/SIBS/SIBS.H

@@ -90,6 +90,7 @@ class SIBS
         (
             const scalar xStart,
             const scalarField& y,
+            const label li,
             const scalarField& dydx,
             const scalarField& dfdx,
             const scalarSquareMatrix& dfdy,
@@ -136,6 +137,7 @@ public:
         (
             scalar& x,
             scalarField& y,
+            const label li,
             scalar& dxTry
         ) const;
 };

src/ODE/ODESolvers/SIBS/SIMPR.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2019 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -31,6 +31,7 @@ void Foam::SIBS::SIMPR
 (
     const scalar xStart,
     const scalarField& y,
+    const label li,
     const scalarField& dydx,
     const scalarField& dfdx,
     const scalarSquareMatrix& dfdy,
@@ -71,7 +72,7 @@ void Foam::SIBS::SIMPR
 
     scalar x = xStart + h;
 
-    odes_.derivatives(x, ytemp, yEnd);
+    odes_.derivatives(x, ytemp, li, yEnd);
 
     for (label nn=2; nn<=nSteps; nn++)
     {
@@ -89,7 +90,7 @@ void Foam::SIBS::SIMPR
 
         x += h;
 
-        odes_.derivatives(x, ytemp, yEnd);
+        odes_.derivatives(x, ytemp, li, yEnd);
     }
     for (label i=0; i<n_; i++)
     {

src/ODE/ODESolvers/Trapezoid/Trapezoid.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2013-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2019 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -68,6 +68,7 @@ Foam::scalar Foam::Trapezoid::solve
 (
     const scalar x0,
     const scalarField& y0,
+    const label li,
     const scalarField& dydx0,
     const scalar dx,
     scalarField& y
@@ -80,7 +81,7 @@ Foam::scalar Foam::Trapezoid::solve
     }
 
     // Evaluate the system for the predicted state
-    odes_.derivatives(x0 + dx, y, err_);
+    odes_.derivatives(x0 + dx, y, li, err_);
 
     // Update the state as the average between the prediction and the correction
     // and estimate the error from the difference
@@ -98,10 +99,11 @@ void Foam::Trapezoid::solve
 (
     scalar& x,
     scalarField& y,
+    const label li,
     scalar& dxTry
 ) const
 {
-    adaptiveSolver::solve(odes_, x, y, dxTry);
+    adaptiveSolver::solve(odes_, x, y, li, dxTry);
 }
 
 

src/ODE/ODESolvers/Trapezoid/Trapezoid.H

@@ -87,6 +87,7 @@ public:
         (
             const scalar x0,
             const scalarField& y0,
+            const label li,
             const scalarField& dydx0,
             const scalar dx,
             scalarField& y
@@ -97,6 +98,7 @@ public:
         (
             scalar& x,
             scalarField& y,
+            const label li,
             scalar& dxTry
         ) const;
 };

src/ODE/ODESolvers/adaptiveSolver/adaptiveSolver.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2019 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -60,20 +60,21 @@ void Foam::adaptiveSolver::solve
     const ODESystem& odes,
     scalar& x,
     scalarField& y,
+    const label li,
     scalar& dxTry
 ) const
 {
     scalar dx = dxTry;
     scalar err = 0.0;
 
-    odes.derivatives(x, y, dydx0_);
+    odes.derivatives(x, y, li, dydx0_);
 
     // Loop over solver and adjust step-size as necessary
     // to achieve desired error
     do
     {
         // Solve step and provide error estimate
-        err = solve(x, y, dydx0_, dx, yTemp_);
+        err = solve(x, y, li, dydx0_, dx, yTemp_);
 
         // If error is large reduce dx
         if (err > 1)

src/ODE/ODESolvers/adaptiveSolver/adaptiveSolver.H

@@ -82,6 +82,7 @@ public:
         (
             const scalar x0,
             const scalarField& y0,
+            const label li,
             const scalarField& dydx0,
             const scalar dx,
             scalarField& y
@@ -93,6 +94,7 @@ public:
             const ODESystem& ode,
             scalar& x,
             scalarField& y,
+            const label li,
             scalar& dxTry
         ) const;
 };

src/ODE/ODESolvers/rodas23/rodas23.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2013-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2019 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -100,12 +100,13 @@ Foam::scalar Foam::rodas23::solve
 (
     const scalar x0,
     const scalarField& y0,
+    const label li,
     const scalarField& dydx0,
     const scalar dx,
     scalarField& y
 ) const
 {
-    odes_.jacobian(x0, y0, dfdx_, dfdy_);
+    odes_.jacobian(x0, y0, li, dfdx_, dfdy_);
 
     for (label i=0; i<n_; i++)
     {
@@ -142,7 +143,7 @@ Foam::scalar Foam::rodas23::solve
         y[i] = y0[i] + dy_[i];
     }
 
-    odes_.derivatives(x0 + dx, y, dydx_);
+    odes_.derivatives(x0 + dx, y, li, dydx_);
 
     forAll(k3_, i)
     {
@@ -158,7 +159,7 @@ Foam::scalar Foam::rodas23::solve
         y[i] = y0[i] + dy_[i];
     }
 
-    odes_.derivatives(x0 + dx, y, dydx_);
+    odes_.derivatives(x0 + dx, y, li, dydx_);
 
     forAll(err_, i)
     {
@@ -180,10 +181,11 @@ void Foam::rodas23::solve
 (
     scalar& x,
     scalarField& y,
+    const label li,
     scalar& dxTry
 ) const
 {
-    adaptiveSolver::solve(odes_, x, y, dxTry);
+    adaptiveSolver::solve(odes_, x, y, li, dxTry);
 }
 
 

src/ODE/ODESolvers/rodas23/rodas23.H

@@ -113,6 +113,7 @@ public:
         (
             const scalar x0,
             const scalarField& y0,
+            const label li,
             const scalarField& dydx0,
             const scalar dx,
             scalarField& y
@@ -123,6 +124,7 @@ public:
         (
             scalar& x,
             scalarField& y,
+            const label li,
             scalar& dxTry
         ) const;
 };

src/ODE/ODESolvers/rodas34/rodas34.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2013-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2019 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -125,12 +125,13 @@ Foam::scalar Foam::rodas34::solve
 (
     const scalar x0,
     const scalarField& y0,
+    const label li,
     const scalarField& dydx0,
     const scalar dx,
     scalarField& y
 ) const
 {
-    odes_.jacobian(x0, y0, dfdx_, dfdy_);
+    odes_.jacobian(x0, y0, li, dfdx_, dfdy_);
 
     for (label i=0; i<n_; i++)
     {
@@ -158,7 +159,7 @@ Foam::scalar Foam::rodas34::solve
         y[i] = y0[i] + a21*k1_[i];
     }
 
-    odes_.derivatives(x0 + c2*dx, y, dydx_);
+    odes_.derivatives(x0 + c2*dx, y, li, dydx_);
 
     forAll(k2_, i)
     {
@@ -173,7 +174,7 @@ Foam::scalar Foam::rodas34::solve
         y[i] = y0[i] + a31*k1_[i] + a32*k2_[i];
     }
 
-    odes_.derivatives(x0 + c3*dx, y, dydx_);
+    odes_.derivatives(x0 + c3*dx, y, li, dydx_);
 
     forAll(k3_, i)
     {
@@ -188,7 +189,7 @@ Foam::scalar Foam::rodas34::solve
         y[i] = y0[i] + a41*k1_[i] + a42*k2_[i] + a43*k3_[i];
     }
 
-    odes_.derivatives(x0 + c4*dx, y, dydx_);
+    odes_.derivatives(x0 + c4*dx, y, li, dydx_);
 
     forAll(k4_, i)
     {
@@ -205,7 +206,7 @@ Foam::scalar Foam::rodas34::solve
         y[i] = y0[i] + dy_[i];
     }
 
-    odes_.derivatives(x0 + dx, y, dydx_);
+    odes_.derivatives(x0 + dx, y, li, dydx_);
 
     forAll(k5_, i)
     {
@@ -222,7 +223,7 @@ Foam::scalar Foam::rodas34::solve
         y[i] = y0[i] + dy_[i];
     }
 
-    odes_.derivatives(x0 + dx, y, dydx_);
+    odes_.derivatives(x0 + dx, y, li, dydx_);
 
     forAll(err_, i)
     {
@@ -245,10 +246,11 @@ void Foam::rodas34::solve
 (
     scalar& x,
     scalarField& y,
+    const label li,
     scalar& dxTry
 ) const
 {
-    adaptiveSolver::solve(odes_, x, y, dxTry);
+    adaptiveSolver::solve(odes_, x, y, li, dxTry);
 }
 
 

src/ODE/ODESolvers/rodas34/rodas34.H

@@ -117,6 +117,7 @@ public:
         (
             const scalar x0,
             const scalarField& y0,
+            const label li,
             const scalarField& dydx0,
             const scalar dx,
             scalarField& y
@@ -127,6 +128,7 @@ public:
         (
             scalar& x,
             scalarField& y,
+            const label li,
             scalar& dxTry
         ) const;
 };

src/ODE/ODESolvers/seulex/seulex.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2013-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2019 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -104,6 +104,7 @@ bool Foam::seulex::seul
 (
     const scalar x0,
     const scalarField& y0,
+    const label li,
     const scalar dxTot,
     const label k,
     scalarField& y,
@@ -126,7 +127,7 @@ bool Foam::seulex::seul
     LUDecompose(a_, pivotIndices_);
 
     scalar xnew = x0 + dx;
-    odes_.derivatives(xnew, y0, dy_);
+    odes_.derivatives(xnew, y0, li, dy_);
     LUBacksubstitute(a_, pivotIndices_, dy_);
 
     yTemp_ = y0;
@@ -145,7 +146,7 @@ bool Foam::seulex::seul
             }
             dy1 = sqrt(dy1);
 
-            odes_.derivatives(x0 + dx, yTemp_, dydx_);
+            odes_.derivatives(x0 + dx, yTemp_, li, dydx_);
             for (label i=0; i<n_; i++)
             {
                 dy_[i] = dydx_[i] - dy_[i]/dx;
@@ -181,7 +182,7 @@ bool Foam::seulex::seul
             }
         }
 
-        odes_.derivatives(xnew, yTemp_, dy_);
+        odes_.derivatives(xnew, yTemp_, li, dy_);
         LUBacksubstitute(a_, pivotIndices_, dy_);
     }
 
@@ -246,6 +247,7 @@ void Foam::seulex::solve
 (
     scalar& x,
     scalarField& y,
+    const label li,
     stepState& step
 ) const
 {
@@ -275,7 +277,7 @@ void Foam::seulex::solve
 
     if (theta_ > jacRedo_)
     {
-        odes_.jacobian(x, y, dfdx_, dfdy_);
+        odes_.jacobian(x, y, li, dfdx_, dfdy_);
         jacUpdated = true;
     }
 
@@ -299,7 +301,7 @@ void Foam::seulex::solve
 
         for (k=0; k<=kTarg_+1; k++)
         {
-            bool success = seul(x, y0_, dx, k, ySequence_, scale_);
+            bool success = seul(x, y0_, li, dx, k, ySequence_, scale_);
 
             if (!success)
             {
@@ -427,7 +429,7 @@ void Foam::seulex::solve
 
                 if (theta_ > jacRedo_ && !jacUpdated)
                 {
-                    odes_.jacobian(x, y, dfdx_, dfdy_);
+                    odes_.jacobian(x, y, li, dfdx_, dfdy_);
                     jacUpdated = true;
                 }
             }

src/ODE/ODESolvers/seulex/seulex.H

@@ -108,6 +108,7 @@ class seulex
         (
             const scalar x0,
             const scalarField& y0,
+            const label li,
             const scalar dxTot,
             const label k,
             scalarField& y,
@@ -150,6 +151,7 @@ public:
         (
             scalar& x,
             scalarField& y,
+            const label li,
             stepState& step
         ) const;
 };

src/ODE/ODESystem/ODESystem.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2019 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -67,19 +67,25 @@ public:
         virtual label nEqns() const = 0;
 
         //- Calculate the derivatives in dydx
+        //  for the current state x and y
+        //  and optional index into the list of systems to solve li
         virtual void derivatives
         (
             const scalar x,
             const scalarField& y,
+            const label li,
             scalarField& dydx
         ) const = 0;
 
         //- Calculate the Jacobian of the system
-        //  Need by the stiff-system solvers
+        //  for the current state x and y
+        //  and optional index into the list of systems to solve li.
+        //  Need by stiff-system solvers
         virtual void jacobian
         (
             const scalar x,
             const scalarField& y,
+            const label li,
             scalarField& dfdx,
             scalarSquareMatrix& dfdy
         ) const = 0;

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

@@ -103,9 +103,10 @@ Foam::StandardChemistryModel<ReactionThermo, ThermoType>::
 template<class ReactionThermo, class ThermoType>
 void Foam::StandardChemistryModel<ReactionThermo, ThermoType>::omega
 (
-    const scalarField& c,
-    const scalar T,
     const scalar p,
+    const scalar T,
+    const scalarField& c,
+    const label li,
     scalarField& dcdt
 ) const
 {
@@ -116,7 +117,7 @@ void Foam::StandardChemistryModel<ReactionThermo, ThermoType>::omega
     {
         const Reaction<ThermoType>& R = reactions_[i];
 
-        R.omega(p, T, c, dcdt);
+        R.omega(p, T, c, li, dcdt);
     }
 }
 
@@ -125,9 +126,10 @@ template<class ReactionThermo, class ThermoType>
 Foam::scalar Foam::StandardChemistryModel<ReactionThermo, ThermoType>::omegaI
 (
     const label index,
-    const scalarField& c,
-    const scalar T,
     const scalar p,
+    const scalar T,
+    const scalarField& c,
+    const label li,
     scalar& pf,
     scalar& cf,
     label& lRef,
@@ -137,7 +139,7 @@ Foam::scalar Foam::StandardChemistryModel<ReactionThermo, ThermoType>::omegaI
 ) const
 {
     const Reaction<ThermoType>& R = reactions_[index];
-    scalar w = R.omega(p, T, c, pf, cf, lRef, pr, cr, rRef);
+    scalar w = R.omega(p, T, c, li, pf, cf, lRef, pr, cr, rRef);
     return(w);
 }
 
@@ -147,6 +149,7 @@ void Foam::StandardChemistryModel<ReactionThermo, ThermoType>::derivatives
 (
     const scalar time,
     const scalarField& c,
+    const label li,
     scalarField& dcdt
 ) const
 {
@@ -158,7 +161,7 @@ void Foam::StandardChemistryModel<ReactionThermo, ThermoType>::derivatives
         c_[i] = max(c[i], 0);
     }
 
-    omega(c_, T, p, dcdt);
+    omega(p, T, c_, li, dcdt);
 
     // Constant pressure
     // dT/dt = ...
@@ -197,6 +200,7 @@ void Foam::StandardChemistryModel<ReactionThermo, ThermoType>::jacobian
 (
     const scalar t,
     const scalarField& c,
+    const label li,
     scalarField& dcdt,
     scalarSquareMatrix& J
 ) const
@@ -227,8 +231,8 @@ void Foam::StandardChemistryModel<ReactionThermo, ThermoType>::jacobian
     {
         const Reaction<ThermoType>& R = reactions_[ri];
         scalar kfwd, kbwd;
-        R.dwdc(p, T, c_, J, dcdt, omegaI, kfwd, kbwd, false, dummy);
+        R.dwdc(p, T, c_, li, J, dcdt, omegaI, kfwd, kbwd, false, dummy);
-        R.dwdT(p, T, c_, omegaI, kfwd, kbwd, J, false, dummy, nSpecie_);
+        R.dwdT(p, T, c_, li, omegaI, kfwd, kbwd, J, false, dummy, nSpecie_);
     }
 
     // The species derivatives of the temperature term are partially computed
@@ -318,7 +322,7 @@ Foam::StandardChemistryModel<ReactionThermo, ThermoType>::tc() const
             {
                 const Reaction<ThermoType>& R = reactions_[i];
 
-                R.omega(pi, Ti, c_, pf, cf, lRef, pr, cr, rRef);
+                R.omega(pi, Ti, c_, celli, pf, cf, lRef, pr, cr, rRef);
 
                 forAll(R.rhs(), s)
                 {
@@ -410,7 +414,10 @@ Foam::StandardChemistryModel<ReactionThermo, ThermoType>::calculateRR
         }
 
         const Reaction<ThermoType>& R = reactions_[ri];
-        const scalar omegai = R.omega(pi, Ti, c_, pf, cf, lRef, pr, cr, rRef);
+        const scalar omegai = R.omega
+        (
+            pi, Ti, c_, celli, pf, cf, lRef, pr, cr, rRef
+        );
 
         forAll(R.lhs(), s)
         {
@@ -461,7 +468,7 @@ void Foam::StandardChemistryModel<ReactionThermo, ThermoType>::calculate()
             c_[i] = rhoi*Yi/specieThermo_[i].W();
         }
 
-        omega(c_, Ti, pi, dcdt_);
+        omega(pi, Ti, c_, celli, dcdt_);
 
         for (label i=0; i<nSpecie_; i++)
         {
@@ -517,7 +524,7 @@ Foam::scalar Foam::StandardChemistryModel<ReactionThermo, ThermoType>::solve
             while (timeLeft > small)
             {
                 scalar dt = timeLeft;
-                this->solve(c_, Ti, pi, dt, this->deltaTChem_[celli]);
+                this->solve(pi, Ti, c_, celli, dt, this->deltaTChem_[celli]);
                 timeLeft -= dt;
             }
 

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

@@ -154,9 +154,10 @@ public:
         //- dc/dt = omega, rate of change in concentration, for each species
         virtual void omega
         (
-            const scalarField& c,
-            const scalar T,
             const scalar p,
+            const scalar T,
+            const scalarField& c,
+            const label li,
             scalarField& dcdt
         ) const;
 
@@ -166,9 +167,10 @@ public:
         virtual scalar omegaI
         (
             label iReaction,
-            const scalarField& c,
-            const scalar T,
             const scalar p,
+            const scalar T,
+            const scalarField& c,
+            const label li,
             scalar& pf,
             scalar& cf,
             label& lRef,
@@ -227,6 +229,7 @@ public:
             (
                 const scalar t,
                 const scalarField& c,
+                const label li,
                 scalarField& dcdt
             ) const;
 
@@ -234,15 +237,17 @@ public:
             (
                 const scalar t,
                 const scalarField& c,
+                const label li,
                 scalarField& dcdt,
                 scalarSquareMatrix& J
             ) const;
 
             virtual void solve
             (
-                scalarField &c,
-                scalar& T,
                 scalar& p,
+                scalar& T,
+                scalarField& c,
+                const label li,
                 scalar& deltaT,
                 scalar& subDeltaT
             ) const = 0;

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

@@ -149,9 +149,10 @@ Foam::TDACChemistryModel<ReactionThermo, ThermoType>::~TDACChemistryModel()
 template<class ReactionThermo, class ThermoType>
 void Foam::TDACChemistryModel<ReactionThermo, ThermoType>::omega
 (
-    const scalarField& c, // Contains all species even when mechRed is active
-    const scalar T,
     const scalar p,
+    const scalar T,
+    const scalarField& c, // Contains all species even when mechRed is active
+    const label li,
     scalarField& dcdt
 ) const
 {
@@ -170,7 +171,7 @@ void Foam::TDACChemistryModel<ReactionThermo, ThermoType>::omega
 
             scalar omegai = R.omega
             (
-                p, T, c, pf, cf, lRef, pr, cr, rRef
+                p, T, c, li, pf, cf, lRef, pr, cr, rRef
             );
 
             forAll(R.lhs(), s)
@@ -204,9 +205,10 @@ template<class ReactionThermo, class ThermoType>
 Foam::scalar Foam::TDACChemistryModel<ReactionThermo, ThermoType>::omega
 (
     const Reaction<ThermoType>& R,
-    const scalarField& c, // Contains all species even when mechRed is active
-    const scalar T,
     const scalar p,
+    const scalar T,
+    const scalarField& c, // Contains all species even when mechRed is active
+    const label li,
     scalar& pf,
     scalar& cf,
     label& lRef,
@@ -215,8 +217,8 @@ Foam::scalar Foam::TDACChemistryModel<ReactionThermo, ThermoType>::omega
     label& rRef
 ) const
 {
-    const scalar kf = R.kf(p, T, c);
+    const scalar kf = R.kf(p, T, c, li);
-    const scalar kr = R.kr(kf, p, T, c);
+    const scalar kr = R.kr(kf, p, T, c, li);
 
     const label Nl = R.lhs().size();
     const label Nr = R.rhs().size();
@@ -314,6 +316,7 @@ void Foam::TDACChemistryModel<ReactionThermo, ThermoType>::derivatives
 (
     const scalar time,
     const scalarField& c,
+    const label li,
     scalarField& dcdt
 ) const
 {
@@ -345,7 +348,7 @@ void Foam::TDACChemistryModel<ReactionThermo, ThermoType>::derivatives
         }
     }
 
-    omega(this->c_, T, p, dcdt);
+    omega(p, T, this->c_, li, dcdt);
 
     // Constant pressure
     // dT/dt = ...
@@ -398,6 +401,7 @@ void Foam::TDACChemistryModel<ReactionThermo, ThermoType>::jacobian
 (
     const scalar t,
     const scalarField& c,
+    const label li,
     scalarField& dcdt,
     scalarSquareMatrix& J
 ) const
@@ -451,6 +455,7 @@ void Foam::TDACChemistryModel<ReactionThermo, ThermoType>::jacobian
                 p,
                 T,
                 this->c_,
+                li,
                 J,
                 dcdt,
                 omegaI,
@@ -464,6 +469,7 @@ void Foam::TDACChemistryModel<ReactionThermo, ThermoType>::jacobian
                 p,
                 T,
                 this->c_,
+                li,
                 omegaI,
                 kfwd,
                 kbwd,
@@ -651,7 +657,7 @@ Foam::scalar Foam::TDACChemistryModel<ReactionThermo, ThermoType>::solve
             if (reduced)
             {
                 // Reduce mechanism change the number of species (only active)
-                mechRed_->reduceMechanism(c, Ti, pi);
+                mechRed_->reduceMechanism(pi, Ti, c, celli);
                 nActiveSpecies += mechRed_->NsSimp();
                 nAvg++;
                 scalar timeIncr = clockTime_.timeIncrement();
@@ -672,7 +678,12 @@ Foam::scalar Foam::TDACChemistryModel<ReactionThermo, ThermoType>::solve
                     // Solve the reduced set of ODE
                     this->solve
                     (
-                        simplifiedC_, Ti, pi, dt, this->deltaTChem_[celli]
+                        pi,
+                        Ti,
+                        simplifiedC_,
+                        celli,
+                        dt,
+                        this->deltaTChem_[celli]
                     );
 
                     for (label i=0; i<NsDAC_; i++)
@@ -682,7 +693,7 @@ Foam::scalar Foam::TDACChemistryModel<ReactionThermo, ThermoType>::solve
                 }
                 else
                 {
-                    this->solve(c, Ti, pi, dt, this->deltaTChem_[celli]);
+                    this->solve(pi, Ti, c, celli, dt, this->deltaTChem_[celli]);
                 }
                 timeLeft -= dt;
             }
@@ -713,7 +724,7 @@ Foam::scalar Foam::TDACChemistryModel<ReactionThermo, ThermoType>::solve
                     Rphiq[Rphiq.size()-1] = pi;
                 }
                 label growOrAdd =
-                    tabulation_->add(phiq, Rphiq, rhoi, deltaT[celli]);
+                    tabulation_->add(phiq, Rphiq, celli, rhoi, deltaT[celli]);
                 if (growOrAdd)
                 {
                     this->setTabulationResultsAdd(celli);

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

@@ -173,9 +173,10 @@ public:
         //- dc/dt = omega, rate of change in concentration, for each species
         virtual void omega
         (
-            const scalarField& c,
-            const scalar T,
             const scalar p,
+            const scalar T,
+            const scalarField& c,
+            const label li,
             scalarField& dcdt
         ) const;
 
@@ -184,9 +185,10 @@ public:
         virtual scalar omega
         (
             const Reaction<ThermoType>& r,
-            const scalarField& c,
-            const scalar T,
             const scalar p,
+            const scalar T,
+            const scalarField& c,
+            const label li,
             scalar& pf,
             scalar& cf,
             label& lRef,
@@ -215,6 +217,7 @@ public:
             (
                 const scalar t,
                 const scalarField& c,
+                const label li,
                 scalarField& dcdt
             ) const;
 
@@ -222,15 +225,17 @@ public:
             (
                 const scalar t,
                 const scalarField& c,
+                const label li,
                 scalarField& dcdt,
                 scalarSquareMatrix& J
             ) const;
 
             virtual void solve
             (
-                scalarField& c,
-                scalar& T,
                 scalar& p,
+                scalar& T,
+                scalarField& c,
+                const label li,
                 scalar& deltaT,
                 scalar& subDeltaT
             ) const = 0;

src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DAC/DAC.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2016-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2016-2019 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -270,9 +270,10 @@ Foam::chemistryReductionMethods::DAC<CompType, ThermoType>::~DAC()
 template<class CompType, class ThermoType>
 void Foam::chemistryReductionMethods::DAC<CompType, ThermoType>::reduceMechanism
 (
-    const scalarField &c,
+    const scalar p,
     const scalar T,
-    const scalar p
+    const scalarField& c,
+    const label li
 )
 {
     scalarField& completeC(this->chemistry_.completeC());
@@ -303,10 +304,11 @@ void Foam::chemistryReductionMethods::DAC<CompType, ThermoType>::reduceMechanism
     forAll(this->chemistry_.reactions(), i)
     {
         const Reaction<ThermoType>& R = this->chemistry_.reactions()[i];
+
         // for each reaction compute omegai
         scalar omegai = this->chemistry_.omega
         (
-            R, c1, T, p, pf, cf, lRef, pr, cr, rRef
+            R, p, T, c1, li, pf, cf, lRef, pr, cr, rRef
         );
 
         // Then for each pair of species composing this reaction,

src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DAC/DAC.H

@@ -147,9 +147,10 @@ public:
         //- Reduce the mechanism
         virtual void reduceMechanism
         (
-            const scalarField &c,
+            const scalar p,
             const scalar T,
-            const scalar p
+            const scalarField& c,
+            const label li
         );
 };
 

src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DRG/DRG.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2016-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2016-2019 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -69,9 +69,10 @@ Foam::chemistryReductionMethods::DRG<CompType, ThermoType>::~DRG()
 template<class CompType, class ThermoType>
 void Foam::chemistryReductionMethods::DRG<CompType, ThermoType>::reduceMechanism
 (
-    const scalarField &c,
+    const scalar p,
     const scalar T,
-    const scalar p
+    const scalarField& c,
+    const label li
 )
 {
     scalarField c1(this->nSpecie_+2, 0.0);
@@ -102,11 +103,12 @@ void Foam::chemistryReductionMethods::DRG<CompType, ThermoType>::reduceMechanism
     forAll(this->chemistry_.reactions(), i)
     {
         const Reaction<ThermoType>& R = this->chemistry_.reactions()[i];
+
         // For each reaction compute omegai
         scalar omegai = this->chemistry_.omega
         (
-         R, c1, T, p, pf, cf, lRef, pr, cr, rRef
+            R, p, T, c1, li, pf, cf, lRef, pr, cr, rRef
-         );
+        );
 
 
         // Then for each pair of species composing this reaction,

src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DRG/DRG.H

@@ -83,9 +83,10 @@ public:
         //- Reduce the mechanism
         virtual void reduceMechanism
         (
-            const scalarField &c,
+            const scalar p,
             const scalar T,
-            const scalar p
+            const scalarField& c,
+            const label li
         );
 };
 

src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DRGEP/DRGEP.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2016-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2016-2019 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -115,9 +115,10 @@ template<class CompType, class ThermoType>
 void Foam::chemistryReductionMethods::DRGEP<CompType, ThermoType>::
 reduceMechanism
 (
-    const scalarField &c,
+    const scalar p,
     const scalar T,
-    const scalar p
+    const scalarField& c,
+    const label li
 )
 {
     scalarField& completeC(this->chemistry_.completeC());
@@ -153,8 +154,8 @@ reduceMechanism
         // for each reaction compute omegai
         scalar omegai = this->chemistry_.omega
         (
-         R, c1, T, p, pf, cf, lRef, pr, cr, rRef
+            R, p,T, c1, li, pf, cf, lRef, pr, cr, rRef
-         );
+        );
         omegaV[i] = omegai;
 
         // then for each pair of species composing this reaction,

src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/DRGEP/DRGEP.H

@@ -154,9 +154,10 @@ public:
         //- Reduce the mechanism
         virtual void reduceMechanism
         (
-            const scalarField &c,
+            const scalar p,
             const scalar T,
-            const scalar p
+            const scalarField& c,
+            const label li
         );
 };
 

src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/EFA/EFA.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2016-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2016-2019 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -94,9 +94,10 @@ Foam::chemistryReductionMethods::EFA<CompType, ThermoType>::~EFA()
 template<class CompType, class ThermoType>
 void Foam::chemistryReductionMethods::EFA<CompType, ThermoType>::reduceMechanism
 (
-    const scalarField &c,
+    const scalar p,
     const scalar T,
-    const scalar p
+    const scalarField& c,
+    const label li
 )
 {
     scalarField& completeC(this->chemistry_.completeC());
@@ -132,10 +133,11 @@ void Foam::chemistryReductionMethods::EFA<CompType, ThermoType>::reduceMechanism
     forAll(this->chemistry_.reactions(), i)
     {
         const Reaction<ThermoType>& R = this->chemistry_.reactions()[i];
+
         // for each reaction compute omegai
         this->chemistry_.omega
         (
-            R, c1, T, p, pf, cf, lRef, pr, cr, rRef
+            R, p, T, c1, li, pf, cf, lRef, pr, cr, rRef
         );
         scalar fr = mag(pf*cf)+mag(pr*cr);
         scalar NCi(0.0),NHi(0.0),NOi(0.0),NNi(0.0);

src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/EFA/EFA.H

@@ -80,9 +80,10 @@ public:
         //- Reduce the mechanism
         virtual void reduceMechanism
         (
-            const scalarField &c,
+            const scalar p,
             const scalar T,
-            const scalar p
+            const scalarField& c,
+            const label li
         );
 };
 

src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/PFA/PFA.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2016-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2016-2019 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -72,9 +72,10 @@ Foam::chemistryReductionMethods::PFA<CompType, ThermoType>::~PFA()
 template<class CompType, class ThermoType>
 void Foam::chemistryReductionMethods::PFA<CompType, ThermoType>::reduceMechanism
 (
-    const scalarField &c,
+    const scalar p,
     const scalar T,
-    const scalar p
+    const scalarField& c,
+    const label li
 )
 {
     scalarField& completeC(this->chemistry_.completeC());
@@ -110,7 +111,7 @@ void Foam::chemistryReductionMethods::PFA<CompType, ThermoType>::reduceMechanism
         // for each reaction compute omegai
         scalar omegai = this->chemistry_.omega
         (
-            R, c1, T, p, pf, cf, lRef, pr, cr, rRef
+            R, p, T, c1, li, pf, cf, lRef, pr, cr, rRef
         );
 
         // then for each pair of species composing this reaction,

src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/PFA/PFA.H

@@ -82,9 +82,10 @@ public:
         //- Reduce the mechanism
         virtual void reduceMechanism
         (
-            const scalarField &c,
+            const scalar p,
             const scalar T,
-            const scalar p
+            const scalarField& c,
+            const label li
         );
 };
 

src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/chemistryReductionMethod/chemistryReductionMethod.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2016-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2016-2019 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -149,9 +149,10 @@ public:
         //- Reduce the mechanism
         virtual void reduceMechanism
         (
-            const scalarField &c,
+            const scalar p,
             const scalar T,
-            const scalar p
+            const scalarField& c,
+            const label li
         ) = 0;
 };
 

src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/noChemistryReduction/noChemistryReduction.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2016-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2016-2019 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -53,9 +53,10 @@ template<class CompType, class ThermoType>
 void Foam::chemistryReductionMethods::none<CompType, ThermoType>::
 reduceMechanism
 (
-    const scalarField &c,
+    const scalar p,
     const scalar T,
-    const scalar p
+    const scalarField& c,
+    const label li
 )
 {
     NotImplemented;

src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/reduction/noChemistryReduction/noChemistryReduction.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2016-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2016-2019 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -76,9 +76,10 @@ public:
         //- Reduce the mechanism
         virtual void reduceMechanism
         (
-            const scalarField &c,
+            const scalar p,
             const scalar T,
-            const scalar p
+            const scalarField& c,
+            const label li
         );
 };
 

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2016-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2016-2019 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -341,6 +341,7 @@ void Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::computeA
 (
     scalarSquareMatrix& A,
     const scalarField& Rphiq,
+    const label li,
     const scalar rhoi,
     const scalar dt
 )
@@ -374,7 +375,7 @@ void Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::computeA
     // A = C(psi0,t0)/(I-dt*J(psi(t0+dt)))
     // where C(psi0,t0) = I
     scalarField dcdt(speciesNumber + 2, Zero);
-    this->chemistry_.jacobian(runTime_.value(), Rcq, dcdt, A);
+    this->chemistry_.jacobian(runTime_.value(), Rcq, li, dcdt, A);
 
     // The jacobian is computed according to the molar concentration
     // the following conversion allows the code to use A with mass fraction
@@ -536,6 +537,7 @@ Foam::label Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::add
 (
     const scalarField& phiq,
     const scalarField& Rphiq,
+    const label li,
     const scalar rho,
     const scalar deltaT
 )
@@ -614,7 +616,7 @@ Foam::label Foam::chemistryTabulationMethods::ISAT<CompType, ThermoType>::add
     // Compute the A matrix needed to store the chemPoint.
     label ASize = this->chemistry_.nEqns() + nAdditionalEqns_ - 2;
     scalarSquareMatrix A(ASize, Zero);
-    computeA(A, Rphiq, rho, deltaT);
+    computeA(A, Rphiq, li, rho, deltaT);
 
     chemisTree().insertNewLeaf
     (

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

@@ -168,6 +168,7 @@ class ISAT
         (
             scalarSquareMatrix& A,
             const scalarField& Rphiq,
+            const label li,
             const scalar rho,
             const scalar dt
         );
@@ -232,6 +233,7 @@ public:
         (
             const scalarField& phiq,
             const scalarField& Rphiq,
+            const label li,
             const scalar rho,
             const scalar deltaT
         );

src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/chemistryTabulationMethod/chemistryTabulationMethod.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2016-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2016-2019 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -161,6 +161,7 @@ public:
         (
             const scalarField& phiQ,
             const scalarField& RphiQ,
+            const label li,
             const scalar rho,
             const scalar deltaT
         ) = 0;

src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/noChemistryTabulation/noChemistryTabulation.H

@@ -105,6 +105,7 @@ public:
         (
             const scalarField& phiq,
             const scalarField& Rphiq,
+            const label li,
             const scalar rho,
             const scalar deltaT
         )

src/thermophysicalModels/chemistryModel/chemistrySolver/EulerImplicit/EulerImplicit.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2019 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -91,9 +91,10 @@ void Foam::EulerImplicit<ChemistryModel>::updateRRInReactionI
 template<class ChemistryModel>
 void Foam::EulerImplicit<ChemistryModel>::solve
 (
-    scalarField& c,
-    scalar& T,
     scalar& p,
+    scalar& T,
+    scalarField& c,
+    const label li,
     scalar& deltaT,
     scalar& subDeltaT
 ) const
@@ -124,8 +125,10 @@ void Foam::EulerImplicit<ChemistryModel>::solve
         scalar pf, cf, pr, cr;
         label lRef, rRef;
 
-        const scalar omegai =
+        const scalar omegai
-            this->omegaI(i, c, T, p, pf, cf, lRef, pr, cr, rRef);
+        (
+            this->omegaI(i, p, T, c, li, pf, cf, lRef, pr, cr, rRef)
+        );
 
         scalar corr = 1;
         if (eqRateLimiter_)

src/thermophysicalModels/chemistryModel/chemistrySolver/EulerImplicit/EulerImplicit.H

@@ -107,9 +107,10 @@ public:
         //- Update the concentrations and return the chemical time
         virtual void solve
         (
-            scalarField& c,
-            scalar& T,
             scalar& p,
+            scalar& T,
+            scalarField& c,
+            const label li,
             scalar& deltaT,
             scalar& subDeltaT
         ) const;

src/thermophysicalModels/chemistryModel/chemistrySolver/chemistrySolver/chemistrySolver.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2019 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -70,9 +70,10 @@ public:
         //- Update the concentrations and return the chemical time
         virtual void solve
         (
-            scalarField &c,
-            scalar& T,
             scalar& p,
+            scalar& T,
+            scalarField& c,
+            const label li,
             scalar& deltaT,
             scalar& subDeltaT
         ) const = 0;

src/thermophysicalModels/chemistryModel/chemistrySolver/noChemistrySolver/noChemistrySolver.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2019 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -50,9 +50,10 @@ Foam::noChemistrySolver<ChemistryModel>::~noChemistrySolver()
 template<class ChemistryModel>
 void Foam::noChemistrySolver<ChemistryModel>::solve
 (
+    scalar&,
+    scalar&,
     scalarField&,
-    scalar&,
+    const label li,
-    scalar&,
     scalar&,
     scalar&
 ) const

src/thermophysicalModels/chemistryModel/chemistrySolver/noChemistrySolver/noChemistrySolver.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2019 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -74,9 +74,10 @@ public:
         //- Update the concentrations and return the chemical time
         virtual void solve
         (
-            scalarField& c,
-            scalar& T,
             scalar& p,
+            scalar& T,
+            scalarField& c,
+            const label li,
             scalar& deltaT,
             scalar& subDeltaT
         ) const;

src/thermophysicalModels/chemistryModel/chemistrySolver/ode/ode.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2019 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -49,9 +49,10 @@ Foam::ode<ChemistryModel>::~ode()
 template<class ChemistryModel>
 void Foam::ode<ChemistryModel>::solve
 (
-    scalarField& c,
-    scalar& T,
     scalar& p,
+    scalar& T,
+    scalarField& c,
+    const label li,
     scalar& deltaT,
     scalar& subDeltaT
 ) const
@@ -73,7 +74,7 @@ void Foam::ode<ChemistryModel>::solve
     cTp_[nSpecie] = T;
     cTp_[nSpecie+1] = p;
 
-    odeSolver_->solve(0, deltaT, cTp_, subDeltaT);
+    odeSolver_->solve(0, deltaT, cTp_, li, subDeltaT);
 
     for (int i=0; i<nSpecie; i++)
     {

src/thermophysicalModels/chemistryModel/chemistrySolver/ode/ode.H

@@ -83,9 +83,10 @@ public:
         //- Update the concentrations and return the chemical time
         virtual void solve
         (
-            scalarField& c,
-            scalar& T,
             scalar& p,
+            scalar& T,
+            scalarField& c,
+            const label li,
             scalar& deltaT,
             scalar& subDeltaT
         ) const;

src/thermophysicalModels/solidChemistryModel/pyrolysisChemistryModel/pyrolysisChemistryModel.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2013-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2019 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -48,8 +48,7 @@ pyrolysisChemistryModel
     nGases_(pyrolisisGases_.size()),
     nSpecie_(this->Ys_.size() + nGases_),
     RRg_(nGases_),
-    Ys0_(this->nSolids_),
+    Ys0_(this->nSolids_)
-    cellCounter_(0)
 {
     // create the fields for the chemistry sources
     forAll(this->RRs_, fieldi)
@@ -184,17 +183,16 @@ template<class CompType, class SolidThermo, class GasThermo>
 Foam::scalarField Foam::
 pyrolysisChemistryModel<CompType, SolidThermo, GasThermo>::omega
 (
-    const scalarField& c,
-    const scalar T,
     const scalar p,
+    const scalar T,
+    const scalarField& c,
+    const label li,
     const bool updateC0
 ) const
 {
     scalar pf, cf, pr, cr;
     label lRef, rRef;
 
-    const label celli = cellCounter_;
-
     scalarField om(nEqns(), 0.0);
 
     forAll(this->reactions_, i)
@@ -207,7 +205,7 @@ pyrolysisChemistryModel<CompType, SolidThermo, GasThermo>::omega
 
         scalar omegai = omega
         (
-            R, c, T, p, pf, cf, lRef, pr, cr, rRef
+            R, p, T, c, li, pf, cf, lRef, pr, cr, rRef
         );
         scalar rhoL = 0.0;
         forAll(R.lhs(), s)
@@ -226,7 +224,7 @@ pyrolysisChemistryModel<CompType, SolidThermo, GasThermo>::omega
 
             if (updateC0)
             {
-                Ys0_[si][celli] += sr*omegai;
+                Ys0_[si][li] += sr*omegai;
             }
         }
         forAll(R.grhs(), g)
@@ -245,9 +243,10 @@ Foam::scalar
 Foam::pyrolysisChemistryModel<CompType, SolidThermo, GasThermo>::omega
 (
     const Reaction<SolidThermo>& R,
-    const scalarField& c,
-    const scalar T,
     const scalar p,
+    const scalar T,
+    const scalarField& c,
+    const label li,
     scalar& pf,
     scalar& cf,
     label& lRef,
@@ -258,14 +257,12 @@ Foam::pyrolysisChemistryModel<CompType, SolidThermo, GasThermo>::omega
 {
     scalarField c1(nSpecie_, 0.0);
 
-    label celli = cellCounter_;
-
     for (label i=0; i<nSpecie_; i++)
     {
         c1[i] = max(0.0, c[i]);
     }
 
-    scalar kf = R.kf(p, T, c1);
+    scalar kf = R.kf(p, T, c1, li);
 
     const label Nl = R.lhs().size();
 
@@ -275,8 +272,8 @@ Foam::pyrolysisChemistryModel<CompType, SolidThermo, GasThermo>::omega
         const scalar exp = R.lhs()[si].exponent;
 
         kf *=
-            pow(c1[si]/Ys0_[si][celli], exp)
+            pow(c1[si]/Ys0_[si][li], exp)
-           *(Ys0_[si][celli]);
+           *(Ys0_[si][li]);
     }
 
     return kf;
@@ -288,9 +285,10 @@ Foam::scalar Foam::pyrolysisChemistryModel<CompType, SolidThermo, GasThermo>::
 omegaI
 (
     const label index,
-    const scalarField& c,
-    const scalar T,
     const scalar p,
+    const scalar T,
+    const scalarField& c,
+    const label li,
     scalar& pf,
     scalar& cf,
     label& lRef,
@@ -301,7 +299,7 @@ omegaI
 {
 
     const Reaction<SolidThermo>& R = this->reactions_[index];
-    scalar w = omega(R, c, T, p, pf, cf, lRef, pr, cr, rRef);
+    scalar w = omega(R, p, T, c, li, pf, cf, lRef, pr, cr, rRef);
     return(w);
 }
 
@@ -311,7 +309,8 @@ void Foam::pyrolysisChemistryModel<CompType, SolidThermo, GasThermo>::
 derivatives
 (
     const scalar time,
-    const scalarField &c,
+    const scalarField& c,
+    const label li,
     scalarField& dcdt
 ) const
 {
@@ -320,7 +319,7 @@ derivatives
 
     dcdt = 0.0;
 
-    dcdt = omega(c, T, p);
+    dcdt = omega(p, T, c, li);
 
     // Total mass concentration
     scalar cTot = 0.0;
@@ -354,6 +353,7 @@ jacobian
 (
     const scalar t,
     const scalarField& c,
+    const label li,
     scalarField& dcdt,
     scalarSquareMatrix& dfdc
 ) const
@@ -377,13 +377,13 @@ jacobian
     }
 
     // length of the first argument must be nSolids
-    dcdt = omega(c2, T, p);
+    dcdt = omega(p, T, c2, li);
 
     for (label ri=0; ri<this->reactions_.size(); ri++)
     {
         const Reaction<SolidThermo>& R = this->reactions_[ri];
 
-        scalar kf0 = R.kf(p, T, c2);
+        scalar kf0 = R.kf(p, T, c2, li);
 
         forAll(R.lhs(), j)
         {
@@ -434,14 +434,13 @@ jacobian
 
     // calculate the dcdT elements numerically
     scalar delta = 1.0e-8;
-    scalarField dcdT0 = omega(c2, T - delta, p);
+    scalarField dcdT0 = omega(p , T - delta, c2, li);
-    scalarField dcdT1 = omega(c2, T + delta, p);
+    scalarField dcdT1 = omega(p, T + delta, c2, li);
 
     for (label i=0; i<nEqns(); i++)
     {
         dfdc[i][nSpecie_] = 0.5*(dcdT1[i] - dcdT0[i])/delta;
     }
-
 }
 
 
@@ -489,8 +488,6 @@ calculate()
 
     forAll(rho, celli)
     {
-        cellCounter_ = celli;
-
         const scalar delta = this->mesh().V()[celli];
 
         if (this->reactingCells_[celli])
@@ -505,7 +502,7 @@ calculate()
                 c[i] = rhoi*this->Ys_[i][celli]*delta;
             }
 
-            const scalarField dcdt = omega(c, Ti, pi, true);
+            const scalarField dcdt = omega(pi, Ti, c, celli, true);
 
             forAll(this->RRs_, i)
             {
@@ -570,8 +567,6 @@ Foam::pyrolysisChemistryModel<CompType, SolidThermo, GasThermo>::solve
     {
         if (this->reactingCells_[celli])
         {
-            cellCounter_ = celli;
-
             scalar rhoi = rho[celli];
             scalar pi = p[celli];
             scalar Ti = T[celli];
@@ -590,7 +585,7 @@ Foam::pyrolysisChemistryModel<CompType, SolidThermo, GasThermo>::solve
             while (timeLeft > small)
             {
                 scalar dt = timeLeft;
-                this->solve(c, Ti, pi, dt, this->deltaTChem_[celli]);
+                this->solve(pi, Ti, c, celli, dt, this->deltaTChem_[celli]);
                 timeLeft -= dt;
             }
 
@@ -608,7 +603,7 @@ Foam::pyrolysisChemistryModel<CompType, SolidThermo, GasThermo>::solve
             }
 
             // Update Ys0_
-            dc = omega(c0, Ti, pi, true);
+            dc = omega(pi, Ti, c0, celli, true);
         }
     }
 
@@ -654,9 +649,10 @@ Foam::pyrolysisChemistryModel<CompType, SolidThermo, GasThermo>::gasHs
 template<class CompType, class SolidThermo, class GasThermo>
 void Foam::pyrolysisChemistryModel<CompType, SolidThermo, GasThermo>::solve
 (
-    scalarField &c,
-    scalar& T,
     scalar& p,
+    scalar& T,
+    scalarField& c,
+    const label li,
     scalar& deltaT,
     scalar& subDeltaT
 ) const

src/thermophysicalModels/solidChemistryModel/pyrolysisChemistryModel/pyrolysisChemistryModel.H

@@ -87,9 +87,6 @@ private:
         //- List of accumulative solid concentrations
         mutable PtrList<volScalarField> Ys0_;
 
-        //- Cell counter
-        label cellCounter_;
-
 
 public:
 
@@ -125,9 +122,10 @@ public:
         //- dc/dt = omega, rate of change in concentration, for each species
         virtual scalarField omega
         (
-            const scalarField& c,
-            const scalar T,
             const scalar p,
+            const scalar T,
+            const scalarField& c,
+            const label li,
             const bool updateC0 = false
         ) const;
 
@@ -137,9 +135,10 @@ public:
         virtual scalar omega
         (
             const Reaction<SolidThermo>& r,
-            const scalarField& c,
-            const scalar T,
             const scalar p,
+            const scalar T,
+            const scalarField& c,
+            const label li,
             scalar& pf,
             scalar& cf,
             label& lRef,
@@ -155,9 +154,10 @@ public:
         virtual scalar omegaI
         (
             label iReaction,
-            const scalarField& c,
-            const scalar T,
             const scalar p,
+            const scalar T,
+            const scalarField& c,
+            const label li,
             scalar& pf,
             scalar& cf,
             label& lRef,
@@ -204,6 +204,7 @@ public:
             (
                 const scalar t,
                 const scalarField& c,
+                const label li,
                 scalarField& dcdt
             ) const;
 
@@ -211,15 +212,17 @@ public:
             (
                 const scalar t,
                 const scalarField& c,
+                const label li,
                 scalarField& dcdt,
                 scalarSquareMatrix& dfdc
             ) const;
 
             virtual void solve
             (
-                scalarField &c,
-                scalar& T,
                 scalar& p,
+                scalar& T,
+                scalarField& c,
+                const label li,
                 scalar& deltaT,
                 scalar& subDeltaT
             ) const;

src/thermophysicalModels/solidChemistryModel/solidChemistryModel/solidChemistryModel.H

@@ -126,9 +126,10 @@ public:
         //- dc/dt = omega, rate of change in concentration, for each species
         virtual scalarField omega
         (
-            const scalarField& c,
-            const scalar T,
             const scalar p,
+            const scalar T,
+            const scalarField& c,
+            const label li,
             const bool updateC0 = false
         ) const = 0;
 
@@ -137,9 +138,10 @@ public:
         virtual scalar omega
         (
             const Reaction<SolidThermo>& r,
-            const scalarField& c,
-            const scalar T,
             const scalar p,
+            const scalar T,
+            const scalarField& c,
+            const label li,
             scalar& pf,
             scalar& cf,
             label& lRef,
@@ -154,9 +156,10 @@ public:
         virtual scalar omegaI
         (
             label iReaction,
-            const scalarField& c,
-            const scalar T,
             const scalar p,
+            const scalar T,
+            const scalarField& c,
+            const label li,
             scalar& pf,
             scalar& cf,
             label& lRef,
@@ -204,6 +207,7 @@ public:
             (
                 const scalar t,
                 const scalarField& c,
+                const label li,
                 scalarField& dcdt
             ) const = 0;
 
@@ -211,15 +215,17 @@ public:
             (
                 const scalar t,
                 const scalarField& c,
+                const label li,
                 scalarField& dcdt,
                 scalarSquareMatrix& dfdc
             ) const = 0;
 
             virtual void solve
             (
-                scalarField &c,
-                scalar& T,
                 scalar& p,
+                scalar& T,
+                scalarField& c,
+                const label li,
                 scalar& deltaT,
                 scalar& subDeltaT
             ) const = 0;

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

@@ -101,14 +101,16 @@ public:
         (
             const scalar p,
             const scalar T,
-            const scalarField& c
+            const scalarField& c,
+            const label li
         ) const;
 
         inline scalar ddT
         (
             const scalar p,
             const scalar T,
-            const scalarField& c
+            const scalarField& c,
+            const label li
         ) const;
 
         //- Third-body efficiencies (beta = 1-alpha)
@@ -123,6 +125,7 @@ public:
             const scalar p,
             const scalar T,
             const scalarField& c,
+            const label li,
             scalarField& dcidc
         ) const;
 
@@ -132,7 +135,8 @@ public:
         (
             const scalar p,
             const scalar T,
-            const scalarField& c
+            const scalarField& c,
+            const label li
         ) const;
 
 

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

@@ -57,7 +57,8 @@ inline Foam::scalar Foam::solidArrheniusReactionRate::operator()
 (
     const scalar,
     const scalar T,
-    const scalarField&
+    const scalarField&,
+    const label li
 ) const
 {
     if (T < Tcrit_)
@@ -75,7 +76,8 @@ inline Foam::scalar Foam::solidArrheniusReactionRate::ddT
 (
     const scalar p,
     const scalar T,
-    const scalarField&
+    const scalarField&,
+    const label li
 ) const
 {
     if (T < Tcrit_)
@@ -101,6 +103,7 @@ inline void Foam::solidArrheniusReactionRate::dcidc
     const scalar p,
     const scalar T,
     const scalarField& c,
+    const label li,
     scalarField& dcidc
 ) const
 {}
@@ -110,7 +113,8 @@ inline Foam::scalar Foam::solidArrheniusReactionRate::dcidT
 (
     const scalar p,
     const scalar T,
-    const scalarField& c
+    const scalarField& c,
+    const label li
 ) 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     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2019 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -99,10 +99,11 @@ Foam::scalar Foam::IrreversibleReaction
 (
     const scalar p,
     const scalar T,
-    const scalarField& c
+    const scalarField& c,
+    const label li
 ) const
 {
-    return k_(p, T, c);
+    return k_(p, T, c, li);
 }
 
 
@@ -122,7 +123,8 @@ Foam::scalar Foam::IrreversibleReaction
     const scalar kfwd,
     const scalar p,
     const scalar T,
-    const scalarField& c
+    const scalarField& c,
+    const label li
 ) const
 {
     return 0;
@@ -144,7 +146,8 @@ Foam::scalar Foam::IrreversibleReaction
 (
     const scalar p,
     const scalar T,
-    const scalarField& c
+    const scalarField& c,
+    const label li
 ) const
 {
     return 0;
@@ -166,10 +169,11 @@ Foam::scalar Foam::IrreversibleReaction
 (
     const scalar p,
     const scalar T,
-    const scalarField& c
+    const scalarField& c,
+    const label li
 ) const
 {
-    return k_.ddT(p, T, c);
+    return k_.ddT(p, T, c, li);
 }
 
 
@@ -189,6 +193,7 @@ Foam::scalar Foam::IrreversibleReaction
     const scalar p,
     const scalar T,
     const scalarField& c,
+    const label li,
     const scalar dkfdT,
     const scalar kr
 ) const
@@ -231,10 +236,11 @@ void Foam::IrreversibleReaction
     const scalar p,
     const scalar T,
     const scalarField& c,
+    const label li,
     scalarField& dcidc
 ) const
 {
-    k_.dcidc(p, T, c, dcidc);
+    k_.dcidc(p, T, c, li, dcidc);
 }
 
 
@@ -253,10 +259,11 @@ Foam::scalar Foam::IrreversibleReaction
 (
     const scalar p,
     const scalar T,
-    const scalarField& c
+    const scalarField& c,
+    const label li
 ) const
 {
-    return k_.dcidT(p, T, c);
+    return k_.dcidT(p, T, c, li);
 }
 
 

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

@@ -146,7 +146,8 @@ public:
             (
                 const scalar p,
                 const scalar T,
-                const scalarField& c
+                const scalarField& c,
+                const label li
             ) const;
 
             //- Reverse rate constant from the given forward rate constant
@@ -156,7 +157,8 @@ public:
                 const scalar kfwd,
                 const scalar p,
                 const scalar T,
-                const scalarField& c
+                const scalarField& c,
+                const label li
             ) const;
 
             //- Reverse rate constant
@@ -165,7 +167,8 @@ public:
             (
                 const scalar p,
                 const scalar T,
-                const scalarField& c
+                const scalarField& c,
+                const label li
             ) const;
 
 
@@ -176,7 +179,8 @@ public:
             (
                 const scalar p,
                 const scalar T,
-                const scalarField& c
+                const scalarField& c,
+                const label li
             ) const;
 
             //- Temperature derivative of reverse rate
@@ -186,6 +190,7 @@ public:
                 const scalar p,
                 const scalar T,
                 const scalarField& c,
+                const label li,
                 const scalar dkfdT,
                 const scalar kr
             ) const;
@@ -202,6 +207,7 @@ public:
                 const scalar p,
                 const scalar T,
                 const scalarField& c,
+                const label li,
                 scalarField& dcidc
             ) const;
 
@@ -211,7 +217,8 @@ public:
             (
                 const scalar p,
                 const scalar T,
-                const scalarField& c
+                const scalarField& c,
+                const label li
             ) const;
 
 

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2019 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -124,10 +124,11 @@ Foam::NonEquilibriumReversibleReaction
 (
     const scalar p,
     const scalar T,
-    const scalarField& c
+    const scalarField& c,
+    const label li
 ) const
 {
-    return fk_(p, T, c);
+    return fk_(p, T, c, li);
 }
 
 
@@ -148,10 +149,11 @@ Foam::NonEquilibriumReversibleReaction
     const scalar,
     const scalar p,
     const scalar T,
-    const scalarField& c
+    const scalarField& c,
+    const label li
 ) const
 {
-    return rk_(p, T, c);
+    return rk_(p, T, c, li);
 }
 
 
@@ -171,10 +173,11 @@ Foam::NonEquilibriumReversibleReaction
 (
     const scalar p,
     const scalar T,
-    const scalarField& c
+    const scalarField& c,
+    const label li
 ) const
 {
-    return rk_(p, T, c);
+    return rk_(p, T, c, li);
 }
 
 
@@ -194,10 +197,11 @@ Foam::NonEquilibriumReversibleReaction
 (
     const scalar p,
     const scalar T,
-    const scalarField& c
+    const scalarField& c,
+    const label li
 ) const
 {
-    return fk_.ddT(p, T, c);
+    return fk_.ddT(p, T, c, li);
 }
 
 
@@ -218,11 +222,12 @@ Foam::NonEquilibriumReversibleReaction
     const scalar p,
     const scalar T,
     const scalarField& c,
+    const label li,
     const scalar dkfdT,
     const scalar kr
 ) const
 {
-    return rk_.ddT(p, T, c);
+    return rk_.ddT(p, T, c, li);
 }
 
 
@@ -260,10 +265,11 @@ void Foam::NonEquilibriumReversibleReaction
     const scalar p,
     const scalar T,
     const scalarField& c,
+    const label li,
     scalarField& dcidc
 ) const
 {
-    fk_.dcidc(p, T, c, dcidc);
+    fk_.dcidc(p, T, c, li, dcidc);
 }
 
 
@@ -282,10 +288,11 @@ Foam::scalar Foam::NonEquilibriumReversibleReaction
 (
     const scalar p,
     const scalar T,
-    const scalarField& c
+    const scalarField& c,
+    const label li
 ) const
 {
-    return fk_.dcidT(p, T, c);
+    return fk_.dcidT(p, T, c, li);
 }
 
 

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

@@ -134,7 +134,8 @@ public:
             (
                 const scalar p,
                 const scalar T,
-                const scalarField& c
+                const scalarField& c,
+                const label li
             ) const;
 
             //- Reverse rate constant from the given formard rate constant
@@ -143,7 +144,8 @@ public:
                 const scalar kfwd,
                 const scalar p,
                 const scalar T,
-                const scalarField& c
+                const scalarField& c,
+                const label li
             ) const;
 
             //- Reverse rate constant.
@@ -153,7 +155,8 @@ public:
             (
                 const scalar p,
                 const scalar T,
-                const scalarField& c
+                const scalarField& c,
+                const label li
             ) const;
 
 
@@ -164,7 +167,8 @@ public:
             (
                 const scalar p,
                 const scalar T,
-                const scalarField& c
+                const scalarField& c,
+                const label li
             ) const;
 
             //- Temperature derivative of backward rate
@@ -173,6 +177,7 @@ public:
                 const scalar p,
                 const scalar T,
                 const scalarField& c,
+                const label li,
                 const scalar dkfdT,
                 const scalar kr
             ) const;
@@ -189,6 +194,7 @@ public:
                 const scalar p,
                 const scalar T,
                 const scalarField& c,
+                const label li,
                 scalarField& dcidc
             ) const;
 
@@ -198,7 +204,8 @@ public:
             (
                 const scalar p,
                 const scalar T,
-                const scalarField& c
+                const scalarField& c,
+                const label li
             ) const;
 
 

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

@@ -205,6 +205,7 @@ void Foam::Reaction<ReactionThermo>::ddot
     const scalar p,
     const scalar T,
     const scalarField& c,
+    const label li,
     scalarField& d
 ) const
 {
@@ -217,6 +218,7 @@ void Foam::Reaction<ReactionThermo>::fdot
     const scalar p,
     const scalar T,
     const scalarField& c,
+    const label li,
     scalarField& f
 ) const
 {
@@ -229,6 +231,7 @@ void Foam::Reaction<ReactionThermo>::omega
     const scalar p,
     const scalar T,
     const scalarField& c,
+    const label li,
     scalarField& dcdt
 ) const
 {
@@ -237,7 +240,7 @@ void Foam::Reaction<ReactionThermo>::omega
 
     scalar omegaI = omega
     (
-        p, T, c, pf, cf, lRef, pr, cr, rRef
+        p, T, c, li, pf, cf, lRef, pr, cr, rRef
     );
 
     forAll(lhs_, i)
@@ -261,6 +264,7 @@ Foam::scalar Foam::Reaction<ReactionThermo>::omega
     const scalar p,
     const scalar T,
     const scalarField& c,
+    const label li,
     scalar& pf,
     scalar& cf,
     label& lRef,
@@ -272,8 +276,8 @@ Foam::scalar Foam::Reaction<ReactionThermo>::omega
 
     scalar clippedT = min(max(T, this->Tlow()), this->Thigh());
 
-    const scalar kf = this->kf(p, clippedT, c);
+    const scalar kf = this->kf(p, clippedT, c, li);
-    const scalar kr = this->kr(kf, p, clippedT, c);
+    const scalar kr = this->kr(kf, p, clippedT, c, li);
 
     pf = 1;
     pr = 1;
@@ -375,6 +379,7 @@ void Foam::Reaction<ReactionThermo>::dwdc
     const scalar p,
     const scalar T,
     const scalarField& c,
+    const label li,
     scalarSquareMatrix& J,
     scalarField& dcdt,
     scalar& omegaI,
@@ -387,7 +392,7 @@ void Foam::Reaction<ReactionThermo>::dwdc
     scalar pf, cf, pr, cr;
     label lRef, rRef;
 
-    omegaI = omega(p, T, c, pf, cf, lRef, pr, cr, rRef);
+    omegaI = omega(p, T, c, li, pf, cf, lRef, pr, cr, rRef);
 
     forAll(lhs_, i)
     {
@@ -402,8 +407,8 @@ void Foam::Reaction<ReactionThermo>::dwdc
         dcdt[si] += sr*omegaI;
     }
 
-    kfwd = this->kf(p, T, c);
+    kfwd = this->kf(p, T, c, li);
-    kbwd = this->kr(kfwd, p, T, c);
+    kbwd = this->kr(kfwd, p, T, c, li);
 
     forAll(lhs_, j)
     {
@@ -504,7 +509,7 @@ void Foam::Reaction<ReactionThermo>::dwdc
     {
         // This temporary array needs to be cached for efficiency
         scalarField dcidc(beta.size());
-        this->dcidc(p, T, c, dcidc);
+        this->dcidc(p, T, c, li, dcidc);
 
         forAll(beta, j)
         {
@@ -538,6 +543,7 @@ void Foam::Reaction<ReactionThermo>::dwdT
     const scalar p,
     const scalar T,
     const scalarField& c,
+    const label li,
     const scalar omegaI,
     const scalar kfwd,
     const scalar kbwd,
@@ -550,8 +556,8 @@ void Foam::Reaction<ReactionThermo>::dwdT
     scalar kf = kfwd;
     scalar kr = kbwd;
 
-    scalar dkfdT = this->dkfdT(p, T, c);
+    scalar dkfdT = this->dkfdT(p, T, c, li);
-    scalar dkrdT = this->dkrdT(p, T, c, dkfdT, kr);
+    scalar dkrdT = this->dkrdT(p, T, c, li, dkfdT, kr);
 
     scalar sumExp = 0.0;
     forAll(lhs_, i)
@@ -582,7 +588,7 @@ void Foam::Reaction<ReactionThermo>::dwdT
 
     // For reactions including third-body efficiencies or pressure dependent
     // reaction, an additional term is needed
-    scalar dcidT = this->dcidT(p, T, c);
+    scalar dcidT = this->dcidT(p, T, c, li);
     dcidT *= omegaI;
 
     // J(i, indexT) = sum_reactions nu_i dqdT

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

@@ -210,6 +210,7 @@ public:
                 const scalar p,
                 const scalar T,
                 const scalarField& c,
+                const label li,
                 scalarField& d
             ) const;
 
@@ -219,6 +220,7 @@ public:
                 const scalar p,
                 const scalar T,
                 const scalarField& c,
+                const label li,
                 scalarField& f
             ) const;
 
@@ -228,6 +230,7 @@ public:
                 const scalar p,
                 const scalar T,
                 const scalarField& c,
+                const label li,
                 scalarField& dcdt
             ) const;
 
@@ -237,6 +240,7 @@ public:
                 const scalar p,
                 const scalar T,
                 const scalarField& c,
+                const label li,
                 scalar& pf,
                 scalar& cf,
                 label& lRef,
@@ -252,7 +256,8 @@ public:
             (
                 const scalar p,
                 const scalar T,
-                const scalarField& c
+                const scalarField& c,
+                const label li
             ) const = 0;
 
             //- Reverse rate constant from the given forward rate constant
@@ -261,7 +266,8 @@ public:
                 const scalar kfwd,
                 const scalar p,
                 const scalar T,
-                const scalarField& c
+                const scalarField& c,
+                const label li
             ) const = 0;
 
             //- Reverse rate constant
@@ -269,7 +275,8 @@ public:
             (
                 const scalar p,
                 const scalar T,
-                const scalarField& c
+                const scalarField& c,
+                const label li
             ) const = 0;
 
 
@@ -282,6 +289,7 @@ public:
                 const scalar p,
                 const scalar T,
                 const scalarField& c,
+                const label li,
                 scalarSquareMatrix& J,
                 scalarField& dcdt,
                 scalar& omegaI,
@@ -298,6 +306,7 @@ public:
                 const scalar p,
                 const scalar T,
                 const scalarField& c,
+                const label li,
                 const scalar omegaI,
                 const scalar kfwd,
                 const scalar kbwd,
@@ -312,7 +321,8 @@ public:
             (
                 const scalar p,
                 const scalar T,
-                const scalarField& c
+                const scalarField& c,
+                const label li
             ) const = 0;
 
             //- Temperature derivative of reverse rate
@@ -321,6 +331,7 @@ public:
                 const scalar p,
                 const scalar T,
                 const scalarField& c,
+                const label li,
                 const scalar dkfdT,
                 const scalar kr
             ) const = 0;
@@ -336,6 +347,7 @@ public:
                 const scalar p,
                 const scalar T,
                 const scalarField& c,
+                const label li,
                 scalarField& dcidc
             ) const = 0;
 
@@ -344,7 +356,8 @@ public:
             (
                 const scalar p,
                 const scalar T,
-                const scalarField& c
+                const scalarField& c,
+                const label li
             ) const = 0;
 
 

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

@@ -106,7 +106,8 @@ Foam::scalar Foam::ReactionProxy<ReactionThermo>::kf
 (
     const scalar p,
     const scalar T,
-    const scalarField& c
+    const scalarField& c,
+    const label li
 ) const
 {
     NotImplemented;
@@ -120,7 +121,8 @@ Foam::scalar Foam::ReactionProxy<ReactionThermo>::kr
     const scalar kfwd,
     const scalar p,
     const scalar T,
-    const scalarField& c
+    const scalarField& c,
+    const label li
 ) const
 {
     NotImplemented;
@@ -133,7 +135,8 @@ Foam::scalar Foam::ReactionProxy<ReactionThermo>::kr
 (
     const scalar p,
     const scalar T,
-    const scalarField& c
+    const scalarField& c,
+    const label li
 ) const
 {
     NotImplemented;
@@ -146,7 +149,8 @@ Foam::scalar Foam::ReactionProxy<ReactionThermo>::dkfdT
 (
     const scalar p,
     const scalar T,
-    const scalarField& c
+    const scalarField& c,
+    const label li
 ) const
 {
     NotImplemented;
@@ -160,6 +164,7 @@ Foam::scalar Foam::ReactionProxy<ReactionThermo>::dkrdT
     const scalar p,
     const scalar T,
     const scalarField& c,
+    const label li,
     const scalar dkfdT,
     const scalar kr
 ) const
@@ -184,6 +189,7 @@ void Foam::ReactionProxy<ReactionThermo>::dcidc
     const scalar p,
     const scalar T,
     const scalarField& c,
+    const label li,
     scalarField& dcidc
 ) const
 {
@@ -196,7 +202,8 @@ Foam::scalar Foam::ReactionProxy<ReactionThermo>::dcidT
 (
     const scalar p,
     const scalar T,
-    const scalarField& c
+    const scalarField& c,
+    const label li
 ) const
 {
     NotImplemented;

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

@@ -107,7 +107,8 @@ public:
             (
                 const scalar p,
                 const scalar T,
-                const scalarField& c
+                const scalarField& c,
+                const label li
             ) const;
 
             //- Reverse rate constant from the given forward rate constant
@@ -116,7 +117,8 @@ public:
                 const scalar kfwd,
                 const scalar p,
                 const scalar T,
-                const scalarField& c
+                const scalarField& c,
+                const label li
             ) const;
 
             //- Reverse rate constant
@@ -124,7 +126,8 @@ public:
             (
                 const scalar p,
                 const scalar T,
-                const scalarField& c
+                const scalarField& c,
+                const label li
             ) const;
 
 
@@ -135,7 +138,8 @@ public:
             (
                 const scalar p,
                 const scalar T,
-                const scalarField& c
+                const scalarField& c,
+                const label li
             ) const;
 
             //- Temperature derivative of reverse rate
@@ -144,6 +148,7 @@ public:
                 const scalar p,
                 const scalar T,
                 const scalarField& c,
+                const label li,
                 const scalar dkfdT,
                 const scalar kr
             ) const;
@@ -159,6 +164,7 @@ public:
                 const scalar p,
                 const scalar T,
                 const scalarField& c,
+                const label li,
                 scalarField& dcidc
             ) const;
 
@@ -167,7 +173,8 @@ public:
             (
                 const scalar p,
                 const scalar T,
-                const scalarField& c
+                const scalarField& c,
+                const label li
             ) const;
 };
 

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2019 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -99,10 +99,11 @@ Foam::scalar Foam::ReversibleReaction
 (
     const scalar p,
     const scalar T,
-    const scalarField& c
+    const scalarField& c,
+    const label li
 ) const
 {
-    return k_(p, T, c);
+    return k_(p, T, c, li);
 }
 
 
@@ -122,7 +123,8 @@ Foam::scalar Foam::ReversibleReaction
     const scalar kfwd,
     const scalar p,
     const scalar T,
-    const scalarField& c
+    const scalarField& c,
+    const label li
 ) const
 {
     return kfwd/max(this->Kc(p, T), rootSmall);
@@ -144,10 +146,11 @@ Foam::scalar Foam::ReversibleReaction
 (
     const scalar p,
     const scalar T,
-    const scalarField& c
+    const scalarField& c,
+    const label li
 ) const
 {
-    return kr(kf(p, T, c), p, T, c);
+    return kr(kf(p, T, c, li), p, T, c, li);
 }
 
 
@@ -166,10 +169,11 @@ Foam::scalar Foam::ReversibleReaction
 (
     const scalar p,
     const scalar T,
-    const scalarField& c
+    const scalarField& c,
+    const label li
 ) const
 {
-    return k_.ddT(p, T, c);
+    return k_.ddT(p, T, c, li);
 }
 
 
@@ -189,6 +193,7 @@ Foam::scalar Foam::ReversibleReaction
     const scalar p,
     const scalar T,
     const scalarField& c,
+    const label li,
     const scalar dkfdT,
     const scalar kr
 ) const
@@ -233,10 +238,11 @@ void Foam::ReversibleReaction
     const scalar p,
     const scalar T,
     const scalarField& c,
+    const label li,
     scalarField& dcidc
 ) const
 {
-    k_.dcidc(p, T, c, dcidc);
+    k_.dcidc(p, T, c, li, dcidc);
 }
 
 
@@ -255,10 +261,11 @@ Foam::scalar Foam::ReversibleReaction
 (
     const scalar p,
     const scalar T,
-    const scalarField& c
+    const scalarField& c,
+    const label li
 ) const
 {
-    return k_.dcidT(p, T, c);
+    return k_.dcidT(p, T, c, li);
 }
 
 

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

@@ -143,7 +143,8 @@ public:
             (
                 const scalar p,
                 const scalar T,
-                const scalarField& c
+                const scalarField& c,
+                const label li
             ) const;
 
             //- Reverse rate constant from the given formard rate constant
@@ -152,7 +153,8 @@ public:
                 const scalar kfwd,
                 const scalar p,
                 const scalar T,
-                const scalarField& c
+                const scalarField& c,
+                const label li
             ) const;
 
             //- Reverse rate constant.
@@ -162,7 +164,8 @@ public:
             (
                 const scalar p,
                 const scalar T,
-                const scalarField& c
+                const scalarField& c,
+                const label li
             ) const;
 
 
@@ -173,7 +176,8 @@ public:
             (
                 const scalar p,
                 const scalar T,
-                const scalarField& c
+                const scalarField& c,
+                const label li
             ) const;
 
             //- Temperature derivative of backward rate
@@ -182,6 +186,7 @@ public:
                 const scalar p,
                 const scalar T,
                 const scalarField& c,
+                const label li,
                 const scalar dkfdT,
                 const scalar kr
             ) const;
@@ -198,6 +203,7 @@ public:
                 const scalar p,
                 const scalar T,
                 const scalarField& c,
+                const label li,
                 scalarField& dcidc
             ) const;
 
@@ -207,7 +213,8 @@ public:
             (
                 const scalar p,
                 const scalar T,
-                const scalarField& c
+                const scalarField& c,
+                const label li
             ) const;
 
 

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

@@ -97,14 +97,16 @@ public:
         (
             const scalar p,
             const scalar T,
-            const scalarField& c
+            const scalarField& c,
+            const label li
         ) const;
 
         inline scalar ddT
         (
             const scalar p,
             const scalar T,
-            const scalarField& c
+            const scalarField& c,
+            const label li
         ) const;
 
         //- Third-body efficiencies (beta = 1-alpha)
@@ -119,6 +121,7 @@ public:
             const scalar p,
             const scalar T,
             const scalarField& c,
+            const label li,
             scalarField& dcidc
         ) const;
 
@@ -128,7 +131,8 @@ public:
         (
             const scalar p,
             const scalar T,
-            const scalarField& c
+            const scalarField& c,
+            const label li
         ) const;
 
         //- Write to stream

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

@@ -56,7 +56,8 @@ inline Foam::scalar Foam::ArrheniusReactionRate::operator()
 (
     const scalar p,
     const scalar T,
-    const scalarField&
+    const scalarField&,
+    const label
 ) const
 {
     scalar ak = A_;
@@ -79,7 +80,8 @@ inline Foam::scalar Foam::ArrheniusReactionRate::ddT
 (
     const scalar p,
     const scalar T,
-    const scalarField&
+    const scalarField&,
+    const label
 ) const
 {
     scalar ak = A_;
@@ -110,6 +112,7 @@ inline void Foam::ArrheniusReactionRate::dcidc
     const scalar p,
     const scalar T,
     const scalarField& c,
+    const label li,
     scalarField& dcidc
 ) const
 {}
@@ -119,7 +122,8 @@ inline Foam::scalar Foam::ArrheniusReactionRate::dcidT
 (
     const scalar p,
     const scalar T,
-    const scalarField& c
+    const scalarField& c,
+    const label li
 ) const
 {
     return 0;

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

@@ -107,14 +107,16 @@ public:
         (
             const scalar p,
             const scalar T,
-            const scalarField& c
+            const scalarField& c,
+            const label li
         ) const;
 
         inline scalar ddT
         (
             const scalar p,
             const scalar T,
-            const scalarField& c
+            const scalarField& c,
+            const label li
         ) const;
 
         //- Third-body efficiencies (beta = 1-alpha)
@@ -132,6 +134,7 @@ public:
             const scalar p,
             const scalar T,
             const scalarField& c,
+            const label li,
             scalarField& dcidc
         ) const;
 
@@ -141,7 +144,8 @@ public:
         (
             const scalar p,
             const scalar T,
-            const scalarField& c
+            const scalarField& c,
+            const label li
         ) const;
 
         //- Write to stream

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2019 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -91,11 +91,12 @@ inline Foam::scalar Foam::ChemicallyActivatedReactionRate
 (
     const scalar p,
     const scalar T,
-    const scalarField& c
+    const scalarField& c,
+    const label li
 ) const
 {
-    const scalar k0 = k0_(p, T, c);
+    const scalar k0 = k0_(p, T, c, li);
-    const scalar kInf = kInf_(p, T, c);
+    const scalar kInf = kInf_(p, T, c, li);
     const scalar Pr = k0*thirdBodyEfficiencies_.M(c)/kInf;
 
     return k0*(1/(1 + Pr))*F_(T, Pr);
@@ -111,14 +112,15 @@ inline Foam::scalar Foam::ChemicallyActivatedReactionRate
 (
     const scalar p,
     const scalar T,
-    const scalarField& c
+    const scalarField& c,
+    const label li
 ) const
 {
-    const scalar k0 = k0_(p, T, c);
+    const scalar k0 = k0_(p, T, c, li);
-    const scalar kInf = kInf_(p, T, c);
+    const scalar kInf = kInf_(p, T, c, li);
     const scalar Pr = k0*thirdBodyEfficiencies_.M(c)/kInf;
 
-    return (1/(1 + Pr))*F_(T, Pr)*k0_.ddT(p, T, c);
+    return (1/(1 + Pr))*F_(T, Pr)*k0_.ddT(p, T, c, li);
 }
 
 
@@ -132,6 +134,7 @@ inline void Foam::ChemicallyActivatedReactionRate
     const scalar p,
     const scalar T,
     const scalarField& c,
+    const label li,
     scalarField& dcidc
 ) const
 {
@@ -139,8 +142,8 @@ inline void Foam::ChemicallyActivatedReactionRate
 
     if (M > small)
     {
-        const scalar k0 = k0_(p, T, c);
+        const scalar k0 = k0_(p, T, c, li);
-        const scalar kInf = kInf_(p, T, c);
+        const scalar kInf = kInf_(p, T, c, li);
 
         const scalar Pr = k0*M/kInf;
         const scalar F = F_(T, Pr);
@@ -171,20 +174,21 @@ inline Foam::scalar Foam::ChemicallyActivatedReactionRate
 (
     const scalar p,
     const scalar T,
-    const scalarField& c
+    const scalarField& c,
+    const label li
 ) const
 {
     const scalar M = thirdBodyEfficiencies_.M(c);
 
     if (M > small)
     {
-        const scalar k0 = k0_(p, T, c);
+        const scalar k0 = k0_(p, T, c, li);
-        const scalar kInf = kInf_(p, T, c);
+        const scalar kInf = kInf_(p, T, c, li);
 
         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);
+            Pr*(k0_.ddT(p, T, c, li)/k0 - kInf_.ddT(p, T, c, li)/kInf - 1/T);
         const scalar dFdT = F_.ddT(Pr, F, dPrdT, T);
 
         return (-dPrdT/(1 + Pr) + dFdT/F);

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

@@ -104,14 +104,16 @@ public:
         (
             const scalar p,
             const scalar T,
-            const scalarField& c
+            const scalarField& c,
+            const label li
         ) const;
 
         inline scalar ddT
         (
             const scalar p,
             const scalar T,
-            const scalarField& c
+            const scalarField& c,
+            const label li
         ) const;
 
         //- Third-body efficiencies (beta = 1-alpha)
@@ -128,6 +130,7 @@ public:
             const scalar p,
             const scalar T,
             const scalarField& c,
+            const label li,
             scalarField& dcidc
         ) const;
 
@@ -136,7 +139,8 @@ public:
         (
             const scalar p,
             const scalar T,
-            const scalarField& c
+            const scalarField& c,
+            const label li
         ) const;
 
         //- Write to stream

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2019 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -82,11 +82,12 @@ Foam::FallOffReactionRate<ReactionRate, FallOffFunction>::operator()
 (
     const scalar p,
     const scalar T,
-    const scalarField& c
+    const scalarField& c,
+    const label li
 ) const
 {
-    const scalar k0 = k0_(p, T, c);
+    const scalar k0 = k0_(p, T, c, li);
-    const scalar kInf = kInf_(p, T, c);
+    const scalar kInf = kInf_(p, T, c, li);
     const scalar Pr = k0*thirdBodyEfficiencies_.M(c)/kInf;
 
     return kInf*(Pr/(1 + Pr))*F_(T, Pr);
@@ -99,14 +100,15 @@ Foam::FallOffReactionRate<ReactionRate, FallOffFunction>::ddT
 (
     const scalar p,
     const scalar T,
-    const scalarField& c
+    const scalarField& c,
+    const label li
 ) const
 {
-    const scalar k0 = k0_(p, T, c);
+    const scalar k0 = k0_(p, T, c, li);
-    const scalar kInf = kInf_(p, T, c);
+    const scalar kInf = kInf_(p, T, c, li);
     const scalar Pr = k0*thirdBodyEfficiencies_.M(c)/kInf;
 
-    return (Pr/(1 + Pr))*F_(T, Pr)*kInf_.ddT(p, T, c);
+    return (Pr/(1 + Pr))*F_(T, Pr)*kInf_.ddT(p, T, c, li);
 }
 
 
@@ -116,6 +118,7 @@ inline void Foam::FallOffReactionRate<ReactionRate, FallOffFunction>::dcidc
     const scalar p,
     const scalar T,
     const scalarField& c,
+    const label li,
     scalarField& dcidc
 ) const
 {
@@ -123,8 +126,8 @@ inline void Foam::FallOffReactionRate<ReactionRate, FallOffFunction>::dcidc
 
     if (M > small)
     {
-        const scalar k0 = k0_(p, T, c);
+        const scalar k0 = k0_(p, T, c, li);
-        const scalar kInf = kInf_(p, T, c);
+        const scalar kInf = kInf_(p, T, c, li);
         const scalar Pr = k0*M/kInf;
         const scalar F = F_(T, Pr);
 
@@ -151,20 +154,21 @@ Foam::FallOffReactionRate<ReactionRate, FallOffFunction>::dcidT
 (
     const scalar p,
     const scalar T,
-    const scalarField& c
+    const scalarField& c,
+    const label li
 ) const
 {
     const scalar M = thirdBodyEfficiencies_.M(c);
 
     if (M > small)
     {
-        const scalar k0 = k0_(p, T, c);
+        const scalar k0 = k0_(p, T, c, li);
-        const scalar kInf = kInf_(p, T, c);
+        const scalar kInf = kInf_(p, T, c, li);
 
         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);
+            Pr*(k0_.ddT(p, T, c, li)/k0 - kInf_.ddT(p, T, c, li)/kInf - 1/T);
         const scalar dFdT = F_.ddT(Pr, F, dPrdT, T);
 
         return (dPrdT/(Pr*(1 + Pr)) + dFdT/F);

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

@@ -100,14 +100,16 @@ public:
         (
             const scalar p,
             const scalar T,
-            const scalarField& c
+            const scalarField& c,
+            const label li
         ) const;
 
         inline scalar ddT
         (
             const scalar p,
             const scalar T,
-            const scalarField& c
+            const scalarField& c,
+            const label li
         ) const;
 
         //- Third-body efficiencies (beta = 1-alpha)
@@ -121,6 +123,7 @@ public:
             const scalar p,
             const scalar T,
             const scalarField& c,
+            const label li,
             scalarField& dcidc
         ) const;
 
@@ -129,7 +132,8 @@ public:
         (
             const scalar p,
             const scalar T,
-            const scalarField& c
+            const scalarField& c,
+            const label li
         ) const;
 
         //- Write to stream

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2019 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -59,7 +59,8 @@ inline Foam::scalar Foam::JanevReactionRate::operator()
 (
     const scalar p,
     const scalar T,
-    const scalarField&
+    const scalarField&,
+    const label
 ) const
 {
     scalar lta = A_;
@@ -93,7 +94,8 @@ inline Foam::scalar Foam::JanevReactionRate::ddT
 (
     const scalar p,
     const scalar T,
-    const scalarField&
+    const scalarField&,
+    const label
 ) const
 {
     scalar lta = A_;
@@ -142,6 +144,7 @@ inline void Foam::JanevReactionRate::dcidc
     const scalar p,
     const scalar T,
     const scalarField& c,
+    const label li,
     scalarField& dcidc
 ) const
 {}
@@ -151,7 +154,8 @@ inline Foam::scalar Foam::JanevReactionRate::dcidT
 (
     const scalar p,
     const scalar T,
-    const scalarField& c
+    const scalarField& c,
+    const label li
 ) const
 {
     return 0;

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

@@ -99,14 +99,16 @@ public:
         (
             const scalar p,
             const scalar T,
-            const scalarField& c
+            const scalarField& c,
+            const label li
         ) const;
 
         inline scalar ddT
         (
             const scalar p,
             const scalar T,
-            const scalarField& c
+            const scalarField& c,
+            const label li
         ) const;
 
         //- Third-body efficiencies (beta = 1-alpha)
@@ -120,6 +122,7 @@ public:
             const scalar p,
             const scalar T,
             const scalarField& c,
+            const label li,
             scalarField& dcidc
         ) const;
 
@@ -128,7 +131,8 @@ public:
         (
             const scalar p,
             const scalar T,
-            const scalarField& c
+            const scalarField& c,
+            const label li
         ) const;
 
         //- Write to stream

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

@@ -62,7 +62,8 @@ inline Foam::scalar Foam::LandauTellerReactionRate::operator()
 (
     const scalar p,
     const scalar T,
-    const scalarField&
+    const scalarField&,
+    const label
 ) const
 {
     scalar lta = A_;
@@ -102,7 +103,8 @@ inline Foam::scalar Foam::LandauTellerReactionRate::ddT
 (
     const scalar p,
     const scalar T,
-    const scalarField&
+    const scalarField&,
+    const label
 ) const
 {
     scalar lta = A_;
@@ -157,6 +159,7 @@ inline void Foam::LandauTellerReactionRate::dcidc
     const scalar p,
     const scalar T,
     const scalarField& c,
+    const label li,
     scalarField& dcidc
 ) const
 {}
@@ -166,7 +169,8 @@ inline Foam::scalar Foam::LandauTellerReactionRate::dcidT
 (
     const scalar p,
     const scalar T,
-    const scalarField& c
+    const scalarField& c,
+    const label li
 ) const
 {
     return 0;

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

@@ -102,14 +102,16 @@ public:
         (
             const scalar p,
             const scalar T,
-            const scalarField& c
+            const scalarField& c,
+            const label li
         ) const;
 
         inline scalar ddT
         (
             const scalar p,
             const scalar T,
-            const scalarField& c
+            const scalarField& c,
+            const label li
         ) const;
 
         //- Third-body efficiencies (beta = 1-alpha)
@@ -123,6 +125,7 @@ public:
             const scalar p,
             const scalar T,
             const scalarField& c,
+            const label li,
             scalarField& dcidc
         ) const;
 
@@ -131,7 +134,8 @@ public:
         (
             const scalar p,
             const scalar T,
-            const scalarField& c
+            const scalarField& c,
+            const label li
         ) const;
 
         //- Write to stream

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

@@ -68,7 +68,8 @@ inline Foam::scalar Foam::LangmuirHinshelwoodReactionRate::operator()
 (
     const scalar p,
     const scalar T,
-    const scalarField& c
+    const scalarField& c,
+    const label li
 ) const
 {
     const scalar c0m = pow(c[r_[0]], m_[0]);
@@ -90,7 +91,8 @@ inline Foam::scalar Foam::LangmuirHinshelwoodReactionRate::ddT
 (
     const scalar p,
     const scalar T,
-    const scalarField& c
+    const scalarField& c,
+    const label li
 ) const
 {
     const scalar c0m = pow(c[r_[0]], m_[0]);
@@ -125,6 +127,7 @@ inline void Foam::LangmuirHinshelwoodReactionRate::dcidc
     const scalar p,
     const scalar T,
     const scalarField& c,
+    const label li,
     scalarField& dcidc
 ) const
 {}
@@ -134,7 +137,8 @@ inline Foam::scalar Foam::LangmuirHinshelwoodReactionRate::dcidT
 (
     const scalar p,
     const scalar T,
-    const scalarField& c
+    const scalarField& c,
+    const label li
 ) const
 {
     return 0;

src/thermophysicalModels/specie/reaction/reactionRate/MichaelisMenten/MichaelisMentenReactionRate.H

@@ -104,14 +104,16 @@ public:
         (
             const scalar p,
             const scalar T,
-            const scalarField& c
+            const scalarField& c,
+            const label li
         ) const;
 
         inline scalar ddT
         (
             const scalar p,
             const scalar T,
-            const scalarField& c
+            const scalarField& c,
+            const label li
         ) const;
 
         //- Third-body efficiencies (beta = 1-alpha)
@@ -125,6 +127,7 @@ public:
             const scalar p,
             const scalar T,
             const scalarField& c,
+            const label li,
             scalarField& dcidc
         ) const;
 
@@ -133,7 +136,8 @@ public:
         (
             const scalar p,
             const scalar T,
-            const scalarField& c
+            const scalarField& c,
+            const label li
         ) const;
 
         //- Write to stream

src/thermophysicalModels/specie/reaction/reactionRate/MichaelisMenten/MichaelisMentenReactionRateI.H

@@ -46,7 +46,8 @@ inline Foam::scalar Foam::MichaelisMentenReactionRate::operator()
 (
     const scalar p,
     const scalar T,
-    const scalarField& c
+    const scalarField& c,
+    const label li
 ) const
 {
     return Vmax_/(Km_ + c[s_]);
@@ -57,7 +58,8 @@ inline Foam::scalar Foam::MichaelisMentenReactionRate::ddT
 (
     const scalar p,
     const scalar T,
-    const scalarField& c
+    const scalarField& c,
+    const label li
 ) const
 {
     return 0;
@@ -76,6 +78,7 @@ inline void Foam::MichaelisMentenReactionRate::dcidc
     const scalar p,
     const scalar T,
     const scalarField& c,
+    const label li,
     scalarField& dcidc
 ) const
 {
@@ -87,7 +90,8 @@ inline Foam::scalar Foam::MichaelisMentenReactionRate::dcidT
 (
     const scalar p,
     const scalar T,
-    const scalarField& c
+    const scalarField& c,
+    const label li
 ) const
 {
     return 0;

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2019 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -84,14 +84,16 @@ public:
         (
             const scalar p,
             const scalar T,
-            const scalarField& c
+            const scalarField& c,
+            const label li
         ) const;
 
         inline scalar ddT
         (
             const scalar p,
             const scalar T,
-            const scalarField& c
+            const scalarField& c,
+            const label li
         ) const;
 
         //- Third-body efficiencies (beta = 1-alpha)
@@ -105,6 +107,7 @@ public:
             const scalar p,
             const scalar T,
             const scalarField& c,
+            const label li,
             scalarField& dcidc
         ) const;
 
@@ -113,7 +116,8 @@ public:
         (
             const scalar p,
             const scalar T,
-            const scalarField& c
+            const scalarField& c,
+            const label li
         ) const;
 
         //- Write to stream

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2019 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -47,7 +47,8 @@ inline Foam::scalar Foam::infiniteReactionRate::operator()
 (
     const scalar p,
     const scalar,
-    const scalarField&
+    const scalarField&,
+    const label
 ) const
 {
     return (1);
@@ -57,7 +58,8 @@ inline Foam::scalar Foam::infiniteReactionRate::ddT
 (
     const scalar p,
     const scalar,
-    const scalarField&
+    const scalarField&,
+    const label
 ) const
 {
     return (0);
@@ -76,6 +78,7 @@ inline void Foam::infiniteReactionRate::dcidc
     const scalar p,
     const scalar T,
     const scalarField& c,
+    const label li,
     scalarField& dcidc
 ) const
 {}
@@ -85,7 +88,8 @@ inline Foam::scalar Foam::infiniteReactionRate::dcidT
 (
     const scalar p,
     const scalar T,
-    const scalarField& c
+    const scalarField& c,
+    const label li
 ) const
 {
     return 0;

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

@@ -100,14 +100,16 @@ public:
         (
             const scalar p,
             const scalar T,
-            const scalarField& c
+            const scalarField& c,
+            const label li
         ) const;
 
         inline scalar ddT
         (
             const scalar p,
             const scalar T,
-            const scalarField& c
+            const scalarField& c,
+            const label li
         ) const;
 
         //- Third-body efficiencies (beta = 1-alpha)
@@ -121,6 +123,7 @@ public:
             const scalar p,
             const scalar T,
             const scalarField& c,
+            const label li,
             scalarField& dcidc
         ) const;
 
@@ -129,7 +132,8 @@ public:
         (
             const scalar p,
             const scalar T,
-            const scalarField& c
+            const scalarField& c,
+            const label li
         ) const;
 
         //- Write to stream

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

@@ -59,7 +59,8 @@ inline Foam::scalar Foam::powerSeriesReactionRate::operator()
 (
     const scalar p,
     const scalar T,
-    const scalarField&
+    const scalarField&,
+    const label
 ) const
 {
     scalar lta = A_;
@@ -86,7 +87,8 @@ inline Foam::scalar Foam::powerSeriesReactionRate::ddT
 (
     const scalar p,
     const scalar T,
-    const scalarField&
+    const scalarField&,
+    const label
 ) const
 {
     scalar lta = A_;
@@ -124,6 +126,7 @@ inline void Foam::powerSeriesReactionRate::dcidc
     const scalar p,
     const scalar T,
     const scalarField& c,
+    const label li,
     scalarField& dcidc
 ) const
 {}
@@ -133,7 +136,8 @@ inline Foam::scalar Foam::powerSeriesReactionRate::dcidT
 (
     const scalar p,
     const scalar T,
-    const scalarField& c
+    const scalarField& c,
+    const label li
 ) const
 {
     return 0;

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

@@ -97,14 +97,16 @@ public:
         (
             const scalar p,
             const scalar T,
-            const scalarField& c
+            const scalarField& c,
+            const label li
         ) const;
 
         inline scalar ddT
         (
             const scalar p,
             const scalar T,
-            const scalarField& c
+            const scalarField& c,
+            const label li
         ) const;
 
         //- Third-body efficiencies (beta = 1-alpha)
@@ -121,6 +123,7 @@ public:
             const scalar p,
             const scalar T,
             const scalarField& c,
+            const label li,
             scalarField& dcidc
         ) const;
 
@@ -129,7 +132,8 @@ public:
         (
             const scalar p,
             const scalar T,
-            const scalarField& c
+            const scalarField& c,
+            const label li
         ) const;
 
         //- Write to stream

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2019 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -76,12 +76,13 @@ inline Foam::scalar Foam::thirdBodyArrheniusReactionRate::operator()
 (
     const scalar p,
     const scalar T,
-    const scalarField& c
+    const scalarField& c,
+    const label li
 ) const
 {
     return
         thirdBodyEfficiencies_.M(c)
-       *ArrheniusReactionRate::operator()(p, T, c);
+       *ArrheniusReactionRate::operator()(p, T, c, li);
 }
 
 
@@ -89,12 +90,13 @@ inline Foam::scalar Foam::thirdBodyArrheniusReactionRate::ddT
 (
     const scalar p,
     const scalar T,
-    const scalarField& c
+    const scalarField& c,
+    const label li
 ) const
 {
     return
         thirdBodyEfficiencies_.M(c)
-       *ArrheniusReactionRate::ddT(p, T, c);
+       *ArrheniusReactionRate::ddT(p, T, c, li);
 }
 
 
@@ -103,6 +105,7 @@ inline void Foam::thirdBodyArrheniusReactionRate::dcidc
     const scalar p,
     const scalar T,
     const scalarField& c,
+    const label li,
     scalarField& dcidc
 ) const
 {
@@ -118,7 +121,8 @@ inline Foam::scalar Foam::thirdBodyArrheniusReactionRate::dcidT
 (
     const scalar p,
     const scalar T,
-    const scalarField& c
+    const scalarField& c,
+    const label li
 ) const
 {
     return -1.0/T;