Merge 5588fea992 into 47517f2ebb

2025-12-08 06:57:46 +00:00 · 2024-11-26 17:13:19 +08:00
parent 47517f2ebb 5588fea992
commit 0b37f85f26
354 changed files with 19277 additions and 69 deletions
--- a/.gitignore
+++ b/.gitignore
@ -84,3 +84,23 @@ doc/Doxygen/DTAGS

 # Ignore the test directory
 /tutorialsTest
+
+#Ignore all tutorial files that are not configs, resources or scripts
+/tutorials
+
+!Allwmake
+!Alltest
+!Allclean
+!Allrun
+!*.yaml
+!*.cti
+!*Properties
+
+!/tutorials/resources/
+
+!/tutorials/**/system/
+!/tutorials/**/0/
+
+!/tutorials/**/constant/*reaction*
+!/tutorials/**/constant/*thermo*
+
--- a/applications/solvers/CTreactingFoam/CTreactingFoam.C
+++ b/applications/solvers/CTreactingFoam/CTreactingFoam.C
@ -0,0 +1,135 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    reactingFoam
+
+Description
+    Solver for combustion with chemical reactions.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "turbulentFluidThermoModel.H"
+#include "psiReactionThermo.H"
+#include "CombustionModel.H"
+#include "multivariateScheme.H"
+#include "pimpleControl.H"
+#include "pressureControl.H"
+#include "fvOptions.H"
+#include "localEulerDdtScheme.H"
+#include "fvcSmooth.H"
+#include "reactingMixture.H"
+#include "thermoPhysicsTypes.H"
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "postProcess.H"
+
+    #include "setRootCaseLists.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "createControl.H"
+    #include "createTimeControls.H"
+    #include "initContinuityErrs.H"
+    #include "createFields.H"
+    #include "createFieldRefs.H"
+
+    turbulence->validate();
+
+    if (!LTS)
+    {
+        #include "compressibleCourantNo.H"
+        #include "setInitialDeltaT.H"
+    }
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.run())
+    {
+        #include "readTimeControls.H"
+
+        if (LTS)
+        {
+            #include "setRDeltaT.H"
+        }
+        else
+        {
+            #include "compressibleCourantNo.H"
+            #include "setDeltaT.H"
+        }
+
+        runTime++;
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        #include "rhoEqn.H"
+
+        while (pimple.loop())
+        {
+            #include "UEqn.H"
+            #include "YEqn.H"
+            #include "EEqn.H"
+
+            // --- Pressure corrector loop
+            while (pimple.correct())
+            {
+                if (pimple.consistent())
+                {
+                    #include "pcEqn.H"
+                }
+                else
+                {
+                    #include "pEqn.H"
+                }
+            }
+
+            if (pimple.turbCorr())
+            {
+                turbulence->correct();
+            }
+        }
+
+        rho = thermo.rho();
+
+        runTime.write();
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+    
+    #include "LewisNumber.H" 
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //
--- a/applications/solvers/CTreactingFoam/EEqn.H
+++ b/applications/solvers/CTreactingFoam/EEqn.H
@ -0,0 +1,47 @@
+{
+    volScalarField tempT = T; 
+    volScalarField& he = thermo.he();
+
+    fvScalarMatrix EEqn
+    (
+        fvm::ddt(rho, he) + mvConvection->fvmDiv(phi, he)
+      + fvc::ddt(rho, K) + fvc::div(phi, K)
+      + (
+            he.name() == "e"
+          ? fvc::div
+            (
+                fvc::absolute(phi/fvc::interpolate(rho), U),
+                p,
+                "div(phiv,p)"
+            )
+          : -dpdt
+        )
+      - fvm::laplacian(turbulence->alphaEff(), he)
+     ==
+        Qdot
+      + fvOptions(rho, he)
+    );
+    
+    forAll(Y, k)
+    {
+        EEqn -= fvc::laplacian(turbulence->alphaEff()*hsi[k], Y[k]);
+        EEqn -= fvc::div(J[k], hsi[k], "div(Ji,hsi)");
+    }
+
+    EEqn.relax();
+
+    fvOptions.constrain(EEqn);
+
+    EEqn.solve();
+
+    fvOptions.correct(he);
+
+    thermo.correct();
+
+    Info << "min/max(T) = "
+         << min(T).value() << ", " << max(T).value() << endl;
+    // update oldT for CoDAC
+    //oldT = tempT;
+}
+
+
--- a/applications/solvers/CTreactingFoam/LewisNumber.H
+++ b/applications/solvers/CTreactingFoam/LewisNumber.H
@ -0,0 +1,26 @@
+/*OFstream LewisNumber
+(
+    "LewisNumber.dat"
+);
+
+LewisNumber
+    << "species" << "      Lewis Number" << endl;
+
+forAll(mesh.cells(), celli)
+{
+    if ( abs( T[celli]-max(T).value() ) < 1.0e-20 )
+    {
+        forAll (Y, k)
+        {
+            scalar Diff =  dynamic_cast<const reactingMixture<gasHThermoPhysics>&>
+                              (composition).D(p, T, Y)[k][celli];
+            if (Diff > 1.0e-20) 
+            {
+                scalar Lewis = thermo.alpha()[celli]/Diff;
+                LewisNumber
+                << k << tab << tab << Lewis << endl;    
+            }
+
+        }
+    }
+}*/
--- a/applications/solvers/CTreactingFoam/Make/files
+++ b/applications/solvers/CTreactingFoam/Make/files
@ -0,0 +1,3 @@
+CTreactingFoam.C
+
+EXE = $(FOAM_APPBIN)/CTreactingFoam
--- a/applications/solvers/CTreactingFoam/Make/options
+++ b/applications/solvers/CTreactingFoam/Make/options
@ -0,0 +1,39 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
+    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
+    -I$(LIB_SRC)/ODE/lnInclude \
+    -I$(LIB_SRC)/combustionModels/lnInclude \
+    -I$(WM_CANTERA_DIR) \
+    -I$(WM_CANTERA_DIR)/cantera \
+    -I$(WM_CANTERA_DIR)/cantera/kinetics \
+    -I$(WM_CANTERA_DIR)/cantera/base \
+    -I$(WM_CANTERA_DIR)/cantera/transport \
+    -I$(WM_CANTERA_DIR)/cantera/thermo \
+    -I$(WM_CANTERA_DIR)/cantera/ext \
+    -I$(WM_CANTERA_DIR)/cantera/ext/fmt
+
+EXE_LIBS = \
+    -L$(FOAM_USER_LIBBIN) \
+    -lfiniteVolume \
+    -lmeshTools \
+    -lsampling \
+    -lturbulenceModels \
+    -lcompressibleTransportModels \
+    -lspecie \
+    -lfluidThermophysicalModels \
+    -lODE \
+    -lchemistryModel \
+    -lreactionThermophysicalModels \
+    -lcompressibleTurbulenceModels \
+    -lfvOptions \
+    -lcombustionModels \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared
--- a/applications/solvers/CTreactingFoam/UEqn.H
+++ b/applications/solvers/CTreactingFoam/UEqn.H
@ -0,0 +1,25 @@
+// Solve the Momentum equation
+
+MRF.correctBoundaryVelocity(U);
+
+tmp<fvVectorMatrix> tUEqn
+(
+    fvm::ddt(rho, U) + fvm::div(phi, U)
+  + MRF.DDt(rho, U)
+  + turbulence->divDevRhoReff(U)
+ ==
+    fvOptions(rho, U)
+);
+fvVectorMatrix& UEqn = tUEqn.ref();
+
+UEqn.relax();
+
+fvOptions.constrain(UEqn);
+
+if (pimple.momentumPredictor())
+{
+    solve(UEqn == -fvc::grad(p));
+
+    fvOptions.correct(U);
+    K = 0.5*magSqr(U);
+}
--- a/applications/solvers/CTreactingFoam/YEqn.H
+++ b/applications/solvers/CTreactingFoam/YEqn.H
@ -0,0 +1,207 @@
+tmp<fv::convectionScheme<scalar>> mvConvection
+(
+    fv::convectionScheme<scalar>::New
+    (
+        mesh,
+        fields,
+        phi,
+        mesh.divScheme("div(phi,Yi_h)")
+    )
+);
+
+
+// Added by Taaresh - Mix Avg Diffusivity (m2/s)
+
+PtrList<volScalarField> D =
+dynamic_cast<const reactingMixture<gasHThermoPhysics>&>
+(composition).D(p, T, Y, rho);
+
+PtrList<volScalarField> hsi(Y.size());
+PtrList<surfaceScalarField> J(Y.size());
+PtrList<volScalarField> tempY(Y.size());
+PtrList<volVectorField> YV(Y.size());
+
+
+forAll(Y, i)
+{
+    hsi.set
+    (
+        i,        
+        new volScalarField
+        (
+            IOobject
+            (
+                "hsi",
+                mesh.time().timeName(),
+                mesh
+            ),
+            mesh,
+            dimensionedScalar
+            (
+                "hsi", 
+                dimEnergy/dimMass, 
+                Zero
+            )
+        )
+    );
+    
+    tempY.set
+    (
+        i,        
+        new volScalarField
+        (
+            IOobject
+            (
+                "tempY",
+                mesh.time().timeName(),
+                mesh
+            ),
+            mesh,
+            scalar(0.0)
+        )
+    );
+
+    // product of mass fraction and diffusion vel of each specie
+
+    YV.set
+    (
+    i,
+    new volVectorField
+    (
+        IOobject
+        (
+        "YV",
+        mesh.time().timeName(),
+        mesh
+        ),
+        mesh,
+        dimensionedVector("YV", dimVelocity, vector(0.0, 0.0, 0.0))
+    )
+    );
+
+
+}
+
+forAll (Y, i)
+{
+    volScalarField& tHsi = hsi[i];    
+    forAll(tHsi, celli)
+    {
+            tHsi[celli] = composition.Hs(i, p[celli], T[celli]);
+    }
+    volScalarField::Boundary& Bf = tHsi.boundaryFieldRef();    
+    forAll(Bf, patchi)
+    {
+        forAll(Bf[patchi], facei)    
+        {
+            Bf[patchi][facei] = 
+            composition.Hs
+            (
+                i,
+                p.boundaryField()[patchi][facei],
+                T.boundaryField()[patchi][facei]
+            );
+        }
+    }
+}
+
+// add turbulence effect on D, assuming constant turbulence Lewis number
+// turbulence thermal diffusivity alphat
+volScalarField alphat = turbulence->alphat();
+//volScalarField nut = turbulence->nut();
+//Info << rho*nut/(dimensionedScalar("SMALL", dimensionSet(1,-1,-1,0,0), Foam::SMALL)+alphat) << endl;
+const scalar Let(thermo.lookupOrDefault<scalar>("Let", 1.0));
+volScalarField Dt = alphat/(rho*Let);
+forAll(Y, i)
+{
+    D[i] = D[i] + Dt;
+}
+
+{
+
+    reaction->correct();
+    Qdot = reaction->Qdot();
+    volScalarField Yt(0.0*Y[0]);
+
+    Vc *= 0.0;
+
+// Compute Correction Velocity
+
+    forAll(Y, i)
+    {
+        Vc += D[i]*fvc::grad(Y[i]) + (D[i]*Y[i]/Wmix)*fvc::grad(Wmix);
+    }
+
+    phiVc = linearInterpolate(rho*Vc) & mesh.Sf();
+
+
+    forAll(Y, i)
+    {
+        if (i != inertIndex && composition.active(i))
+        {
+            volScalarField& Yi = Y[i];
+
+            fvScalarMatrix YiEqn
+            (
+                fvm::ddt(rho, Yi)
+              + mvConvection->fvmDiv(phi, Yi)
+         == fvm::laplacian(rho*D[i], Yi) + fvc::laplacian(rho*D[i]*Yi/Wmix, Wmix)
+        - fvm::div(phiVc, Yi, "div(phiVc,Yi_h)") + 
+                reaction->R(Yi)
+              + fvOptions(rho, Yi)
+            );
+
+            YiEqn.relax();
+
+            fvOptions.constrain(YiEqn);
+
+            YiEqn.solve(mesh.solver("Yi"));
+
+            fvOptions.correct(Yi);
+
+            Yi.max(0.0);
+            Yt += Yi;
+           
+        }
+    }   
+
+    Y[inertIndex] = scalar(1) - Yt;
+    Y[inertIndex].max(0.0);
+
+    volVectorField YV_balance(0.0*YV[0]);
+
+    forAll(YV, i)
+    {
+    // Calculate diff velocity for each species (including inertSpecies) for use in EEqn.H
+
+    YV[i] = (-D[i]*fvc::grad(Y[i])) - (D[i]*Y[i]/Wmix)*fvc::grad(Wmix);
+    
+        if (i != inertIndex && composition.active(i))
+        {
+        YV_balance += YV[i]; 
+        }
+    }
+    
+    YV[inertIndex] = -YV_balance; // since sum(Y[i]*V[i]) over all species = 0 for global mass conservation
+
+    forAll(Y, i)
+    {
+
+    // Species Diffusion Flux term used in EEqn.H
+
+        J.set
+        (
+            i,         
+            linearInterpolate(YV[i]*rho) & mesh.Sf() 
+        );    
+    }
+
+}
+
+
+D_O2 = D[O2_index];
+D_N2 = D[N2_index];
+D_O = D[O_index];
+
+Info << "Max O2 diff = " << max(D[O2_index]) << endl; 
+
--- a/applications/solvers/CTreactingFoam/YEqn_WuSplitting.H
+++ b/applications/solvers/CTreactingFoam/YEqn_WuSplitting.H
@ -0,0 +1,265 @@
+tmp<fv::convectionScheme<scalar>> mvConvection
+(
+    fv::convectionScheme<scalar>::New
+    (
+        mesh,
+        fields,
+        phi,
+        mesh.divScheme("div(phi,Yi_h)")
+    )
+);
+
+
+PtrList<volVectorField> Dm =
+dynamic_cast<const reactingMixture<gasHThermoPhysics>&>
+(composition).Dm(p, T, Y, rho);
+
+PtrList<volScalarField> hsi(Y.size());
+PtrList<surfaceScalarField> J(Y.size());
+PtrList<volScalarField> D(Y.size());
+PtrList<volScalarField> tempY(Y.size()), YSave(Y.size());
+forAll(Y, i)
+{
+    hsi.set
+    (
+        i,        
+        new volScalarField
+        (
+            IOobject
+            (
+                "hsi",
+                mesh.time().timeName(),
+                mesh
+            ),
+            mesh,
+            dimensionedScalar
+            (
+                "hsi", 
+                dimEnergy/dimMass, 
+                Zero
+            )
+        )
+    );
+    
+    tempY.set
+    (
+        i,        
+        new volScalarField
+        (
+            IOobject
+            (
+                "tempY",
+                mesh.time().timeName(),
+                mesh
+            ),
+            mesh,
+            scalar(0.0)
+        )
+    );
+    
+    YSave.set(i, Y[i]);
+    
+    D.set
+    (
+        i, 
+        mag(Dm[i])/(dimensionedScalar("SMALL", dimensionSet(0,-1,0,0,0), Foam::SMALL)
+        +mag(fvc::grad(Y[i]).ref()))        
+    );
+
+    J.set
+    (
+        i,         
+        linearInterpolate(Dm[i]*rho) & mesh.Sf()
+    );    
+}
+
+forAll (Y, i)
+{
+    volScalarField& tHsi = hsi[i];    
+    forAll(tHsi, celli)
+    {
+        tHsi[celli] = composition.Hs(i, p[celli], T[celli]);
+    }
+    volScalarField::Boundary& Bf = tHsi.boundaryFieldRef();    
+    forAll(Bf, patchi)
+    {
+        forAll(Bf[patchi], facei)    
+        {
+            Bf[patchi][facei] = 
+            composition.Hs
+            (
+                i,
+                p.boundaryField()[patchi][facei],
+                T.boundaryField()[patchi][facei]
+            );
+        }
+    }
+}
+
+
+{
+    tempY = Y;
+    volScalarField Yt(0.0*Y[0]);
+    
+    // do simple balanced splitting
+    // reference: Computer Physics Communications 243 (2019) 81–96
+    // evaluate cn term, which will be constant during the splitting integration 
+
+    // -----------------------------step 1, integrate R-cn by dt------------------------------------
+    // cn is included in the chemistry model
+    // but we do not necessarily need cn to be integrated in the chemistry model 
+    // because we can use explicit for Transport (Cn) and then S, it will be no splitting error   
+    scalar dtSave = runTime.deltaT().value();
+    reaction->correct();
+    Qdot = reaction->Qdot();
+     
+    // --------Substep 1a, integrate the transport term explicitly by dt--------
+    forAll(Y, i)
+    {
+        if (i != inertIndex && composition.active(i))
+        {
+            volScalarField& Yi = Y[i];
+
+            fvScalarMatrix YiEqn
+            (
+                fvm::ddt(rho, Yi)
+              + mvConvection->fvcDiv(phi, Yi)
+              - fvc::laplacian(rho*D[i], Yi)
+              == fvOptions(rho, Yi)
+            );
+
+            fvOptions.constrain(YiEqn);
+
+            YiEqn.solve(mesh.solver("Yi"));
+
+            fvOptions.correct(Yi);
+
+            Yi.max(0.0);
+            Yt += Yi;
+
+        }
+    }
+
+    Y[inertIndex] = scalar(1) - Yt;
+    Y[inertIndex].max(0.0);
+    //YSave[inertIndex] = Y[inertIndex];
+
+
+    // update Rho oldTime field, this could be done with a smarter way: write 
+    //    fvm::ddt(rho, Yi)-fvc::ddt(rho, Yi) == reaction->R(Yi)
+    // for the transport Equation, here we still use the "update oldTime field" way
+    volScalarField& tRho = rho.oldTime(); 
+    volScalarField::Boundary& Bf = tRho.boundaryFieldRef();    
+    forAll (tRho, celli)
+    {
+        tRho[celli] = rho[celli];
+    }
+ 
+    forAll(Bf, patchi)
+    {
+        forAll(Bf[patchi], facei)   
+        {
+            Bf[patchi][facei] = rho.boundaryField()[patchi][facei];
+        }
+    } 
+    
+    //update Y.oldTime(), which will be used in fvm::ddt
+    forAll (Y, i)
+    {
+        volScalarField& tYi = Y[i].oldTime();
+        
+        forAll(tYi, celli)
+        {
+            tYi[celli] = Y[i][celli];
+        }
+        volScalarField::Boundary& Bf = tYi.boundaryFieldRef(); 
+        forAll(Bf, patchi)
+        {
+            forAll(Bf[patchi], facei)   
+            {
+                Bf[patchi][facei] = Y[i].boundaryField()[patchi][facei];
+            }
+        }
+    } 
+    // --------Substep 1b, integrate the chemical source term by dt--------    
+    forAll(Y, i)
+    {
+        if (i != inertIndex && composition.active(i))
+        {
+            volScalarField& Yi = Y[i];
+            fvScalarMatrix YiEqn
+            (
+               fvm::ddt(rho, Yi)
+             ==
+               reaction->R(Yi)
+            );
+            YiEqn.relax();
+
+            YiEqn.solve(mesh.solver("Yi"));
+
+            Yi.max(0.0);
+            Yt += Yi;           
+        }
+    }
+    Y[inertIndex] = scalar(1) - Yt;
+    Y[inertIndex].max(0.0);     
+    
+    // ----------------------- step 2, integrate T+cn by dt/2 ----------------------------
+    //update Y.oldTime(), which will be used in fvm::ddt
+    forAll (Y, i)
+    {
+        volScalarField& tYi = Y[i].oldTime();
+        
+        forAll(tYi, celli)
+        {
+            tYi[celli] = Y[i][celli];
+        }
+        volScalarField::Boundary& Bf = tYi.boundaryFieldRef(); 
+        forAll(Bf, patchi)
+        {
+            forAll(Bf[patchi], facei)   
+            {
+                Bf[patchi][facei] = Y[i].boundaryField()[patchi][facei];
+            }
+        }
+    }     
+    Yt = 0.0;
+    runTime.setDeltaT(dtSave/2.0);
+    forAll(Y, i)
+    {
+        if (i != inertIndex && composition.active(i))
+        {
+            volScalarField& Yi = Y[i];
+
+            fvScalarMatrix YiEqn
+            (
+                fvm::ddt(rho, Yi)
+                + mvConvection->fvmDiv(phi, Yi)
+                - fvm::laplacian(rho*D[i], Yi)
+             ==
+                mvConvection->fvcDiv(phi, YSave[i])
+                - fvc::laplacian(rho*D[i], YSave[i])
+              + fvOptions(rho, Yi)
+            );
+
+            YiEqn.relax();
+
+            fvOptions.constrain(YiEqn);
+
+            YiEqn.solve(mesh.solver("Yi"));
+
+            fvOptions.correct(Yi);
+
+            Yi.max(0.0);
+            Yt += Yi;
+        }
+    }
+    Y[inertIndex] = scalar(1) - Yt;
+    Y[inertIndex].max(0.0);     
+
+    runTime.setDeltaT(dtSave);   
+    // update oldY here for CoDAC
+    // oldY = tempY;
+}
+
+
--- a/applications/solvers/CTreactingFoam/createFieldRefs.H
+++ b/applications/solvers/CTreactingFoam/createFieldRefs.H
@ -0,0 +1,3 @@
+const volScalarField& psi = thermo.psi();
+const volScalarField& T = thermo.T();
+const label inertIndex(composition.species()[inertSpecie]);
--- a/applications/solvers/CTreactingFoam/createFields.H
+++ b/applications/solvers/CTreactingFoam/createFields.H
@ -0,0 +1,286 @@
+#include "createRDeltaT.H"
+
+Info<< "Reading thermophysical properties\n" << endl;
+autoPtr<psiReactionThermo> pThermo(psiReactionThermo::New(mesh));
+
+psiReactionThermo& thermo = pThermo();
+thermo.validate(args.executable(), "h", "e");
+
+basicSpecieMixture& composition = thermo.composition();
+PtrList<volScalarField>& Y = composition.Y();
+
+// initilize omega, diff, omegaT and diffT for use in CEMA
+//PtrList<volScalarField>& omega = composition.omega();
+//PtrList<volScalarField>& diff = composition.diff();
+//volScalarField& omegaT = composition.omegaT();
+//volScalarField& diffT = composition.diffT();
+
+// initilize oldY for CoDAC, it will be large numbers to aviod similarity initially
+//PtrList<volScalarField>& oldY = composition.oldY(); 
+//forAll(Y, i)
+//{
+//    oldY[i] = 1000.0;
+//}
+const word inertSpecie(thermo.lookup("inertSpecie"));
+if (!composition.species().found(inertSpecie))
+{
+    FatalIOErrorIn(args.executable().c_str(), thermo)
+        << "Inert specie " << inertSpecie << " not found in available species "
+        << composition.species() << exit(FatalIOError);
+}
+
+volScalarField rho
+(
+    IOobject
+    (
+        "rho",
+        runTime.timeName(),
+        mesh
+    ),
+    thermo.rho()
+);
+
+Info<< "Reading field U\n" << endl;
+volVectorField U
+(
+    IOobject
+    (
+        "U",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);
+
+
+volVectorField Vc
+(
+    IOobject
+    (
+        "Vc",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh,
+    dimensionedVector("Vc", dimVelocity, vector(0,0,0))
+);
+
+volScalarField& p = thermo.p();
+
+#include "compressibleCreatePhi.H"
+
+pressureControl pressureControl(p, rho, pimple.dict(), false);
+
+mesh.setFluxRequired(p.name());
+
+Info << "Creating turbulence model.\n" << nl;
+autoPtr<compressible::turbulenceModel> turbulence
+(
+    compressible::turbulenceModel::New
+    (
+        rho,
+        U,
+        phi,
+        thermo
+    )
+);
+
+Info<< "Creating reaction model\n" << endl;
+autoPtr<CombustionModel<psiReactionThermo>> reaction
+(
+    CombustionModel<psiReactionThermo>::New(thermo, turbulence())
+);
+
+Info<< "Creating field dpdt\n" << endl;
+volScalarField dpdt
+(
+    IOobject
+    (
+        "dpdt",
+        runTime.timeName(),
+        mesh
+    ),
+    mesh,
+    dimensionedScalar("dpdt", p.dimensions()/dimTime, 0)
+);
+
+Info<< "Creating field kinetic energy K\n" << endl;
+volScalarField K("K", 0.5*magSqr(U));
+
+
+multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
+
+forAll(Y, i)
+{
+    fields.add(Y[i]);
+}
+fields.add(thermo.he());
+
+volScalarField Qdot
+(
+    IOobject
+    (
+        "Qdot",
+        runTime.timeName(),
+        mesh,
+        IOobject::READ_IF_PRESENT,
+        IOobject::AUTO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("Qdot", dimEnergy/dimVolume/dimTime, 0.0)
+);
+
+// add mixture fraction field
+volScalarField Zmix 
+(
+    IOobject
+    (
+        "Zmix",
+        runTime.timeName(),
+        mesh,
+        IOobject::READ_IF_PRESENT,
+        IOobject::AUTO_WRITE
+    ),
+    mesh,
+    scalar(0.0)
+);
+
+// add progress variable field
+volScalarField Prog 
+(
+    IOobject
+    (
+        "Prog",
+        runTime.timeName(),
+        mesh,
+        IOobject::READ_IF_PRESENT,
+        IOobject::AUTO_WRITE
+    ),
+    mesh,
+    scalar(0.0)
+);
+#include "createMRF.H"
+#include "createFvOptions.H"
+
+label O_index, O2_index, N2_index;
+
+forAll(Y,i)
+{
+    if (Y[i].member() == "O")
+    {
+        O_index = i;
+        Info << "O_index = " << O_index << endl;
+        Info << Y[O_index].member() << endl;
+    }
+    else if (Y[i].member() == "O2")
+    {
+        O2_index = i;
+        Info << "O2_index = " << O2_index << endl;
+        Info << Y[O2_index].member() << endl;
+    }
+    else if (Y[i].member() == "N2")
+    {
+        N2_index = i;
+        Info << "N2 index = " << N2_index <<  endl;
+    }
+
+}
+
+volScalarField D_O2
+{
+    IOobject
+    (
+        "D_O2",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("D_O2", dimArea/dimTime, 0.0)
+};
+
+volScalarField D_N2
+{
+    IOobject
+    (
+        "D_N2",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("D_N2", dimArea/dimTime, 0.0)
+};
+
+
+volScalarField D_O
+{
+    IOobject
+    (
+        "D_O",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("D_O", dimArea/dimTime, 0.0)
+};
+
+
+
+volScalarField gradY_O2
+{
+    IOobject
+    (
+        "gradY_O2",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("gradY_O2", dimless/dimLength, 0.0)
+};
+
+OFstream maxDataFile
+(
+        "max_data_myrf.dat"
+);
+
+volScalarField Wmix
+{
+    IOobject
+    (
+        "Wmix",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    thermo.W()
+};
+
+
+surfaceScalarField phiVc
+(
+    IOobject
+    (
+        "phiVc",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    linearInterpolate(rho*Vc) & mesh.Sf()
+);
+
+
+
+
--- a/applications/solvers/CTreactingFoam/pEqn.H
+++ b/applications/solvers/CTreactingFoam/pEqn.H
@ -0,0 +1,104 @@
+rho = thermo.rho();
+
+volScalarField rAU(1.0/UEqn.A());
+surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
+volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
+
+if (pimple.nCorrPISO() <= 1)
+{
+    tUEqn.clear();
+}
+
+if (pimple.transonic())
+{
+    surfaceScalarField phid
+    (
+        "phid",
+        fvc::interpolate(psi)
+       *(
+            fvc::flux(HbyA)
+          + MRF.zeroFilter
+            (
+                rhorAUf*fvc::ddtCorr(rho, U, phi)/fvc::interpolate(rho)
+            )
+        )
+    );
+
+    MRF.makeRelative(fvc::interpolate(psi), phid);
+
+    while (pimple.correctNonOrthogonal())
+    {
+        fvScalarMatrix pEqn
+        (
+            fvm::ddt(psi, p)
+          + fvm::div(phid, p)
+          - fvm::laplacian(rhorAUf, p)
+         ==
+            fvOptions(psi, p, rho.name())
+        );
+
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
+
+        if (pimple.finalNonOrthogonalIter())
+        {
+            phi == pEqn.flux();
+        }
+    }
+}
+else
+{
+    surfaceScalarField phiHbyA
+    (
+        "phiHbyA",
+        (
+            fvc::flux(rho*HbyA)
+          + MRF.zeroFilter(rhorAUf*fvc::ddtCorr(rho, U, phi))
+        )
+    );
+
+    MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
+
+    // Update the pressure BCs to ensure flux consistency
+    constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF);
+
+    while (pimple.correctNonOrthogonal())
+    {
+        fvScalarMatrix pEqn
+        (
+            fvm::ddt(psi, p)
+          + fvc::div(phiHbyA)
+          - fvm::laplacian(rhorAUf, p)
+         ==
+            fvOptions(psi, p, rho.name())
+        );
+
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
+
+        if (pimple.finalNonOrthogonalIter())
+        {
+            phi = phiHbyA + pEqn.flux();
+        }
+    }
+}
+
+#include "rhoEqn.H"
+#include "compressibleContinuityErrs.H"
+
+// Explicitly relax pressure for momentum corrector
+p.relax();
+
+U = HbyA - rAU*fvc::grad(p);
+U.correctBoundaryConditions();
+fvOptions.correct(U);
+K = 0.5*magSqr(U);
+
+if (pressureControl.limit(p))
+{
+    p.correctBoundaryConditions();
+    rho = thermo.rho();
+}
+
+if (thermo.dpdt())
+{
+    dpdt = fvc::ddt(p);
+}
--- a/applications/solvers/CTreactingFoam/pcEqn.H
+++ b/applications/solvers/CTreactingFoam/pcEqn.H
@ -0,0 +1,121 @@
+rho = thermo.rho();
+
+volScalarField rAU(1.0/UEqn.A());
+volScalarField rAtU(1.0/(1.0/rAU - UEqn.H1()));
+volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
+
+if (pimple.nCorrPISO() <= 1)
+{
+    tUEqn.clear();
+}
+
+if (pimple.transonic())
+{
+    surfaceScalarField phid
+    (
+        "phid",
+        fvc::interpolate(psi)
+       *(
+            fvc::flux(HbyA)
+          + MRF.zeroFilter
+            (
+                fvc::interpolate(rho*rAU)*fvc::ddtCorr(rho, U, phi)
+               /fvc::interpolate(rho)
+            )
+        )
+    );
+
+    MRF.makeRelative(fvc::interpolate(psi), phid);
+
+    surfaceScalarField phic
+    (
+        "phic",
+        fvc::interpolate(rho*(rAtU - rAU))*fvc::snGrad(p)*mesh.magSf()
+    );
+
+    HbyA -= (rAU - rAtU)*fvc::grad(p);
+
+    volScalarField rhorAtU("rhorAtU", rho*rAtU);
+
+    while (pimple.correctNonOrthogonal())
+    {
+        fvScalarMatrix pEqn
+        (
+            fvm::ddt(psi, p)
+          + fvm::div(phid, p)
+          + fvc::div(phic)
+          - fvm::laplacian(rhorAtU, p)
+         ==
+            fvOptions(psi, p, rho.name())
+        );
+
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
+
+        if (pimple.finalNonOrthogonalIter())
+        {
+            phi == phic + pEqn.flux();
+        }
+    }
+}
+else
+{
+    surfaceScalarField phiHbyA
+    (
+        "phiHbyA",
+        (
+            fvc::flux(rho*HbyA)
+          + MRF.zeroFilter(fvc::interpolate(rho*rAU)*fvc::ddtCorr(rho, U, phi))
+        )
+    );
+
+    MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
+
+    phiHbyA += fvc::interpolate(rho*(rAtU - rAU))*fvc::snGrad(p)*mesh.magSf();
+    HbyA -= (rAU - rAtU)*fvc::grad(p);
+
+    volScalarField rhorAtU("rhorAtU", rho*rAtU);
+
+    // Update the pressure BCs to ensure flux consistency
+    constrainPressure(p, rho, U, phiHbyA, rhorAtU, MRF);
+
+    while (pimple.correctNonOrthogonal())
+    {
+        fvScalarMatrix pEqn
+        (
+            fvm::ddt(psi, p)
+          + fvc::div(phiHbyA)
+          - fvm::laplacian(rhorAtU, p)
+         ==
+            fvOptions(psi, p, rho.name())
+        );
+
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
+
+        if (pimple.finalNonOrthogonalIter())
+        {
+            phi = phiHbyA + pEqn.flux();
+        }
+    }
+}
+
+#include "rhoEqn.H"
+#include "compressibleContinuityErrs.H"
+
+// Explicitly relax pressure for momentum corrector
+p.relax();
+
+U = HbyA - rAtU*fvc::grad(p);
+U.correctBoundaryConditions();
+fvOptions.correct(U);
+K = 0.5*magSqr(U);
+
+if (pressureControl.limit(p))
+{
+    p.correctBoundaryConditions();
+    rho = thermo.rho();
+}
+
+if (thermo.dpdt())
+{
+    dpdt = fvc::ddt(p);
+}
--- a/applications/solvers/CTreactingFoam/setRDeltaT.H
+++ b/applications/solvers/CTreactingFoam/setRDeltaT.H
@ -0,0 +1,186 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2013-2018 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+{
+    volScalarField& rDeltaT = trDeltaT.ref();
+
+    const dictionary& pimpleDict = pimple.dict();
+
+    // Maximum flow Courant number
+    scalar maxCo(readScalar(pimpleDict.lookup("maxCo")));
+
+    // Maximum time scale
+    scalar maxDeltaT(pimpleDict.lookupOrDefault<scalar>("maxDeltaT", great));
+
+    // Smoothing parameter (0-1) when smoothing iterations > 0
+    scalar rDeltaTSmoothingCoeff
+    (
+        pimpleDict.lookupOrDefault<scalar>("rDeltaTSmoothingCoeff", 0.1)
+    );
+
+    // Damping coefficient (1-0)
+    scalar rDeltaTDampingCoeff
+    (
+        pimpleDict.lookupOrDefault<scalar>("rDeltaTDampingCoeff", 1.0)
+    );
+
+    // Maximum change in cell temperature per iteration
+    // (relative to previous value)
+    scalar alphaTemp(pimpleDict.lookupOrDefault("alphaTemp", 0.05));
+
+    // Maximum change in cell concentration per iteration
+    // (relative to reference value)
+    scalar alphaY(pimpleDict.lookupOrDefault("alphaY", 1.0));
+
+    Info<< "Time scales min/max:" << endl;
+
+    // Cache old reciprocal time scale field
+    volScalarField rDeltaT0("rDeltaT0", rDeltaT);
+
+    // Flow time scale
+    {
+        rDeltaT.ref() =
+        (
+            fvc::surfaceSum(mag(phi))()()
+           /((2*maxCo)*mesh.V()*rho())
+        );
+
+        // Limit the largest time scale
+        rDeltaT.max(1/maxDeltaT);
+
+        Info<< "    Flow        = "
+            << 1/gMax(rDeltaT.primitiveField()) << ", "
+            << 1/gMin(rDeltaT.primitiveField()) << endl;
+    }
+
+    // Heat release rate time scale
+    if (alphaTemp < 1)
+    {
+        volScalarField::Internal rDeltaTT
+        (
+            mag(Qdot)/(alphaTemp*rho*thermo.Cp()*T)
+        );
+
+        Info<< "    Temperature = "
+            << 1/(gMax(rDeltaTT.field()) + vSmall) << ", "
+            << 1/(gMin(rDeltaTT.field()) + vSmall) << endl;
+
+        rDeltaT.ref() = max(rDeltaT(), rDeltaTT);
+    }
+
+    // Reaction rate time scale
+    if (alphaY < 1)
+    {
+        dictionary Yref(pimpleDict.subDict("Yref"));
+
+        volScalarField::Internal rDeltaTY
+        (
+            IOobject
+            (
+                "rDeltaTY",
+                runTime.timeName(),
+                mesh
+            ),
+            mesh,
+            dimensionedScalar("rDeltaTY", rDeltaT.dimensions(), 0)
+        );
+
+        bool foundY = false;
+
+        forAll(Y, i)
+        {
+            if (i != inertIndex && composition.active(i))
+            {
+                volScalarField& Yi = Y[i];
+
+                if (Yref.found(Yi.name()))
+                {
+                    foundY = true;
+                    scalar Yrefi = readScalar(Yref.lookup(Yi.name()));
+
+                    rDeltaTY.field() = max
+                    (
+                        mag
+                        (
+                            reaction->R(Yi)().source()
+                           /((Yrefi*alphaY)*(rho*mesh.V()))
+                        ),
+                        rDeltaTY
+                    );
+                }
+            }
+        }
+
+        if (foundY)
+        {
+            Info<< "    Composition = "
+                << 1/(gMax(rDeltaTY.field()) + vSmall) << ", "
+                << 1/(gMin(rDeltaTY.field()) + vSmall) << endl;
+
+            rDeltaT.ref() = max(rDeltaT(), rDeltaTY);
+        }
+        else
+        {
+            IOWarningIn(args.executable().c_str(), Yref)
+                << "Cannot find any active species in Yref " << Yref
+                << endl;
+        }
+    }
+
+    // Update tho boundary values of the reciprocal time-step
+    rDeltaT.correctBoundaryConditions();
+
+    // Spatially smooth the time scale field
+    if (rDeltaTSmoothingCoeff < 1)
+    {
+        fvc::smooth(rDeltaT, rDeltaTSmoothingCoeff);
+    }
+
+    // Limit rate of change of time scale
+    // - reduce as much as required
+    // - only increase at a fraction of old time scale
+    if
+    (
+        rDeltaTDampingCoeff < 1
+     && runTime.timeIndex() > runTime.startTimeIndex() + 1
+    )
+    {
+        rDeltaT = max
+        (
+            rDeltaT,
+            (scalar(1) - rDeltaTDampingCoeff)*rDeltaT0
+        );
+    }
+
+    // Update tho boundary values of the reciprocal time-step
+    rDeltaT.correctBoundaryConditions();
+
+    Info<< "    Overall     = "
+        << 1/gMax(rDeltaT.primitiveField())
+        << ", " << 1/gMin(rDeltaT.primitiveField()) << endl;
+}
+
+
+// ************************************************************************* //
--- a/applications/solvers/DNS/dnsFoam/Make/options
+++ b/applications/solvers/DNS/dnsFoam/Make/options
@ -4,6 +4,8 @@ EXE_INC = \
    -I$(LIB_SRC)/meshTools/lnInclude

 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -lrandomProcesses \
    -lfiniteVolume \
    -lmeshTools
--- a/applications/solvers/basic/laplacianFoam/Make/options
+++ b/applications/solvers/basic/laplacianFoam/Make/options
@ -5,4 +5,6 @@ EXE_INC = \
 EXE_LIBS = \
    -lfiniteVolume \
    -lfvOptions \
-    -lmeshTools
+    -lmeshTools \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared
--- a/applications/solvers/basic/potentialFoam/Make/options
+++ b/applications/solvers/basic/potentialFoam/Make/options
@ -4,6 +4,8 @@ EXE_INC = \
    -I$(LIB_SRC)/sampling/lnInclude

 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -lfiniteVolume \
    -lmeshTools \
    -lsampling
--- a/applications/solvers/basic/scalarTransportFoam/Make/options
+++ b/applications/solvers/basic/scalarTransportFoam/Make/options
@ -4,6 +4,8 @@ EXE_INC = \
    -I$(LIB_SRC)/sampling/lnInclude

 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -lfiniteVolume \
    -lfvOptions \
    -lmeshTools \
--- a/applications/solvers/combustion/PDRFoam/Make/options
+++ b/applications/solvers/combustion/PDRFoam/Make/options
@ -34,6 +34,8 @@ EXE_INC = \
    -I$(LIB_SRC)/triSurface/lnInclude \

 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -lengine \
    -lmeshTools \
    -lturbulenceModels \
--- a/applications/solvers/combustion/XiFoam/Make/options
+++ b/applications/solvers/combustion/XiFoam/Make/options
@ -12,6 +12,8 @@ EXE_INC = \
    -I$(LIB_SRC)/meshTools/lnInclude

 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -lengine \
    -lturbulenceModels \
    -lcompressibleTurbulenceModels \
--- a/applications/solvers/combustion/XiFoam/XiEngineFoam/Make/options
+++ b/applications/solvers/combustion/XiFoam/XiEngineFoam/Make/options
@ -13,6 +13,8 @@ EXE_INC = \
    -I$(LIB_SRC)/thermophysicalModels/laminarFlameSpeed/lnInclude

 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -lfiniteVolume \
    -lfvOptions \
    -lsampling \
--- a/applications/solvers/combustion/chemFoam/Make/options
+++ b/applications/solvers/combustion/chemFoam/Make/options
@ -6,9 +6,19 @@ EXE_INC = \
    -I$(LIB_SRC)/ODE/lnInclude\
    -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(WM_CANTERA_DIR) \
+    -I$(WM_CANTERA_DIR)/cantera \
+    -I$(WM_CANTERA_DIR)/cantera/kinetics \
+    -I$(WM_CANTERA_DIR)/cantera/base \
+    -I$(WM_CANTERA_DIR)/cantera/transport \
+    -I$(WM_CANTERA_DIR)/cantera/thermo \
+    -I$(WM_CANTERA_DIR)/cantera/ext \
+    -I$(WM_CANTERA_DIR)/cantera/ext/fmt

 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -lchemistryModel \
    -lfiniteVolume \
-    -lmeshTools
+    -lmeshTools 
--- a/applications/solvers/combustion/coldEngineFoam/Make/options
+++ b/applications/solvers/combustion/coldEngineFoam/Make/options
@ -15,6 +15,8 @@ EXE_INC = \
    -I$(LIB_SRC)/meshTools/lnInclude \

 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -lengine \
    -lturbulenceModels \
    -lcompressibleTurbulenceModels \
--- a/applications/solvers/combustion/fireFoam/Make/options
+++ b/applications/solvers/combustion/fireFoam/Make/options
@ -22,9 +22,19 @@ EXE_INC = \
    -I$(LIB_SRC)/regionModels/pyrolysisModels/lnInclude \
    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
    -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
-    -I$(LIB_SRC)/ODE/lnInclude
+    -I$(LIB_SRC)/ODE/lnInclude \
+    -I$(WM_CANTERA_DIR) \
+    -I$(WM_CANTERA_DIR)/cantera \
+    -I$(WM_CANTERA_DIR)/cantera/kinetics \
+    -I$(WM_CANTERA_DIR)/cantera/base \
+    -I$(WM_CANTERA_DIR)/cantera/transport \
+    -I$(WM_CANTERA_DIR)/cantera/thermo \
+    -I$(WM_CANTERA_DIR)/cantera/ext \
+    -I$(WM_CANTERA_DIR)/cantera/ext/fmt

 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -lfiniteVolume \
    -lfvOptions \
    -lmeshTools \
@ -47,4 +57,4 @@ EXE_LIBS = \
    -lregionCoupling \
    -llagrangianIntermediate \
    -llagrangianTurbulence \
-    -lODE
+    -lODE 
--- a/applications/solvers/combustion/reactingFoam/Make/options
+++ b/applications/solvers/combustion/reactingFoam/Make/options
@ -10,9 +10,19 @@ EXE_INC = \
    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
    -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
    -I$(LIB_SRC)/ODE/lnInclude \
-    -I$(LIB_SRC)/combustionModels/lnInclude
+    -I$(LIB_SRC)/combustionModels/lnInclude \
+    -I$(WM_CANTERA_DIR) \
+    -I$(WM_CANTERA_DIR)/cantera \
+    -I$(WM_CANTERA_DIR)/cantera/kinetics \
+    -I$(WM_CANTERA_DIR)/cantera/base \
+    -I$(WM_CANTERA_DIR)/cantera/transport \
+    -I$(WM_CANTERA_DIR)/cantera/thermo \
+    -I$(WM_CANTERA_DIR)/cantera/ext \
+    -I$(WM_CANTERA_DIR)/cantera/ext/fmt

 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -lfiniteVolume \
    -lfvOptions \
    -lmeshTools \
--- a/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/Make/options
+++ b/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/Make/options
@ -16,6 +16,8 @@ EXE_INC = \


 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -lfiniteVolume \
    -lfvOptions \
    -lmeshTools \
--- a/applications/solvers/combustion/reactingFoam/rhoReactingFoam/Make/options
+++ b/applications/solvers/combustion/reactingFoam/rhoReactingFoam/Make/options
@ -16,6 +16,8 @@ EXE_INC = \


 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -lfiniteVolume \
    -lfvOptions \
    -lmeshTools \
--- a/applications/solvers/compressible/rhoCentralFoam/Make/options
+++ b/applications/solvers/compressible/rhoCentralFoam/Make/options
@ -10,6 +10,8 @@ EXE_INC = \
    -I$(LIB_SRC)/meshTools/lnInclude

 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -lfiniteVolume \
    -lfvOptions \
    -lcompressibleTransportModels \
--- a/applications/solvers/compressible/rhoCentralFoam/rhoCentralDyMFoam/Make/options
+++ b/applications/solvers/compressible/rhoCentralFoam/rhoCentralDyMFoam/Make/options
@ -12,6 +12,8 @@ EXE_INC = \
    -I$(LIB_SRC)/meshTools/lnInclude

 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -lfiniteVolume \
    -lfvOptions \
    -lcompressibleTransportModels \
--- a/applications/solvers/compressible/rhoPimpleFoam/Make/options
+++ b/applications/solvers/compressible/rhoPimpleFoam/Make/options
@ -11,6 +11,8 @@ EXE_INC = \
    -I$(LIB_SRC)/dynamicMesh/lnInclude

 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -lcompressibleTransportModels \
    -lfluidThermophysicalModels \
    -lspecie \
--- a/applications/solvers/compressible/rhoSimpleFoam/Make/options
+++ b/applications/solvers/compressible/rhoSimpleFoam/Make/options
@ -9,6 +9,8 @@ EXE_INC = \
    -I$(LIB_SRC)/meshTools/lnInclude \

 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -lcompressibleTransportModels \
    -lfluidThermophysicalModels \
    -lspecie \
--- a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousSimpleFoam/Make/options
+++ b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousSimpleFoam/Make/options
@ -11,6 +11,8 @@ EXE_INC = \
    -I$(LIB_SRC)/meshTools/lnInclude \

 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -lcompressibleTransportModels \
    -lfluidThermophysicalModels \
    -lspecie \
--- a/applications/solvers/compressible/sonicFoam/Make/options
+++ b/applications/solvers/compressible/sonicFoam/Make/options
@ -8,6 +8,8 @@ EXE_INC = \
    -I$(LIB_SRC)/finiteVolume/lnInclude

 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -lcompressibleTransportModels \
    -lfluidThermophysicalModels \
    -lspecie \
--- a/applications/solvers/compressible/sonicFoam/sonicDyMFoam/Make/options
+++ b/applications/solvers/compressible/sonicFoam/sonicDyMFoam/Make/options
@ -13,6 +13,8 @@ EXE_INC = \
    -I$(LIB_SRC)/meshTools/lnInclude \

 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -lcompressibleTransportModels \
    -lfluidThermophysicalModels \
    -lspecie \
--- a/applications/solvers/compressible/sonicFoam/sonicLiquidFoam/Make/options
+++ b/applications/solvers/compressible/sonicFoam/sonicLiquidFoam/Make/options
@ -3,5 +3,7 @@ EXE_INC = \
    -I$(LIB_SRC)/meshTools/lnInclude

 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -lfiniteVolume \
    -lmeshTools
--- a/applications/solvers/discreteMethods/dsmc/dsmcFoam/Make/options
+++ b/applications/solvers/discreteMethods/dsmc/dsmcFoam/Make/options
@ -5,6 +5,8 @@ EXE_INC = \
    -I$(LIB_SRC)/meshTools/lnInclude

 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -lmeshTools \
    -lfiniteVolume \
    -llagrangian \
--- a/applications/solvers/discreteMethods/molecularDynamics/mdEquilibrationFoam/Make/options
+++ b/applications/solvers/discreteMethods/molecularDynamics/mdEquilibrationFoam/Make/options
@ -7,6 +7,8 @@ EXE_INC = \
    -I$(LIB_SRC)/meshTools/lnInclude

 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -lmeshTools \
    -lfiniteVolume \
    -llagrangian \
--- a/applications/solvers/discreteMethods/molecularDynamics/mdFoam/Make/options
+++ b/applications/solvers/discreteMethods/molecularDynamics/mdFoam/Make/options
@ -7,6 +7,8 @@ EXE_INC = \
    -I$(LIB_SRC)/meshTools/lnInclude

 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -lmeshTools \
    -lfiniteVolume \
    -llagrangian \
--- a/applications/solvers/electromagnetics/electrostaticFoam/Make/options
+++ b/applications/solvers/electromagnetics/electrostaticFoam/Make/options
@ -3,5 +3,7 @@ EXE_INC = \
    -I$(LIB_SRC)/meshTools/lnInclude

 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -lfiniteVolume \
    -lmeshTools
--- a/applications/solvers/electromagnetics/magneticFoam/Make/options
+++ b/applications/solvers/electromagnetics/magneticFoam/Make/options
@ -3,5 +3,7 @@ EXE_INC = \
    -I$(LIB_SRC)/meshTools/lnInclude

 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -lfiniteVolume \
    -lmeshTools
--- a/applications/solvers/electromagnetics/mhdFoam/Make/options
+++ b/applications/solvers/electromagnetics/mhdFoam/Make/options
@ -3,5 +3,7 @@ EXE_INC = \
    -I$(LIB_SRC)/meshTools/lnInclude

 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -lfiniteVolume \
    -lmeshTools
--- a/applications/solvers/financial/financialFoam/Make/options
+++ b/applications/solvers/financial/financialFoam/Make/options
@ -4,6 +4,8 @@ EXE_INC = \
    -I$(LIB_SRC)/sampling/lnInclude

 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -lfiniteVolume \
    -lmeshTools \
    -lsampling
--- a/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/Make/options
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/Make/options
@ -11,6 +11,8 @@ EXE_INC = \
    -I$(LIB_SRC)/finiteVolume/lnInclude

 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -lturbulenceModels \
    -lincompressibleTurbulenceModels \
    -lincompressibleTransportModels \
--- a/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/Make/options
+++ b/applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/Make/options
@ -9,6 +9,8 @@ EXE_INC = \
    -I$(LIB_SRC)/meshTools/lnInclude \

 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -lturbulenceModels \
    -lincompressibleTurbulenceModels \
    -lincompressibleTransportModels \
--- a/applications/solvers/heatTransfer/buoyantPimpleFoam/Make/options
+++ b/applications/solvers/heatTransfer/buoyantPimpleFoam/Make/options
@ -9,6 +9,8 @@ EXE_INC = \
    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude

 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -lfiniteVolume \
    -lsampling \
    -lmeshTools \
--- a/applications/solvers/heatTransfer/buoyantSimpleFoam/Make/options
+++ b/applications/solvers/heatTransfer/buoyantSimpleFoam/Make/options
@ -9,6 +9,8 @@ EXE_INC = \
    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude

 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -lfiniteVolume \
    -lfvOptions \
    -lsampling \
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/Make/options
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/Make/options
@ -24,6 +24,8 @@ EXE_INC = \


 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -lcompressibleTransportModels \
    -lfluidThermophysicalModels \
    -lspecie \
--- a/applications/solvers/heatTransfer/thermoFoam/Make/options
+++ b/applications/solvers/heatTransfer/thermoFoam/Make/options
@ -9,6 +9,8 @@ EXE_INC = \
    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude

 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -lfiniteVolume \
    -lsampling \
    -lmeshTools \
--- a/applications/solvers/incompressible/adjointShapeOptimizationFoam/Make/options
+++ b/applications/solvers/incompressible/adjointShapeOptimizationFoam/Make/options
@ -10,6 +10,8 @@ EXE_INC = \
    -I$(LIB_SRC)/meshTools/lnInclude \

 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -lturbulenceModels \
    -lincompressibleTurbulenceModels \
    -lincompressibleTransportModels \
--- a/applications/solvers/incompressible/boundaryFoam/Make/options
+++ b/applications/solvers/incompressible/boundaryFoam/Make/options
@ -8,6 +8,8 @@ EXE_INC = \
    -I$(LIB_SRC)/sampling/lnInclude

 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -lturbulenceModels \
    -lincompressibleTurbulenceModels \
    -lincompressibleTransportModels \
--- a/applications/solvers/incompressible/icoFoam/Make/options
+++ b/applications/solvers/incompressible/icoFoam/Make/options
@ -3,5 +3,7 @@ EXE_INC = \
    -I$(LIB_SRC)/meshTools/lnInclude

 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -lfiniteVolume \
    -lmeshTools
--- a/applications/solvers/incompressible/nonNewtonianIcoFoam/Make/options
+++ b/applications/solvers/incompressible/nonNewtonianIcoFoam/Make/options
@ -5,6 +5,8 @@ EXE_INC = \
    -I$(LIB_SRC)/meshTools/lnInclude

 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -lincompressibleTransportModels \
    -lfiniteVolume \
    -lmeshTools
--- a/applications/solvers/incompressible/pimpleFoam/Make/options
+++ b/applications/solvers/incompressible/pimpleFoam/Make/options
@ -10,6 +10,8 @@ EXE_INC = \
    -I$(LIB_SRC)/meshTools/lnInclude \

 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -lturbulenceModels \
    -lincompressibleTurbulenceModels \
    -lincompressibleTransportModels \
--- a/applications/solvers/incompressible/pimpleFoam/SRFPimpleFoam/Make/options
+++ b/applications/solvers/incompressible/pimpleFoam/SRFPimpleFoam/Make/options
@ -10,6 +10,8 @@ EXE_INC = \
    -I$(LIB_SRC)/sampling/lnInclude

 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -lturbulenceModels \
    -lincompressibleTurbulenceModels \
    -lincompressibleTransportModels \
--- a/applications/solvers/incompressible/pisoFoam/Make/options
+++ b/applications/solvers/incompressible/pisoFoam/Make/options
@ -8,6 +8,8 @@ EXE_INC = \
    -I$(LIB_SRC)/sampling/lnInclude

 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -lturbulenceModels \
    -lincompressibleTurbulenceModels \
    -lincompressibleTransportModels \
--- a/applications/solvers/incompressible/shallowWaterFoam/Make/options
+++ b/applications/solvers/incompressible/shallowWaterFoam/Make/options
@ -3,5 +3,7 @@ EXE_INC = \
    -I$(LIB_SRC)/meshTools/lnInclude

 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -lfiniteVolume \
    -lmeshTools
--- a/applications/solvers/incompressible/simpleFoam/Make/options
+++ b/applications/solvers/incompressible/simpleFoam/Make/options
@ -9,6 +9,8 @@ EXE_INC = \


 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -lturbulenceModels \
    -lincompressibleTurbulenceModels \
    -lincompressibleTransportModels \
--- a/applications/solvers/incompressible/simpleFoam/SRFSimpleFoam/Make/options
+++ b/applications/solvers/incompressible/simpleFoam/SRFSimpleFoam/Make/options
@ -10,6 +10,8 @@ EXE_INC = \


 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -lturbulenceModels \
    -lincompressibleTurbulenceModels \
    -lincompressibleTransportModels \
--- a/applications/solvers/incompressible/simpleFoam/porousSimpleFoam/Make/options
+++ b/applications/solvers/incompressible/simpleFoam/porousSimpleFoam/Make/options
@ -10,6 +10,8 @@ EXE_INC = \


 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -lturbulenceModels \
    -lincompressibleTurbulenceModels \
    -lincompressibleTransportModels \
--- a/applications/solvers/lagrangian/DPMFoam/DPMDyMFoam/MPPICDyMFoam/Make/options
+++ b/applications/solvers/lagrangian/DPMFoam/DPMDyMFoam/MPPICDyMFoam/Make/options
@ -22,6 +22,8 @@ EXE_INC = \
    -I$(LIB_SRC)/dynamicMesh/lnInclude

 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -lfiniteVolume \
    -lfvOptions \
    -lmeshTools \
--- a/applications/solvers/lagrangian/DPMFoam/DPMDyMFoam/Make/options
+++ b/applications/solvers/lagrangian/DPMFoam/DPMDyMFoam/Make/options
@ -21,6 +21,8 @@ EXE_INC = \
    -I$(LIB_SRC)/dynamicMesh/lnInclude

 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -llagrangian \
    -llagrangianIntermediate \
    -llagrangianTurbulence \
--- a/applications/solvers/lagrangian/DPMFoam/MPPICFoam/Make/options
+++ b/applications/solvers/lagrangian/DPMFoam/MPPICFoam/Make/options
@ -19,6 +19,8 @@ EXE_INC = \
    -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \

 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -lfiniteVolume \
    -lfvOptions \
    -lmeshTools \
--- a/applications/solvers/lagrangian/DPMFoam/Make/options
+++ b/applications/solvers/lagrangian/DPMFoam/Make/options
@ -18,6 +18,8 @@ EXE_INC = \
    -I$(LIB_SRC)/meshTools/lnInclude

 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -llagrangian \
    -llagrangianIntermediate \
    -llagrangianTurbulence \
--- a/applications/solvers/lagrangian/coalChemistryFoam/Make/options
+++ b/applications/solvers/lagrangian/coalChemistryFoam/Make/options
@ -23,6 +23,8 @@ EXE_INC = \
    -I$(LIB_SRC)/sampling/lnInclude

 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -lfiniteVolume \
    -lmeshTools \
    -lturbulenceModels \
--- a/applications/solvers/lagrangian/icoUncoupledKinematicParcelFoam/Make/options
+++ b/applications/solvers/lagrangian/icoUncoupledKinematicParcelFoam/Make/options
@ -16,6 +16,8 @@ EXE_INC = \
    -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude

 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -llagrangian \
    -llagrangianIntermediate \
    -llagrangianTurbulence \
--- a/applications/solvers/lagrangian/icoUncoupledKinematicParcelFoam/icoUncoupledKinematicParcelDyMFoam/Make/options
+++ b/applications/solvers/lagrangian/icoUncoupledKinematicParcelFoam/icoUncoupledKinematicParcelDyMFoam/Make/options
@ -20,6 +20,8 @@ EXE_INC = \


 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -llagrangian \
    -llagrangianIntermediate \
    -llagrangianTurbulence \
--- a/applications/solvers/lagrangian/reactingParcelFoam/Make/options
+++ b/applications/solvers/lagrangian/reactingParcelFoam/Make/options
@ -24,6 +24,8 @@ EXE_INC = \
    -I$(FOAM_SOLVERS)/combustion/reactingFoam

 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -lfiniteVolume \
    -lfvOptions \
    -lsampling \
--- a/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/Make/options
+++ b/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/Make/options
@ -25,6 +25,8 @@ EXE_INC = \


 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -lfiniteVolume \
    -lmeshTools \
    -lturbulenceModels \
--- a/applications/solvers/lagrangian/sprayFoam/Make/options
+++ b/applications/solvers/lagrangian/sprayFoam/Make/options
@ -24,6 +24,8 @@ EXE_INC = \
    -I$(LIB_SRC)/combustionModels/lnInclude

 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -lturbulenceModels \
    -lcompressibleTurbulenceModels \
    -llagrangian \
--- a/applications/solvers/lagrangian/sprayFoam/engineFoam/Make/options
+++ b/applications/solvers/lagrangian/sprayFoam/engineFoam/Make/options
@ -28,6 +28,8 @@ EXE_INC = \


 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -lfiniteVolume \
    -lmeshTools \
    -lsampling \
--- a/applications/solvers/lagrangian/sprayFoam/sprayDyMFoam/Make/options
+++ b/applications/solvers/lagrangian/sprayFoam/sprayDyMFoam/Make/options
@ -28,6 +28,8 @@ EXE_INC = \
    -I$(LIB_SRC)/meshTools/lnInclude

 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -lturbulenceModels \
    -lcompressibleTurbulenceModels \
    -llagrangian \
--- a/applications/solvers/lagrangian/uncoupledKinematicParcelFoam/Make/options
+++ b/applications/solvers/lagrangian/uncoupledKinematicParcelFoam/Make/options
@ -14,6 +14,8 @@ EXE_INC = \
    -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude

 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -llagrangian \
    -llagrangianIntermediate \
    -llagrangianTurbulence \
--- a/applications/solvers/lagrangian/uncoupledKinematicParcelFoam/uncoupledKinematicParcelDyMFoam/Make/options
+++ b/applications/solvers/lagrangian/uncoupledKinematicParcelFoam/uncoupledKinematicParcelDyMFoam/Make/options
@ -17,6 +17,8 @@ EXE_INC = \
    -I$(LIB_SRC)/dynamicFvMesh/lnInclude

 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -llagrangian \
    -llagrangianIntermediate \
    -llagrangianTurbulence \
--- a/applications/solvers/multiphase/cavitatingFoam/Make/options
+++ b/applications/solvers/multiphase/cavitatingFoam/Make/options
@ -10,6 +10,8 @@ EXE_INC = \
    -I$(LIB_SRC)/meshTools/lnInclude

 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -lturbulenceModels \
    -lincompressibleTurbulenceModels \
    -lbarotropicCompressibilityModel \
--- a/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/Make/options
+++ b/applications/solvers/multiphase/cavitatingFoam/cavitatingDyMFoam/Make/options
@ -13,6 +13,8 @@ EXE_INC = \
    -I$(LIB_SRC)/dynamicFvMesh/lnInclude

 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -lturbulenceModels \
    -lincompressibleTurbulenceModels \
    -lbarotropicCompressibilityModel \
--- a/applications/solvers/multiphase/compressibleInterFoam/Make/options
+++ b/applications/solvers/multiphase/compressibleInterFoam/Make/options
@ -14,6 +14,8 @@ EXE_INC = \
    -I$(LIB_SRC)/meshTools/lnInclude

 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -ltwoPhaseMixtureThermo \
    -ltwoPhaseSurfaceTension \
    -lcompressibleTransportModels \
--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/Make/options
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/Make/options
@ -17,6 +17,8 @@ EXE_INC = \
    -I$(LIB_SRC)/finiteVolume/lnInclude

 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -ltwoPhaseMixtureThermo \
    -ltwoPhaseSurfaceTension \
    -lcompressibleTransportModels \
--- a/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/Make/options
+++ b/applications/solvers/multiphase/compressibleInterFoam/compressibleInterFilmFoam/Make/options
@ -27,6 +27,8 @@ EXE_INC = \
    -I$(LIB_SRC)/fvOptions/lnInclude

 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -ltwoPhaseMixtureThermo \
    -ltwoPhaseSurfaceTension \
    -lcompressibleTransportModels \
--- a/applications/solvers/multiphase/compressibleMultiphaseInterFoam/Make/options
+++ b/applications/solvers/multiphase/compressibleMultiphaseInterFoam/Make/options
@ -13,6 +13,8 @@ EXE_INC = \
    -I$(LIB_SRC)/meshTools/lnInclude

 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -lmultiphaseMixtureThermo \
    -lcompressibleTransportModels \
    -lfluidThermophysicalModels \
--- a/applications/solvers/multiphase/driftFluxFoam/Make/options
+++ b/applications/solvers/multiphase/driftFluxFoam/Make/options
@ -12,6 +12,8 @@ EXE_INC = \
    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude

 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -ldriftFluxTransportModels \
    -ldriftFluxRelativeVelocityModels \
    -lfiniteVolume \
--- a/applications/solvers/multiphase/interFoam/Make/options
+++ b/applications/solvers/multiphase/interFoam/Make/options
@ -13,6 +13,8 @@ EXE_INC = \
    -I$(LIB_SRC)/sampling/lnInclude

 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -limmiscibleIncompressibleTwoPhaseMixture \
    -lturbulenceModels \
    -lincompressibleTurbulenceModels \
--- a/applications/solvers/multiphase/interFoam/interMixingFoam/Make/options
+++ b/applications/solvers/multiphase/interFoam/interMixingFoam/Make/options
@ -17,6 +17,8 @@ EXE_INC = \
    -I$(LIB_SRC)/sampling/lnInclude

 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -ltwoPhaseMixture \
    -ltwoPhaseProperties \
    -lincompressibleTransportModels \
--- a/applications/solvers/multiphase/interPhaseChangeFoam/Make/options
+++ b/applications/solvers/multiphase/interPhaseChangeFoam/Make/options
@ -11,6 +11,8 @@ EXE_INC = \
    -I$(LIB_SRC)/sampling/lnInclude

 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -lphaseChangeTwoPhaseMixtures \
    -ltwoPhaseMixture \
    -linterfaceProperties \
--- a/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/Make/options
+++ b/applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeDyMFoam/Make/options
@ -14,6 +14,8 @@ EXE_INC = \
    -I$(LIB_SRC)/sampling/lnInclude

 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -lphaseChangeTwoPhaseMixtures \
    -ltwoPhaseMixture \
    -linterfaceProperties \
--- a/applications/solvers/multiphase/multiphaseEulerFoam/Make/options
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/Make/options
@ -13,6 +13,8 @@ EXE_INC = \


 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -lmultiphaseSystem \
    -linterfaceProperties \
    -lincompressibleTransportModels \
--- a/applications/solvers/multiphase/multiphaseInterFoam/Make/options
+++ b/applications/solvers/multiphase/multiphaseInterFoam/Make/options
@ -14,6 +14,8 @@ EXE_INC = \
    -I$(LIB_SRC)/sampling/lnInclude

 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -lmultiphaseInterFoam \
    -linterfaceProperties \
    -lincompressibleTransportModels \
--- a/applications/solvers/multiphase/potentialFreeSurfaceFoam/Make/options
+++ b/applications/solvers/multiphase/potentialFreeSurfaceFoam/Make/options
@ -9,6 +9,8 @@ EXE_INC = \


 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -lincompressibleTransportModels \
    -lturbulenceModels \
    -lincompressibleTurbulenceModels \
--- a/applications/solvers/multiphase/potentialFreeSurfaceFoam/potentialFreeSurfaceDyMFoam/Make/options
+++ b/applications/solvers/multiphase/potentialFreeSurfaceFoam/potentialFreeSurfaceDyMFoam/Make/options
@ -13,6 +13,8 @@ EXE_INC = \


 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -lincompressibleTransportModels \
    -lturbulenceModels \
    -lincompressibleTurbulenceModels \
--- a/applications/solvers/multiphase/reactingEulerFoam/interfacialCompositionModels/Make/options
+++ b/applications/solvers/multiphase/reactingEulerFoam/interfacialCompositionModels/Make/options
@ -12,9 +12,19 @@ EXE_INC = \
    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/phaseCompressible/lnInclude \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(WM_CANTERA_DIR) \
+    -I$(WM_CANTERA_DIR)/cantera \
+    -I$(WM_CANTERA_DIR)/cantera/kinetics \
+    -I$(WM_CANTERA_DIR)/cantera/base \
+    -I$(WM_CANTERA_DIR)/cantera/transport \
+    -I$(WM_CANTERA_DIR)/cantera/thermo \
+    -I$(WM_CANTERA_DIR)/cantera/ext \
+    -I$(WM_CANTERA_DIR)/cantera/ext/fmt

 LIB_LIBS = \
    -lfluidThermophysicalModels \
    -lreactionThermophysicalModels \
-    -lspecie
+    -lspecie \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared
--- a/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/Make/options
+++ b/applications/solvers/multiphase/reactingEulerFoam/interfacialModels/Make/options
@ -7,9 +7,19 @@ EXE_INC = \
    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
    -I$(LIB_SRC)/TurbulenceModels/phaseCompressible/lnInclude \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(WM_CANTERA_DIR) \
+    -I$(WM_CANTERA_DIR)/cantera \
+    -I$(WM_CANTERA_DIR)/cantera/kinetics \
+    -I$(WM_CANTERA_DIR)/cantera/base \
+    -I$(WM_CANTERA_DIR)/cantera/transport \
+    -I$(WM_CANTERA_DIR)/cantera/thermo \
+    -I$(WM_CANTERA_DIR)/cantera/ext \
+    -I$(WM_CANTERA_DIR)/cantera/ext/fmt

 LIB_LIBS = \
    -lcompressibleTransportModels \
    -lfluidThermophysicalModels \
-    -lspecie
+    -lspecie \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared
--- a/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/Make/options
+++ b/applications/solvers/multiphase/reactingEulerFoam/phaseSystems/Make/options
@ -12,7 +12,18 @@ EXE_INC = \
    -I$(LIB_SRC)/combustionModels/lnInclude \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
    -I$(LIB_SRC)/meshTools/lnInclude \
-    -I$(LIB_SRC)/sampling/lnInclude
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(WM_CANTERA_DIR) \
+    -I$(WM_CANTERA_DIR)/cantera \
+    -I$(WM_CANTERA_DIR)/cantera/kinetics \
+    -I$(WM_CANTERA_DIR)/cantera/base \
+    -I$(WM_CANTERA_DIR)/cantera/transport \
+    -I$(WM_CANTERA_DIR)/cantera/thermo \
+    -I$(WM_CANTERA_DIR)/cantera/ext \
+    -I$(WM_CANTERA_DIR)/cantera/ext/fmt

 LIB_LIBS = \
-    -lcombustionModels
+    -L$(FOAM_USER_LIBBIN) \
+    -lcombustionModels \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/Make/options
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingMultiphaseEulerFoam/Make/options
@ -14,6 +14,8 @@ EXE_INC = \
    -I$(LIB_SRC)/sampling/lnInclude

 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -lreactingPhaseSystem \
    -lreactingMultiphaseSystem \
    -lreactingEulerianInterfacialModels \
--- a/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/Make/options
+++ b/applications/solvers/multiphase/reactingEulerFoam/reactingTwoPhaseEulerFoam/Make/options
@ -14,6 +14,8 @@ EXE_INC = \
    -I$(LIB_SRC)/sampling/lnInclude

 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -lreactingPhaseSystem \
    -lreactingTwoPhaseSystem \
    -lreactingEulerianInterfacialModels \
--- a/applications/solvers/multiphase/twoLiquidMixingFoam/Make/options
+++ b/applications/solvers/multiphase/twoLiquidMixingFoam/Make/options
@ -11,6 +11,8 @@ EXE_INC = \
    -I$(LIB_SRC)/meshTools/lnInclude

 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -ltwoPhaseMixture \
    -lturbulenceModels \
    -lincompressibleTurbulenceModels \
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/Make/options
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/Make/options
@ -13,6 +13,8 @@ EXE_INC = \
    -I$(LIB_SRC)/sampling/lnInclude

 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -lcompressibleTransportModels \
    -lfluidThermophysicalModels \
    -lspecie \
--- a/applications/solvers/stressAnalysis/solidDisplacementFoam/Make/options
+++ b/applications/solvers/stressAnalysis/solidDisplacementFoam/Make/options
@ -5,5 +5,7 @@ EXE_INC = \
    -I$(LIB_SRC)/meshTools/lnInclude

 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -lfiniteVolume \
    -lmeshTools
--- a/applications/solvers/stressAnalysis/solidEquilibriumDisplacementFoam/Make/options
+++ b/applications/solvers/stressAnalysis/solidEquilibriumDisplacementFoam/Make/options
@ -6,5 +6,7 @@ EXE_INC = \
    -I$(LIB_SRC)/meshTools/lnInclude

 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -lfiniteVolume \
    -lmeshTools
--- a/applications/test/Function1/Make/options
+++ b/applications/test/Function1/Make/options
@ -6,6 +6,8 @@ EXE_INC = \
    -I$(LIB_SRC)/meshTools/lnInclude

 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -llagrangianIntermediate \
    -lradiationModels \
    -lregionModels \
--- a/applications/test/GAMGAgglomeration/Make/options
+++ b/applications/test/GAMGAgglomeration/Make/options
@ -3,5 +3,7 @@ EXE_INC = \
    -I$(LIB_SRC)/meshTools/lnInclude

 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -lfiniteVolume \
    -lmeshTools
--- a/applications/test/PatchEdgeFaceWave/Make/options
+++ b/applications/test/PatchEdgeFaceWave/Make/options
@ -3,5 +3,7 @@ EXE_INC = \
    -I$(LIB_SRC)/meshTools/lnInclude

 EXE_LIBS = \
+    -L$(CANTERA_LIB_PATH) \
+    -lcantera_shared \
    -lfiniteVolume \
    -lmeshTools
--- a/Show More
+++ b/Show More