diff --git a/src/callbacks_step/stepsize_dg1d.jl b/src/callbacks_step/stepsize_dg1d.jl
index edc25ec78f..cdb7ed21de 100644
--- a/src/callbacks_step/stepsize_dg1d.jl
+++ b/src/callbacks_step/stepsize_dg1d.jl
@@ -11,7 +11,7 @@ function max_dt(u, t, mesh::TreeMesh{1},
     # e.g. for steady-state linear advection
     max_scaled_speed = nextfloat(zero(t))
 
-    for element in eachelement(dg, cache)
+    @batch reduction=(max, max_scaled_speed) for element in eachelement(dg, cache)
         max_lambda1 = zero(max_scaled_speed)
         for i in eachnode(dg)
             u_node = get_node_vars(u, equations, dg, i, element)
@@ -31,7 +31,7 @@ function max_dt(u, t, mesh::TreeMesh{1},
     # e.g. for steady-state linear advection
     max_scaled_speed = nextfloat(zero(t))
 
-    for element in eachelement(dg, cache)
+    @batch reduction=(max, max_scaled_speed) for element in eachelement(dg, cache)
         max_lambda1, = max_abs_speeds(equations)
         inv_jacobian = cache.elements.inverse_jacobian[element]
         max_scaled_speed = max(max_scaled_speed, inv_jacobian * max_lambda1)
@@ -46,7 +46,7 @@ function max_dt(u, t, mesh::StructuredMesh{1},
     # e.g. for steady-state linear advection
     max_scaled_speed = nextfloat(zero(t))
 
-    for element in eachelement(dg, cache)
+    @batch reduction=(max, max_scaled_speed) for element in eachelement(dg, cache)
         max_lambda1 = zero(max_scaled_speed)
 
         for i in eachnode(dg)
@@ -70,7 +70,7 @@ function max_dt(u, t, mesh::StructuredMesh{1},
     # e.g. for steady-state linear advection
     max_scaled_speed = nextfloat(zero(t))
 
-    for element in eachelement(dg, cache)
+    @batch reduction=(max, max_scaled_speed) for element in eachelement(dg, cache)
         max_lambda1, = max_abs_speeds(equations)
 
         for i in eachnode(dg)
diff --git a/src/callbacks_step/stepsize_dg2d.jl b/src/callbacks_step/stepsize_dg2d.jl
index c7922cecc6..b895dcd5b5 100644
--- a/src/callbacks_step/stepsize_dg2d.jl
+++ b/src/callbacks_step/stepsize_dg2d.jl
@@ -11,7 +11,7 @@ function max_dt(u, t, mesh::TreeMesh{2},
     # e.g. for steady-state linear advection
     max_scaled_speed = nextfloat(zero(t))
 
-    for element in eachelement(dg, cache)
+    @batch reduction=(max, max_scaled_speed) for element in eachelement(dg, cache)
         max_lambda1 = max_lambda2 = zero(max_scaled_speed)
         for j in eachnode(dg), i in eachnode(dg)
             u_node = get_node_vars(u, equations, dg, i, j, element)
@@ -33,7 +33,7 @@ function max_dt(u, t, mesh::TreeMesh{2},
     # e.g. for steady-state linear advection
     max_scaled_speed = nextfloat(zero(t))
 
-    for element in eachelement(dg, cache)
+    @batch reduction=(max, max_scaled_speed) for element in eachelement(dg, cache)
         max_lambda1, max_lambda2 = max_abs_speeds(equations)
         inv_jacobian = cache.elements.inverse_jacobian[element]
         max_scaled_speed = max(max_scaled_speed,
@@ -87,7 +87,7 @@ function max_dt(u, t,
 
     @unpack contravariant_vectors, inverse_jacobian = cache.elements
 
-    for element in eachelement(dg, cache)
+    @batch reduction=(max, max_scaled_speed) for element in eachelement(dg, cache)
         max_lambda1 = max_lambda2 = zero(max_scaled_speed)
         for j in eachnode(dg), i in eachnode(dg)
             u_node = get_node_vars(u, equations, dg, i, j, element)
@@ -125,7 +125,7 @@ function max_dt(u, t,
 
     max_lambda1, max_lambda2 = max_abs_speeds(equations)
 
-    for element in eachelement(dg, cache)
+    @batch reduction=(max, max_scaled_speed) for element in eachelement(dg, cache)
         for j in eachnode(dg), i in eachnode(dg)
             # Local speeds transformed to the reference element
             Ja11, Ja12 = get_contravariant_vector(1, contravariant_vectors, i, j,
diff --git a/src/callbacks_step/stepsize_dg3d.jl b/src/callbacks_step/stepsize_dg3d.jl
index 49976de650..3324e819ce 100644
--- a/src/callbacks_step/stepsize_dg3d.jl
+++ b/src/callbacks_step/stepsize_dg3d.jl
@@ -11,7 +11,7 @@ function max_dt(u, t, mesh::TreeMesh{3},
     # e.g. for steady-state linear advection
     max_scaled_speed = nextfloat(zero(t))
 
-    for element in eachelement(dg, cache)
+    @batch reduction=(max, max_scaled_speed) 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)
@@ -34,7 +34,7 @@ function max_dt(u, t, mesh::TreeMesh{3},
     # e.g. for steady-state linear advection
     max_scaled_speed = nextfloat(zero(t))
 
-    for element in eachelement(dg, cache)
+    @batch reduction=(max, max_scaled_speed) for element in eachelement(dg, cache)
         max_lambda1, max_lambda2, max_lambda3 = max_abs_speeds(equations)
         inv_jacobian = cache.elements.inverse_jacobian[element]
         max_scaled_speed = max(max_scaled_speed,
@@ -52,7 +52,7 @@ function max_dt(u, t, mesh::Union{StructuredMesh{3}, P4estMesh{3}, T8codeMesh{3}
 
     @unpack contravariant_vectors = cache.elements
 
-    for element in eachelement(dg, cache)
+    @batch reduction=(max, max_scaled_speed) 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)
@@ -92,7 +92,7 @@ function max_dt(u, t, mesh::Union{StructuredMesh{3}, P4estMesh{3}, T8codeMesh{3}
 
     max_lambda1, max_lambda2, max_lambda3 = max_abs_speeds(equations)
 
-    for element in eachelement(dg, cache)
+    @batch reduction=(max, max_scaled_speed) for element in eachelement(dg, cache)
         for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
             Ja11, Ja12, Ja13 = get_contravariant_vector(1, contravariant_vectors, i, j,
                                                         k, element)