Merge branch 'master' into GPU

Author: maorz1998

.github/workflows/CPU_inferencce_validation.yml

@@ -50,8 +50,8 @@ jobs:
         && ls $PWD
         && ls flareFGM_Table_Download
         && cp -r flareFGM_Table_Download/SandiaD/flare.tbl examples/dfLowMachFoam/2DSandiaD_flareFGM/
-        && git clone https://github.com/deepcombustion/deepcombustion.git 
-        && cp -r deepcombustion/DeePCK/Model/HE04_Hydrogen_ESH2_GMS_sub_20221101/ mechanisms/ && source ~/miniconda3/etc/profile.d/conda.sh && conda activate libcantera && source /opt/openfoam7/etc/bashrc 
+        && git clone https://github.com/intelligent-algorithm-team/intelligent-combustion.git 
+        && cp -r intelligent-combustion/DeePCK/Model/HE04_Hydrogen_ESH2_GMS_sub_20221101/ mechanisms/ && source ~/miniconda3/etc/profile.d/conda.sh && conda activate libcantera && source /opt/openfoam7/etc/bashrc 
         && . configure.sh --use_pytorch && source ./bashrc && . install.sh 
         && cd test && ./Allrun && conda deactivate "
 

.gitmodules

@@ -0,0 +1,6 @@
+[submodule "submodules/libROUNDSchemes"]
+	path = submodules/libROUNDSchemes
+	url = https://github.com/advanCFD/libROUNDSchemes.git
+[submodule "submodules/WENOEXT"]
+	path = submodules/WENOEXT
+	url = https://github.com/WENO-OF/WENOEXT.git

README.md

