force reading combustionProperties, update dfHighSpeedFoam and examples

applications/solvers/dfHighSpeedFoam/Make/options

@@ -17,6 +17,7 @@ EXE_INC = -std=c++14 \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(DF_SRC)/CanteraMixture/lnInclude \
     -I$(DF_SRC)/dfChemistryModel/lnInclude \
+    -I$(DF_SRC)/dfCombustionModels/lnInclude \
     -I$(CANTERA_ROOT)/include \
     $(if $(LIBTORCH_ROOT),-I$(LIBTORCH_ROOT)/include,) \
     $(if $(LIBTORCH_ROOT),-I$(LIBTORCH_ROOT)/include/torch/csrc/api/include,) \
@@ -34,6 +35,7 @@ EXE_LIBS = \
     -ldfCompressibleTurbulenceModels \
     -lCanteraMixture \
     -ldfChemistryModel \
+    -ldfCombustionModels  \
     $(CANTERA_ROOT)/lib/libcantera.so \
     $(if $(LIBTORCH_ROOT),$(LIBTORCH_ROOT)/lib/libtorch.so,) \
     $(if $(LIBTORCH_ROOT),$(LIBTORCH_ROOT)/lib/libc10.so,) \

applications/solvers/dfHighSpeedFoam/createFields.H

@@ -125,12 +125,18 @@ autoPtr<compressible::turbulenceModel> turbulence
 
 multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
 
-dfChemistryModel<basicThermo> chemistry(thermo);
-PtrList<volScalarField>& Y = chemistry.Y();
-const word inertSpecie(chemistry.lookup("inertSpecie"));
-const label inertIndex(chemistry.species()[inertSpecie]);
+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();
+chemistry->correctThermo();
 Info<< "At initial time, min/max(T) = " << min(T).value() << ", " << max(T).value() << endl;
 
 forAll(Y, i)
@@ -139,7 +145,7 @@ forAll(Y, i)
 }
 fields.add(thermo.he());
 
-const scalar Sct = chemistry.lookupOrDefault("Sct", 1.);
+const scalar Sct = chemistry->lookupOrDefault("Sct", 1.);
 volScalarField diffAlphaD
 (
     IOobject

applications/solvers/dfHighSpeedFoam/dfHighSpeedFoam.C

@@ -50,6 +50,7 @@ Description
 #include "localEulerDdtScheme.H"
 #include "fvcSmooth.H"
 #include "PstreamGlobals.H"
+#include "CombustionModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

applications/solvers/dfHighSpeedFoam/rhoEEqn.H

@@ -19,7 +19,7 @@ e = rhoE/rho - 0.5*magSqr(U);
 e.correctBoundaryConditions();
 
 ha = e + p/rho;
-chemistry.correctThermo(); // before this, we must update ha = e + p/rho
+chemistry->correctThermo(); // before this, we must update ha = e + p/rho
 
 rhoE.boundaryFieldRef() == rho.boundaryField()*(e.boundaryField() + 0.5*magSqr(U.boundaryField()));
 
@@ -38,7 +38,7 @@ if (!inviscid)
     eEqn.solve("e");
 
     ha = e + p/rho;
-    chemistry.correctThermo();
+    chemistry->correctThermo();
     rhoE = rho*(e + 0.5*magSqr(U));
 }
 

applications/solvers/dfHighSpeedFoam/rhoYEqn.H

@@ -6,7 +6,7 @@ if (!inviscid)
 
     forAll(Y, i)
     {
-        sumYDiffError += chemistry.rhoD(i)*fvc::grad(Y[i]);
+        sumYDiffError += chemistry->rhoD(i)*fvc::grad(Y[i]);
     }
 }
 
@@ -23,14 +23,14 @@ tmp<fv::convectionScheme<scalar>> mvConvection
 
 {
     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);
     end = std::clock();
     time_monitor_chem += double(end - start) / double(CLOCKS_PER_SEC);
 
-    volScalarField Yt(0*Y[0]);
+    volScalarField Yt(0.0*Y[0]);
 
     start = std::clock();
     forAll(Y, i)
@@ -44,16 +44,16 @@ tmp<fv::convectionScheme<scalar>> mvConvection
                  fvm::ddt(rho, Yi)
                + mvConvection->fvmDiv(phi, Yi)
               ==
-                 chemistry.RR(i)
+                combustion->R(Yi)
             );
 
             if (!inviscid)
             {
                 const surfaceScalarField phiUc = linearInterpolate(sumYDiffError) & mesh.Sf();
 
-                hDiffCorrFlux += chemistry.hai(i)*(chemistry.rhoD(i)*fvc::grad(Yi) - Yi*sumYDiffError);
-                diffAlphaD += fvc::laplacian(thermo.alpha()*chemistry.hai(i), Yi);
-                tmp<volScalarField> DEff = chemistry.rhoD(i) + turbulence->mut()/Sct;
+                hDiffCorrFlux += chemistry->hai(i)*(chemistry->rhoD(i)*fvc::grad(Yi) - Yi*sumYDiffError);
+                diffAlphaD += fvc::laplacian(thermo.alpha()*chemistry->hai(i), Yi);
+                tmp<volScalarField> DEff = chemistry->rhoD(i) + turbulence->mut()/Sct;
 
                 YiEqn -= fvm::laplacian(DEff(), Yi) - mvConvection->fvmDiv(phiUc, Yi);
             }

examples/dfHighSpeedFoam/oneD_detonationH2/constant/combustionProperties

@@ -0,0 +1,20 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  7
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      combustionProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+combustionModel  laminar;
+
+// ************************************************************************* //

examples/dfHighSpeedFoam/oneD_reactiveShockTube/0_orig/T

@@ -28,8 +28,7 @@ boundaryField
 
     Right
     {
-        type            inletOutlet;
-        inletValue      uniform 748.472;
+        type            fixedValue;
         value           uniform 748.472;
     }
 

examples/dfHighSpeedFoam/oneD_reactiveShockTube/constant/combustionProperties

@@ -0,0 +1,20 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  7
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      combustionProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+combustionModel  laminar;
+
+// ************************************************************************* //

examples/dfHighSpeedFoam/twoD_detonationH2/constant/combustionProperties

@@ -0,0 +1,20 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  7
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      combustionProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+combustionModel  laminar;
+
+// ************************************************************************* //

src/dfCombustionModels/combustionModel/combustionModelTemplates.C

@@ -46,17 +46,7 @@ Foam::autoPtr<CombustionModel> Foam::combustionModel::New
         )
     );
 
-    word combModelName("none");
-    if (combIO.typeHeaderOk<IOdictionary>(false))
-    {
-        IOdictionary(combIO).lookup("combustionModel") >> combModelName;
-    }
-    else
-    {
-        Info<< "Combustion model not active: "
-            << thermo.phasePropertyName(combustionProperties)
-            << " not found" << endl;
-    }
+    word combModelName = IOdictionary(combIO).lookup("combustionModel");
 
     Info<< "Selecting combustion model " << combModelName << endl;