Gpu

.gitignore

@@ -70,3 +70,6 @@ __pycache__/
 lib/
 bin/
 .vscode/
+result/
+*result*
+*profile*

Allwclean

@@ -25,3 +25,4 @@ wclean ./applications/solvers/dfHighSpeedFoam
 rm -rf src_orig/
 rm -rf bin/
 rm -rf lib/
+rm -rf src_gpu/build

applications/solvers/dfLowMachFoam/CMakeLists.txt

@@ -1,6 +1,8 @@
 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})
@@ -26,6 +28,8 @@ 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
@@ -83,6 +87,9 @@ include_directories(
     ${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
@@ -98,6 +105,9 @@ target_link_libraries(${PROJECT_NAME}
     ${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})

applications/solvers/dfLowMachFoam/EEqn.H

@@ -1,8 +1,113 @@
 {
     volScalarField& he = thermo.he();
+#ifdef GPUSolver_
+    start1 = std::clock();
+    UEqn_GPU.updatePsi(&U[0][0]);
+    UEqn_GPU.correctBoundaryConditions();
+    U.correctBoundaryConditions();
+    K = 0.5*magSqr(U);
+    end1 = std::clock();
+    time_monitor_UEqn += double(end1 - start1) / double(CLOCKS_PER_SEC);
+    time_monitor_UEqn_correctBC += double(end1 - start1) / double(CLOCKS_PER_SEC);
 
+    // prepare data on CPU
+    start1 = std::clock();
+    start2 = std::clock();
+    // const tmp<volScalarField> alphaEff_tmp(thermo.alpha());
+    // const volScalarField& alphaEff = alphaEff_tmp();
+    double *alphaEff = nullptr; // tmp
+    end2 = std::clock();
+    int eeqn_offset = 0;
+    int patchNum = 0;
+
+    forAll(he.boundaryField(), patchi)
+    {
+        patchNum++;
+        const fvsPatchScalarField& pw = mesh.surfaceInterpolation::weights().boundaryField()[patchi];
+        int patchSize = pw.size();
+
+        // construct gradient manually
+        const fvPatchScalarField& hew = he.boundaryField()[patchi];
+        const basicThermo& bThermo = basicThermo::lookupThermo(hew);
+        const scalarField& ppw = bThermo.p().boundaryField()[patchi];
+        fvPatchScalarField& Tw =
+            const_cast<fvPatchScalarField&>(bThermo.T().boundaryField()[patchi]);
+        scalarField& Tw_v = Tw;
+
+        Tw.evaluate();
+        const scalarField& patchDeltaCoeff = mesh.boundary()[patchi].deltaCoeffs();
+        const scalarField heInternal = bThermo.he(ppw, Tw, patchi)();
+        const scalarField heBoundary = bThermo.he(ppw, Tw, mesh.boundary()[patchi].faceCells())();
+        const scalarField patchGradMau = patchDeltaCoeff * (heInternal - heBoundary);
+
+        const scalarField& patchK = K.boundaryField()[patchi];
+        // const scalarField& patchAlphaEff = alphaEff.boundaryField()[patchi]; // not H2Dcopy when use UnityLewis
+        // const scalarField& patchGrad = he.boundaryField()[patchi].gradientBoundaryCoeffs(); // gradient_
+
+        // const DimensionedField<scalar, volMesh>& patchHa_ = he.boundaryField()[patchi];
+        // const gradientEnergyFvPatchScalarField patchHa(mesh.boundary()[patchi], patchHa_);
+        // const scalarField& patchGrad = patchHa.gradient(); // gradient_
+        memcpy(boundary_K + eeqn_offset, &patchK[0], patchSize*sizeof(double));
+        // memcpy(boundary_alphaEff + eeqn_offset, &patchAlphaEff[0], patchSize*sizeof(double));
+        memcpy(boundary_gradient + eeqn_offset, &patchGradMau[0], patchSize*sizeof(double));
+
+        eeqn_offset += patchSize;
+    }
+    end1 = std::clock();
+    time_monitor_EEqn += double(end1 - start1) / double(CLOCKS_PER_SEC);
+    time_monitor_EEqn_mtxAssembly += double(end1 - start1) / double(CLOCKS_PER_SEC);
+    time_monitor_EEqn_mtxAssembly_CPU_prepare += double(end1 - start1) / double(CLOCKS_PER_SEC);
+    fprintf(stderr, "time_monitor_EEqn_mtxAssembly_CPU_prepare: %lf, build alphaEff time: %lf, patchNum: %d\n",
+            time_monitor_EEqn_mtxAssembly_CPU_prepare,
+            double(end2 - start2) / double(CLOCKS_PER_SEC), patchNum);
+
+    // prepare data on GPU
+    start1 = std::clock();
+    he.oldTime();
+    K.oldTime();
+    EEqn_GPU.prepare_data(&he.oldTime()[0], &K[0], &K.oldTime()[0], alphaEff,
+            &dpdt[0], boundary_K, boundary_alphaEff, boundary_gradient);
+    EEqn_GPU.sync();
+    end1 = std::clock();
+    time_monitor_EEqn += double(end1 - start1) / double(CLOCKS_PER_SEC);
+    time_monitor_EEqn_mtxAssembly += double(end1 - start1) / double(CLOCKS_PER_SEC);
+    time_monitor_EEqn_mtxAssembly_GPU_prepare += double(end1 - start1) / double(CLOCKS_PER_SEC);
+
+    start1 = std::clock();
+    EEqn_GPU.initializeTimeStep();
+    EEqn_GPU.fvm_ddt();
+    EEqn_GPU.fvm_div();
+    EEqn_GPU.fvm_laplacian();
+    EEqn_GPU.fvc_ddt();
+    EEqn_GPU.fvc_div_phi_scalar();
+    EEqn_GPU.fvc_div_vector();
+    EEqn_GPU.add_to_source();
+    EEqn_GPU.sync();
+    end1 = std::clock();
+    time_monitor_EEqn += double(end1 - start1) / double(CLOCKS_PER_SEC);
+    time_monitor_EEqn_mtxAssembly += double(end1 - start1) / double(CLOCKS_PER_SEC);
+    time_monitor_EEqn_mtxAssembly_GPU_run += double(end1 - start1) / double(CLOCKS_PER_SEC);
+
+    // check value of mtxAssembly, no time monitor
+    // EEqn_GPU.checkValue(true);
+
+    start1 = std::clock();
+    EEqn_GPU.solve();
+    end1 = std::clock();
+    time_monitor_EEqn += double(end1 - start1) / double(CLOCKS_PER_SEC);
+    time_monitor_EEqn_solve += double(end1 - start1) / double(CLOCKS_PER_SEC);
+
+    start1 = std::clock();
+    EEqn_GPU.updatePsi(&he[0]);
+    he.correctBoundaryConditions();
+    he.write();
+    end1 = std::clock();
+    time_monitor_EEqn += double(end1 - start1) / double(CLOCKS_PER_SEC);
+    time_monitor_EEqn_correctBC += double(end1 - start1) / double(CLOCKS_PER_SEC);
+#else
+    start1 = std::clock();
     fvScalarMatrix EEqn
-        (
+    (
 
             fvm::ddt(rho, he) + mvConvection->fvmDiv(phi, he)
         +   fvc::ddt(rho, K) + fvc::div(phi, K)
@@ -22,8 +127,15 @@
                 )
             )
         );
+    end1 = std::clock();
+    time_monitor_EEqn += double(end1 - start1) / double(CLOCKS_PER_SEC);
+    time_monitor_EEqn_mtxAssembly += double(end1 - start1) / double(CLOCKS_PER_SEC);
 
-        EEqn.relax();
-
-        EEqn.solve("ha");
+    EEqn.relax();
+    start1 = std::clock();
+    EEqn.solve("ha");
+    end1 = std::clock();
+    time_monitor_EEqn += double(end1 - start1) / double(CLOCKS_PER_SEC);
+    time_monitor_EEqn_solve += double(end1 - start1) / double(CLOCKS_PER_SEC);
+#endif
 }

applications/solvers/dfLowMachFoam/Make/options

@@ -1,12 +1,15 @@
 -include $(GENERAL_RULES)/mplibType
 
 EXE_INC = -std=c++14 \
+    -g \
+    -fopenmp \
     -Wno-unused-variable \
     -Wno-unused-but-set-variable \
     -Wno-old-style-cast \
     $(PFLAGS) $(PINC) \
     $(if $(LIBTORCH_ROOT),-DUSE_LIBTORCH,) \
     $(if $(PYTHON_INC_DIR),-DUSE_PYTORCH,) \
+    $(if $(AMGX_DIR),-DGPUSolver_,) \
     -I$(LIB_SRC)/transportModels/compressible/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
@@ -23,7 +26,10 @@ EXE_INC = -std=c++14 \
     -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)
