Skip to content

Commit de4ab97

Browse files
MarcoArtianoJoshuaLampertranocha
authored
Add GPU kernel for calc_sources! (#3012)
* add source term gpu kernel * add dispatch for no source terms * add 1d and 3d source gpu kernels * specify signature * add tests * Update test/test_amdgpu_2d.jl * Update test/test_amdgpu_2d.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> * refactoring * format * refactor also 3D kernels * format * fix tests * fix tests * fix elixirs * fix tests * fix tests * Apply suggestions from code review Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com> * rename tests, delete 1D P4est * add kernel abstraction tests * add cuda tests * fix tests * fix allocations * fix allocations 3D * fix allocations * Apply suggestions from code review Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> * Update test/test_amdgpu_3d.jl * Apply suggestions from code review Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> * add comments about storage_type and real_type choices Co-authored-by: Marco Artiano <57838732+MarcoArtiano@users.noreply.github.com> * Apply suggestions from code review Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com> --------- Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com>
1 parent 7ff42d4 commit de4ab97

18 files changed

Lines changed: 661 additions & 88 deletions

‎examples/p4est_2d_dgsem/elixir_euler_source_terms.jl‎

Lines changed: 70 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,70 @@
1+
# The same setup as tree_2d_dgsem/elixir_euler_source_terms.jl
2+
# to verify the P4estMesh implementation against TreeMesh
3+
4+
using OrdinaryDiffEqLowStorageRK
5+
using Trixi
6+
7+
###############################################################################
8+
# semidiscretization of the compressible Euler equations
9+
gamma = 1.4
10+
equations = CompressibleEulerEquations2D(gamma)
11+
12+
initial_condition = initial_condition_convergence_test
13+
14+
# Up to version 0.13.0, `max_abs_speed_naive` was used as the default wave speed estimate of
15+
# `const flux_lax_friedrichs = FluxLaxFriedrichs(), i.e., `FluxLaxFriedrichs(max_abs_speed = max_abs_speed_naive)`.
16+
# In the `StepsizeCallback`, though, the less diffusive `max_abs_speeds` is employed which is consistent with `max_abs_speed`.
17+
# Thus, we exchanged in PR#2458 the default wave speed used in the LLF flux to `max_abs_speed`.
18+
# To ensure that every example still runs we specify explicitly `FluxLaxFriedrichs(max_abs_speed_naive)`.
19+
# We remark, however, that the now default `max_abs_speed` is in general recommended due to compliance with the
20+
# `StepsizeCallback` (CFL-Condition) and less diffusion.
21+
solver = DGSEM(polydeg = 3, surface_flux = FluxLaxFriedrichs(max_abs_speed_naive))
22+
23+
coordinates_min = (0.0, 0.0)
24+
coordinates_max = (2.0, 2.0)
25+
26+
trees_per_dimension = (16, 16)
27+
mesh = P4estMesh(trees_per_dimension,
28+
polydeg = 3, initial_refinement_level = 0,
29+
coordinates_min = coordinates_min, coordinates_max = coordinates_max,
30+
periodicity = true)
31+
32+
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver;
33+
source_terms = source_terms_convergence_test,
34+
boundary_conditions = boundary_condition_periodic)
35+
36+
###############################################################################
37+
# ODE solvers, callbacks etc.
38+
39+
tspan = (0.0, 2.0)
40+
# Create ODE problem with time span from 0.0 to 1.0
41+
# Setting `real_type` allows to change the real number type, e.g., to `Float32`.
42+
# This is particularly useful when changing the `storage_type` to a GPU array
43+
# type such as `ROCArray` (AMD) or `CuArray` (NVIDIA CUDA).
44+
ode = semidiscretize(semi, tspan; real_type = nothing, storage_type = nothing)
45+
46+
summary_callback = SummaryCallback()
47+
48+
analysis_interval = 100
49+
analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
50+
51+
alive_callback = AliveCallback(analysis_interval = analysis_interval)
52+
53+
save_solution = SaveSolutionCallback(interval = 100,
54+
save_initial_solution = true,
55+
save_final_solution = true,
56+
solution_variables = cons2prim)
57+
58+
stepsize_callback = StepsizeCallback(cfl = 1.0)
59+
60+
callbacks = CallbackSet(summary_callback,
61+
analysis_callback, alive_callback,
62+
save_solution,
63+
stepsize_callback)
64+
65+
###############################################################################
66+
# run the simulation
67+
68+
sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false);
69+
dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
70+
ode_default_options()..., callback = callbacks);

‎examples/p4est_3d_dgsem/elixir_euler_source_terms.jl‎