@@ -18,7 +18,7 @@
 
 DeepFlame is a deep learning empowered computational fluid dynamics package for single or multiphase, laminar or turbulent, reacting flows at all speeds. It aims to provide an open-source platform to combine the individual strengths of [OpenFOAM](https://openfoam.org), [Cantera](https://cantera.org), and [PyTorch](https://pytorch.org/) libraries for deep learning assisted reacting flow simulations. It also has the scope to leverage the next-generation heterogenous supercomputing and AI acceleration infrastructures such as GPU and FPGA.
 
-The deep learning algorithms and models used in the DeepFlame tutorial examples are developed and trained independently by our collaborators team – [DeepCombustion](https://github.com/deepcombustion/deepcombustion). Please refer to their website for detailed information.
+The deep learning algorithms and models used in the DeepFlame tutorial examples are developed and trained independently by our collaborators team – [Intelligent Combustion](https://github.com/intelligent-algorithm-team/intelligent-combustion.git). Please refer to their website for detailed information.
 
 ## Documentation
 Detailed guide for installation and tutorials is available on [our documentation website](https://deepflame.deepmodeling.com).

applications/solvers/df0DFoam/createFields.H

@@ -2,10 +2,11 @@
 
 Info<< "Reading thermophysical properties\n" << endl;
 
-fluidThermo* pThermo = new hePsiThermo<psiThermo, CanteraMixture>(mesh, word::null);
+// fluidThermo* pThermo = new hePsiThermo<psiThermo, CanteraMixture>(mesh, word::null);
+fluidThermo* pThermo = new heRhoThermo<rhoThermo, CanteraMixture>(mesh, word::null);
 fluidThermo& thermo = *pThermo;
 
-thermo.validate(args.executable(), "ha");
+// thermo.validate(args.executable(), "ha");
 
 volScalarField& p = thermo.p();
 
@@ -93,6 +94,8 @@ dfChemistryModel<basicThermo> chemistry(thermo);
 PtrList<volScalarField>& Y = chemistry.Y();
 const word inertSpecie(chemistry.lookup("inertSpecie"));
 const label inertIndex(chemistry.species()[inertSpecie]);
+chemistry.setEnergyName("ha");
+chemistry.updateEnergy();
 
 forAll(Y, i)
 {

applications/solvers/df0DFoam/df0DFoam.C

@@ -27,7 +27,8 @@ Description
 \*---------------------------------------------------------------------------*/
 #include "dfChemistryModel.H"
 #include "CanteraMixture.H"
-#include "hePsiThermo.H"
+// #include "hePsiThermo.H"
+#include "heRhoThermo.H"
 
 #ifdef USE_PYTORCH
 #include <pybind11/embed.h>

applications/solvers/dfHighSpeedFoam/createFields.H

@@ -2,8 +2,10 @@
 
 Info<< "Reading thermophysical properties\n" << endl;
 
-psiThermo* pThermo = new hePsiThermo<psiThermo, CanteraMixture>(mesh, word::null);
-psiThermo& thermo = *pThermo;
+// psiThermo* pThermo = new hePsiThermo<psiThermo, CanteraMixture>(mesh, word::null);
+// psiThermo& thermo = *pThermo;
+rhoThermo* pThermo = new heRhoThermo<rhoThermo, CanteraMixture>(mesh, word::null);
+rhoThermo& thermo = *pThermo;
 
 //move from creatFieldRefs.H to createFields.H
 //p needed to be created before e
@@ -27,7 +29,7 @@ dictionary thermoDict
     )
 );
 
-bool inviscid(thermoDict.lookup("inviscid"));
+bool inviscid(thermoDict.lookupOrDefault("inviscid",false));
 
 Info<< "Reading field U\n" << endl;
 volVectorField U
@@ -56,9 +58,6 @@ volScalarField rho
     thermo.rho()
 );
 
-volScalarField e = thermo.he() - p/rho;
-volScalarField& ha = thermo.he();
-
 volVectorField rhoU
 (
     IOobject
@@ -72,19 +71,6 @@ volVectorField rhoU
     rho*U
 );
 
-volScalarField rhoE
-(
-    IOobject
-    (
-        "rhoE",
-        runTime.timeName(),
-        mesh,
-        IOobject::NO_READ,
-        IOobject::NO_WRITE
-    ),
-    rho*(e + 0.5*magSqr(U))
-);
-
 surfaceScalarField pos
 (
     IOobject
@@ -135,6 +121,24 @@ dfChemistryModel<basicThermo>* chemistry = combustion->chemistry();
 PtrList<volScalarField>& Y = chemistry->Y();
 const word inertSpecie(chemistry->lookup("inertSpecie"));
 const label inertIndex(chemistry->species()[inertSpecie]);
+chemistry->setEnergyName("ea");
+chemistry->updateEnergy();
+
+volScalarField& ea = thermo.he();
+volScalarField ha = ea + p/rho;
+
+volScalarField rhoE
+(
+    IOobject
+    (
+        "rhoE",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    rho*(ea + 0.5*magSqr(U))
+);
 
 chemistry->correctThermo();
 Info<< "At initial time, min/max(T) = " << min(T).value() << ", " << max(T).value() << endl;
@@ -145,6 +149,28 @@ forAll(Y, i)
 }
 fields.add(thermo.he());
 
+const label nspecies(chemistry->species().size());
+PtrList<volScalarField> rhoYi(nspecies);
+forAll(rhoYi,i)
+{
+    rhoYi.set
+    (
+        i,
+        new volScalarField
+        (
+            IOobject
+            (
+                "rhoYi" + Y[i].name(),
+                runTime.timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            rho*Y[i]
+        )
+    );
+}
+
 const scalar Sct = chemistry->lookupOrDefault("Sct", 1.);
 volScalarField diffAlphaD
 (

applications/solvers/dfHighSpeedFoam/dfHighSpeedFoam.C

@@ -32,7 +32,8 @@ Description
 
 #include "dfChemistryModel.H"
 #include "CanteraMixture.H"
-#include "hePsiThermo.H"
+// #include "hePsiThermo.H"
+#include "heRhoThermo.H"
 
 #ifdef USE_PYTORCH
 #include <pybind11/embed.h>
@@ -47,7 +48,8 @@ Description
 
 #include "fvCFD.H"
 #include "dynamicFvMesh.H"
-#include "psiThermo.H"
+// #include "psiThermo.H"
+#include "rhoThermo.H"
 #include "turbulentFluidThermoModel.H"
 #include "fixedRhoFvPatchScalarField.H"
 #include "directionInterpolate.H"
@@ -130,15 +132,57 @@ int main(int argc, char *argv[])
         surfaceScalarField rho_pos(interpolate(rho, pos));
         surfaceScalarField rho_neg(interpolate(rho, neg));
 
+        PtrList<surfaceScalarField> rhoYi_pos(nspecies);
+        PtrList<surfaceScalarField> rhoYi_neg(nspecies);
+        forAll(rhoYi_pos,i)
+        {
+            rhoYi_pos.set
+            (
+                i,
+                new surfaceScalarField
+                (
+                    IOobject
+                    (
+                        "rhoYi_pos" + Y[i].name(),
+                        runTime.timeName(),
+                        mesh,
+                        IOobject::NO_READ,
+                        IOobject::NO_WRITE
+                    ),
+                    interpolate(rhoYi[i], pos,"Yi")
+                )
+            );
+        }
+
+        forAll(rhoYi_neg,i)
+        {
+            rhoYi_neg.set
+            (
+                i,
+                new surfaceScalarField
+                (
+                    IOobject
+                    (
+                        "rhoYi_neg" + Y[i].name(),
+                        runTime.timeName(),
+                        mesh,
+                        IOobject::NO_READ,
+                        IOobject::NO_WRITE
+                    ),
+                    interpolate(rhoYi[i], neg,"Yi")
+                )
+            );
+        }
+
         surfaceVectorField rhoU_pos(interpolate(rhoU, pos, U.name()));
         surfaceVectorField rhoU_neg(interpolate(rhoU, neg, U.name()));
 
         volScalarField rPsi("rPsi", 1.0/psi);
         surfaceScalarField rPsi_pos(interpolate(rPsi, pos, T.name()));
         surfaceScalarField rPsi_neg(interpolate(rPsi, neg, T.name()));
 
-        surfaceScalarField e_pos(interpolate(e, pos, T.name()));
-        surfaceScalarField e_neg(interpolate(e, neg, T.name()));
+        surfaceScalarField ea_pos(interpolate(ea, pos, T.name()));
+        surfaceScalarField ea_neg(interpolate(ea, neg, T.name()));
 
         surfaceVectorField U_pos("U_pos", rhoU_pos/rho_pos);
         surfaceVectorField U_neg("U_neg", rhoU_neg/rho_neg);
@@ -215,6 +259,27 @@ int main(int argc, char *argv[])
 
         phi = aphiv_pos*rho_pos + aphiv_neg*rho_neg;
 
+        PtrList<surfaceScalarField> phiYi(nspecies);
+        forAll(phiYi,i)
+        {
+            phiYi.set
+            (
+                i,
+                new surfaceScalarField
+                (
+                    IOobject
+                    (
+                        "phiYi_" + Y[i].name(),
+                        runTime.timeName(),
+                        mesh,
+                        IOobject::NO_READ,
+                        IOobject::NO_WRITE
+                    ),
+                    aphiv_pos*rhoYi_pos[i] + aphiv_neg*rhoYi_neg[i]
+                )
+            );
+        }
+
         surfaceVectorField phiUp
         (
             (aphiv_pos*rhoU_pos + aphiv_neg*rhoU_neg)
@@ -224,8 +289,8 @@ int main(int argc, char *argv[])
         surfaceScalarField phiEp
         (
             "phiEp",
-            aphiv_pos*(rho_pos*(e_pos + 0.5*magSqr(U_pos)) + p_pos)
-          + aphiv_neg*(rho_neg*(e_neg + 0.5*magSqr(U_neg)) + p_neg)
+            aphiv_pos*(rho_pos*(ea_pos + 0.5*magSqr(U_pos)) + p_pos)
+          + aphiv_neg*(rho_neg*(ea_neg + 0.5*magSqr(U_neg)) + p_neg)
           + aSf*p_pos - aSf*p_neg
         );
 

applications/solvers/dfHighSpeedFoam/rhoEEqn.H

@@ -15,31 +15,31 @@ solve
     - fvc::div(sigmaDotU)
 );
 
-e = rhoE/rho - 0.5*magSqr(U);
-e.correctBoundaryConditions();
+ea = rhoE/rho - 0.5*magSqr(U);
+ea.correctBoundaryConditions();
 
-ha = e + p/rho;
+ha = ea + p/rho;
 chemistry->correctThermo(); // before this, we must update ha = e + p/rho
 
-rhoE.boundaryFieldRef() == rho.boundaryField()*(e.boundaryField() + 0.5*magSqr(U.boundaryField()));
+rhoE.boundaryFieldRef() == rho.boundaryField()*(ea.boundaryField() + 0.5*magSqr(U.boundaryField()));
 
 if (!inviscid)
 {
     fvScalarMatrix eEqn
     (
-           fvm::ddt(rho, e) - fvc::ddt(rho, e)
+           fvm::ddt(rho, ea) - fvc::ddt(rho, ea)
         // alpha in deepflame is considered to calculate h by default (kappa/Cp), so multiply gamma to correct alpha
-         - fvm::laplacian(turbulence->alphaEff()*thermo.gamma(), e)
+         - fvm::laplacian(turbulence->alphaEff()*thermo.gamma(), ea)
          + diffAlphaD
       ==
            fvc::div(hDiffCorrFlux)
     );
 
-    eEqn.solve("e");
+    eEqn.solve("ea");
 
-    ha = e + p/rho;
+    ha = ea + p/rho;
     chemistry->correctThermo();
-    rhoE = rho*(e + 0.5*magSqr(U));
+    rhoE = rho*(ea + 0.5*magSqr(U));
 }
 
 Info<< "min/max(T) = "