diff --git a/NEWS.md b/NEWS.md
index bca7d03e95..9371cafa07 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -4,6 +4,14 @@ Trixi.jl follows the interpretation of [semantic versioning (semver)](https://ju
 used in the Julia ecosystem. Notable changes will be documented in this file
 for human readability.
 
+## Changes in the v0.8 lifecycle
+
+#### Changed
+
+- The AMR routines for `P4estMesh` and `T8codeMesh` were changed to work on the product
+  of the Jacobian and the conserved variables instead of the conserved variables only
+  to make AMR fully conservative ([#2028]). This may change AMR results slightly.
+
 ## Changes when updating to v0.8 from v0.7.x
 
 #### Added
diff --git a/Project.toml b/Project.toml
index d5627651d6..923cf5248a 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "Trixi"
 uuid = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb"
 authors = ["Michael Schlottke-Lakemper <m.schlottke-lakemper@hlrs.de>", "Gregor Gassner <ggassner@uni-koeln.de>", "Hendrik Ranocha <mail@ranocha.de>", "Andrew R. Winters <andrew.ross.winters@liu.se>", "Jesse Chan <jesse.chan@rice.edu>"]
-version = "0.8.7-DEV"
+version = "0.8.8-DEV"
 
 [deps]
 CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2"
diff --git a/examples/p4est_2d_dgsem/elixir_euler_weak_blast_wave_amr.jl b/examples/p4est_2d_dgsem/elixir_euler_weak_blast_wave_amr.jl
new file mode 100644
index 0000000000..6868067371
--- /dev/null
+++ b/examples/p4est_2d_dgsem/elixir_euler_weak_blast_wave_amr.jl
@@ -0,0 +1,117 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+
+equations = CompressibleEulerEquations2D(1.4)
+
+function initial_condition_weak_blast_wave(x, t, equations::CompressibleEulerEquations2D)
+    # Set up polar coordinates
+    inicenter = SVector(0.0, 0.0)
+    x_norm = x[1] - inicenter[1]
+    y_norm = x[2] - inicenter[2]
+    r = sqrt(x_norm^2 + y_norm^2)
+
+    r0 = 0.2
+    E = 1
+    p0_inner = 3
+    p0_outer = 1
+
+    # Calculate primitive variables
+    rho = 1.1
+    v1 = 0.0
+    v2 = 0.0
+    p = r > r0 ? p0_outer : p0_inner
+
+    return prim2cons(SVector(rho, v1, v2, p), equations)
+end
+
+initial_condition = initial_condition_weak_blast_wave
+
+# Get the DG approximation space
+
+# Activate the shock capturing + flux differencing
+surface_flux = flux_lax_friedrichs
+volume_flux = flux_ranocha
+polydeg = 4
+basis = LobattoLegendreBasis(polydeg)
+indicator_sc = IndicatorHennemannGassner(equations, basis,
+                                         alpha_max = 0.5,
+                                         alpha_min = 0.001,
+                                         alpha_smooth = true,
+                                         variable = density_pressure)
+volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
+                                                 volume_flux_dg = volume_flux,
+                                                 volume_flux_fv = surface_flux)
+
+solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux,
+               volume_integral = volume_integral)
+
+###############################################################################
+
+# Affine type mapping to take the [-1,1]^2 domain
+# and warp it as described in https://arxiv.org/abs/2012.12040
+# Warping with the coefficient 0.2 is even more extreme.
+function mapping_twist(xi, eta)
+    y = eta + 0.125 * cos(1.5 * pi * xi) * cos(0.5 * pi * eta)
+    x = xi + 0.125 * cos(0.5 * pi * xi) * cos(2.0 * pi * y)
+    return SVector(x, y)
+end
+
+# The mesh below can be made periodic
+# Create P4estMesh with 8 x 8 trees
+trees_per_dimension = (8, 8)
+mesh = P4estMesh(trees_per_dimension, polydeg = 4,
+                 mapping = mapping_twist,
+                 initial_refinement_level = 0,
+                 periodicity = true)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 1.2)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 400
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
+                                     save_analysis = true,
+                                     extra_analysis_errors = (:conservation_error,))
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(dt = 0.2,
+                                     save_initial_solution = true,
+                                     save_final_solution = true)
+
+amr_indicator = IndicatorLöhner(semi, variable = Trixi.density)
+amr_controller = ControllerThreeLevel(semi, amr_indicator,
+                                      base_level = 0,
+                                      med_level = 1, med_threshold = 0.05,
+                                      max_level = 2, max_threshold = 0.1)
+amr_callback = AMRCallback(semi, amr_controller,
+                           interval = 5,
+                           adapt_initial_condition = true,
+                           adapt_initial_condition_only_refine = true)
+
+stepsize_callback = StepsizeCallback(cfl = 0.5)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback,
+                        alive_callback,
+                        save_solution,
+                        amr_callback,
+                        stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
+            dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep = false, callback = callbacks);#, maxiters=4);
+summary_callback() # print the timer summary
diff --git a/examples/p4est_3d_dgsem/elixir_euler_weak_blast_wave_amr.jl b/examples/p4est_3d_dgsem/elixir_euler_weak_blast_wave_amr.jl
new file mode 100644
index 0000000000..b5b5622000
--- /dev/null
+++ b/examples/p4est_3d_dgsem/elixir_euler_weak_blast_wave_amr.jl
@@ -0,0 +1,114 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+
+equations = CompressibleEulerEquations3D(1.4)
+
+function initial_condition_weak_blast_wave(x, t,
+                                           equations::CompressibleEulerEquations3D)
+    # Set up polar coordinates
+    inicenter = SVector(0.0, 0.0, 0.0)
+    x_norm = x[1] - inicenter[1]
+    y_norm = x[2] - inicenter[2]
+    z_norm = x[3] - inicenter[3]
+    r = sqrt(x_norm^2 + y_norm^2 + z_norm^2)
+
+    r0 = 0.2
+    E = 1.0
+    p0_inner = 3
+    p0_outer = 1
+
+    # Calculate primitive variables
+    rho = 1.1
+    v1 = 0.0
+    v2 = 0.0
+    v3 = 0.0
+    p = r > r0 ? p0_outer : p0_inner
+
+    return prim2cons(SVector(rho, v1, v2, v3, p), equations)
+end
+
+initial_condition = initial_condition_weak_blast_wave
+
+surface_flux = flux_lax_friedrichs
+volume_flux = flux_ranocha
+polydeg = 4
+basis = LobattoLegendreBasis(polydeg)
+indicator_sc = IndicatorHennemannGassner(equations, basis,
+                                         alpha_max = 1.0,
+                                         alpha_min = 0.001,
+                                         alpha_smooth = true,
+                                         variable = density_pressure)
+volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
+                                                 volume_flux_dg = volume_flux,
+                                                 volume_flux_fv = surface_flux)
+
+solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux,
+               volume_integral = volume_integral)
+
+# Setup a periodic mesh with 4 x 4 x 4 trees and 8 x 8 x 8 elements
+trees_per_dimension = (4, 4, 4)
+
+# Affine type mapping to take the [-1,1]^3 domain
+# and warp it as described in https://arxiv.org/abs/2012.12040
+function mapping_twist(xi, eta, zeta)
+    y = eta + 1 / 6 * (cos(1.5 * pi * xi) * cos(0.5 * pi * eta) * cos(0.5 * pi * zeta))
+
+    x = xi + 1 / 6 * (cos(0.5 * pi * xi) * cos(2 * pi * y) * cos(0.5 * pi * zeta))
+
+    z = zeta + 1 / 6 * (cos(0.5 * pi * x) * cos(pi * y) * cos(0.5 * pi * zeta))
+
+    return SVector(x, y, z)
+end
+
+mesh = P4estMesh(trees_per_dimension,
+                 polydeg = 2,
+                 mapping = mapping_twist,
+                 initial_refinement_level = 1,
+                 periodicity = true)
+
+# Create the semidiscretization object
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 1.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
+                                     extra_analysis_errors = (:conservation_error,))
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+amr_indicator = IndicatorLöhner(semi, variable = Trixi.density)
+amr_controller = ControllerThreeLevel(semi, amr_indicator,
+                                      base_level = 1,
+                                      med_level = 2, med_threshold = 0.05,
+                                      max_level = 3, max_threshold = 0.15)
+amr_callback = AMRCallback(semi, amr_controller,
+                           interval = 1,
+                           adapt_initial_condition = false,
+                           adapt_initial_condition_only_refine = false)
+
+stepsize_callback = StepsizeCallback(cfl = 0.5)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback,
+                        alive_callback,
+                        amr_callback,
+                        stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
+            dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep = false, callback = callbacks);
+summary_callback() # print the timer summary
diff --git a/examples/t8code_2d_dgsem/elixir_advection_amr_unstructured_flag.jl b/examples/t8code_2d_dgsem/elixir_advection_amr_unstructured_flag.jl
index 9138586ccc..d889385481 100644
--- a/examples/t8code_2d_dgsem/elixir_advection_amr_unstructured_flag.jl
+++ b/examples/t8code_2d_dgsem/elixir_advection_amr_unstructured_flag.jl
@@ -61,7 +61,8 @@ amr_controller = ControllerThreeLevel(semi, IndicatorMax(semi, variable = first)
 amr_callback = AMRCallback(semi, amr_controller,
                            interval = 5,
                            adapt_initial_condition = true,
-                           adapt_initial_condition_only_refine = true)
+                           adapt_initial_condition_only_refine = true,
+                           dynamic_load_balancing = true)
 
 stepsize_callback = StepsizeCallback(cfl = 0.7)
 
diff --git a/examples/t8code_2d_dgsem/elixir_euler_weak_blast_wave_amr.jl b/examples/t8code_2d_dgsem/elixir_euler_weak_blast_wave_amr.jl
new file mode 100644
index 0000000000..a3366caa31
--- /dev/null
+++ b/examples/t8code_2d_dgsem/elixir_euler_weak_blast_wave_amr.jl
@@ -0,0 +1,114 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+
+equations = CompressibleEulerEquations2D(1.4)
+
+function initial_condition_weak_blast_wave(x, t, equations::CompressibleEulerEquations2D)
+    # Set up polar coordinates
+    inicenter = SVector(0.0, 0.0)
+    x_norm = x[1] - inicenter[1]
+    y_norm = x[2] - inicenter[2]
+    r = sqrt(x_norm^2 + y_norm^2)
+
+    r0 = 0.2
+    E = 1
+    p0_inner = 3
+    p0_outer = 1
+
+    # Calculate primitive variables
+    rho = 1.1
+    v1 = 0.0
+    v2 = 0.0
+    p = r > r0 ? p0_outer : p0_inner
+
+    return prim2cons(SVector(rho, v1, v2, p), equations)
+end
+
+initial_condition = initial_condition_weak_blast_wave
+
+# Get the DG approximation space
+
+# Activate the shock capturing + flux differencing
+surface_flux = flux_lax_friedrichs
+volume_flux = flux_ranocha
+polydeg = 4
+basis = LobattoLegendreBasis(polydeg)
+indicator_sc = IndicatorHennemannGassner(equations, basis,
+                                         alpha_max = 0.5,
+                                         alpha_min = 0.001,
+                                         alpha_smooth = true,
+                                         variable = density_pressure)
+volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
+                                                 volume_flux_dg = volume_flux,
+                                                 volume_flux_fv = surface_flux)
+
+solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux,
+               volume_integral = volume_integral)
+
+# Affine type mapping to take the [-1,1]^2 domain
+# and warp it as described in https://arxiv.org/abs/2012.12040
+# Warping with the coefficient 0.2 is even more extreme.
+function mapping_twist(xi, eta)
+    y = eta + 0.125 * cos(1.5 * pi * xi) * cos(0.5 * pi * eta)
+    x = xi + 0.125 * cos(0.5 * pi * xi) * cos(2.0 * pi * y)
+    return SVector(x, y)
+end
+
+# The mesh below can be made periodic
+# Create T8codeMesh with 8 x 8 trees
+trees_per_dimension = (8, 8)
+mesh = T8codeMesh(trees_per_dimension, polydeg = 4,
+                  mapping = mapping_twist,
+                  initial_refinement_level = 0,
+                  periodicity = true)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 1.2)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 400
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
+                                     save_analysis = true,
+                                     extra_analysis_errors = (:conservation_error,))
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+amr_indicator = IndicatorLöhner(semi, variable = Trixi.density)
+amr_controller = ControllerThreeLevel(semi, amr_indicator,
+                                      base_level = 0,
+                                      med_level = 1, med_threshold = 0.05,
+                                      max_level = 2, max_threshold = 0.1)
+amr_callback = AMRCallback(semi, amr_controller,
+                           interval = 5,
+                           adapt_initial_condition = true,
+                           adapt_initial_condition_only_refine = true)
+
+stepsize_callback = StepsizeCallback(cfl = 0.5)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback,
+                        alive_callback,
+                        amr_callback,
+                        stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
+            dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep = false, callback = callbacks);#, maxiters=4);
+summary_callback() # print the timer summary
+
+# Finalize `T8codeMesh` to make sure MPI related objects in t8code are
+# released before `MPI` finalizes.
+!isinteractive() && finalize(mesh)
diff --git a/examples/t8code_3d_dgsem/elixir_euler_weak_blast_wave_amr.jl b/examples/t8code_3d_dgsem/elixir_euler_weak_blast_wave_amr.jl
new file mode 100644
index 0000000000..106b482114
--- /dev/null
+++ b/examples/t8code_3d_dgsem/elixir_euler_weak_blast_wave_amr.jl
@@ -0,0 +1,116 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+
+equations = CompressibleEulerEquations3D(1.4)
+
+function initial_condition_weak_blast_wave(x, t,
+                                           equations::CompressibleEulerEquations3D)
+    # Set up polar coordinates
+    inicenter = SVector(0.0, 0.0, 0.0)
+    x_norm = x[1] - inicenter[1]
+    y_norm = x[2] - inicenter[2]
+    z_norm = x[3] - inicenter[3]
+    r = sqrt(x_norm^2 + y_norm^2 + z_norm^2)
+
+    r0 = 0.2
+    E = 1.0
+    p0_inner = 3
+    p0_outer = 1
+
+    # Calculate primitive variables
+    rho = 1.1
+    v1 = 0.0
+    v2 = 0.0
+    v3 = 0.0
+    p = r > r0 ? p0_outer : p0_inner
+
+    return prim2cons(SVector(rho, v1, v2, v3, p), equations)
+end
+
+initial_condition = initial_condition_weak_blast_wave
+
+surface_flux = flux_lax_friedrichs
+volume_flux = flux_ranocha
+polydeg = 4
+basis = LobattoLegendreBasis(polydeg)
+indicator_sc = IndicatorHennemannGassner(equations, basis,
+                                         alpha_max = 1.0,
+                                         alpha_min = 0.001,
+                                         alpha_smooth = true,
+                                         variable = density_pressure)
+volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
+                                                 volume_flux_dg = volume_flux,
+                                                 volume_flux_fv = surface_flux)
+
+solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux,
+               volume_integral = volume_integral)
+
+# Setup a periodic mesh with 4 x 4 x 4 trees and 8 x 8 x 8 elements
+trees_per_dimension = (4, 4, 4)
+
+function mapping_twist(xi, eta, zeta)
+    y = eta + 1 / 6 * (cos(1.5 * pi * xi) * cos(0.5 * pi * eta) * cos(0.5 * pi * zeta))
+
+    x = xi + 1 / 6 * (cos(0.5 * pi * xi) * cos(2 * pi * y) * cos(0.5 * pi * zeta))
+
+    z = zeta + 1 / 6 * (cos(0.5 * pi * x) * cos(pi * y) * cos(0.5 * pi * zeta))
+
+    return SVector(x, y, z)
+end
+
+mesh = T8codeMesh(trees_per_dimension,
+                  polydeg = 2,
+                  mapping = mapping_twist,
+                  initial_refinement_level = 1,
+                  periodicity = true)
+
+# Create the semidiscretization object
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 1.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
+                                     extra_analysis_errors = (:conservation_error,))
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+amr_indicator = IndicatorLöhner(semi, variable = Trixi.density)
+amr_controller = ControllerThreeLevel(semi, amr_indicator,
+                                      base_level = 1,
+                                      med_level = 2, med_threshold = 0.05,
+                                      max_level = 3, max_threshold = 0.15)
+amr_callback = AMRCallback(semi, amr_controller,
+                           interval = 1,
+                           adapt_initial_condition = false,
+                           adapt_initial_condition_only_refine = false)
+
+stepsize_callback = StepsizeCallback(cfl = 0.5)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback,
+                        alive_callback,
+                        amr_callback,
+                        stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
+            dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep = false, callback = callbacks);
+summary_callback() # print the timer summary
+
+# Finalize `T8codeMesh` to make sure MPI related objects in t8code are
+# released before `MPI` finalizes.
+!isinteractive() && finalize(mesh)
diff --git a/examples/tree_2d_dgsem/elixir_shallowwater_wall.jl b/examples/tree_2d_dgsem/elixir_shallowwater_wall.jl
index b8dbad5068..f8f601d412 100644
--- a/examples/tree_2d_dgsem/elixir_shallowwater_wall.jl
+++ b/examples/tree_2d_dgsem/elixir_shallowwater_wall.jl
@@ -29,8 +29,8 @@ boundary_condition = boundary_condition_slip_wall
 ###############################################################################
 # Get the DG approximation space
 
