Skip to content

Commit 3fec553

Browse files
tim-pearsontpadioleau
tim-pearson
and
tpadioleau
authored
Vec godunov (#23)
* running test prim_to_cons, (pre cleanup) * gtest for small N=16 testcase prim_to_cons * PrimToConsVectorized benchmark * PrimToConsVectorized speedup against scalar over plateau * PrimToConsWorstRem benchmark added * kokkos version 4.7.1 -> 5.0.0 * Seprate global utils for benchmark and tests * include utils directory * pass cons/prim arrays to prim_to_cons kernel, updated plots * simd compatible store and load functions * to_cons optimized, negligible slowdown comapred to fully inlined version * forcw inline for load, to_cons, store * inlined vectorized prim_to_cons passing tests * prim_to_cons_vec passing tests with remainders * constexp for scalar remainder, cons_to_prim implementation * cons_to_prim passing tests * neglible difference cons_to_prim vect<->scalar * prim_to_cons + cons_to_prim vectorized * first implementation time_step_vec, passing tests for skx * optimization: max reduction + division outside kernal loop * use reduction_idenitiy time step cuda compatibilty * vectorized godunov with benchmark + test * hllc select for conditionals * remove constexpr for remainders, add godunov remainder test * update setup scripts a100 skx, vectorized simulation, global utils * remove constexpr for remainder, to create PR * only use OMP_PLACES=numa_domains * Update euler_operators/cons_to_prim.hpp Co-authored-by: Thomas Padioleau <thomas.padioleau@cea.fr> * remove test/ dir duplicate, clang tidy + format * shellcheck, pylint * compare helper for cons and prim arrays in tests * (pylint):pip install pandas, (clang-tidy): widening cast fix * Udpate adastra setup scripts and plot.py * MI300 4871378 benchmark with current run_bench.sh * Add run_test.sh for mi300 * Optimized hllc and godunov implementations * Updated Godunov tests for opti * Formatting * Remove casts in hllc * Update bm save paths, consistency adastra<->ruche * Use common mapping + data_handle in godunov * Flatten loop for compare() in test_godunov * Uncomment prepare scripts sections, format hllc headers * Use APU Kokkos Arch flag for mi300a * linting and iwyu * kokkos simd include bm euler_simulation fro iwyu * Try Kokkos_SIMD.hpp for iwyu * test_godunov.cpp update inlcude for iwyu * benchmark_godunov.cpp update inlcude for iwyu * godunov.hpp update inlcude for iwyu * benchmark_euler_simulation.cpp update include for iwyu * benchmark_godunov.cpp update include for iwyu * test_godunov.cpp update include for iwyu * hllc.hpp update include for iwyu * Specify 8 to --parallel flag, remove old modules, fix inlcudes in hllc * Remove inti setup scripts --------- Co-authored-by: Thomas Padioleau <thomas.padioleau@cea.fr>
1 parent 9d5e028 commit 3fec553

17 files changed

Lines changed: 734 additions & 35 deletions

‎benchmarks/benchmark_euler_simulation.cpp‎

Lines changed: 49 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -70,6 +70,55 @@ void EulerSimulation(benchmark::State& state)
7070
set_constant_cells_processed(state, size(cons_arrays));
7171
}
7272

73+
void EulerSimulationVectorized(benchmark::State& state)
74+
{
75+
auto const nx = int_cast<index_t>(state.range());
76+
real_t const cfl_factor = 0.49;
77+
real_t const gamma = 1.4;
78+
auto nn = static_cast<std::size_t>(nx) + 2;
79+
std::size_t const n3 = nn * nn * nn;
80+
81+
real_t const dx = 1. / static_cast<real_t>(nx);
82+
PerfectGas<real_t> const eos(gamma);
83+
UniformMesh3d<real_t> const mesh(dx, dx, dx);
84+
hllc const riemann_solver;
85+
Kokkos::DefaultExecutionSpace const exec_space;
86+
EulerPrimArrays const prims_alloc = create_prim_arrays_1d<real_t>(exec_space, n3);
87+
EulerPrimArrays const prim_arrays = to_mdspan<Kokkos::mdspan<
88+
real_t,
89+
Kokkos::dextents<index_t, 3>,
90+
Kokkos::layout_left>>(prims_alloc, nx + 2, nx + 2, nx + 2);
91+
EulerConsArrays const cons_alloc = create_cons_arrays_1d<real_t>(exec_space, n3);
92+
EulerConsArrays const cons_arrays = to_mdspan<Kokkos::mdspan<
93+
real_t,
94+
Kokkos::dextents<index_t, 3>,
95+
Kokkos::layout_left>>(cons_alloc, nx + 2, nx + 2, nx + 2);
96+
97+
init_implode(exec_space, prim_arrays, mesh);
98+
prim_to_cons_vec(exec_space, as_const(prim_arrays), cons_arrays, eos);
99+
exec_space.fence();
100+
101+
for ([[maybe_unused]] auto _ : state) {
102+
real_t const dt = time_step_vec(exec_space, as_const(prim_arrays), eos, mesh);
103+
104+
godunov_vec(
105+
exec_space,
106+
as_const(prim_arrays),
107+
cons_arrays,
108+
eos,
109+
mesh,
110+
riemann_solver,
111+
cfl_factor * dt);
112+
113+
boundary_conditions_periodic(exec_space, cons_arrays, 1);
114+
115+
cons_to_prim_vec(exec_space, as_const(cons_arrays), prim_arrays, eos);
116+
exec_space.fence();
117+
}
118+
set_constant_cells_processed(state, size(cons_arrays));
119+
}
120+
73121
} // namespace
74122

