Gpu dev

applications/solvers/dfLowMachFoam/UEqn.H

@@ -11,7 +11,6 @@
         fvm::div(phi, U)
         +  
         turbulence->divDevRhoReff(U)
-        == -fvc::grad(p)
     );
     fvVectorMatrix& UEqn = tUEqn.ref();
     TICK_STOP(CPU assembly time);
@@ -64,7 +63,7 @@
 #if defined DEBUG_
     // UEqn.relax();
     TICK_START;
-    UEqn.solve();
+    solve(UEqn == -fvc::grad(p));
     K.oldTime();
     K = 0.5*magSqr(U);
     TICK_STOP(CPU solve time);
@@ -118,7 +117,7 @@
         //     h_internal_coeffs.data(), h_boundary_coeffs.data(), 
         //     // &gradU[0][0], h_boundary_gradU,
         //     printFlag);
-        //UEqn_GPU.compareU(&U[0][0], h_boundary_u_tmp, printFlag);
+        // UEqn_GPU.compareU(&U[0][0], h_boundary_u_tmp, printFlag);
     }
     DEBUG_TRACE;
 #endif
@@ -129,7 +128,6 @@
     (
         fvm::ddt(rho, U) + fvm::div(phi, U)
     + turbulence->divDevRhoReff(U) 
-    == -fvc::grad(p)
     );
     fvVectorMatrix& UEqn = tUEqn.ref();
     end1 = std::clock();
@@ -140,7 +138,7 @@
     start1 = std::clock();
     if (pimple.momentumPredictor())
     {
-        solve(UEqn);
+        solve(UEqn == -fvc::grad(p));
         K.oldTime();
         K = 0.5*magSqr(U);
     }

applications/solvers/dfLowMachFoam/YEqn.H

@@ -22,6 +22,15 @@
     //randomInitField<volScalarField>(const_cast<volScalarField&>(chemistry->hai(5)));
     //randomInitField<volScalarField>(const_cast<volScalarField&>(chemistry->hai(6)));
 
+    // auto& mgcs = dynamic_cast<fv::multivariateGaussConvectionScheme<scalar>&>(mvConvection.ref());
+    // tmp<surfaceInterpolationScheme<scalar>> tinterpScheme_ = mgcs.interpolationScheme()()(Y[0]);
+    // tmp<surfaceScalarField> tweights = tinterpScheme_().weights(Y[0]);
+    // const surfaceScalarField& weights = tweights();
+    // Info << "CPU weights\n" << weights << endl;
+
+    // auto& limitedScheme_ = dynamic_cast<const limitedSurfaceInterpolationScheme<scalar>&>(tinterpScheme_());
+    // Info << "CPU limiter\n" << limitedScheme_.limiter(Y[0]) << endl;
+
     forAll(Y, i)
     {
         sumYDiffError += chemistry->rhoD(i)*fvc::grad(Y[i]);
@@ -44,6 +53,12 @@
         )
     );
 