-volume_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal)
-surface_flux = (flux_lax_friedrichs, flux_nonconservative_ersing_etal)
+volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal)
+surface_flux = (flux_lax_friedrichs, flux_nonconservative_wintermeyer_etal)
 solver = DGSEM(polydeg = 3, surface_flux = surface_flux,
                volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
 
diff --git a/src/Trixi.jl b/src/Trixi.jl
index 23a8cfe1d0..90fa590c50 100644
--- a/src/Trixi.jl
+++ b/src/Trixi.jl
@@ -179,7 +179,6 @@ export flux, flux_central, flux_lax_friedrichs, flux_hll, flux_hllc, flux_hlle,
        flux_kennedy_gruber, flux_shima_etal, flux_ec,
        flux_fjordholm_etal, flux_nonconservative_fjordholm_etal,
        flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal,
-       flux_nonconservative_ersing_etal,
        flux_chan_etal, flux_nonconservative_chan_etal, flux_winters_etal,
        hydrostatic_reconstruction_audusse_etal, flux_nonconservative_audusse_etal,
        FluxPlusDissipation, DissipationGlobalLaxFriedrichs, DissipationLocalLaxFriedrichs,
diff --git a/src/callbacks_step/amr_dg2d.jl b/src/callbacks_step/amr_dg2d.jl
index 94524b23a3..062020b6d4 100644
--- a/src/callbacks_step/amr_dg2d.jl
+++ b/src/callbacks_step/amr_dg2d.jl
@@ -76,7 +76,10 @@ function rebalance_solver!(u_ode::AbstractVector, mesh::TreeMesh{2}, equations,
     end # GC.@preserve old_u_ode
 end
 
-# Refine elements in the DG solver based on a list of cell_ids that should be refined
+# Refine elements in the DG solver based on a list of cell_ids that should be refined.
+# On `P4estMesh`, if an element refines the solution scaled by the Jacobian `J*u` is interpolated
+# from the parent element into the four children elements. The solution on each child
+# element is then recovered by dividing by the new element Jacobians.
 function refine!(u_ode::AbstractVector, adaptor, mesh::Union{TreeMesh{2}, P4estMesh{2}},
                  equations, dg::DGSEM, cache, elements_to_refine)
     # Return early if there is nothing to do
@@ -96,9 +99,23 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::Union{TreeMesh{2}, P4estM
     # Retain current solution data
     old_n_elements = nelements(dg, cache)
     old_u_ode = copy(u_ode)
-    GC.@preserve old_u_ode begin # OBS! If we don't GC.@preserve old_u_ode, it might be GC'ed
+    old_inverse_jacobian = copy(cache.elements.inverse_jacobian)
+    # OBS! If we don't GC.@preserve old_u_ode and old_inverse_jacobian, they might be GC'ed
+    GC.@preserve old_u_ode old_inverse_jacobian begin
         old_u = wrap_array(old_u_ode, mesh, equations, dg, cache)
 
+        if mesh isa P4estMesh
+            # Loop over all elements in old container and scale the old solution by the Jacobian
+            # prior to projection
+            for old_element_id in 1:old_n_elements
+                for v in eachvariable(equations)
+                    old_u[v, .., old_element_id] .= (old_u[v, .., old_element_id] ./
+                                                     old_inverse_jacobian[..,
+                                                                          old_element_id])
+                end
+            end
+        end
+
         reinitialize_containers!(mesh, equations, dg, cache)
 
         resize!(u_ode,
@@ -112,10 +129,34 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::Union{TreeMesh{2}, P4estM
                 # Refine element and store solution directly in new data structure
                 refine_element!(u, element_id, old_u, old_element_id,
                                 adaptor, equations, dg)
+
+                if mesh isa P4estMesh
+                    # Before `element_id` is incremented, divide by the new Jacobians on each
+                    # child element and save the result
+                    for m in 0:3 # loop over the children
+                        for v in eachvariable(equations)
+                            u[v, .., element_id + m] .*= (0.25f0 .*
+                                                          cache.elements.inverse_jacobian[..,
+                                                                                          element_id + m])
+                        end
+                    end
+                end
+
+                # Increment `element_id` on the refined mesh with the number
+                # of children, i.e., 4 in 2D
                 element_id += 2^ndims(mesh)
             else
-                # Copy old element data to new element container
-                @views u[:, .., element_id] .= old_u[:, .., old_element_id]
+                if mesh isa P4estMesh
+                    # Copy old element data to new element container and remove Jacobian scaling
+                    for v in eachvariable(equations)
+                        u[v, .., element_id] .= (old_u[v, .., old_element_id] .*
+                                                 old_inverse_jacobian[..,
+                                                                      old_element_id])
+                    end
+                else # isa TreeMesh
+                    @views u[:, .., element_id] .= old_u[:, .., old_element_id]
+                end
+                # No refinement occurred, so increment `element_id` on the new mesh by one
                 element_id += 1
             end
         end
@@ -125,7 +166,7 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::Union{TreeMesh{2}, P4estM
         @assert element_id ==
                 nelements(dg, cache) +
                 1||element_id == nelements(dg, cache) + 2^ndims(mesh) "element_id = $element_id, nelements(dg, cache) = $(nelements(dg, cache))"
-    end # GC.@preserve old_u_ode
+    end # GC.@preserve old_u_ode old_inverse_jacobian
 
     # Sanity check
     if mesh isa TreeMesh && isperiodic(mesh.tree) && nmortars(cache.mortars) == 0 &&
@@ -228,7 +269,10 @@ function refine_element!(u::AbstractArray{<:Any, 4}, element_id,
     return nothing
 end
 
-# Coarsen elements in the DG solver based on a list of cell_ids that should be removed
+# Coarsen elements in the DG solver based on a list of cell_ids that should be removed.
+# On `P4estMesh`, if an element coarsens the solution scaled by the Jacobian `J*u` is projected
+# from the four children elements back onto the parent element. The solution on the parent
+# element is then recovered by dividing by the new element Jacobian.
 function coarsen!(u_ode::AbstractVector, adaptor,
                   mesh::Union{TreeMesh{2}, P4estMesh{2}},
                   equations, dg::DGSEM, cache, elements_to_remove)
@@ -246,12 +290,26 @@ function coarsen!(u_ode::AbstractVector, adaptor,
     to_be_removed = falses(nelements(dg, cache))
     to_be_removed[elements_to_remove] .= true
 
-    # Retain current solution data
+    # Retain current solution data and Jacobians
     old_n_elements = nelements(dg, cache)
     old_u_ode = copy(u_ode)
-    GC.@preserve old_u_ode begin # OBS! If we don't GC.@preserve old_u_ode, it might be GC'ed
+    old_inverse_jacobian = copy(cache.elements.inverse_jacobian)
+    # OBS! If we don't GC.@preserve old_u_ode and old_inverse_jacobian, they might be GC'ed
+    GC.@preserve old_u_ode old_inverse_jacobian begin
         old_u = wrap_array(old_u_ode, mesh, equations, dg, cache)
 
+        if mesh isa P4estMesh
+            # Loop over all elements in old container and scale the old solution by the Jacobian
+            # prior to projection
+            for old_element_id in 1:old_n_elements
+                for v in eachvariable(equations)
+                    old_u[v, .., old_element_id] .= (old_u[v, .., old_element_id] ./
+                                                     old_inverse_jacobian[..,
+                                                                          old_element_id])
+                end
+            end
+        end
+
         reinitialize_containers!(mesh, equations, dg, cache)
 
         resize!(u_ode,
@@ -277,17 +335,39 @@ function coarsen!(u_ode::AbstractVector, adaptor,
                 # Coarsen elements and store solution directly in new data structure
                 coarsen_elements!(u, element_id, old_u, old_element_id,
                                   adaptor, equations, dg)
+
+                if mesh isa P4estMesh
+                    # Before `element_id` is incremented, divide by the new Jacobian and save
+                    # the result in the parent element
+                    for v in eachvariable(equations)
+                        u[v, .., element_id] .*= (4 .*
+                                                  cache.elements.inverse_jacobian[..,
+                                                                                  element_id])
+                    end
+                end
+
+                # Increment `element_id` on the coarsened mesh by one and `skip` = 3 in 2D
+                # because 4 children elements become 1 parent element
                 element_id += 1
                 skip = 2^ndims(mesh) - 1
             else
-                # Copy old element data to new element container
-                @views u[:, .., element_id] .= old_u[:, .., old_element_id]
+                if mesh isa P4estMesh
+                    # Copy old element data to new element container and remove Jacobian scaling
+                    for v in eachvariable(equations)
+                        u[v, .., element_id] .= (old_u[v, .., old_element_id] .*
+                                                 old_inverse_jacobian[..,
+                                                                      old_element_id])
+                    end
+                else # isa TreeMesh
+                    @views u[:, .., element_id] .= old_u[:, .., old_element_id]
+                end
+                # No coarsening occurred, so increment `element_id` on the new mesh by one
                 element_id += 1
             end
         end
         # If everything is correct, we should have processed all elements.
         @assert element_id==nelements(dg, cache) + 1 "element_id = $element_id, nelements(dg, cache) = $(nelements(dg, cache))"
-    end # GC.@preserve old_u_ode
+    end # GC.@preserve old_u_ode old_inverse_jacobian
 
     # Sanity check
     if mesh isa TreeMesh && isperiodic(mesh.tree) && nmortars(cache.mortars) == 0 &&
@@ -404,16 +484,27 @@ function adapt!(u_ode::AbstractVector, adaptor, mesh::T8codeMesh{2}, equations,
     # Note: This is true for `quads`.
     T8_CHILDREN = 4
 
-    # Retain current solution data.
+    # Retain current solution and inverse Jacobian data.
     old_u_ode = copy(u_ode)
+    old_inverse_jacobian = copy(cache.elements.inverse_jacobian)
 
+    # OBS! If we don't GC.@preserve old_u_ode and old_inverse_jacobian, they might be GC'ed
     GC.@preserve old_u_ode begin
         old_u = wrap_array(old_u_ode, mesh, equations, dg, cache)
 
+        # Loop over all elements in old container and scale the old solution by the Jacobian
+        # prior to interpolation or projection
+        for old_element_id in 1:old_nelems
+            for v in eachvariable(equations)
+                old_u[v, .., old_element_id] .= (old_u[v, .., old_element_id] ./
+                                                 old_inverse_jacobian[..,
+                                                                      old_element_id])
+            end
+        end
+
         reinitialize_containers!(mesh, equations, dg, cache)
 
-        resize!(u_ode,
-                nvariables(equations) * nnodes(dg)^ndims(mesh) * nelements(dg, cache))
+        resize!(u_ode, nvariables(equations) * ndofs(mesh, dg, cache))
         u = wrap_array(u_ode, mesh, equations, dg, cache)
 
         while old_index <= old_nelems && new_index <= new_nelems
@@ -422,6 +513,18 @@ function adapt!(u_ode::AbstractVector, adaptor, mesh::T8codeMesh{2}, equations,
                 # Refine element and store solution directly in new data structure.
                 refine_element!(u, new_index, old_u, old_index, adaptor, equations, dg)
 
+                # Before indices are incremented divide by the new Jacobians on each
+                # child element and save the result
+                for m in 0:3 # loop over the children
+                    for v in eachvariable(equations)
+                        u[v, .., new_index + m] .*= (0.25f0 .*
+                                                     cache.elements.inverse_jacobian[..,
+                                                                                     new_index + m])
+                    end
+                end
+
+                # Increment `old_index` on the original mesh and the `new_index`
+                # on the refined mesh with the number of children, i.e., T8_CHILDREN = 4
                 old_index += 1
                 new_index += T8_CHILDREN
 
@@ -436,19 +539,34 @@ function adapt!(u_ode::AbstractVector, adaptor, mesh::T8codeMesh{2}, equations,
                 coarsen_elements!(u, new_index, old_u, old_index, adaptor, equations,
                                   dg)
 
+                # Before the indices are incremented divide by the new Jacobian and save
+                # the result again in the parent element
+                for v in eachvariable(equations)
+                    u[v, .., new_index] .*= (4 .* cache.elements.inverse_jacobian[..,
+                                                                             new_index])
+                end
+
+                # Increment `old_index` on the original mesh with the number of children
+                # (T8_CHILDREN = 4 in 2D) and the `new_index` by one for the single
+                # coarsened element
                 old_index += T8_CHILDREN
                 new_index += 1
 
             else # No changes.
 
-                # Copy old element data to new element container.
-                @views u[:, .., new_index] .= old_u[:, .., old_index]
+                # Copy old element data to new element container and remove Jacobian scaling
+                for v in eachvariable(equations)
+                    u[v, .., new_index] .= (old_u[v, .., old_index] .*
+                                            old_inverse_jacobian[.., old_index])
+                end
 
+                # No refinement / coarsening occurred, so increment element index
+                # on each mesh by one
                 old_index += 1
                 new_index += 1
             end
         end # while
-    end # GC.@preserve old_u_ode
+    end # GC.@preserve old_u_ode old_inverse_jacobian
 
     return nothing
 end
diff --git a/src/callbacks_step/amr_dg3d.jl b/src/callbacks_step/amr_dg3d.jl
index 3f67951baf..594c30dcca 100644
--- a/src/callbacks_step/amr_dg3d.jl
+++ b/src/callbacks_step/amr_dg3d.jl
@@ -5,9 +5,11 @@
 @muladd begin
 #! format: noindent
 
-# Refine elements in the DG solver based on a list of cell_ids that should be refined
-function refine!(u_ode::AbstractVector, adaptor,
-                 mesh::Union{TreeMesh{3}, P4estMesh{3}},
+# Refine elements in the DG solver based on a list of cell_ids that should be refined.
+# On `P4estMesh`, if an element refines the solution scaled by the Jacobian `J*u` is interpolated
+# from the parent element into the eight children elements. The solution on each child
+# element is then recovered by dividing by the new element Jacobians.
+function refine!(u_ode::AbstractVector, adaptor, mesh::Union{TreeMesh{3}, P4estMesh{3}},
                  equations, dg::DGSEM, cache, elements_to_refine)
     # Return early if there is nothing to do
     if isempty(elements_to_refine)
@@ -26,9 +28,23 @@ function refine!(u_ode::AbstractVector, adaptor,
     # Retain current solution data
     old_n_elements = nelements(dg, cache)
     old_u_ode = copy(u_ode)
-    GC.@preserve old_u_ode begin # OBS! If we don't GC.@preserve old_u_ode, it might be GC'ed
+    old_inverse_jacobian = copy(cache.elements.inverse_jacobian)
+    # OBS! If we don't GC.@preserve old_u_ode and old_inverse_jacobian, they might be GC'ed
+    GC.@preserve old_u_ode old_inverse_jacobian begin
         old_u = wrap_array(old_u_ode, mesh, equations, dg, cache)
 
+        if mesh isa P4estMesh
+            # Loop over all elements in old container and scale the old solution by the Jacobian
+            # prior to projection
+            for old_element_id in 1:old_n_elements
+                for v in eachvariable(equations)
+                    old_u[v, .., old_element_id] .= (old_u[v, .., old_element_id] ./
+                                                     old_inverse_jacobian[..,
+                                                                          old_element_id])
+                end
+            end
+        end
+
         reinitialize_containers!(mesh, equations, dg, cache)
 
         resize!(u_ode,
@@ -44,12 +60,36 @@ function refine!(u_ode::AbstractVector, adaptor,
         for old_element_id in 1:old_n_elements
             if needs_refinement[old_element_id]
                 # Refine element and store solution directly in new data structure
-                refine_element!(u, element_id, old_u, old_element_id,
-                                adaptor, equations, dg, u_tmp1, u_tmp2)
+                refine_element!(u, element_id, old_u, old_element_id, adaptor,
+                                equations, dg, u_tmp1, u_tmp2)
+
+                if mesh isa P4estMesh
+                    # Before `element_id` is incremented, divide by the new Jacobians on each
+                    # child element and save the result
+                    for m in 0:7 # loop over the children
+                        for v in eachvariable(equations)
+                            u[v, .., element_id + m] .*= (0.125f0 .*
+                                                          cache.elements.inverse_jacobian[..,
+                                                                                          element_id + m])
+                        end
+                    end
+                end
+
+                # Increment `element_id` on the refined mesh with the number
+                # of children, i.e., 8 in 3D
                 element_id += 2^ndims(mesh)
             else
-                # Copy old element data to new element container
-                @views u[:, .., element_id] .= old_u[:, .., old_element_id]
+                if mesh isa P4estMesh
+                    # Copy old element data to new element container and remove Jacobian scaling
+                    for v in eachvariable(equations)
+                        u[v, .., element_id] .= (old_u[v, .., old_element_id] .*
+                                                 old_inverse_jacobian[..,
+                                                                      old_element_id])
+                    end
+                else # isa TreeMesh
+                    @views u[:, .., element_id] .= old_u[:, .., old_element_id]
+                end
+                # No refinement occurred, so increment `element_id` on the new mesh by one
                 element_id += 1
             end
         end
@@ -59,7 +99,7 @@ function refine!(u_ode::AbstractVector, adaptor,
         @assert element_id ==
                 nelements(dg, cache) +
                 1||element_id == nelements(dg, cache) + 2^ndims(mesh) "element_id = $element_id, nelements(dg, cache) = $(nelements(dg, cache))"
-    end # GC.@preserve old_u_ode
+    end # GC.@preserve old_u_ode old_inverse_jacobian
 
     # Sanity check
     if mesh isa TreeMesh && isperiodic(mesh.tree) && nmortars(cache.mortars) == 0
@@ -145,7 +185,10 @@ function refine_element!(u::AbstractArray{<:Any, 5}, element_id,
     return nothing
 end
 
-# Coarsen elements in the DG solver based on a list of cell_ids that should be removed
+# Coarsen elements in the DG solver based on a list of cell_ids that should be removed.
+# On `P4estMesh`, if an element coarsens the solution scaled by the Jacobian `J*u` is projected
+# from the eight children elements back onto the parent element. The solution on the parent
+# element is then recovered by dividing by the new element Jacobian.
 function coarsen!(u_ode::AbstractVector, adaptor,
                   mesh::Union{TreeMesh{3}, P4estMesh{3}},
                   equations, dg::DGSEM, cache, elements_to_remove)
@@ -163,12 +206,26 @@ function coarsen!(u_ode::AbstractVector, adaptor,
     to_be_removed = falses(nelements(dg, cache))
     to_be_removed[elements_to_remove] .= true
 
-    # Retain current solution data
+    # Retain current solution data and Jacobians
     old_n_elements = nelements(dg, cache)
     old_u_ode = copy(u_ode)
-    GC.@preserve old_u_ode begin # OBS! If we don't GC.@preserve old_u_ode, it might be GC'ed
+    old_inverse_jacobian = copy(cache.elements.inverse_jacobian)
+    # OBS! If we don't GC.@preserve old_u_ode and old_inverse_jacobian, they might be GC'ed
+    GC.@preserve old_u_ode old_inverse_jacobian begin
         old_u = wrap_array(old_u_ode, mesh, equations, dg, cache)
 
+        if mesh isa P4estMesh
+            # Loop over all elements in old container and scale the old solution by the Jacobian
+            # prior to projection
+            for old_element_id in 1:old_n_elements
+                for v in eachvariable(equations)
+                    old_u[v, .., old_element_id] .= (old_u[v, .., old_element_id] ./
+                                                     old_inverse_jacobian[..,
+                                                                          old_element_id])
+                end
+            end
+        end
+
         reinitialize_containers!(mesh, equations, dg, cache)
 
         resize!(u_ode,
@@ -196,19 +253,41 @@ function coarsen!(u_ode::AbstractVector, adaptor,
                 @assert all(to_be_removed[old_element_id:(old_element_id + 2^ndims(mesh) - 1)]) "bad cell/element order"
 
                 # Coarsen elements and store solution directly in new data structure
-                coarsen_elements!(u, element_id, old_u, old_element_id,
-                                  adaptor, equations, dg, u_tmp1, u_tmp2)
+                coarsen_elements!(u, element_id, old_u, old_element_id, adaptor,
+                                  equations, dg, u_tmp1, u_tmp2)
+
+                if mesh isa P4estMesh
+                    # Before `element_id` is incremented, divide by the new Jacobian and save
+                    # the result in the parent element
+                    for v in eachvariable(equations)
+                        u[v, .., element_id] .*= (8 .*
+                                                  cache.elements.inverse_jacobian[..,
+                                                                                  element_id])
+                    end
+                end
+
+                # Increment `element_id` on the coarsened mesh by one and `skip` = 7 in 3D
+                # because 8 children elements become 1 parent element
                 element_id += 1
                 skip = 2^ndims(mesh) - 1
             else
-                # Copy old element data to new element container
-                @views u[:, .., element_id] .= old_u[:, .., old_element_id]
+                if mesh isa P4estMesh
+                    # Copy old element data to new element container and remove Jacobian scaling
+                    for v in eachvariable(equations)
+                        u[v, .., element_id] .= (old_u[v, .., old_element_id] .*
+                                                 old_inverse_jacobian[..,
+                                                                      old_element_id])
+                    end
+                else # isa TreeMesh
+                    @views u[:, .., element_id] .= old_u[:, .., old_element_id]
+                end
+                # No coarsening occurred, so increment `element_id` on the new mesh by one
                 element_id += 1
             end
         end
         # If everything is correct, we should have processed all elements.
         @assert element_id==nelements(dg, cache) + 1 "element_id = $element_id, nelements(dg, cache) = $(nelements(dg, cache))"
-    end # GC.@preserve old_u_ode
+    end # GC.@preserve old_u_ode old_inverse_jacobian
 
     # Sanity check
     if mesh isa TreeMesh && isperiodic(mesh.tree) && nmortars(cache.mortars) == 0
@@ -335,16 +414,27 @@ function adapt!(u_ode::AbstractVector, adaptor, mesh::T8codeMesh{3}, equations,
     # Note: This is only true for `hexs`.
     T8_CHILDREN = 8
 
-    # Retain current solution data.
+    # Retain current solution and inverse Jacobian data.
     old_u_ode = copy(u_ode)
+    old_inverse_jacobian = copy(cache.elements.inverse_jacobian)
 
+    # OBS! If we don't GC.@preserve old_u_ode and old_inverse_jacobian, they might be GC'ed
     GC.@preserve old_u_ode begin
         old_u = wrap_array(old_u_ode, mesh, equations, dg, cache)
 
+        # Loop over all elements in old container and scale the old solution by the Jacobian
+        # prior to interpolation or projection
+        for old_element_id in 1:old_nelems
+            for v in eachvariable(equations)
+                old_u[v, .., old_element_id] .= (old_u[v, .., old_element_id] ./
+                                                 old_inverse_jacobian[..,
+                                                                      old_element_id])
+            end
+        end
+
         reinitialize_containers!(mesh, equations, dg, cache)
 
-        resize!(u_ode,
-                nvariables(equations) * ndofs(mesh, dg, cache))
+        resize!(u_ode, nvariables(equations) * ndofs(mesh, dg, cache))
         u = wrap_array(u_ode, mesh, equations, dg, cache)
 
         u_tmp1 = Array{eltype(u), 4}(undef, nvariables(equations), nnodes(dg),
@@ -359,6 +449,18 @@ function adapt!(u_ode::AbstractVector, adaptor, mesh::T8codeMesh{3}, equations,
                 refine_element!(u, new_index, old_u, old_index, adaptor, equations, dg,
                                 u_tmp1, u_tmp2)
 
+                # Before indices are incremented divide by the new Jacobians on each
+                # child element and save the result
+                for m in 0:7 # loop over the children
+                    for v in eachvariable(equations)
+                        u[v, .., new_index + m] .*= (0.125f0 .*
+                                                     cache.elements.inverse_jacobian[..,
+                                                                                     new_index + m])
+                    end
+                end
+
+                # Increment `old_index` on the original mesh and the `new_index`
+                # on the refined mesh with the number of children, i.e., T8_CHILDREN = 8
                 old_index += 1
                 new_index += T8_CHILDREN
 
@@ -373,19 +475,34 @@ function adapt!(u_ode::AbstractVector, adaptor, mesh::T8codeMesh{3}, equations,
                 coarsen_elements!(u, new_index, old_u, old_index, adaptor, equations,
                                   dg, u_tmp1, u_tmp2)
 
+                # Before the indices are incremented divide by the new Jacobian and save
+                # the result again in the parent element
+                for v in eachvariable(equations)
+                    u[v, .., new_index] .*= (8 .* cache.elements.inverse_jacobian[..,
+                                                                             new_index])
+                end
+
+                # Increment `old_index` on the original mesh with the number of children
+                # (T8_CHILDREN = 8 in 3D) and the `new_index` by one for the single
+                # coarsened element
                 old_index += T8_CHILDREN
                 new_index += 1
 
             else # No changes.
 
-                # Copy old element data to new element container.
-                @views u[:, .., new_index] .= old_u[:, .., old_index]
+                # Copy old element data to new element container and remove Jacobian scaling
+                for v in eachvariable(equations)
+                    u[v, .., new_index] .= (old_u[v, .., old_index] .*
+                                            old_inverse_jacobian[.., old_index])
+                end
 
+                # No refinement / coarsening occurred, so increment element index
+                # on each mesh by one
                 old_index += 1
                 new_index += 1
             end
         end # while
-    end # GC.@preserve old_u_ode
+    end # GC.@preserve old_u_ode old_inverse_jacobian
 
     return nothing
 end
diff --git a/src/equations/polytropic_euler_2d.jl b/src/equations/polytropic_euler_2d.jl
index e900fd6407..571d4723f6 100644
--- a/src/equations/polytropic_euler_2d.jl
+++ b/src/equations/polytropic_euler_2d.jl
@@ -62,7 +62,7 @@ in combination with [`source_terms_convergence_test`](@ref).
 function initial_condition_convergence_test(x, t, equations::PolytropicEulerEquations2D)
     # manufactured initial condition from Winters (2019) [0.1007/s10543-019-00789-w]
     # domain must be set to [0, 1] x [0, 1]
-    h = 8 + cos(2 * pi * x[1]) * sin(2 * pi * x[2]) * cos(2 * pi * t)
+    h = 8 + cospi(2 * x[1]) * sinpi(2 * x[2]) * cospi(2 * t)
 
     return SVector(h, h / 2, 3 * h / 2)
 end
@@ -78,19 +78,20 @@ Source terms used for convergence tests in combination with
     rho, v1, v2 = cons2prim(u, equations)
 
     # Residual from Winters (2019) [0.1007/s10543-019-00789-w] eq. (5.2).
-    h = 8 + cos(2 * pi * x[1]) * sin(2 * pi * x[2]) * cos(2 * pi * t)
-    h_t = -2 * pi * cos(2 * pi * x[1]) * sin(2 * pi * x[2]) * sin(2 * pi * t)
-    h_x = -2 * pi * sin(2 * pi * x[1]) * sin(2 * pi * x[2]) * cos(2 * pi * t)
-    h_y = 2 * pi * cos(2 * pi * x[1]) * cos(2 * pi * x[2]) * cos(2 * pi * t)
+    RealT = eltype(u)
+    h = 8 + cospi(2 * x[1]) * sinpi(2 * x[2]) * cospi(2 * t)
+    h_t = -2 * convert(RealT, pi) * cospi(2 * x[1]) * sinpi(2 * x[2]) * sinpi(2 * t)
+    h_x = -2 * convert(RealT, pi) * sinpi(2 * x[1]) * sinpi(2 * x[2]) * cospi(2 * t)
+    h_y = 2 * convert(RealT, pi) * cospi(2 * x[1]) * cospi(2 * x[2]) * cospi(2 * t)
 
     rho_x = h_x
     rho_y = h_y
 
     b = equations.kappa * equations.gamma * h^(equations.gamma - 1)
 
-    r_1 = h_t + h_x / 2 + 3 / 2 * h_y
-    r_2 = h_t / 2 + h_x / 4 + b * rho_x + 3 / 4 * h_y
-    r_3 = 3 / 2 * h_t + 3 / 4 * h_x + 9 / 4 * h_y + b * rho_y
+    r_1 = h_t + h_x / 2 + 3 * h_y / 2
+    r_2 = h_t / 2 + h_x / 4 + b * rho_x + 3 * h_y / 4
+    r_3 = 3 * h_t / 2 + 3 * h_x / 4 + 9 * h_y / 4 + b * rho_y
 
     return SVector(r_1, r_2, r_3)
 end
@@ -113,9 +114,10 @@ function initial_condition_weak_blast_wave(x, t, equations::PolytropicEulerEquat
     phi = atan(y_norm, x_norm)
 
     # Calculate primitive variables
-    rho = r > 0.5 ? 1.0 : 1.1691
-    v1 = r > 0.5 ? 0.0 : 0.1882 * cos(phi)
-    v2 = r > 0.5 ? 0.0 : 0.1882 * sin(phi)
+    RealT = eltype(x)
+    rho = r > 0.5f0 ? one(RealT) : convert(RealT, 1.1691)
+    v1 = r > 0.5f0 ? zero(RealT) : convert(RealT, 0.1882) * cos(phi)
+    v2 = r > 0.5f0 ? zero(RealT) : convert(RealT, 0.1882) * sin(phi)
 
     return prim2cons(SVector(rho, v1, v2), equations)
 end
@@ -181,17 +183,17 @@ For details see Section 3.2 of the following reference
     v_dot_n_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2]
 
     # Compute the necessary mean values
-    if equations.gamma == 1.0 # isothermal gas
+    if equations.gamma == 1 # isothermal gas
         rho_mean = ln_mean(rho_ll, rho_rr)
     else # equations.gamma > 1 # polytropic gas
         rho_mean = stolarsky_mean(rho_ll, rho_rr, equations.gamma)
     end
-    v1_avg = 0.5 * (v1_ll + v1_rr)
-    v2_avg = 0.5 * (v2_ll + v2_rr)
-    p_avg = 0.5 * (p_ll + p_rr)
+    v1_avg = 0.5f0 * (v1_ll + v1_rr)
+    v2_avg = 0.5f0 * (v2_ll + v2_rr)
+    p_avg = 0.5f0 * (p_ll + p_rr)
 
     # Calculate fluxes depending on normal_direction
-    f1 = rho_mean * 0.5 * (v_dot_n_ll + v_dot_n_rr)
+    f1 = rho_mean * 0.5f0 * (v_dot_n_ll + v_dot_n_rr)
     f2 = f1 * v1_avg + p_avg * normal_direction[1]
     f3 = f1 * v2_avg + p_avg * normal_direction[2]
 
@@ -207,21 +209,21 @@ end
     p_rr = equations.kappa * rho_rr^equations.gamma
 
     # Compute the necessary mean values
-    if equations.gamma == 1.0 # isothermal gas
+    if equations.gamma == 1 # isothermal gas
         rho_mean = ln_mean(rho_ll, rho_rr)
     else # equations.gamma > 1 # polytropic gas
         rho_mean = stolarsky_mean(rho_ll, rho_rr, equations.gamma)
     end
-    v1_avg = 0.5 * (v1_ll + v1_rr)
-    v2_avg = 0.5 * (v2_ll + v2_rr)
-    p_avg = 0.5 * (p_ll + p_rr)
+    v1_avg = 0.5f0 * (v1_ll + v1_rr)
+    v2_avg = 0.5f0 * (v2_ll + v2_rr)
+    p_avg = 0.5f0 * (p_ll + p_rr)
 
     if orientation == 1 # x-direction
-        f1 = rho_mean * 0.5 * (v1_ll + v1_rr)
+        f1 = rho_mean * 0.5f0 * (v1_ll + v1_rr)
         f2 = f1 * v1_avg + p_avg
         f3 = f1 * v2_avg
     else # y-direction
-        f1 = rho_mean * 0.5 * (v2_ll + v2_rr)
+        f1 = rho_mean * 0.5f0 * (v2_ll + v2_rr)
         f2 = f1 * v1_avg
         f3 = f1 * v2_avg + p_avg
     end
@@ -360,14 +362,14 @@ end
     v_square = v1^2 + v2^2
     p = pressure(u, equations)
     # Form of the internal energy depends on gas type
-    if equations.gamma == 1.0 # isothermal gas
+    if equations.gamma == 1 # isothermal gas
         internal_energy = equations.kappa * log(rho)
     else # equations.gamma > 1 # polytropic gas
         internal_energy = equations.kappa * rho^(equations.gamma - 1) /
-                          (equations.gamma - 1.0)
+                          (equations.gamma - 1)
     end
 
-    w1 = internal_energy + p / rho - 0.5 * v_square
+    w1 = internal_energy + p / rho - 0.5f0 * v_square
     w2 = v1
     w3 = v2
 
diff --git a/src/equations/shallow_water_1d.jl b/src/equations/shallow_water_1d.jl
index 998deb04de..3c218eee9d 100644
--- a/src/equations/shallow_water_1d.jl
+++ b/src/equations/shallow_water_1d.jl
@@ -201,20 +201,27 @@ end
 Non-symmetric two-point volume flux discretizing the nonconservative (source) term
 that contains the gradient of the bottom topography [`ShallowWaterEquations1D`](@ref).
 
-Further details are available in the paper:
+Gives entropy conservation and well-balancedness on both the volume and surface when combined with 
+[`flux_wintermeyer_etal`](@ref).
+
+Further details are available in the papers:
 - Niklas Wintermeyer, Andrew R. Winters, Gregor J. Gassner and David A. Kopriva (2017)
   An entropy stable nodal discontinuous Galerkin method for the two dimensional
   shallow water equations on unstructured curvilinear meshes with discontinuous bathymetry
   [DOI: 10.1016/j.jcp.2017.03.036](https://doi.org/10.1016/j.jcp.2017.03.036)
+- Patrick Ersing, Andrew R. Winters (2023)
+  An entropy stable discontinuous Galerkin method for the two-layer shallow water equations on 
+  curvilinear meshes
+  [DOI: 10.48550/arXiv.2306.12699](https://doi.org/10.48550/arXiv.2306.12699)
 """
 @inline function flux_nonconservative_wintermeyer_etal(u_ll, u_rr, orientation::Integer,
                                                        equations::ShallowWaterEquations1D)
     # Pull the necessary left and right state information
     h_ll = waterheight(u_ll, equations)
-    b_rr = u_rr[3]
+    b_jump = u_rr[3] - u_ll[3]
 
     # Bottom gradient nonconservative term: (0, g h b_x, 0)
-    f = SVector(0, equations.gravity * h_ll * b_rr, 0)
+    f = SVector(0, equations.gravity * h_ll * b_jump, 0)
 
     return f
 end
@@ -226,11 +233,8 @@ end
 Non-symmetric two-point surface flux discretizing the nonconservative (source) term of
 that contains the gradient of the bottom topography [`ShallowWaterEquations1D`](@ref).
 
-This contains additional terms compared to [`flux_nonconservative_wintermeyer_etal`](@ref)
-that account for possible discontinuities in the bottom topography function.
-Thus, this flux should be used in general at interfaces. For flux differencing volume terms,
-[`flux_nonconservative_wintermeyer_etal`](@ref) is analytically equivalent but slightly
-cheaper.
+This flux can be used together with [`flux_fjordholm_etal`](@ref) at interfaces to ensure entropy
+conservation and well-balancedness.
 
 Further details for the original finite volume formulation are available in
 - Ulrik S. Fjordholm, Siddhartha Mishr and Eitan Tadmor (2011)
@@ -256,7 +260,6 @@ and for curvilinear 2D case in the paper:
     #       cross-averaging across a discontinuous bottom topography
     #  (ii) True surface part that uses `h_average` and `b_jump` to handle discontinuous bathymetry
     f = SVector(0,
-                equations.gravity * h_ll * b_ll +
                 equations.gravity * h_average * b_jump,
                 0)
 
@@ -285,7 +288,7 @@ Further details on the hydrostatic reconstruction and its motivation can be foun
                                                    orientation::Integer,
                                                    equations::ShallowWaterEquations1D)
     # Pull the water height and bottom topography on the left
-    h_ll, _, b_ll = u_ll
+    h_ll, _, _ = u_ll
 
     # Create the hydrostatic reconstruction for the left solution state
     u_ll_star, _ = hydrostatic_reconstruction_audusse_etal(u_ll, u_rr, equations)
@@ -293,52 +296,11 @@ Further details on the hydrostatic reconstruction and its motivation can be foun
     # Copy the reconstructed water height for easier to read code
     h_ll_star = u_ll_star[1]
 
-    # Includes two parts:
-    #   (i)  Diagonal (consistent) term from the volume flux that uses `b_ll` to avoid
-    #        cross-averaging across a discontinuous bottom topography
-    #   (ii) True surface part that uses `h_ll` and `h_ll_star` to handle discontinuous bathymetry
     return SVector(0,
-                   equations.gravity * h_ll * b_ll +
                    equations.gravity * (h_ll^2 - h_ll_star^2),
                    0)
 end
 
-"""
-    flux_nonconservative_ersing_etal(u_ll, u_rr, orientation::Integer,
-                                     equations::ShallowWaterEquations1D)
-
-!!! warning "Experimental code"
-    This numerical flux is experimental and may change in any future release.
-
-Non-symmetric path-conservative two-point volume flux discretizing the nonconservative (source) term
-that contains the gradient of the bottom topography [`ShallowWaterEquations1D`](@ref).
-
-This is a modified version of [`flux_nonconservative_wintermeyer_etal`](@ref) that gives entropy 
-conservation and well-balancedness in both the volume and surface when combined with 
-[`flux_wintermeyer_etal`](@ref).
-
-For further details see:
-- Patrick Ersing, Andrew R. Winters (2023)
-  An entropy stable discontinuous Galerkin method for the two-layer shallow water equations on 
-  curvilinear meshes
-  [DOI: 10.48550/arXiv.2306.12699](https://doi.org/10.48550/arXiv.2306.12699)
-"""
-@inline function flux_nonconservative_ersing_etal(u_ll, u_rr, orientation::Integer,
-                                                  equations::ShallowWaterEquations1D)
-    # Pull the necessary left and right state information
-    h_ll = waterheight(u_ll, equations)
-    b_rr = u_rr[3]
-    b_ll = u_ll[3]
-
-    # Calculate jump
-    b_jump = b_rr - b_ll
-
-    # Bottom gradient nonconservative term: (0, g h b_x, 0)
-    f = SVector(0, equations.gravity * h_ll * b_jump, 0)
-
-    return f
-end
-
 """
     flux_fjordholm_etal(u_ll, u_rr, orientation,
                         equations::ShallowWaterEquations1D)
@@ -378,7 +340,8 @@ end
 
 Total energy conservative (mathematical entropy for shallow water equations) split form.
 When the bottom topography is nonzero this scheme will be well-balanced when used as a `volume_flux`.
-The `surface_flux` should still use, e.g., [`flux_fjordholm_etal`](@ref).
+For the `surface_flux` either [`flux_wintermeyer_etal`](@ref) or [`flux_fjordholm_etal`](@ref) can
+be used to ensure well-balancedness and entropy conservation.
 
 Further details are available in Theorem 1 of the paper:
 - Niklas Wintermeyer, Andrew R. Winters, Gregor J. Gassner and David A. Kopriva (2017)
@@ -571,7 +534,7 @@ end
 # Note, only the first two are the entropy variables, the third entry still
 # just carries the bottom topography values for convenience
 @inline function cons2entropy(u, equations::ShallowWaterEquations1D)
-    h, h_v, b = u
+    h, _, b = u
 
     v = velocity(u, equations)
 
diff --git a/src/equations/shallow_water_2d.jl b/src/equations/shallow_water_2d.jl
index db8f00fc15..4ecaf3b6e1 100644
--- a/src/equations/shallow_water_2d.jl
+++ b/src/equations/shallow_water_2d.jl
@@ -274,28 +274,30 @@ end
 Non-symmetric two-point volume flux discretizing the nonconservative (source) term
 that contains the gradient of the bottom topography [`ShallowWaterEquations2D`](@ref).
 
-On curvilinear meshes, this nonconservative flux depends on both the
-contravariant vector (normal direction) at the current node and the averaged
-one. This is different from numerical fluxes used to discretize conservative
-terms.
+For the `surface_flux` either [`flux_wintermeyer_etal`](@ref) or [`flux_fjordholm_etal`](@ref) can
+be used to ensure well-balancedness and entropy conservation.
 
-Further details are available in the paper:
+Further details are available in the papers:
 - Niklas Wintermeyer, Andrew R. Winters, Gregor J. Gassner and David A. Kopriva (2017)
   An entropy stable nodal discontinuous Galerkin method for the two dimensional
   shallow water equations on unstructured curvilinear meshes with discontinuous bathymetry
   [DOI: 10.1016/j.jcp.2017.03.036](https://doi.org/10.1016/j.jcp.2017.03.036)
+- Patrick Ersing, Andrew R. Winters (2023)
+  An entropy stable discontinuous Galerkin method for the two-layer shallow water equations on 
+  curvilinear meshes
+  [DOI: 10.48550/arXiv.2306.12699](https://doi.org/10.48550/arXiv.2306.12699)
 """
 @inline function flux_nonconservative_wintermeyer_etal(u_ll, u_rr, orientation::Integer,
                                                        equations::ShallowWaterEquations2D)
     # Pull the necessary left and right state information
     h_ll = waterheight(u_ll, equations)
-    b_rr = u_rr[4]
+    b_jump = u_rr[4] - u_ll[4]
 
     # Bottom gradient nonconservative term: (0, g h b_x, g h b_y, 0)
     if orientation == 1
-        f = SVector(0, equations.gravity * h_ll * b_rr, 0, 0)
+        f = SVector(0, equations.gravity * h_ll * b_jump, 0, 0)
     else # orientation == 2
-        f = SVector(0, 0, equations.gravity * h_ll * b_rr, 0)
+        f = SVector(0, 0, equations.gravity * h_ll * b_jump, 0)
     end
     return f
 end
@@ -306,12 +308,12 @@ end
                                                        equations::ShallowWaterEquations2D)
     # Pull the necessary left and right state information
     h_ll = waterheight(u_ll, equations)
-    b_rr = u_rr[4]
-    # Note this routine only uses the `normal_direction_average` and the average of the
-    # bottom topography to get a quadratic split form DG gradient on curved elements
+    b_jump = u_rr[4] - u_ll[4]
+
+    # Bottom gradient nonconservative term: (0, g h b_x, g h b_y, 0)
     return SVector(0,
-                   normal_direction_average[1] * equations.gravity * h_ll * b_rr,
-                   normal_direction_average[2] * equations.gravity * h_ll * b_rr,
+                   normal_direction_average[1] * equations.gravity * h_ll * b_jump,
+                   normal_direction_average[2] * equations.gravity * h_ll * b_jump,
                    0)
 end
 
@@ -326,16 +328,8 @@ end
 Non-symmetric two-point surface flux discretizing the nonconservative (source) term of
 that contains the gradient of the bottom topography [`ShallowWaterEquations2D`](@ref).
 
-On curvilinear meshes, this nonconservative flux depends on both the
-contravariant vector (normal direction) at the current node and the averaged
-one. This is different from numerical fluxes used to discretize conservative
-terms.
-
-This contains additional terms compared to [`flux_nonconservative_wintermeyer_etal`](@ref)
-that account for possible discontinuities in the bottom topography function.
-Thus, this flux should be used in general at interfaces. For flux differencing volume terms,
-[`flux_nonconservative_wintermeyer_etal`](@ref) is analytically equivalent but slightly
-cheaper.
+This flux can be used together with [`flux_fjordholm_etal`](@ref) at interfaces to ensure entropy
+conservation and well-balancedness.
 
 Further details for the original finite volume formulation are available in
 - Ulrik S. Fjordholm, Siddhartha Mishr and Eitan Tadmor (2011)
@@ -356,18 +350,13 @@ and for curvilinear 2D case in the paper:
     h_average = 0.5f0 * (h_ll + h_rr)
     b_jump = b_rr - b_ll
 
-    # Includes two parts:
-    #   (i)  Diagonal (consistent) term from the volume flux that uses `b_ll` to avoid
-    #        cross-averaging across a discontinuous bottom topography
-    #   (ii) True surface part that uses `h_average` and `b_jump` to handle discontinuous bathymetry
+    # Bottom gradient nonconservative term: (0, g h b_x, g h b_y, 0)
     if orientation == 1
         f = SVector(0,
-                    equations.gravity * h_ll * b_ll +
                     equations.gravity * h_average * b_jump,
                     0, 0)
     else # orientation == 2
         f = SVector(0, 0,
-                    equations.gravity * h_ll * b_ll +
                     equations.gravity * h_average * b_jump,
                     0)
     end
@@ -383,20 +372,12 @@ end
     h_ll, _, _, b_ll = u_ll
     h_rr, _, _, b_rr = u_rr
 
-    # Comes in two parts:
-    #   (i)  Diagonal (consistent) term from the volume flux that uses `normal_direction_average`
-    #        but we use `b_ll` to avoid cross-averaging across a discontinuous bottom topography
-
-    f2 = normal_direction_average[1] * equations.gravity * h_ll * b_ll
-    f3 = normal_direction_average[2] * equations.gravity * h_ll * b_ll
-
-    #   (ii) True surface part that uses `normal_direction_ll`, `h_average` and `b_jump`
-    #        to handle discontinuous bathymetry
     h_average = 0.5f0 * (h_ll + h_rr)
     b_jump = b_rr - b_ll
 
-    f2 += normal_direction_ll[1] * equations.gravity * h_average * b_jump
-    f3 += normal_direction_ll[2] * equations.gravity * h_average * b_jump
+    # Bottom gradient nonconservative term: (0, g h b_x, g h b_y, 0)
+    f2 = normal_direction_average[1] * equations.gravity * h_average * b_jump
+    f3 = normal_direction_average[2] * equations.gravity * h_average * b_jump
 
     # First and last equations do not have a nonconservative flux
     f1 = f4 = 0
@@ -472,18 +453,12 @@ Further details for the hydrostatic reconstruction and its motivation can be fou
     # Copy the reconstructed water height for easier to read code
     h_ll_star = u_ll_star[1]
 
-    # Includes two parts:
-    #   (i)  Diagonal (consistent) term from the volume flux that uses `b_ll` to avoid
-    #        cross-averaging across a discontinuous bottom topography
-    #   (ii) True surface part that uses `h_ll` and `h_ll_star` to handle discontinuous bathymetry
     if orientation == 1
         f = SVector(0,
-                    equations.gravity * h_ll * b_ll +
                     equations.gravity * (h_ll^2 - h_ll_star^2),
                     0, 0)
     else # orientation == 2
         f = SVector(0, 0,
-                    equations.gravity * h_ll * b_ll +
                     equations.gravity * (h_ll^2 - h_ll_star^2),
                     0)
     end
@@ -504,18 +479,8 @@ end
     # Copy the reconstructed water height for easier to read code
     h_ll_star = u_ll_star[1]
 
-    # Comes in two parts:
-    #   (i)  Diagonal (consistent) term from the volume flux that uses `normal_direction_average`
-    #        but we use `b_ll` to avoid cross-averaging across a discontinuous bottom topography
-
-    f2 = normal_direction_average[1] * equations.gravity * h_ll * b_ll
-    f3 = normal_direction_average[2] * equations.gravity * h_ll * b_ll
-
-    #   (ii) True surface part that uses `normal_direction_ll`, `h_ll` and `h_ll_star`
-    #        to handle discontinuous bathymetry
-
-    f2 += normal_direction_ll[1] * equations.gravity * (h_ll^2 - h_ll_star^2)
-    f3 += normal_direction_ll[2] * equations.gravity * (h_ll^2 - h_ll_star^2)
+    f2 = normal_direction_average[1] * equations.gravity * (h_ll^2 - h_ll_star^2)
+    f3 = normal_direction_average[2] * equations.gravity * (h_ll^2 - h_ll_star^2)
 
     # First and last equations do not have a nonconservative flux
     f1 = f4 = 0
@@ -523,73 +488,6 @@ end
     return SVector(f1, f2, f3, f4)
 end
 
-"""
-    flux_nonconservative_ersing_etal(u_ll, u_rr, orientation::Integer,
-                                     equations::ShallowWaterEquations2D)
-    flux_nonconservative_ersing_etal(u_ll, u_rr,
-                                     normal_direction_ll::AbstractVector,
-                                     normal_direction_average::AbstractVector,
-                                     equations::ShallowWaterEquations2D)
-
-!!! warning "Experimental code"
-    This numerical flux is experimental and may change in any future release.
-
-Non-symmetric path-conservative two-point volume flux discretizing the nonconservative (source) term
-that contains the gradient of the bottom topography [`ShallowWaterEquations2D`](@ref).
-
-On curvilinear meshes, this nonconservative flux depends on both the
-contravariant vector (normal direction) at the current node and the averaged
-one. This is different from numerical fluxes used to discretize conservative
-terms.
-
-This is a modified version of [`flux_nonconservative_wintermeyer_etal`](@ref) that gives entropy 
-conservation and well-balancedness in both the volume and surface when combined with 
-[`flux_wintermeyer_etal`](@ref).
-
-For further details see:
-- Patrick Ersing, Andrew R. Winters (2023)
-  An entropy stable discontinuous Galerkin method for the two-layer shallow water equations on 
-  curvilinear meshes
-  [DOI: 10.48550/arXiv.2306.12699](https://doi.org/10.48550/arXiv.2306.12699)
-"""
-@inline function flux_nonconservative_ersing_etal(u_ll, u_rr, orientation::Integer,
-                                                  equations::ShallowWaterEquations2D)
-    # Pull the necessary left and right state information
-    h_ll = waterheight(u_ll, equations)
-    b_rr = u_rr[4]
-    b_ll = u_ll[4]
-
-    # Calculate jump
-    b_jump = b_rr - b_ll
-
-    # Bottom gradient nonconservative term: (0, g h b_x, g h b_y, 0)
-    if orientation == 1
-        f = SVector(0, equations.gravity * h_ll * b_jump, 0, 0)
-    else # orientation == 2
-        f = SVector(0, 0, equations.gravity * h_ll * b_jump, 0)
-    end
-    return f
-end
-
-@inline function flux_nonconservative_ersing_etal(u_ll, u_rr,
-                                                  normal_direction_ll::AbstractVector,
-                                                  normal_direction_average::AbstractVector,
-                                                  equations::ShallowWaterEquations2D)
-    # Pull the necessary left and right state information
-    h_ll = waterheight(u_ll, equations)
-    b_rr = u_rr[4]
-    b_ll = u_ll[4]
-
-    # Calculate jump
-    b_jump = b_rr - b_ll
-    # Note this routine only uses the `normal_direction_average` and the average of the
-    # bottom topography to get a quadratic split form DG gradient on curved elements
-    return SVector(0,
-                   normal_direction_average[1] * equations.gravity * h_ll * b_jump,
-                   normal_direction_average[2] * equations.gravity * h_ll * b_jump,
-                   0)
-end
-
 """
     flux_fjordholm_etal(u_ll, u_rr, orientation_or_normal_direction,
                         equations::ShallowWaterEquations2D)
@@ -664,7 +562,8 @@ end
 
 Total energy conservative (mathematical entropy for shallow water equations) split form.
 When the bottom topography is nonzero this scheme will be well-balanced when used as a `volume_flux`.
-The `surface_flux` should still use, e.g., [`flux_fjordholm_etal`](@ref).
+For the `surface_flux` either [`flux_wintermeyer_etal`](@ref) or [`flux_fjordholm_etal`](@ref) can
+be used to ensure well-balancedness and entropy conservation.
 
 Further details are available in Theorem 1 of the paper:
 - Niklas Wintermeyer, Andrew R. Winters, Gregor J. Gassner and David A. Kopriva (2017)
diff --git a/src/meshes/structured_mesh.jl b/src/meshes/structured_mesh.jl
index 553aabbbc2..9e5c55197d 100644
--- a/src/meshes/structured_mesh.jl
+++ b/src/meshes/structured_mesh.jl
@@ -255,13 +255,15 @@ function correction_term_3d(x, y, z, faces)
                              linear_interpolate(x, faces[1](1, z), faces[2](1, z)) +
                              linear_interpolate(z, faces[5](x, 1), faces[6](x, 1)))
 
-    # Correction for x-terms
+    # Correction for z-terms
     c_z = linear_interpolate(z,
                              linear_interpolate(x, faces[1](y, -1), faces[2](y, -1)) +
                              linear_interpolate(y, faces[3](x, -1), faces[4](x, -1)),
                              linear_interpolate(x, faces[1](y, 1), faces[2](y, 1)) +
                              linear_interpolate(y, faces[3](x, 1), faces[4](x, 1)))
 
+    # Each of the 12 edges are counted twice above
+    # so we divide the correction term by two
     return 0.5 * (c_x + c_y + c_z)
 end
 
diff --git a/src/meshes/surface_interpolant.jl b/src/meshes/surface_interpolant.jl
index c3ee330a39..dec62d6a97 100644
--- a/src/meshes/surface_interpolant.jl
+++ b/src/meshes/surface_interpolant.jl
@@ -52,7 +52,7 @@ function derivative_at(s, boundary_curve::CurvedSurface)
            y_coordinate_at_s_on_boundary_curve_prime
 end
 
-# Chebysehv-Gauss-Lobatto nodes and weights for use with curved boundaries
+# Chebyshev-Gauss-Lobatto nodes and weights for use with curved boundaries
 function chebyshev_gauss_lobatto_nodes_weights(n_nodes::Integer)
 
     # Initialize output
diff --git a/test/test_mpi_p4est_2d.jl b/test/test_mpi_p4est_2d.jl
index 6d66bc68a2..29de4efc6d 100644
--- a/test/test_mpi_p4est_2d.jl
+++ b/test/test_mpi_p4est_2d.jl
@@ -97,8 +97,8 @@ const EXAMPLES_DIR = pkgdir(Trixi, "examples", "p4est_2d_dgsem")
     @trixi_testset "elixir_advection_amr_unstructured_flag.jl" begin
         @test_trixi_include(joinpath(EXAMPLES_DIR,
                                      "elixir_advection_amr_unstructured_flag.jl"),
-                            l2=[0.0012766060609964525],
-                            linf=[0.01750280631586159],
+                            l2=[0.0012808538770535593],
+                            linf=[0.01752690016659812],
                             coverage_override=(maxiters = 6,))
 
         # Ensure that we do not have excessive memory allocations
diff --git a/test/test_mpi_p4est_3d.jl b/test/test_mpi_p4est_3d.jl
index cca9093ec5..4f9465b85d 100644
--- a/test/test_mpi_p4est_3d.jl
+++ b/test/test_mpi_p4est_3d.jl
@@ -69,8 +69,8 @@ const EXAMPLES_DIR = pkgdir(Trixi, "examples", "p4est_3d_dgsem")
     @trixi_testset "elixir_advection_amr_unstructured_curved.jl" begin
         @test_trixi_include(joinpath(EXAMPLES_DIR,
                                      "elixir_advection_amr_unstructured_curved.jl"),
-                            l2=[1.6236411810065552e-5],
-                            linf=[0.0010554006923731395],
+                            l2=[1.6163120948209677e-5],
+                            linf=[0.0010572201890564834],
                             tspan=(0.0, 1.0),
                             coverage_override=(maxiters = 6,
                                                initial_refinement_level = 0,
diff --git a/test/test_mpi_t8code_2d.jl b/test/test_mpi_t8code_2d.jl
index 7c7fc03898..75e65c8c38 100644
--- a/test/test_mpi_t8code_2d.jl
+++ b/test/test_mpi_t8code_2d.jl
@@ -97,8 +97,9 @@ const EXAMPLES_DIR = pkgdir(Trixi, "examples", "t8code_2d_dgsem")
     @trixi_testset "elixir_advection_amr_unstructured_flag.jl" begin
         @test_trixi_include(joinpath(EXAMPLES_DIR,
                                      "elixir_advection_amr_unstructured_flag.jl"),
-                            l2=[0.001980652042312077],
-                            linf=[0.0328882442132265],
+                            l2=[0.002019623611753929],
+                            linf=[0.03542375961299987],
+                            dynamic_load_balancing=false,
                             coverage_override=(maxiters = 6,))
 
         # Ensure that we do not have excessive memory allocations
diff --git a/test/test_mpi_t8code_3d.jl b/test/test_mpi_t8code_3d.jl
index a15690a762..2e837f79ad 100644
--- a/test/test_mpi_t8code_3d.jl
+++ b/test/test_mpi_t8code_3d.jl
@@ -69,8 +69,8 @@ const EXAMPLES_DIR = pkgdir(Trixi, "examples", "t8code_3d_dgsem")
     @trixi_testset "elixir_advection_amr_unstructured_curved.jl" begin
         @test_trixi_include(joinpath(EXAMPLES_DIR,
                                      "elixir_advection_amr_unstructured_curved.jl"),
-                            l2=[2.0556575425846923e-5],
-                            linf=[0.00105682693484822],
+                            l2=[2.0535121347526814e-5],
+                            linf=[0.0010586603797777504],
                             tspan=(0.0, 1.0),
                             coverage_override=(maxiters = 6,
                                                initial_refinement_level = 0,
diff --git a/test/test_p4est_2d.jl b/test/test_p4est_2d.jl
index ef5ac2c9de..f56dc9cd76 100644
--- a/test/test_p4est_2d.jl
+++ b/test/test_p4est_2d.jl
@@ -78,8 +78,8 @@ end
 @trixi_testset "elixir_advection_amr_unstructured_flag.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR,
                                  "elixir_advection_amr_unstructured_flag.jl"),
-                        l2=[0.0012766060609964525],
-                        linf=[0.01750280631586159],
+                        l2=[0.0012808538770535593],
+                        linf=[0.01752690016659812],
                         coverage_override=(maxiters = 6,))
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
@@ -301,16 +301,16 @@ end
 @trixi_testset "elixir_euler_blast_wave_amr.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_blast_wave_amr.jl"),
                         l2=[
-                            6.32183914e-01,
-                            3.86914231e-01,
-                            3.86869171e-01,
-                            1.06575688e+00,
+                            0.6321850210104147,
+                            0.38691446170269167,
+                            0.3868695626809587,
+                            1.0657553825683956,
                         ],
                         linf=[
-                            2.76020890e+00,
-                            2.32659890e+00,
-                            2.32580837e+00,
-                            2.15778188e+00,
+                            2.7602280007469666,
+                            2.3265993814913672,
+                            2.3258078438689673,
+                            2.1577683028925416,
                         ],
                         tspan=(0.0, 0.3),
                         coverage_override=(maxiters = 6,))
@@ -327,16 +327,16 @@ end
 @trixi_testset "elixir_euler_wall_bc_amr.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_wall_bc_amr.jl"),
                         l2=[
-                            0.020291447969983396,
-                            0.017479614254319948,
-                            0.011387644425613437,
-                            0.0514420126021293,
+                            0.02026685991647352,
+                            0.017467584076280237,
+                            0.011378371604813321,
+                            0.05138942558296091,
                         ],
                         linf=[
-                            0.3582779022370579,
-                            0.32073537890751663,
-                            0.221818049107692,
-                            0.9209559420400415,
+                            0.35924402060711524,
+                            0.32068389566068806,
+                            0.2361141752119986,
+                            0.9289840057748628,
                         ],
                         tspan=(0.0, 0.15))
     # Ensure that we do not have excessive memory allocations
@@ -352,16 +352,16 @@ end
 @trixi_testset "elixir_euler_forward_step_amr.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_forward_step_amr.jl"),
                         l2=[
-                            0.004194875320833303,
-                            0.003785140699353966,
-                            0.0013696609105790351,
-                            0.03265268616046424,
+                            0.004191480950848891,
+                            0.003781298410569231,
+                            0.0013470418422981045,
+                            0.03262817609394949,
                         ],
                         linf=[
-                            2.0585399781442852,
-                            2.213428805506876,
-                            3.862362410419163,
-                            17.75187237459251,
+                            2.0581500751947113,
+                            2.2051301367971288,
+                            3.8502467979250254,
+                            17.750333649853616,
                         ],
                         tspan=(0.0, 0.0001),
                         rtol=1.0e-7,
@@ -409,16 +409,16 @@ end
 @trixi_testset "elixir_euler_supersonic_cylinder.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_supersonic_cylinder.jl"),
                         l2=[
-                            0.026798021911954406,
-                            0.05118546368109259,
-                            0.03206703583774831,
-                            0.19680026567208672,
+                            0.02676082999794676,
+                            0.05110830068968181,
+                            0.03205164257040607,
+                            0.1965981012724311,
                         ],
                         linf=[
-                            3.653905721692421,
-                            4.285035711361009,
-                            6.8544353186357645,
-                            31.748244912257533,
+                            3.6830683476364476,
+                            4.284442685012427,
+                            6.857777546171545,
+                            31.749285097390576,
                         ],
                         tspan=(0.0, 0.001),
                         skip_coverage=true)
@@ -529,16 +529,17 @@ end
 
 @trixi_testset "elixir_mhd_rotor.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_rotor.jl"),
-                        l2=[0.4552084651735862, 0.8918048264575757, 0.832471223081887,
+                        l2=[0.4552094211937344, 0.8918052934760807, 0.8324715234680768,
                             0.0,
-                            0.9801264164951583, 0.10475690769435382, 0.1555132409149897,
+                            0.9801268321975978, 0.10475722739111007,
+                            0.15551326369033164,
                             0.0,
-                            2.0597079362763556e-5],
-                        linf=[10.194181233788775, 18.25472397868819, 10.031307436191334,
+                            2.0602990858239632e-5],
+                        linf=[10.19421969147307, 18.254409357804683, 10.031954811332596,
                             0.0,
-                            19.647239392277378, 1.3938810140985936, 1.8724965294853084,
+                            19.646870938371492, 1.3938679692894465, 1.8725058401937984,
                             0.0,
-                            0.0016290067532561904],
+                            0.0016201762010257296],
                         tspan=(0.0, 0.02))
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
@@ -617,12 +618,12 @@ end
     @test_trixi_include(joinpath(EXAMPLES_DIR,
                                  "elixir_euler_NACA0012airfoil_mach085.jl"),
                         l2=[
-                            5.634402680811982e-7, 6.748066107517321e-6,
-                            1.091879472416885e-5, 0.0006686372064029146,
+                            5.56114097044427e-7, 6.62284247153255e-6,
+                            1.0823259724601275e-5, 0.000659804574787503,
                         ],
                         linf=[
-                            0.0021456247890772823, 0.03957142889488085,
-                            0.03832024233032798, 2.6628739573358495,
+                            0.002157589754528455, 0.039163189253511164,
+                            0.038386804399707625, 2.6685831417913914,
                         ],
                         amr_interval=1,
                         base_level=0, med_level=1, max_level=1,
@@ -652,8 +653,8 @@ end
     lift = Trixi.analyze(lift_coefficient, du, u, tspan[2], mesh, equations, solver,
                          semi.cache, semi)
 
-    @test isapprox(lift, 0.029076443678087403, atol = 1e-13)
-    @test isapprox(drag, 0.13564720009197903, atol = 1e-13)
+    @test isapprox(lift, 0.029094009322876882, atol = 1e-13)
+    @test isapprox(drag, 0.13579200776643238, atol = 1e-13)
 end
 
 @trixi_testset "elixir_euler_blast_wave_pure_fv.jl" begin
@@ -684,6 +685,40 @@ end
         @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
     end
 end
+
+@trixi_testset "elixir_euler_weak_blast_wave_amr.jl" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_weak_blast_wave_amr.jl"),
+                        l2=[
+                            0.11134260363848127,
+                            0.11752357091804219,
+                            0.11829112104640764,
+                            0.7557891142955036,
+                        ],
+                        linf=[
+                            0.5728647031475109,
+                            0.8353132977670252,
+                            0.8266797080712205,
+                            3.9792506230548317,
+                        ],
+                        tspan=(0.0, 0.1),
+                        coverage_override=(maxiters = 6,))
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
+    let
+        t = sol.t[end]
+        u_ode = sol.u[end]
+        du_ode = similar(u_ode)
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+    end
+    # Check for conservation
+    state_integrals = Trixi.integrate(sol.u[2], semi)
+    initial_state_integrals = analysis_callback.affect!.initial_state_integrals
+
+    @test isapprox(state_integrals[1], initial_state_integrals[1], atol = 1e-13)
+    @test isapprox(state_integrals[2], initial_state_integrals[2], atol = 1e-13)
+    @test isapprox(state_integrals[3], initial_state_integrals[3], atol = 1e-13)
+    @test isapprox(state_integrals[4], initial_state_integrals[4], atol = 1e-13)
+end
 end
 
 # Clean up afterwards: delete Trixi.jl output directory
diff --git a/test/test_p4est_3d.jl b/test/test_p4est_3d.jl
index 3c3bc1b2ab..3432bd69b2 100644
--- a/test/test_p4est_3d.jl
+++ b/test/test_p4est_3d.jl
@@ -78,8 +78,8 @@ end
 @trixi_testset "elixir_advection_amr_unstructured_curved.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR,
                                  "elixir_advection_amr_unstructured_curved.jl"),
-                        l2=[1.6236411810065552e-5],
-                        linf=[0.0010554006923731395],
+                        l2=[1.6163120948209677e-5],
+                        linf=[0.0010572201890564834],
                         tspan=(0.0, 1.0),
                         coverage_override=(maxiters = 6, initial_refinement_level = 0,
                                            base_level = 0, med_level = 1, max_level = 2))
@@ -571,16 +571,16 @@ end
 
 @trixi_testset "elixir_mhd_shockcapturing_amr.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_shockcapturing_amr.jl"),
-                        l2=[0.006298541670176575, 0.0064368506652601265,
-                            0.007108729762852636, 0.006530420607206385,
-                            0.02061185869237284, 0.005562033787605515,
-                            0.007571716276627825, 0.005571862660453231,
-                            3.909755063709152e-6],
-                        linf=[0.20904054009050665, 0.18622917151105936,
-                            0.2347957890323218, 0.19432508025509926,
-                            0.6858860133405615, 0.15172116633332622,
-                            0.22432820727833747, 0.16805989780225183,
-                            0.000535219040687628],
+                        l2=[0.006297229188299052, 0.0064363477630573936,
+                            0.007109134822960387, 0.0065295379843073945,
+                            0.02061487028361094, 0.005561406556868266,
+                            0.007570747563219415, 0.005571060186624124,
+                            3.910359570546058e-6],
+                        linf=[0.20904050617411984, 0.18630026905465372,
+                            0.23476537952044518, 0.19430178061639747,
+                            0.6858488631108304, 0.15169972134884624,
+                            0.22431157069631724, 0.16823638724229162,
+                            0.0005352202836463904],
                         tspan=(0.0, 0.04),
                         coverage_override=(maxiters = 6, initial_refinement_level = 1,
                                            base_level = 1, max_level = 2))
@@ -615,6 +615,43 @@ end
         @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
     end
 end
+
+@trixi_testset "elixir_euler_weak_blast_wave_amr.jl" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_weak_blast_wave_amr.jl"),
+                        l2=[
+                            0.011345993108796831,
+                            0.018525073963833696,
+                            0.019102348105917946,
+                            0.01920515438943838,
+                            0.15060493968460148,
+                        ],
+                        linf=[
+                            0.2994949779783401,
+                            0.5530175050084679,
+                            0.5335803757792128,
+                            0.5647252867336123,
+                            3.6462732329242566,
+                        ],
+                        tspan=(0.0, 0.025),
+                        coverage_override=(maxiters = 6,))
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
+    let
+        t = sol.t[end]
+        u_ode = sol.u[end]
+        du_ode = similar(u_ode)
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+    end
+    # Check for conservation
+    state_integrals = Trixi.integrate(sol.u[2], semi)
+    initial_state_integrals = analysis_callback.affect!.initial_state_integrals
+
+    @test isapprox(state_integrals[1], initial_state_integrals[1], atol = 1e-13)
+    @test isapprox(state_integrals[2], initial_state_integrals[2], atol = 1e-13)
+    @test isapprox(state_integrals[3], initial_state_integrals[3], atol = 1e-13)
+    @test isapprox(state_integrals[4], initial_state_integrals[4], atol = 1e-13)
+    @test isapprox(state_integrals[5], initial_state_integrals[5], atol = 1e-13)
+end
 end
 
 # Clean up afterwards: delete Trixi.jl output directory
diff --git a/test/test_t8code_2d.jl b/test/test_t8code_2d.jl
index b63d2a105a..c1fcc35521 100644
--- a/test/test_t8code_2d.jl
+++ b/test/test_t8code_2d.jl
@@ -121,8 +121,8 @@ end
 @trixi_testset "elixir_advection_amr_unstructured_flag.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR,
                                  "elixir_advection_amr_unstructured_flag.jl"),
-                        l2=[0.001993165013217687],
-                        linf=[0.032891018571625796],
+                        l2=[0.002019623611753929],
+                        linf=[0.03542375961299987],
                         coverage_override=(maxiters = 6,))
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
@@ -305,16 +305,16 @@ end
     # This test is identical to the one in `test_p4est_2d.jl` besides minor
     # deviations in the expected error norms.
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_rotor.jl"),
-                        l2=[0.44211360369891683, 0.8805178316216257, 0.8262710688468049,
+                        l2=[0.44207324634847545, 0.8804644301177857, 0.8262542320669426,
                             0.0,
-                            0.9616090460973586, 0.10386643568745411,
-                            0.15403457366543802, 0.0,
-                            2.8399715649715473e-5],
-                        linf=[10.04369305341599, 17.995640564998403, 9.576041548174265,
+                            0.9615023124189027, 0.10386709616755131, 0.1540308191628843,
                             0.0,
-                            19.429658884314534, 1.3821395681242314, 1.818559351543182,
+                            2.8350276854372125e-5],
+                        linf=[10.04548675437385, 17.998696852394836, 9.575802136190026,
                             0.0,
-                            0.002261930217575465],
+                            19.431290746184473, 1.3821685018474321, 1.8186235976551453,
+                            0.0,
+                            0.002309422702635547],
                         tspan=(0.0, 0.02))
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
@@ -325,6 +325,40 @@ end
         @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
     end
 end
+
+@trixi_testset "elixir_euler_weak_blast_wave_amr.jl" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_weak_blast_wave_amr.jl"),
+                        l2=[
+                            0.10823279736983638,
+                            0.1158152939803735,
+                            0.11633970342992006,
+                            0.751152651902375,
+                        ],
+                        linf=[
+                            0.5581611332828653,
+                            0.8354026029724041,
+                            0.834485181423738,
+                            3.923553028014343,
+                        ],
+                        tspan=(0.0, 0.1),
+                        coverage_override=(maxiters = 6,))
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
+    let
+        t = sol.t[end]
+        u_ode = sol.u[end]
+        du_ode = similar(u_ode)
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+    end
+    # Check for conservation
+    state_integrals = Trixi.integrate(sol.u[2], semi)
+    initial_state_integrals = analysis_callback.affect!.initial_state_integrals
+
+    @test isapprox(state_integrals[1], initial_state_integrals[1], atol = 1e-13)
+    @test isapprox(state_integrals[2], initial_state_integrals[2], atol = 1e-13)
+    @test isapprox(state_integrals[3], initial_state_integrals[3], atol = 1e-13)
+    @test isapprox(state_integrals[4], initial_state_integrals[4], atol = 1e-13)
+end
 end
 
 # Clean up afterwards: delete Trixi.jl output directory
diff --git a/test/test_t8code_3d.jl b/test/test_t8code_3d.jl
index 5c452444be..940d2c4337 100644
--- a/test/test_t8code_3d.jl
+++ b/test/test_t8code_3d.jl
@@ -61,7 +61,7 @@ mkdir(outdir)
         @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_nonconforming.jl"),
                             l2=[0.00253595715323843],
                             linf=[0.016486952252155795])
-        # Ensure that we do not have excessive memory allocations 
+        # Ensure that we do not have excessive memory allocations
         # (e.g., from type instabilities)
         let
             t = sol.t[end]
@@ -80,7 +80,7 @@ mkdir(outdir)
                             linf=[0.0007889950196294793],
                             coverage_override=(maxiters = 6, initial_refinement_level = 1,
                                                base_level = 1, med_level = 2, max_level = 3))
-        # Ensure that we do not have excessive memory allocations 
+        # Ensure that we do not have excessive memory allocations
         # (e.g., from type instabilities)
         let
             t = sol.t[end]
@@ -95,12 +95,12 @@ mkdir(outdir)
     @trixi_testset "elixir_advection_amr_unstructured_curved.jl" begin
         @test_trixi_include(joinpath(EXAMPLES_DIR,
                                      "elixir_advection_amr_unstructured_curved.jl"),
-                            l2=[2.0556575425846923e-5],
-                            linf=[0.00105682693484822],
+                            l2=[2.0535121347526814e-5],
+                            linf=[0.0010586603797777504],
                             tspan=(0.0, 1.0),
                             coverage_override=(maxiters = 6, initial_refinement_level = 0,
                                                base_level = 0, med_level = 1, max_level = 2))
-        # Ensure that we do not have excessive memory allocations 
+        # Ensure that we do not have excessive memory allocations
         # (e.g., from type instabilities)
         let
             t = sol.t[end]
@@ -129,7 +129,7 @@ mkdir(outdir)
                                 0.008526972236273522,
                             ],
                             tspan=(0.0, 0.01))
-        # Ensure that we do not have excessive memory allocations 
+        # Ensure that we do not have excessive memory allocations
         # (e.g., from type instabilities)
         let
             t = sol.t[end]
@@ -158,7 +158,7 @@ mkdir(outdir)
                                 0.01562861968368434,
                             ],
                             tspan=(0.0, 1.0))
-        # Ensure that we do not have excessive memory allocations 
+        # Ensure that we do not have excessive memory allocations
         # (e.g., from type instabilities)
         let
             t = sol.t[end]
@@ -186,7 +186,7 @@ mkdir(outdir)
                                 9.412914891981927e-12,
                             ],
                             tspan=(0.0, 0.03))
-        # Ensure that we do not have excessive memory allocations 
+        # Ensure that we do not have excessive memory allocations
         # (e.g., from type instabilities)
         let
             t = sol.t[end]
@@ -214,7 +214,7 @@ mkdir(outdir)
                                 9.592326932761353e-13,
                             ],
                             tspan=(0.0, 0.1), atol=5.0e-13,)
-        # Ensure that we do not have excessive memory allocations 
+        # Ensure that we do not have excessive memory allocations
         # (e.g., from type instabilities)
         let
             t = sol.t[end]
@@ -243,7 +243,7 @@ mkdir(outdir)
                             ],
                             tspan=(0.0, 0.2),
                             coverage_override=(polydeg = 3,)) # Prevent long compile time in CI
-        # Ensure that we do not have excessive memory allocations 
+        # Ensure that we do not have excessive memory allocations
         # (e.g., from type instabilities)
         let
             t = sol.t[end]
@@ -273,7 +273,7 @@ mkdir(outdir)
                             ],
                             tspan=(0.0, 0.3),
                             coverage_override=(polydeg = 3,)) # Prevent long compile time in CI
-        # Ensure that we do not have excessive memory allocations 
+        # Ensure that we do not have excessive memory allocations
         # (e.g., from type instabilities)
         let
             t = sol.t[end]
@@ -316,6 +316,43 @@ mkdir(outdir)
             @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
         end
     end
+
+    @trixi_testset "elixir_euler_weak_blast_wave_amr.jl" begin
+        @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_weak_blast_wave_amr.jl"),
+                            l2=[
+                                0.010014531529951328,
+                                0.0176268986746271,
+                                0.01817514447099777,
+                                0.018271085903740675,
+                                0.15193033077438198,
+                            ],
+                            linf=[
+                                0.2898958869606375,
+                                0.529717119064458,
+                                0.5567193302705906,
+                                0.570663236219957,
+                                3.5496520808512027,
+                            ],
+                            tspan=(0.0, 0.025),
+                            coverage_override=(maxiters = 6,))
+        # Ensure that we do not have excessive memory allocations
+        # (e.g., from type instabilities)
+        let
+            t = sol.t[end]
+            u_ode = sol.u[end]
+            du_ode = similar(u_ode)
+            @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+        end
+        # Check for conservation
+        state_integrals = Trixi.integrate(sol.u[2], semi)
+        initial_state_integrals = analysis_callback.affect!.initial_state_integrals
+
+        @test isapprox(state_integrals[1], initial_state_integrals[1], atol = 1e-13)
+        @test isapprox(state_integrals[2], initial_state_integrals[2], atol = 1e-13)
+        @test isapprox(state_integrals[3], initial_state_integrals[3], atol = 1e-13)
+        @test isapprox(state_integrals[4], initial_state_integrals[4], atol = 1e-13)
+        @test isapprox(state_integrals[5], initial_state_integrals[5], atol = 1e-13)
+    end
 end
 
 # Clean up afterwards: delete Trixi.jl output directory
