diff --git a/examples/dgmulti_2d/elixir_navierstokes_convergence.jl b/examples/dgmulti_2d/elixir_navierstokes_convergence.jl
index b529b43a41d..d3647d102fb 100644
--- a/examples/dgmulti_2d/elixir_navierstokes_convergence.jl
+++ b/examples/dgmulti_2d/elixir_navierstokes_convergence.jl
@@ -178,9 +178,11 @@ end
 initial_condition = initial_condition_navier_stokes_convergence_test
 
 # BC types
-velocity_bc_top_bottom = NoSlip((x, t, equations) -> initial_condition_navier_stokes_convergence_test(x,
-                                                                                                      t,
-                                                                                                      equations)[2:3])
+velocity_bc_top_bottom = NoSlip() do x, t, equations
+    u_cons = initial_condition_navier_stokes_convergence_test(x, t, equations)
+    return SVector(u_cons[2] / u_cons[1], u_cons[3] / u_cons[1])
+end
+
 heat_bc_top_bottom = Adiabatic((x, t, equations) -> 0.0)
 boundary_condition_top_bottom = BoundaryConditionNavierStokesWall(velocity_bc_top_bottom,
                                                                   heat_bc_top_bottom)
diff --git a/examples/dgmulti_2d/elixir_navierstokes_convergence_curved.jl b/examples/dgmulti_2d/elixir_navierstokes_convergence_curved.jl
index 92216bf6c29..e969e8c7553 100644
--- a/examples/dgmulti_2d/elixir_navierstokes_convergence_curved.jl
+++ b/examples/dgmulti_2d/elixir_navierstokes_convergence_curved.jl
@@ -186,9 +186,11 @@ end
 initial_condition = initial_condition_navier_stokes_convergence_test
 
 # BC types
-velocity_bc_top_bottom = NoSlip((x, t, equations) -> initial_condition_navier_stokes_convergence_test(x,
-                                                                                                      t,
-                                                                                                      equations)[2:3])
+velocity_bc_top_bottom = NoSlip() do x, t, equations
+    u_cons = initial_condition_navier_stokes_convergence_test(x, t, equations)
+    return SVector(u_cons[2] / u_cons[1], u_cons[3] / u_cons[1])
+end
+
 heat_bc_top_bottom = Adiabatic((x, t, equations) -> 0.0)
 boundary_condition_top_bottom = BoundaryConditionNavierStokesWall(velocity_bc_top_bottom,
                                                                   heat_bc_top_bottom)
diff --git a/examples/dgmulti_3d/elixir_navierstokes_convergence.jl b/examples/dgmulti_3d/elixir_navierstokes_convergence.jl
index 5fa0ad7ce60..d7deb0a8bd4 100644
--- a/examples/dgmulti_3d/elixir_navierstokes_convergence.jl
+++ b/examples/dgmulti_3d/elixir_navierstokes_convergence.jl
@@ -221,9 +221,11 @@ end
 initial_condition = initial_condition_navier_stokes_convergence_test
 
 # BC types
-velocity_bc_top_bottom = NoSlip((x, t, equations) -> initial_condition_navier_stokes_convergence_test(x,
-                                                                                                      t,
-                                                                                                      equations)[2:4])
+velocity_bc_top_bottom = NoSlip() do x, t, equations
+    u_cons = initial_condition_navier_stokes_convergence_test(x, t, equations)
+    return SVector(u_cons[2] / u_cons[1], u_cons[3] / u_cons[1], u_cons[4] / u_cons[1])
+end
+
 heat_bc_top_bottom = Adiabatic((x, t, equations) -> 0.0)
 boundary_condition_top_bottom = BoundaryConditionNavierStokesWall(velocity_bc_top_bottom,
                                                                   heat_bc_top_bottom)
diff --git a/examples/dgmulti_3d/elixir_navierstokes_convergence_curved.jl b/examples/dgmulti_3d/elixir_navierstokes_convergence_curved.jl
index c58d78d2581..0860fbc4285 100644
--- a/examples/dgmulti_3d/elixir_navierstokes_convergence_curved.jl
+++ b/examples/dgmulti_3d/elixir_navierstokes_convergence_curved.jl
@@ -229,9 +229,11 @@ end
 initial_condition = initial_condition_navier_stokes_convergence_test
 
 # BC types