+    $(PYTHON_INC_DIR) \
+    $(if $(AMGX_DIR), -I$(DF_ROOT)/src_gpu,) \
+    $(if $(AMGX_DIR), -I/usr/local/cuda-11.6/include,) \
+    $(if $(AMGX_DIR), -I$(AMGX_DIR)/include,)
 
 EXE_LIBS = \
     -lcompressibleTransportModels \
@@ -44,4 +50,8 @@ EXE_LIBS = \
     $(if $(LIBTORCH_ROOT),-lpthread,) \
     $(if $(LIBTORCH_ROOT),$(DF_SRC)/dfChemistryModel/DNNInferencer/build/libDNNInferencer.so,) \
     $(if $(PYTHON_LIB_DIR),-L$(PYTHON_LIB_DIR),) \
-    $(if $(PYTHON_LIB_DIR),-lpython3.8,)
+    $(if $(PYTHON_LIB_DIR),-lpython3.8,) \
+    $(if $(AMGX_DIR), /usr/local/cuda-11.6/lib64/libcudart.so,) \
+    $(if $(AMGX_DIR), $(DF_ROOT)/src_gpu/build/libdfMatrix.so,) \
+    $(if $(AMGX_DIR), $(AMGX_DIR)/build/libamgxsh.so,)
+