diff --git a/test/test_tree_1d_shallowwater.jl b/test/test_tree_1d_shallowwater.jl
index 8fe3291a93..42a91e578e 100644
--- a/test/test_tree_1d_shallowwater.jl
+++ b/test/test_tree_1d_shallowwater.jl
@@ -98,7 +98,7 @@ end
     end
 end
 
-@trixi_testset "elixir_shallowwater_well_balanced.jl with flux_nonconservative_ersing_etal" begin
+@trixi_testset "elixir_shallowwater_well_balanced.jl with flux_nonconservative_wintermeyer_etal" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_well_balanced.jl"),
                         l2=[
                             0.10416666834254838,
@@ -107,9 +107,7 @@ end
                         ],
                         linf=[2.0000000000000004, 3.0610625110157164e-14, 2.0],
                         surface_flux=(flux_wintermeyer_etal,
-                                      flux_nonconservative_ersing_etal),
-                        volume_flux=(flux_wintermeyer_etal,
-                                     flux_nonconservative_ersing_etal),
+                                      flux_nonconservative_wintermeyer_etal),
                         tspan=(0.0, 0.25))
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
@@ -169,7 +167,7 @@ end
     end
 end
 
-@trixi_testset "elixir_shallowwater_source_terms.jl with flux_nonconservative_ersing_etal" begin
+@trixi_testset "elixir_shallowwater_source_terms.jl with flux_nonconservative_wintermeyer_etal" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"),
                         l2=[
                             0.005774284062933275,
@@ -182,9 +180,7 @@ end
                             9.098379777450205e-5,
                         ],
                         surface_flux=(flux_wintermeyer_etal,
-                                      flux_nonconservative_ersing_etal),
-                        volume_flux=(flux_wintermeyer_etal,
-                                     flux_nonconservative_ersing_etal),
+                                      flux_nonconservative_wintermeyer_etal),
                         tspan=(0.0, 0.025))
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
@@ -200,13 +196,13 @@ end
     @test_trixi_include(joinpath(EXAMPLES_DIR,
                                  "elixir_shallowwater_source_terms_dirichlet.jl"),
                         l2=[
-                            0.0022851099219788917,
-                            0.01560453773635554,
-                            4.43649172558535e-5,
+                            0.0022667320585353927,
+                            0.01571629729279524,
+                            4.4364917255842716e-5,
                         ],
                         linf=[
-                            0.008934615705174398,
-                            0.059403169140869405,
+                            0.008945234652224965,
+                            0.059403165802872415,
                             9.098379777405796e-5,
                         ],
                         tspan=(0.0, 0.025))
