Skip to content

Commit

Permalink
3D version with z=1
Browse files Browse the repository at this point in the history
  • Loading branch information
@mehmetyusufoglu
mehmetyusufoglu committed Aug 1, 2024
1 parent f282ecf commit 77633bc
Show file tree
Hide file tree
Showing 2 changed files with 217 additions and 156 deletions.
145 changes: 90 additions & 55 deletions benchmarks/cloverleaf/src/cloverLeafKernels.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,11 +5,12 @@
#include <experimental/mdspan>

using Data = float;
using Dim2 = alpaka::DimInt<2u>;
using Dim3 = alpaka::DimInt<3u>;
using Idx = std::uint32_t;

const Idx nx = 960; // Number of cells in x direction
const Idx nx = 1920; // Number of cells in x direction
const Idx ny = 960; // Number of cells in y direction
const Idx nz = 1; // Number of cells in z direction

// Kernel to update the halo regions
struct UpdateHaloKernel
Expand All @@ -19,20 +20,26 @@ struct UpdateHaloKernel
{
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)
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) = density(nx - 2, j);
density(i, j, k) = density(nx - 2, j, k);
if(i == nx - 1)
density(i, j) = density(1, j);
density(i, j, k) = density(1, j, k);

if(j == 0)
density(i, j) = density(i, ny - 2);
density(i, j, k) = density(i, ny - 2, k);
if(j == ny - 1)
density(i, j) = density(i, 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);
}
}
};
Expand All @@ -51,11 +58,12 @@ struct IdealGasKernel
{
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)
if(i < nx && j < ny && k < nz)
{
pressure(i, j) = (gamma - 1.0f) * density(i, j) * energy(i, j);
soundspeed(i, j) = sqrt(gamma * pressure(i, j) / density(i, j));
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));
}
}
};
Expand All @@ -70,28 +78,28 @@ struct InitializerKernel
MdSpan energy,
MdSpan pressure,
MdSpan velocityX,
MdSpan velocityY) const -> void
MdSpan velocityY,
MdSpan velocityZ) 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)
if(i < nx && j < ny && k < nz)
{
if(i >= static_cast<Idx>(0.0 * nx) && i < static_cast<Idx>(5.0 * nx / 10.0)
&& j >= static_cast<Idx>(0.0 * ny) && j < static_cast<Idx>(2.0 * ny / 10.0))
{
density(i, j) = 1.0f;
energy(i, j) = 2.5f;
}
else
density(i, j, k) = 0.2f; // State 1 density
energy(i, j, k) = 1.0f; // State 1 energy

if(i > 0 && i < 5.0f * nx / 10.0f && j > 0 && j < 2.0f * ny / 10.0f)
{
density(i, j) = 0.2f;
energy(i, j) = 1.0f;
density(i, j, k) = 1.0f; // State 2 density
energy(i, j, k) = 2.5f; // State 2 energy
}

pressure(i, j) = 1.0f;
velocityX(i, j) = 0.0f;
velocityY(i, j) = 0.0f;
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
}
}
};
Expand All @@ -107,19 +115,22 @@ struct EOSKernel
MdSpan pressure,
MdSpan velocityX,
MdSpan velocityY,
MdSpan velocityZ,
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)
if(i < nx && j < ny && k < nz)
{
// Compute pressure using ideal gas law: P = (gamma - 1) * density * energy
pressure(i, j) = (gamma - 1.0f) * density(i, j) * energy(i, j);
pressure(i, j, k) = (gamma - 1.0f) * density(i, j, k) * energy(i, j, k);

// Additional calculations to update velocities (this is a simplified example)
velocityX(i, j) += pressure(i, j) * 0.1f;
velocityY(i, j) += pressure(i, j) * 0.1f;
velocityX(i, j, k) += pressure(i, j, k) * 0.1f;
velocityY(i, j, k) += pressure(i, j, k) * 0.1f;
velocityZ(i, j, k) += pressure(i, j, k) * 0.1f;
}
}
};
Expand All @@ -135,21 +146,25 @@ struct FluxKernel
MdSpan pressure,
MdSpan velocityX,
MdSpan velocityY,
MdSpan velocityZ,
MdSpan fluxDensity,
MdSpan fluxEnergy,
MdSpan fluxVelocityX,
MdSpan fluxVelocityY) const -> void
MdSpan fluxVelocityY,
MdSpan fluxVelocityZ) 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)
if(i < nx && j < ny && k < nz)
{
// Compute fluxes (this is a simplified example)
fluxDensity(i, j) = density(i, j) * velocityX(i, j);
fluxEnergy(i, j) = energy(i, j) * velocityX(i, j);
fluxVelocityX(i, j) = velocityX(i, j) * velocityX(i, j) + pressure(i, j);
fluxVelocityY(i, j) = velocityY(i, j) * velocityX(i, j);
fluxDensity(i, j, k) = density(i, j, k) * velocityX(i, j, k);
fluxEnergy(i, j, k) = energy(i, j, k) * velocityX(i, j, k);
fluxVelocityX(i, j, k) = velocityX(i, j, k) * velocityX(i, j, k) + pressure(i, j, k);
fluxVelocityY(i, j, k) = velocityY(i, j, k) * velocityX(i, j, k);
fluxVelocityZ(i, j, k) = velocityZ(i, j, k) * velocityX(i, j, k);
}
}
};
Expand All @@ -158,15 +173,21 @@ struct FluxKernel
struct AdvectionKernel
{
template<typename TAcc, typename MdSpan>
ALPAKA_FN_ACC auto operator()(TAcc const& acc, MdSpan density, MdSpan velocityX, MdSpan velocityY) const -> void
ALPAKA_FN_ACC auto operator()(
TAcc const& acc,
MdSpan density,
MdSpan velocityX,
MdSpan velocityY,
MdSpan velocityZ) 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)
if(i < nx && j < ny && k < nz)
{
// Simple advection calculation (this is a simplified example)
density(i, j) += (velocityX(i, j) + velocityY(i, j)) * 0.01f;
density(i, j, k) += (velocityX(i, j, k) + velocityY(i, j, k) + velocityZ(i, j, k)) * 0.01f;
}
}
};
Expand All @@ -181,21 +202,25 @@ struct LagrangianKernel
MdSpan energy,
MdSpan velocityX,
MdSpan velocityY,
MdSpan velocityZ,
MdSpan fluxDensity,
MdSpan fluxEnergy,
MdSpan fluxVelocityX,
MdSpan fluxVelocityY) const -> void
MdSpan fluxVelocityY,
MdSpan fluxVelocityZ) 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)
if(i < nx && j < ny && k < nz)
{
// Update the cell-centered variables based on flux calculations
density(i, j) -= fluxDensity(i, j) * 0.1f;
energy(i, j) -= fluxEnergy(i, j) * 0.1f;
velocityX(i, j) -= fluxVelocityX(i, j) * 0.1f;
velocityY(i, j) -= fluxVelocityY(i, j) * 0.1f;
density(i, j, k) -= fluxDensity(i, j, k) * 0.1f;
energy(i, j, k) -= fluxEnergy(i, j, k) * 0.1f;
velocityX(i, j, k) -= fluxVelocityX(i, j, k) * 0.1f;
velocityY(i, j, k) -= fluxVelocityY(i, j, k) * 0.1f;
velocityZ(i, j, k) -= fluxVelocityZ(i, j, k) * 0.1f;
}
}
};
Expand All @@ -209,21 +234,24 @@ struct ViscosityKernel
MdSpan density,
MdSpan velocityX,
MdSpan velocityY,
MdSpan velocityZ,
MdSpan pressure,
MdSpan viscosity) 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 > 0 && i < nx - 1 && j > 0 && j < ny - 1)
if(i > 0 && i < nx - 1 && j > 0 && j < ny - 1 && k > 0 && k < nz - 1)
{
float gradVx = (velocityX(i + 1, j) - velocityX(i - 1, j)) * 0.5f;
float gradVy = (velocityY(i, j + 1) - velocityY(i, j - 1)) * 0.5f;
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;

viscosity(i, j) = density(i, j) * (gradVx * gradVx + gradVy * gradVy) * 0.01f;
viscosity(i, j, k) = density(i, j, k) * (gradVx * gradVx + gradVy * gradVy + gradVz * gradVz) * 0.01f;

// Apply viscosity to pressure
pressure(i, j) += viscosity(i, j);
pressure(i, j, k) += viscosity(i, j, k);
}
}
};
Expand All @@ -232,16 +260,23 @@ struct ViscosityKernel
struct MaxVelocityKernel
{
template<typename TAcc, typename MdSpan>
ALPAKA_FN_ACC auto operator()(TAcc const& acc, MdSpan velocityX, MdSpan velocityY, Data* maxVelocity) const -> void
ALPAKA_FN_ACC auto operator()(
TAcc const& acc,
MdSpan velocityX,
MdSpan velocityY,
MdSpan velocityZ,
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];
auto const k = alpaka::getIdx<alpaka::Grid, alpaka::Threads>(acc)[2];

if(i < nx && j < ny)
if(i < nx && j < ny && k < nz)
{
float vx = velocityX(i, j);
float vy = velocityY(i, j);
float v = alpaka::math::sqrt(acc, (vx * vx + vy * vy));
float vx = velocityX(i, j, k);
float vy = velocityY(i, j, k);
float vz = velocityZ(i, j, k);
float v = alpaka::math::sqrt(acc, (vx * vx + vy * vy + vz * vz));

// Atomic operation to find the maximum velocity
float val = alpaka::atomicMax(acc, maxVelocity, v);
Expand All @@ -250,10 +285,10 @@ struct MaxVelocityKernel
}
};

[[maybe_unused]] static float calculateTimeStep(float dx, float dy, float maxVelocity, float cMax)
[[maybe_unused]] static float calculateTimeStep(float dx, float dy, float dz, float maxVelocity, float cMax)
{
// Compute the smallest grid spacing
float minDx = std::min({dx, dy});
float minDx = std::min({dx, dy, dz});

// Calculate the time step based on the CFL condition
float dt = cMax * minDx / maxVelocity;
Expand Down
Loading

0 comments on commit 77633bc

Please sign in to comment.