applications/solvers/dfLowMachFoam/UEqn.H

@@ -1,17 +1,132 @@
 // Solve the Momentum equation
+#ifdef GPUSolver_
+    start1 = std::clock();
+    int offset = 0;
+    const tmp<volScalarField> nuEff_tmp(turbulence->nuEff());
+    const volScalarField& nuEff = nuEff_tmp();
+    forAll(U.boundaryField(), patchi)
+    {
+        const scalarField& patchP = p.boundaryField()[patchi];
+        const vectorField& patchU = U.boundaryField()[patchi];
+        const scalarField& patchRho = rho.boundaryField()[patchi];
+        const scalarField& patchNuEff = nuEff.boundaryField()[patchi];
 
-tmp<fvVectorMatrix> tUEqn
-(
-    fvm::ddt(rho, U) + fvm::div(phi, U)
-  + turbulence->divDevRhoReff(U)
-);
-fvVectorMatrix& UEqn = tUEqn.ref();
+        int patchSize = patchP.size();
 
-UEqn.relax();
+        // boundary pressure
+        memcpy(boundary_pressure_init+offset, &patchP[0], patchSize*sizeof(double));
+        // boundary velocity
+        memcpy(boundary_velocity_init+3*offset, &patchU[0][0], 3*patchSize*sizeof(double));
+        // boundary nuEff
+        memcpy(boundary_nuEff_init+offset, &patchNuEff[0], patchSize*sizeof(double));
+        // boundary rho
+        memcpy(boundary_rho_init+offset, &patchRho[0], patchSize*sizeof(double));
+        offset += patchSize;
+    }
+    end1 = std::clock();
+    time_monitor_UEqn += double(end1 - start1) / double(CLOCKS_PER_SEC);
+    time_monitor_UEqn_mtxAssembly += double(end1 - start1) / double(CLOCKS_PER_SEC);
+    time_monitor_UEqn_mtxAssembly_CPU_prepare += double(end1 - start1) / double(CLOCKS_PER_SEC);
 