@@ -224,13 +220,13 @@ end
     @test_trixi_include(joinpath(EXAMPLES_DIR,
                                  "elixir_shallowwater_source_terms_dirichlet.jl"),
                         l2=[
-                            0.0022956052733432287,
-                            0.015540053559855601,
-                            4.43649172558535e-5,
+                            0.0022774071143995952,
+                            0.01566214422689219,
+                            4.4364917255842716e-5,
                         ],
                         linf=[
-                            0.008460440313118323,
-                            0.05720939349382359,
+                            0.008451721489057373,
+                            0.05720939055279217,
                             9.098379777405796e-5,
                         ],
                         surface_flux=(FluxHydrostaticReconstruction(FluxHLL(min_max_speed_naive),
diff --git a/test/test_tree_2d_shallowwater.jl b/test/test_tree_2d_shallowwater.jl
index 9a3ba36c7d..bcad663008 100644
--- a/test/test_tree_2d_shallowwater.jl
+++ b/test/test_tree_2d_shallowwater.jl
@@ -114,7 +114,7 @@ end
     end
 end
 
-@trixi_testset "elixir_shallowwater_well_balanced.jl with flux_nonconservative_ersing_etal" begin
+@trixi_testset "elixir_shallowwater_well_balanced.jl with flux_nonconservative_wintermeyer_etal" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_well_balanced.jl"),
                         l2=[
                             0.9130579602987146,
@@ -129,9 +129,7 @@ end
                             2.1130620376156584,
                         ],
                         surface_flux=(flux_wintermeyer_etal,
-                                      flux_nonconservative_ersing_etal),
-                        volume_flux=(flux_wintermeyer_etal,
-                                     flux_nonconservative_ersing_etal),
+                                      flux_nonconservative_wintermeyer_etal),
                         tspan=(0.0, 0.25))
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
@@ -172,15 +170,15 @@ end
     @test_trixi_include(joinpath(EXAMPLES_DIR,
                                  "elixir_shallowwater_source_terms_dirichlet.jl"),
                         l2=[
-                            0.0018746929418489125,
-                            0.017332321628469628,
-                            0.01634953679145536,
-                            6.274146767717023e-5,
+                            0.0018596727473552813,
+                            0.017306217777629147,
+                            0.016367646997420396,
+                            6.274146767723934e-5,
                         ],
                         linf=[
-                            0.016262353691956388,
-                            0.08726160620859424,
-                            0.09043621801418844,
+                            0.016548007102923368,
+                            0.08726160568822783,
+                            0.09043621622245013,
                             0.0001819675955490041,
                         ],
                         tspan=(0.0, 0.025))
@@ -248,7 +246,7 @@ end
     end
 end
 
-@trixi_testset "elixir_shallowwater_source_terms.jl with flux_nonconservative_ersing_etal" begin
+@trixi_testset "elixir_shallowwater_source_terms.jl with flux_nonconservative_wintermeyer_etal" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"),
                         l2=[
                             0.002471853426064005,
@@ -263,9 +261,7 @@ end
                             0.0001819675955490041,
                         ],
                         surface_flux=(flux_wintermeyer_etal,
-                                      flux_nonconservative_ersing_etal),
-                        volume_flux=(flux_wintermeyer_etal,
-                                     flux_nonconservative_ersing_etal),
+                                      flux_nonconservative_wintermeyer_etal),
                         tspan=(0.0, 0.25))
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
diff --git a/test/test_type.jl b/test/test_type.jl
index 9b8c3366cf..d13b626b06 100644
--- a/test/test_type.jl
+++ b/test/test_type.jl
@@ -1459,6 +1459,69 @@ isdir(outdir) && rm(outdir, recursive = true)
         end
     end
 
