adapt solvers with dfCombustionModel

applications/solvers/dfLowMachFoam/Make/options

@@ -14,6 +14,7 @@ EXE_INC = -std=c++14 \
     -I$(LIB_SRC)/Pstream/mpi \
     -I$(DF_SRC)/CanteraMixture/lnInclude \
     -I$(DF_SRC)/dfChemistryModel/lnInclude \
+    -I$(DF_SRC)/dfCombustionModels/lnInclude \
     -I$(CANTERA_ROOT)/include \
     -I$(TORCH_ROOT)/include \
     -I$(TORCH_ROOT)/include/torch/csrc/api/include
@@ -29,6 +30,7 @@ EXE_LIBS = \
     -ldfCompressibleTurbulenceModels \
     -lCanteraMixture \
     -ldfChemistryModel \
+    -ldfCombustionModels  \
     $(CANTERA_ROOT)/lib/libcantera.so \
     $(TORCH_ROOT)/lib/libtorch.so \
     $(TORCH_ROOT)/lib/libc10.so \

applications/solvers/dfLowMachFoam/YEqn.H

@@ -15,13 +15,13 @@ tmp<fv::convectionScheme<scalar>> mvConvection
 
 forAll(Y, i)
 {
-    Derror += chemistry.rhoD(i)*fvc::grad(Y[i]);
+    Derror += chemistry->rhoD(i)*fvc::grad(Y[i]);
 }
 const surfaceScalarField phiUc = linearInterpolate(Derror) & mesh.Sf();
 
 {
     start = std::clock();
-    chemistry.solve(mesh.time().deltaTValue());
+    combustion->correct();
     label flag_mpi_init;
     MPI_Initialized(&flag_mpi_init);
     if(flag_mpi_init) MPI_Barrier(PstreamGlobals::MPI_COMM_FOAM);
@@ -34,20 +34,20 @@ const surfaceScalarField phiUc = linearInterpolate(Derror) & mesh.Sf();
     forAll(Y, i)
     {
         volScalarField& Yi = Y[i];
-        vk += chemistry.hai(i)*(chemistry.rhoD(i)*fvc::grad(Yi) - Yi*Derror);
-        diffsh += fvc::laplacian(thermo.alpha()*chemistry.hai(i), Yi);
+        vk += chemistry->hai(i)*(chemistry->rhoD(i)*fvc::grad(Yi) - Yi*Derror);
+        diffsh += fvc::laplacian(thermo.alpha()*chemistry->hai(i), Yi);
 
         if (i != inertIndex)
         {
-            tmp<volScalarField> DEff = chemistry.rhoD(i) + turbulence->mut()/Sct;
+            tmp<volScalarField> DEff = chemistry->rhoD(i) + turbulence->mut()/Sct;
             fvScalarMatrix YiEqn
             (
                 fvm::ddt(rho, Yi)
               + mvConvection->fvmDiv(phi, Yi)
               - fvm::laplacian(DEff(), Yi)
               + mvConvection->fvmDiv(phiUc, Yi)
              ==
-                chemistry.RR(i)
+                combustion->R(Yi)
             );
 
             YiEqn.relax();

applications/solvers/dfLowMachFoam/createFields.H

@@ -22,14 +22,6 @@ volScalarField rho
     thermo.rho()
 );
 
-dfChemistryModel<basicThermo> chemistry(thermo);
-PtrList<volScalarField>& Y = chemistry.Y();
-const word inertSpecie(chemistry.lookup("inertSpecie"));
-const label inertIndex(chemistry.species()[inertSpecie]);
-
-chemistry.correctThermo();
-Info<< "At initial time, min/max(T) = " << min(T).value() << ", " << max(T).value() << endl;
-
 Info<< "Reading field U\n" << endl;
 volVectorField U
 (
@@ -62,6 +54,21 @@ autoPtr<compressible::turbulenceModel> turbulence
     )
 );
 
+Info<< "Creating reaction model\n" << endl;
+autoPtr<CombustionModel<basicThermo>> combustion
+(
+    CombustionModel<basicThermo>::New(thermo, turbulence())
+);
+Info<< "end Creating reaction model\n" << endl;
+dfChemistryModel<basicThermo>* chemistry = combustion->chemistry();
+PtrList<volScalarField>& Y = chemistry->Y();
+const word inertSpecie(chemistry->lookup("inertSpecie"));
+const label inertIndex(chemistry->species()[inertSpecie]);
+
+
+chemistry->correctThermo();
+Info<< "At initial time, min/max(T) = " << min(T).value() << ", " << max(T).value() << endl;
+
 Info<< "Creating field dpdt\n" << endl;
 volScalarField dpdt
 (
@@ -85,7 +92,7 @@ forAll(Y, i)
 }
 fields.add(thermo.he());
 
-const scalar Sct = chemistry.lookupOrDefault("Sct", 1.);
+const scalar Sct = chemistry->lookupOrDefault("Sct", 1.);
 volScalarField diffsh
 (
     IOobject

applications/solvers/dfLowMachFoam/dfLowMachFoam.C

@@ -45,6 +45,8 @@ Description
 #include "localEulerDdtScheme.H"
 #include "fvcSmooth.H"
 #include "PstreamGlobals.H"
+#include "basicThermo.H"
+#include "CombustionModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -120,7 +122,7 @@ int main(int argc, char *argv[])
 
             #include "YEqn.H"//added by zhangyan 20220226
             #include "EEqn.H"
-            chemistry.correctThermo();
+            chemistry->correctThermo();
             Info<< "min/max(T) = " << min(T).value() << ", " << max(T).value() << endl;
 
             // --- Pressure corrector loop

applications/solvers/dfSprayFoam/Make/options

@@ -24,6 +24,7 @@ EXE_INC = -std=c++14 \
     -I$(DF_SRC)/thermophysicalModels/SLGThermo/lnInclude \
     -I$(DF_SRC)/CanteraMixture/lnInclude \
     -I$(DF_SRC)/dfChemistryModel/lnInclude \
+    -I$(DF_SRC)/dfCombustionModels/lnInclude \
     -I$(CANTERA_ROOT)/include \
     -I$(TORCH_ROOT)/include \
     -I$(TORCH_ROOT)/include/torch/csrc/api/include
@@ -49,6 +50,7 @@ EXE_LIBS = \
     -ldfLagrangianSpray \
     -lCanteraMixture \
     -ldfChemistryModel \
+    -ldfCombustionModels  \
     $(CANTERA_ROOT)/lib/libcantera.so \
     $(TORCH_ROOT)/lib/libtorch.so \
     $(TORCH_ROOT)/lib/libc10.so \

applications/solvers/dfSprayFoam/YEqn.H

@@ -10,8 +10,7 @@ tmp<fv::convectionScheme<scalar>> mvConvection
 );
 
 {
-    //combustion->correct();
-    chemistry.solve(mesh.time().deltaTValue());
+    combustion->correct();
     volScalarField Yt(0.0*Y[0]);
 
     forAll(Y, i)
@@ -27,7 +26,7 @@ tmp<fv::convectionScheme<scalar>> mvConvection
               - fvm::laplacian(turbulence->muEff(), Yi)
              ==
                 parcels.SYi(i, Yi)
-              + chemistry.RR(i)
+              + combustion->R(Yi)
               //+ fvOptions(rho, Yi)
             );
 

applications/solvers/dfSprayFoam/createFields.H

@@ -23,14 +23,6 @@ volScalarField rho
     thermo.rho()
 );
 
