diff --git a/src/callbacks_step/analysis_dg1d.jl b/src/callbacks_step/analysis_dg1d.jl
index 82064c7cb8..444d218f11 100644
--- a/src/callbacks_step/analysis_dg1d.jl
+++ b/src/callbacks_step/analysis_dg1d.jl
@@ -40,16 +40,15 @@ function create_cache_analysis(analyzer, mesh::StructuredMesh{1},
 end
 
 function calc_error_norms(func, u, t, analyzer,
-                          mesh::StructuredMesh{1}, equations, initial_condition,
+                          mesh::TreeMesh{1}, equations, initial_condition,
                           dg::DGSEM, cache, cache_analysis)
     @unpack vandermonde, weights = analyzer
-    @unpack node_coordinates, inverse_jacobian = cache.elements
-    @unpack u_local, x_local, jacobian_local = cache_analysis
+    @unpack node_coordinates = cache.elements
+    @unpack u_local, x_local = cache_analysis
 
     # Set up data structures
     l2_error = zero(func(get_node_vars(u, equations, dg, 1, 1), equations))
     linf_error = copy(l2_error)
-    total_volume = zero(real(mesh))
 
     # Iterate over all elements for error calculations
     for element in eachelement(dg, cache)
@@ -57,40 +56,38 @@ function calc_error_norms(func, u, t, analyzer,
         multiply_dimensionwise!(u_local, vandermonde, view(u, :, :, element))
         multiply_dimensionwise!(x_local, vandermonde,
                                 view(node_coordinates, :, :, element))
-        multiply_scalar_dimensionwise!(jacobian_local, vandermonde,
-                                       inv.(view(inverse_jacobian, :, element)))
 
         # Calculate errors at each analysis node
+        volume_jacobian_ = volume_jacobian(element, mesh, cache)
+
         for i in eachnode(analyzer)
             u_exact = initial_condition(get_node_coords(x_local, equations, dg, i), t,
                                         equations)
             diff = func(u_exact, equations) -
                    func(get_node_vars(u_local, equations, dg, i), equations)
-            # We take absolute value as we need the Jacobian here for the volume calculation
-            abs_jacobian_local_i = abs(jacobian_local[i])
-
-            l2_error += diff .^ 2 * (weights[i] * abs_jacobian_local_i)
+            l2_error += diff .^ 2 * (weights[i] * volume_jacobian_)
             linf_error = @. max(linf_error, abs(diff))
-            total_volume += weights[i] * abs_jacobian_local_i
         end
     end
 
     # For L2 error, divide by total volume
-    l2_error = @. sqrt(l2_error / total_volume)
+    total_volume_ = total_volume(mesh)
+    l2_error = @. sqrt(l2_error / total_volume_)
 
     return l2_error, linf_error
 end
 
 function calc_error_norms(func, u, t, analyzer,
-                          mesh::TreeMesh{1}, equations, initial_condition,
+                          mesh::StructuredMesh{1}, equations, initial_condition,
                           dg::DGSEM, cache, cache_analysis)
     @unpack vandermonde, weights = analyzer
-    @unpack node_coordinates = cache.elements
-    @unpack u_local, x_local = cache_analysis
+    @unpack node_coordinates, inverse_jacobian = cache.elements
+    @unpack u_local, x_local, jacobian_local = cache_analysis
 
     # Set up data structures
     l2_error = zero(func(get_node_vars(u, equations, dg, 1, 1), equations))
     linf_error = copy(l2_error)
+    total_volume = zero(real(mesh))
 
     # Iterate over all elements for error calculations
     for element in eachelement(dg, cache)
@@ -98,23 +95,26 @@ function calc_error_norms(func, u, t, analyzer,
         multiply_dimensionwise!(u_local, vandermonde, view(u, :, :, element))
         multiply_dimensionwise!(x_local, vandermonde,
                                 view(node_coordinates, :, :, element))
+        multiply_scalar_dimensionwise!(jacobian_local, vandermonde,
+                                       inv.(view(inverse_jacobian, :, element)))
 
         # Calculate errors at each analysis node
-        volume_jacobian_ = volume_jacobian(element, mesh, cache)
-
         for i in eachnode(analyzer)
             u_exact = initial_condition(get_node_coords(x_local, equations, dg, i), t,
                                         equations)
             diff = func(u_exact, equations) -
                    func(get_node_vars(u_local, equations, dg, i), equations)
-            l2_error += diff .^ 2 * (weights[i] * volume_jacobian_)
+            # We take absolute value as we need the Jacobian here for the volume calculation
+            abs_jacobian_local_i = abs(jacobian_local[i])
+
+            l2_error += diff .^ 2 * (weights[i] * abs_jacobian_local_i)
             linf_error = @. max(linf_error, abs(diff))
+            total_volume += weights[i] * abs_jacobian_local_i
         end
     end
 
     # For L2 error, divide by total volume
-    total_volume_ = total_volume(mesh)
-    l2_error = @. sqrt(l2_error / total_volume_)
+    l2_error = @. sqrt(l2_error / total_volume)
 
     return l2_error, linf_error
 end
@@ -214,7 +214,7 @@ function analyze(::Val{:linf_divb}, du, u, t,
 
     # integrate over all elements to get the divergence-free condition errors
     linf_divb = zero(eltype(u))
-    for element in eachelement(dg, cache)
+    @batch reduction=(max, linf_divb) for element in eachelement(dg, cache)
         for i in eachnode(dg)
             divb = zero(eltype(u))
             for k in eachnode(dg)
diff --git a/src/callbacks_step/analysis_dg2d.jl b/src/callbacks_step/analysis_dg2d.jl
index d8da0ec7f7..658778828b 100644
--- a/src/callbacks_step/analysis_dg2d.jl
+++ b/src/callbacks_step/analysis_dg2d.jl
@@ -338,7 +338,7 @@ function analyze(::Val{:linf_divb}, du, u, t,
 
     # integrate over all elements to get the divergence-free condition errors
     linf_divb = zero(eltype(u))
-    for element in eachelement(dg, cache)
+    @batch reduction=(max, linf_divb) for element in eachelement(dg, cache)
         for j in eachnode(dg), i in eachnode(dg)
             divb = zero(eltype(u))
             for k in eachnode(dg)
@@ -368,7 +368,7 @@ function analyze(::Val{:linf_divb}, du, u, t,
 
     # integrate over all elements to get the divergence-free condition errors
     linf_divb = zero(eltype(u))
-    for element in eachelement(dg, cache)
+    @batch reduction=(max, linf_divb) for element in eachelement(dg, cache)
         for j in eachnode(dg), i in eachnode(dg)
             divb = zero(eltype(u))
             # Get the contravariant vectors Ja^1 and Ja^2
diff --git a/src/callbacks_step/analysis_dg3d.jl b/src/callbacks_step/analysis_dg3d.jl
index 8bb6c63016..c5011d558d 100644
--- a/src/callbacks_step/analysis_dg3d.jl
+++ b/src/callbacks_step/analysis_dg3d.jl
@@ -378,7 +378,7 @@ function analyze(::Val{:linf_divb}, du, u, t,
 
     # integrate over all elements to get the divergence-free condition errors
     linf_divb = zero(eltype(u))
-    for element in eachelement(dg, cache)
+    @batch reduction=(max, linf_divb) for element in eachelement(dg, cache)
         for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
             divb = zero(eltype(u))
             for l in eachnode(dg)
@@ -416,7 +416,7 @@ function analyze(::Val{:linf_divb}, du, u, t,
 
     # integrate over all elements to get the divergence-free condition errors
     linf_divb = zero(eltype(u))
-    for element in eachelement(dg, cache)
+    @batch reduction=(max, linf_divb) for element in eachelement(dg, cache)
         for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
             divb = zero(eltype(u))
             # Get the contravariant vectors Ja^1, Ja^2, and Ja^3
diff --git a/src/callbacks_step/analysis_dgmulti.jl b/src/callbacks_step/analysis_dgmulti.jl
index 8d43bfa90e..b3223340ae 100644
--- a/src/callbacks_step/analysis_dgmulti.jl
+++ b/src/callbacks_step/analysis_dgmulti.jl
@@ -125,7 +125,7 @@ function analyze(::Val{:linf_divb}, du, u, t,
     uEltype = eltype(B1)
     linf_divB = zero(uEltype)
     local_divB = zeros(uEltype, size(B1, 1))
-    for e in eachelement(mesh, dg, cache)
+    @batch reduction=(max, linf_divb) for e in eachelement(mesh, dg, cache)
         compute_local_divergence!(local_divB, e, view.(B, :, e), mesh, dg, cache)
 
         # compute maximum norm