Lines changed: 72 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,72 @@
1+
# The same setup as tree_3d_dgsem/elixir_euler_source_terms.jl
2+
# to verify the StructuredMesh implementation against TreeMesh
3+
4+
using OrdinaryDiffEqLowStorageRK
5+
using Trixi
6+
7+
###############################################################################
8+
# semidiscretization of the compressible Euler equations
9+
gamma = 1.4
10+
equations = CompressibleEulerEquations3D(gamma)
11+
12+
initial_condition = initial_condition_convergence_test
13+
14+
# Up to version 0.13.0, `max_abs_speed_naive` was used as the default wave speed estimate of
15+
# `const flux_lax_friedrichs = FluxLaxFriedrichs(), i.e., `FluxLaxFriedrichs(max_abs_speed = max_abs_speed_naive)`.
16+
# In the `StepsizeCallback`, though, the less diffusive `max_abs_speeds` is employed which is consistent with `max_abs_speed`.
17+
# Thus, we exchanged in PR#2458 the default wave speed used in the LLF flux to `max_abs_speed`.
18+
# To ensure that every example still runs we specify explicitly `FluxLaxFriedrichs(max_abs_speed_naive)`.
19+
# We remark, however, that the now default `max_abs_speed` is in general recommended due to compliance with the
20+
# `StepsizeCallback` (CFL-Condition) and less diffusion.
21+
solver = DGSEM(polydeg = 3, surface_flux = FluxLaxFriedrichs(max_abs_speed_naive),
22+
volume_integral = VolumeIntegralWeakForm())
23+
24+
coordinates_min = (0.0, 0.0, 0.0)
25+
coordinates_max = (2.0, 2.0, 2.0)
26+
27+
trees_per_dimension = (4, 4, 4)
28+
29+
mesh = P4estMesh(trees_per_dimension, polydeg = 3,
30+
coordinates_min = coordinates_min, coordinates_max = coordinates_max,
31+
initial_refinement_level = 1,
32+
periodicity = true)
33+
34+
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver;
35+
source_terms = source_terms_convergence_test,
36+
boundary_conditions = boundary_condition_periodic)
37+
38+
###############################################################################
39+
# ODE solvers, callbacks etc.
40+
41+
tspan = (0.0, 5.0)
42+
# Create ODE problem with time span from 0.0 to 1.0
43+
# Setting `real_type` allows to change the real number type, e.g., to `Float32`.
44+
# This is particularly useful when changing the `storage_type` to a GPU array
45+
# type such as `ROCArray` (AMD) or `CuArray` (NVIDIA CUDA).
46+
ode = semidiscretize(semi, tspan; real_type = nothing, storage_type = nothing)
47+
48+
summary_callback = SummaryCallback()
49+
50+
analysis_interval = 100
51+
analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
52+
53+
alive_callback = AliveCallback(analysis_interval = analysis_interval)
54+
55+
save_solution = SaveSolutionCallback(interval = 100,
56+
save_initial_solution = true,
57+
save_final_solution = true,
58+
solution_variables = cons2prim)
59+
60+
stepsize_callback = StepsizeCallback(cfl = 0.6)
61+
62+
callbacks = CallbackSet(summary_callback,
63+
analysis_callback, alive_callback,
64+
save_solution,
65+
stepsize_callback)
66+
67+
###############################################################################
68+
# run the simulation
69+
70+
sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false);
71+
dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
72+
ode_default_options()..., callback = callbacks);

‎src/solvers/dgsem_p4est/dg.jl‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -93,6 +93,7 @@ include("dg_2d_gpu.jl")
9393
include("dg_3d.jl")
9494
include("dg_3d_parabolic.jl")
9595
include("dg_parallel.jl")
96+
include("dg_3d_gpu.jl")
9697

9798
# Subcell limiters
9899
include("subcell_limiters.jl")

