diff --git a/src/callbacks_step/analysis_dg2d.jl b/src/callbacks_step/analysis_dg2d.jl
index e089133fa17..b66a4176317 100644
--- a/src/callbacks_step/analysis_dg2d.jl
+++ b/src/callbacks_step/analysis_dg2d.jl
@@ -285,8 +285,11 @@ function analyze(::Val{:l2_divb}, du, u, t,
                                                          dg, cache, derivative_matrix
         divb = zero(eltype(u))
         for k in eachnode(dg)
-            B1_kj, _, _ = magnetic_field(u[:, k, j, element], equations)
-            _, B2_ik, _ = magnetic_field(u[:, i, k, element], equations)
+            u_kj = get_node_vars(u, equations, dg, k, j, element)
+            u_ik = get_node_vars(u, equations, dg, i, k, element)
+
+            B1_kj, _, _ = magnetic_field(u_kj, equations)
+            _, B2_ik, _ = magnetic_field(u_ik, equations)
 
             divb += (derivative_matrix[i, k] * B1_kj +
                      derivative_matrix[j, k] * B2_ik)
@@ -310,8 +313,11 @@ function analyze(::Val{:l2_divb}, du, u, t,
         Ja21, Ja22 = get_contravariant_vector(2, contravariant_vectors, i, j, element)
         # Compute the transformed divergence
         for k in eachnode(dg)
-            B1_kj, B2_kj, _ = magnetic_field(u[:, k, j, element], equations)
-            B1_ik, B2_ik, _ = magnetic_field(u[:, i, k, element], equations)
+            u_kj = get_node_vars(u, equations, dg, k, j, element)
+            u_ik = get_node_vars(u, equations, dg, i, k, element)
+
+            B1_kj, B2_kj, _ = magnetic_field(u_kj, equations)
+            B1_ik, B2_ik, _ = magnetic_field(u_ik, equations)
 
             divb += (derivative_matrix[i, k] *
                      (Ja11 * B1_kj + Ja12 * B2_kj) +
@@ -334,8 +340,11 @@ function analyze(::Val{:linf_divb}, du, u, t,
         for j in eachnode(dg), i in eachnode(dg)
             divb = zero(eltype(u))
             for k in eachnode(dg)
-                B1_kj, _, _ = magnetic_field(u[:, k, j, element], equations)
-                _, B2_ik, _ = magnetic_field(u[:, i, k, element], equations)
+                u_kj = get_node_vars(u, equations, dg, k, j, element)
+                u_ik = get_node_vars(u, equations, dg, i, k, element)
+
+                B1_kj, _, _ = magnetic_field(u_kj, equations)
+                _, B2_ik, _ = magnetic_field(u_ik, equations)
 
                 divb += (derivative_matrix[i, k] * B1_kj +
                          derivative_matrix[j, k] * B2_ik)
@@ -367,8 +376,11 @@ function analyze(::Val{:linf_divb}, du, u, t,
                                                   element)
             # Compute the transformed divergence
             for k in eachnode(dg)
-                B1_kj, B2_kj, _ = magnetic_field(u[:, k, j, element], equations)
-                B1_ik, B2_ik, _ = magnetic_field(u[:, i, k, element], equations)
+                u_kj = get_node_vars(u, equations, dg, k, j, element)
+                u_ik = get_node_vars(u, equations, dg, i, k, element)
+
+                B1_kj, B2_kj, _ = magnetic_field(u_kj, equations)
+                B1_ik, B2_ik, _ = magnetic_field(u_ik, equations)
 
                 divb += (derivative_matrix[i, k] *
                          (Ja11 * B1_kj + Ja12 * B2_kj) +
diff --git a/src/callbacks_step/analysis_dg3d.jl b/src/callbacks_step/analysis_dg3d.jl
index 5d9e38e5f53..0857d2b4c01 100644
--- a/src/callbacks_step/analysis_dg3d.jl
+++ b/src/callbacks_step/analysis_dg3d.jl
@@ -349,9 +349,13 @@ function analyze(::Val{:l2_divb}, du, u, t,
                                                     element)
         # Compute the transformed divergence
         for l in eachnode(dg)
-            B_ljk = magnetic_field(u[:, l, j, k, element], equations)
-            B_ilk = magnetic_field(u[:, i, l, k, element], equations)
-            B_ijl = magnetic_field(u[:, i, j, l, element], equations)
+            u_ljk = get_node_vars(u, equations, dg, l, j, k, element)
+            u_ilk = get_node_vars(u, equations, dg, i, l, k, element)
+            u_ijl = get_node_vars(u, equations, dg, i, j, l, element)
+
+            B_ljk = magnetic_field(u_ljk, equations)
+            B_ilk = magnetic_field(u_ilk, equations)
+            B_ijl = magnetic_field(u_ijl, equations)
 
             divb += (derivative_matrix[i, l] *
                      (Ja11 * B_ljk[1] + Ja12 * B_ljk[2] + Ja13 * B_ljk[3]) +
@@ -422,9 +426,13 @@ function analyze(::Val{:linf_divb}, du, u, t,
                                                         k, element)
             # Compute the transformed divergence
             for l in eachnode(dg)
-                B_ljk = magnetic_field(u[:, l, j, k, element], equations)
-                B_ilk = magnetic_field(u[:, i, l, k, element], equations)
-                B_ijl = magnetic_field(u[:, i, j, l, element], equations)
+                u_ljk = get_node_vars(u, equations, dg, l, j, k, element)
+                u_ilk = get_node_vars(u, equations, dg, i, l, k, element)
+                u_ijl = get_node_vars(u, equations, dg, i, j, l, element)
+
+                B_ljk = magnetic_field(u_ljk, equations)
+                B_ilk = magnetic_field(u_ilk, equations)
+                B_ijl = magnetic_field(u_ijl, equations)
 
                 divb += (derivative_matrix[i, l] * (Ja11 * B_ljk[1] +
                           Ja12 * B_ljk[2] + Ja13 * B_ljk[3]) +