+    // auto& mgcs = dynamic_cast<fv::multivariateGaussConvectionScheme<scalar>&>(mvConvection.ref());
+    // tmp<surfaceInterpolationScheme<scalar>> tinterpScheme_ = mgcs.interpolationScheme()()(Y[0]);
+    // tmp<surfaceScalarField> tweights = tinterpScheme_().weights(Y[0]);
+    // const surfaceScalarField& weights = tweights();
+    // Info << "CPU weights\n" << weights << endl;
+
     start1 = std::clock();
     forAll(Y, i)
     {
@@ -62,17 +77,6 @@ MPI_Initialized(&flag_mpi_init);
 if(flag_mpi_init) MPI_Barrier(PstreamGlobals::MPI_COMM_FOAM);
 
 {
-    if (!splitting)
-    {
-        std::chrono::steady_clock::time_point start = std::chrono::steady_clock::now();
-        combustion->correct();
-        //label flag_mpi_init;
-        //MPI_Initialized(&flag_mpi_init);
-        if(flag_mpi_init) MPI_Barrier(PstreamGlobals::MPI_COMM_FOAM);
-        std::chrono::steady_clock::time_point stop = std::chrono::steady_clock::now();
-        std::chrono::duration<double> processingTime = std::chrono::duration_cast<std::chrono::duration<double>>(stop - start);
-        time_monitor_chem += processingTime.count();
-    }
 
 #ifdef GPUSolverNew_
 #if defined DEBUG_
@@ -265,6 +269,18 @@ if(flag_mpi_init) MPI_Barrier(PstreamGlobals::MPI_COMM_FOAM);
 
     fflush(stderr);
 #else
+    if (!splitting)
+    {
+        std::chrono::steady_clock::time_point start = std::chrono::steady_clock::now();
+        combustion->correct();
+        //label flag_mpi_init;
+        //MPI_Initialized(&flag_mpi_init);
+        if(flag_mpi_init) MPI_Barrier(PstreamGlobals::MPI_COMM_FOAM);
+        std::chrono::steady_clock::time_point stop = std::chrono::steady_clock::now();
+        std::chrono::duration<double> processingTime = std::chrono::duration_cast<std::chrono::duration<double>>(stop - start);
+        time_monitor_chem += processingTime.count();
+    }
+
     start2 = std::clock();
     volScalarField Yt(0.0*Y[0]);
     int speciesIndex = 0;

applications/solvers/dfLowMachFoam/createGPUSolver.H

@@ -6,9 +6,10 @@ int nRanks, myRank, localRank, mpi_init_flag = 0;
 
 dfMatrixDataBase dfDataBase;
 dfThermo thermo_GPU(dfDataBase);
+dfChemistrySolver chemistrySolver_GPU(dfDataBase);
 dfRhoEqn rhoEqn_GPU(dfDataBase);
 dfUEqn UEqn_GPU(dfDataBase);
-dfYEqn YEqn_GPU(dfDataBase);
+dfYEqn YEqn_GPU(dfDataBase, chemistrySolver_GPU);
 dfEEqn EEqn_GPU(dfDataBase, thermo_GPU);
 dfpEqn pEqn_GPU(dfDataBase);
 
@@ -325,7 +326,17 @@ void createGPUBase(const IOdictionary& CanteraTorchProperties, fvMesh& mesh, Ptr
 
     dfDataBase.createConstantFieldsInternal();
     dfDataBase.createConstantFieldsBoundary();
-    dfDataBase.initConstantFieldsInternal(&mesh.Sf()[0][0], &mesh.magSf()[0], &mesh.surfaceInterpolation::weights()[0], &mesh.deltaCoeffs()[0], &mesh.V()[0]);
+
+    // construct mesh distance for limitedLinear scheme
+    vectorField meshDistance = mesh.Sf();
+    forAll(meshDistance, facei) {
+        label own = owner[facei];
+        label nei = neighbour[facei];
+        meshDistance[facei] = mesh.C()[nei] - mesh.C()[own];
+    }
+
+    dfDataBase.initConstantFieldsInternal(&mesh.Sf()[0][0], &mesh.magSf()[0], &mesh.surfaceInterpolation::weights()[0], 
+            &mesh.deltaCoeffs()[0], &mesh.V()[0], &meshDistance[0][0]);
     dfDataBase.initConstantFieldsBoundary(boundary_sf, boundary_mag_sf, boundary_delta_coeffs, boundary_weights, boundary_face_cell, 
             patch_type_calculated, patch_type_extropolated);
 
@@ -586,7 +597,7 @@ void createGPUpEqn(const IOdictionary& CanteraTorchProperties, volScalarField& p
 
 void createGPUThermo(const IOdictionary& CanteraTorchProperties, volScalarField& T, volScalarField& he, 
         const volScalarField& psi, const volScalarField& alpha, const volScalarField& mu,
-        const volScalarField& K, const volScalarField& dpdt) {
+        const volScalarField& K, const volScalarField& dpdt, dfChemistryModel<basicThermo>* chemistry) {
     DEBUG_TRACE;
     // initialize dfThermo
     string mechanismFile;
@@ -604,6 +615,8 @@ void createGPUThermo(const IOdictionary& CanteraTorchProperties, volScalarField&
     double *h_boundary_thermo_alpha = new double[dfDataBase.num_boundary_surfaces];
     double *h_boundary_mu = new double[dfDataBase.num_boundary_surfaces];
     double *h_boundary_k = new double[dfDataBase.num_boundary_surfaces];
+    double *h_boundary_thermo_rhoD = new double[dfDataBase.num_boundary_surfaces * dfDataBase.num_species];
+    double *h_thermo_rhoD = new double[dfDataBase.num_cells * dfDataBase.num_species];
 
     // initialize thermo boundary
     std::vector<int> patch_type_T(dfDataBase.num_patches);
@@ -651,6 +664,14 @@ void createGPUThermo(const IOdictionary& CanteraTorchProperties, volScalarField&
                     dynamic_cast<const processorFvPatchField<scalar>&>(patchK).patchInternalField()();
             memcpy(h_boundary_k + offset + patchsize, &patchKInternal[0], patchsize * sizeof(double));
 
+            for (int i = 0; i < dfDataBase.num_species; i++) {
+                const fvPatchScalarField& patchRhoD = chemistry->rhoD(i).boundaryField()[patchi];
+                memcpy(h_boundary_thermo_rhoD + i * dfDataBase.num_boundary_surfaces + offset, &patchRhoD[0], patchsize * sizeof(double));
+                scalarField patchRhoDInternal = 
+                        dynamic_cast<const processorFvPatchField<scalar>&>(patchRhoD).patchInternalField()();
+                memcpy(h_boundary_thermo_rhoD + i * dfDataBase.num_boundary_surfaces + offset + patchsize, &patchRhoDInternal[0], patchsize * sizeof(double));
+            }
+
             offset += patchsize * 2;
         } else {
             memcpy(h_boundary_T + offset, &patchT[0], patchsize * sizeof(double));
@@ -660,13 +681,29 @@ void createGPUThermo(const IOdictionary& CanteraTorchProperties, volScalarField&
             memcpy(h_boundary_mu + offset, &patchMu[0], patchsize * sizeof(double));
             memcpy(h_boundary_k + offset, &patchK[0], patchsize * sizeof(double));
 
+            for (int i = 0; i < dfDataBase.num_species; i++) {
+                const fvPatchScalarField& patchRhoD = chemistry->rhoD(i).boundaryField()[patchi];
+                memcpy(h_boundary_thermo_rhoD + i * dfDataBase.num_boundary_surfaces + offset, &patchRhoD[0], patchsize * sizeof(double));
+            }
             offset += patchsize;
         }
     }
+    for (int i = 0; i < dfDataBase.num_species; i++) {
+        memcpy(h_thermo_rhoD + i * dfDataBase.num_cells, &chemistry->rhoD(i)[0], dfDataBase.num_cells * sizeof(double));
+    }
+    double *h_T = dfDataBase.getFieldPointer("T", location::cpu, position::internal);
+    memcpy(h_T, &T[0], dfDataBase.cell_value_bytes);
+
     thermo_GPU.setConstantFields(patch_type_T);
-    thermo_GPU.initNonConstantFields(&T[0], &he[0], &psi[0], &alpha[0], &mu[0], &K[0], &dpdt[0],
-            h_boundary_T, h_boundary_he, h_boundary_thermo_psi, h_boundary_thermo_alpha, h_boundary_mu, h_boundary_k);
+    thermo_GPU.initNonConstantFields(h_T, &he[0], &psi[0], &alpha[0], &mu[0], &K[0], &dpdt[0], h_thermo_rhoD,
+            h_boundary_T, h_boundary_he, h_boundary_thermo_psi, h_boundary_thermo_alpha, h_boundary_mu, h_boundary_k, h_boundary_thermo_rhoD);
 
     delete h_boundary_T;
     delete h_boundary_he;
+    delete h_boundary_thermo_psi;
+    delete h_boundary_thermo_alpha;
+    delete h_boundary_mu;
+    delete h_boundary_k;
+    delete h_boundary_thermo_rhoD;
+    delete h_thermo_rhoD;
 }

applications/solvers/dfLowMachFoam/dfLowMachFoam.C

@@ -76,6 +76,7 @@ Description
     #include "dfMatrixOpBase.H"
     #include "dfNcclBase.H"
     #include "dfThermo.H"
+    #include "dfChemistrySolver.H"
     #include <cuda_runtime.h>
     #include <thread>
 
@@ -87,9 +88,13 @@ Description
     #include "upwind.H"
     #include "GenFvMatrix.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
@@ -216,7 +221,11 @@ int main(int argc, char *argv[])
 
     const volScalarField& mu = thermo.mu();
     const volScalarField& alpha = thermo.alpha();
-    createGPUThermo(CanteraTorchProperties, T, thermo.he(), psi, alpha, mu, K, dpdt);
+    createGPUThermo(CanteraTorchProperties, T, thermo.he(), psi, alpha, mu, K, dpdt, chemistry);
+    if (chemistry->ifChemstry())
+    {
+        chemistrySolver_GPU.setConstantValue(dfDataBase.num_cells, dfDataBase.num_species, 4096);
+    }
 #endif
 
     end1 = std::clock();
@@ -302,13 +311,15 @@ int main(int argc, char *argv[])
                 thermo_GPU.sync();
             #if defined DEBUG_
                 // check correctThermo
+                int speciesIndex = 6;
                 chemistry->correctThermo(); // reference
                 double *h_boundary_T_tmp = new double[dfDataBase.num_boundary_surfaces];
                 double *h_boundary_he_tmp = new double[dfDataBase.num_boundary_surfaces];
                 double *h_boundary_mu_tmp = new double[dfDataBase.num_boundary_surfaces];
                 double *h_boundary_rho_tmp = new double[dfDataBase.num_boundary_surfaces];
                 double *h_boundary_thermo_alpha_tmp = new double[dfDataBase.num_boundary_surfaces];    
                 double *h_boundary_thermo_psi_tmp = new double[dfDataBase.num_boundary_surfaces];
+                double *h_boundary_thermo_rhoD_tmp = new double[dfDataBase.num_boundary_surfaces];
                 offset = 0;
                 forAll(T.boundaryField(), patchi)
                 {
@@ -318,6 +329,7 @@ int main(int argc, char *argv[])
                     const fvPatchScalarField& patchAlpha = alpha.boundaryField()[patchi];
                     const fvPatchScalarField& patchRho = thermo.rho()().boundaryField()[patchi];
                     const fvPatchScalarField& patchT = T.boundaryField()[patchi];
+                    const fvPatchScalarField& patchRhoD = chemistry->rhoD(speciesIndex).boundaryField()[patchi];
 
                     int patchsize = patchT.size();
                     if (patchT.type() == "processor") {
@@ -351,6 +363,11 @@ int main(int argc, char *argv[])
                                 dynamic_cast<const processorFvPatchField<scalar>&>(patchRho).patchInternalField()();
                         memcpy(h_boundary_rho_tmp + offset + patchsize, &patchRhoInternal[0], patchsize * sizeof(double));
 
+                        memcpy(h_boundary_thermo_rhoD_tmp + offset, &patchRhoD[0], patchsize * sizeof(double));
+                        scalarField patchRhoDInternal = 
+                                dynamic_cast<const processorFvPatchField<scalar>&>(patchRhoD).patchInternalField()();
+                        memcpy(h_boundary_thermo_rhoD_tmp + offset + patchsize, &patchRhoDInternal[0], patchsize * sizeof(double));
+
                         offset += patchsize * 2;
                     } else {
                         memcpy(h_boundary_T_tmp + offset, &patchT[0], patchsize * sizeof(double));
@@ -359,6 +376,7 @@ int main(int argc, char *argv[])
                         memcpy(h_boundary_thermo_alpha_tmp + offset, &patchAlpha[0], patchsize * sizeof(double));
                         memcpy(h_boundary_mu_tmp + offset, &patchMu[0], patchsize * sizeof(double));
                         memcpy(h_boundary_rho_tmp + offset, &patchRho[0], patchsize * sizeof(double));
+                        memcpy(h_boundary_thermo_rhoD_tmp + offset, &patchRhoD[0], patchsize * sizeof(double));
 
                         offset += patchsize;
                     }
@@ -369,12 +387,13 @@ int main(int argc, char *argv[])
                     MPI_Comm_rank(MPI_COMM_WORLD, &rank);
                 }
                 if (!mpi_init_flag || rank == 0) {
-                    thermo_GPU.compareT(&T[0], h_boundary_T_tmp, printFlag);
+                    // thermo_GPU.compareT(&T[0], h_boundary_T_tmp, printFlag);
                     // thermo_GPU.compareHe(&thermo.he()[0], h_boundary_he_tmp, printFlag);
                     // thermo_GPU.comparePsi(&psi[0], h_boundary_thermo_psi_tmp, printFlag);
                     // thermo_GPU.compareAlpha(&alpha[0], h_boundary_thermo_alpha_tmp, printFlag);
                     // thermo_GPU.compareMu(&mu[0], h_boundary_mu_tmp, printFlag);
                     // thermo_GPU.compareRho(&thermo.rho()()[0], h_boundary_rho_tmp, printFlag);
+                    // thermo_GPU.compareRhoD(&chemistry->rhoD(speciesIndex)[0], h_boundary_thermo_rhoD_tmp, speciesIndex, printFlag);
                 }
 
                 delete h_boundary_T_tmp;
@@ -396,11 +415,10 @@ int main(int argc, char *argv[])
             }
             // update T for debug
             #ifdef GPUSolverNew_
-            double *h_T_tmp = new double[dfDataBase.num_cells];
+            double *h_T = dfDataBase.getFieldPointer("T", location::cpu, position::internal);
             double *h_boundary_T_tmp = new double[dfDataBase.num_boundary_surfaces];
-            thermo_GPU.updateCPUT(h_T_tmp, h_boundary_T_tmp);
-
-            memcpy(&T[0], h_T_tmp, T.size() * sizeof(double));
+            thermo_GPU.updateCPUT(h_T, h_boundary_T_tmp);
+            memcpy(&T[0], h_T, T.size() * sizeof(double));
             offset = 0;
             forAll(T.boundaryField(), patchi) {
                 const fvPatchScalarField& const_patchT = T.boundaryField()[patchi];
@@ -414,7 +432,6 @@ int main(int argc, char *argv[])
                     offset += patchsize;
                 }
             }
-            delete h_T_tmp;
             delete h_boundary_T_tmp;
             #endif
 
@@ -469,7 +486,7 @@ int main(int argc, char *argv[])
 #endif
 
 #ifdef GPUSolverNew_
-        // write U for
+        // write U
         UEqn_GPU.postProcess();
         memcpy(&U[0][0], dfDataBase.h_u, dfDataBase.cell_value_vec_bytes);
         U.correctBoundaryConditions();

applications/solvers/dfLowMachFoam/pEqn_CPU.H

@@ -68,8 +68,7 @@ else
     fvScalarMatrix pDDtEqn
     (
         fvc::ddt(rho) + 
-        // psi*correction(fvm::ddt(p)) // TODO: bugs in correction fvMatrix
-        psi*fvm::ddt(p)
+        psi*correction(fvm::ddt(p))
         + fvc::div(phiHbyA)
     );
 

applications/solvers/dfLowMachFoam/pEqn_GPU.H

@@ -98,8 +98,7 @@ UEqn_GPU.sync();
     fvScalarMatrix pDDtEqn
     (
         fvc::ddt(rho) + 
-        // psi*correction(fvm::ddt(p)) // TODO: bugs in correction fvMatrix
-        psi*fvm::ddt(p)
+        psi*correction(fvm::ddt(p))
         + fvc::div(phiHbyA)
     );
     while (pimple.correctNonOrthogonal())

src_gpu/AmgXSolver.H

@@ -30,32 +30,6 @@
 // mpi
 #include <mpi.h>
 
-
-/** \brief A macro to check the returned CUDA error code.
- *
- * \param call [in] Function call to CUDA API.
- */
-# define CHECK(call)                                                        \
-do                                                    \
-{                                                     \
-    const cudaError_t error_code = call;              \
-    if (error_code != cudaSuccess)                    \
-    {                                                 \
-        printf("CUDA Error:\n");                      \
-        printf("    File:       %s\n", __FILE__);     \
-        printf("    Line:       %d\n", __LINE__);     \
-        printf("    Error code: %d\n", error_code);   \
-        printf("    Error text: %s\n",                \
-            cudaGetErrorString(error_code));          \
-        exit(1);                                      \
-    }                                                 \
-} while (0)
-
-
-
-
-
-
 /** \brief A wrapper class for coupling PETSc and AmgX.
  *
  * This class is a wrapper of AmgX library for PETSc. PETSc users only need to