diff --git a/NEWS.md b/NEWS.md
index bca7d03e954..9371cafa07f 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/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 00000000000..68680673712
--- /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 00000000000..b5b56220004
--- /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 9138586cccf..d8893854811 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 00000000000..a3366caa317
--- /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 00000000000..106b4821144
--- /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/src/callbacks_step/amr_dg2d.jl b/src/callbacks_step/amr_dg2d.jl
index 94524b23a3a..062020b6d42 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 3f67951bafe..594c30dcca5 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/meshes/structured_mesh.jl b/src/meshes/structured_mesh.jl
index 553aabbbc20..9e5c55197d8 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 22d14e38c5c..6a7c5c7834b 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 6d66bc68a26..29de4efc6d0 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 cca9093ec51..4f9465b85dc 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 7c7fc03898c..75e65c8c380 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 a15690a7629..2e837f79ad8 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 ef5ac2c9de3..f56dc9cd769 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 7483cde2752..7b69d1c0cf2 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))
@@ -542,16 +542,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))
@@ -586,6 +586,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 b63d2a105ac..c1fcc355218 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 5c452444be0..940d2c43372 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