From d50532f4a87af72ea976cbb62cd56901607f1883 Mon Sep 17 00:00:00 2001
From: Marco Artiano <artianomarco97@gmail.com>
Date: Sun, 17 May 2026 22:33:31 +0200
Subject: [PATCH 01/29] add source term gpu kernel

---
 src/solvers/dgsem_tree/dg_2d.jl | 29 ++++++++++++++++++++++++++---
 1 file changed, 26 insertions(+), 3 deletions(-)

diff --git a/src/solvers/dgsem_tree/dg_2d.jl b/src/solvers/dgsem_tree/dg_2d.jl
index 815e0c0df13..74b93dd3f76 100644
--- a/src/solvers/dgsem_tree/dg_2d.jl
+++ b/src/solvers/dgsem_tree/dg_2d.jl
@@ -178,7 +178,7 @@ function rhs!(du, u, t,
 
     # Calculate source terms
     @trixi_timeit_ext backend timer() "source terms" begin
-        calc_sources!(du, u, t, source_terms, equations, dg, cache)
+        calc_sources!(backend, du, u, t, source_terms, equations, dg, cache)
     end
 
     return nothing
@@ -1359,12 +1359,12 @@ function apply_jacobian!(backend::Nothing, du, mesh::TreeMesh{2},
 end
 
 # Need dimension specific version to avoid error at dispatching
-function calc_sources!(du, u, t, source_terms::Nothing,
+function calc_sources!(backend::Nothing, du, u, t, source_terms::Nothing,
                        equations::AbstractEquations{2}, dg::DG, cache)
     return nothing
 end
 
-function calc_sources!(du, u, t, source_terms,
+function calc_sources!(backend::Nothing, du, u, t, source_terms,
                        equations::AbstractEquations{2}, dg::DG, cache)
     @unpack node_coordinates = cache.elements
 
@@ -1380,4 +1380,27 @@ function calc_sources!(du, u, t, source_terms,
 
     return nothing
 end
+
+@kernel function calc_sources_KAkernel!(du, u, t, source_terms,
+                                        node_coordinates, equations, dg, cache)
+    i, j, element = @index(Global, NTuple)
+    u_local = get_node_vars(u, equations, dg, i, j, element)
+    x_local = get_node_coords(node_coordinates, equations, dg, i, j, element)
+
+    du_local = source_terms(u_local, x_local, t, equations)
+
+    add_to_node_vars!(du, du_local, equations, dg, i, j, element)
+end
+
+function calc_sources!(backend::Backend, du, u, t, source_terms,
+                       equations::AbstractEquations{2}, dg::DG, cache)
+    nelements(dg, cache) == 0 && return nothing
+    @unpack node_coordinates = cache.elements
+    kernel_cache = kernel_filter_cache(cache)
+    kernel! = calc_sources_KAkernel!(backend)
+    kernel!(du, u, t, source_terms, node_coordinates, equations, dg, kernel_cache,
+            ndrange = (nnodes(dg), nnodes(dg), nelements(dg, cache)))
+
+    return nothing
+end
 end # @muladd

From 453f04dc12da94e7845dd61831b633b02d7aa255 Mon Sep 17 00:00:00 2001
From: Marco Artiano <artianomarco97@gmail.com>
Date: Sun, 17 May 2026 22:46:07 +0200
Subject: [PATCH 02/29] add dispatch for no source terms

---
 src/solvers/dgsem_tree/dg_2d.jl | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/src/solvers/dgsem_tree/dg_2d.jl b/src/solvers/dgsem_tree/dg_2d.jl
index 74b93dd3f76..7b66250abdc 100644
--- a/src/solvers/dgsem_tree/dg_2d.jl
+++ b/src/solvers/dgsem_tree/dg_2d.jl
@@ -1403,4 +1403,9 @@ function calc_sources!(backend::Backend, du, u, t, source_terms,
 
     return nothing
 end
+
+function calc_sources!(backend::Backend, du, u, t, source_terms::Nothing,
+                       equations::AbstractEquations{2}, dg::DG, cache)
+    return nothing
+end
 end # @muladd

From 5b7e0e8014789c9f2fd920385f4394492ccd79f4 Mon Sep 17 00:00:00 2001
From: Marco Artiano <artianomarco97@gmail.com>
Date: Sun, 17 May 2026 23:25:12 +0200
Subject: [PATCH 03/29] add 1d and 3d source gpu kernels

---
 src/solvers/dgsem_p4est/dg_2d.jl          |  2 +-
 src/solvers/dgsem_p4est/dg_3d_parallel.jl |  2 +-
 src/solvers/dgsem_structured/dg.jl        |  2 +-
 src/solvers/dgsem_tree/dg_1d.jl           | 35 +++++++++++++++++++++--
 src/solvers/dgsem_tree/dg_2d_parallel.jl  |  2 +-
 src/solvers/dgsem_tree/dg_3d.jl           | 33 +++++++++++++++++++--
 6 files changed, 67 insertions(+), 9 deletions(-)

diff --git a/src/solvers/dgsem_p4est/dg_2d.jl b/src/solvers/dgsem_p4est/dg_2d.jl
index f60ece112c1..77eec3a57e1 100644
--- a/src/solvers/dgsem_p4est/dg_2d.jl
+++ b/src/solvers/dgsem_p4est/dg_2d.jl
@@ -1338,7 +1338,7 @@ function rhs!(du, u, t, u_parent, semis,
 
     # Calculate source terms
     @trixi_timeit timer() "source terms" begin
-        calc_sources!(du, u, t, source_terms, equations, dg, cache)
+        calc_sources!(backend, du, u, t, source_terms, equations, dg, cache)
     end
 
     return nothing
diff --git a/src/solvers/dgsem_p4est/dg_3d_parallel.jl b/src/solvers/dgsem_p4est/dg_3d_parallel.jl
index 91a700a847c..6fc4e2009b0 100644
--- a/src/solvers/dgsem_p4est/dg_3d_parallel.jl
+++ b/src/solvers/dgsem_p4est/dg_3d_parallel.jl
@@ -107,7 +107,7 @@ function rhs!(du, u, t,
 
     # Calculate source terms
     @trixi_timeit timer() "source terms" begin
-        calc_sources!(du, u, t, source_terms, equations, dg, cache)
+        calc_sources!(backend, du, u, t, source_terms, equations, dg, cache)
     end
 
     # Finish to send MPI data
diff --git a/src/solvers/dgsem_structured/dg.jl b/src/solvers/dgsem_structured/dg.jl
index f3e8365b6f6..5ac27049151 100644
--- a/src/solvers/dgsem_structured/dg.jl
+++ b/src/solvers/dgsem_structured/dg.jl
@@ -88,7 +88,7 @@ function rhs!(du, u, t,
 
     # Calculate source terms
     @trixi_timeit timer() "source terms" begin
-        calc_sources!(du, u, t, source_terms, equations, dg, cache)
+        calc_sources!(backend, du, u, t, source_terms, equations, dg, cache)
     end
 
     return nothing
diff --git a/src/solvers/dgsem_tree/dg_1d.jl b/src/solvers/dgsem_tree/dg_1d.jl
index 65383168835..f7f5f88ccd0 100644
--- a/src/solvers/dgsem_tree/dg_1d.jl
+++ b/src/solvers/dgsem_tree/dg_1d.jl
@@ -111,7 +111,7 @@ function rhs!(du, u, t,
 
     # Calculate source terms
     @trixi_timeit timer() "source terms" begin
-        calc_sources!(du, u, t, source_terms, equations, dg, cache)
+        calc_sources!(backend, du, u, t, source_terms, equations, dg, cache)
     end
 
     return nothing
@@ -765,12 +765,12 @@ function apply_jacobian!(backend::Nothing, du, mesh::TreeMesh{1},
 end
 
 # Need dimension specific version to avoid error at dispatching
-function calc_sources!(du, u, t, source_terms::Nothing,
+function calc_sources!(backend::Nothing, du, u, t, source_terms::Nothing,
                        equations::AbstractEquations{1}, dg::DG, cache)
     return nothing
 end
 
-function calc_sources!(du, u, t, source_terms,
+function calc_sources!(backend::Nothing, du, u, t, source_terms,
                        equations::AbstractEquations{1}, dg::DG, cache)
     @unpack node_coordinates = cache.elements
 
@@ -786,4 +786,33 @@ function calc_sources!(du, u, t, source_terms,
 
     return nothing
 end
+
+@kernel function calc_sources_KAkernel!(du, u, t, source_terms,
+                                        node_coordinates,
+                                        equations::AbstractEquations{1}, dg, cache)
+    i, element = @index(Global, NTuple)
+    u_local = get_node_vars(u, equations, dg, i, element)
+    x_local = get_node_coords(node_coordinates, equations, dg, i, element)
+
+    du_local = source_terms(u_local, x_local, t, equations)
+
+    add_to_node_vars!(du, du_local, equations, dg, i, element)
+end
+
+function calc_sources!(backend::Backend, du, u, t, source_terms,
+                       equations::AbstractEquations{1}, dg::DG, cache)
+    nelements(dg, cache) == 0 && return nothing
+    @unpack node_coordinates = cache.elements
+    kernel_cache = kernel_filter_cache(cache)
+    kernel! = calc_sources_KAkernel!(backend)
+    kernel!(du, u, t, source_terms, node_coordinates, equations, dg, kernel_cache,
+            ndrange = (nnodes(dg), nelements(dg, cache)))
+
+    return nothing
+end
+
+function calc_sources!(backend::Backend, du, u, t, source_terms::Nothing,
+                       equations::AbstractEquations{1}, dg::DG, cache)
+    return nothing
+end
 end # @muladd
diff --git a/src/solvers/dgsem_tree/dg_2d_parallel.jl b/src/solvers/dgsem_tree/dg_2d_parallel.jl
index d4b6192a4df..822ee9b7d26 100644
--- a/src/solvers/dgsem_tree/dg_2d_parallel.jl
+++ b/src/solvers/dgsem_tree/dg_2d_parallel.jl
@@ -554,7 +554,7 @@ function rhs!(du, u, t,
 
     # Calculate source terms
     @trixi_timeit timer() "source terms" begin
-        calc_sources!(du, u, t, source_terms, equations, dg, cache)
+        calc_sources!(backend, du, u, t, source_terms, equations, dg, cache)
     end
 
     # Finish to send MPI data
diff --git a/src/solvers/dgsem_tree/dg_3d.jl b/src/solvers/dgsem_tree/dg_3d.jl
index 3fe1bc8b0a3..06632a77edc 100644
--- a/src/solvers/dgsem_tree/dg_3d.jl
+++ b/src/solvers/dgsem_tree/dg_3d.jl
@@ -1422,12 +1422,12 @@ function apply_jacobian!(backend::Nothing, du, mesh::TreeMesh{3},
 end
 
 # Need dimension specific version to avoid error at dispatching
-function calc_sources!(du, u, t, source_terms::Nothing,
+function calc_sources!(backend::Nothing, du, u, t, source_terms::Nothing,
                        equations::AbstractEquations{3}, dg::DG, cache)
     return nothing
 end
 
-function calc_sources!(du, u, t, source_terms,
+function calc_sources!(backend::Nothing, du, u, t, source_terms,
                        equations::AbstractEquations{3}, dg::DG, cache)
     @unpack node_coordinates = cache.elements
 
@@ -1443,4 +1443,33 @@ function calc_sources!(du, u, t, source_terms,
 
     return nothing
 end
+
+@kernel function calc_sources_KAkernel!(du, u, t, source_terms,
+                                        node_coordinates,
+                                        equations::AbstractEquations{3}, dg, cache)
+    i, j, k, element = @index(Global, NTuple)
+    u_local = get_node_vars(u, equations, dg, i, j, k, element)
+    x_local = get_node_coords(node_coordinates, equations, dg, i, j, k, element)
+
+    du_local = source_terms(u_local, x_local, t, equations)
+
+    add_to_node_vars!(du, du_local, equations, dg, i, j, k, element)
+end
+
+function calc_sources!(backend::Backend, du, u, t, source_terms,
+                       equations::AbstractEquations{3}, dg::DG, cache)
+    nelements(dg, cache) == 0 && return nothing
+    @unpack node_coordinates = cache.elements
+    kernel_cache = kernel_filter_cache(cache)
+    kernel! = calc_sources_KAkernel!(backend)
+    kernel!(du, u, t, source_terms, node_coordinates, equations, dg, kernel_cache,
+            ndrange = (nnodes(dg), nnodes(dg), nnodes(dg), nelements(dg, cache)))
+
+    return nothing
+end
+
+function calc_sources!(backend::Backend, du, u, t, source_terms::Nothing,
+                       equations::AbstractEquations{3}, dg::DG, cache)
+    return nothing
+end
 end # @muladd

From 158e4818b04b3edb0d41b48255fd1f248d69435c Mon Sep 17 00:00:00 2001
From: Marco Artiano <artianomarco97@gmail.com>
Date: Sun, 17 May 2026 23:27:37 +0200
Subject: [PATCH 04/29] specify signature

---
 src/solvers/dgsem_tree/dg_2d.jl | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/solvers/dgsem_tree/dg_2d.jl b/src/solvers/dgsem_tree/dg_2d.jl
index 7b66250abdc..373938976bd 100644
--- a/src/solvers/dgsem_tree/dg_2d.jl
+++ b/src/solvers/dgsem_tree/dg_2d.jl
@@ -1382,7 +1382,8 @@ function calc_sources!(backend::Nothing, du, u, t, source_terms,
 end
 
 @kernel function calc_sources_KAkernel!(du, u, t, source_terms,
-                                        node_coordinates, equations, dg, cache)
+                                        node_coordinates,
+                                        equations::AbstractEquations{2}, dg, cache)
     i, j, element = @index(Global, NTuple)
     u_local = get_node_vars(u, equations, dg, i, j, element)
     x_local = get_node_coords(node_coordinates, equations, dg, i, j, element)

From 8f8b0efe1d0c67a193a27b131927b191104e558c Mon Sep 17 00:00:00 2001
From: Marco Artiano <artianomarco97@gmail.com>
Date: Mon, 18 May 2026 12:03:37 +0200
Subject: [PATCH 05/29] add tests

---
 .../elixir_euler_source_terms_gpu.jl          | 66 ++++++++++++++++++
 .../elixir_euler_source_terms_gpu.jl          | 68 +++++++++++++++++++
 test/test_amdgpu_2d.jl                        | 66 ++++++++++++++++++
 test/test_amdgpu_3d.jl                        | 67 ++++++++++++++++++
 4 files changed, 267 insertions(+)
 create mode 100644 examples/p4est_2d_dgsem/elixir_euler_source_terms_gpu.jl
 create mode 100644 examples/p4est_3d_dgsem/elixir_euler_source_terms_gpu.jl

diff --git a/examples/p4est_2d_dgsem/elixir_euler_source_terms_gpu.jl b/examples/p4est_2d_dgsem/elixir_euler_source_terms_gpu.jl
new file mode 100644
index 00000000000..755227e5985
--- /dev/null
+++ b/examples/p4est_2d_dgsem/elixir_euler_source_terms_gpu.jl
@@ -0,0 +1,66 @@
+# The same setup as tree_2d_dgsem/elixir_euler_source_terms.jl
+# to verify the StructuredMesh implementation against TreeMesh
+
+using OrdinaryDiffEqLowStorageRK
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+gamma = 1.4
+equations = CompressibleEulerEquations2D(gamma)
+
+initial_condition = initial_condition_convergence_test
+
+# Up to version 0.13.0, `max_abs_speed_naive` was used as the default wave speed estimate of
+# `const flux_lax_friedrichs = FluxLaxFriedrichs(), i.e., `FluxLaxFriedrichs(max_abs_speed = max_abs_speed_naive)`.
+# In the `StepsizeCallback`, though, the less diffusive `max_abs_speeds` is employed which is consistent with `max_abs_speed`.
+# Thus, we exchanged in PR#2458 the default wave speed used in the LLF flux to `max_abs_speed`.
+# To ensure that every example still runs we specify explicitly `FluxLaxFriedrichs(max_abs_speed_naive)`.
+# We remark, however, that the now default `max_abs_speed` is in general recommended due to compliance with the
+# `StepsizeCallback` (CFL-Condition) and less diffusion.
+solver = DGSEM(polydeg = 3, surface_flux = FluxLaxFriedrichs(max_abs_speed_naive))
+
+coordinates_min = (0.0, 0.0)
+coordinates_max = (2.0, 2.0)
+
+trees_per_dimension = (16, 16)
+mesh = P4estMesh(trees_per_dimension,
+                 polydeg = 3, initial_refinement_level = 0,
+                 coordinates_min = coordinates_min, coordinates_max = coordinates_max,
+                 periodicity = true)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver;
+                                    source_terms = source_terms_convergence_test,
+                                    boundary_conditions = boundary_condition_periodic)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 2.0)
+ode = semidiscretize(semi, tspan; real_type = real_type, storage_type = storage_type)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(interval = 100,
+                                     save_initial_solution = true,
+                                     save_final_solution = true,
+                                     solution_variables = cons2prim)
+
+stepsize_callback = StepsizeCallback(cfl = 1.0)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback,
+                        save_solution,
+                        stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false);
+            dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            ode_default_options()..., callback = callbacks);
diff --git a/examples/p4est_3d_dgsem/elixir_euler_source_terms_gpu.jl b/examples/p4est_3d_dgsem/elixir_euler_source_terms_gpu.jl
new file mode 100644
index 00000000000..6dda9dcc8fb
--- /dev/null
+++ b/examples/p4est_3d_dgsem/elixir_euler_source_terms_gpu.jl
@@ -0,0 +1,68 @@
+# The same setup as tree_3d_dgsem/elixir_euler_source_terms.jl
+# to verify the StructuredMesh implementation against TreeMesh
+
+using OrdinaryDiffEqLowStorageRK
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+gamma = 1.4
+equations = CompressibleEulerEquations3D(gamma)
+
+initial_condition = initial_condition_convergence_test
+
+# Up to version 0.13.0, `max_abs_speed_naive` was used as the default wave speed estimate of
+# `const flux_lax_friedrichs = FluxLaxFriedrichs(), i.e., `FluxLaxFriedrichs(max_abs_speed = max_abs_speed_naive)`.
+# In the `StepsizeCallback`, though, the less diffusive `max_abs_speeds` is employed which is consistent with `max_abs_speed`.
+# Thus, we exchanged in PR#2458 the default wave speed used in the LLF flux to `max_abs_speed`.
+# To ensure that every example still runs we specify explicitly `FluxLaxFriedrichs(max_abs_speed_naive)`.
+# We remark, however, that the now default `max_abs_speed` is in general recommended due to compliance with the
+# `StepsizeCallback` (CFL-Condition) and less diffusion.
+solver = DGSEM(polydeg = 3, surface_flux = FluxLaxFriedrichs(max_abs_speed_naive),
+               volume_integral = VolumeIntegralWeakForm())
+
+coordinates_min = (0.0, 0.0, 0.0)
+coordinates_max = (2.0, 2.0, 2.0)
+
+trees_per_dimension = (4, 4, 4)
+
+mesh = P4estMesh(trees_per_dimension, polydeg = 3,
+                 coordinates_min = coordinates_min, coordinates_max = coordinates_max,
+                 initial_refinement_level = 1,
+                 periodicity = true)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver;
+                                    source_terms = source_terms_convergence_test,
+                                    boundary_conditions = boundary_condition_periodic)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 5.0)
+ode = semidiscretize(semi, tspan; real_type = real_type, storage_type = storage_type)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(interval = 100,
+                                     save_initial_solution = true,
+                                     save_final_solution = true,
+                                     solution_variables = cons2prim)
+
+stepsize_callback = StepsizeCallback(cfl = 0.6)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback,
+                        save_solution,
+                        stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false);
+            dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            ode_default_options()..., callback = callbacks);
diff --git a/test/test_amdgpu_2d.jl b/test/test_amdgpu_2d.jl
index 13f98e1b958..0512f217ea9 100644
--- a/test/test_amdgpu_2d.jl
+++ b/test/test_amdgpu_2d.jl
@@ -67,6 +67,72 @@ end
     @test Trixi.storage_type(ode.p.cache.mortars) === ROCArray
 end
 
+@trixi_testset "elixir_euler_source_terms_gpu.jl native" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_source_terms_gpu.jl"),
+                        # Expected errors are exactly the same as with TreeMesh!
+                        l2=[
+                            4.893619139889976e-5,
+                            5.3526950567182756e-5,
+                            5.35269505672133e-5,
+                            5.352695056735998e-5,
+                            0.00015172095200428318
+                        ],
+                        linf=[
+                            0.00031179856625374036,
+                            0.0003368725355339386,
+                            0.0003368725355383795,
+                            0.00033687253560787944,
+                            0.0013193387520935573
+                        ])
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
+    semi = ode.p # `semidiscretize` adapts the semi, so we need to obtain it from the ODE problem.
+    @test_allocations(Trixi.rhs!, semi, sol, 1000)
+    @test real(ode.p.solver) == Float64
+    @test real(ode.p.solver.basis) == Float64
+    @test real(ode.p.solver.mortar) == Float64
+    # TODO: `mesh` is currently not `adapt`ed correctly
+    @test real(ode.p.mesh) == Float64
+
+    @test ode.u0 isa Array
+    @test ode.p.solver.basis.derivative_matrix isa Array
+
+    @test Trixi.storage_type(ode.p.cache.elements) === Array
+    @test Trixi.storage_type(ode.p.cache.interfaces) === Array
+    @test Trixi.storage_type(ode.p.cache.boundaries) === Array
+    @test Trixi.storage_type(ode.p.cache.mortars) === Array
+end
+
+@trixi_testset "elixir_source_terms_gpu.jl Float32 / AMDGPU" begin
+    # Using AMDGPU inside the testset since otherwise the bindings are hiddend by the anonymous modules
+    using AMDGPU
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_source_terms_gpu.jl"),
+                        # Expected errors are exactly the same as with TreeMesh!
+                        l2=Float32[2.4917018095933837e-6, 2.7148269885239423e-6,
+                                   2.695290306860358e-6, 6.243861976167833e-6],
+                        linf=Float32[1.6489475493930428e-5, 1.7499923706143505e-5, 1.893043518075288e-5, 6.214141845717336e-5]
+                        RealT_for_test_tolerances = Float32,
+                        real_type=Float32,
+                        storage_type=ROCArray, gamma=Float32(1.4))
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
+    semi = ode.p # `semidiscretize` adapts the semi, so we need to obtain it from the ODE problem.
+    @test_allocations(Trixi.rhs!, semi, sol, 100_000)
+    @test real(ode.p.solver) == Float32
+    @test real(ode.p.solver.basis) == Float32
+    @test real(ode.p.solver.mortar) == Float32
+    # TODO: `mesh` is currently not `adapt`ed correctly
+    @test real(ode.p.mesh) == Float64
+
+    @test ode.u0 isa ROCArray
+    @test ode.p.solver.basis.derivative_matrix isa ROCArray
+
+    @test Trixi.storage_type(ode.p.cache.elements) === ROCArray
+    @test Trixi.storage_type(ode.p.cache.interfaces) === ROCArray
+    @test Trixi.storage_type(ode.p.cache.boundaries) === ROCArray
+    @test Trixi.storage_type(ode.p.cache.mortars) === ROCArray
+end
+
 # Clean up afterwards: delete Trixi.jl output directory
 @test_nowarn isdir(outdir) && rm(outdir, recursive = true)
 end
diff --git a/test/test_amdgpu_3d.jl b/test/test_amdgpu_3d.jl
index b70ca55b962..f0548b056fc 100644
--- a/test/test_amdgpu_3d.jl
+++ b/test/test_amdgpu_3d.jl
@@ -67,6 +67,73 @@ end
     @test Trixi.storage_type(ode.p.cache.mortars) === ROCArray
 end
 
+@trixi_testset "elixir_euler_source_terms_gpu.jl native" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_source_terms_gpu.jl"),
+                        # Expected errors are exactly the same as with TreeMesh!
+                        l2=[
+                            4.893619139889976e-5,
+                            5.3526950567182756e-5,
+                            5.35269505672133e-5,
+                            5.352695056735998e-5,
+                            0.00015172095200428318
+                        ],
+                        linf=[
+                            0.00031179856625374036,
+                            0.0003368725355339386,
+                            0.0003368725355383795,
+                            0.00033687253560787944,
+                            0.0013193387520935573
+                        ])
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
+    semi = ode.p # `semidiscretize` adapts the semi, so we need to obtain it from the ODE problem.
+    @test_allocations(Trixi.rhs!, semi, sol, 1000)
+    @test real(ode.p.solver) == Float64
+    @test real(ode.p.solver.basis) == Float64
+    @test real(ode.p.solver.mortar) == Float64
+    # TODO: `mesh` is currently not `adapt`ed correctly
+    @test real(ode.p.mesh) == Float64
+
+    @test ode.u0 isa Array
+    @test ode.p.solver.basis.derivative_matrix isa Array
+
+    @test Trixi.storage_type(ode.p.cache.elements) === Array
+    @test Trixi.storage_type(ode.p.cache.interfaces) === Array
+    @test Trixi.storage_type(ode.p.cache.boundaries) === Array
+    @test Trixi.storage_type(ode.p.cache.mortars) === Array
+end
+
+@trixi_testset "elixir_source_terms_gpu.jl Float32 / AMDGPU" begin
+    # Using AMDGPU inside the testset since otherwise the bindings are hiddend by the anonymous modules
+    using AMDGPU
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_source_terms_gpu.jl"),
+                        # Expected errors are exactly the same as with TreeMesh!
+                        l2=Float32[4.912578089985958e-5, 5.3683407014580115e-5,
+                                   5.368099834769191e-5, 5.371664525206341e-5,
+                                   0.00015186256300882088],
+                        linf=Float32[0.00032772542853032327, 0.00035144807715092874, 0.0003549051465479014, 0.00035573961157475686, 0.0013591384887696734]
+                        RealT_for_test_tolerances = Float32,
+                        real_type=Float32,
+                        storage_type=ROCArray, gamma=Float32(1.4))
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
+    semi = ode.p # `semidiscretize` adapts the semi, so we need to obtain it from the ODE problem.
+    @test_allocations(Trixi.rhs!, semi, sol, 100_000)
+    @test real(ode.p.solver) == Float32
+    @test real(ode.p.solver.basis) == Float32
+    @test real(ode.p.solver.mortar) == Float32
+    # TODO: `mesh` is currently not `adapt`ed correctly
+    @test real(ode.p.mesh) == Float64
+
+    @test ode.u0 isa ROCArray
+    @test ode.p.solver.basis.derivative_matrix isa ROCArray
+
+    @test Trixi.storage_type(ode.p.cache.elements) === ROCArray
+    @test Trixi.storage_type(ode.p.cache.interfaces) === ROCArray
+    @test Trixi.storage_type(ode.p.cache.boundaries) === ROCArray
+    @test Trixi.storage_type(ode.p.cache.mortars) === ROCArray
+end
+
 # Clean up afterwards: delete Trixi.jl output directory
 @test_nowarn isdir(outdir) && rm(outdir, recursive = true)
 end

From 5add2610890427de5c3153832f3783c5e778328d Mon Sep 17 00:00:00 2001
From: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com>
Date: Mon, 18 May 2026 17:17:30 +0200
Subject: [PATCH 06/29] Update test/test_amdgpu_2d.jl

---
 test/test_amdgpu_2d.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/test_amdgpu_2d.jl b/test/test_amdgpu_2d.jl
index 0512f217ea9..45bf641b4e9 100644
--- a/test/test_amdgpu_2d.jl
+++ b/test/test_amdgpu_2d.jl
@@ -110,7 +110,7 @@ end
                         # Expected errors are exactly the same as with TreeMesh!
                         l2=Float32[2.4917018095933837e-6, 2.7148269885239423e-6,
                                    2.695290306860358e-6, 6.243861976167833e-6],
-                        linf=Float32[1.6489475493930428e-5, 1.7499923706143505e-5, 1.893043518075288e-5, 6.214141845717336e-5]
+                        linf=Float32[1.6489475493930428e-5, 1.7499923706143505e-5, 1.893043518075288e-5, 6.214141845717336e-5],
                         RealT_for_test_tolerances = Float32,
                         real_type=Float32,
                         storage_type=ROCArray, gamma=Float32(1.4))

From 0a5599208ddeedad9bc18d136594a7b3a92a4f40 Mon Sep 17 00:00:00 2001
From: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com>
Date: Mon, 18 May 2026 17:28:18 +0200
Subject: [PATCH 07/29] Update test/test_amdgpu_2d.jl

Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com>
---
 test/test_amdgpu_2d.jl | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/test/test_amdgpu_2d.jl b/test/test_amdgpu_2d.jl
index 45bf641b4e9..2d641906cd9 100644
--- a/test/test_amdgpu_2d.jl
+++ b/test/test_amdgpu_2d.jl
@@ -110,8 +110,9 @@ end
                         # Expected errors are exactly the same as with TreeMesh!
                         l2=Float32[2.4917018095933837e-6, 2.7148269885239423e-6,
                                    2.695290306860358e-6, 6.243861976167833e-6],
-                        linf=Float32[1.6489475493930428e-5, 1.7499923706143505e-5, 1.893043518075288e-5, 6.214141845717336e-5],
-                        RealT_for_test_tolerances = Float32,
+                        linf=Float32[1.6489475493930428e-5, 1.7499923706143505e-5,
+                                     1.893043518075288e-5, 6.214141845717336e-5],
+                        RealT_for_test_tolerances=Float32,
                         real_type=Float32,
                         storage_type=ROCArray, gamma=Float32(1.4))
     # Ensure that we do not have excessive memory allocations

From 244504d2f3b6b086cb3f4fc499a171825654c8f1 Mon Sep 17 00:00:00 2001
From: Marco Artiano <artianomarco97@gmail.com>
Date: Mon, 18 May 2026 18:54:50 +0200
Subject: [PATCH 08/29] refactoring

---
 src/solvers/dgsem_p4est/dg.jl        |  3 +++
 src/solvers/dgsem_p4est/dg_1d_gpu.jl | 36 ++++++++++++++++++++++++++++
 src/solvers/dgsem_p4est/dg_2d_gpu.jl | 36 ++++++++++++++++++++++++++++
 src/solvers/dgsem_p4est/dg_3d_gpu.jl | 36 ++++++++++++++++++++++++++++
 src/solvers/dgsem_tree/dg_1d.jl      | 28 ----------------------
 src/solvers/dgsem_tree/dg_2d.jl      | 28 ----------------------
 src/solvers/dgsem_tree/dg_3d.jl      | 29 ----------------------
 7 files changed, 111 insertions(+), 85 deletions(-)
 create mode 100644 src/solvers/dgsem_p4est/dg_1d_gpu.jl
 create mode 100644 src/solvers/dgsem_p4est/dg_2d_gpu.jl
 create mode 100644 src/solvers/dgsem_p4est/dg_3d_gpu.jl

diff --git a/src/solvers/dgsem_p4est/dg.jl b/src/solvers/dgsem_p4est/dg.jl
index 12e13f0d239..455d5d2ed32 100644
--- a/src/solvers/dgsem_p4est/dg.jl
+++ b/src/solvers/dgsem_p4est/dg.jl
@@ -85,13 +85,16 @@ end
 end
 
 include("containers.jl")
+include("dg_1d_gpu.jl")
 
 include("dg_2d.jl")
 include("dg_2d_parabolic.jl")
+include("dg_2d_gpu.jl")
 
 include("dg_3d.jl")
 include("dg_3d_parabolic.jl")
 include("dg_parallel.jl")
+include("dg_3d_gpu.jl")
 
 # Subcell limiters
 include("subcell_limiters.jl")
diff --git a/src/solvers/dgsem_p4est/dg_1d_gpu.jl b/src/solvers/dgsem_p4est/dg_1d_gpu.jl
new file mode 100644
index 00000000000..b9bdba9dfa9
--- /dev/null
+++ b/src/solvers/dgsem_p4est/dg_1d_gpu.jl
@@ -0,0 +1,36 @@
+# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
+# Since these FMAs can increase the performance of many numerical algorithms,
+# we need to opt-in explicitly.
+# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
+@muladd begin
+#! format: noindent
+
+@kernel function calc_sources_KAkernel!(du, u, t, source_terms,
+                                        node_coordinates,
+                                        equations::AbstractEquations{1}, dg, cache)
+    i, element = @index(Global, NTuple)
+    u_local = get_node_vars(u, equations, dg, i, element)
+    x_local = get_node_coords(node_coordinates, equations, dg, i, element)
+
+    du_local = source_terms(u_local, x_local, t, equations)
+
+    add_to_node_vars!(du, du_local, equations, dg, i, element)
+end
+
+function calc_sources!(backend::Backend, du, u, t, source_terms,
+                       equations::AbstractEquations{1}, dg::DG, cache)
+    nelements(dg, cache) == 0 && return nothing
+    @unpack node_coordinates = cache.elements
+    kernel_cache = kernel_filter_cache(cache)
+    kernel! = calc_sources_KAkernel!(backend)
+    kernel!(du, u, t, source_terms, node_coordinates, equations, dg, kernel_cache,
+            ndrange = (nnodes(dg), nelements(dg, cache)))
+
+    return nothing
+end
+
+function calc_sources!(backend::Backend, du, u, t, source_terms::Nothing,
+                       equations::AbstractEquations{1}, dg::DG, cache)
+    return nothing
+end
+end #muladd
diff --git a/src/solvers/dgsem_p4est/dg_2d_gpu.jl b/src/solvers/dgsem_p4est/dg_2d_gpu.jl
new file mode 100644
index 00000000000..93871032a2b
--- /dev/null
+++ b/src/solvers/dgsem_p4est/dg_2d_gpu.jl
@@ -0,0 +1,36 @@
+# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
+# Since these FMAs can increase the performance of many numerical algorithms,
+# we need to opt-in explicitly.
+# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
+@muladd begin
+#! format: noindent
+
+@kernel function calc_sources_KAkernel!(du, u, t, source_terms,
+                                        node_coordinates,
+                                        equations::AbstractEquations{2}, dg, cache)
+    i, j, element = @index(Global, NTuple)
+    u_local = get_node_vars(u, equations, dg, i, j, element)
+    x_local = get_node_coords(node_coordinates, equations, dg, i, j, element)
+
+    du_local = source_terms(u_local, x_local, t, equations)
+
+    add_to_node_vars!(du, du_local, equations, dg, i, j, element)
+end
+
+function calc_sources!(backend::Backend, du, u, t, source_terms,
+                       equations::AbstractEquations{2}, dg::DG, cache)
+    nelements(dg, cache) == 0 && return nothing
+    @unpack node_coordinates = cache.elements
+    kernel_cache = kernel_filter_cache(cache)
+    kernel! = calc_sources_KAkernel!(backend)
+    kernel!(du, u, t, source_terms, node_coordinates, equations, dg, kernel_cache,
+            ndrange = (nnodes(dg), nnodes(dg), nelements(dg, cache)))
+
+    return nothing
+end
+
+function calc_sources!(backend::Backend, du, u, t, source_terms::Nothing,
+                       equations::AbstractEquations{2}, dg::DG, cache)
+    return nothing
+end
+end #muladd
diff --git a/src/solvers/dgsem_p4est/dg_3d_gpu.jl b/src/solvers/dgsem_p4est/dg_3d_gpu.jl
new file mode 100644
index 00000000000..79e0488995e
--- /dev/null
+++ b/src/solvers/dgsem_p4est/dg_3d_gpu.jl
@@ -0,0 +1,36 @@
+# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
+# Since these FMAs can increase the performance of many numerical algorithms,
+# we need to opt-in explicitly.
+# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
+@muladd begin
+#! format: noindent
+
+@kernel function calc_sources_KAkernel!(du, u, t, source_terms,
+                                        node_coordinates,
+                                        equations::AbstractEquations{3}, dg, cache)
+    i, j, k, element = @index(Global, NTuple)
+    u_local = get_node_vars(u, equations, dg, i, j, k, element)
+    x_local = get_node_coords(node_coordinates, equations, dg, i, j, k, element)
+
+    du_local = source_terms(u_local, x_local, t, equations)
+
+    add_to_node_vars!(du, du_local, equations, dg, i, j, k, element)
+end
+
+function calc_sources!(backend::Backend, du, u, t, source_terms,
+                       equations::AbstractEquations{3}, dg::DG, cache)
+    nelements(dg, cache) == 0 && return nothing
+    @unpack node_coordinates = cache.elements
+    kernel_cache = kernel_filter_cache(cache)
+    kernel! = calc_sources_KAkernel!(backend)
+    kernel!(du, u, t, source_terms, node_coordinates, equations, dg, kernel_cache,
+            ndrange = (nnodes(dg), nnodes(dg), nnodes(dg), nelements(dg, cache)))
+
+    return nothing
+end
+
+function calc_sources!(backend::Backend, du, u, t, source_terms::Nothing,
+                       equations::AbstractEquations{3}, dg::DG, cache)
+    return nothing
+end
+end #muladd
diff --git a/src/solvers/dgsem_tree/dg_1d.jl b/src/solvers/dgsem_tree/dg_1d.jl
index f7f5f88ccd0..edd77bc7e2a 100644
--- a/src/solvers/dgsem_tree/dg_1d.jl
+++ b/src/solvers/dgsem_tree/dg_1d.jl
@@ -787,32 +787,4 @@ function calc_sources!(backend::Nothing, du, u, t, source_terms,
     return nothing
 end
 
-@kernel function calc_sources_KAkernel!(du, u, t, source_terms,
-                                        node_coordinates,
-                                        equations::AbstractEquations{1}, dg, cache)
-    i, element = @index(Global, NTuple)
-    u_local = get_node_vars(u, equations, dg, i, element)
-    x_local = get_node_coords(node_coordinates, equations, dg, i, element)
-
-    du_local = source_terms(u_local, x_local, t, equations)
-
-    add_to_node_vars!(du, du_local, equations, dg, i, element)
-end
-
-function calc_sources!(backend::Backend, du, u, t, source_terms,
-                       equations::AbstractEquations{1}, dg::DG, cache)
-    nelements(dg, cache) == 0 && return nothing
-    @unpack node_coordinates = cache.elements
-    kernel_cache = kernel_filter_cache(cache)
-    kernel! = calc_sources_KAkernel!(backend)
-    kernel!(du, u, t, source_terms, node_coordinates, equations, dg, kernel_cache,
-            ndrange = (nnodes(dg), nelements(dg, cache)))
-
-    return nothing
-end
-
-function calc_sources!(backend::Backend, du, u, t, source_terms::Nothing,
-                       equations::AbstractEquations{1}, dg::DG, cache)
-    return nothing
-end
 end # @muladd
diff --git a/src/solvers/dgsem_tree/dg_2d.jl b/src/solvers/dgsem_tree/dg_2d.jl
index 373938976bd..68d760fd502 100644
--- a/src/solvers/dgsem_tree/dg_2d.jl
+++ b/src/solvers/dgsem_tree/dg_2d.jl
@@ -1381,32 +1381,4 @@ function calc_sources!(backend::Nothing, du, u, t, source_terms,
     return nothing
 end
 
-@kernel function calc_sources_KAkernel!(du, u, t, source_terms,
-                                        node_coordinates,
-                                        equations::AbstractEquations{2}, dg, cache)
-    i, j, element = @index(Global, NTuple)
-    u_local = get_node_vars(u, equations, dg, i, j, element)
-    x_local = get_node_coords(node_coordinates, equations, dg, i, j, element)
-
-    du_local = source_terms(u_local, x_local, t, equations)
-
-    add_to_node_vars!(du, du_local, equations, dg, i, j, element)
-end
-
-function calc_sources!(backend::Backend, du, u, t, source_terms,
-                       equations::AbstractEquations{2}, dg::DG, cache)
-    nelements(dg, cache) == 0 && return nothing
-    @unpack node_coordinates = cache.elements
-    kernel_cache = kernel_filter_cache(cache)
-    kernel! = calc_sources_KAkernel!(backend)
-    kernel!(du, u, t, source_terms, node_coordinates, equations, dg, kernel_cache,
-            ndrange = (nnodes(dg), nnodes(dg), nelements(dg, cache)))
-
-    return nothing
-end
-
-function calc_sources!(backend::Backend, du, u, t, source_terms::Nothing,
-                       equations::AbstractEquations{2}, dg::DG, cache)
-    return nothing
-end
 end # @muladd
diff --git a/src/solvers/dgsem_tree/dg_3d.jl b/src/solvers/dgsem_tree/dg_3d.jl
index 06632a77edc..177dd5953a6 100644
--- a/src/solvers/dgsem_tree/dg_3d.jl
+++ b/src/solvers/dgsem_tree/dg_3d.jl
@@ -1443,33 +1443,4 @@ function calc_sources!(backend::Nothing, du, u, t, source_terms,
 
     return nothing
 end
-
-@kernel function calc_sources_KAkernel!(du, u, t, source_terms,
-                                        node_coordinates,
-                                        equations::AbstractEquations{3}, dg, cache)
-    i, j, k, element = @index(Global, NTuple)
-    u_local = get_node_vars(u, equations, dg, i, j, k, element)
-    x_local = get_node_coords(node_coordinates, equations, dg, i, j, k, element)
-
-    du_local = source_terms(u_local, x_local, t, equations)
-
-    add_to_node_vars!(du, du_local, equations, dg, i, j, k, element)
-end
-
-function calc_sources!(backend::Backend, du, u, t, source_terms,
-                       equations::AbstractEquations{3}, dg::DG, cache)
-    nelements(dg, cache) == 0 && return nothing
-    @unpack node_coordinates = cache.elements
-    kernel_cache = kernel_filter_cache(cache)
-    kernel! = calc_sources_KAkernel!(backend)
-    kernel!(du, u, t, source_terms, node_coordinates, equations, dg, kernel_cache,
-            ndrange = (nnodes(dg), nnodes(dg), nnodes(dg), nelements(dg, cache)))
-
-    return nothing
-end
-
-function calc_sources!(backend::Backend, du, u, t, source_terms::Nothing,
-                       equations::AbstractEquations{3}, dg::DG, cache)
-    return nothing
-end
 end # @muladd

From d897a73a17dd4a617eca436c5ac96bb4b4f0f470 Mon Sep 17 00:00:00 2001
From: Marco Artiano <artianomarco97@gmail.com>
Date: Mon, 18 May 2026 18:55:14 +0200
Subject: [PATCH 09/29] format

---
 src/solvers/dgsem_tree/dg_1d.jl | 1 -
 src/solvers/dgsem_tree/dg_2d.jl | 1 -
 2 files changed, 2 deletions(-)

diff --git a/src/solvers/dgsem_tree/dg_1d.jl b/src/solvers/dgsem_tree/dg_1d.jl
index edd77bc7e2a..8f195490767 100644
--- a/src/solvers/dgsem_tree/dg_1d.jl
+++ b/src/solvers/dgsem_tree/dg_1d.jl
@@ -786,5 +786,4 @@ function calc_sources!(backend::Nothing, du, u, t, source_terms,
 
     return nothing
 end
-
 end # @muladd
diff --git a/src/solvers/dgsem_tree/dg_2d.jl b/src/solvers/dgsem_tree/dg_2d.jl
index 68d760fd502..29c881943ae 100644
--- a/src/solvers/dgsem_tree/dg_2d.jl
+++ b/src/solvers/dgsem_tree/dg_2d.jl
@@ -1380,5 +1380,4 @@ function calc_sources!(backend::Nothing, du, u, t, source_terms,
 
     return nothing
 end
-
 end # @muladd

From 74e140da6fcbcf15d7ba46b898ac4c239ac8d37d Mon Sep 17 00:00:00 2001
From: Marco Artiano <artianomarco97@gmail.com>
Date: Mon, 18 May 2026 19:28:59 +0200
Subject: [PATCH 10/29] refactor also 3D kernels

---
 src/solvers/dgsem_p4est/dg_3d.jl     | 76 ----------------------------
 src/solvers/dgsem_p4est/dg_3d_gpu.jl | 76 ++++++++++++++++++++++++++++
 2 files changed, 76 insertions(+), 76 deletions(-)

diff --git a/src/solvers/dgsem_p4est/dg_3d.jl b/src/solvers/dgsem_p4est/dg_3d.jl
index faac5fa46f1..ab3b910925c 100644
--- a/src/solvers/dgsem_p4est/dg_3d.jl
+++ b/src/solvers/dgsem_p4est/dg_3d.jl
@@ -106,27 +106,6 @@ function prolong2interfaces!(backend::Nothing, cache, u,
     return nothing
 end
 
-function prolong2interfaces!(backend::Backend, cache, u,
-                             mesh::Union{P4estMesh{3}, T8codeMesh{3}},
-                             equations, dg::DG)
-    @unpack interfaces = cache
-    @unpack neighbor_ids, node_indices = cache.interfaces
-    index_range = eachnode(dg)
-
-    kernel! = prolong2interfaces_KAkernel!(backend)
-    kernel!(interfaces.u, u, typeof(mesh), equations, neighbor_ids, node_indices,
-            index_range,
-            ndrange = ninterfaces(interfaces))
-    return nothing
-end
-
-@kernel function prolong2interfaces_KAkernel!(interface_u, u, MeshT, equations,
-                                              neighbor_ids, node_indices, index_range)
-    interface = @index(Global)
-    prolong2interfaces_per_interface!(interface_u, u, MeshT, equations, neighbor_ids,
-                                      node_indices, index_range, interface)
-end
-
 @inline function prolong2interfaces_per_interface!(u_interface, u,
                                                    ::Type{<:Union{P4estMesh{3},
                                                                   T8codeMesh{3}}},
@@ -223,38 +202,6 @@ function calc_interface_flux!(backend::Nothing, surface_flux_values,
     return nothing
 end
 
-function calc_interface_flux!(backend::Backend, surface_flux_values,
-                              mesh::Union{P4estMesh{3}, T8codeMesh{3}},
-                              have_nonconservative_terms,
-                              equations, surface_integral, dg::DG, cache)
-    @unpack neighbor_ids, node_indices = cache.interfaces
-    @unpack contravariant_vectors = cache.elements
-    index_range = eachnode(dg)
-
-    kernel! = calc_interface_flux_KAkernel!(backend)
-    kernel!(surface_flux_values, typeof(mesh), have_nonconservative_terms, equations,
-            surface_integral, typeof(dg), cache.interfaces.u,
-            neighbor_ids, node_indices, contravariant_vectors, index_range,
-            ndrange = ninterfaces(cache.interfaces))
-    return nothing
-end
-
-@kernel function calc_interface_flux_KAkernel!(surface_flux_values, MeshT,
-                                               have_nonconservative_terms, equations,
-                                               surface_integral, SolverT, u_interface,
-                                               neighbor_ids, node_indices,
-                                               contravariant_vectors, index_range)
-    interface = @index(Global)
-    calc_interface_flux_per_interface!(surface_flux_values,
-                                       MeshT,
-                                       have_nonconservative_terms,
-                                       equations, surface_integral, SolverT,
-                                       u_interface,
-                                       neighbor_ids, node_indices,
-                                       contravariant_vectors,
-                                       index_range, interface)
-end
-
 @inline function calc_interface_flux_per_interface!(surface_flux_values,
                                                     MeshT::Type{<:Union{P4estMesh{3},
                                                                         T8codeMesh{3}}},
@@ -1026,29 +973,6 @@ function calc_surface_integral!(backend::Nothing, du, u,
     return nothing
 end
 
-function calc_surface_integral!(backend::Backend, du, u,
-                                mesh::Union{P4estMesh{3}, T8codeMesh{3}},
-                                equations,
-                                surface_integral::SurfaceIntegralWeakForm,
-                                dg::DGSEM, cache)
-    @unpack inverse_weights = dg.basis
-    @unpack surface_flux_values = cache.elements
-
-    kernel! = calc_surface_integral_KAkernel!(backend)
-    kernel!(du, typeof(mesh), equations, surface_integral, dg, inverse_weights[1],
-            surface_flux_values, ndrange = nelements(cache.elements))
-    return nothing
-end
-
-@kernel function calc_surface_integral_KAkernel!(du, MeshT, equations,
-                                                 surface_integral, dg, factor,
-                                                 surface_flux_values)
-    element = @index(Global)
-    calc_surface_integral_per_element!(du, MeshT,
-                                       equations, surface_integral, dg, factor,
-                                       surface_flux_values, element)
-end
-
 @inline function calc_surface_integral_per_element!(du,
                                                     ::Type{<:Union{P4estMesh{3},
                                                                    T8codeMesh{3}}},
diff --git a/src/solvers/dgsem_p4est/dg_3d_gpu.jl b/src/solvers/dgsem_p4est/dg_3d_gpu.jl
index 79e0488995e..046cf498dac 100644
--- a/src/solvers/dgsem_p4est/dg_3d_gpu.jl
+++ b/src/solvers/dgsem_p4est/dg_3d_gpu.jl
@@ -5,6 +5,82 @@
 @muladd begin
 #! format: noindent
 
+function prolong2interfaces!(backend::Backend, cache, u,
+                             mesh::Union{P4estMesh{3}, T8codeMesh{3}},
+                             equations, dg::DG)
+    @unpack interfaces = cache
+    @unpack neighbor_ids, node_indices = cache.interfaces
+    index_range = eachnode(dg)
+
+    kernel! = prolong2interfaces_KAkernel!(backend)
+    kernel!(interfaces.u, u, typeof(mesh), equations, neighbor_ids, node_indices,
+            index_range,
+            ndrange = ninterfaces(interfaces))
+    return nothing
+end
+
+@kernel function prolong2interfaces_KAkernel!(interface_u, u, MeshT, equations,
+                                              neighbor_ids, node_indices, index_range)
+    interface = @index(Global)
+    prolong2interfaces_per_interface!(interface_u, u, MeshT, equations, neighbor_ids,
+                                      node_indices, index_range, interface)
+end
+
+function calc_surface_integral!(backend::Backend, du, u,
+                                mesh::Union{P4estMesh{3}, T8codeMesh{3}},
+                                equations,
+                                surface_integral::SurfaceIntegralWeakForm,
+                                dg::DGSEM, cache)
+    @unpack inverse_weights = dg.basis
+    @unpack surface_flux_values = cache.elements
+
+    kernel! = calc_surface_integral_KAkernel!(backend)
+    kernel!(du, typeof(mesh), equations, surface_integral, dg, inverse_weights[1],
+            surface_flux_values, ndrange = nelements(cache.elements))
+    return nothing
+end
+
+@kernel function calc_surface_integral_KAkernel!(du, MeshT, equations,
+                                                 surface_integral, dg, factor,
+                                                 surface_flux_values)
+    element = @index(Global)
+    calc_surface_integral_per_element!(du, MeshT,
+                                       equations, surface_integral, dg, factor,
+                                       surface_flux_values, element)
+end
+
+function calc_interface_flux!(backend::Backend, surface_flux_values,
+                              mesh::Union{P4estMesh{3}, T8codeMesh{3}},
+                              have_nonconservative_terms,
+                              equations, surface_integral, dg::DG, cache)
+    @unpack neighbor_ids, node_indices = cache.interfaces
+    @unpack contravariant_vectors = cache.elements
+    index_range = eachnode(dg)
+
+    kernel! = calc_interface_flux_KAkernel!(backend)
+    kernel!(surface_flux_values, typeof(mesh), have_nonconservative_terms, equations,
+            surface_integral, typeof(dg), cache.interfaces.u,
+            neighbor_ids, node_indices, contravariant_vectors, index_range,
+            ndrange = ninterfaces(cache.interfaces))
+    return nothing
+end
+
+@kernel function calc_interface_flux_KAkernel!(surface_flux_values, MeshT,
+                                               have_nonconservative_terms, equations,
+                                               surface_integral, SolverT, u_interface,
+                                               neighbor_ids, node_indices,
+                                               contravariant_vectors, index_range)
+    interface = @index(Global)
+    calc_interface_flux_per_interface!(surface_flux_values,
+                                       MeshT,
+                                       have_nonconservative_terms,
+                                       equations, surface_integral, SolverT,
+                                       u_interface,
+                                       neighbor_ids, node_indices,
+                                       contravariant_vectors,
+                                       index_range, interface)
+end
+
 @kernel function calc_sources_KAkernel!(du, u, t, source_terms,
                                         node_coordinates,
                                         equations::AbstractEquations{3}, dg, cache)

From f1209fe297696b836baeedf36867b5c0a50bfdd6 Mon Sep 17 00:00:00 2001
From: Marco Artiano <artianomarco97@gmail.com>
Date: Mon, 18 May 2026 22:18:05 +0200
Subject: [PATCH 11/29] format

---
 src/solvers/dgsem_p4est/dg_2d_gpu.jl | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/solvers/dgsem_p4est/dg_2d_gpu.jl b/src/solvers/dgsem_p4est/dg_2d_gpu.jl
index bb64111589e..17c5d3c2b9d 100644
--- a/src/solvers/dgsem_p4est/dg_2d_gpu.jl
+++ b/src/solvers/dgsem_p4est/dg_2d_gpu.jl
@@ -118,8 +118,8 @@ end
     apply_jacobian_per_quadrature_node!(du, MeshT, equations, dg, inverse_jacobian,
                                         i, j, element)
 end
-  
- @kernel function calc_sources_KAkernel!(du, u, t, source_terms,
+
+@kernel function calc_sources_KAkernel!(du, u, t, source_terms,
                                         node_coordinates,
                                         equations::AbstractEquations{2}, dg, cache)
     i, j, element = @index(Global, NTuple)
@@ -147,4 +147,4 @@ function calc_sources!(backend::Backend, du, u, t, source_terms::Nothing,
                        equations::AbstractEquations{2}, dg::DG, cache)
     return nothing
 end
-end #muladd
\ No newline at end of file
+end #muladd

From c9780365e6d4d7ab21af3547bd5e7b1d30c594f6 Mon Sep 17 00:00:00 2001
From: Marco Artiano <artianomarco97@gmail.com>
Date: Mon, 18 May 2026 22:58:18 +0200
Subject: [PATCH 12/29] fix tests

---
 test/test_amdgpu_2d.jl | 12 ++++++++----
 test/test_amdgpu_3d.jl | 14 ++++++++++----
 2 files changed, 18 insertions(+), 8 deletions(-)

diff --git a/test/test_amdgpu_2d.jl b/test/test_amdgpu_2d.jl
index 2d641906cd9..fba3828380a 100644
--- a/test/test_amdgpu_2d.jl
+++ b/test/test_amdgpu_2d.jl
@@ -108,10 +108,14 @@ end
     using AMDGPU
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_source_terms_gpu.jl"),
                         # Expected errors are exactly the same as with TreeMesh!
-                        l2=Float32[2.4917018095933837e-6, 2.7148269885239423e-6,
-                                   2.695290306860358e-6, 6.243861976167833e-6],
-                        linf=Float32[1.6489475493930428e-5, 1.7499923706143505e-5,
-                                     1.893043518075288e-5, 6.214141845717336e-5],
+                        l2=Float32[2.4917018095933837e-6,
+                                   2.7148269885239423e-6,
+                                   2.695290306860358e-6,
+                                   6.243861976167833e-6],
+                        linf=Float32[1.6489475493930428e-5,
+                                     1.7499923706143505e-5,
+                                     1.893043518075288e-5,
+                                     6.214141845717336e-5],
                         RealT_for_test_tolerances=Float32,
                         real_type=Float32,
                         storage_type=ROCArray, gamma=Float32(1.4))
diff --git a/test/test_amdgpu_3d.jl b/test/test_amdgpu_3d.jl
index f0548b056fc..f6d734e1f7f 100644
--- a/test/test_amdgpu_3d.jl
+++ b/test/test_amdgpu_3d.jl
@@ -108,11 +108,17 @@ end
     using AMDGPU
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_source_terms_gpu.jl"),
                         # Expected errors are exactly the same as with TreeMesh!
-                        l2=Float32[4.912578089985958e-5, 5.3683407014580115e-5,
-                                   5.368099834769191e-5, 5.371664525206341e-5,
+                        l2=Float32[4.912578089985958e-5,
+                                   5.3683407014580115e-5,
+                                   5.368099834769191e-5,
+                                   5.371664525206341e-5,
                                    0.00015186256300882088],
-                        linf=Float32[0.00032772542853032327, 0.00035144807715092874, 0.0003549051465479014, 0.00035573961157475686, 0.0013591384887696734]
-                        RealT_for_test_tolerances = Float32,
+                        linf=Float32[0.00032772542853032327,
+                                     0.00035144807715092874,
+                                     0.0003549051465479014,
+                                     0.00035573961157475686,
+                                     0.0013591384887696734],
+                        RealT_for_test_tolerances=Float32,
                         real_type=Float32,
                         storage_type=ROCArray, gamma=Float32(1.4))
     # Ensure that we do not have excessive memory allocations

From 5a5948f568423e356d1ee17923ae10f247fa1d2a Mon Sep 17 00:00:00 2001
From: Marco Artiano <artianomarco97@gmail.com>
Date: Tue, 19 May 2026 07:29:52 +0200
Subject: [PATCH 13/29] fix tests

---
 test/test_amdgpu_2d.jl | 4 ++--
 test/test_amdgpu_3d.jl | 4 ++--
 2 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/test/test_amdgpu_2d.jl b/test/test_amdgpu_2d.jl
index fba3828380a..f2cf860de9c 100644
--- a/test/test_amdgpu_2d.jl
+++ b/test/test_amdgpu_2d.jl
@@ -68,7 +68,7 @@ end
 end
 
 @trixi_testset "elixir_euler_source_terms_gpu.jl native" begin
-    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_source_terms_gpu.jl"),
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms_gpu.jl"),
                         # Expected errors are exactly the same as with TreeMesh!
                         l2=[
                             4.893619139889976e-5,
@@ -106,7 +106,7 @@ end
 @trixi_testset "elixir_source_terms_gpu.jl Float32 / AMDGPU" begin
     # Using AMDGPU inside the testset since otherwise the bindings are hiddend by the anonymous modules
     using AMDGPU
-    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_source_terms_gpu.jl"),
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms_gpu.jl"),
                         # Expected errors are exactly the same as with TreeMesh!
                         l2=Float32[2.4917018095933837e-6,
                                    2.7148269885239423e-6,
diff --git a/test/test_amdgpu_3d.jl b/test/test_amdgpu_3d.jl
index f6d734e1f7f..91d39ae253d 100644
--- a/test/test_amdgpu_3d.jl
+++ b/test/test_amdgpu_3d.jl
@@ -68,7 +68,7 @@ end
 end
 
 @trixi_testset "elixir_euler_source_terms_gpu.jl native" begin
-    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_source_terms_gpu.jl"),
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms_gpu.jl"),
                         # Expected errors are exactly the same as with TreeMesh!
                         l2=[
                             4.893619139889976e-5,
@@ -106,7 +106,7 @@ end
 @trixi_testset "elixir_source_terms_gpu.jl Float32 / AMDGPU" begin
     # Using AMDGPU inside the testset since otherwise the bindings are hiddend by the anonymous modules
     using AMDGPU
-    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_source_terms_gpu.jl"),
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms_gpu.jl"),
                         # Expected errors are exactly the same as with TreeMesh!
                         l2=Float32[4.912578089985958e-5,
                                    5.3683407014580115e-5,

From d66c664f762188adf3b8a6f5f2483e300b6db03c Mon Sep 17 00:00:00 2001
From: Marco Artiano <artianomarco97@gmail.com>
Date: Tue, 19 May 2026 08:16:26 +0200
Subject: [PATCH 14/29] fix elixirs

---
 examples/p4est_2d_dgsem/elixir_euler_source_terms_gpu.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/examples/p4est_2d_dgsem/elixir_euler_source_terms_gpu.jl b/examples/p4est_2d_dgsem/elixir_euler_source_terms_gpu.jl
index 755227e5985..a4de88b34ed 100644
--- a/examples/p4est_2d_dgsem/elixir_euler_source_terms_gpu.jl
+++ b/examples/p4est_2d_dgsem/elixir_euler_source_terms_gpu.jl
@@ -37,7 +37,7 @@ semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver;
 # ODE solvers, callbacks etc.
 
 tspan = (0.0, 2.0)
-ode = semidiscretize(semi, tspan; real_type = real_type, storage_type = storage_type)
+ode = semidiscretize(semi, tspan; real_type = nothing, storage_type = nothing)
 
 summary_callback = SummaryCallback()
 

From b9bbea5dbc4764a6c38c21cbd50be83492230a2b Mon Sep 17 00:00:00 2001
From: Marco Artiano <artianomarco97@gmail.com>
Date: Tue, 19 May 2026 08:54:06 +0200
Subject: [PATCH 15/29] fix tests

---
 test/test_amdgpu_2d.jl | 24 +++++++++---------------
 test/test_amdgpu_3d.jl |  2 +-
 2 files changed, 10 insertions(+), 16 deletions(-)

diff --git a/test/test_amdgpu_2d.jl b/test/test_amdgpu_2d.jl
index f2cf860de9c..17e73d2e052 100644
--- a/test/test_amdgpu_2d.jl
+++ b/test/test_amdgpu_2d.jl
@@ -70,20 +70,14 @@ end
 @trixi_testset "elixir_euler_source_terms_gpu.jl native" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms_gpu.jl"),
                         # Expected errors are exactly the same as with TreeMesh!
-                        l2=[
-                            4.893619139889976e-5,
-                            5.3526950567182756e-5,
-                            5.35269505672133e-5,
-                            5.352695056735998e-5,
-                            0.00015172095200428318
-                        ],
-                        linf=[
-                            0.00031179856625374036,
-                            0.0003368725355339386,
-                            0.0003368725355383795,
-                            0.00033687253560787944,
-                            0.0013193387520935573
-                        ])
+                        l2=[9.321181254378498e-7,
+                            1.418121074369651e-6,
+                            1.4181210743821669e-6,
+                            4.824553091168877e-6],
+                        linf=[9.577246532499473e-6,
+                            1.1707525985116263e-5,
+                            1.1707525982673772e-5,
+                            4.886961559069647e-5])
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
     semi = ode.p # `semidiscretize` adapts the semi, so we need to obtain it from the ODE problem.
@@ -103,7 +97,7 @@ end
     @test Trixi.storage_type(ode.p.cache.mortars) === Array
 end
 
-@trixi_testset "elixir_source_terms_gpu.jl Float32 / AMDGPU" begin
+@trixi_testset "elixir_euler_source_terms_gpu.jl Float32 / AMDGPU" begin
     # Using AMDGPU inside the testset since otherwise the bindings are hiddend by the anonymous modules
     using AMDGPU
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms_gpu.jl"),
diff --git a/test/test_amdgpu_3d.jl b/test/test_amdgpu_3d.jl
index 91d39ae253d..e8d9b7f9e61 100644
--- a/test/test_amdgpu_3d.jl
+++ b/test/test_amdgpu_3d.jl
@@ -103,7 +103,7 @@ end
     @test Trixi.storage_type(ode.p.cache.mortars) === Array
 end
 
-@trixi_testset "elixir_source_terms_gpu.jl Float32 / AMDGPU" begin
+@trixi_testset "elixir_euler_source_terms_gpu.jl Float32 / AMDGPU" begin
     # Using AMDGPU inside the testset since otherwise the bindings are hiddend by the anonymous modules
     using AMDGPU
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms_gpu.jl"),

From 55ded1f0375813de9e0967b33c3d664bb8303baf Mon Sep 17 00:00:00 2001
From: Marco Artiano <artianomarco97@gmail.com>
Date: Tue, 19 May 2026 09:08:30 +0200
Subject: [PATCH 16/29] fix tests

---
 examples/p4est_3d_dgsem/elixir_euler_source_terms_gpu.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/examples/p4est_3d_dgsem/elixir_euler_source_terms_gpu.jl b/examples/p4est_3d_dgsem/elixir_euler_source_terms_gpu.jl
index 6dda9dcc8fb..3e1aee8be47 100644
--- a/examples/p4est_3d_dgsem/elixir_euler_source_terms_gpu.jl
+++ b/examples/p4est_3d_dgsem/elixir_euler_source_terms_gpu.jl
@@ -39,7 +39,7 @@ semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver;
 # ODE solvers, callbacks etc.
 
 tspan = (0.0, 5.0)
-ode = semidiscretize(semi, tspan; real_type = real_type, storage_type = storage_type)
+ode = semidiscretize(semi, tspan; real_type = nothing, storage_type = nothing)
 
 summary_callback = SummaryCallback()
 

From 303cd35097fdf165949f36960e55118321c398fc Mon Sep 17 00:00:00 2001
From: Marco Artiano <57838732+MarcoArtiano@users.noreply.github.com>
Date: Tue, 19 May 2026 15:38:22 +0200
Subject: [PATCH 17/29] Apply suggestions from code review

Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com>
---
 examples/p4est_2d_dgsem/elixir_euler_source_terms_gpu.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/examples/p4est_2d_dgsem/elixir_euler_source_terms_gpu.jl b/examples/p4est_2d_dgsem/elixir_euler_source_terms_gpu.jl
index a4de88b34ed..cdc8d2bcc86 100644
--- a/examples/p4est_2d_dgsem/elixir_euler_source_terms_gpu.jl
+++ b/examples/p4est_2d_dgsem/elixir_euler_source_terms_gpu.jl
@@ -1,5 +1,5 @@
 # The same setup as tree_2d_dgsem/elixir_euler_source_terms.jl
-# to verify the StructuredMesh implementation against TreeMesh
+# to verify the P4estMesh implementation against TreeMesh
 
 using OrdinaryDiffEqLowStorageRK
 using Trixi

From 8733fddcb13a3b9509b7902a9211c322e1b50756 Mon Sep 17 00:00:00 2001
From: MarcoArtiano <artianomarco97@gmail.com>
Date: Tue, 19 May 2026 15:45:11 +0200
Subject: [PATCH 18/29] rename tests, delete 1D P4est

---
 ...ms_gpu.jl => elixir_euler_source_terms.jl} |  0
 ...ms_gpu.jl => elixir_euler_source_terms.jl} |  0
 src/solvers/dgsem_p4est/dg.jl                 |  1 -
 src/solvers/dgsem_p4est/dg_1d_gpu.jl          | 36 -------------------
 test/test_amdgpu_2d.jl                        |  8 ++---
 test/test_amdgpu_3d.jl                        |  8 ++---
 6 files changed, 8 insertions(+), 45 deletions(-)
 rename examples/p4est_2d_dgsem/{elixir_euler_source_terms_gpu.jl => elixir_euler_source_terms.jl} (100%)
 rename examples/p4est_3d_dgsem/{elixir_euler_source_terms_gpu.jl => elixir_euler_source_terms.jl} (100%)
 delete mode 100644 src/solvers/dgsem_p4est/dg_1d_gpu.jl

diff --git a/examples/p4est_2d_dgsem/elixir_euler_source_terms_gpu.jl b/examples/p4est_2d_dgsem/elixir_euler_source_terms.jl
similarity index 100%
rename from examples/p4est_2d_dgsem/elixir_euler_source_terms_gpu.jl
rename to examples/p4est_2d_dgsem/elixir_euler_source_terms.jl
diff --git a/examples/p4est_3d_dgsem/elixir_euler_source_terms_gpu.jl b/examples/p4est_3d_dgsem/elixir_euler_source_terms.jl
similarity index 100%
rename from examples/p4est_3d_dgsem/elixir_euler_source_terms_gpu.jl
rename to examples/p4est_3d_dgsem/elixir_euler_source_terms.jl
diff --git a/src/solvers/dgsem_p4est/dg.jl b/src/solvers/dgsem_p4est/dg.jl
index 455d5d2ed32..713cb491f03 100644
--- a/src/solvers/dgsem_p4est/dg.jl
+++ b/src/solvers/dgsem_p4est/dg.jl
@@ -85,7 +85,6 @@ end
 end
 
 include("containers.jl")
-include("dg_1d_gpu.jl")
 
 include("dg_2d.jl")
 include("dg_2d_parabolic.jl")
diff --git a/src/solvers/dgsem_p4est/dg_1d_gpu.jl b/src/solvers/dgsem_p4est/dg_1d_gpu.jl
deleted file mode 100644
index b9bdba9dfa9..00000000000
--- a/src/solvers/dgsem_p4est/dg_1d_gpu.jl
+++ /dev/null
@@ -1,36 +0,0 @@
-# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
-# Since these FMAs can increase the performance of many numerical algorithms,
-# we need to opt-in explicitly.
-# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
-@muladd begin
-#! format: noindent
-
-@kernel function calc_sources_KAkernel!(du, u, t, source_terms,
-                                        node_coordinates,
-                                        equations::AbstractEquations{1}, dg, cache)
-    i, element = @index(Global, NTuple)
-    u_local = get_node_vars(u, equations, dg, i, element)
-    x_local = get_node_coords(node_coordinates, equations, dg, i, element)
-
-    du_local = source_terms(u_local, x_local, t, equations)
-
-    add_to_node_vars!(du, du_local, equations, dg, i, element)
-end
-
-function calc_sources!(backend::Backend, du, u, t, source_terms,
-                       equations::AbstractEquations{1}, dg::DG, cache)
-    nelements(dg, cache) == 0 && return nothing
-    @unpack node_coordinates = cache.elements
-    kernel_cache = kernel_filter_cache(cache)
-    kernel! = calc_sources_KAkernel!(backend)
-    kernel!(du, u, t, source_terms, node_coordinates, equations, dg, kernel_cache,
-            ndrange = (nnodes(dg), nelements(dg, cache)))
-
-    return nothing
-end
-
-function calc_sources!(backend::Backend, du, u, t, source_terms::Nothing,
-                       equations::AbstractEquations{1}, dg::DG, cache)
-    return nothing
-end
-end #muladd
diff --git a/test/test_amdgpu_2d.jl b/test/test_amdgpu_2d.jl
index 17e73d2e052..4738a208a58 100644
--- a/test/test_amdgpu_2d.jl
+++ b/test/test_amdgpu_2d.jl
@@ -67,8 +67,8 @@ end
     @test Trixi.storage_type(ode.p.cache.mortars) === ROCArray
 end
 
-@trixi_testset "elixir_euler_source_terms_gpu.jl native" begin
-    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms_gpu.jl"),
+@trixi_testset "elixir_euler_source_terms.jl native" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms.jl"),
                         # Expected errors are exactly the same as with TreeMesh!
                         l2=[9.321181254378498e-7,
                             1.418121074369651e-6,
@@ -97,10 +97,10 @@ end
     @test Trixi.storage_type(ode.p.cache.mortars) === Array
 end
 
-@trixi_testset "elixir_euler_source_terms_gpu.jl Float32 / AMDGPU" begin
+@trixi_testset "elixir_euler_source_terms.jl Float32 / AMDGPU" begin
     # Using AMDGPU inside the testset since otherwise the bindings are hiddend by the anonymous modules
     using AMDGPU
-    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms_gpu.jl"),
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms.jl"),
                         # Expected errors are exactly the same as with TreeMesh!
                         l2=Float32[2.4917018095933837e-6,
                                    2.7148269885239423e-6,
diff --git a/test/test_amdgpu_3d.jl b/test/test_amdgpu_3d.jl
index e8d9b7f9e61..828c2f5efbe 100644
--- a/test/test_amdgpu_3d.jl
+++ b/test/test_amdgpu_3d.jl
@@ -67,8 +67,8 @@ end
     @test Trixi.storage_type(ode.p.cache.mortars) === ROCArray
 end
 
-@trixi_testset "elixir_euler_source_terms_gpu.jl native" begin
-    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms_gpu.jl"),
+@trixi_testset "elixir_euler_source_terms.jl native" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms.jl"),
                         # Expected errors are exactly the same as with TreeMesh!
                         l2=[
                             4.893619139889976e-5,
@@ -103,10 +103,10 @@ end
     @test Trixi.storage_type(ode.p.cache.mortars) === Array
 end
 
-@trixi_testset "elixir_euler_source_terms_gpu.jl Float32 / AMDGPU" begin
+@trixi_testset "elixir_euler_source_terms.jl Float32 / AMDGPU" begin
     # Using AMDGPU inside the testset since otherwise the bindings are hiddend by the anonymous modules
     using AMDGPU
-    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms_gpu.jl"),
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms.jl"),
                         # Expected errors are exactly the same as with TreeMesh!
                         l2=Float32[4.912578089985958e-5,
                                    5.3683407014580115e-5,

From 01c2764ce15929bfa30f963a4d214412213af37f Mon Sep 17 00:00:00 2001
From: MarcoArtiano <artianomarco97@gmail.com>
Date: Tue, 19 May 2026 15:49:17 +0200
Subject: [PATCH 19/29] add kernel abstraction tests

---
 test/test_kernelabstractions.jl | 86 +++++++++++++++++++++++++++++++++
 1 file changed, 86 insertions(+)

diff --git a/test/test_kernelabstractions.jl b/test/test_kernelabstractions.jl
index a1a771ee402..7a714b7a96b 100644
--- a/test/test_kernelabstractions.jl
+++ b/test/test_kernelabstractions.jl
@@ -44,6 +44,45 @@ end
     semi = ode.p # `semidiscretize` adapts the semi, so we need to obtain it from the ODE problem.
     @test_allocations(Trixi.rhs!, ode.p, sol, 60_000)
 end
+
+@trixi_testset "elixir_euler_source_terms.jl native" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms.jl"),
+                        # Expected errors are exactly the same as with TreeMesh!
+                        l2=[9.321181254378498e-7,
+                            1.418121074369651e-6,
+                            1.4181210743821669e-6,
+                            4.824553091168877e-6],
+                        linf=[9.577246532499473e-6,
+                            1.1707525985116263e-5,
+                            1.1707525982673772e-5,
+                            4.886961559069647e-5])
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
+    semi = ode.p # `semidiscretize` adapts the semi, so we need to obtain it from the ODE problem.
+    @test_allocations(Trixi.rhs!, semi, sol, 1000)
+end
+
+@trixi_testset "elixir_euler_source_terms.jl Float32 / AMDGPU" begin
+    # Using AMDGPU inside the testset since otherwise the bindings are hiddend by the anonymous modules
+    using AMDGPU
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms.jl"),
+                        # Expected errors are exactly the same as with TreeMesh!
+                        l2=Float32[2.4917018095933837e-6,
+                                   2.7148269885239423e-6,
+                                   2.695290306860358e-6,
+                                   6.243861976167833e-6],
+                        linf=Float32[1.6489475493930428e-5,
+                                     1.7499923706143505e-5,
+                                     1.893043518075288e-5,
+                                     6.214141845717336e-5],
+                        RealT_for_test_tolerances=Float32,
+                        real_type=Float32,
+                        gamma=Float32(1.4))
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
+    semi = ode.p # `semidiscretize` adapts the semi, so we need to obtain it from the ODE problem.
+    @test_allocations(Trixi.rhs!, semi, sol, 100_000)
+end
 end
 
 @testset "KernelAbstractions CPU 3D" begin
@@ -74,6 +113,53 @@ end
     semi = ode.p # `semidiscretize` adapts the semi, so we need to obtain it from the ODE problem.
     @test_allocations(Trixi.rhs!, semi, sol, 370_000)
 end
+
+@trixi_testset "elixir_euler_source_terms.jl native" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms.jl"),
+                        # Expected errors are exactly the same as with TreeMesh!
+                        l2=[
+                            4.893619139889976e-5,
+                            5.3526950567182756e-5,
+                            5.35269505672133e-5,
+                            5.352695056735998e-5,
+                            0.00015172095200428318
+                        ],
+                        linf=[
+                            0.00031179856625374036,
+                            0.0003368725355339386,
+                            0.0003368725355383795,
+                            0.00033687253560787944,
+                            0.0013193387520935573
+                        ])
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
+    semi = ode.p # `semidiscretize` adapts the semi, so we need to obtain it from the ODE problem.
+    @test_allocations(Trixi.rhs!, semi, sol, 1000)
+end
+
+@trixi_testset "elixir_euler_source_terms.jl Float32 / AMDGPU" begin
+    # Using AMDGPU inside the testset since otherwise the bindings are hiddend by the anonymous modules
+    using AMDGPU
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms.jl"),
+                        # Expected errors are exactly the same as with TreeMesh!
+                        l2=Float32[4.912578089985958e-5,
+                                   5.3683407014580115e-5,
+                                   5.368099834769191e-5,
+                                   5.371664525206341e-5,
+                                   0.00015186256300882088],
+                        linf=Float32[0.00032772542853032327,
+                                     0.00035144807715092874,
+                                     0.0003549051465479014,
+                                     0.00035573961157475686,
+                                     0.0013591384887696734],
+                        RealT_for_test_tolerances=Float32,
+                        real_type=Float32,
+                        gamma=Float32(1.4))
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
+    semi = ode.p # `semidiscretize` adapts the semi, so we need to obtain it from the ODE problem.
+    @test_allocations(Trixi.rhs!, semi, sol, 100_000)
+end
 end
 
 # Clean up afterwards: delete Trixi.jl output directory

From cf75145c6f6a65aa6b87866905544b43b5f48acb Mon Sep 17 00:00:00 2001
From: MarcoArtiano <artianomarco97@gmail.com>
Date: Tue, 19 May 2026 15:53:27 +0200
Subject: [PATCH 20/29] add cuda tests

---
 test/test_cuda_2d.jl | 65 +++++++++++++++++++++++++++++++++++++++
 test/test_cuda_3d.jl | 73 ++++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 138 insertions(+)

diff --git a/test/test_cuda_2d.jl b/test/test_cuda_2d.jl
index 85fcb75139f..02fb12a0ba0 100644
--- a/test/test_cuda_2d.jl
+++ b/test/test_cuda_2d.jl
@@ -67,6 +67,71 @@ end
     @test Trixi.storage_type(ode.p.cache.mortars) === CuArray
 end
 
+@trixi_testset "elixir_euler_source_terms.jl native" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms.jl"),
+                        # Expected errors are exactly the same as with TreeMesh!
+                        l2=[9.321181254378498e-7,
+                            1.418121074369651e-6,
+                            1.4181210743821669e-6,
+                            4.824553091168877e-6],
+                        linf=[9.577246532499473e-6,
+                            1.1707525985116263e-5,
+                            1.1707525982673772e-5,
+                            4.886961559069647e-5])
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
+    semi = ode.p # `semidiscretize` adapts the semi, so we need to obtain it from the ODE problem.
+    @test_allocations(Trixi.rhs!, semi, sol, 1000)
+    @test real(ode.p.solver) == Float64
+    @test real(ode.p.solver.basis) == Float64
+    @test real(ode.p.solver.mortar) == Float64
+    # TODO: `mesh` is currently not `adapt`ed correctly
+    @test real(ode.p.mesh) == Float64
+
+    @test ode.u0 isa Array
+    @test ode.p.solver.basis.derivative_matrix isa Array
+
+    @test Trixi.storage_type(ode.p.cache.elements) === Array
+    @test Trixi.storage_type(ode.p.cache.interfaces) === Array
+    @test Trixi.storage_type(ode.p.cache.boundaries) === Array
+    @test Trixi.storage_type(ode.p.cache.mortars) === Array
+end
+
+@trixi_testset "elixir_euler_source_terms.jl Float32 / CUDA" begin
+    # Using CUDA inside the testset since otherwise the bindings are hiddend by the anonymous modules
+    using CUDA
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms.jl"),
+                        # Expected errors are exactly the same as with TreeMesh!
+                        l2=Float32[2.4917018095933837e-6,
+                                   2.7148269885239423e-6,
+                                   2.695290306860358e-6,
+                                   6.243861976167833e-6],
+                        linf=Float32[1.6489475493930428e-5,
+                                     1.7499923706143505e-5,
+                                     1.893043518075288e-5,
+                                     6.214141845717336e-5],
+                        RealT_for_test_tolerances=Float32,
+                        real_type=Float32,
+                        storage_type=CuArray, gamma=Float32(1.4))
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
+    semi = ode.p # `semidiscretize` adapts the semi, so we need to obtain it from the ODE problem.
+    @test_allocations(Trixi.rhs!, semi, sol, 100_000)
+    @test real(ode.p.solver) == Float32
+    @test real(ode.p.solver.basis) == Float32
+    @test real(ode.p.solver.mortar) == Float32
+    # TODO: `mesh` is currently not `adapt`ed correctly
+    @test real(ode.p.mesh) == Float64
+
+    @test ode.u0 isa CuArray
+    @test ode.p.solver.basis.derivative_matrix isa CuArray
+
+    @test Trixi.storage_type(ode.p.cache.elements) === CuArray
+    @test Trixi.storage_type(ode.p.cache.interfaces) === CuArray
+    @test Trixi.storage_type(ode.p.cache.boundaries) === CuArray
+    @test Trixi.storage_type(ode.p.cache.mortars) === CuArray
+end
+
 # Clean up afterwards: delete Trixi.jl output directory
 @test_nowarn isdir(outdir) && rm(outdir, recursive = true)
 end
diff --git a/test/test_cuda_3d.jl b/test/test_cuda_3d.jl
index 6c590332555..7a9174d79f4 100644
--- a/test/test_cuda_3d.jl
+++ b/test/test_cuda_3d.jl
@@ -67,6 +67,79 @@ end
     @test Trixi.storage_type(ode.p.cache.mortars) === CuArray
 end
 
+@trixi_testset "elixir_euler_source_terms.jl native" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms.jl"),
+                        # Expected errors are exactly the same as with TreeMesh!
+                        l2=[
+                            4.893619139889976e-5,
+                            5.3526950567182756e-5,
+                            5.35269505672133e-5,
+                            5.352695056735998e-5,
+                            0.00015172095200428318
+                        ],
+                        linf=[
+                            0.00031179856625374036,
+                            0.0003368725355339386,
+                            0.0003368725355383795,
+                            0.00033687253560787944,
+                            0.0013193387520935573
+                        ])
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
+    semi = ode.p # `semidiscretize` adapts the semi, so we need to obtain it from the ODE problem.
+    @test_allocations(Trixi.rhs!, semi, sol, 1000)
+    @test real(ode.p.solver) == Float64
+    @test real(ode.p.solver.basis) == Float64
+    @test real(ode.p.solver.mortar) == Float64
+    # TODO: `mesh` is currently not `adapt`ed correctly
+    @test real(ode.p.mesh) == Float64
+
+    @test ode.u0 isa Array
+    @test ode.p.solver.basis.derivative_matrix isa Array
+
+    @test Trixi.storage_type(ode.p.cache.elements) === Array
+    @test Trixi.storage_type(ode.p.cache.interfaces) === Array
+    @test Trixi.storage_type(ode.p.cache.boundaries) === Array
+    @test Trixi.storage_type(ode.p.cache.mortars) === Array
+end
+
+@trixi_testset "elixir_euler_source_terms.jl Float32 / CUDA" begin
+    # Using CUDA inside the testset since otherwise the bindings are hiddend by the anonymous modules
+    using CUDA
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms.jl"),
+                        # Expected errors are exactly the same as with TreeMesh!
+                        l2=Float32[4.912578089985958e-5,
+                                   5.3683407014580115e-5,
+                                   5.368099834769191e-5,
+                                   5.371664525206341e-5,
+                                   0.00015186256300882088],
+                        linf=Float32[0.00032772542853032327,
+                                     0.00035144807715092874,
+                                     0.0003549051465479014,
+                                     0.00035573961157475686,
+                                     0.0013591384887696734],
+                        RealT_for_test_tolerances=Float32,
+                        real_type=Float32,
+                        storage_type=CuArray, gamma=Float32(1.4))
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
+    semi = ode.p # `semidiscretize` adapts the semi, so we need to obtain it from the ODE problem.
+    @test_allocations(Trixi.rhs!, semi, sol, 100_000)
+    @test real(ode.p.solver) == Float32
+    @test real(ode.p.solver.basis) == Float32
+    @test real(ode.p.solver.mortar) == Float32
+    # TODO: `mesh` is currently not `adapt`ed correctly
+    @test real(ode.p.mesh) == Float64
+
+    @test ode.u0 isa CuArray
+    @test ode.p.solver.basis.derivative_matrix isa CuArray
+
+    @test Trixi.storage_type(ode.p.cache.elements) === CuArray
+    @test Trixi.storage_type(ode.p.cache.interfaces) === CuArray
+    @test Trixi.storage_type(ode.p.cache.boundaries) === CuArray
+    @test Trixi.storage_type(ode.p.cache.mortars) === CuArray
+end
+
 # Clean up afterwards: delete Trixi.jl output directory
 @test_nowarn isdir(outdir) && rm(outdir, recursive = true)
 end

From 20e890d1323924e2f8cdad02b27c9bc5e9315ead Mon Sep 17 00:00:00 2001
From: MarcoArtiano <artianomarco97@gmail.com>
Date: Tue, 19 May 2026 18:35:24 +0200
Subject: [PATCH 21/29] fix tests

---
 test/test_kernelabstractions.jl | 12 ++++++++----
 1 file changed, 8 insertions(+), 4 deletions(-)

diff --git a/test/test_kernelabstractions.jl b/test/test_kernelabstractions.jl
index 7a714b7a96b..dbe5cc7e585 100644
--- a/test/test_kernelabstractions.jl
+++ b/test/test_kernelabstractions.jl
@@ -46,7 +46,8 @@ end
 end
 
 @trixi_testset "elixir_euler_source_terms.jl native" begin
-    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms.jl"),
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "p4est_2d_dgsem",
+                                 "elixir_euler_source_terms.jl"),
                         # Expected errors are exactly the same as with TreeMesh!
                         l2=[9.321181254378498e-7,
                             1.418121074369651e-6,
@@ -65,7 +66,8 @@ end
 @trixi_testset "elixir_euler_source_terms.jl Float32 / AMDGPU" begin
     # Using AMDGPU inside the testset since otherwise the bindings are hiddend by the anonymous modules
     using AMDGPU
-    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms.jl"),
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "p4est_2d_dgsem",
+                                 "elixir_euler_source_terms.jl"),
                         # Expected errors are exactly the same as with TreeMesh!
                         l2=Float32[2.4917018095933837e-6,
                                    2.7148269885239423e-6,
@@ -115,7 +117,8 @@ end
 end
 
 @trixi_testset "elixir_euler_source_terms.jl native" begin
-    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms.jl"),
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "p4est_3d_dgsem",
+                                 "elixir_euler_source_terms.jl"),
                         # Expected errors are exactly the same as with TreeMesh!
                         l2=[
                             4.893619139889976e-5,
@@ -140,7 +143,8 @@ end
 @trixi_testset "elixir_euler_source_terms.jl Float32 / AMDGPU" begin
     # Using AMDGPU inside the testset since otherwise the bindings are hiddend by the anonymous modules
     using AMDGPU
-    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms.jl"),
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "p4est_3d_dgsem",
+                                 "elixir_euler_source_terms.jl"),
                         # Expected errors are exactly the same as with TreeMesh!
                         l2=Float32[4.912578089985958e-5,
                                    5.3683407014580115e-5,

From 031528dbdc5f5335f2fec03cbd53b7d832f21e77 Mon Sep 17 00:00:00 2001
From: MarcoArtiano <artianomarco97@gmail.com>
Date: Tue, 19 May 2026 21:21:34 +0200
Subject: [PATCH 22/29] fix allocations

---
 test/test_kernelabstractions.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/test/test_kernelabstractions.jl b/test/test_kernelabstractions.jl
index dbe5cc7e585..95bae4e2f37 100644
--- a/test/test_kernelabstractions.jl
+++ b/test/test_kernelabstractions.jl
@@ -60,7 +60,7 @@ end
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
     semi = ode.p # `semidiscretize` adapts the semi, so we need to obtain it from the ODE problem.
-    @test_allocations(Trixi.rhs!, semi, sol, 1000)
+    @test_allocations(Trixi.rhs!, semi, sol, 100_000)
 end
 
 @trixi_testset "elixir_euler_source_terms.jl Float32 / AMDGPU" begin
@@ -137,7 +137,7 @@ end
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
     semi = ode.p # `semidiscretize` adapts the semi, so we need to obtain it from the ODE problem.
-    @test_allocations(Trixi.rhs!, semi, sol, 1000)
+    @test_allocations(Trixi.rhs!, semi, sol, 100_000)
 end
 
 @trixi_testset "elixir_euler_source_terms.jl Float32 / AMDGPU" begin

From 6044e49b7d56c56bccb4c53c501cfd8ec2f32a71 Mon Sep 17 00:00:00 2001
From: MarcoArtiano <artianomarco97@gmail.com>
Date: Tue, 19 May 2026 21:44:04 +0200
Subject: [PATCH 23/29] fix allocations 3D

---
 test/test_kernelabstractions.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/test_kernelabstractions.jl b/test/test_kernelabstractions.jl
index 95bae4e2f37..6c3bbc6943c 100644
--- a/test/test_kernelabstractions.jl
+++ b/test/test_kernelabstractions.jl
@@ -137,7 +137,7 @@ end
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
     semi = ode.p # `semidiscretize` adapts the semi, so we need to obtain it from the ODE problem.
-    @test_allocations(Trixi.rhs!, semi, sol, 100_000)
+    @test_allocations(Trixi.rhs!, semi, sol, 400_000)
 end
 
 @trixi_testset "elixir_euler_source_terms.jl Float32 / AMDGPU" begin

From e19fe5f44a63f1547845ddc8d136c9cf951fbe84 Mon Sep 17 00:00:00 2001
From: MarcoArtiano <artianomarco97@gmail.com>
Date: Tue, 19 May 2026 21:59:49 +0200
Subject: [PATCH 24/29] fix allocations

---
 test/test_kernelabstractions.jl | 10 +++-------
 1 file changed, 3 insertions(+), 7 deletions(-)

diff --git a/test/test_kernelabstractions.jl b/test/test_kernelabstractions.jl
index 6c3bbc6943c..c8ceb266e26 100644
--- a/test/test_kernelabstractions.jl
+++ b/test/test_kernelabstractions.jl
@@ -63,9 +63,7 @@ end
     @test_allocations(Trixi.rhs!, semi, sol, 100_000)
 end
 
-@trixi_testset "elixir_euler_source_terms.jl Float32 / AMDGPU" begin
-    # Using AMDGPU inside the testset since otherwise the bindings are hiddend by the anonymous modules
-    using AMDGPU
+@trixi_testset "elixir_euler_source_terms.jl Float32" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "p4est_2d_dgsem",
                                  "elixir_euler_source_terms.jl"),
                         # Expected errors are exactly the same as with TreeMesh!
@@ -140,9 +138,7 @@ end
     @test_allocations(Trixi.rhs!, semi, sol, 400_000)
 end
 
-@trixi_testset "elixir_euler_source_terms.jl Float32 / AMDGPU" begin
-    # Using AMDGPU inside the testset since otherwise the bindings are hiddend by the anonymous modules
-    using AMDGPU
+@trixi_testset "elixir_euler_source_terms.jl Float32" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "p4est_3d_dgsem",
                                  "elixir_euler_source_terms.jl"),
                         # Expected errors are exactly the same as with TreeMesh!
@@ -162,7 +158,7 @@ end
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
     semi = ode.p # `semidiscretize` adapts the semi, so we need to obtain it from the ODE problem.
-    @test_allocations(Trixi.rhs!, semi, sol, 100_000)
+    @test_allocations(Trixi.rhs!, semi, sol, 400_000)
 end
 end
 

From f3e650d82f38927610cca72a6b96a453a0cf6725 Mon Sep 17 00:00:00 2001
From: Marco Artiano <57838732+MarcoArtiano@users.noreply.github.com>
Date: Thu, 21 May 2026 14:01:29 +0200
Subject: [PATCH 25/29] Apply suggestions from code review

Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com>
---
 test/test_amdgpu_2d.jl          | 1 -
 test/test_amdgpu_3d.jl          | 1 -
 test/test_cuda_2d.jl            | 1 -
 test/test_cuda_3d.jl            | 1 -
 test/test_kernelabstractions.jl | 2 --
 5 files changed, 6 deletions(-)

diff --git a/test/test_amdgpu_2d.jl b/test/test_amdgpu_2d.jl
index 4738a208a58..595a7781f54 100644
--- a/test/test_amdgpu_2d.jl
+++ b/test/test_amdgpu_2d.jl
@@ -101,7 +101,6 @@ end
     # Using AMDGPU inside the testset since otherwise the bindings are hiddend by the anonymous modules
     using AMDGPU
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms.jl"),
-                        # Expected errors are exactly the same as with TreeMesh!
                         l2=Float32[2.4917018095933837e-6,
                                    2.7148269885239423e-6,
                                    2.695290306860358e-6,
diff --git a/test/test_amdgpu_3d.jl b/test/test_amdgpu_3d.jl
index 828c2f5efbe..8a35178dbf8 100644
--- a/test/test_amdgpu_3d.jl
+++ b/test/test_amdgpu_3d.jl
@@ -107,7 +107,6 @@ end
     # Using AMDGPU inside the testset since otherwise the bindings are hiddend by the anonymous modules
     using AMDGPU
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms.jl"),
-                        # Expected errors are exactly the same as with TreeMesh!
                         l2=Float32[4.912578089985958e-5,
                                    5.3683407014580115e-5,
                                    5.368099834769191e-5,
diff --git a/test/test_cuda_2d.jl b/test/test_cuda_2d.jl
index 02fb12a0ba0..1623fb8f961 100644
--- a/test/test_cuda_2d.jl
+++ b/test/test_cuda_2d.jl
@@ -101,7 +101,6 @@ end
     # Using CUDA inside the testset since otherwise the bindings are hiddend by the anonymous modules
     using CUDA
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms.jl"),
-                        # Expected errors are exactly the same as with TreeMesh!
                         l2=Float32[2.4917018095933837e-6,
                                    2.7148269885239423e-6,
                                    2.695290306860358e-6,
diff --git a/test/test_cuda_3d.jl b/test/test_cuda_3d.jl
index 7a9174d79f4..ee1a68e854a 100644
--- a/test/test_cuda_3d.jl
+++ b/test/test_cuda_3d.jl
@@ -107,7 +107,6 @@ end
     # Using CUDA inside the testset since otherwise the bindings are hiddend by the anonymous modules
     using CUDA
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms.jl"),
-                        # Expected errors are exactly the same as with TreeMesh!
                         l2=Float32[4.912578089985958e-5,
                                    5.3683407014580115e-5,
                                    5.368099834769191e-5,
diff --git a/test/test_kernelabstractions.jl b/test/test_kernelabstractions.jl
index c8ceb266e26..5adf971ef70 100644
--- a/test/test_kernelabstractions.jl
+++ b/test/test_kernelabstractions.jl
@@ -66,7 +66,6 @@ end
 @trixi_testset "elixir_euler_source_terms.jl Float32" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "p4est_2d_dgsem",
                                  "elixir_euler_source_terms.jl"),
-                        # Expected errors are exactly the same as with TreeMesh!
                         l2=Float32[2.4917018095933837e-6,
                                    2.7148269885239423e-6,
                                    2.695290306860358e-6,
@@ -141,7 +140,6 @@ end
 @trixi_testset "elixir_euler_source_terms.jl Float32" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "p4est_3d_dgsem",
                                  "elixir_euler_source_terms.jl"),
-                        # Expected errors are exactly the same as with TreeMesh!
                         l2=Float32[4.912578089985958e-5,
                                    5.3683407014580115e-5,
                                    5.368099834769191e-5,

From b15439e05627d706d61fc84eb9bac235123878af Mon Sep 17 00:00:00 2001
From: Marco Artiano <57838732+MarcoArtiano@users.noreply.github.com>
Date: Thu, 21 May 2026 14:02:05 +0200
Subject: [PATCH 26/29] Update test/test_amdgpu_3d.jl

---
 test/test_amdgpu_3d.jl | 1 -
 1 file changed, 1 deletion(-)

diff --git a/test/test_amdgpu_3d.jl b/test/test_amdgpu_3d.jl
index 8a35178dbf8..5802aa3ae68 100644
--- a/test/test_amdgpu_3d.jl
+++ b/test/test_amdgpu_3d.jl
@@ -69,7 +69,6 @@ end
 
 @trixi_testset "elixir_euler_source_terms.jl native" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms.jl"),
-                        # Expected errors are exactly the same as with TreeMesh!
                         l2=[
                             4.893619139889976e-5,
                             5.3526950567182756e-5,

From 7f62e6e76ce122b92a65b4d4aef8b4055e39e711 Mon Sep 17 00:00:00 2001
From: Marco Artiano <57838732+MarcoArtiano@users.noreply.github.com>
Date: Thu, 21 May 2026 14:27:10 +0200
Subject: [PATCH 27/29] Apply suggestions from code review

Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com>
---
 test/test_cuda_3d.jl            | 1 -
 test/test_kernelabstractions.jl | 1 -
 2 files changed, 2 deletions(-)

diff --git a/test/test_cuda_3d.jl b/test/test_cuda_3d.jl
index ee1a68e854a..1d573b11881 100644
--- a/test/test_cuda_3d.jl
+++ b/test/test_cuda_3d.jl
@@ -69,7 +69,6 @@ end
 
 @trixi_testset "elixir_euler_source_terms.jl native" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms.jl"),
-                        # Expected errors are exactly the same as with TreeMesh!
                         l2=[
                             4.893619139889976e-5,
                             5.3526950567182756e-5,
diff --git a/test/test_kernelabstractions.jl b/test/test_kernelabstractions.jl
index 5adf971ef70..5887cab4aca 100644
--- a/test/test_kernelabstractions.jl
+++ b/test/test_kernelabstractions.jl
@@ -116,7 +116,6 @@ end
 @trixi_testset "elixir_euler_source_terms.jl native" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "p4est_3d_dgsem",
                                  "elixir_euler_source_terms.jl"),
-                        # Expected errors are exactly the same as with TreeMesh!
                         l2=[
                             4.893619139889976e-5,
                             5.3526950567182756e-5,

From d4f6ff9179f565b95d109149631c90e2ac8afb73 Mon Sep 17 00:00:00 2001
From: Marco Artiano <57838732+MarcoArtiano@users.noreply.github.com>
Date: Thu, 21 May 2026 16:59:23 +0200
Subject: [PATCH 28/29] add comments about storage_type and real_type choices

Co-authored-by: Marco Artiano <57838732+MarcoArtiano@users.noreply.github.com>
---
 examples/p4est_2d_dgsem/elixir_euler_source_terms.jl | 4 ++++
 examples/p4est_3d_dgsem/elixir_euler_source_terms.jl | 4 ++++
 2 files changed, 8 insertions(+)

diff --git a/examples/p4est_2d_dgsem/elixir_euler_source_terms.jl b/examples/p4est_2d_dgsem/elixir_euler_source_terms.jl
index cdc8d2bcc86..f85658c5423 100644
--- a/examples/p4est_2d_dgsem/elixir_euler_source_terms.jl
+++ b/examples/p4est_2d_dgsem/elixir_euler_source_terms.jl
@@ -37,6 +37,10 @@ semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver;
 # ODE solvers, callbacks etc.
 
 tspan = (0.0, 2.0)
+# Create ODE problem with time span from 0.0 to 1.0
+# Setting `real_type` allows to change the real number type, e.g., to `Float32`.
+# This is particularly useful when changing the `storage_type` to a GPU array
+# type such as `ROCArray` (AMD) or `CuArray` (NVIDIA CUDA).
 ode = semidiscretize(semi, tspan; real_type = nothing, storage_type = nothing)
 
 summary_callback = SummaryCallback()
diff --git a/examples/p4est_3d_dgsem/elixir_euler_source_terms.jl b/examples/p4est_3d_dgsem/elixir_euler_source_terms.jl
index 3e1aee8be47..3fce3e54fb2 100644
--- a/examples/p4est_3d_dgsem/elixir_euler_source_terms.jl
+++ b/examples/p4est_3d_dgsem/elixir_euler_source_terms.jl
@@ -39,6 +39,10 @@ semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver;
 # ODE solvers, callbacks etc.
 
 tspan = (0.0, 5.0)
+# Create ODE problem with time span from 0.0 to 1.0
+# Setting `real_type` allows to change the real number type, e.g., to `Float32`.
+# This is particularly useful when changing the `storage_type` to a GPU array
+# type such as `ROCArray` (AMD) or `CuArray` (NVIDIA CUDA).
 ode = semidiscretize(semi, tspan; real_type = nothing, storage_type = nothing)
 
 summary_callback = SummaryCallback()

From 34215b6167cc9be658d35b99be96c72f19525ef3 Mon Sep 17 00:00:00 2001
From: Marco Artiano <57838732+MarcoArtiano@users.noreply.github.com>
Date: Thu, 21 May 2026 17:08:23 +0200
Subject: [PATCH 29/29] Apply suggestions from code review

Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com>
---
 test/test_amdgpu_2d.jl          | 5 ++++-
 test/test_amdgpu_3d.jl          | 5 ++++-
 test/test_cuda_2d.jl            | 5 ++++-
 test/test_cuda_3d.jl            | 5 ++++-
 test/test_kernelabstractions.jl | 4 ++--
 5 files changed, 18 insertions(+), 6 deletions(-)

diff --git a/test/test_amdgpu_2d.jl b/test/test_amdgpu_2d.jl
index 595a7781f54..145786eac24 100644
--- a/test/test_amdgpu_2d.jl
+++ b/test/test_amdgpu_2d.jl
@@ -87,6 +87,7 @@ end
     @test real(ode.p.solver.mortar) == Float64
     # TODO: `mesh` is currently not `adapt`ed correctly
     @test real(ode.p.mesh) == Float64
+    @test typeof(equations.gamma) == Float64
 
     @test ode.u0 isa Array
     @test ode.p.solver.basis.derivative_matrix isa Array
@@ -111,7 +112,8 @@ end
                                      6.214141845717336e-5],
                         RealT_for_test_tolerances=Float32,
                         real_type=Float32,
-                        storage_type=ROCArray, gamma=Float32(1.4))
+                        storage_type=ROCArray,
+                        gamma=Float32(1.4)) # TODO: This should not be required
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
     semi = ode.p # `semidiscretize` adapts the semi, so we need to obtain it from the ODE problem.
@@ -121,6 +123,7 @@ end
     @test real(ode.p.solver.mortar) == Float32
     # TODO: `mesh` is currently not `adapt`ed correctly
     @test real(ode.p.mesh) == Float64
+    @test typeof(equations.gamma) == Float32
 
     @test ode.u0 isa ROCArray
     @test ode.p.solver.basis.derivative_matrix isa ROCArray
diff --git a/test/test_amdgpu_3d.jl b/test/test_amdgpu_3d.jl
index 5802aa3ae68..328861732a2 100644
--- a/test/test_amdgpu_3d.jl
+++ b/test/test_amdgpu_3d.jl
@@ -92,6 +92,7 @@ end
     @test real(ode.p.solver.mortar) == Float64
     # TODO: `mesh` is currently not `adapt`ed correctly
     @test real(ode.p.mesh) == Float64
+    @test typeof(equations.gamma) == Float64
 
     @test ode.u0 isa Array
     @test ode.p.solver.basis.derivative_matrix isa Array
@@ -118,7 +119,8 @@ end
                                      0.0013591384887696734],
                         RealT_for_test_tolerances=Float32,
                         real_type=Float32,
-                        storage_type=ROCArray, gamma=Float32(1.4))
+                        storage_type=ROCArray,
+                        gamma=Float32(1.4)) # TODO: This should not be required
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
     semi = ode.p # `semidiscretize` adapts the semi, so we need to obtain it from the ODE problem.
@@ -128,6 +130,7 @@ end
     @test real(ode.p.solver.mortar) == Float32
     # TODO: `mesh` is currently not `adapt`ed correctly
     @test real(ode.p.mesh) == Float64
+    @test typeof(equations.gamma) == Float32
 
     @test ode.u0 isa ROCArray
     @test ode.p.solver.basis.derivative_matrix isa ROCArray
diff --git a/test/test_cuda_2d.jl b/test/test_cuda_2d.jl
index 1623fb8f961..70c822814c6 100644
--- a/test/test_cuda_2d.jl
+++ b/test/test_cuda_2d.jl
@@ -87,6 +87,7 @@ end
     @test real(ode.p.solver.mortar) == Float64
     # TODO: `mesh` is currently not `adapt`ed correctly
     @test real(ode.p.mesh) == Float64
+    @test typeof(equations.gamma) == Float64
 
     @test ode.u0 isa Array
     @test ode.p.solver.basis.derivative_matrix isa Array
@@ -111,7 +112,8 @@ end
                                      6.214141845717336e-5],
                         RealT_for_test_tolerances=Float32,
                         real_type=Float32,
-                        storage_type=CuArray, gamma=Float32(1.4))
+                        storage_type=CuArray,
+                        gamma=Float32(1.4)) # TODO: This should not be required
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
     semi = ode.p # `semidiscretize` adapts the semi, so we need to obtain it from the ODE problem.
@@ -121,6 +123,7 @@ end
     @test real(ode.p.solver.mortar) == Float32
     # TODO: `mesh` is currently not `adapt`ed correctly
     @test real(ode.p.mesh) == Float64
+    @test typeof(equations.gamma) == Float32
 
     @test ode.u0 isa CuArray
     @test ode.p.solver.basis.derivative_matrix isa CuArray
diff --git a/test/test_cuda_3d.jl b/test/test_cuda_3d.jl
index 1d573b11881..95809d22146 100644
--- a/test/test_cuda_3d.jl
+++ b/test/test_cuda_3d.jl
@@ -92,6 +92,7 @@ end
     @test real(ode.p.solver.mortar) == Float64
     # TODO: `mesh` is currently not `adapt`ed correctly
     @test real(ode.p.mesh) == Float64
+    @test typeof(equations.gamma) == Float64
 
     @test ode.u0 isa Array
     @test ode.p.solver.basis.derivative_matrix isa Array
@@ -118,7 +119,8 @@ end
                                      0.0013591384887696734],
                         RealT_for_test_tolerances=Float32,
                         real_type=Float32,
-                        storage_type=CuArray, gamma=Float32(1.4))
+                        storage_type=CuArray,
+                        gamma=Float32(1.4)) # TODO: This should not be required
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
     semi = ode.p # `semidiscretize` adapts the semi, so we need to obtain it from the ODE problem.
@@ -128,6 +130,7 @@ end
     @test real(ode.p.solver.mortar) == Float32
     # TODO: `mesh` is currently not `adapt`ed correctly
     @test real(ode.p.mesh) == Float64
+    @test typeof(equations.gamma) == Float64
 
     @test ode.u0 isa CuArray
     @test ode.p.solver.basis.derivative_matrix isa CuArray
diff --git a/test/test_kernelabstractions.jl b/test/test_kernelabstractions.jl
index 5887cab4aca..1e39e058131 100644
--- a/test/test_kernelabstractions.jl
+++ b/test/test_kernelabstractions.jl
@@ -76,7 +76,7 @@ end
                                      6.214141845717336e-5],
                         RealT_for_test_tolerances=Float32,
                         real_type=Float32,
-                        gamma=Float32(1.4))
+                        gamma=Float32(1.4)) # TODO: This should not be required
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
     semi = ode.p # `semidiscretize` adapts the semi, so we need to obtain it from the ODE problem.
@@ -151,7 +151,7 @@ end
                                      0.0013591384887696734],
                         RealT_for_test_tolerances=Float32,
                         real_type=Float32,
-                        gamma=Float32(1.4))
+                        gamma=Float32(1.4)) # TODO: This should not be required
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
     semi = ode.p # `semidiscretize` adapts the semi, so we need to obtain it from the ODE problem.