+    @timed_testset "Laplace Diffusion 1D" begin
+        for RealT in (Float32, Float64)
+            equations = @inferred LinearScalarAdvectionEquation1D(RealT(1))
+            equations_parabolic = @inferred LaplaceDiffusion1D(RealT(0.1), equations)
+
+            x = SVector(zero(RealT))
+            t = zero(RealT)
+            u = gradients = SVector(one(RealT))
+            orientation = 1
+
+            @test eltype(@inferred flux(u, gradients, orientation, equations_parabolic)) ==
+                  RealT
+
+            # TODO: BC tests for BoundaryConditionDirichlet
+            # TODO: BC tests for BoundaryConditionNeumann
+        end
+    end
+
+    @timed_testset "Laplace Diffusion 2D" begin
+        for RealT in (Float32, Float64)
+            equations = LinearScalarAdvectionEquation2D(RealT(1), RealT(1))
+            equations_parabolic = LaplaceDiffusion2D(RealT(0.1), equations)
+
+            x = SVector(zero(RealT), zero(RealT))
+            t = zero(RealT)
+            u = u_inner = u_outer = inv_h = gradients = SVector(one(RealT), one(RealT))
+            orientations = [1, 2]
+
+            for orientation in orientations
+                @test eltype(@inferred flux(u, gradients, orientation, equations_parabolic)) ==
+                      RealT
+            end
+
+            parabolic_solver = ViscousFormulationLocalDG(RealT(0.1))
+            @test eltype(@inferred Trixi.penalty(u_outer, u_inner, inv_h,
+                                                 equations_parabolic, parabolic_solver)) ==
+                  RealT
+        end
+    end
+
+    @timed_testset "Laplace Diffusion 3D" begin
+        for RealT in (Float32, Float64)
+            equations = LinearScalarAdvectionEquation3D(RealT(1), RealT(1), RealT(1))
+            equations_parabolic = LaplaceDiffusion3D(RealT(0.1), equations)
+
+            x = SVector(zero(RealT), zero(RealT), zero(RealT))
+            t = zero(RealT)
+            u = u_inner = u_outer = inv_h = gradients = SVector(one(RealT), one(RealT),
+                                                                one(RealT))
+            orientations = [1, 2, 3]
+
+            for orientation in orientations
+                @test eltype(@inferred flux(u, gradients, orientation, equations_parabolic)) ==
+                      RealT
+            end
+
+            parabolic_solver = ViscousFormulationLocalDG(RealT(0.1))
+            @test eltype(@inferred Trixi.penalty(u_outer, u_inner, inv_h,
+                                                 equations_parabolic, parabolic_solver)) ==
+                  RealT
+        end
+    end
+
     @timed_testset "Lattice Boltzmann 2D" begin
         for RealT in (Float32, Float64)
             equations = @inferred LatticeBoltzmannEquations2D(Ma = RealT(0.1), Re = 1000)
