refactor a little

benchmarks/cloverleaf/CMakeLists.txt

@@ -19,6 +19,10 @@ if(NOT TARGET alpaka::alpaka)
     endif()
 endif()
 
+if (alpaka_USE_MDSPAN STREQUAL "OFF")
+    message(STATUS "The conv2DWithMdspan example requires mdspan. Please set alpaka_USE_MDSPAN accordingly. Example disabled.")
+    return()
+endif ()
 
 set(_TARGET_NAME "cloverleaf")
 append_recursive_files_add_to_src_group("src/" "src/" "cpp" _FILES_SOURCE)

benchmarks/cloverleaf/src/cloverLeafKernels.hpp

@@ -12,6 +12,62 @@ const Idx nx = 512; // Number of cells in x direction
 const Idx ny = 512; // Number of cells in y direction
 const Idx nz = 512; // Number of cells in z direction
 
+// Kernel to update the halo regions
+struct UpdateHaloKernel
+{
+    template<typename TAcc, typename MdSpan>
+    ALPAKA_FN_ACC auto operator()(TAcc const& acc, MdSpan density) const -> void
+    {
+        auto const i = alpaka::getIdx<alpaka::Grid, alpaka::Threads>(acc)[0];
+        auto const j = alpaka::getIdx<alpaka::Grid, alpaka::Threads>(acc)[1];
+        auto const k = alpaka::getIdx<alpaka::Grid, alpaka::Threads>(acc)[2];
+
+        if(i < nx && j < ny && k < nz)
+        {
+            // Update halo cells for density (simplified example)
+            // Assuming a single layer halo, and periodic boundary conditions
+            if(i == 0)
+                density(i, j, k) = density(nx - 2, j, k);
+            if(i == nx - 1)
+                density(i, j, k) = density(1, j, k);
+
+            if(j == 0)
+                density(i, j, k) = density(i, ny - 2, k);
+            if(j == ny - 1)
+                density(i, j, k) = density(i, 1, k);
+
+            if(k == 0)
+                density(i, j, k) = density(i, j, nz - 2);
+            if(k == nz - 1)
+                density(i, j, k) = density(i, j, 1);
+        }
+    }
+};
+
+// Kernel to compute the ideal gas equation of state
+struct IdealGasKernel
+{
+    template<typename TAcc, typename MdSpan>
+    ALPAKA_FN_ACC auto operator()(
+        TAcc const& acc,
+        MdSpan density,
+        MdSpan energy,
+        MdSpan pressure,
+        MdSpan soundspeed,
+        float gamma) const -> void
+    {
+        auto const i = alpaka::getIdx<alpaka::Grid, alpaka::Threads>(acc)[0];
+        auto const j = alpaka::getIdx<alpaka::Grid, alpaka::Threads>(acc)[1];
+        auto const k = alpaka::getIdx<alpaka::Grid, alpaka::Threads>(acc)[2];
+
+        if(i < nx && j < ny && k < nz)
+        {
+            pressure(i, j, k) = (gamma - 1.0f) * density(i, j, k) * energy(i, j, k);
+            soundspeed(i, j, k) = sqrt(gamma * pressure(i, j, k) / density(i, j, k));
+        }
+    }
+};
+
 // Kernel to initialize the simulation variables
 struct InitializerKernel
 {
@@ -25,19 +81,22 @@ struct InitializerKernel
         MdSpan velocityY,
         MdSpan velocityZ) const -> void
     {
-        // Get thread index, the center of filter-matrix is positioned to the item on this index.
         auto const i = alpaka::getIdx<alpaka::Grid, alpaka::Threads>(acc)[0];
         auto const j = alpaka::getIdx<alpaka::Grid, alpaka::Threads>(acc)[1];
         auto const k = alpaka::getIdx<alpaka::Grid, alpaka::Threads>(acc)[2];
 
         if(i < nx && j < ny && k < nz)
         {
             density(i, j, k) = 1.0f; // Initial density
-            energy(i, j, k) = 1.0f; // Initial energy
+            energy(i, j, k) = 2.5f; // Initial energy
             pressure(i, j, k) = 1.0f; // Initial pressure
-            velocityX(i, j, k) = 0.0f; // Initial velocity in x direction
-            velocityY(i, j, k) = 0.0f; // Initial velocity in y direction
-            velocityZ(i, j, k) = 0.0f; // Initial velocity in z direction
+            velocityX(i, j, k) = 0.0f; // Initial velocity
+
+            if(i < nx && j < ny && k < nz)
+            {
+                // Simple advection calculation (this is a simplified example)
+                density(i, j, k) += (velocityX(i, j, k) + velocityY(i, j, k) + velocityZ(i, j, k)) * 0.01f;
+            }
         }
     }
 };
@@ -130,6 +189,7 @@ struct AdvectionKernel
     }
 };
 
+// Kernel for the Lagrangian step
 struct LagrangianKernel
 {
     template<typename TAcc, typename MdSpan>
@@ -162,6 +222,7 @@ struct LagrangianKernel
     }
 };
 
+// Kernel for viscosity calculations
 struct ViscosityKernel
 {
     template<typename TAcc, typename MdSpan>
@@ -178,9 +239,8 @@ struct ViscosityKernel
         auto const j = alpaka::getIdx<alpaka::Grid, alpaka::Threads>(acc)[1];
         auto const k = alpaka::getIdx<alpaka::Grid, alpaka::Threads>(acc)[2];
 
-        if(i < nx && j < ny && k < nz)
+        if(i > 0 && i < nx - 1 && j > 0 && j < ny - 1 && k > 0 && k < nz - 1)
         {
-            // Calculate artificial viscosity (this is a simplified example)
             float gradVx = (velocityX(i + 1, j, k) - velocityX(i - 1, j, k)) * 0.5f;
             float gradVy = (velocityY(i, j + 1, k) - velocityY(i, j - 1, k)) * 0.5f;
             float gradVz = (velocityZ(i, j, k + 1) - velocityZ(i, j, k - 1)) * 0.5f;
@@ -193,6 +253,7 @@ struct ViscosityKernel
     }
 };
 
+// Kernel to find the maximum velocity
 struct MaxVelocityKernel
 {
     template<typename TAcc, typename MdSpan>
@@ -201,7 +262,7 @@ struct MaxVelocityKernel
         MdSpan velocityX,
         MdSpan velocityY,
         MdSpan velocityZ,
-        float* maxVelocity) const -> void
+        Data* maxVelocity) const -> void
     {
         auto const i = alpaka::getIdx<alpaka::Grid, alpaka::Threads>(acc)[0];
         auto const j = alpaka::getIdx<alpaka::Grid, alpaka::Threads>(acc)[1];
@@ -215,7 +276,8 @@ struct MaxVelocityKernel
             float v = alpaka::math::sqrt(acc, (vx * vx + vy * vy + vz * vz));
 
             // Atomic operation to find the maximum velocity
-            alpaka::atomicMax(acc, maxVelocity, v);
+            float val = alpaka::atomicMax(acc, maxVelocity, v);
+            maxVelocity[0] = val;
         }
     }
 };