diff --git a/examples/p4est_3d_dgsem/elixir_advection_basic_fd.jl b/examples/p4est_3d_dgsem/elixir_advection_basic_fd.jl
index e4ebe1692b..c38a4edd75 100644
--- a/examples/p4est_3d_dgsem/elixir_advection_basic_fd.jl
+++ b/examples/p4est_3d_dgsem/elixir_advection_basic_fd.jl
@@ -12,7 +12,7 @@ equations = LinearScalarAdvectionEquation3D(advection_velocity)
 
 # Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux
 solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs,
-               volume_integral=VolumeIntegralFluxDifferencing(flux_lax_friedrichs), backend=backendcl)
+               volume_integral=VolumeIntegralFluxDifferencing(flux_lax_friedrichs), backend=backend)
 
 coordinates_min = (-1.0, -1.0, -1.0) # minimum coordinates (min(x), min(y), min(z))
 coordinates_max = ( 1.0,  1.0,  1.0) # maximum coordinates (max(x), max(y), max(z))
diff --git a/examples/p4est_3d_dgsem/elixir_euler_taylor_green_vortex.jl b/examples/p4est_3d_dgsem/elixir_euler_taylor_green_vortex.jl
index 182b60ae41..b8fba2d5b8 100644
--- a/examples/p4est_3d_dgsem/elixir_euler_taylor_green_vortex.jl
+++ b/examples/p4est_3d_dgsem/elixir_euler_taylor_green_vortex.jl
@@ -1,5 +1,6 @@
 using OrdinaryDiffEq
 using Trixi
+using KernelAbstractions
 
 ###############################################################################
 # semidiscretization of the compressible Euler equations
@@ -24,10 +25,13 @@ function initial_condition_taylor_green_vortex(x, t, equations::CompressibleEule
 
   return prim2cons(SVector(rho, v1, v2, v3, p), equations)
 end
+
+backend = CPU()
+
 initial_condition = initial_condition_taylor_green_vortex
 
 solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs,
-               volume_integral=VolumeIntegralFluxDifferencing(flux_lax_friedrichs))
+               volume_integral=VolumeIntegralFluxDifferencing(flux_lax_friedrichs), backend=backend)
 
 coordinates_min = (-1.0, -1.0, -1.0) .* pi
 coordinates_max = ( 1.0,  1.0,  1.0) .* pi
@@ -38,14 +42,14 @@ mesh = P4estMesh(trees_per_dimension, polydeg=1,
                  coordinates_min=coordinates_min, coordinates_max=coordinates_max,
                  initial_refinement_level=1)
 
-semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver; backend=backend)
 
 
 ###############################################################################
 # ODE solvers, callbacks etc.
 
 tspan = (0.0, 5.0)
-ode = semidiscretize(semi, tspan)
+ode = semidiscretize(semi, tspan; offload=true, backend=backend)
 
 summary_callback = SummaryCallback()
 
diff --git a/src/callbacks_step/stepsize_dg3d.jl b/src/callbacks_step/stepsize_dg3d.jl
index 9079974caa..5ef7010e33 100644
--- a/src/callbacks_step/stepsize_dg3d.jl
+++ b/src/callbacks_step/stepsize_dg3d.jl
@@ -50,25 +50,35 @@ function max_dt(u, t, mesh::Union{StructuredMesh{3}, P4estMesh{3}},
     # e.g. for steady-state linear advection
     max_scaled_speed = nextfloat(zero(t))
 
-    @unpack contravariant_vectors = cache.elements
+    @unpack contravariant_vectors, inverse_jacobian = cache.elements
+
+    tmp_u = copyto!(CPU(),
+                    allocate(CPU(), eltype(u), size(u)),
+                    u)
+    tmp_contravariant_vectors = copyto!(CPU(),
+                                   allocate(CPU(), eltype(contravariant_vectors), size(contravariant_vectors)),
+                                   contravariant_vectors)
+    tmp_inverse_jacobian = copyto!(CPU(),
+                                   allocate(CPU(), eltype(inverse_jacobian), size(inverse_jacobian)),
+                                   inverse_jacobian)
 
     for element in eachelement(dg, cache)
         max_lambda1 = max_lambda2 = max_lambda3 = zero(max_scaled_speed)
         for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
-            u_node = get_node_vars(u, equations, dg, i, j, k, element)
+            u_node = get_node_vars(tmp_u, equations, dg, i, j, k, element)
             lambda1, lambda2, lambda3 = max_abs_speeds(u_node, equations)
 
-            Ja11, Ja12, Ja13 = get_contravariant_vector(1, contravariant_vectors, i, j,
+            Ja11, Ja12, Ja13 = get_contravariant_vector(1, tmp_contravariant_vectors, i, j,
                                                         k, element)
             lambda1_transformed = abs(Ja11 * lambda1 + Ja12 * lambda2 + Ja13 * lambda3)
-            Ja21, Ja22, Ja23 = get_contravariant_vector(2, contravariant_vectors, i, j,
+            Ja21, Ja22, Ja23 = get_contravariant_vector(2, tmp_contravariant_vectors, i, j,
                                                         k, element)
             lambda2_transformed = abs(Ja21 * lambda1 + Ja22 * lambda2 + Ja23 * lambda3)
-            Ja31, Ja32, Ja33 = get_contravariant_vector(3, contravariant_vectors, i, j,
+            Ja31, Ja32, Ja33 = get_contravariant_vector(3, tmp_contravariant_vectors, i, j,
                                                         k, element)
             lambda3_transformed = abs(Ja31 * lambda1 + Ja32 * lambda2 + Ja33 * lambda3)
 
-            inv_jacobian = abs(cache.elements.inverse_jacobian[i, j, k, element])
+            inv_jacobian = abs(tmp_inverse_jacobian[i, j, k, element])
 
             max_lambda1 = max(max_lambda1, inv_jacobian * lambda1_transformed)
             max_lambda2 = max(max_lambda2, inv_jacobian * lambda2_transformed)