@@ -1926,6 +1989,80 @@ isdir(outdir) && rm(outdir, recursive = true)
         end
     end
 
+    @timed_testset "Polytropic Euler 2D" begin
+        for RealT in (Float32, Float64)
+            equations1 = @inferred PolytropicEulerEquations2D(RealT(1),
+                                                              RealT(1)) # equations.gamma == 1
+            equations2 = @inferred PolytropicEulerEquations2D(RealT(1.4), RealT(0.5))  # equations.gamma > 1
+
+            for equations in (equations1, equations2)
+                x = SVector(zero(RealT), zero(RealT))
+                t = zero(RealT)
+                u = u_ll = u_rr = SVector(one(RealT), one(RealT), one(RealT))
+                orientations = [1, 2]
+                directions = [1, 2, 3, 4]
+                normal_direction = SVector(one(RealT), zero(RealT))
+
+                @test eltype(@inferred initial_condition_convergence_test(x, t, equations)) ==
+                      RealT
+                @test eltype(@inferred source_terms_convergence_test(u, x, t, equations)) ==
+                      RealT
+                @test eltype(@inferred initial_condition_weak_blast_wave(x, t, equations)) ==
+                      RealT
+
+                @test eltype(@inferred flux(u, normal_direction, equations)) == RealT
+                if RealT == Float32
+                    # check `ln_mean` and `stolarsky_mean` (test broken)
+                    @test_broken eltype(@inferred flux_winters_etal(u_ll, u_rr,
+                                                                    normal_direction,
+                                                                    equations)) ==
+                                 RealT
+                else
+                    @test eltype(@inferred flux_winters_etal(u_ll, u_rr, normal_direction,
+                                                             equations)) ==
+                          RealT
+                end
+                @test eltype(@inferred min_max_speed_naive(u_ll, u_rr, normal_direction,
+                                                           equations)) ==
+                      RealT
+                @test eltype(@inferred min_max_speed_davis(u_ll, u_rr, normal_direction,
+                                                           equations)) ==
+                      RealT
+                @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, normal_direction,
+                                                           equations)) ==
+                      RealT
+
+                for orientation in orientations
+                    @test eltype(@inferred flux(u, orientation, equations)) == RealT
+                    if RealT == Float32
+                        # check `ln_mean` and `stolarsky_mean` (test broken)
+                        @test_broken eltype(@inferred flux_winters_etal(u_ll, u_rr,
+                                                                        orientation,
+                                                                        equations)) ==
+                                     RealT
+                    else
+                        @test eltype(@inferred flux_winters_etal(u_ll, u_rr, orientation,
+                                                                 equations)) ==
+                              RealT
+                    end
+                    @test eltype(@inferred min_max_speed_davis(u_ll, u_rr, orientation,
+                                                               equations)) ==
+                          RealT
+                    @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, orientation,
+                                                               equations)) ==
+                          RealT
+                end
+
+                @test eltype(@inferred Trixi.max_abs_speeds(u, equations)) == RealT
+                @test eltype(@inferred cons2prim(u, equations)) == RealT
+                @test eltype(@inferred prim2cons(u, equations)) == RealT
+                @test eltype(@inferred cons2entropy(u, equations)) == RealT
+                @test typeof(@inferred density(u, equations)) == RealT
+                @test typeof(@inferred pressure(u, equations)) == RealT
+            end
+        end
+    end
+
     @timed_testset "Shallow Water 1D" begin
         for RealT in (Float32, Float64)
             equations = @inferred ShallowWaterEquations1D(gravity_constant = RealT(9.81))
