diff --git a/Allwmake b/Allwmake
index 574729e3..dd7ff807 100755
--- a/Allwmake
+++ b/Allwmake
@@ -25,5 +25,6 @@ wmake applications/solvers/dfLowMachFoam
 wmake applications/solvers/dfHighSpeedFoam
 wmake applications/solvers/dfSprayFoam
 wmake applications/solvers/dfBuoyancyFoam
+wmake applications/solvers/dfSteadyFoam
 
 wmake applications/utilities/flameSpeed
diff --git a/applications/solvers/dfSteadyFoam/CMakeLists.txt b/applications/solvers/dfSteadyFoam/CMakeLists.txt
new file mode 100644
index 00000000..68168e71
--- /dev/null
+++ b/applications/solvers/dfSteadyFoam/CMakeLists.txt
@@ -0,0 +1,122 @@
+cmake_minimum_required(VERSION 3.5)
+project(dfLowMachFoam LANGUAGES CXX)
+FIND_PACKAGE(MPI REQUIRED)
+FIND_PACKAGE(OpenMP REQUIRED)
+FIND_PACKAGE(CUDA REQUIRED)
+
+# Check valid thirdParty
+if(DEFINED ENV{WM_PROJECT_DIR})
+	MESSAGE(STATUS "OpenFOAM: " $ENV{WM_PROJECT_DIR})
+else()
+	message(FATAL_ERROR "OpenFOAM is not sourced")
+endif(DEFINED ENV{WM_PROJECT_DIR})
+
+if(DEFINED ENV{CANTERA_ROOT})
+	MESSAGE(STATUS "libcantera: " $ENV{CANTERA_ROOT})
+  SET(CANTERA_ROOT $ENV{CANTERA_ROOT})
+else()
+	message(FATAL_ERROR "libcantera directory is not specified")
+endif(DEFINED ENV{CANTERA_ROOT})
+
+# define variables
+SET(OpenFOAM_LIB_DIR $ENV{FOAM_LIBBIN})
+SET(OpenFOAM_SRC $ENV{FOAM_SRC})
+
+SET(DF_ROOT $ENV{DF_ROOT})
+SET(DF_SRC $ENV{DF_SRC})
+SET(SRC_ORIG $ENV{SRC_ORIG})
+
+# set compilation options
+SET(CMAKE_EXE_LINKER_FLAGS "-fuse-ld=bfd -Xlinker --add-needed -Xlinker --no-as-needed")
+SET (CMAKE_C_FLAGS ${CMAKE_C_FLAGS} ${OpenMP_C_FLAGS})
+SET (CMAKE_CXX_FLAGS ${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS})
+
+SET(CMAKE_C_COMPILER g++)
+SET(PATH_LIB_OPENMPI "openmpi-system")  # Foundation version
+SET(EXE_COMPILE_OPTION "-std=c++14 -m64 -Dlinux64 -DWM_ARCH_OPTION=64
+-DWM_DP -DWM_LABEL_SIZE=32 -Wall -Wextra -Wold-style-cast -Wnon-virtual-dtor
+-Wno-unused-parameter -Wno-invalid-offsetof -Wno-attributes -O3
+-DNoRepository -ftemplate-depth-100
+-Wno-unused-variable -Wno-unused-but-set-variable -Wno-old-style-cast -DOMPI_SKIP_MPICXX
+-pthread -fPIC")
+add_definitions("${EXE_COMPILE_OPTION}")
+
+# add header files
+FUNCTION(R_SEARCH search_path return_list)
+    FILE(GLOB_RECURSE new_list ${search_path}/*.H)
+    SET(dir_list "")
+    FOREACH(file_path ${new_list})
+        GET_FILENAME_COMPONENT(dir_path ${file_path} PATH)
+        SET(dir_list ${dir_list} ${dir_path})
+    ENDFOREACH()
+    LIST(REMOVE_DUPLICATES dir_list)
+    SET(${return_list} ${dir_list} PARENT_SCOPE)
+ENDFUNCTION(R_SEARCH)
+
+R_SEARCH(${DF_SRC}/dfCombustionModels dfcombustion_inc)
+R_SEARCH(${DF_SRC}/dfCanteraMixture dfcantera_inc)
+R_SEARCH(${DF_SRC}/lagrangian/intermediate dflagrangianinter_inc)
+R_SEARCH(${DF_SRC}/lagrangian/spray dflagrangianspray_inc)
+R_SEARCH(${DF_SRC}/lagrangian/turbulence dflagrangianturb_inc)
+R_SEARCH(${DF_SRC}/dfChemistryModel dfchemistry_inc)
+R_SEARCH(${DF_SRC}/thermophysicalModels/thermophysicalProperties dfthermophysicalprop_inc)
+R_SEARCH(${DF_SRC}/thermophysicalModels/thermophysicalProperties dfthermophysicalprop_inc)
+R_SEARCH(${DF_SRC}/thermophysicalModels/basic dfthermophysicalbasic_inc)
+R_SEARCH(${DF_SRC}/thermophysicalModels/SLGThermo dfthermophysicalslg_inc)
+R_SEARCH(${DF_SRC}/TurbulenceModels dfturbulence_inc)
+R_SEARCH(${DF_SRC}/dynamicMesh dfnewdynamic_inc)
+R_SEARCH(${DF_SRC}/dynamicFvMesh dffvdynamic_inc)
+
+include_directories(
+    ${OpenFOAM_SRC}/finiteVolume/lnInclude
+    ${OpenFOAM_SRC}/OSspecific/POSIX/lnInclude
+    ${OpenFOAM_SRC}/OpenFOAM/lnInclude
+    ${OpenFOAM_SRC}/transportModels/compressible/lnInclude
+    ${OpenFOAM_SRC}/thermophysicalModels/basic/lnInclude
+    ${OpenFOAM_SRC}/TurbulenceModels/turbulenceModels/lnInclude
+    ${OpenFOAM_SRC}/TurbulenceModels/compressible/lnInclude
+    ${OpenFOAM_SRC}/finiteVolume/cfdTools
+    ${OpenFOAM_SRC}/finiteVolume/lnInclude
+    ${OpenFOAM_SRC}/meshTools/lnInclude
+    ${OpenFOAM_SRC}/sampling/lnInclude
+    ${OpenFOAM_SRC}/dynamicFvMesh/lnInclude
+    ${OpenFOAM_SRC}/Pstream/mpi
+    ${dfcantera_inc}
+    ${dfchemistry_inc}
+    ${dfcombustion_inc}
+    ${CANTERA_ROOT}/include
+    ${MPI_INCLUDE_PATH}
+    ${PROJECT_SOURCE_DIR}
+    ${CUDA_INCLUDE_DIRS}
+    /home/runze/AmgX/AMGX/include
+    /home/runze/deepflame-dev/src_gpu
+)
+
+# add execution
+add_executable(${PROJECT_NAME} ${PROJECT_SOURCE_DIR}/dfLowMachFoam.C)
+
+target_link_libraries(${PROJECT_NAME}
+    $ENV{FOAM_LIBBIN}/libfiniteVolume.so libmeshTools.so libcompressibleTransportModels.so
+    libturbulenceModels.so libsampling.so libOpenFOAM.so
+    ${CANTERA_ROOT}/lib/libcantera_shared.so.2
+    ${DF_ROOT}/lib/libdfChemistryModel.so
+    ${DF_ROOT}/lib/libdfCanteraMixture.so
+    ${DF_ROOT}/lib/libdfFluidThermophysicalModels.so
+    ${DF_ROOT}/lib/libdfCombustionModels.so
+    $ENV{FOAM_LIBBIN}/openmpi-system/libPstream.so
+    ${MPI_LIBRARIES}
+    ${CUDA_LIBRARIES}
+    /home/runze/AmgX/AMGX/build/libamgxsh.so
+    /home/runze/deepflame-dev/src_gpu/build/libdfMatrix.so
+)
+
+if(DEFINED ENV{PYTHON_INC_DIR})
+  add_definitions(-DUSE_PYTORCH)
+  # https://pybind11.readthedocs.io/en/stable/advanced/embedding.html
+  find_package(pybind11)
+  target_link_libraries(${PROJECT_NAME} pybind11::embed)
+endif()
+
+# install
+set(CMAKE_INSTALL_PREFIX ${DF_ROOT})
+install(TARGETS ${PROJECT_NAME} DESTINATION bin)
diff --git a/applications/solvers/dfSteadyFoam/EEqn.H b/applications/solvers/dfSteadyFoam/EEqn.H
new file mode 100644
index 00000000..e14cd68f
--- /dev/null
+++ b/applications/solvers/dfSteadyFoam/EEqn.H
@@ -0,0 +1,27 @@
+{
+    volScalarField& he = thermo.he();
+    
+    fvScalarMatrix EEqn
+    (
+
+         fvm::div(phi, he)
+      + (
+            he.name() == "e"
+          ? fvc::div(phi, volScalarField("Ekp", 0.5*magSqr(U) + p/rho))
+          : fvc::div(phi, volScalarField("K", 0.5*magSqr(U)))
+        )
+      - fvm::laplacian(turbulence->alphaEff(), he)
+     ==
+          combustion->Qdot()
+        + fvOptions(rho, he)
+    );
+
+    EEqn.relax();
+
+    fvOptions.constrain(EEqn);
+
+    EEqn.solve("hs");
+    // EEqn.solve("ha");
+
+    fvOptions.correct(he);
+}
diff --git a/applications/solvers/dfSteadyFoam/Make/files b/applications/solvers/dfSteadyFoam/Make/files
new file mode 100644
index 00000000..4cce0ad5
--- /dev/null
+++ b/applications/solvers/dfSteadyFoam/Make/files
@@ -0,0 +1,3 @@
+dfSteadyFoam.C
+
+EXE = $(DF_APPBIN)/dfSteadyFoam
diff --git a/applications/solvers/dfSteadyFoam/Make/options b/applications/solvers/dfSteadyFoam/Make/options
new file mode 100644
index 00000000..28feb963
--- /dev/null
+++ b/applications/solvers/dfSteadyFoam/Make/options
@@ -0,0 +1,82 @@
+-include $(GENERAL_RULES)/mplibType
+
+EXE_INC = -std=c++14 \
+    -g \
+    -fopenmp \
+    -Wno-unused-variable \
+    -Wno-unused-but-set-variable \
+    -Wno-old-style-cast \
+    -I. \
+    $(PFLAGS) $(PINC) \
+    $(if $(LIBTORCH_ROOT),-DUSE_LIBTORCH,) \
+    $(if $(PYTHON_INC_DIR),-DUSE_PYTORCH,) \
+    -I$(FOAM_APP)/solvers/lagrangian/reactingParcelFoam \
+    -I$(FOAM_APP)/solvers/compressible/rhoSimpleFoam \
+    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
+    -I$(FOAM_APP)/solvers/lagrangian/reactingParcelFoam \
+    -I$(LIB_SRC)/finiteVolume/cfdTools \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
+    -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
+    -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
+    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
+    -I$(DF_SRC)/lagrangian/intermediate/lnInclude \
+    -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
+    -I$(DF_SRC)/lagrangian/spray/lnInclude \
+    -I$(LIB_SRC)/lagrangian/spray/lnInclude \
+    -I$(LIB_SRC)/lagrangian/distributionModels/lnInclude \
+    -I$(DF_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \
+    -I$(DF_SRC)/thermophysicalModels/SLGThermo/lnInclude \
+    -I$(LIB_SRC)/Pstream/mpi \
+    -I$(DF_SRC)/dfCanteraMixture/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,) \
+    $(PYTHON_INC_DIR) \
+    $(if $(AMGX_DIR), -I$(DF_ROOT)/src_gpu,) \
+    $(if $(AMGX_DIR), -I/usr/local/cuda/include,) \
+    $(if $(AMGX_DIR), -I$(AMGX_DIR)/include,) \
+    $(if $(ODE_GPU_SOLVER), -I$(OPENCC_PATH)/include,) \
+    $(if $(ODE_GPU_SOLVER), -DODE_GPU_SOLVER,)
+
+EXE_LIBS = \
+    -lcompressibleTransportModels \
+    -lturbulenceModels \
+    -llagrangian \
+    -lregionModels \
+    -ldfSurfaceFilmModels \
+    -lfiniteVolume \
+    -ltopoChangerFvMesh \
+    -lmeshTools \
+    -lsampling \
+    -L$(DF_LIBBIN) \
+    -ldfFluidThermophysicalModels \
+    -ldfCompressibleTurbulenceModels \
+    -ldfThermophysicalProperties \
+    -ldfSLGThermo \
+    -ldfLagrangianIntermediate \
+    -ldfLagrangianTurbulence \
+    -ldfLagrangianSpray \
+    -ldfCanteraMixture \
+    -ldfChemistryModel \
+    -ldfCombustionModels  \
+    $(CANTERA_ROOT)/lib/libcantera.so \
+    $(if $(LIBTORCH_ROOT),$(LIBTORCH_ROOT)/lib/libtorch.so,) \
+    $(if $(LIBTORCH_ROOT),$(LIBTORCH_ROOT)/lib/libc10.so,) \
+    $(if $(LIBTORCH_ROOT),-rdynamic,) \
+    $(if $(LIBTORCH_ROOT),-lpthread,) \
+    $(if $(LIBTORCH_ROOT),$(DF_SRC)/dfChemistryModel/DNNInferencer/build/libDNNInferencer.so,) \
+    $(if $(PYTHON_LIB_DIR),$(PYTHON_LIB_DIR),) \
+    $(if $(AMGX_DIR), /usr/local/cuda/lib64/libcudart.so,) \
+	$(if $(AMGX_DIR), /usr/local/cuda/lib64/libnccl.so,) \
+    $(if $(AMGX_DIR), $(DF_ROOT)/src_gpu/build/libdfMatrix.so,) \
+    $(if $(AMGX_DIR), $(AMGX_DIR)/build/libamgxsh.so,) \
+    $(if $(ODE_GPU_SOLVER), $(ODE_GPU_SOLVER)/lib/libopencc.so,)
diff --git a/applications/solvers/dfSteadyFoam/UEqn.H b/applications/solvers/dfSteadyFoam/UEqn.H
new file mode 100644
index 00000000..57d5dd38
--- /dev/null
+++ b/applications/solvers/dfSteadyFoam/UEqn.H
@@ -0,0 +1,24 @@
+// Solve the Momentum equation
+
+    MRF.correctBoundaryVelocity(U);
+
+    tmp<fvVectorMatrix> tUEqn
+    (
+          fvm::div(phi, U)
+      +   MRF.DDt(rho, U)
+      +   turbulence->divDevRhoReff(U)
+    ==
+        fvOptions(rho, U)
+    );
+    fvVectorMatrix& UEqn = tUEqn.ref();
+
+    UEqn.relax();
+
+    fvOptions.constrain(UEqn);
+
+    solve(UEqn == -fvc::grad(p));
+
+    fvOptions.correct(U);
+
+
+
diff --git a/applications/solvers/dfSteadyFoam/YEqn.H b/applications/solvers/dfSteadyFoam/YEqn.H
new file mode 100644
index 00000000..60e0bae9
--- /dev/null
+++ b/applications/solvers/dfSteadyFoam/YEqn.H
@@ -0,0 +1,59 @@
+     // should only for CPUSolver
+
+    tmp<fv::convectionScheme<scalar>> mvConvection
+    (
+        fv::convectionScheme<scalar>::New
+        (
+            mesh,
+            fields,
+            phi,
+            mesh.divScheme("div(phi,Yi_h)")
+        )
+    );
+
+
+    label flag_mpi_init;
+    MPI_Initialized(&flag_mpi_init);
+    if(flag_mpi_init) MPI_Barrier(PstreamGlobals::MPI_COMM_FOAM);
+
+    {
+               
+    combustion->correct();
+
+    //label flag_mpi_init;
+    //MPI_Initialized(&flag_mpi_init);
+    if(flag_mpi_init) MPI_Barrier(PstreamGlobals::MPI_COMM_FOAM);
+
+
+    volScalarField Yt(0.0*Y[0]);
+    forAll(Y, i)
+    {
+        if (i != inertIndex)
+        {
+            volScalarField& Yi = Y[i];
+            
+            fvScalarMatrix YiEqn
+            (
+                mvConvection->fvmDiv(phi, Yi)
+              - fvm::laplacian(turbulence->muEff(), Yi)
+             ==
+                combustion->R(Yi)
+              + fvOptions(rho, Yi)
+            );
+
+            YiEqn.relax();
+
+            fvOptions.constrain(YiEqn);
+
+            YiEqn.solve("Yi");
+
+            fvOptions.correct(Yi);
+
+            Yi.max(0.0);
+            Yt += Yi;
+        }
+    }
+
+    Y[inertIndex] = scalar(1) - Yt;
+    Y[inertIndex].max(0.0);
+}
diff --git a/applications/solvers/dfSteadyFoam/createFields.H b/applications/solvers/dfSteadyFoam/createFields.H
new file mode 100644
index 00000000..5ec3d786
--- /dev/null
+++ b/applications/solvers/dfSteadyFoam/createFields.H
@@ -0,0 +1,209 @@
+Info<< "Reading thermophysical properties\n" << endl;
+
+CanteraMixture::setEnergyName("hs");
+// 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");
+
+// SLGThermo slgThermo(mesh, thermo);
+
+const volScalarField& psi = thermo.psi();
+volScalarField& p = thermo.p();
+volScalarField& T = thermo.T();
+T.correctBoundaryConditions();
+volScalarField rho
+(
+    IOobject
+    (
+        "rho",
+        runTime.timeName(),
+        mesh,
+        IOobject::READ_IF_PRESENT,
+        IOobject::AUTO_WRITE
+    ),
+    thermo.rho()
+);
+
+
+Info<< "Reading field U\n" << endl;
+volVectorField U
+(
+    IOobject
+    (
+        "U",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);
+
+// #include "compressibleCreatePhi.H"
+surfaceScalarField phi
+(
+    IOobject
+    (
+        "phi",
+        runTime.timeName(),
+        mesh,
+        IOobject::READ_IF_PRESENT,
+        IOobject::AUTO_WRITE
+    ),
+    linearInterpolate(rho*U) & mesh.Sf()
+);
+
+U.correctBoundaryConditions();
+p.correctBoundaryConditions();
+
+pressureControl pressureControl(p, rho, simple.dict(), false);
+
+mesh.setFluxRequired(p.name());
+
+Info<< "Creating turbulence model\n" << endl;
+autoPtr<compressible::turbulenceModel> turbulence
+(
+    compressible::turbulenceModel::New
+    (
+        rho,
+        U,
+        phi,
+        thermo
+    )
+);
+
+
+Info<< "Creating reaction model\n" << endl;
+autoPtr<CombustionModel<basicThermo>> combustion
+(
+    CombustionModel<basicThermo>::New(thermo, turbulence())
+);
+Info<< "end Creating reaction model\n" << endl;
+
+
+const word combModelName(mesh.objectRegistry::lookupObject<IOdictionary>("combustionProperties").lookup("combustionModel"));
+Info << "Combustion Model Name is confirmed as "<< combModelName << endl;
+
+const word turbName(mesh.objectRegistry::lookupObject<IOdictionary>("turbulenceProperties").lookup("simulationType"));
+
+dfChemistryModel<basicThermo>* chemistry = combustion->chemistry();
+PtrList<volScalarField>& Y = chemistry->Y();
+forAll(Y, i)
+{
+    Y[i].correctBoundaryConditions();
+}
+const word inertSpecie(chemistry->lookup("inertSpecie"));
+const label inertIndex(chemistry->species()[inertSpecie]);
+// chemistry->setEnergyName("hs");
+// chemistry->setEnergyName("hs");
+// Info << "energy name is  "<< thermo.he().name() << endl;
+chemistry->updateEnergy();
+
+
+if(combModelName!="ESF" && combModelName!="flareFGM"  && combModelName!="DeePFGM" && combModelName!="FSD")
+{
+    chemistry->correctThermo();
+    Info<< "At initial time, min/max(T) = " << min(T).value() << ", " << max(T).value() << endl;
+}
+rho = thermo.rho();
+
+//for dpdt
+
+Info<< "Creating field kinetic energy K\n" << endl;
+volScalarField K("K", 0.5*magSqr(U));
+
+multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
+
+dimensionedScalar initialMass = fvc::domainIntegrate(rho);
+
+#include "createMRF.H"
+#include "createFvOptions.H"
+
+if(combModelName!="ESF" && combModelName!="flareFGM"  && combModelName!="DeePFGM" && combModelName!="FSD")
+{
+forAll(Y, i)
+{
+    fields.add(Y[i]);
+}
+fields.add(thermo.he());
+}
+
+
+const scalar Sct = chemistry->lookupOrDefault("Sct", 1.);
+// volScalarField diffAlphaD
+// (
+//     IOobject
+//     (
+//         "diffAlphaD",
+//         runTime.timeName(),
+//         mesh,
+//         IOobject::NO_READ,
+//         IOobject::NO_WRITE
+//     ),
+//     mesh,
+//     dimensionedScalar(dimEnergy/dimTime/dimVolume, 0)
+// );
+// volVectorField hDiffCorrFlux
+// (
+//     IOobject
+//     (
+//         "hDiffCorrFlux",
+//         runTime.timeName(),
+//         mesh,
+//         IOobject::NO_READ,
+//         IOobject::NO_WRITE
+//     ),
+//     mesh,
+//     dimensionedVector(dimensionSet(1,0,-3,0,0,0,0), Zero)
+// );
+// volVectorField sumYDiffError
+// (
+//     IOobject
+//     (
+//         "sumYDiffError",
+//         runTime.timeName(),
+//         mesh,
+//         IOobject::NO_READ,
+//         IOobject::NO_WRITE
+//     ),
+//     mesh,
+//     dimensionedVector("sumYDiffError", dimDynamicViscosity/dimLength, Zero)
+// );
+
+IOdictionary CanteraTorchProperties
+(
+    IOobject
+    (
+        "CanteraTorchProperties",
+        runTime.constant(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::NO_WRITE
+    )
+);
+
+volScalarField UEqn_A
+(
+    IOobject
+    (
+        "A("+U.name()+')',
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    mesh,
+    dimensionedScalar(dimensionSet(1,-3,-1,0,0,0,0), Zero),
+    extrapolatedCalculatedFvPatchScalarField::typeName
+);
+
+const Switch splitting = CanteraTorchProperties.lookupOrDefault("splittingStrategy", false);
+#ifdef USE_PYTORCH
+    const Switch log_ = CanteraTorchProperties.subDict("TorchSettings").lookupOrDefault("log", false);
+    const Switch torch_ = CanteraTorchProperties.subDict("TorchSettings").lookupOrDefault("torch", false);
+#endif
+#ifdef USE_LIBTORCH
+    const Switch log_ = CanteraTorchProperties.subDict("TorchSettings").lookupOrDefault("log", false);
+    const Switch torch_ = CanteraTorchProperties.subDict("TorchSettings").lookupOrDefault("torch", false);
+#endif
diff --git a/applications/solvers/dfSteadyFoam/dfSteadyFoam.C b/applications/solvers/dfSteadyFoam/dfSteadyFoam.C
new file mode 100644
index 00000000..ec4b047a
--- /dev/null
+++ b/applications/solvers/dfSteadyFoam/dfSteadyFoam.C
@@ -0,0 +1,206 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2011-2019 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
+    rhoPimpleFoam
+
+Description
+    Transient solver for turbulent flow of compressible fluids for HVAC and
+    similar applications, with optional mesh motion and mesh topology changes.
+
+    Uses the flexible PIMPLE (PISO-SIMPLE) solution for time-resolved and
+    pseudo-transient simulations.
+
+\*---------------------------------------------------------------------------*/
+#include "stdlib.h"
+#include "dfChemistryModel.H"
+#include "CanteraMixture.H"
+// #include "hePsiThermo.H"
+#include "heRhoThermo.H"
+
+#ifdef USE_PYTORCH
+#include <pybind11/embed.h>
+#include <pybind11/numpy.h>
+#include <pybind11/stl.h> //used to convert
+#endif
+
+#ifdef USE_LIBTORCH
+#include <torch/script.h>
+#include "DNNInferencer.H"
+#endif
+
+#ifdef ODE_GPU_SOLVER
+#include "opencc.h"
+#endif
+
+#include "fvCFD.H"
+#include "fluidThermo.H"
+#include "turbulentFluidThermoModel.H"
+
+// #include "psiReactionThermo.H"
+#include "basicThermo.H"
+#include "CombustionModel.H"
+#include "multivariateScheme.H"
+
+#include "simpleControl.H"
+#include "pressureControl.H"
+#include "fvOptions.H"
+#include "fvcSmooth.H"
+
+// #include "PstreamGlobals.H"
+// #include "CombustionModel.H"
+
+// to be decided
+// #include "basicSprayCloud.H"
+// #include "SLGThermo.H"
+
+//#define GPUSolver_
+// #define TIME
+// #define DEBUG_
+// #define SHOW_MEMINFO
+
+    #ifdef GPUSolver_
+    #include "dfMatrixDataBase.H"
+    #include "AmgXSolver.H"
+    #include "dfUEqn.H"
+    #include "dfYEqn.H"
+    #include "dfRhoEqn.H"
+    #include "dfEEqn.H"
+    #include "dfpEqn.H"
+    #include "dfMatrixOpBase.H"
+    #include "dfNcclBase.H"
+    #include "dfThermo.H"
+    #include "dfChemistrySolver.H"
+    #include <cuda_runtime.h>
+    #include <thread>
+
+    #include "processorFvPatchField.H"
+    #include "cyclicFvPatchField.H"
+    #include "processorCyclicFvPatchField.H"
+    #include "createGPUSolver.H"
+
+    #include "upwind.H"
+    #include "CanteraMixture.H"
+    #include "multivariateGaussConvectionScheme.H"
+    #include "limitedSurfaceInterpolationScheme.H"
+#else
+    // #include "processorFvPatchField.H"
+    // #include "cyclicFvPatchField.H"
+    // #include "multivariateGaussConvectionScheme.H"
+    // #include "limitedSurfaceInterpolationScheme.H"
+    int myRank = -1;
+    int mpi_init_flag = 0;
+#endif
+
+int offset;
+
+#ifdef TIME
+    #define TICK_START \
+        start_new = std::clock(); 
+    #define TICK_STOP(prefix) \
+        stop_new = std::clock(); \
+        Foam::Info << #prefix << " time = " << double(stop_new - start_new) / double(CLOCKS_PER_SEC) << " s" << Foam::endl;
+#else
+    #define TICK_START
+    #define TICK_STOP(prefix)
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+#ifdef USE_PYTORCH
+    pybind11::scoped_interpreter guard{};//start python interpreter
+#endif
+    #include "postProcess.H"
+    #include "setRootCaseLists.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "createControl.H"
+    #include "createFields.H"
+    #include "initContinuityErrs.H"
+
+    #ifdef ODE_GPU_SOLVER
+    #include "createFields_GPU.H"
+    #endif
+
+    turbulence->validate();
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+    
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (simple.loop(runTime))
+    {
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        // #include "rhoEqn.H"
+
+        #include "UEqn.H"
+
+        if(combModelName!="ESF" && combModelName!="flareFGM"  && combModelName!="DeePFGM" && combModelName!="FSD")
+        {
+            #include "YEqn.H"
+
+            #include "EEqn.H"
+            
+            chemistry->correctThermo();
+
+            Info<< "min/max(T) = "
+                << min(T).value() << ", " << max(T).value() << endl;
+        }
+        else
+        {
+            combustion->correct();
+        }
+
+        // update T for debug
+
+        if (simple.consistent())
+        {
+            #include "pcEqn.H"
+        }
+        else
+        {
+            #include "pEqn.H"
+        }
+
+        Info<< "min/max(p) = " << min(p).value() << ", " << max(p).value() << endl;
+
+        turbulence->correct();
+        // rho = thermo.rho();
+
+        runTime.write();
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/dfSteadyFoam/pEqn.H b/applications/solvers/dfSteadyFoam/pEqn.H
new file mode 100644
index 00000000..7011a1ad
--- /dev/null
+++ b/applications/solvers/dfSteadyFoam/pEqn.H
@@ -0,0 +1,106 @@
+volScalarField rAU(1.0/UEqn.A());
+surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
+volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
+tUEqn.clear();
+
+bool closedVolume = false;
+
+surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(rho)*fvc::flux(HbyA));
+MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
+
+// Update the pressure BCs to ensure flux consistency
+constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF);
+
+if (simple.transonic())
+{
+    surfaceScalarField phid
+    (
+        "phid",
+        (fvc::interpolate(psi)/fvc::interpolate(rho))*phiHbyA
+    );
+
+    phiHbyA -= fvc::interpolate(psi*p)*phiHbyA/fvc::interpolate(rho);
+
+    while (simple.correctNonOrthogonal())
+    {
+        fvScalarMatrix pEqn
+        (
+            fvc::div(phiHbyA)
+          + fvm::div(phid, p)
+          - fvm::laplacian(rhorAUf, p)
+          ==
+            fvOptions(psi, p, rho.name())
+        );
+
+        // Relax the pressure equation to ensure diagonal-dominance
+        pEqn.relax();
+
+        pEqn.setReference
+        (
+            pressureControl.refCell(),
+            pressureControl.refValue()
+        );
+
+        pEqn.solve();
+
+        if (simple.finalNonOrthogonalIter())
+        {
+            phi = phiHbyA + pEqn.flux();
+        }
+    }
+}
+else
+{
+    closedVolume = adjustPhi(phiHbyA, U, p);
+
+    while (simple.correctNonOrthogonal())
+    {
+        fvScalarMatrix pEqn
+        (
+            fvc::div(phiHbyA)
+          - fvm::laplacian(rhorAUf, p)
+          ==
+            fvOptions(psi, p, rho.name())
+        );
+
+        pEqn.setReference
+        (
+            pressureControl.refCell(),
+            pressureControl.refValue()
+        );
+
+        pEqn.solve();
+
+        if (simple.finalNonOrthogonalIter())
+        {
+            phi = phiHbyA + pEqn.flux();
+        }
+    }
+}
+
+#include "incompressible/continuityErrs.H"
+
+// Explicitly relax pressure for momentum corrector
+p.relax();
+
+U = HbyA - rAU*fvc::grad(p);
+U.correctBoundaryConditions();
+fvOptions.correct(U);
+
+pressureControl.limit(p);
+
+// For closed-volume cases adjust the pressure and density levels
+// to obey overall mass continuity
+if (closedVolume)
+{
+    p += (initialMass - fvc::domainIntegrate(psi*p))
+        /fvc::domainIntegrate(psi);
+    p.correctBoundaryConditions();
+}
+
+rho = thermo.rho();
+
+if (!simple.transonic())
+{
+    rho.relax();
+}
diff --git a/applications/solvers/dfSteadyFoam/pcEqn.H b/applications/solvers/dfSteadyFoam/pcEqn.H
new file mode 100644
index 00000000..3020d711
--- /dev/null
+++ b/applications/solvers/dfSteadyFoam/pcEqn.H
@@ -0,0 +1,119 @@
+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));
+tUEqn.clear();
+
+bool closedVolume = false;
+
+surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(rho)*fvc::flux(HbyA));
+MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
+
+volScalarField rhorAtU("rhorAtU", rho*rAtU);
+
+// Update the pressure BCs to ensure flux consistency
+constrainPressure(p, rho, U, phiHbyA, rhorAtU, MRF);
+
+if (simple.transonic())
+{
+    surfaceScalarField phid
+    (
+        "phid",
+        (fvc::interpolate(psi)/fvc::interpolate(rho))*phiHbyA
+    );
+
+    phiHbyA +=
+        fvc::interpolate(rho*(rAtU - rAU))*fvc::snGrad(p)*mesh.magSf()
+      - fvc::interpolate(psi*p)*phiHbyA/fvc::interpolate(rho);
+
+    HbyA -= (rAU - rAtU)*fvc::grad(p);
+
+    while (simple.correctNonOrthogonal())
+    {
+        fvScalarMatrix pEqn
+        (
+            fvc::div(phiHbyA)
+          + fvm::div(phid, p)
+          - fvm::laplacian(rhorAtU, p)
+         ==
+            fvOptions(psi, p, rho.name())
+        );
+
+        // Relax the pressure equation to maintain diagonal dominance
+        pEqn.relax();
+
+        pEqn.setReference
+        (
+            pressureControl.refCell(),
+            pressureControl.refValue()
+        );
+
+        pEqn.solve();
+
+        if (simple.finalNonOrthogonalIter())
+        {
+            phi = phiHbyA + pEqn.flux();
+        }
+    }
+}
+else
+{
+    closedVolume = adjustPhi(phiHbyA, U, p);
+
+    phiHbyA += fvc::interpolate(rho*(rAtU - rAU))*fvc::snGrad(p)*mesh.magSf();
+    HbyA -= (rAU - rAtU)*fvc::grad(p);
+
+    while (simple.correctNonOrthogonal())
+    {
+        fvScalarMatrix pEqn
+        (
+            fvc::div(phiHbyA)
+          - fvm::laplacian(rhorAtU, p)
+          ==
+            fvOptions(psi, p, rho.name())
+        );
+
+        pEqn.setReference
+        (
+            pressureControl.refCell(),
+            pressureControl.refValue()
+        );
+
+        pEqn.solve();
+
+        if (simple.finalNonOrthogonalIter())
+        {
+            phi = phiHbyA + pEqn.flux();
+        }
+    }
+}
+
+// The incompressibe form of the continuity error check is appropriate for
+// steady-state compressible also.
+#include "incompressible/continuityErrs.H"
+
+// Explicitly relax pressure for momentum corrector
+p.relax();
+
+U = HbyA - rAtU*fvc::grad(p);
+U.correctBoundaryConditions();
+fvOptions.correct(U);
+
+pressureControl.limit(p);
+
+// For closed-volume cases adjust the pressure and density levels
+// to obey overall mass continuity
+if (closedVolume)
+{
+    p += (initialMass - fvc::domainIntegrate(psi*p))
+        /fvc::domainIntegrate(psi);
+    p.correctBoundaryConditions();
+}
+
+rho = thermo.rho();
+
+if (!simple.transonic())
+{
+    rho.relax();
+}
diff --git a/applications/solvers/dfSteadyFoam/rhoEqn.H b/applications/solvers/dfSteadyFoam/rhoEqn.H
new file mode 100644
index 00000000..3987109b
--- /dev/null
+++ b/applications/solvers/dfSteadyFoam/rhoEqn.H
@@ -0,0 +1,45 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+Global
+    rhoEqn
+
+Description
+    Solve the continuity for density.
+
+\*---------------------------------------------------------------------------*/
+
+    fvScalarMatrix rhoEqn
+    (
+        fvm::ddt(rho)
+      + fvc::div(phi)
+    );
+    
+    // fvOptions.constrain(rhoEqn);
+
+    rhoEqn.solve();
+
+    // fvOptions.correct(rho);
+    
+    
+// ************************************************************************* //