multi-gpu base

Allwmake

@@ -19,8 +19,8 @@ wmake src/dynamicMesh
 wmake src/dynamicFvMesh
 
 wmake applications/solvers/df0DFoam
-wmake applications/solvers/dfLowMachFoam
-wmake applications/solvers/dfHighSpeedFoam
-wmake applications/solvers/dfSprayFoam
+# wmake applications/solvers/dfLowMachFoam
+# wmake applications/solvers/dfHighSpeedFoam
+# wmake applications/solvers/dfSprayFoam
 
 wmake applications/utilities/flameSpeed

applications/solvers/df0DFoam/Make/options

@@ -1,4 +1,12 @@
+-include $(GENERAL_RULES)/mplibType
+PYTHON_INC_DIR := $(shell python3 -m pybind11 --includes)
+PYTHON_LIB_DIR := $(shell python3 -c "from distutils import sysconfig; \
+    import os.path as op; v = sysconfig.get_config_vars(); \
+    fpaths = [op.join(v[pv], v['LDLIBRARY']) for pv in ('LIBDIR', 'LIBPL')]; \
+    print(list(filter(op.exists, fpaths))[0])" | xargs dirname)
+
 EXE_INC = -std=c++14 \
+    $(PFLAGS) $(PINC) \
     -I$(LIB_SRC)/transportModels/compressible/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
@@ -10,9 +18,9 @@ EXE_INC = -std=c++14 \
     -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
     -I$(DF_SRC)/CanteraMixture/lnInclude \
     -I$(DF_SRC)/dfChemistryModel/lnInclude \
+    -I$(LIB_SRC)/Pstream/mpi \
     -I$(CANTERA_ROOT)/include \
-    -I$(TORCH_ROOT)/include \
-    -I$(TORCH_ROOT)/include/torch/csrc/api/include
+    $(PYTHON_INC_DIR)
 
 EXE_LIBS = \
     -lcompressibleTransportModels \
@@ -23,12 +31,13 @@ EXE_LIBS = \
     -lmeshTools \
     -lsampling \
     -L$(FOAM_USER_LIBBIN) \
+    -L$(PYTHON_LIB_DIR) \
+    -L$(DF_SRC)/dfChemistryModel \
     -ldfFluidThermophysicalModels \
     -ldfCompressibleTurbulenceModels \
     -lCanteraMixture \
     -ldfChemistryModel \
     $(CANTERA_ROOT)/lib/libcantera.so \
-    $(TORCH_ROOT)/lib/libtorch.so \
-    $(TORCH_ROOT)/lib/libc10.so \
     -rdynamic \
-    -lpthread
+    -lpthread \
+    -lpython3.8

applications/solvers/df0DFoam/df0DFoam.C

@@ -7,27 +7,37 @@
 -------------------------------------------------------------------------------
 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 "dfChemistryModel.H"
 #include "CanteraMixture.H"
 #include "hePsiThermo.H"
+#include <pybind11/embed.h>
+#include <pybind11/numpy.h>
+#include <pybind11/stl.h> //used to convert
 
 #include "fvCFD.H"
 #include "dynamicFvMesh.H"
@@ -44,9 +54,15 @@ Description
 
 int main(int argc, char *argv[])
 {
+    // initialize pybind
+    pybind11::scoped_interpreter guard{}; //start python interpreter
+
     #include "postProcess.H"
 
-    #include "setRootCaseLists.H"
+    // #include "setRootCaseLists.H"
+    #include "listOptions.H"
+    #include "setRootCase2.H"
+    #include "listOutput.H"
     #include "createTime.H"
     #include "createDynamicFvMesh.H"
     #include "createDyMControls.H"
@@ -99,4 +115,4 @@ int main(int argc, char *argv[])
 }
 
 
-// ************************************************************************* //
+// ************************************************************************* //

applications/solvers/df0DFoam/setRootCase2.H