@@ -1971,8 +2108,6 @@ isdir(outdir) && rm(outdir, recursive = true)
             @test eltype(@inferred flux_nonconservative_audusse_etal(u_ll, u_rr,
                                                                      orientation,
                                                                      equations)) == RealT
-            @test eltype(@inferred flux_nonconservative_ersing_etal(u_ll, u_rr, orientation,
-                                                                    equations)) == RealT
             @test eltype(@inferred flux_fjordholm_etal(u_ll, u_rr, orientation,
                                                        equations)) == RealT
             @test eltype(@inferred flux_wintermeyer_etal(u_ll, u_rr, orientation,
@@ -2059,10 +2194,6 @@ isdir(outdir) && rm(outdir, recursive = true)
                                                                      normal_direction_ll,
                                                                      normal_direction_average,
                                                                      equations)) == RealT
-            @test eltype(@inferred flux_nonconservative_ersing_etal(u_ll, u_rr,
-                                                                    normal_direction_ll,
-                                                                    normal_direction_average,
-                                                                    equations)) == RealT
             @test eltype(@inferred flux_fjordholm_etal(u_ll, u_rr, normal_direction,
                                                        equations)) == RealT
             @test eltype(@inferred flux_wintermeyer_etal(u_ll, u_rr, normal_direction,
@@ -2107,10 +2238,6 @@ isdir(outdir) && rm(outdir, recursive = true)
                                                                          orientation,
                                                                          equations)) ==
                       RealT