‎src/solvers/dgsem_p4est/dg_2d.jl‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1245,7 +1245,7 @@ function rhs!(du, u, t, u_parent, semis,
12451245

12461246
# Calculate source terms
12471247
@trixi_timeit timer() "source terms" begin
1248-
calc_sources!(du, u, t, source_terms, equations, dg, cache)
1248+
calc_sources!(backend, du, u, t, source_terms, equations, dg, cache)
12491249
end
12501250

12511251
return nothing

‎src/solvers/dgsem_p4est/dg_2d_gpu.jl‎

Lines changed: 29 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -118,4 +118,33 @@ end
118118
apply_jacobian_per_quadrature_node!(du, MeshT, equations, dg, inverse_jacobian,
119119
i, j, element)
120120
end
121+
122+
@kernel function calc_sources_KAkernel!(du, u, t, source_terms,
123+
node_coordinates,
124+
equations::AbstractEquations{2}, dg, cache)
125+
i, j, element = @index(Global, NTuple)
126+
u_local = get_node_vars(u, equations, dg, i, j, element)
127+
x_local = get_node_coords(node_coordinates, equations, dg, i, j, element)
128+
129+
du_local = source_terms(u_local, x_local, t, equations)
130+
131+
add_to_node_vars!(du, du_local, equations, dg, i, j, element)
132+
end
133+
134+
function calc_sources!(backend::Backend, du, u, t, source_terms,
135+
equations::AbstractEquations{2}, dg::DG, cache)
136+
nelements(dg, cache) == 0 && return nothing
137+
@unpack node_coordinates = cache.elements
138+
kernel_cache = kernel_filter_cache(cache)
139+
kernel! = calc_sources_KAkernel!(backend)
140+
kernel!(du, u, t, source_terms, node_coordinates, equations, dg, kernel_cache,
141+
ndrange = (nnodes(dg), nnodes(dg), nelements(dg, cache)))
142+
143+
return nothing
144+
end
145+
146+
function calc_sources!(backend::Backend, du, u, t, source_terms::Nothing,
147+
equations::AbstractEquations{2}, dg::DG, cache)
148+
return nothing
121149
end
150+
end #muladd

‎src/solvers/dgsem_p4est/dg_3d.jl‎

Lines changed: 0 additions & 76 deletions
Original file line numberDiff line numberDiff line change
@@ -106,27 +106,6 @@ function prolong2interfaces!(backend::Nothing, cache, u,
106106
return nothing
107107
end
108108

109-
function prolong2interfaces!(backend::Backend, cache, u,
110-
mesh::Union{P4estMesh{3}, T8codeMesh{3}},
111-
equations, dg::DG)
112-
@unpack interfaces = cache
113-
@unpack neighbor_ids, node_indices = cache.interfaces
114-
index_range = eachnode(dg)
115-
116-
kernel! = prolong2interfaces_KAkernel!(backend)
117-
kernel!(interfaces.u, u, typeof(mesh), equations, neighbor_ids, node_indices,
118-
index_range,
119-
ndrange = ninterfaces(interfaces))
120-
return nothing
121-
end
122-
123-
@kernel function prolong2interfaces_KAkernel!(interface_u, u, MeshT, equations,
124-
neighbor_ids, node_indices, index_range)
125-
interface = @index(Global)
126-
prolong2interfaces_per_interface!(interface_u, u, MeshT, equations, neighbor_ids,
127-
node_indices, index_range, interface)
128-
end
129-
130109
@inline function prolong2interfaces_per_interface!(u_interface, u,
131110
::Type{<:Union{P4estMesh{3},
132111
T8codeMesh{3}}},
@@ -223,38 +202,6 @@ function calc_interface_flux!(backend::Nothing, surface_flux_values,
223202
return nothing
224203
end
225204

226-
function calc_interface_flux!(backend::Backend, surface_flux_values,
227-
mesh::Union{P4estMesh{3}, T8codeMesh{3}},
228-
have_nonconservative_terms,
229-
equations, surface_integral, dg::DG, cache)
230-
@unpack neighbor_ids, node_indices = cache.interfaces
231-
@unpack contravariant_vectors = cache.elements
232-
index_range = eachnode(dg)
233-
234-
kernel! = calc_interface_flux_KAkernel!(backend)
235-
kernel!(surface_flux_values, typeof(mesh), have_nonconservative_terms, equations,
236-
surface_integral, typeof(dg), cache.interfaces.u,
237-
neighbor_ids, node_indices, contravariant_vectors, index_range,
238-
ndrange = ninterfaces(cache.interfaces))
239-
return nothing
240-
end
241-
242-
@kernel function calc_interface_flux_KAkernel!(surface_flux_values, MeshT,
243-
have_nonconservative_terms, equations,
244-
surface_integral, SolverT, u_interface,
245-
neighbor_ids, node_indices,
246-
contravariant_vectors, index_range)
247-
interface = @index(Global)
248-
calc_interface_flux_per_interface!(surface_flux_values,
249-
MeshT,
250-
have_nonconservative_terms,
251-
equations, surface_integral, SolverT,
252-
u_interface,
253-
neighbor_ids, node_indices,
254-
contravariant_vectors,
255-
index_range, interface)
256-
end
257-
258205
@inline function calc_interface_flux_per_interface!(surface_flux_values,
259206
MeshT::Type{<:Union{P4estMesh{3},
260207
T8codeMesh{3}}},
@@ -1026,29 +973,6 @@ function calc_surface_integral!(backend::Nothing, du, u,
1026973
return nothing
1027974
end
1028975

1029-
function calc_surface_integral!(backend::Backend, du, u,
1030-
mesh::Union{P4estMesh{3}, T8codeMesh{3}},
1031-
equations,
1032-
surface_integral::SurfaceIntegralWeakForm,
1033-
dg::DGSEM, cache)
1034-
@unpack inverse_weights = dg.basis
1035-
@unpack surface_flux_values = cache.elements
1036-
1037-
kernel! = calc_surface_integral_KAkernel!(backend)
1038-
kernel!(du, typeof(mesh), equations, surface_integral, dg, inverse_weights[1],
1039-
surface_flux_values, ndrange = nelements(cache.elements))
1040-
return nothing
1041-
end
1042-
1043-
@kernel function calc_surface_integral_KAkernel!(du, MeshT, equations,
1044-
surface_integral, dg, factor,
1045-
surface_flux_values)
1046-
element = @index(Global)
1047-
calc_surface_integral_per_element!(du, MeshT,
1048-
equations, surface_integral, dg, factor,
1049-
surface_flux_values, element)
1050-
end
1051-
1052976
@inline function calc_surface_integral_per_element!(du,
1053977
::Type{<:Union{P4estMesh{3},
1054978
T8codeMesh{3}}},

0 commit comments

Comments
 (0)