GPU version of UEqn.A and UEqn.H

applications/solvers/dfLowMachFoam/UEqn.H

@@ -1,33 +1,4 @@
 // Solve the Momentum equation
-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();
-#ifdef CPUSolver_
-    time_monitor_UEqn += double(end1 - start1) / double(CLOCKS_PER_SEC);
-    time_monitor_UEqn_mtxAssembly += double(end1 - start1) / double(CLOCKS_PER_SEC);
-#endif
-
-#ifdef CPUSolver_
-    // 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
-
 #ifdef GPUSolver_
     start1 = std::clock();
     UEqn_GPU.fvm_ddt(&U.oldTime()[0][0]);
@@ -111,6 +82,31 @@ end1 = std::clock();
     //     K = 0.5*magSqr(U);
     // }
     // UEqn_GPU.checkValue(false);
+#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();

applications/solvers/dfLowMachFoam/YEqn.H

@@ -107,7 +107,10 @@ time_monitor_corrDiff += double(end - start) / double(CLOCKS_PER_SEC);
 
         if (i != inertIndex)
         {
-            #ifdef CPUSolver_
+            #ifdef GPUSolver_
+                YEqn_GPU.updatePsi(&Yi[0], speciesIndex);
+                Yi.correctBoundaryConditions();
+            #else
                 tmp<volScalarField> DEff = chemistry->rhoD(i) + turbulence->mut()/Sct;
                 fvScalarMatrix YiEqn
                 (
@@ -127,11 +130,6 @@ time_monitor_corrDiff += double(end - start) / double(CLOCKS_PER_SEC);
                 YiEqn.solve("Yi");
             #endif
 
-            #ifdef GPUSolver_
-                YEqn_GPU.updatePsi(&Yi[0], speciesIndex);
-                Yi.correctBoundaryConditions();
-            #endif
-
             Yi.max(0.0);
             Yt += Yi;
             ++speciesIndex;

applications/solvers/dfLowMachFoam/pEqn.H

@@ -7,9 +7,53 @@ if (!pimple.simpleRho())
 // pressure solution
 const volScalarField psip0(psi*p);
 
-volScalarField rAU(1.0/UEqn.A());
-surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
-volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
+#ifdef GPUSolver_
+    // UEqn.H()
+    volVectorField UEqn_H
+    (
+        IOobject
+        (
+            "H("+U.name()+')',
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedVector(dimensionSet(1,-2,-2,0,0,0,0), Zero),
+        extrapolatedCalculatedFvPatchScalarField::typeName
+    );
+    UEqn_GPU.H(&UEqn_H[0][0]);
+    UEqn_H.correctBoundaryConditions();
+
+    // UEqn.A()
+    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
+    );
+    UEqn_GPU.A(&UEqn_A[0]);
+    UEqn_A.correctBoundaryConditions();
+#endif
+
+#ifdef GPUSolver_
+    volScalarField rAU(1.0/UEqn_A);
+    surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
+    volVectorField HbyA(constrainHbyA(rAU*UEqn_H, U, p));
+#else
+    volScalarField rAU(1.0/UEqn.A());
+    surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
+    volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
+#endif
 
 if (pimple.nCorrPiso() <= 1)
 {
@@ -99,6 +143,7 @@ p.relax();
 U = HbyA - rAU*fvc::grad(p);
 U.correctBoundaryConditions();
 K = 0.5*magSqr(U);
+UEqn_GPU.correctPsi(&U[0][0]);
 
 if (pimple.simpleRho())
 {

applications/solvers/dfLowMachFoam/rhoEqn.H

@@ -28,18 +28,6 @@ Description
     Solve the continuity for density.
 
 \*---------------------------------------------------------------------------*/
-
-#ifdef CPUSolver_
-{
-    fvScalarMatrix rhoEqn
-    (
-        fvm::ddt(rho)
-      + fvc::div(phi)
-    );
-
-    rhoEqn.solve();
-}
-#endif
 #ifdef GPUSolver_
 {
     rho.oldTime();
@@ -57,6 +45,16 @@ Description
     rhoEqn_GPU.updatePsi(&rho.primitiveFieldRef()[0]);
     rho.correctBoundaryConditions();
 }
+#else
+{
+    fvScalarMatrix rhoEqn
+    (
+        fvm::ddt(rho)
+      + fvc::div(phi)
+    );
+
+    rhoEqn.solve();
+}
 #endif
 
 // ************************************************************************* //

src_gpu/dfUEqn.H

@@ -19,8 +19,8 @@ private:
     int *d_A_csr_row_index, *d_A_csr_diag_index, *d_A_csr_col_index;
 
     // Matrix variables
-    double *d_A_csr, *d_b, *d_psi = nullptr;
-    double *h_A_csr, *h_b, *h_psi = nullptr;
+    double *d_A_csr, *d_b, *d_psi, *d_H, *d_A, *d_ueqn_internal_coeffs, *d_ueqn_boundary_coeffs= nullptr;
+    double *h_A_csr, *h_b, *h_psi, *h_H, *h_A = nullptr;
 
 public:
     dfUEqn();
@@ -48,6 +48,13 @@ public:
 
     void add_fvMatrix(double *turbSrc_low, double *turbSrc_diag, double *turbSrc_upp, double *turbSrc_source);
 
+    void A(double *Psi);
+
+    void H(double *Psi);
+
     void solve();
+
     void updatePsi(double *Psi);
+
+    void correctPsi(double *Psi);
 };