-                @test eltype(eltype(@inferred flux_nonconservative_ersing_etal(u_ll, u_rr,
-                                                                               orientation,
-                                                                               equations))) ==
-                      RealT
                 @test eltype(@inferred flux_fjordholm_etal(u_ll, u_rr, orientation,
                                                            equations)) == RealT
                 @test eltype(@inferred flux_wintermeyer_etal(u_ll, u_rr, orientation,
diff --git a/test/test_unstructured_2d.jl b/test/test_unstructured_2d.jl
index 6682fc2375..4a2e727c1d 100644
--- a/test/test_unstructured_2d.jl
+++ b/test/test_unstructured_2d.jl
@@ -378,16 +378,16 @@ end
 @trixi_testset "elixir_shallowwater_well_balanced.jl with FluxHydrostaticReconstruction" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_well_balanced.jl"),
                         l2=[
-                            1.2164292510839085,
-                            1.2643106818778908e-12,
-                            1.269230436589819e-12,
-                            1.2164292510839079,
+                            1.2164292510839063,
+                            1.2676379081600215e-12,
+                            1.255855785593831e-12,
+                            1.2164292510839074,
                         ],
                         linf=[
-                            1.513851228231562,
-                            1.6670644673575802e-11,
-                            1.8426585188623954e-11,
-                            1.513851228231574,
+                            1.5138512282315604,
+                            1.658245722058109e-11,
+                            1.8665562182185795e-11,
+                            1.5138512282315737,
                         ],
                         surface_flux=(FluxHydrostaticReconstruction(flux_lax_friedrichs,
                                                                     hydrostatic_reconstruction_audusse_etal),
@@ -403,7 +403,7 @@ end
     end
 end
 
-@trixi_testset "elixir_shallowwater_well_balanced.jl with flux_nonconservative_ersing_etal" begin
+@trixi_testset "elixir_shallowwater_well_balanced.jl with flux_nonconservative_wintermeyer_etal" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_well_balanced.jl"),
                         l2=[
                             1.2164292510839083,
@@ -418,9 +418,7 @@ end
                             1.513851228231574,
                         ],
                         surface_flux=(flux_wintermeyer_etal,
-                                      flux_nonconservative_ersing_etal),
-                        volume_flux=(flux_wintermeyer_etal,
-                                     flux_nonconservative_ersing_etal),
+                                      flux_nonconservative_wintermeyer_etal),
                         tspan=(0.0, 0.25))
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
@@ -485,7 +483,7 @@ end
     end
 end
 
-@trixi_testset "elixir_shallowwater_source_terms.jl with flux_nonconservative_ersing_etal" begin
+@trixi_testset "elixir_shallowwater_source_terms.jl with flux_nonconservative_wintermeyer_etal" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"),
                         l2=[
                             0.001118046975499805,
@@ -500,9 +498,7 @@ end
                             2.6407324614341476e-5,
                         ],
                         surface_flux=(flux_wintermeyer_etal,
-                                      flux_nonconservative_ersing_etal),
-                        volume_flux=(flux_wintermeyer_etal,
-                                     flux_nonconservative_ersing_etal),
+                                      flux_nonconservative_wintermeyer_etal),
                         tspan=(0.0, 0.025))
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
@@ -544,12 +540,16 @@ end
 @trixi_testset "elixir_shallowwater_dirichlet.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_dirichlet.jl"),
                         l2=[
-                            1.1577518608938916e-5, 4.859252379740366e-13,
-                            4.639600837197925e-13, 1.1577518608952174e-5,
+                            1.1577518608950964e-5,
+                            4.761947272222427e-13,
+                            4.546045873135486e-13,
+                            1.157751860893347e-5,
                         ],
                         linf=[
-                            8.3940638787805e-5, 1.1446362498574484e-10,
-                            1.1124515748367981e-10, 8.39406387962427e-5,
+                            8.394063879002545e-5,
+                            1.1211566736150389e-10,
+                            1.0890426250906834e-10,
+                            8.394063879602065e-5,
                         ],
                         tspan=(0.0, 2.0))
     # Ensure that we do not have excessive memory allocations