-velocity_bc_top_bottom = NoSlip((x, t, equations) -> initial_condition_navier_stokes_convergence_test(x,
-                                                                                                      t,
-                                                                                                      equations)[2:4])
+velocity_bc_top_bottom = NoSlip() do x, t, equations
+    u_cons = initial_condition_navier_stokes_convergence_test(x, t, equations)
+    return SVector(u_cons[2] / u_cons[1], u_cons[3] / u_cons[1], u_cons[4] / u_cons[1])
+end
+
 heat_bc_top_bottom = Adiabatic((x, t, equations) -> 0.0)
 boundary_condition_top_bottom = BoundaryConditionNavierStokesWall(velocity_bc_top_bottom,
                                                                   heat_bc_top_bottom)
diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl
index 43ce7ac1ed6..27b58b3a3f4 100644
--- a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl
+++ b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl
@@ -83,12 +83,14 @@ heat_airfoil = Adiabatic((x, t, equations) -> 0.0)
 boundary_conditions_airfoil = BoundaryConditionNavierStokesWall(velocity_airfoil,
                                                                 heat_airfoil)
 
-function momenta_initial_condition_mach08_flow(x, t, equations)
-    u = initial_condition_mach08_flow(x, t, equations)
-    momenta = SVector(u[2], u[3])
+function velocities_initial_condition_mach08_flow(x, t, equations)
+    u_cons = initial_condition_mach08_flow(x, t, equations)
+    return SVector(u_cons[2] / u_cons[1], u_cons[3] / u_cons[1])
 end