75123
BENCHMARK(EulerSimulation)->UseRealTime()->DenseRange(8, 31, 8)->DenseRange(32, 320, 32);
124+
BENCHMARK(EulerSimulationVectorized)->UseRealTime()->DenseRange(8, 31, 8)->DenseRange(32, 320, 32);

‎benchmarks/benchmark_godunov.cpp‎

Lines changed: 120 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -53,6 +53,125 @@ void Godunov(benchmark::State& state)
5353
set_constant_bytes_processed(state, size_bytes(prim_arrays) + (2 * size_bytes(cons_arrays)));
5454
}
5555

56+
void GodunovWorstRem(benchmark::State& state)
57+
{
58+
auto const n = int_cast<index_t>(state.range() + 2);
59+
PerfectGas<real_t> const eos(1.4);
60+
UniformMesh3d<real_t> const mesh(1., 1., 1.);
61+
real_t const dt = 1E-9;
62+
auto nn = static_cast<std::size_t>(n);
63+
std::size_t const n3 = nn * nn * nn;
64+
65+
Kokkos::DefaultExecutionSpace const exec_space;
66+
EulerPrimArrays const prims_alloc = create_prim_arrays_1d<real_t>(exec_space, n3);
67+
EulerPrimArrays const prim_arrays = to_mdspan<Kokkos::mdspan<
68+
real_t,
69+
Kokkos::dextents<index_t, 3>,
70+
Kokkos::layout_left>>(prims_alloc, n, n, n);
71+
EulerConsArrays const cons_alloc = create_cons_arrays_1d<real_t>(exec_space, n3);
72+
EulerConsArrays const cons_arrays = to_mdspan<Kokkos::mdspan<
73+
real_t,
74+
Kokkos::dextents<index_t, 3>,
75+
Kokkos::layout_left>>(cons_alloc, n, n, n);
76+
EulerPrim<real_t> const prim {.d = 1, .p = 1, .ux0 = 0, .ux1 = 0, .ux2 = 0};
77+
init_from_state(exec_space, prim_arrays, prim);
78+
init_from_state(exec_space, cons_arrays, to_cons(prim, eos.internal_energy(prim.d, prim.p)));
79+
exec_space.fence();
80+
81+
for ([[maybe_unused]] auto _ : state) {
82+
godunov(exec_space, as_const(prim_arrays), cons_arrays, eos, mesh, hllc(), dt);
83+
exec_space.fence();
84+
benchmark::ClobberMemory();
85+
}
86+
87+
set_constant_cells_processed(state, size(cons_arrays));
88+
set_constant_bytes_processed(state, size_bytes(prim_arrays) + (2 * size_bytes(cons_arrays)));
89+
}
90+
91+
void GodunovVectorized(benchmark::State& state)
92+
{
93+
auto const n = int_cast<index_t>(state.range() + 2);
94+
PerfectGas<real_t> const eos(1.4);
95+
UniformMesh3d<real_t> const mesh(1., 1., 1.);
96+
real_t const dt = 1E-9;
97+
auto nn = static_cast<std::size_t>(n);
98+
std::size_t const n3 = nn * nn * nn;
99+
100+
Kokkos::DefaultExecutionSpace const exec_space;
101+
EulerPrimArrays const prims_alloc = create_prim_arrays_1d<real_t>(exec_space, n3);
102+
EulerPrimArrays const prim_arrays = to_mdspan<Kokkos::mdspan<
103+
real_t,
104+
Kokkos::dextents<index_t, 3>,
105+
Kokkos::layout_left>>(prims_alloc, n, n, n);
106+
EulerConsArrays const cons_alloc = create_cons_arrays_1d<real_t>(exec_space, n3);
107+
EulerConsArrays const cons_arrays = to_mdspan<Kokkos::mdspan<
108+
real_t,
109+
Kokkos::dextents<index_t, 3>,
110+
Kokkos::layout_left>>(cons_alloc, n, n, n);
111+
EulerPrim<real_t> const prim {.d = 1, .p = 1, .ux0 = 0, .ux1 = 0, .ux2 = 0};
112+
init_from_state(exec_space, prim_arrays, prim);
113+
init_from_state(exec_space, cons_arrays, to_cons(prim, eos.internal_energy(prim.d, prim.p)));
114+
exec_space.fence();
115+
116+
for ([[maybe_unused]] auto _ : state) {
117+
godunov_vec(exec_space, as_const(prim_arrays), cons_arrays, eos, mesh, hllc(), dt);
118+
exec_space.fence();
119+
benchmark::ClobberMemory();
120+
}
121+
122+
set_constant_cells_processed(state, size(cons_arrays));
123+
set_constant_bytes_processed(state, size_bytes(prim_arrays) + (2 * size_bytes(cons_arrays)));
124+
}
125+
126+
void GodunovVectorizedWorstRem(benchmark::State& state)
127+
{
128+
auto const n = int_cast<index_t>(state.range() + 2);
129+
PerfectGas<real_t> const eos(1.4);
130+
UniformMesh3d<real_t> const mesh(1., 1., 1.);
131+
real_t const dt = 1E-9;
132+
133+
auto nn = static_cast<std::size_t>(n);
134+
std::size_t const n3 = nn * nn * nn;
135+
136+
Kokkos::DefaultExecutionSpace const exec_space;
137+
EulerPrimArrays const prims_alloc = create_prim_arrays_1d<real_t>(exec_space, n3);
138+
EulerPrimArrays const prim_arrays = to_mdspan<Kokkos::mdspan<
139+
real_t,
140+
Kokkos::dextents<index_t, 3>,
141+
Kokkos::layout_left>>(prims_alloc, n, n, n);
142+
EulerConsArrays const cons_alloc = create_cons_arrays_1d<real_t>(exec_space, n3);
143+
EulerConsArrays const cons_arrays = to_mdspan<Kokkos::mdspan<
144+
real_t,
145+
Kokkos::dextents<index_t, 3>,
146+
Kokkos::layout_left>>(cons_alloc, n, n, n);
147+
EulerPrim<real_t> const prim {.d = 1, .p = 1, .ux0 = 0, .ux1 = 0, .ux2 = 0};
148+
init_from_state(exec_space, prim_arrays, prim);
149+
init_from_state(exec_space, cons_arrays, to_cons(prim, eos.internal_energy(prim.d, prim.p)));
150+
exec_space.fence();
151+
152+
for ([[maybe_unused]] auto _ : state) {
153+
godunov_vec(exec_space, as_const(prim_arrays), cons_arrays, eos, mesh, hllc(), dt);
154+
exec_space.fence();
155+
benchmark::ClobberMemory();
156+
}
157+
158+
set_constant_cells_processed(state, size(cons_arrays));
159+
set_constant_bytes_processed(state, size_bytes(prim_arrays) + (2 * size_bytes(cons_arrays)));
160+
}
161+
56162
} // namespace
57163