@@ -0,0 +1,5 @@
+Foam::argList args(argc,argv,true,true,/*initialise=*/false);
+if (!args.checkRootCase())
+{
+    Foam::FatalError.exit();
+}

applications/solvers/dfSprayFoam/EEqn.H

@@ -1,14 +1,6 @@
 {
     volScalarField& he = thermo.he();
 
-
-    tmp<fvScalarMatrix> thcSource(new fvScalarMatrix(he, dimEnergy/dimTime));
-    fvScalarMatrix& hcSource = thcSource.ref();
-    forAll(Y, i)
-    {
-        hcSource.source() -= parcels.rhoTrans(i)*chemistry->mixture().Hc(i)/runTime.deltaT();
-    }
-    //hSource.source() -= parcels.hsTrans()/runTime.deltaT();// equivalent to parcels.Sh(he)
     fvScalarMatrix EEqn
     (
         fvm::ddt(rho, he) + mvConvection->fvmDiv(phi, he)
@@ -19,7 +11,6 @@
      ==
         rho*(U&g)
       + parcels.Sh(he)
-      + hcSource
       + fvc::div(hDiffCorrFlux)
       //+ radiation->Sh(thermo, he)
       //+ fvOptions(rho, he)

bashrc.raw

@@ -1,9 +1,7 @@
 export DF_ROOT=pwd
 export DF_SRC=pwd/src
 export SRC_ORIG=pwd/src_orig
-export TORCH_ROOT=TORCH_DIR
 export CANTERA_ROOT=CONDA_PREFIX
-export TORCH_ROOT=TORCH_DIR
 export CANTERA_DATA=$CANTERA_ROOT/share/cantera/data
 export LD_LIBRARY_PATH=$CANTERA_ROOT/lib:$LD_LIBRARY_PATH
-export LD_LIBRARY_PATH=$TORCH_ROOT/lib:$LD_LIBRARY_PATH
+export LD_LIBRARY_PATH=$TORCH_ROOT/lib:$LD_LIBRARY_PATH

examples/df0DFoam/zeroD_cubicReactor/CH4/torchSolver/constant/CanteraTorchProperties

@@ -31,7 +31,7 @@ zeroDReactor
 
 torch on;
 torchModel "new_modle.pt"; 
-torchParameters  
+torchParameters1  
 {
     Xmu  
     (
@@ -140,6 +140,18 @@ torchParameters
         3.357557946945996e-05
         0.0
     );
+    Tact 0  ;
+    Qdotact  0 ;
+}
+torchParameters2
+{
+    Tact 1700;
+    Qdotact  1e5;
+}
+torchParameters3
+{
+    Tact 2500;
+    Qdotact  5e5;
 }
 loadbalancing
 {

examples/df0DFoam/zeroD_cubicReactor/CH4/torchSolver/inference2.py

@@ -0,0 +1,173 @@
+from builtins import Exception, print
+from calendar import prcal
+import torch
+import numpy as np
+import math
+import time
+import json
+import os
+from easydict import EasyDict as edict
+import torch.profiler
+import os
+
+torch.set_printoptions(precision=10)
+print('position 0 in inference.py')
+device = torch.device("cuda")
+device_ids = range(torch.cuda.device_count())
+
+
+
+class MyGELU(torch.nn.Module):
+    def __init__(self):
+        super(MyGELU, self).__init__()
+        self.torch_PI = 3.1415926536
+
+    def forward(self, x):
+        return 0.5 * x * (1 + torch.tanh(
+            math.sqrt(2 / self.torch_PI) * (x + 0.044715 * torch.pow(x, 3))))
+
+
+def json2Parser(json_path):
+    """load json and return parser-like object"""
+    with open(json_path, 'r') as f:
+        args = json.load(f)
+    return edict(args)
+
+
+class Net(torch.nn.Module):
+    def __init__(self):
+        super(Net, self).__init__()
+        neurons = layers
+        self.depth = len(neurons) - 1
+        self.actfun = MyGELU()
+        self.layers = []
+        for i in range(self.depth - 1):
+            self.layers.append(torch.nn.Linear(neurons[i], neurons[i + 1]))
+            self.layers.append(self.actfun)
+        self.layers.append(torch.nn.Linear(neurons[-2], neurons[-1]))  # last layer
+        self.fc = torch.nn.Sequential(*self.layers)
+
+    def forward(self, x):
+        x = self.fc(x)
+        return x
+try:
+    #glbal variable will only init once when called interperter
+    #load parameters from json
+    setting0 = json2Parser('settingsdrm19_0.json')
+    setting1 = json2Parser('settingsdrm19_1.json')
+    setting2 = json2Parser('settingsdrm19_2.json')
+
+    lamda = setting0.power_transform
+    delta_t = setting0.delta_t
+    dim = setting0.dim
+    layers = setting0.layers
+
+    Xmu0 = torch.tensor(setting0.Xmu).unsqueeze(0).to(device)
+    Xstd0 = torch.tensor(setting0.Xstd).unsqueeze(0).to(device=device)
+    Ymu0 = torch.tensor(setting0.Ymu).unsqueeze(0).to(device=device)
+    Ystd0 = torch.tensor(setting0.Ystd).unsqueeze(0).to(device=device)
+
+    Xmu1 = torch.tensor(setting1.Xmu).unsqueeze(0).to(device=device)
+    Xstd1 = torch.tensor(setting1.Xstd).unsqueeze(0).to(device=device)
+    Ymu1 = torch.tensor(setting1.Ymu).unsqueeze(0).to(device=device)
+    Ystd1 = torch.tensor(setting1.Ystd).unsqueeze(0).to(device=device)
+
+    Xmu2 = torch.tensor(setting2.Xmu).unsqueeze(0).to(device=device)
+    Xstd2 = torch.tensor(setting2.Xstd).unsqueeze(0).to(device=device)
+    Ymu2 = torch.tensor(setting2.Ymu).unsqueeze(0).to(device=device)
+    Ystd2 = torch.tensor(setting2.Ystd).unsqueeze(0).to(device=device)
+    print('position 1 in inference.py')
+
+    #load module  
+    model0 = Net()
+    model1 = Net()
+    model2 = Net()
+    check_point0 = torch.load('modeldrm19_0.pt')
+    check_point1 = torch.load('modeldrm19_1.pt')
+    check_point2 = torch.load('modeldrm19_2.pt')
+    model0.load_state_dict(check_point0)
+    model1.load_state_dict(check_point1)
+    model2.load_state_dict(check_point2)
+    model0.to(device=device)
+    model1.to(device=device)
+    model2.to(device=device)
+    if len(device_ids) > 1:
+        model0 = torch.nn.DataParallel(model0, device_ids=device_ids)
+        model1 = torch.nn.DataParallel(model1, device_ids=device_ids)
+        model2 = torch.nn.DataParallel(model2, device_ids=device_ids)
+    print('call init')
+except Exception as e:
+    print(e.args)
+
+
+def inference(vec0, vec1, vec2):
+    '''
+    use model to inference
+    '''
+    #args = np.reshape(args, (-1, 9)) #reshape to formed size
+    vec0 = np.reshape(vec0, (-1, 24))
+    vec1 = np.reshape(vec1, (-1, 24))
+    vec2 = np.reshape(vec2, (-1, 24))
+
+    try:
+        with torch.no_grad():
+            input0_ = torch.from_numpy(vec0).double().to(device=device) #cast ndarray to torch tensor
+            input1_ = torch.from_numpy(vec1).double().to(device=device) #cast ndarray to torch tensor
+            input2_ = torch.from_numpy(vec2).double().to(device=device) #cast ndarray to torch tensor
+
+            # pre_processing
+            rho0 = input0_[:, 0].unsqueeze(1)
+            input0_Y = input0_[:, 3:].clone()
+            input0_bct = input0_[:, 1:]
+            input0_bct[:, 2:] = (input0_bct[:, 2:]**(lamda) - 1) / lamda #BCT
+            input0_normalized = (input0_bct - Xmu0) / Xstd0
+            input0_normalized[:, -1] = 0 #set Y_AR to 0
+            input0_normalized = input0_normalized.float()
+
+            rho1 = input1_[:, 0].unsqueeze(1)
+            input1_Y = input1_[:, 3:].clone()
+            input1_bct = input1_[:, 1:]
+            input1_bct[:, 2:] = (input1_bct[:, 2:]**(lamda) - 1) / lamda #BCT
+            input1_normalized = (input1_bct - Xmu1) / Xstd1
+            input1_normalized[:, -1] = 0 #set Y_AR to 0
+            input1_normalized = input1_normalized.float()
+
+
+            rho2 = input2_[:, 0].unsqueeze(1)
+            input2_Y = input2_[:, 3:].clone()
+            input2_bct = input2_[:, 1:]
+            input2_bct[:, 2:] = (input2_bct[:, 2:]**(lamda) - 1) / lamda #BCT
+            input2_normalized = (input2_bct - Xmu2) / Xstd2
+            input2_normalized[:, -1] = 0 #set Y_AR to 0
+            input2_normalized = input2_normalized.float()
+
+            #inference
+            output0_normalized = model0(input0_normalized)
+            output1_normalized = model1(input1_normalized)
+            output2_normalized = model2(input2_normalized)
+
+
+            # post_processing
+            output0_bct = (output0_normalized * Ystd0 + Ymu0) * delta_t + input0_bct
+            output0_Y = (lamda * output0_bct[:, 2:] + 1)**(1 / lamda)
+            output0_Y = output0_Y / torch.sum(input=output0_Y, dim=1, keepdim=True)
+            output0 = (output0_Y - input0_Y) * rho0 / delta_t   
+            output0 = output0.cpu().numpy()
+
+            output1_bct = (output1_normalized * Ystd1 + Ymu1) * delta_t + input1_bct
+            output1_Y = (lamda * output1_bct[:, 2:] + 1)**(1 / lamda)
+            output1_Y = output1_Y / torch.sum(input=output1_Y, dim=1, keepdim=True)
+            output1 = (output1_Y - input1_Y) * rho1 / delta_t
+            output1 = output1.cpu().numpy()
+
+            output2_bct = (output2_normalized * Ystd2 + Ymu2) * delta_t + input2_bct
+            output2_Y = (lamda * output2_bct[:, 2:] + 1)**(1 / lamda)
+            output2_Y = output2_Y / torch.sum(input=output2_Y, dim=1, keepdim=True)
+            output2 = (output2_Y - input2_Y) * rho2 / delta_t
+            output2 = output2.cpu().numpy()
+
+            result = np.append(output0, output1, axis=0)
+            result = np.append(result, output2, axis=0)
+            return result
+    except Exception as e:
+        print(e.args)

examples/df0DFoam/zeroD_cubicReactor/CH4/torchSolver/system/blockMeshDict

@@ -31,7 +31,7 @@ vertices
 
 blocks
 (
-    hex (0 1 2 3 4 5 6 7) (1 1 1) simpleGrading (1 1 1)
+    hex (0 1 2 3 4 5 6 7) (50 50 50) simpleGrading (1 1 1)
 );
 
 edges

examples/df0DFoam/zeroD_cubicReactor/CH4/torchSolver/system/controlDict

@@ -23,7 +23,7 @@ startTime       0;
 
 stopAt          endTime;
 
-endTime         0.001;
+endTime         1e-3;
 
 deltaT          1e-06;
 
@@ -33,7 +33,7 @@ adjustTimeStep  off;
 
 writeControl    runTime;
 
-writeInterval   0.01
+writeInterval   0.01;
 
 purgeWrite      0;