-velocity_bc_square = NoSlip((x, t, equations) -> momenta_initial_condition_mach08_flow(x, t,
-                                                                                       equations))
+
+velocity_bc_square = NoSlip((x, t, equations) -> velocities_initial_condition_mach08_flow(x,
+                                                                                          t,
+                                                                                          equations))
 
 heat_bc_square = Adiabatic((x, t, equations) -> 0.0)
 boundary_condition_square = BoundaryConditionNavierStokesWall(velocity_bc_square,
diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_convergence.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_convergence.jl
index 54ec09d2be8..d3f7a01a04c 100644
--- a/examples/p4est_2d_dgsem/elixir_navierstokes_convergence.jl
+++ b/examples/p4est_2d_dgsem/elixir_navierstokes_convergence.jl
@@ -179,9 +179,10 @@ initial_condition = initial_condition_navier_stokes_convergence_test
 
 # BC types
 velocity_bc_top_bottom = NoSlip() do x, t, equations
-    u = initial_condition_navier_stokes_convergence_test(x, t, equations)
-    return SVector(u[2], u[3])
+    u_cons = initial_condition_navier_stokes_convergence_test(x, t, equations)
+    return SVector(u_cons[2] / u_cons[1], u_cons[3] / u_cons[1])
 end
+
 heat_bc_top_bottom = Adiabatic((x, t, equations) -> 0.0)
 boundary_condition_top_bottom = BoundaryConditionNavierStokesWall(velocity_bc_top_bottom,
                                                                   heat_bc_top_bottom)
diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_convergence_nonperiodic.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_convergence_nonperiodic.jl
index b4177fe8538..9125d2e03bf 100644
--- a/examples/p4est_2d_dgsem/elixir_navierstokes_convergence_nonperiodic.jl
+++ b/examples/p4est_2d_dgsem/elixir_navierstokes_convergence_nonperiodic.jl
@@ -178,9 +178,11 @@ end
 initial_condition = initial_condition_navier_stokes_convergence_test
 
 # BC types
-velocity_bc_top_bottom = NoSlip((x, t, equations) -> initial_condition_navier_stokes_convergence_test(x,
-                                                                                                      t,
-                                                                                                      equations)[2:3])
+velocity_bc_top_bottom = NoSlip() do x, t, equations
+    u_cons = initial_condition_navier_stokes_convergence_test(x, t, equations)
+    return SVector(u_cons[2] / u_cons[1], u_cons[3] / u_cons[1])
+end
+
 heat_bc_top_bottom = Adiabatic((x, t, equations) -> 0.0)
 boundary_condition_top_bottom = BoundaryConditionNavierStokesWall(velocity_bc_top_bottom,
                                                                   heat_bc_top_bottom)
diff --git a/examples/p4est_3d_dgsem/elixir_navierstokes_convergence.jl b/examples/p4est_3d_dgsem/elixir_navierstokes_convergence.jl
index c640b255b05..06f19432ee9 100644
--- a/examples/p4est_3d_dgsem/elixir_navierstokes_convergence.jl
+++ b/examples/p4est_3d_dgsem/elixir_navierstokes_convergence.jl
@@ -223,9 +223,10 @@ initial_condition = initial_condition_navier_stokes_convergence_test
 
 # BC types
 velocity_bc_top_bottom = NoSlip() do x, t, equations
-    u = initial_condition_navier_stokes_convergence_test(x, t, equations)
-    return SVector(u[2], u[3], u[4])
+    u_cons = initial_condition_navier_stokes_convergence_test(x, t, equations)
+    return SVector(u_cons[2] / u_cons[1], u_cons[3] / u_cons[1], u_cons[4] / u_cons[1])
 end
+
 heat_bc_top_bottom = Adiabatic((x, t, equations) -> 0.0)
 boundary_condition_top_bottom = BoundaryConditionNavierStokesWall(velocity_bc_top_bottom,
                                                                   heat_bc_top_bottom)
@@ -264,3 +265,5 @@ time_int_tol = 1e-8
 sol = solve(ode, RDPK3SpFSAL49(); abstol = time_int_tol, reltol = time_int_tol, dt = 1e-5,
             ode_default_options()..., callback = callbacks)
 summary_callback() # print the timer summary
+
+display(analysis_callback(sol))
diff --git a/examples/tree_1d_dgsem/elixir_navierstokes_convergence_walls.jl b/examples/tree_1d_dgsem/elixir_navierstokes_convergence_walls.jl
index d534b7b189b..83364c5f356 100644
--- a/examples/tree_1d_dgsem/elixir_navierstokes_convergence_walls.jl
+++ b/examples/tree_1d_dgsem/elixir_navierstokes_convergence_walls.jl
@@ -121,9 +121,10 @@ end
 initial_condition = initial_condition_navier_stokes_convergence_test
 
 # BC types
-velocity_bc_left_right = NoSlip((x, t, equations) -> initial_condition_navier_stokes_convergence_test(x,
-                                                                                                      t,
-                                                                                                      equations)[2])
+velocity_bc_left_right = NoSlip() do x, t, equations
+    u_cons = initial_condition_navier_stokes_convergence_test(x, t, equations)
+    return u_cons[2] / u_cons[1]
+end
 
 heat_bc_left = Isothermal((x, t, equations) -> Trixi.temperature(initial_condition_navier_stokes_convergence_test(x,
                                                                                                                   t,
diff --git a/examples/tree_1d_dgsem/elixir_navierstokes_convergence_walls_amr.jl b/examples/tree_1d_dgsem/elixir_navierstokes_convergence_walls_amr.jl
index 4f9cd6365fa..b114e84aa5f 100644
--- a/examples/tree_1d_dgsem/elixir_navierstokes_convergence_walls_amr.jl
+++ b/examples/tree_1d_dgsem/elixir_navierstokes_convergence_walls_amr.jl
@@ -121,9 +121,10 @@ end
 initial_condition = initial_condition_navier_stokes_convergence_test
 
 # BC types
-velocity_bc_left_right = NoSlip((x, t, equations) -> initial_condition_navier_stokes_convergence_test(x,
-                                                                                                      t,
-                                                                                                      equations)[2])
+velocity_bc_left_right = NoSlip() do x, t, equations
+    u_cons = initial_condition_navier_stokes_convergence_test(x, t, equations)
+    return u_cons[2] / u_cons[1]
+end
 
 heat_bc_left = Isothermal((x, t, equations) -> Trixi.temperature(initial_condition_navier_stokes_convergence_test(x,
                                                                                                                   t,
diff --git a/examples/tree_2d_dgsem/elixir_navierstokes_convergence.jl b/examples/tree_2d_dgsem/elixir_navierstokes_convergence.jl
index b0c8678baad..7b4a3f00e19 100644
--- a/examples/tree_2d_dgsem/elixir_navierstokes_convergence.jl
+++ b/examples/tree_2d_dgsem/elixir_navierstokes_convergence.jl
@@ -179,8 +179,8 @@ initial_condition = initial_condition_navier_stokes_convergence_test
 
 # BC types
 velocity_bc_top_bottom = NoSlip() do x, t, equations
-    u = initial_condition_navier_stokes_convergence_test(x, t, equations)
-    return SVector(u[2], u[3])
+    u_cons = initial_condition_navier_stokes_convergence_test(x, t, equations)
+    return SVector(u_cons[2] / u_cons[1], u_cons[3] / u_cons[1])
 end
 heat_bc_top_bottom = Adiabatic((x, t, equations) -> 0.0)
 boundary_condition_top_bottom = BoundaryConditionNavierStokesWall(velocity_bc_top_bottom,
diff --git a/examples/tree_3d_dgsem/elixir_navierstokes_convergence.jl b/examples/tree_3d_dgsem/elixir_navierstokes_convergence.jl
index 3ada9503c6a..cf16ebff332 100644
--- a/examples/tree_3d_dgsem/elixir_navierstokes_convergence.jl
+++ b/examples/tree_3d_dgsem/elixir_navierstokes_convergence.jl
@@ -223,8 +223,8 @@ initial_condition = initial_condition_navier_stokes_convergence_test
 
 # BC types
 velocity_bc_top_bottom = NoSlip() do x, t, equations
-    u = initial_condition_navier_stokes_convergence_test(x, t, equations)
-    return SVector(u[2], u[3], u[4])
+    u_cons = initial_condition_navier_stokes_convergence_test(x, t, equations)
+    return SVector(u_cons[2] / u_cons[1], u_cons[3] / u_cons[1], u_cons[4] / u_cons[1])
 end
 heat_bc_top_bottom = Adiabatic((x, t, equations) -> 0.0)
 boundary_condition_top_bottom = BoundaryConditionNavierStokesWall(velocity_bc_top_bottom,
diff --git a/test/test_parabolic_1d.jl b/test/test_parabolic_1d.jl
index 48b2b9ce4d6..e13c1ed31e3 100644
--- a/test/test_parabolic_1d.jl
+++ b/test/test_parabolic_1d.jl
@@ -130,14 +130,14 @@ end
     @test_trixi_include(joinpath(examples_dir(), "tree_1d_dgsem",
                                  "elixir_navierstokes_convergence_walls.jl"),
                         l2=[
-                            0.00047023310868269237,
-                            0.00032181736027057234,
-                            0.0014966266486095025
+                            0.0004702331100298379,
+                            0.0003218173539588441,
+                            0.001496626616191212
                         ],
                         linf=[
-                            0.002996375101363302,
-                            0.0028639041695096433,
-                            0.012691132694550689
+                            0.0029963751636357117,
+                            0.002863904080244726,
+                            0.012691132434907715
                         ])
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
@@ -157,14 +157,14 @@ end
                                                                                 Prandtl = prandtl_number(),
                                                                                 gradient_variables = GradientVariablesEntropy()),
                         l2=[
-                            0.0004608500483647771,
-                            0.00032431091222851285,
-                            0.0015159733360626845
+                            0.00046085004909354776,
+                            0.0003243109084492897,
+                            0.0015159733164383632
                         ],
                         linf=[
-                            0.002754803146635787,
-                            0.0028567713744625124,
-                            0.012941793784197131
+                            0.0027548031865172184,
+                            0.0028567713569609024,
+                            0.012941793735691931
                         ])
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
@@ -183,14 +183,14 @@ end
                                                                                 mu = mu(),
                                                                                 Prandtl = prandtl_number()),
                         l2=[
-                            2.5278845598681636e-5,
-                            2.5540145802666872e-5,
-                            0.0001211867535580826
+                            2.5279985418616284e-5,
+                            2.5543804380243317e-5,
+                            0.000121197807780003
                         ],
                         linf=[
-                            0.0001466387202588848,
-                            0.00019422419092429135,
-                            0.0009556449835592673
+                            0.00014663941942760772,
+                            0.00019421961959378687,
+                            0.0009556612427097377
                         ])
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
@@ -210,14 +210,14 @@ end
                                                                                 Prandtl = prandtl_number(),
                                                                                 gradient_variables = GradientVariablesEntropy()),
                         l2=[
-                            2.4593521887223632e-5,
-                            2.3928212900127102e-5,
-                            0.00011252332663824173
+                            2.459436737779634e-5,
+                            2.393023570132773e-5,
+                            0.00011253258031437205
                         ],
                         linf=[
-                            0.00011850494672183132,
-                            0.00018987676556476442,
-                            0.0009597461727750556
+                            0.00011850784723010932,
+                            0.00018988145760703384,
+                            0.0009597937547916047
                         ])
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl
index ceefb65e99b..dc3e776154f 100644
--- a/test/test_parabolic_2d.jl
+++ b/test/test_parabolic_2d.jl
@@ -124,16 +124,16 @@ end
                                  "elixir_navierstokes_convergence.jl"),
                         cells_per_dimension=(4, 4), tspan=(0.0, 0.1),
                         l2=[
-                            0.0015355076812510957,
-                            0.0033843168272696756,
-                            0.0036531858107443434,
-                            0.009948436427519214
+                            0.0015355076812512347,
+                            0.0033843168272690953,
+                            0.0036531858107444093,
+                            0.009948436427520873
                         ],
                         linf=[
-                            0.005522560467190019,
-                            0.013425258500730508,
-                            0.013962115643482154,
-                            0.027483102120502423
+                            0.005522560467186466,
+                            0.013425258500736614,
+                            0.013962115643462503,
+                            0.0274831021205042
                         ])
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
@@ -150,16 +150,16 @@ end
                                  "elixir_navierstokes_convergence_curved.jl"),
                         cells_per_dimension=(4, 4), tspan=(0.0, 0.1),
                         l2=[
-                            0.004255101916146187,
-                            0.011118488923215765,
-                            0.011281831283462686,
-                            0.03573656447388509
+                            0.004255101916146402,
+                            0.011118488923215394,
+                            0.011281831283462784,
+                            0.035736564473886
                         ],
                         linf=[
-                            0.015071710669706473,
-                            0.04103132025858458,
-                            0.03990424085750277,
-                            0.1309401718598764
+                            0.015071710669707805,
+                            0.04103132025860057,
+                            0.03990424085748012,
+                            0.13094017185987106
                         ],)
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
@@ -273,16 +273,16 @@ end
                                                                                        energy_internal,
                                                                                        enstrophy)),
                         l2=[
-                            0.002111672530658797,
-                            0.0034322351490857846,
-                            0.0038742528195910416,
-                            0.012469246082568561
+                            0.0021116725306624543,
+                            0.003432235149083229,
+                            0.003874252819605527,
+                            0.012469246082535005
                         ],
                         linf=[
-                            0.012006418939223495,
-                            0.035520871209746126,
-                            0.024512747492231427,
-                            0.11191122588756564
+                            0.012006418939279007,
+                            0.03552087120962882,
+                            0.02451274749189282,
+                            0.11191122588626357
                         ])
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
@@ -303,16 +303,16 @@ end
                                                                                                                                               equations),
                                                                                              equations)),
                         l2=[
-                            0.002103629650383915,
-                            0.003435843933396454,
-                            0.00386735987813341,
-                            0.012670355349235728
+                            0.0021036296503840883,
+                            0.003435843933397192,
+                            0.003867359878114748,
+                            0.012670355349293195
                         ],
                         linf=[
-                            0.012006261793147788,
-                            0.03550212518982032,
-                            0.025107947319661185,
-                            0.11647078036571124
+                            0.01200626179308184,
+                            0.03550212518997239,
+                            0.025107947320178275,
+                            0.11647078036751068
                         ])
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
@@ -330,16 +330,16 @@ end
                         initial_refinement_level=2, tspan=(0.0, 0.1),
                         gradient_variables=GradientVariablesEntropy(),
                         l2=[
-                            0.0021403742517389513,
-                            0.0034258287094908572,
-                            0.0038915122886898517,
-                            0.012506862343013842
+                            0.002140374251729127,
+                            0.003425828709496601,
+                            0.0038915122887358097,
+                            0.012506862342858291
                         ],
                         linf=[
-                            0.012244412004628336,
-                            0.03507559186162224,
-                            0.024580892345558894,
-                            0.11425600758350107
+                            0.012244412004772665,
+                            0.03507559186131113,
+                            0.02458089234472249,
+                            0.11425600758024679
                         ])
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
@@ -361,16 +361,16 @@ end
                                                                                                                                               equations),
                                                                                              equations)),
                         l2=[
-                            0.0021349737347844907,
-                            0.0034301388278203033,
-                            0.0038928324474291572,
-                            0.012693611436230873
+                            0.0021349737347923716,
+                            0.0034301388278178365,
+                            0.0038928324473968836,
+                            0.012693611436338
                         ],
                         linf=[
-                            0.01224423627586213,
-                            0.035054066314102905,
-                            0.025099598504931965,
-                            0.11795616324751634
+                            0.012244236275761766,
+                            0.03505406631430898,
+                            0.025099598505644406,
+                            0.11795616324985403
                         ])
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
@@ -388,16 +388,16 @@ end
                         initial_refinement_level=2, tspan=(0.0, 0.1),
                         volume_integral=VolumeIntegralFluxDifferencing(flux_central),
                         l2=[
-                            0.0021116725306633594,
-                            0.0034322351490827557,
-                            0.0038742528196093542,
-                            0.012469246082526909
+                            0.0021116725306612075,
+                            0.0034322351490838703,
+                            0.0038742528196011594,
+                            0.012469246082545557
                         ],
                         linf=[
-                            0.012006418939291663,
-                            0.035520871209594115,
-                            0.024512747491801577,
-                            0.11191122588591007
+                            0.012006418939262131,
+                            0.0355208712096602,
+                            0.024512747491999436,
+                            0.11191122588669522
                         ])
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
@@ -431,11 +431,15 @@ end
                 ode_default_options()..., callback = callbacks)
     l2_error, linf_error = analysis_callback(sol)
     @test l2_error ≈
