first commit

GPUTest/GPUTestRefBase.H

@@ -0,0 +1,63 @@
+
+// unittest of fvc::grad(U)
+void test_fvc_grad_vector_orig(dfMatrixDataBase& dfDataBase, Foam::fvMesh& mesh, volVectorField& U, initType type,
+        dfMatrixDataBaseOrig* dfDataBaseOrig) 
+{
+    if (type == initType::randomInit) {
+        U.oldTime();
+        randomInitField<volVectorField>(U);
+    }
+
+    // run CPU
+    // volTensorField fvc_ouput_tensor = fvc::grad(U);
+    volTensorField fvc_ouput_tensor = gaussGradSchemeGrad(U);
+
+    // prepare for run GPU
+    // prepare U on GPU
+    uploadRegisteredField<volVectorField>(dfDataBase, U, "u");
+
+    double *d_fvc_ouput_tensor = nullptr, *d_fvc_ouput_boundary_tensor = nullptr, *d_fvc_ouput_boundary_tensor_init = nullptr;
+    checkCudaErrors(cudaMalloc((void**)&d_fvc_ouput_tensor, dfDataBase.cell_value_tsr_bytes));
+    checkCudaErrors(cudaMalloc((void**)&d_fvc_ouput_boundary_tensor, dfDataBase.boundary_surface_value_tsr_bytes));
+    checkCudaErrors(cudaMalloc((void**)&d_fvc_ouput_boundary_tensor_init, dfDataBase.boundary_surface_value_tsr_bytes));
+    checkCudaErrors(cudaMemset(d_fvc_ouput_tensor, 0, dfDataBase.cell_value_tsr_bytes));
+    checkCudaErrors(cudaMemset(d_fvc_ouput_boundary_tensor, 0, dfDataBase.boundary_surface_value_tsr_bytes));
+    checkCudaErrors(cudaMemset(d_fvc_ouput_boundary_tensor_init, 0, dfDataBase.boundary_surface_value_tsr_bytes));
+
+    fvc_grad_vector_orig(dfDataBase.stream, dfDataBaseOrig, dfDataBase, d_fvc_ouput_tensor, d_fvc_ouput_boundary_tensor_init, d_fvc_ouput_boundary_tensor);
+
+    // compare result
+    bool printFlag = false;
+    std::vector<double> h_fvc_ouput_tensor(dfDataBase.num_cells * 9);
+    checkCudaErrors(cudaMemcpy(h_fvc_ouput_tensor.data(), d_fvc_ouput_tensor, dfDataBase.cell_value_tsr_bytes, cudaMemcpyDeviceToHost));
+    checkVectorEqual(dfDataBase.num_cells * 9, &fvc_ouput_tensor[0][0], h_fvc_ouput_tensor.data(), 1e-14, printFlag);
+}
+
+void test_fvc_grad_scalar_orig(dfMatrixDataBase& dfDataBase, Foam::fvMesh& mesh, volScalarField& p, initType type,
+        dfMatrixDataBaseOrig* dfDataBaseOrig) 
+{
+    if (type == initType::randomInit) {
+        p.oldTime();
+        randomInitField<volScalarField>(p);
+    }
+
+    // run CPU
+    // volVectorField fvc_ouput_vector = fvc::grad(p);
+    volVectorField fvc_ouput_vector = gaussGradSchemeGrad(p);
+
+    // prepare for run GPU
+    // prepare p on GPU
+    uploadRegisteredField<volScalarField>(dfDataBase, p, "p");
+
+    double *d_fvc_ouput_vector = nullptr;
+    checkCudaErrors(cudaMalloc((void**)&d_fvc_ouput_vector, dfDataBase.cell_value_vec_bytes));
+    checkCudaErrors(cudaMemset(d_fvc_ouput_vector, 0, dfDataBase.cell_value_vec_bytes));
+
+    fvc_grad_scalar_orig(dfDataBase.stream, dfDataBaseOrig, dfDataBase, d_fvc_ouput_vector);
+
+    // compare result
+    bool printFlag = false;
+    std::vector<double> h_fvc_ouput_vector(dfDataBase.num_cells * 3);
+    checkCudaErrors(cudaMemcpy(h_fvc_ouput_vector.data(), d_fvc_ouput_vector, dfDataBase.cell_value_vec_bytes, cudaMemcpyDeviceToHost));
+    checkVectorEqual(dfDataBase.num_cells * 3, &fvc_ouput_vector[0][0], h_fvc_ouput_vector.data(), 1e-14, printFlag);
+}

GPUTest/Make/files

@@ -0,0 +1,4 @@
+unittest.C
+
+EXE = $(DF_APPBIN)/unitTest
+

GPUTest/Make/options

@@ -0,0 +1,50 @@
+-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,) \
+    -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$(LIB_SRC)/finiteVolume/cfdTools \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/dynamicFvMesh/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 \
+	-I$(DF_ROOT)/src_gpu \
+	-I$(DF_ROOT)/src_gpu_orig \
+	-I$(DF_ROOT)/GPUTestRef/lnInclude \
+	-I/usr/local/cuda-11.6/include \
+	-I$(AMGX_DIR)/include
+
+EXE_LIBS = \
+    -lcompressibleTransportModels \
+    -lturbulenceModels \
+    -lfiniteVolume \
+    -lmeshTools \
+    -lsampling \
+    -L$(DF_LIBBIN) \
+    -ldfFluidThermophysicalModels \
+    -ldfCompressibleTurbulenceModels \
+    -ldfCanteraMixture \
+    -ldfChemistryModel \
+    -ldfCombustionModels \
+	-ldfGenMatrix \
+    $(CANTERA_ROOT)/lib/libcantera.so \
+	/usr/local/cuda-11.6/lib64/libcudart.so \
+	$(AMGX_DIR)/build/libamgxsh.so \
+    $(DF_ROOT)/src_gpu/build/libdfMatrix.so \
+	$(DF_ROOT)/src_gpu_orig/build/libdfMatrixOrig.so 
+

GPUTest/correctPhi.H

@@ -0,0 +1,12 @@
+CorrectPhi
+(
+    U,
+    phi,
+    p,
+    rho,
+    psi,
+    dimensionedScalar("rAUf", dimTime, 1),
+    divrhoU(),
+    pimple,
+    true
+);

GPUTest/createFields.H

@@ -0,0 +1,176 @@
+#include "createRDeltaT.H"
+
+Info<< "Reading thermophysical properties\n" << endl;
+
+// 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");
+
+const volScalarField& psi = thermo.psi();
+volScalarField& p = thermo.p();
+volScalarField& T = thermo.T();
+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"
+
+pressureControl pressureControl(p, rho, pimple.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 field dpdt\n" << endl;
+volScalarField dpdt
+(
+    IOobject
+    (
+        "dpdt",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    mesh,
+    dimensionedScalar("dpdt",p.dimensions()/dimTime, 0)
+);
+
+
+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();
+const word inertSpecie(chemistry->lookup("inertSpecie"));
+const label inertIndex(chemistry->species()[inertSpecie]);
+chemistry->setEnergyName("ha");
+chemistry->updateEnergy();
+
+
+chemistry->correctThermo();
+Info<< "At initial time, min/max(T) = " << min(T).value() << ", " << max(T).value() << endl;
+
+//for dpdt
+
+Info<< "Creating field kinetic energy K\n" << endl;
+volScalarField K("K", 0.5*magSqr(U));
+
+multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
+
+if(combModelName!="flareFGM")
+{
+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
+    )
+);
+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