58-
BENCHMARK(Godunov)->UseRealTime()->DenseRange(8, 31, 8)->DenseRange(32, 320, 32);
164+
BENCHMARK(Godunov)
165+
->UseRealTime()
166+
->DenseRange(8, 31, 8)
167+
->DenseRange(32, 127, 16)
168+
->DenseRange(128, 320, 32);
169+
170+
BENCHMARK(GodunovVectorized)
171+
->UseRealTime()
172+
->DenseRange(8, 31, 8)
173+
->DenseRange(32, 127, 16)
174+
->DenseRange(128, 320, 32);
175+
176+
BENCHMARK(GodunovWorstRem)->UseRealTime()->DenseRange(7, 31, 8)->DenseRange(31, 320, 32);
177+
BENCHMARK(GodunovVectorizedWorstRem)->UseRealTime()->DenseRange(7, 31, 8)->DenseRange(31, 320, 32);

‎euler_operators/godunov.hpp‎

Lines changed: 147 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -4,6 +4,7 @@
44
#include <type_traits>
55

66
#include <Kokkos_Core.hpp>
7+
#include <Kokkos_SIMD.hpp>
78
#include <euler_arrays.hpp>
89
#include <hllc.hpp>
910
#include <perfect_gas.hpp>
@@ -28,6 +29,9 @@ void godunov(
2829
hllc const& riemann_solver,
2930
T const dt)
3031
{
32+
namespace KE = Kokkos::Experimental;
33+
using simd_t = KE::simd<T>;
34+
3135
Kokkos::Array<T, 3> const ds = {mesh.ds0(), mesh.ds1(), mesh.ds2()};
3236
T const dtodv = dt / mesh.dv();
3337

@@ -90,3 +94,146 @@ void godunov(
9094
store(cons, cons_arrays, i, j, k);
9195
});
9296
}
97+
98+
template <class SimdType, class T, class IndexType, std::size_t E0, std::size_t E1, std::size_t E2>
99+
void godunov_kernel(
100+
Kokkos::DefaultExecutionSpace const& exec_space,
101+
EulerPrimArrays<Kokkos::mdspan<
102+
T const,
103+
Kokkos::extents<IndexType, E0, E1, E2>,
104+
Kokkos::layout_left>> const& prim_arrays,
105+
EulerConsArrays<Kokkos::mdspan<
106+
T,
107+
Kokkos::extents<IndexType, E0, E1, E2>,
108+
Kokkos::layout_left>> const& cons_arrays,
109+
IndexType nx_begin,
110+
IndexType nx_end,
111+
PerfectGas<T> const& eos,
112+
UniformMesh3d<T> const& mesh,
113+
hllc const& riemann_solver,
114+
T const dt)
115+
{
116+
constexpr IndexType width = SimdType::size();
117+
IndexType const nx_blocks = (nx_end - nx_begin) / width;
118+
IndexType const ny = prim_arrays.d.extent(1);
119+
IndexType const nz = prim_arrays.d.extent(2);
120+
121+
// layout_left strides: stride in y = extent(0), stride in z = extent(0)*extent(1)
122+
IndexType const stride_y = prim_arrays.d.extent(0);
123+
IndexType const stride_z = prim_arrays.d.extent(0) * prim_arrays.d.extent(1);
124+
125+
Kokkos::Array<T, 3> const ds = {mesh.ds0(), mesh.ds1(), mesh.ds2()};
126+
T const dtodv = dt / mesh.dv();
127+
128+
Kokkos::layout_left::mapping const common_mapping = prim_arrays.d.mapping();
129+
EulerPrimArrays const prim_ptrs = data_handle(prim_arrays);
130+
EulerConsArrays const cons_ptrs = data_handle(cons_arrays);
131+
132+
Kokkos::parallel_for(
133+
"godunov_kernel",
134+
Kokkos::MDRangePolicy<
135+
Kokkos::Rank<3, Kokkos::Iterate::Left, Kokkos::Iterate::Left>,
136+
Kokkos::IndexType<IndexType>>(
137+
exec_space,
138+
{0, 1, 1},
139+
{nx_blocks, ny - 1, nz - 1}), // nx_begin already acouting for ghost cells
140+
KOKKOS_LAMBDA(IndexType const bi, IndexType const j, IndexType const k) {
141+
IndexType const base = common_mapping(nx_begin + (bi * width), j, k);
142+
EulerPrim<SimdType> const prim = load<SimdType>(prim_arrays, base);
143+
EulerFlux<SimdType> flux {};
144+
145+
{
146+
EulerPrim const prim_L = load<SimdType>(prim_arrays, base - 1);
147+
EulerPrim const prim_R = load<SimdType>(prim_arrays, base + 1);
148+
EulerFlux const flux_L = riemann_solver(dir_t<0>(), eos, prim_L, prim);
149+
EulerFlux const flux_R = riemann_solver(dir_t<0>(), eos, prim, prim_R);
150+
flux.d += ds[0] * (flux_R.d - flux_L.d);
151+
flux.e += ds[0] * (flux_R.e - flux_L.e);
152+
flux.mx0 += ds[0] * (flux_R.mx0 - flux_L.mx0);
153+
flux.mx1 += ds[0] * (flux_R.mx1 - flux_L.mx1);
154+
flux.mx2 += ds[0] * (flux_R.mx2 - flux_L.mx2);
155+
}
156+
{
157+
EulerPrim const prim_L = load<SimdType>(prim_arrays, base - stride_y);
158+
EulerPrim const prim_R = load<SimdType>(prim_arrays, base + stride_y);
159+
EulerFlux const flux_L = riemann_solver(dir_t<1>(), eos, prim_L, prim);
160+
EulerFlux const flux_R = riemann_solver(dir_t<1>(), eos, prim, prim_R);
161+
flux.d += ds[1] * (flux_R.d - flux_L.d);
162+
flux.e += ds[1] * (flux_R.e - flux_L.e);
163+
flux.mx0 += ds[1] * (flux_R.mx0 - flux_L.mx0);
164+
flux.mx1 += ds[1] * (flux_R.mx1 - flux_L.mx1);
165+
flux.mx2 += ds[1] * (flux_R.mx2 - flux_L.mx2);
166+
}
167+
{
168+
EulerPrim const prim_L = load<SimdType>(prim_arrays, base - stride_z);
169+
EulerPrim const prim_R = load<SimdType>(prim_arrays, base + stride_z);
170+
EulerFlux const flux_L = riemann_solver(dir_t<2>(), eos, prim_L, prim);
171+
EulerFlux const flux_R = riemann_solver(dir_t<2>(), eos, prim, prim_R);
172+
flux.d += ds[2] * (flux_R.d - flux_L.d);
173+
flux.e += ds[2] * (flux_R.e - flux_L.e);
174+
flux.mx0 += ds[2] * (flux_R.mx0 - flux_L.mx0);
175+
flux.mx1 += ds[2] * (flux_R.mx1 - flux_L.mx1);
176+
flux.mx2 += ds[2] * (flux_R.mx2 - flux_L.mx2);
177+
}
178+
179+
EulerCons cons = load<SimdType>(cons_ptrs, base);
180+
cons.d -= dtodv * flux.d;
181+
cons.e -= dtodv * flux.e;
182+
cons.mx0 -= dtodv * flux.mx0;
183+
cons.mx1 -= dtodv * flux.mx1;
184+
cons.mx2 -= dtodv * flux.mx2;
185+
store<SimdType>(cons, cons_ptrs, base);
186+
});
187+
}
188+
189+
template <class T, class IndexType, std::size_t E0, std::size_t E1, std::size_t E2>
190+
void godunov_vec(
191+
Kokkos::DefaultExecutionSpace const& exec_space,
192+
EulerPrimArrays<Kokkos::mdspan<
193+
T const,
194+
Kokkos::extents<IndexType, E0, E1, E2>,
195+
Kokkos::layout_left>> const& prim_arrays,
196+
EulerConsArrays<Kokkos::mdspan<
197+
T,
198+
Kokkos::extents<IndexType, E0, E1, E2>,
199+
Kokkos::layout_left>> const& cons_arrays,
200+
PerfectGas<T> const& eos,
201+
UniformMesh3d<T> const& mesh,
202+
hllc const& riemann_solver,
203+
T const dt)
204+
{
205+
namespace KE = Kokkos::Experimental;
206+
using simd_t = KE::simd<T>;
207+
using simd_scalar_t = KE::basic_simd<T, KE::simd_abi::scalar>;
208+
209+
// interior x-range is [1, nx-1)
210+
IndexType const nx = prim_arrays.d.extent(0);
211+
IndexType const nx_begin = 1;
212+
IndexType const nx_inner = nx - 2; // number of interior cells
213+
IndexType const vec_end = nx_begin + ((nx_inner / simd_t::size()) * simd_t::size());
214+
IndexType const nx_end = nx - 1;
215+
216+
godunov_kernel<simd_t>(
217+
exec_space,
218+
prim_arrays,
219+
cons_arrays,
220+
nx_begin,
221+
vec_end,
222+
eos,
223+
mesh,
224+
riemann_solver,
225+
dt);
226+
227+
if (vec_end < nx_end) {
228+
godunov_kernel<simd_scalar_t>(
229+
exec_space,
230+
prim_arrays,
231+
cons_arrays,
232+
vec_end,
233+
nx_end,
234+
eos,
235+
mesh,
236+
riemann_solver,
237+
dt);
238+
}
239+
}

0 commit comments

Comments
 (0)