-          [0.00024296959173852447; 0.0002093263158670915; 0.0005390572390977262;
-           0.00026753561392341537]
+          [0.00024296959174050973;
+           0.00020932631586399853;
+           0.0005390572390981241;
+           0.00026753561391316933]
     @test linf_error ≈
-          [0.0016210102053424436; 0.002593287648655501; 0.002953907343823712;
-           0.002077119120180271]
+          [0.0016210102053486608;
+           0.0025932876486537016;
+           0.0029539073438284817;
+           0.0020771191202548778]
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
     let
@@ -613,16 +617,16 @@ end
                                  "elixir_navierstokes_convergence.jl"),
                         initial_refinement_level=1, tspan=(0.0, 0.2),
                         l2=[
-                            0.0003811978985836709,
-                            0.0005874314969169538,
-                            0.0009142898787923481,
-                            0.0011613918899727263
+                            0.0003811978986531135,
+                            0.0005874314969137914,
+                            0.0009142898787681551,
+                            0.0011613918893790497
                         ],
                         linf=[
-                            0.0021633623982135752,
-                            0.009484348274135372,
-                            0.004231572066492217,
-                            0.011661660275365193
+                            0.0021633623985426453,
+                            0.009484348273965089,
+                            0.0042315720663082534,
+                            0.011661660264076446
                         ])
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
@@ -639,16 +643,16 @@ end
                                  "elixir_navierstokes_convergence_nonperiodic.jl"),
                         initial_refinement_level=1, tspan=(0.0, 0.2),
                         l2=[
-                            0.00040364962558511795,
-                            0.0005869762481506936,
-                            0.00091488537427274,
-                            0.0011984191566376762
+                            0.0004036496258545996,
+                            0.0005869762480189079,
+                            0.0009148853742181908,
+                            0.0011984191532764543
                         ],
                         linf=[
-                            0.0024993634941723464,
-                            0.009487866203944725,
-                            0.004505829506628117,
-                            0.011634902776245681
+                            0.0024993634989209923,
+                            0.009487866203496731,
+                            0.004505829506103787,
+                            0.011634902753554499
                         ])
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
diff --git a/test/test_parabolic_3d.jl b/test/test_parabolic_3d.jl
index 2690a08cbb9..a5e63c7486c 100644
--- a/test/test_parabolic_3d.jl
+++ b/test/test_parabolic_3d.jl
@@ -17,18 +17,18 @@ isdir(outdir) && rm(outdir, recursive = true)
                                  "elixir_navierstokes_convergence.jl"),
                         cells_per_dimension=(4, 4, 4), tspan=(0.0, 0.1),
                         l2=[
-                            0.0005532847115849239,
-                            0.000659263490965341,
-                            0.0007776436127362806,
-                            0.0006592634909662951,
-                            0.0038073628897809185
+                            0.000553284711585015,
+                            0.0006592634909666629,
+                            0.0007776436127373607,
+                            0.0006592634909664286,
+                            0.0038073628897819524
                         ],
                         linf=[
-                            0.0017039861523615585,
-                            0.002628561703560073,
-                            0.003531057425112172,
-                            0.0026285617036090336,
-                            0.015587829540351095
+                            0.0017039861523628907,
+                            0.002628561703550747,
+                            0.0035310574250866367,
+                            0.002628561703585053,
+                            0.015587829540340437
                         ])
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
@@ -45,18 +45,18 @@ end
                                  "elixir_navierstokes_convergence_curved.jl"),
                         cells_per_dimension=(4, 4, 4), tspan=(0.0, 0.1),
                         l2=[
-                            0.0014027227251207474,
-                            0.0021322235533273513,
-                            0.0027873741447455194,
-                            0.0024587473070627423,
-                            0.00997836818019202
+                            0.001402722725120774,
+                            0.0021322235533272546,
+                            0.002787374144745514,
+                            0.002458747307062751,
+                            0.009978368180191861
                         ],
                         linf=[
-                            0.006341750402837576,
-                            0.010306014252246865,
-                            0.01520740250924979,
-                            0.010968264045485565,
-                            0.047454389831591115
+                            0.0063417504028402405,
+                            0.010306014252245865,
+                            0.015207402509253565,
+                            0.010968264045485343,
+                            0.04745438983155026
                         ])
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
@@ -101,18 +101,18 @@ end
                                  "elixir_navierstokes_convergence.jl"),
                         initial_refinement_level=2, tspan=(0.0, 0.1),
                         l2=[
-                            0.0019582188528512257,
-                            0.002653449504302844,
-                            0.002898264205184629,
-                            0.002653449504302853,
-                            0.009511572365085706
+                            0.0019582188528520267,
+                            0.002653449504302849,
+                            0.002898264205184317,
+                            0.0026534495043028534,
+                            0.009511572365092744
                         ],
                         linf=[
-                            0.013680656759085918,
-                            0.0356910450154318,
-                            0.023526343547736236,
-                            0.035691045015431855,
-                            0.11482570604041165
+                            0.013680656759089693,
+                            0.03569104501543785,
+                            0.023526343547761893,
+                            0.03569104501543733,
+                            0.11482570604049513
                         ])
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
@@ -133,18 +133,18 @@ end
                                                                                                                                               equations),
                                                                                              equations)),
                         l2=[
-                            0.00195468651965362,
+                            0.001954686519653731,
                             0.0026554367897028506,
-                            0.002892730402724066,
-                            0.002655436789702817,
-                            0.009596351796609566
+                            0.0028927304027240026,
+                            0.0026554367897028437,
+                            0.00959635179660988
                         ],
                         linf=[
-                            0.013680508110645473,
-                            0.035673446359424356,
-                            0.024024936779729028,
-                            0.03567344635942474,
-                            0.11839497110809383
+                            0.013680508110646583,
+                            0.03567344635942522,
+                            0.024024936779738822,
+                            0.035673446359425674,
+                            0.11839497110814179
                         ])
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
@@ -162,18 +162,18 @@ end
                         initial_refinement_level=2, tspan=(0.0, 0.1),
                         gradient_variables=GradientVariablesEntropy(),
                         l2=[
-                            0.0019770444875099307,
-                            0.0026524750946399327,
-                            0.00290860030832445,
-                            0.0026524750946399396,
-                            0.009509568981439294
+                            0.0019770444875097737,
+                            0.002652475094640119,
+                            0.0029086003083239236,
+                            0.002652475094640097,
+                            0.009509568981441823
                         ],
                         linf=[
-                            0.01387936112914212,
-                            0.03526260609304053,
-                            0.023554197097368997,
-                            0.035262606093040896,
-                            0.11719963716509518
+                            0.013879361129145007,
+                            0.035262606093049195,
+                            0.02355419709739138,
+                            0.03526260609304984,
+                            0.11719963716518933
                         ])
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
@@ -195,18 +195,18 @@ end
                                                                                                                                               equations),
                                                                                              equations)),
                         l2=[
-                            0.001974631423398113,
-                            0.002654768259143932,
-                            0.002907031063651286,
-                            0.002654768259143901,
-                            0.009587792882971452
+                            0.0019746314233993435,
+                            0.0026547682591448896,
+                            0.0029070310636460494,
+                            0.0026547682591448922,
+                            0.00958779288300152
                         ],
                         linf=[
-                            0.01387919380137137,
-                            0.035244084526358944,
-                            0.02398614622061363,
-                            0.03524408452635828,
-                            0.12005056512506407
+                            0.013879193801400458,
+                            0.03524408452641245,
+                            0.023986146220843566,
+                            0.035244084526412915,
+                            0.1200505651257302
                         ])
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
@@ -224,18 +224,18 @@ end
                         initial_refinement_level=2, tspan=(0.0, 0.1),
                         volume_integral=VolumeIntegralFluxDifferencing(flux_central),
                         l2=[
-                            0.0019582188528180213,
-                            0.002653449504301736,
-                            0.0028982642051960006,
-                            0.0026534495043017384,
-                            0.009511572364811033
+                            0.0019582188528208754,
+                            0.0026534495043017935,
+                            0.002898264205195059,
+                            0.0026534495043017917,
+                            0.009511572364832972
                         ],
                         linf=[
-                            0.013680656758949583,
-                            0.035691045015224444,
-                            0.02352634354676752,
-                            0.035691045015223424,
-                            0.11482570603751441
+                            0.013680656758958687,
+                            0.03569104501523916,
+                            0.02352634354684648,
+                            0.03569104501523987,
+                            0.11482570603774533
                         ])
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
@@ -268,20 +268,16 @@ end
                 dt = 1e-5,
                 ode_default_options()..., callback = callbacks)
     l2_error, linf_error = analysis_callback(sol)