-dfChemistryModel<basicThermo> chemistry(thermo);
-PtrList<volScalarField>& Y = chemistry.Y();
-const word inertSpecie(chemistry.lookup("inertSpecie"));
-const label inertIndex(chemistry.species()[inertSpecie]);
-
-chemistry.correctThermo();
-Info<< "At initial time, min/max(T) = " << min(T).value() << ", " << max(T).value() << endl;
-
 Info<< "\nReading field U\n" << endl;
 volVectorField U
 (
@@ -83,6 +75,22 @@ autoPtr<compressible::turbulenceModel> turbulence
     )
 );
 
+Info<< "Creating reaction model\n" << endl;
+autoPtr<CombustionModel<basicThermo>> combustion
+(
+    CombustionModel<basicThermo>::New(thermo, turbulence())
+);
+Info<< "end Creating reaction model\n" << endl;
+dfChemistryModel<basicThermo>* chemistry = combustion->chemistry();
+PtrList<volScalarField>& Y = chemistry->Y();
+const word inertSpecie(chemistry->lookup("inertSpecie"));
+const label inertIndex(chemistry->species()[inertSpecie]);
+
+
+chemistry->correctThermo();
+Info<< "At initial time, min/max(T) = " << min(T).value() << ", " << max(T).value() << endl;
+
+
 Info<< "Creating field dpdt\n" << endl;
 volScalarField dpdt
 (

applications/solvers/dfSprayFoam/dfSprayFoam.C

@@ -43,6 +43,8 @@ Description
 #include "pimpleControl.H"
 #include "CorrectPhi.H"
 //#include "fvOptions.H"
+#include "basicThermo.H"
+#include "CombustionModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -134,7 +136,7 @@ int main(int argc, char *argv[])
             #include "UEqn.H"
             #include "YEqn.H"
             #include "EEqn.H"
-            chemistry.correctThermo();
+            chemistry->correctThermo();
             Info<< "T gas min/max   " << min(T).value() << ", "
                 << max(T).value() << endl;