-if (pimple.momentumPredictor())
-{
-    solve(UEqn == -fvc::grad(p));
+    start1 = std::clock();
+    UEqn_GPU.initializeTimeStep();
+    U.oldTime();
+    UEqn_GPU.fvm_ddt(&U.oldTime()[0][0]);
+    UEqn_GPU.fvm_div(boundary_pressure_init, boundary_velocity_init, boundary_nuEff_init, boundary_rho_init);
+    UEqn_GPU.fvc_grad(&p[0]);
+    UEqn_GPU.fvc_grad_vector();
+    UEqn_GPU.dev2T();
+    UEqn_GPU.fvc_div_tensor(&nuEff[0]);
+    UEqn_GPU.fvm_laplacian();
+    UEqn_GPU.sync();
+    end1 = std::clock();
+    time_monitor_UEqn += double(end1 - start1) / double(CLOCKS_PER_SEC);
+    time_monitor_UEqn_mtxAssembly += double(end1 - start1) / double(CLOCKS_PER_SEC);
+    time_monitor_UEqn_mtxAssembly_GPU_run += double(end1 - start1) / double(CLOCKS_PER_SEC);
+
+    // start2 = std::clock();
+    // fvVectorMatrix turb_source
+    // (
+    //     turbulence->divDevRhoReff(U)
+    // );
+    // end2 = std::clock();
+    // time_monitor_CPU += double(end2 - start2) / double(CLOCKS_PER_SEC);
+
+    // UEqn_GPU.add_fvMatrix(&turb_source.lower()[0], &turb_source.diag()[0], &turb_source.upper()[0], &turb_source.source()[0][0]);
+    // end1 = std::clock();
+    // time_monitor_UEqn += double(end1 - start1) / double(CLOCKS_PER_SEC);
+    // time_monitor_UEqn_mtxAssembly += double(end1 - start1) / double(CLOCKS_PER_SEC);
+
+    // check value
+    // U.oldTime();
+    // tmp<fvVectorMatrix> tUEqn
+    // (
+        // fvm::ddt(rho, U) 
+        // + 
+        // fvm::div(phi, U)
+        // + 
+        // turbulence->divDevRhoReff(U) 
+        // == -fvc::grad(p)
+    // );
+    // fvVectorMatrix& UEqn = tUEqn.ref();
+    // printf("b_cpu = %e\n", UEqn.source()[1][1]);
+    // forAll(U.boundaryField(), patchi){
+        // labelUList sub_boundary = mesh.boundary()[patchi].faceCells();
+        // forAll(sub_boundary, i){
+        //     if (sub_boundary[i] == 1){
+        //         printf("b_cpu_bou = %e\n", UEqn.boundaryCoeffs()[patchi][i][1]);
+        //         printf("patchi = %d, i = %d\n", patchi, i);
+        //     }
+        // }
+    // }
+    // if (pimple.momentumPredictor())
+    // {
+    //     solve(UEqn);
+    //     Info << "U_CPU\n" << U << endl;
+    //     K = 0.5*magSqr(U);
+    // }
+    // UEqn_GPU.checkValue(true);
+#else
+    start1 = std::clock();
+    tmp<fvVectorMatrix> tUEqn
+    (
+        fvm::ddt(rho, U) + fvm::div(phi, U)
+    + turbulence->divDevRhoReff(U) 
+    == -fvc::grad(p)
+    );
+    fvVectorMatrix& UEqn = tUEqn.ref();
+
+    end1 = std::clock();
+    time_monitor_UEqn += double(end1 - start1) / double(CLOCKS_PER_SEC);
+    time_monitor_UEqn_mtxAssembly += double(end1 - start1) / double(CLOCKS_PER_SEC);
+
+    UEqn.relax();
+    start1 = std::clock();
+    if (pimple.momentumPredictor())
+    {
+        solve(UEqn);
+
+        K = 0.5*magSqr(U);
+    }
+    end1 = std::clock();
+    time_monitor_UEqn += double(end1 - start1) / double(CLOCKS_PER_SEC);
+    time_monitor_UEqn_solve += double(end1 - start1) / double(CLOCKS_PER_SEC);
+#endif
+
+// start1 = std::clock();
+// // // std::thread t(&dfMatrix::solve, &UEqn_GPU);
+// UEqn_GPU.solve();
+// end1 = std::clock();
+// time_monitor_UEqn += double(end1 - start1) / double(CLOCKS_PER_SEC);
+// time_monitor_UEqn_solve += double(end1 - start1) / double(CLOCKS_PER_SEC);
+
+// start1 = std::clock();
+// // // t.join();
+// // UEqn_GPU.updatePsi(&U[0][0]);
+// K = 0.5*magSqr(U);
+// end1 = std::clock();
+// time_monitor_UEqn += double(end1 - start1) / double(CLOCKS_PER_SEC);
+// time_monitor_UEqn_mtxAssembly += double(end1 - start1) / double(CLOCKS_PER_SEC);
+// time_monitor_CPU += double(end1 - start1) / double(CLOCKS_PER_SEC);
+// // Info << "U_amgx = " << U << endl;
 
-    K = 0.5*magSqr(U);
-}