-    @test l2_error ≈ [
-        0.0003109336253407314,
-        0.0006473493036803503,
-        0.0007705277238213672,
-        0.0006280517917198335,
-        0.000903927789884075
-    ]
-    @test linf_error ≈ [
-        0.0023694155365339142,
-        0.010634932622402863,
-        0.006772070862236412,
-        0.010640551561726901,
-        0.019256819038719897
-    ]
+    @test l2_error ≈ [0.00031093362536287433;
+           0.0006473493036800964;
+           0.0007705277238221976;
+           0.0006280517917194624;
+           0.0009039277899421355]
+    @test linf_error ≈ [0.0023694155363713776;
+           0.01063493262248095;
+           0.006772070862041679;
+           0.010640551561807883;
+           0.019256819037817507]
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
     let
@@ -373,18 +369,18 @@ end
                                  "elixir_navierstokes_convergence.jl"),
                         initial_refinement_level=2, tspan=(0.0, 0.1),
                         l2=[
-                            0.00026599105554982194,
-                            0.000461877794472316,
-                            0.0005424899076052261,
-                            0.0004618777944723191,
-                            0.0015846392581126832
+                            0.00026599105557723507,
+                            0.00046187779448444603,
+                            0.0005424899076194272,
+                            0.00046187779448445546,
+                            0.0015846392584275121
                         ],
                         linf=[
-                            0.0025241668929956163,
-                            0.006308461681816373,
-                            0.004334939663169113,
-                            0.006308461681804009,
-                            0.03176343480493493
+                            0.0025241668964857134,
+                            0.006308461684409397,
+                            0.004334939668473314,
+                            0.006308461684396753,
+                            0.03176343483364796
                         ])
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)