@@ -12,6 +12,62 @@ const Idx nx = 512; // Number of cells in x direction
12
12
const Idx ny = 512 ; // Number of cells in y direction
13
13
const Idx nz = 512 ; // Number of cells in z direction
14
14
15
+ // Kernel to update the halo regions
16
+ struct UpdateHaloKernel
17
+ {
18
+ template <typename TAcc, typename MdSpan>
19
+ ALPAKA_FN_ACC auto operator ()(TAcc const & acc, MdSpan density) const -> void
20
+ {
21
+ auto const i = alpaka::getIdx<alpaka::Grid, alpaka::Threads>(acc)[0 ];
22
+ auto const j = alpaka::getIdx<alpaka::Grid, alpaka::Threads>(acc)[1 ];
23
+ auto const k = alpaka::getIdx<alpaka::Grid, alpaka::Threads>(acc)[2 ];
24
+
25
+ if (i < nx && j < ny && k < nz)
26
+ {
27
+ // Update halo cells for density (simplified example)
28
+ // Assuming a single layer halo, and periodic boundary conditions
29
+ if (i == 0 )
30
+ density (i, j, k) = density (nx - 2 , j, k);
31
+ if (i == nx - 1 )
32
+ density (i, j, k) = density (1 , j, k);
33
+
34
+ if (j == 0 )
35
+ density (i, j, k) = density (i, ny - 2 , k);
36
+ if (j == ny - 1 )
37
+ density (i, j, k) = density (i, 1 , k);
38
+
39
+ if (k == 0 )
40
+ density (i, j, k) = density (i, j, nz - 2 );
41
+ if (k == nz - 1 )
42
+ density (i, j, k) = density (i, j, 1 );
43
+ }
44
+ }
45
+ };
46
+
47
+ // Kernel to compute the ideal gas equation of state
48
+ struct IdealGasKernel
49
+ {
50
+ template <typename TAcc, typename MdSpan>
51
+ ALPAKA_FN_ACC auto operator ()(
52
+ TAcc const & acc,
53
+ MdSpan density,
54
+ MdSpan energy,
55
+ MdSpan pressure,
56
+ MdSpan soundspeed,
57
+ float gamma) const -> void
58
+ {
59
+ auto const i = alpaka::getIdx<alpaka::Grid, alpaka::Threads>(acc)[0 ];
60
+ auto const j = alpaka::getIdx<alpaka::Grid, alpaka::Threads>(acc)[1 ];
61
+ auto const k = alpaka::getIdx<alpaka::Grid, alpaka::Threads>(acc)[2 ];
62
+
63
+ if (i < nx && j < ny && k < nz)
64
+ {
65
+ pressure (i, j, k) = (gamma - 1 .0f ) * density (i, j, k) * energy (i, j, k);
66
+ soundspeed (i, j, k) = sqrt (gamma * pressure (i, j, k) / density (i, j, k));
67
+ }
68
+ }
69
+ };
70
+
15
71
// Kernel to initialize the simulation variables
16
72
struct InitializerKernel
17
73
{
@@ -25,19 +81,22 @@ struct InitializerKernel
25
81
MdSpan velocityY,
26
82
MdSpan velocityZ) const -> void
27
83
{
28
- // Get thread index, the center of filter-matrix is positioned to the item on this index.
29
84
auto const i = alpaka::getIdx<alpaka::Grid, alpaka::Threads>(acc)[0 ];
30
85
auto const j = alpaka::getIdx<alpaka::Grid, alpaka::Threads>(acc)[1 ];
31
86
auto const k = alpaka::getIdx<alpaka::Grid, alpaka::Threads>(acc)[2 ];
32
87
33
88
if (i < nx && j < ny && k < nz)
34
89
{
35
90
density (i, j, k) = 1 .0f ; // Initial density
36
- energy (i, j, k) = 1 . 0f ; // Initial energy
91
+ energy (i, j, k) = 2 . 5f ; // Initial energy
37
92
pressure (i, j, k) = 1 .0f ; // Initial pressure
38
- velocityX (i, j, k) = 0 .0f ; // Initial velocity in x direction
39
- velocityY (i, j, k) = 0 .0f ; // Initial velocity in y direction
40
- velocityZ (i, j, k) = 0 .0f ; // Initial velocity in z direction
93
+ velocityX (i, j, k) = 0 .0f ; // Initial velocity
94
+
95
+ if (i < nx && j < ny && k < nz)
96
+ {
97
+ // Simple advection calculation (this is a simplified example)
98
+ density (i, j, k) += (velocityX (i, j, k) + velocityY (i, j, k) + velocityZ (i, j, k)) * 0 .01f ;
99
+ }
41
100
}
42
101
}
43
102
};
@@ -130,6 +189,7 @@ struct AdvectionKernel
130
189
}
131
190
};
132
191
192
+ // Kernel for the Lagrangian step
133
193
struct LagrangianKernel
134
194
{
135
195
template <typename TAcc, typename MdSpan>
@@ -162,6 +222,7 @@ struct LagrangianKernel
162
222
}
163
223
};
164
224
225
+ // Kernel for viscosity calculations
165
226
struct ViscosityKernel
166
227
{
167
228
template <typename TAcc, typename MdSpan>
@@ -178,9 +239,8 @@ struct ViscosityKernel
178
239
auto const j = alpaka::getIdx<alpaka::Grid, alpaka::Threads>(acc)[1 ];
179
240
auto const k = alpaka::getIdx<alpaka::Grid, alpaka::Threads>(acc)[2 ];
180
241
181
- if (i < nx && j < ny && k < nz)
242
+ if (i > 0 && i < nx - 1 && j > 0 && j < ny - 1 && k > 0 && k < nz - 1 )
182
243
{
183
- // Calculate artificial viscosity (this is a simplified example)
184
244
float gradVx = (velocityX (i + 1 , j, k) - velocityX (i - 1 , j, k)) * 0 .5f ;
185
245
float gradVy = (velocityY (i, j + 1 , k) - velocityY (i, j - 1 , k)) * 0 .5f ;
186
246
float gradVz = (velocityZ (i, j, k + 1 ) - velocityZ (i, j, k - 1 )) * 0 .5f ;
@@ -193,6 +253,7 @@ struct ViscosityKernel
193
253
}
194
254
};
195
255
256
+ // Kernel to find the maximum velocity
196
257
struct MaxVelocityKernel
197
258
{
198
259
template <typename TAcc, typename MdSpan>
@@ -201,7 +262,7 @@ struct MaxVelocityKernel
201
262
MdSpan velocityX,
202
263
MdSpan velocityY,
203
264
MdSpan velocityZ,
204
- float * maxVelocity) const -> void
265
+ Data * maxVelocity) const -> void
205
266
{
206
267
auto const i = alpaka::getIdx<alpaka::Grid, alpaka::Threads>(acc)[0 ];
207
268
auto const j = alpaka::getIdx<alpaka::Grid, alpaka::Threads>(acc)[1 ];
@@ -215,7 +276,8 @@ struct MaxVelocityKernel
215
276
float v = alpaka::math::sqrt (acc, (vx * vx + vy * vy + vz * vz));
216
277
217
278
// Atomic operation to find the maximum velocity
218
- alpaka::atomicMax (acc, maxVelocity, v);
279
+ float val = alpaka::atomicMax (acc, maxVelocity, v);
280
+ maxVelocity[0 ] = val;
219
281
}
220
282
}
221
283
};
0 commit comments