From 642da1af9f9cc390f1d3d2a47a5fd07628a632f0 Mon Sep 17 00:00:00 2001
From: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com>
Date: Mon, 19 Jun 2023 09:07:15 +0200
Subject: [PATCH 1/3] Use `Pointerwrapper`s in `P4estMesh` (#1434)

* use PointerWrappers

* fix typo

* fix typo

* fixes

* fixes

* fixes

* unsafe_load_sc returns PointerWrapper

* bug fixes

* bug fixes

* bug fixes

* add comment about

* rename unsafe_load_* to load_pointerwrapper_*

* format

* fix merge conflicts

* fix bad format again

* fix

* add unsafe_wrap_sc again

* fix

* fix

* Update src/auxiliary/p4est.jl

Co-authored-by: Michael Schlottke-Lakemper <michael@sloede.com>

* introduce generic type PointerOrWrapper{T}

* format

---------

Co-authored-by: Michael Schlottke-Lakemper <michael@sloede.com>
Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com>
---
 src/auxiliary/p4est.jl                        |  94 +++++++-----
 src/callbacks_step/amr.jl                     |  30 ++--
 src/callbacks_step/amr_dg.jl                  |   3 +-
 src/callbacks_step/analysis.jl                |   2 +-
 src/meshes/p4est_mesh.jl                      | 117 +++++++-------
 src/solvers/dgsem_p4est/containers.jl         | 122 +++++++--------
 src/solvers/dgsem_p4est/containers_2d.jl      |   2 +-
 src/solvers/dgsem_p4est/containers_3d.jl      |   2 +-
 .../dgsem_p4est/containers_parallel.jl        | 145 +++++++++---------
 src/solvers/dgsem_p4est/dg_parallel.jl        |  80 +++++-----
 10 files changed, 311 insertions(+), 286 deletions(-)

diff --git a/src/auxiliary/p4est.jl b/src/auxiliary/p4est.jl
index 93b5166cd81..968af339cbd 100644
--- a/src/auxiliary/p4est.jl
+++ b/src/auxiliary/p4est.jl
@@ -24,35 +24,38 @@ function init_p4est()
     return nothing
 end
 
+# for convenience to either pass a Ptr or a PointerWrapper
+const PointerOrWrapper = Union{Ptr{T}, PointerWrapper{T}} where {T}
+
 # Convert sc_array of type T to Julia array
-function unsafe_wrap_sc(::Type{T}, sc_array::Ptr{sc_array}) where {T}
-    sc_array_obj = unsafe_load(sc_array)
+function unsafe_wrap_sc(::Type{T}, sc_array_ptr::Ptr{sc_array}) where {T}
+    sc_array_obj = unsafe_load(sc_array_ptr)
     return unsafe_wrap_sc(T, sc_array_obj)
 end
 
 function unsafe_wrap_sc(::Type{T}, sc_array_obj::sc_array) where {T}
     elem_count = sc_array_obj.elem_count
     array = sc_array_obj.array
-
     return unsafe_wrap(Array, Ptr{T}(array), elem_count)
 end
 
-# Load the ith element (1-indexed) of an sc array of type T
-function unsafe_load_sc(::Type{T}, sc_array::Ptr{sc_array}, i = 1) where {T}
-    sc_array_obj = unsafe_load(sc_array)
-    return unsafe_load_sc(T, sc_array_obj, i)
-end
+function unsafe_wrap_sc(::Type{T}, sc_array_pw::PointerWrapper{sc_array}) where {T}
+    elem_count = sc_array_pw.elem_count[]
+    array = sc_array_pw.array
 
-function unsafe_load_sc(::Type{T}, sc_array_obj::sc_array, i = 1) where {T}
-    element_size = sc_array_obj.elem_size
-    @assert element_size == sizeof(T)
+    return unsafe_wrap(Array, Ptr{T}(pointer(array)), elem_count)
+end
 
-    return unsafe_load(Ptr{T}(sc_array_obj.array), i)
+# Load the ith element (1-indexed) of an sc array of type T as PointerWrapper
+function load_pointerwrapper_sc(::Type{T}, sc_array::PointerWrapper{sc_array},
+                                i::Integer = 1) where {T}
+    return PointerWrapper(T, pointer(sc_array.array) + (i - 1) * sizeof(T))
 end
 
 # Create new `p4est` from a p4est_connectivity
 # 2D
-function new_p4est(connectivity::Ptr{p4est_connectivity_t}, initial_refinement_level)
+function new_p4est(connectivity::PointerOrWrapper{p4est_connectivity_t},
+                   initial_refinement_level)
     comm = P4est.uses_mpi() ? mpi_comm() : 0 # Use Trixi.jl's MPI communicator if p4est supports MPI
     p4est_new_ext(comm,
                   connectivity,
@@ -65,7 +68,8 @@ function new_p4est(connectivity::Ptr{p4est_connectivity_t}, initial_refinement_l
 end
 
 # 3D
-function new_p4est(connectivity::Ptr{p8est_connectivity_t}, initial_refinement_level)
+function new_p4est(connectivity::PointerOrWrapper{p8est_connectivity_t},
+                   initial_refinement_level)
     comm = P4est.uses_mpi() ? mpi_comm() : 0 # Use Trixi.jl's MPI communicator if p4est supports MPI
     p8est_new_ext(comm, connectivity, 0, initial_refinement_level, true,
                   2 * sizeof(Int), C_NULL, C_NULL)
@@ -73,13 +77,13 @@ end
 
 # Save `p4est` data to file
 # 2D
-function save_p4est!(file, p4est::Ptr{p4est_t})
+function save_p4est!(file, p4est::PointerOrWrapper{p4est_t})
     # Don't save user data of the quads
     p4est_save(file, p4est, false)
 end
 
 # 3D
-function save_p4est!(file, p8est::Ptr{p8est_t})
+function save_p4est!(file, p8est::PointerOrWrapper{p8est_t})
     # Don't save user data of the quads
     p8est_save(file, p8est, false)
 end
@@ -107,27 +111,33 @@ read_inp_p4est(meshfile, ::Val{3}) = p8est_connectivity_read_inp(meshfile)
 
 # Refine `p4est` if refine_fn_c returns 1
 # 2D
-function refine_p4est!(p4est::Ptr{p4est_t}, recursive, refine_fn_c, init_fn_c)
+function refine_p4est!(p4est::PointerOrWrapper{p4est_t}, recursive, refine_fn_c,
+                       init_fn_c)
     p4est_refine(p4est, recursive, refine_fn_c, init_fn_c)
 end
 # 3D
-function refine_p4est!(p8est::Ptr{p8est_t}, recursive, refine_fn_c, init_fn_c)
+function refine_p4est!(p8est::PointerOrWrapper{p8est_t}, recursive, refine_fn_c,
+                       init_fn_c)
     p8est_refine(p8est, recursive, refine_fn_c, init_fn_c)
 end
 
 # Refine `p4est` if coarsen_fn_c returns 1
 # 2D
-function coarsen_p4est!(p4est::Ptr{p4est_t}, recursive, coarsen_fn_c, init_fn_c)
+function coarsen_p4est!(p4est::PointerOrWrapper{p4est_t}, recursive, coarsen_fn_c,
+                        init_fn_c)
     p4est_coarsen(p4est, recursive, coarsen_fn_c, init_fn_c)
 end
 # 3D
-function coarsen_p4est!(p8est::Ptr{p8est_t}, recursive, coarsen_fn_c, init_fn_c)
+function coarsen_p4est!(p8est::PointerOrWrapper{p8est_t}, recursive, coarsen_fn_c,
+                        init_fn_c)
     p8est_coarsen(p8est, recursive, coarsen_fn_c, init_fn_c)
 end
 
 # Create new ghost layer from p4est, only connections via faces are relevant
 # 2D
-ghost_new_p4est(p4est::Ptr{p4est_t}) = p4est_ghost_new(p4est, P4est.P4EST_CONNECT_FACE)
+function ghost_new_p4est(p4est::PointerOrWrapper{p4est_t})
+    p4est_ghost_new(p4est, P4est.P4EST_CONNECT_FACE)
+end
 # 3D
 # In 3D it is not sufficient to use `P8EST_CONNECT_FACE`. Consider the neighbor elements of a mortar
 # in 3D. We have to determine which MPI ranks are involved in this mortar.
@@ -147,28 +157,37 @@ ghost_new_p4est(p4est::Ptr{p4est_t}) = p4est_ghost_new(p4est, P4est.P4EST_CONNEC
 # `P8EST_CONNECT_FACE`. But if it is not in the ghost layer, it will not be available in
 # `iterate_p4est` and thus we cannot determine its MPI rank
 # (see https://github.com/cburstedde/p4est/blob/439bc9aae849555256ddfe4b03d1f9fe8d18ff0e/src/p8est_iterate.h#L66-L72).
-ghost_new_p4est(p8est::Ptr{p8est_t}) = p8est_ghost_new(p8est, P4est.P8EST_CONNECT_FULL)
+function ghost_new_p4est(p8est::PointerOrWrapper{p8est_t})
+    p8est_ghost_new(p8est, P4est.P8EST_CONNECT_FULL)
+end
 
 # Check if ghost layer is valid
 # 2D
-function ghost_is_valid_p4est(p4est::Ptr{p4est_t}, ghost_layer::Ptr{p4est_ghost_t})
+function ghost_is_valid_p4est(p4est::PointerOrWrapper{p4est_t},
+                              ghost_layer::Ptr{p4est_ghost_t})
     return p4est_ghost_is_valid(p4est, ghost_layer)
 end
 # 3D
-function ghost_is_valid_p4est(p4est::Ptr{p8est_t}, ghost_layer::Ptr{p8est_ghost_t})
+function ghost_is_valid_p4est(p4est::PointerOrWrapper{p8est_t},
+                              ghost_layer::Ptr{p8est_ghost_t})
     return p8est_ghost_is_valid(p4est, ghost_layer)
 end
 
 # Destroy ghost layer
 # 2D
-ghost_destroy_p4est(ghost_layer::Ptr{p4est_ghost_t}) = p4est_ghost_destroy(ghost_layer)
+function ghost_destroy_p4est(ghost_layer::PointerOrWrapper{p4est_ghost_t})
+    p4est_ghost_destroy(ghost_layer)
+end
 # 3D
-ghost_destroy_p4est(ghost_layer::Ptr{p8est_ghost_t}) = p8est_ghost_destroy(ghost_layer)
+function ghost_destroy_p4est(ghost_layer::PointerOrWrapper{p8est_ghost_t})
+    p8est_ghost_destroy(ghost_layer)
+end
 
 # Let `p4est` iterate over each cell volume and cell face.
 # Call iter_volume_c for each cell and iter_face_c for each face.
 # 2D
-function iterate_p4est(p4est::Ptr{p4est_t}, user_data; ghost_layer = C_NULL,
+function iterate_p4est(p4est::PointerOrWrapper{p4est_t}, user_data;
+                       ghost_layer = C_NULL,
                        iter_volume_c = C_NULL, iter_face_c = C_NULL)
     if user_data === C_NULL
         user_data_ptr = user_data
@@ -191,7 +210,8 @@ function iterate_p4est(p4est::Ptr{p4est_t}, user_data; ghost_layer = C_NULL,
 end
 
 # 3D
-function iterate_p4est(p8est::Ptr{p8est_t}, user_data; ghost_layer = C_NULL,
+function iterate_p4est(p8est::PointerOrWrapper{p8est_t}, user_data;
+                       ghost_layer = C_NULL,
                        iter_volume_c = C_NULL, iter_face_c = C_NULL)
     if user_data === C_NULL
         user_data_ptr = user_data
@@ -216,23 +236,25 @@ end
 
 # Load i-th element of the sc_array info.sides of the type p[48]est_iter_face_side_t
 # 2D version
-function unsafe_load_side(info::Ptr{p4est_iter_face_info_t}, i = 1)
-    return unsafe_load_sc(p4est_iter_face_side_t, unsafe_load(info).sides, i)
+function load_pointerwrapper_side(info::PointerWrapper{p4est_iter_face_info_t},
+                                  i::Integer = 1)
+    return load_pointerwrapper_sc(p4est_iter_face_side_t, info.sides, i)
 end
 
 # 3D version
-function unsafe_load_side(info::Ptr{p8est_iter_face_info_t}, i = 1)
-    return unsafe_load_sc(p8est_iter_face_side_t, unsafe_load(info).sides, i)
+function load_pointerwrapper_side(info::PointerWrapper{p8est_iter_face_info_t},
+                                  i::Integer = 1)
+    return load_pointerwrapper_sc(p8est_iter_face_side_t, info.sides, i)
 end
 
 # Load i-th element of the sc_array p4est.trees of the type p[48]est_tree_t
 # 2D version
-function unsafe_load_tree(p4est::Ptr{p4est_t}, i = 1)
-    return unsafe_load_sc(p4est_tree_t, unsafe_load(p4est).trees, i)
+function load_pointerwrapper_tree(p4est::PointerWrapper{p4est_t}, i::Integer = 1)
+    return load_pointerwrapper_sc(p4est_tree_t, p4est.trees, i)
 end
 
 # 3D version
-function unsafe_load_tree(p8est::Ptr{p8est_t}, i = 1)
-    return unsafe_load_sc(p8est_tree_t, unsafe_load(p8est).trees, i)
+function load_pointerwrapper_tree(p8est::PointerWrapper{p8est_t}, i::Integer = 1)
+    return load_pointerwrapper_sc(p8est_tree_t, p8est.trees, i)
 end
 end # @muladd
diff --git a/src/callbacks_step/amr.jl b/src/callbacks_step/amr.jl
index d6e19b79886..bef49b4c482 100644
--- a/src/callbacks_step/amr.jl
+++ b/src/callbacks_step/amr.jl
@@ -348,24 +348,24 @@ end
 
 # Copy controller values to quad user data storage, will be called below
 function copy_to_quad_iter_volume(info, user_data)
-    info_obj = unsafe_load(info)
+    info_pw = PointerWrapper(info)
 
     # Load tree from global trees array, one-based indexing
-    tree = unsafe_load_tree(info_obj.p4est, info_obj.treeid + 1)
+    tree_pw = load_pointerwrapper_tree(info_pw.p4est, info_pw.treeid[] + 1)
     # Quadrant numbering offset of this quadrant
-    offset = tree.quadrants_offset
+    offset = tree_pw.quadrants_offset[]
     # Global quad ID
-    quad_id = offset + info_obj.quadid
+    quad_id = offset + info_pw.quadid[]
 
     # Access user_data = lambda
-    user_data_ptr = Ptr{Int}(user_data)
+    user_data_pw = PointerWrapper(Int, user_data)
     # Load controller_value = lambda[quad_id + 1]
-    controller_value = unsafe_load(user_data_ptr, quad_id + 1)
+    controller_value = user_data_pw[quad_id + 1]
 
     # Access quadrant's user data ([global quad ID, controller_value])
-    quad_data_ptr = Ptr{Int}(unsafe_load(info_obj.quad.p.user_data))
+    quad_data_pw = PointerWrapper(Int, info_pw.quad.p.user_data[])
     # Save controller value to quadrant's user data.
-    unsafe_store!(quad_data_ptr, controller_value, 2)
+    quad_data_pw[2] = controller_value
 
     return nothing
 end
@@ -599,22 +599,22 @@ function current_element_levels(mesh::TreeMesh, solver, cache)
 end
 
 function extract_levels_iter_volume(info, user_data)
-    info_obj = unsafe_load(info)
+    info_pw = PointerWrapper(info)
 
     # Load tree from global trees array, one-based indexing
-    tree = unsafe_load_tree(info_obj.p4est, info_obj.treeid + 1)
+    tree_pw = load_pointerwrapper_tree(info_pw.p4est, info_pw.treeid[] + 1)
     # Quadrant numbering offset of this quadrant
-    offset = tree.quadrants_offset
+    offset = tree_pw.quadrants_offset[]
     # Global quad ID
-    quad_id = offset + info_obj.quadid
+    quad_id = offset + info_pw.quadid[]
     # Julia element ID
     element_id = quad_id + 1
 
-    current_level = unsafe_load(info_obj.quad.level)
+    current_level = info_pw.quad.level[]
 
     # Unpack user_data = current_levels and save current element level
-    ptr = Ptr{Int}(user_data)
-    unsafe_store!(ptr, current_level, element_id)
+    pw = PointerWrapper(Int, user_data)
+    pw[element_id] = current_level
 
     return nothing
 end
diff --git a/src/callbacks_step/amr_dg.jl b/src/callbacks_step/amr_dg.jl
index 19bbebd9254..1dcfdccdea8 100644
--- a/src/callbacks_step/amr_dg.jl
+++ b/src/callbacks_step/amr_dg.jl
@@ -9,8 +9,7 @@
 function rebalance_solver!(u_ode::AbstractVector, mesh::ParallelP4estMesh, equations,
                            dg::DGSEM, cache, old_global_first_quadrant)
     # mpi ranks are 0-based, this array uses 1-based indices
-    global_first_quadrant = unsafe_wrap(Array,
-                                        unsafe_load(mesh.p4est).global_first_quadrant,
+    global_first_quadrant = unsafe_wrap(Array, mesh.p4est.global_first_quadrant,
                                         mpi_nranks() + 1)
     if global_first_quadrant[mpi_rank() + 1] ==
        old_global_first_quadrant[mpi_rank() + 1] &&
diff --git a/src/callbacks_step/analysis.jl b/src/callbacks_step/analysis.jl
index 7fa2e21a244..8cf43a1d15e 100644
--- a/src/callbacks_step/analysis.jl
+++ b/src/callbacks_step/analysis.jl
@@ -508,7 +508,7 @@ function print_amr_information(callbacks, mesh::P4estMesh, solver, cache)
 
     elements_per_level = zeros(P4EST_MAXLEVEL + 1)
 
-    for tree in unsafe_wrap_sc(p4est_tree_t, unsafe_load(mesh.p4est).trees)
+    for tree in unsafe_wrap_sc(p4est_tree_t, mesh.p4est.trees)
         elements_per_level .+= tree.quadrants_per_level
     end
 
diff --git a/src/meshes/p4est_mesh.jl b/src/meshes/p4est_mesh.jl
index ddd6cf473e4..60db285e04f 100644
--- a/src/meshes/p4est_mesh.jl
+++ b/src/meshes/p4est_mesh.jl
@@ -13,9 +13,9 @@ to manage trees and mesh refinement.
 """
 mutable struct P4estMesh{NDIMS, RealT <: Real, IsParallel, P, Ghost, NDIMSP2, NNODES} <:
                AbstractMesh{NDIMS}
-    p4est::P # Either Ptr{p4est_t} or Ptr{p8est_t}
-    is_parallel::IsParallel
-    ghost::Ghost # Either Ptr{p4est_ghost_t} or Ptr{p8est_ghost_t}
+    p4est       :: P # Either PointerWrapper{p4est_t} or PointerWrapper{p8est_t}
+    is_parallel :: IsParallel
+    ghost       :: Ghost # Either PointerWrapper{p4est_ghost_t} or PointerWrapper{p8est_ghost_t}
     # Coordinates at the nodes specified by the tensor product of `nodes` (NDIMS times).
     # This specifies the geometry interpolation for each tree.
     tree_node_coordinates::Array{RealT, NDIMSP2} # [dimension, i, j, k, tree]
@@ -43,18 +43,21 @@ mutable struct P4estMesh{NDIMS, RealT <: Real, IsParallel, P, Ghost, NDIMSP2, NN
             is_parallel = False()
         end
 
+        p4est_pw = PointerWrapper(p4est)
+
         ghost = ghost_new_p4est(p4est)
+        ghost_pw = PointerWrapper(ghost)
 
         mesh = new{NDIMS, eltype(tree_node_coordinates), typeof(is_parallel),
-                   typeof(p4est), typeof(ghost), NDIMS + 2, length(nodes)}(p4est,
-                                                                           is_parallel,
-                                                                           ghost,
-                                                                           tree_node_coordinates,
-                                                                           nodes,
-                                                                           boundary_names,
-                                                                           current_filename,
-                                                                           unsaved_changes,
-                                                                           p4est_partition_allow_for_coarsening)
+                   typeof(p4est_pw), typeof(ghost_pw), NDIMS + 2, length(nodes)}(p4est_pw,
+                                                                                 is_parallel,
+                                                                                 ghost_pw,
+                                                                                 tree_node_coordinates,
+                                                                                 nodes,
+                                                                                 boundary_names,
+                                                                                 current_filename,
+                                                                                 unsaved_changes,
+                                                                                 p4est_partition_allow_for_coarsening)
 
         # Destroy `p4est` structs when the mesh is garbage collected
         finalizer(destroy_mesh, mesh)
@@ -70,14 +73,14 @@ const ParallelP4estMesh{NDIMS} = P4estMesh{NDIMS, <:Real, <:True}
 @inline mpi_parallel(mesh::ParallelP4estMesh) = True()
 
 function destroy_mesh(mesh::P4estMesh{2})
-    connectivity = unsafe_load(mesh.p4est).connectivity
+    connectivity = mesh.p4est.connectivity
     p4est_ghost_destroy(mesh.ghost)
     p4est_destroy(mesh.p4est)
     p4est_connectivity_destroy(connectivity)
 end
 
 function destroy_mesh(mesh::P4estMesh{3})
-    connectivity = unsafe_load(mesh.p4est).connectivity
+    connectivity = mesh.p4est.connectivity
     p8est_ghost_destroy(mesh.ghost)
     p8est_destroy(mesh.p4est)
     p8est_connectivity_destroy(connectivity)
@@ -87,11 +90,10 @@ end
 @inline Base.real(::P4estMesh{NDIMS, RealT}) where {NDIMS, RealT} = RealT
 
 @inline function ntrees(mesh::P4estMesh)
-    trees = unsafe_load(mesh.p4est).trees
-    return unsafe_load(trees).elem_count
+    return mesh.p4est.trees.elem_count[]
 end
 # returns Int32 by default which causes a weird method error when creating the cache
-@inline ncells(mesh::P4estMesh) = Int(unsafe_load(mesh.p4est).local_num_quadrants)
+@inline ncells(mesh::P4estMesh) = Int(mesh.p4est.local_num_quadrants[])
 
 function Base.show(io::IO, mesh::P4estMesh)
     print(io, "P4estMesh{", ndims(mesh), ", ", real(mesh), "}")
@@ -387,14 +389,14 @@ function p4est_mesh_from_hohqmesh_abaqus(meshfile, initial_refinement_level,
                                          n_dimensions, RealT)
     # Create the mesh connectivity using `p4est`
     connectivity = read_inp_p4est(meshfile, Val(n_dimensions))
-    connectivity_obj = unsafe_load(connectivity)
+    connectivity_pw = PointerWrapper(connectivity)
 
     # These need to be of the type Int for unsafe_wrap below to work
-    n_trees::Int = connectivity_obj.num_trees
-    n_vertices::Int = connectivity_obj.num_vertices
+    n_trees::Int = connectivity_pw.num_trees[]
+    n_vertices::Int = connectivity_pw.num_vertices[]
 
     # Extract a copy of the element vertices to compute the tree node coordinates
-    vertices = unsafe_wrap(Array, connectivity_obj.vertices, (3, n_vertices))
+    vertices = unsafe_wrap(Array, connectivity_pw.vertices, (3, n_vertices))
 
     # Readin all the information from the mesh file into a string array
     file_lines = readlines(open(meshfile))
@@ -445,14 +447,14 @@ function p4est_mesh_from_standard_abaqus(meshfile, mapping, polydeg,
                                          initial_refinement_level, n_dimensions, RealT)
     # Create the mesh connectivity using `p4est`
     connectivity = read_inp_p4est(meshfile, Val(n_dimensions))
-    connectivity_obj = unsafe_load(connectivity)
+    connectivity_pw = PointerWrapper(connectivity)
 
     # These need to be of the type Int for unsafe_wrap below to work
-    n_trees::Int = connectivity_obj.num_trees
-    n_vertices::Int = connectivity_obj.num_vertices
+    n_trees::Int = connectivity_pw.num_trees[]
+    n_vertices::Int = connectivity_pw.num_vertices[]
 
-    vertices = unsafe_wrap(Array, connectivity_obj.vertices, (3, n_vertices))
-    tree_to_vertex = unsafe_wrap(Array, connectivity_obj.tree_to_vertex,
+    vertices = unsafe_wrap(Array, connectivity_pw.vertices, (3, n_vertices))
+    tree_to_vertex = unsafe_wrap(Array, connectivity_pw.tree_to_vertex,
                                  (2^n_dimensions, n_trees))
 
     basis = LobattoLegendreBasis(RealT, polydeg)
@@ -1511,17 +1513,18 @@ end
 
 function update_ghost_layer!(mesh::P4estMesh)
     ghost_destroy_p4est(mesh.ghost)
-    mesh.ghost = ghost_new_p4est(mesh.p4est)
+    mesh.ghost = PointerWrapper(ghost_new_p4est(mesh.p4est))
 end
 
 function init_fn(p4est, which_tree, quadrant)
     # Unpack quadrant's user data ([global quad ID, controller_value])
-    ptr = Ptr{Int}(unsafe_load(quadrant.p.user_data))
+    # Use `unsafe_load` here since `quadrant.p.user_data isa Ptr{Ptr{Nothing}}`
+    # and we only need the first (only!) entry
+    pw = PointerWrapper(Int, unsafe_load(quadrant.p.user_data))
 
     # Initialize quad ID as -1 and controller_value as 0 (don't refine or coarsen)
-    unsafe_store!(ptr, -1, 1)
-    unsafe_store!(ptr, 0, 2)
-
+    pw[1] = -1
+    pw[2] = 0
     return nothing
 end
 
@@ -1539,8 +1542,10 @@ end
 function refine_fn(p4est, which_tree, quadrant)
     # Controller value has been copied to the quadrant's user data storage before.
     # Unpack quadrant's user data ([global quad ID, controller_value]).
-    ptr = Ptr{Int}(unsafe_load(quadrant.p.user_data))
-    controller_value = unsafe_load(ptr, 2)
+    # Use `unsafe_load` here since `quadrant.p.user_data isa Ptr{Ptr{Nothing}}`
+    # and we only need the first (only!) entry
+    pw = PointerWrapper(Int, unsafe_load(quadrant.p.user_data))
+    controller_value = pw[2]
 
     if controller_value > 0
         # return true (refine)
@@ -1586,9 +1591,9 @@ function coarsen_fn(p4est, which_tree, quadrants_ptr)
 
     # Controller value has been copied to the quadrant's user data storage before.
     # Load controller value from quadrant's user data ([global quad ID, controller_value]).
-    function controller_value(i)
-        unsafe_load(Ptr{Int}(unsafe_load(quadrants[i].p.user_data)), 2)
-    end
+    # Use `unsafe_load` here since `quadrant.p.user_data isa Ptr{Ptr{Nothing}}`
+    # and we only need the first (only!) entry
+    controller_value(i) = PointerWrapper(Int, unsafe_load(quadrants[i].p.user_data))[2]
 
     # `p4est` calls this function for each 2^ndims quads that could be coarsened to a single one.
     # Only coarsen if all these 2^ndims quads have been marked for coarsening.
@@ -1671,20 +1676,19 @@ end
 
 # Copy global quad ID to quad's user data storage, will be called below
 function save_original_id_iter_volume(info, user_data)
-    info_obj = unsafe_load(info)
+    info_pw = PointerWrapper(info)
 
     # Load tree from global trees array, one-based indexing
-    tree = unsafe_load_tree(info_obj.p4est, info_obj.treeid + 1)
+    tree_pw = load_pointerwrapper_tree(info_pw.p4est, info_pw.treeid[] + 1)
     # Quadrant numbering offset of this quadrant
-    offset = tree.quadrants_offset
+    offset = tree_pw.quadrants_offset[]
     # Global quad ID
-    quad_id = offset + info_obj.quadid
+    quad_id = offset + info_pw.quadid[]
 
     # Unpack quadrant's user data ([global quad ID, controller_value])
-    ptr = Ptr{Int}(unsafe_load(info_obj.quad.p.user_data))
+    pw = PointerWrapper(Int, info_pw.quad.p.user_data[])
     # Save global quad ID
-    unsafe_store!(ptr, quad_id, 1)
-
+    pw[1] = quad_id
     return nothing
 end
 
@@ -1708,24 +1712,23 @@ end
 
 # Extract information about which cells have been changed
 function collect_changed_iter_volume(info, user_data)
-    info_obj = unsafe_load(info)
+    info_pw = PointerWrapper(info)
 
     # The original element ID has been saved to user_data before.
     # Load original quad ID from quad's user data ([global quad ID, controller_value]).
-    quad_data_ptr = Ptr{Int}(unsafe_load(info_obj.quad.p.user_data))
-    original_id = unsafe_load(quad_data_ptr, 1)
+    quad_data_pw = PointerWrapper(Int, info_pw.quad.p.user_data[])
+    original_id = quad_data_pw[1]
 
     # original_id of cells that have been newly created is -1
     if original_id >= 0
         # Unpack user_data = original_cells
-        user_data_ptr = Ptr{Int}(user_data)
+        user_data_pw = PointerWrapper(Int, user_data)
 
         # If quad has an original_id, it existed before refinement/coarsening,
         # and therefore wasn't changed.
         # Mark original_id as "not changed during refinement/coarsening" in original_cells
-        unsafe_store!(user_data_ptr, 0, original_id + 1)
+        user_data_pw[original_id + 1] = 0
     end
-
     return nothing
 end
 
@@ -1756,29 +1759,27 @@ end
 
 # Extract newly created cells
 function collect_new_iter_volume(info, user_data)
-    info_obj = unsafe_load(info)
+    info_pw = PointerWrapper(info)
 
     # The original element ID has been saved to user_data before.
     # Unpack quadrant's user data ([global quad ID, controller_value]).
-    quad_data_ptr = Ptr{Int}(unsafe_load(info_obj.quad.p.user_data))
-    original_id = unsafe_load(quad_data_ptr, 1)
+    original_id = PointerWrapper(Int, info_pw.quad.p.user_data[])[1]
 
     # original_id of cells that have been newly created is -1
     if original_id < 0
         # Load tree from global trees array, one-based indexing
-        tree = unsafe_load_tree(info_obj.p4est, info_obj.treeid + 1)
+        tree_pw = load_pointerwrapper_tree(info_pw.p4est, info_pw.treeid[] + 1)
         # Quadrant numbering offset of this quadrant
-        offset = tree.quadrants_offset
+        offset = tree_pw.quadrants_offset[]
         # Global quad ID
-        quad_id = offset + info_obj.quadid
+        quad_id = offset + info_pw.quadid[]
 
         # Unpack user_data = original_cells
-        user_data_ptr = Ptr{Int}(user_data)
+        user_data_pw = PointerWrapper(Int, user_data)
 
         # Mark cell as "newly created during refinement/coarsening/balancing"
-        unsafe_store!(user_data_ptr, 1, quad_id + 1)
+        user_data_pw[quad_id + 1] = 1
     end
-
     return nothing
 end
 
diff --git a/src/solvers/dgsem_p4est/containers.jl b/src/solvers/dgsem_p4est/containers.jl
index 9b87de777a6..2b9c6987d24 100644
--- a/src/solvers/dgsem_p4est/containers.jl
+++ b/src/solvers/dgsem_p4est/containers.jl
@@ -276,18 +276,18 @@ function init_boundaries!(boundaries, mesh::P4estMesh)
 end
 
 # Function barrier for type stability
-function init_boundaries_iter_face_inner(info, boundaries, boundary_id, mesh)
+function init_boundaries_iter_face_inner(info_pw, boundaries, boundary_id, mesh)
     # Extract boundary data
-    side = unsafe_load_side(info)
+    side_pw = load_pointerwrapper_side(info_pw)
     # Get local tree, one-based indexing
-    tree = unsafe_load_tree(mesh.p4est, side.treeid + 1)
+    tree_pw = load_pointerwrapper_tree(mesh.p4est, side_pw.treeid[] + 1)
     # Quadrant numbering offset of this quadrant
-    offset = tree.quadrants_offset
+    offset = tree_pw.quadrants_offset[]
 
     # Verify before accessing is.full, but this should never happen
-    @assert side.is_hanging == false
+    @assert side_pw.is_hanging[] == false
 
-    local_quad_id = side.is.full.quadid
+    local_quad_id = side_pw.is.full.quadid[]
     # Global ID of this quad
     quad_id = offset + local_quad_id
 
@@ -296,13 +296,13 @@ function init_boundaries_iter_face_inner(info, boundaries, boundary_id, mesh)
     boundaries.neighbor_ids[boundary_id] = quad_id + 1
 
     # Face at which the boundary lies
-    face = side.face
+    face = side_pw.face[]
 
     # Save boundaries.node_indices dimension specific in containers_[23]d.jl
     init_boundary_node_indices!(boundaries, face, boundary_id)
 
     # One-based indexing
-    boundaries.name[boundary_id] = mesh.boundary_names[face + 1, side.treeid + 1]
+    boundaries.name[boundary_id] = mesh.boundary_names[face + 1, side_pw.treeid[] + 1]
 
     return nothing
 end
@@ -479,32 +479,33 @@ end
 # Function barrier for type stability
 function init_surfaces_iter_face_inner(info, user_data)
     @unpack interfaces, mortars, boundaries = user_data
-    elem_count = unsafe_load(info).sides.elem_count
+    info_pw = PointerWrapper(info)
+    elem_count = info_pw.sides.elem_count[]
 
     if elem_count == 2
         # Two neighboring elements => Interface or mortar
 
         # Extract surface data
-        sides = (unsafe_load_side(info, 1), unsafe_load_side(info, 2))
+        sides_pw = (load_pointerwrapper_side(info_pw, 1),
+                    load_pointerwrapper_side(info_pw, 2))
 
-        if sides[1].is_hanging == false && sides[2].is_hanging == false
+        if sides_pw[1].is_hanging[] == false && sides_pw[2].is_hanging[] == false
             # No hanging nodes => normal interface
             if interfaces !== nothing
-                init_interfaces_iter_face_inner(info, sides, user_data)
+                init_interfaces_iter_face_inner(info_pw, sides_pw, user_data)
             end
         else
             # Hanging nodes => mortar
             if mortars !== nothing
-                init_mortars_iter_face_inner(info, sides, user_data)
+                init_mortars_iter_face_inner(info_pw, sides_pw, user_data)
             end
         end
     elseif elem_count == 1
         # One neighboring elements => boundary
         if boundaries !== nothing
-            init_boundaries_iter_face_inner(info, user_data)
+            init_boundaries_iter_face_inner(info_pw, user_data)
         end
     end
-
     return nothing
 end
 
@@ -519,18 +520,18 @@ function init_surfaces!(interfaces, mortars, boundaries, mesh::P4estMesh)
 end
 
 # Initialization of interfaces after the function barrier
-function init_interfaces_iter_face_inner(info, sides, user_data)
+function init_interfaces_iter_face_inner(info_pw, sides_pw, user_data)
     @unpack interfaces, interface_id, mesh = user_data
     user_data.interface_id += 1
 
     # Get Tuple of local trees, one-based indexing
-    trees = (unsafe_load_tree(mesh.p4est, sides[1].treeid + 1),
-             unsafe_load_tree(mesh.p4est, sides[2].treeid + 1))
+    trees_pw = (load_pointerwrapper_tree(mesh.p4est, sides_pw[1].treeid[] + 1),
+                load_pointerwrapper_tree(mesh.p4est, sides_pw[2].treeid[] + 1))
     # Quadrant numbering offsets of the quadrants at this interface
-    offsets = SVector(trees[1].quadrants_offset,
-                      trees[2].quadrants_offset)
+    offsets = SVector(trees_pw[1].quadrants_offset[],
+                      trees_pw[2].quadrants_offset[])
 
-    local_quad_ids = SVector(sides[1].is.full.quadid, sides[2].is.full.quadid)
+    local_quad_ids = SVector(sides_pw[1].is.full.quadid[], sides_pw[2].is.full.quadid[])
     # Global IDs of the neighboring quads
     quad_ids = offsets + local_quad_ids
 
@@ -540,31 +541,30 @@ function init_interfaces_iter_face_inner(info, sides, user_data)
     interfaces.neighbor_ids[2, interface_id] = quad_ids[2] + 1
 
     # Face at which the interface lies
-    faces = (sides[1].face, sides[2].face)
+    faces = (sides_pw[1].face[], sides_pw[2].face[])
 
     # Save interfaces.node_indices dimension specific in containers_[23]d.jl
-    init_interface_node_indices!(interfaces, faces,
-                                 unsafe_load(info).orientation, interface_id)
+    init_interface_node_indices!(interfaces, faces, info_pw.orientation[], interface_id)
 
     return nothing
 end
 
 # Initialization of boundaries after the function barrier
-function init_boundaries_iter_face_inner(info, user_data)
+function init_boundaries_iter_face_inner(info_pw, user_data)
     @unpack boundaries, boundary_id, mesh = user_data
     user_data.boundary_id += 1
 
     # Extract boundary data
-    side = unsafe_load_side(info)
+    side_pw = load_pointerwrapper_side(info_pw)
     # Get local tree, one-based indexing
-    tree = unsafe_load_tree(mesh.p4est, side.treeid + 1)
+    tree_pw = load_pointerwrapper_tree(mesh.p4est, side_pw.treeid[] + 1)
     # Quadrant numbering offset of this quadrant
-    offset = tree.quadrants_offset
+    offset = tree_pw.quadrants_offset[]
 
     # Verify before accessing is.full, but this should never happen
-    @assert side.is_hanging == false
+    @assert side_pw.is_hanging[] == false
 
-    local_quad_id = side.is.full.quadid
+    local_quad_id = side_pw.is.full.quadid[]
     # Global ID of this quad
     quad_id = offset + local_quad_id
 
@@ -573,52 +573,52 @@ function init_boundaries_iter_face_inner(info, user_data)
     boundaries.neighbor_ids[boundary_id] = quad_id + 1
 
     # Face at which the boundary lies
-    face = side.face
+    face = side_pw.face[]
 
     # Save boundaries.node_indices dimension specific in containers_[23]d.jl
     init_boundary_node_indices!(boundaries, face, boundary_id)
 
     # One-based indexing
-    boundaries.name[boundary_id] = mesh.boundary_names[face + 1, side.treeid + 1]
+    boundaries.name[boundary_id] = mesh.boundary_names[face + 1, side_pw.treeid[] + 1]
 
     return nothing
 end
 
 # Initialization of mortars after the function barrier
-function init_mortars_iter_face_inner(info, sides, user_data)
+function init_mortars_iter_face_inner(info_pw, sides_pw, user_data)
     @unpack mortars, mortar_id, mesh = user_data
     user_data.mortar_id += 1
 
     # Get Tuple of local trees, one-based indexing
-    trees = (unsafe_load_tree(mesh.p4est, sides[1].treeid + 1),
-             unsafe_load_tree(mesh.p4est, sides[2].treeid + 1))
+    trees_pw = (load_pointerwrapper_tree(mesh.p4est, sides_pw[1].treeid[] + 1),
+                load_pointerwrapper_tree(mesh.p4est, sides_pw[2].treeid[] + 1))
     # Quadrant numbering offsets of the quadrants at this interface
-    offsets = SVector(trees[1].quadrants_offset,
-                      trees[2].quadrants_offset)
+    offsets = SVector(trees_pw[1].quadrants_offset[],
+                      trees_pw[2].quadrants_offset[])
 
-    if sides[1].is_hanging == true
+    if sides_pw[1].is_hanging[] == true
         # Left is small, right is large
-        faces = (sides[1].face, sides[2].face)
+        faces = (sides_pw[1].face[], sides_pw[2].face[])
 
-        local_small_quad_ids = sides[1].is.hanging.quadid
+        local_small_quad_ids = sides_pw[1].is.hanging.quadid[]
         # Global IDs of the two small quads
         small_quad_ids = offsets[1] .+ local_small_quad_ids
 
         # Just be sure before accessing is.full
-        @assert sides[2].is_hanging == false
-        large_quad_id = offsets[2] + sides[2].is.full.quadid
-    else # sides[2].is_hanging == true
+        @assert sides_pw[2].is_hanging[] == false
+        large_quad_id = offsets[2] + sides_pw[2].is.full.quadid[]
+    else # sides_pw[2].is_hanging[] == true
         # Right is small, left is large.
         # init_mortar_node_indices! below expects side 1 to contain the small elements.
-        faces = (sides[2].face, sides[1].face)
+        faces = (sides_pw[2].face[], sides_pw[1].face[])
 
-        local_small_quad_ids = sides[2].is.hanging.quadid
+        local_small_quad_ids = sides_pw[2].is.hanging.quadid[]
         # Global IDs of the two small quads
         small_quad_ids = offsets[2] .+ local_small_quad_ids
 
         # Just be sure before accessing is.full
-        @assert sides[1].is_hanging == false
-        large_quad_id = offsets[1] + sides[1].is.full.quadid
+        @assert sides_pw[1].is_hanging[] == false
+        large_quad_id = offsets[1] + sides_pw[1].is.full.quadid[]
     end
 
     # Write data to mortar container, 1 and 2 are the small elements
@@ -627,7 +627,7 @@ function init_mortars_iter_face_inner(info, sides, user_data)
     # Last entry is the large element
     mortars.neighbor_ids[end, mortar_id] = large_quad_id + 1
 
-    init_mortar_node_indices!(mortars, faces, unsafe_load(info).orientation, mortar_id)
+    init_mortar_node_indices!(mortars, faces, info_pw.orientation[], mortar_id)
 
     return nothing
 end
@@ -638,34 +638,36 @@ end
 # - boundaries
 # and collect the numbers in `user_data` in this order.
 function count_surfaces_iter_face(info, user_data)
-    elem_count = unsafe_load(info).sides.elem_count
+    info_pw = PointerWrapper(info)
+    elem_count = info_pw.sides.elem_count[]
 
     if elem_count == 2
         # Two neighboring elements => Interface or mortar
 
         # Extract surface data
-        sides = (unsafe_load_side(info, 1), unsafe_load_side(info, 2))
+        sides_pw = (load_pointerwrapper_side(info_pw, 1),
+                    load_pointerwrapper_side(info_pw, 2))
 
-        if sides[1].is_hanging == false && sides[2].is_hanging == false
+        if sides_pw[1].is_hanging[] == false && sides_pw[2].is_hanging[] == false
             # No hanging nodes => normal interface
             # Unpack user_data = [interface_count] and increment interface_count
-            ptr = Ptr{Int}(user_data)
-            id = unsafe_load(ptr, 1)
-            unsafe_store!(ptr, id + 1, 1)
+            pw = PointerWrapper(Int, user_data)
+            id = pw[1]
+            pw[1] = id + 1
         else
             # Hanging nodes => mortar
             # Unpack user_data = [mortar_count] and increment mortar_count
-            ptr = Ptr{Int}(user_data)
-            id = unsafe_load(ptr, 2)
-            unsafe_store!(ptr, id + 1, 2)
+            pw = PointerWrapper(Int, user_data)
+            id = pw[2]
+            pw[2] = id + 1
         end
     elseif elem_count == 1
         # One neighboring elements => boundary
 
         # Unpack user_data = [boundary_count] and increment boundary_count
-        ptr = Ptr{Int}(user_data)
-        id = unsafe_load(ptr, 3)
-        unsafe_store!(ptr, id + 1, 3)
+        pw = PointerWrapper(Int, user_data)
+        id = pw[3]
+        pw[3] = id + 1
     end
 
     return nothing
diff --git a/src/solvers/dgsem_p4est/containers_2d.jl b/src/solvers/dgsem_p4est/containers_2d.jl
index 4f7d903897a..11747f1f175 100644
--- a/src/solvers/dgsem_p4est/containers_2d.jl
+++ b/src/solvers/dgsem_p4est/containers_2d.jl
@@ -52,7 +52,7 @@ function calc_node_coordinates!(node_coordinates,
     p4est_root_len = 1 << P4EST_MAXLEVEL
     p4est_quadrant_len(l) = 1 << (P4EST_MAXLEVEL - l)
 
-    trees = unsafe_wrap_sc(p4est_tree_t, unsafe_load(mesh.p4est).trees)
+    trees = unsafe_wrap_sc(p4est_tree_t, mesh.p4est.trees)
 
     for tree in eachindex(trees)
         offset = trees[tree].quadrants_offset
diff --git a/src/solvers/dgsem_p4est/containers_3d.jl b/src/solvers/dgsem_p4est/containers_3d.jl
index 6cdc2cf9611..e9994fe4569 100644
--- a/src/solvers/dgsem_p4est/containers_3d.jl
+++ b/src/solvers/dgsem_p4est/containers_3d.jl
@@ -43,7 +43,7 @@ function calc_node_coordinates!(node_coordinates,
     p4est_root_len = 1 << P4EST_MAXLEVEL
     p4est_quadrant_len(l) = 1 << (P4EST_MAXLEVEL - l)
 
-    trees = unsafe_wrap_sc(p8est_tree_t, unsafe_load(mesh.p4est).trees)
+    trees = unsafe_wrap_sc(p8est_tree_t, mesh.p4est.trees)
 
     for tree in eachindex(trees)
         offset = trees[tree].quadrants_offset
diff --git a/src/solvers/dgsem_p4est/containers_parallel.jl b/src/solvers/dgsem_p4est/containers_parallel.jl
index 42d6ea44c5e..e7ee1f81478 100644
--- a/src/solvers/dgsem_p4est/containers_parallel.jl
+++ b/src/solvers/dgsem_p4est/containers_parallel.jl
@@ -311,21 +311,24 @@ function init_surfaces_iter_face_inner(info,
     # surfaces at once or any subset of them, some of the unpacked values above may be `nothing` if
     # they're not supposed to be initialized during this call. That is why we need additional
     # `!== nothing` checks below before initializing individual faces.
-    if unsafe_load(info).sides.elem_count == 2
+    info_pw = PointerWrapper(info)
+    if info_pw.sides.elem_count[] == 2
         # Two neighboring elements => Interface or mortar
 
         # Extract surface data
-        sides = (unsafe_load_side(info, 1), unsafe_load_side(info, 2))
+        sides_pw = (load_pointerwrapper_side(info_pw, 1),
+                    load_pointerwrapper_side(info_pw, 2))
 
-        if sides[1].is_hanging == false && sides[2].is_hanging == false
+        if sides_pw[1].is_hanging[] == false && sides_pw[2].is_hanging[] == false
             # No hanging nodes => normal interface or MPI interface
-            if sides[1].is.full.is_ghost == true || sides[2].is.full.is_ghost == true # remote side => MPI interface
+            if sides_pw[1].is.full.is_ghost[] == true ||
+               sides_pw[2].is.full.is_ghost[] == true # remote side => MPI interface
                 if mpi_interfaces !== nothing
-                    init_mpi_interfaces_iter_face_inner(info, sides, user_data)
+                    init_mpi_interfaces_iter_face_inner(info_pw, sides_pw, user_data)
                 end
             else
                 if interfaces !== nothing
-                    init_interfaces_iter_face_inner(info, sides, user_data)
+                    init_interfaces_iter_face_inner(info_pw, sides_pw, user_data)
                 end
             end
         else
@@ -333,18 +336,18 @@ function init_surfaces_iter_face_inner(info,
             # First, we check which side is hanging, i.e., on which side we have the refined cells.
             # Then we check if any of the refined cells or the coarse cell are "ghost" cells, i.e., they
             # belong to another rank. That way we can determine if this is a regular mortar or MPI mortar
-            if sides[1].is_hanging == true
-                @assert sides[2].is_hanging == false
-                if any(sides[1].is.hanging.is_ghost .== true) ||
-                   sides[2].is.full.is_ghost == true
+            if sides_pw[1].is_hanging[] == true
+                @assert sides_pw[2].is_hanging[] == false
+                if any(sides_pw[1].is.hanging.is_ghost[] .== true) ||
+                   sides_pw[2].is.full.is_ghost[] == true
                     face_has_ghost_side = true
                 else
                     face_has_ghost_side = false
                 end
-            else # sides[2].is_hanging == true
-                @assert sides[1].is_hanging == false
-                if sides[1].is.full.is_ghost == true ||
-                   any(sides[2].is.hanging.is_ghost .== true)
+            else # sides_pw[2].is_hanging[] == true
+                @assert sides_pw[1].is_hanging[] == false
+                if sides_pw[1].is.full.is_ghost[] == true ||
+                   any(sides_pw[2].is.hanging.is_ghost[] .== true)
                     face_has_ghost_side = true
                 else
                     face_has_ghost_side = false
@@ -352,15 +355,15 @@ function init_surfaces_iter_face_inner(info,
             end
             # Initialize mortar or MPI mortar
             if face_has_ghost_side && mpi_mortars !== nothing
-                init_mpi_mortars_iter_face_inner(info, sides, user_data)
+                init_mpi_mortars_iter_face_inner(info_pw, sides_pw, user_data)
             elseif !face_has_ghost_side && mortars !== nothing
-                init_mortars_iter_face_inner(info, sides, user_data)
+                init_mortars_iter_face_inner(info_pw, sides_pw, user_data)
             end
         end
-    elseif unsafe_load(info).sides.elem_count == 1
+    elseif info_pw.sides.elem_count[] == 1
         # One neighboring elements => boundary
         if boundaries !== nothing
-            init_boundaries_iter_face_inner(info, user_data)
+            init_boundaries_iter_face_inner(info_pw, user_data)
         end
     end
 
@@ -381,23 +384,23 @@ function init_surfaces!(interfaces, mortars, boundaries, mpi_interfaces, mpi_mor
 end
 
 # Initialization of MPI interfaces after the function barrier
-function init_mpi_interfaces_iter_face_inner(info, sides, user_data)
+function init_mpi_interfaces_iter_face_inner(info_pw, sides_pw, user_data)
     @unpack mpi_interfaces, mpi_interface_id, mesh = user_data
     user_data.mpi_interface_id += 1
 
-    if sides[1].is.full.is_ghost == true
+    if sides_pw[1].is.full.is_ghost[] == true
         local_side = 2
-    elseif sides[2].is.full.is_ghost == true
+    elseif sides_pw[2].is.full.is_ghost[] == true
         local_side = 1
     else
         error("should not happen")
     end
 
     # Get local tree, one-based indexing
-    tree = unsafe_load_tree(mesh.p4est, sides[local_side].treeid + 1)
+    tree_pw = load_pointerwrapper_tree(mesh.p4est, sides_pw[local_side].treeid[] + 1)
     # Quadrant numbering offset of the local quadrant at this interface
-    offset = tree.quadrants_offset
-    tree_quad_id = sides[local_side].is.full.quadid # quadid in the local tree
+    offset = tree_pw.quadrants_offset[]
+    tree_quad_id = sides_pw[local_side].is.full.quadid[] # quadid in the local tree
     # ID of the local neighboring quad, cumulative over local trees
     local_quad_id = offset + tree_quad_id
 
@@ -406,52 +409,52 @@ function init_mpi_interfaces_iter_face_inner(info, sides, user_data)
     mpi_interfaces.local_sides[mpi_interface_id] = local_side
 
     # Face at which the interface lies
-    faces = (sides[1].face, sides[2].face)
+    faces = (sides_pw[1].face[], sides_pw[2].face[])
 
     # Save mpi_interfaces.node_indices dimension specific in containers_[23]d_parallel.jl
     init_mpi_interface_node_indices!(mpi_interfaces, faces, local_side,
-                                     unsafe_load(info).orientation,
+                                     info_pw.orientation[],
                                      mpi_interface_id)
 
     return nothing
 end
 
 # Initialization of MPI mortars after the function barrier
-function init_mpi_mortars_iter_face_inner(info, sides, user_data)
+function init_mpi_mortars_iter_face_inner(info_pw, sides_pw, user_data)
     @unpack mpi_mortars, mpi_mortar_id, mesh = user_data
     user_data.mpi_mortar_id += 1
 
     # Get Tuple of adjacent trees, one-based indexing
-    trees = (unsafe_load_tree(mesh.p4est, sides[1].treeid + 1),
-             unsafe_load_tree(mesh.p4est, sides[2].treeid + 1))
+    trees_pw = (load_pointerwrapper_tree(mesh.p4est, sides_pw[1].treeid[] + 1),
+                load_pointerwrapper_tree(mesh.p4est, sides_pw[2].treeid[] + 1))
     # Quadrant numbering offsets of the quadrants at this mortar
-    offsets = SVector(trees[1].quadrants_offset,
-                      trees[2].quadrants_offset)
+    offsets = SVector(trees_pw[1].quadrants_offset[],
+                      trees_pw[2].quadrants_offset[])
 
-    if sides[1].is_hanging == true
+    if sides_pw[1].is_hanging[] == true
         hanging_side = 1
         full_side = 2
-    else # sides[2].is_hanging == true
+    else # sides_pw[2].is_hanging[] == true
         hanging_side = 2
         full_side = 1
     end
     # Just be sure before accessing is.full or is.hanging later
-    @assert sides[full_side].is_hanging == false
-    @assert sides[hanging_side].is_hanging == true
+    @assert sides_pw[full_side].is_hanging[] == false
+    @assert sides_pw[hanging_side].is_hanging[] == true
 
     # Find small quads that are locally available
-    local_small_quad_positions = findall(sides[hanging_side].is.hanging.is_ghost .==
+    local_small_quad_positions = findall(sides_pw[hanging_side].is.hanging.is_ghost[] .==
                                          false)
 
     # Get id of local small quadrants within their tree
     # Indexing CBinding.Caccessor via a Vector does not work here -> use map instead
-    tree_small_quad_ids = map(p -> sides[hanging_side].is.hanging.quadid[p],
+    tree_small_quad_ids = map(p -> sides_pw[hanging_side].is.hanging.quadid[][p],
                               local_small_quad_positions)
     local_small_quad_ids = offsets[hanging_side] .+ tree_small_quad_ids # ids cumulative over local trees
 
     # Determine if large quadrant is available and if yes, determine its id
-    if sides[full_side].is.full.is_ghost == false
-        local_large_quad_id = offsets[full_side] + sides[full_side].is.full.quadid
+    if sides_pw[full_side].is.full.is_ghost[] == false
+        local_large_quad_id = offsets[full_side] + sides_pw[full_side].is.full.quadid[]
     else
         local_large_quad_id = -1 # large quad is ghost
     end
@@ -470,9 +473,8 @@ function init_mpi_mortars_iter_face_inner(info, sides, user_data)
     mpi_mortars.local_neighbor_positions[mpi_mortar_id] = local_neighbor_positions
 
     # init_mortar_node_indices! expects side 1 to contain small elements
-    faces = (sides[hanging_side].face, sides[full_side].face)
-    init_mortar_node_indices!(mpi_mortars, faces, unsafe_load(info).orientation,
-                              mpi_mortar_id)
+    faces = (sides_pw[hanging_side].face[], sides_pw[full_side].face[])
+    init_mortar_node_indices!(mpi_mortars, faces, info_pw.orientation[], mpi_mortar_id)
 
     return nothing
 end
@@ -485,42 +487,45 @@ end
 # - (MPI) mortars at subdomain boundaries
 # and collect the numbers in `user_data` in this order.
 function count_surfaces_iter_face_parallel(info, user_data)
-    if unsafe_load(info).sides.elem_count == 2
+    info_pw = PointerWrapper(info)
+    if info_pw.sides.elem_count[] == 2
         # Two neighboring elements => Interface or mortar
 
         # Extract surface data
-        sides = (unsafe_load_side(info, 1), unsafe_load_side(info, 2))
+        sides_pw = (load_pointerwrapper_side(info_pw, 1),
+                    load_pointerwrapper_side(info_pw, 2))
 
-        if sides[1].is_hanging == false && sides[2].is_hanging == false
+        if sides_pw[1].is_hanging[] == false && sides_pw[2].is_hanging[] == false
             # No hanging nodes => normal interface or MPI interface
-            if sides[1].is.full.is_ghost == true || sides[2].is.full.is_ghost == true # remote side => MPI interface
+            if sides_pw[1].is.full.is_ghost[] == true ||
+               sides_pw[2].is.full.is_ghost[] == true # remote side => MPI interface
                 # Unpack user_data = [mpi_interface_count] and increment mpi_interface_count
-                ptr = Ptr{Int}(user_data)
-                id = unsafe_load(ptr, 4)
-                unsafe_store!(ptr, id + 1, 4)
+                pw = PointerWrapper(Int, user_data)
+                id = pw[4]
+                pw[4] = id + 1
             else
                 # Unpack user_data = [interface_count] and increment interface_count
-                ptr = Ptr{Int}(user_data)
-                id = unsafe_load(ptr, 1)
-                unsafe_store!(ptr, id + 1, 1)
+                pw = PointerWrapper(Int, user_data)
+                id = pw[1]
+                pw[1] = id + 1
             end
         else
             # Hanging nodes => mortar or MPI mortar
             # First, we check which side is hanging, i.e., on which side we have the refined cells.
             # Then we check if any of the refined cells or the coarse cell are "ghost" cells, i.e., they
             # belong to another rank. That way we can determine if this is a regular mortar or MPI mortar
-            if sides[1].is_hanging == true
-                @assert sides[2].is_hanging == false
-                if any(sides[1].is.hanging.is_ghost .== true) ||
-                   sides[2].is.full.is_ghost == true
+            if sides_pw[1].is_hanging[] == true
+                @assert sides_pw[2].is_hanging[] == false
+                if any(sides_pw[1].is.hanging.is_ghost[] .== true) ||
+                   sides_pw[2].is.full.is_ghost[] == true
                     face_has_ghost_side = true
                 else
                     face_has_ghost_side = false
                 end
-            else # sides[2].is_hanging == true
-                @assert sides[1].is_hanging == false
-                if sides[1].is.full.is_ghost == true ||
-                   any(sides[2].is.hanging.is_ghost .== true)
+            else # sides_pw[2].is_hanging[] == true
+                @assert sides_pw[1].is_hanging[] == false
+                if sides_pw[1].is.full.is_ghost[] == true ||
+                   any(sides_pw[2].is.hanging.is_ghost[] .== true)
                     face_has_ghost_side = true
                 else
                     face_has_ghost_side = false
@@ -528,23 +533,23 @@ function count_surfaces_iter_face_parallel(info, user_data)
             end
             if face_has_ghost_side
                 # Unpack user_data = [mpi_mortar_count] and increment mpi_mortar_count
-                ptr = Ptr{Int}(user_data)
-                id = unsafe_load(ptr, 5)
-                unsafe_store!(ptr, id + 1, 5)
+                pw = PointerWrapper(Int, user_data)
+                id = pw[5]
+                pw[5] = id + 1
             else
                 # Unpack user_data = [mortar_count] and increment mortar_count
-                ptr = Ptr{Int}(user_data)
-                id = unsafe_load(ptr, 2)
-                unsafe_store!(ptr, id + 1, 2)
+                pw = PointerWrapper(Int, user_data)
+                id = pw[2]
+                pw[2] = id + 1
             end
         end
-    elseif unsafe_load(info).sides.elem_count == 1
+    elseif info_pw.sides.elem_count[] == 1
         # One neighboring elements => boundary
 
         # Unpack user_data = [boundary_count] and increment boundary_count
-        ptr = Ptr{Int}(user_data)
-        id = unsafe_load(ptr, 3)
-        unsafe_store!(ptr, id + 1, 3)
+        pw = PointerWrapper(Int, user_data)
+        id = pw[3]
+        pw[3] = id + 1
     end
 
     return nothing
diff --git a/src/solvers/dgsem_p4est/dg_parallel.jl b/src/solvers/dgsem_p4est/dg_parallel.jl
index ac122d048c1..324bc7f3cd6 100644
--- a/src/solvers/dgsem_p4est/dg_parallel.jl
+++ b/src/solvers/dgsem_p4est/dg_parallel.jl
@@ -263,15 +263,13 @@ function init_mpi_cache!(mpi_cache::P4estMPICache, mesh::ParallelP4estMesh,
                                                                                                         uEltype)
 
     # Determine local and total number of elements
-    n_elements_global = Int(unsafe_load(mesh.p4est).global_num_quadrants)
-    n_elements_by_rank = vcat(Int.(unsafe_wrap(Array,
-                                               unsafe_load(mesh.p4est).global_first_quadrant,
+    n_elements_global = Int(mesh.p4est.global_num_quadrants[])
+    n_elements_by_rank = vcat(Int.(unsafe_wrap(Array, mesh.p4est.global_first_quadrant,
                                                mpi_nranks())),
                               n_elements_global) |> diff # diff sufficient due to 0-based quad indices
     n_elements_by_rank = OffsetArray(n_elements_by_rank, 0:(mpi_nranks() - 1))
     # Account for 1-based indexing in Julia
-    first_element_global_id = Int(unsafe_load(unsafe_load(mesh.p4est).global_first_quadrant,
-                                              mpi_rank() + 1)) + 1
+    first_element_global_id = Int(mesh.p4est.global_first_quadrant[mpi_rank() + 1]) + 1
     @assert n_elements_global==sum(n_elements_by_rank) "error in total number of elements"
 
     # TODO reuse existing structures
@@ -379,17 +377,19 @@ function init_neighbor_rank_connectivity_iter_face_inner(info, user_data)
     @unpack interfaces, interface_id, global_interface_ids, neighbor_ranks_interface,
     mortars, mortar_id, global_mortar_ids, neighbor_ranks_mortar, mesh = user_data
 
+    info_pw = PointerWrapper(info)
     # Get the global interface/mortar ids and neighbor rank if current face belongs to an MPI
     # interface/mortar
-    if unsafe_load(info).sides.elem_count == 2 # MPI interfaces/mortars have two neighboring elements
+    if info_pw.sides.elem_count[] == 2 # MPI interfaces/mortars have two neighboring elements
         # Extract surface data
-        sides = (unsafe_load_side(info, 1), unsafe_load_side(info, 2))
+        sides_pw = (load_pointerwrapper_side(info_pw, 1),
+                    load_pointerwrapper_side(info_pw, 2))
 
-        if sides[1].is_hanging == false && sides[2].is_hanging == false # No hanging nodes for MPI interfaces
-            if sides[1].is.full.is_ghost == true
+        if sides_pw[1].is_hanging[] == false && sides_pw[2].is_hanging[] == false # No hanging nodes for MPI interfaces
+            if sides_pw[1].is.full.is_ghost[] == true
                 remote_side = 1
                 local_side = 2
-            elseif sides[2].is.full.is_ghost == true
+            elseif sides_pw[2].is.full.is_ghost[] == true
                 remote_side = 2
                 local_side = 1
             else # both sides are on this rank -> skip since it's a regular interface
@@ -397,16 +397,17 @@ function init_neighbor_rank_connectivity_iter_face_inner(info, user_data)
             end
 
             # Sanity check, current face should belong to current MPI interface
-            local_tree = unsafe_load_tree(mesh.p4est, sides[local_side].treeid + 1) # one-based indexing
-            local_quad_id = local_tree.quadrants_offset +
-                            sides[local_side].is.full.quadid
+            local_tree_pw = load_pointerwrapper_tree(mesh.p4est,
+                                                     sides_pw[local_side].treeid[] + 1) # one-based indexing
+            local_quad_id = local_tree_pw.quadrants_offset[] +
+                            sides_pw[local_side].is.full.quadid[]
             @assert interfaces.local_neighbor_ids[interface_id] == local_quad_id + 1 # one-based indexing
 
             # Get neighbor ID from ghost layer
             proc_offsets = unsafe_wrap(Array,
-                                       unsafe_load(unsafe_load(info).ghost_layer).proc_offsets,
+                                       info_pw.ghost_layer.proc_offsets,
                                        mpi_nranks() + 1)
-            ghost_id = sides[remote_side].is.full.quadid # indexes the ghost layer, 0-based
+            ghost_id = sides_pw[remote_side].is.full.quadid[] # indexes the ghost layer, 0-based
             neighbor_rank = findfirst(r -> proc_offsets[r] <= ghost_id <
                                            proc_offsets[r + 1],
                                       1:mpi_nranks()) - 1 # MPI ranks are 0-based
@@ -415,21 +416,18 @@ function init_neighbor_rank_connectivity_iter_face_inner(info, user_data)
             # Global interface id is the globally unique quadrant id of the quadrant on the primary
             # side (1) multiplied by the number of faces per quadrant plus face
             if local_side == 1
-                offset = unsafe_load(unsafe_load(mesh.p4est).global_first_quadrant,
-                                     mpi_rank() + 1) # one-based indexing
+                offset = mesh.p4est.global_first_quadrant[mpi_rank() + 1] # one-based indexing
                 primary_quad_id = offset + local_quad_id
             else
-                offset = unsafe_load(unsafe_load(mesh.p4est).global_first_quadrant,
-                                     neighbor_rank + 1) # one-based indexing
-                primary_quad_id = offset +
-                                  unsafe_load(sides[1].is.full.quad.p.piggy3.local_num)
+                offset = mesh.p4est.global_first_quadrant[neighbor_rank + 1] # one-based indexing
+                primary_quad_id = offset + sides_pw[1].is.full.quad.p.piggy3.local_num[]
             end
-            global_interface_id = 2 * ndims(mesh) * primary_quad_id + sides[1].face
+            global_interface_id = 2 * ndims(mesh) * primary_quad_id + sides_pw[1].face[]
             global_interface_ids[interface_id] = global_interface_id
 
             user_data.interface_id += 1
         else # hanging node
-            if sides[1].is_hanging == true
+            if sides_pw[1].is_hanging[] == true
                 hanging_side = 1
                 full_side = 2
             else
@@ -437,26 +435,26 @@ function init_neighbor_rank_connectivity_iter_face_inner(info, user_data)
                 full_side = 1
             end
             # Verify before accessing is.full / is.hanging
-            @assert sides[hanging_side].is_hanging == true &&
-                    sides[full_side].is_hanging == false
+            @assert sides_pw[hanging_side].is_hanging[] == true &&
+                    sides_pw[full_side].is_hanging[] == false
 
             # If all quadrants are locally available, this is a regular mortar -> skip
-            if sides[full_side].is.full.is_ghost == false &&
-               all(sides[hanging_side].is.hanging.is_ghost .== false)
+            if sides_pw[full_side].is.full.is_ghost[] == false &&
+               all(sides_pw[hanging_side].is.hanging.is_ghost[] .== false)
                 return nothing
             end
 
-            trees = (unsafe_load_tree(mesh.p4est, sides[1].treeid + 1),
-                     unsafe_load_tree(mesh.p4est, sides[2].treeid + 1))
+            trees_pw = (load_pointerwrapper_tree(mesh.p4est, sides_pw[1].treeid[] + 1),
+                        load_pointerwrapper_tree(mesh.p4est, sides_pw[2].treeid[] + 1))
 
             # Find small quads that are remote and determine which rank owns them
-            remote_small_quad_positions = findall(sides[hanging_side].is.hanging.is_ghost .==
+            remote_small_quad_positions = findall(sides_pw[hanging_side].is.hanging.is_ghost[] .==
                                                   true)
             proc_offsets = unsafe_wrap(Array,
-                                       unsafe_load(unsafe_load(info).ghost_layer).proc_offsets,
+                                       info_pw.ghost_layer.proc_offsets,
                                        mpi_nranks() + 1)
             # indices of small remote quads inside the ghost layer, 0-based
-            ghost_ids = map(pos -> sides[hanging_side].is.hanging.quadid[pos],
+            ghost_ids = map(pos -> sides_pw[hanging_side].is.hanging.quadid[][pos],
                             remote_small_quad_positions)
             neighbor_ranks = map(ghost_ids) do ghost_id
                 return findfirst(r -> proc_offsets[r] <= ghost_id < proc_offsets[r + 1],
@@ -464,28 +462,26 @@ function init_neighbor_rank_connectivity_iter_face_inner(info, user_data)
             end
             # Determine global quad id of large element to determine global MPI mortar id
             # Furthermore, if large element is ghost, add its owner rank to neighbor_ranks
-            if sides[full_side].is.full.is_ghost == true
-                ghost_id = sides[full_side].is.full.quadid
+            if sides_pw[full_side].is.full.is_ghost[] == true
+                ghost_id = sides_pw[full_side].is.full.quadid[]
                 large_quad_owner_rank = findfirst(r -> proc_offsets[r] <= ghost_id <
                                                        proc_offsets[r + 1],
                                                   1:mpi_nranks()) - 1 # MPI ranks are 0-based
                 push!(neighbor_ranks, large_quad_owner_rank)
 
-                offset = unsafe_load(unsafe_load(mesh.p4est).global_first_quadrant,
-                                     large_quad_owner_rank + 1) # one-based indexing
+                offset = mesh.p4est.global_first_quadrant[large_quad_owner_rank + 1] # one-based indexing
                 large_quad_id = offset +
-                                unsafe_load(sides[full_side].is.full.quad.p.piggy3.local_num)
+                                sides_pw[full_side].is.full.quad.p.piggy3.local_num[]
             else
-                offset = unsafe_load(unsafe_load(mesh.p4est).global_first_quadrant,
-                                     mpi_rank() + 1) # one-based indexing
-                large_quad_id = offset + trees[full_side].quadrants_offset +
-                                sides[full_side].is.full.quadid
+                offset = mesh.p4est.global_first_quadrant[mpi_rank() + 1] # one-based indexing
+                large_quad_id = offset + trees_pw[full_side].quadrants_offset[] +
+                                sides_pw[full_side].is.full.quadid[]
             end
             neighbor_ranks_mortar[mortar_id] = neighbor_ranks
             # Global mortar id is the globally unique quadrant id of the large quadrant multiplied by the
             # number of faces per quadrant plus face
             global_mortar_ids[mortar_id] = 2 * ndims(mesh) * large_quad_id +
-                                           sides[full_side].face
+                                           sides_pw[full_side].face[]
 
             user_data.mortar_id += 1
         end

From eb4c91a6ba3c1fe1c1f0bbb9c332ab0068171113 Mon Sep 17 00:00:00 2001
From: David Knapp <david.knapp@dlr.de>
Date: Mon, 19 Jun 2023 10:48:56 +0200
Subject: [PATCH 2/3] Added a dispatchable constructur for DG with
 TensorProductWedges (#1389)

* Added a dispatchable constructur for DG with TensorProductWedges

Adapted the timestepping for DGMulti to use the minimal polynomial degree of the tensorproduct

* Use the maximal poly-degree

* Update src/solvers/dgmulti/types.jl

Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com>

* Change DGMulti-call for Tensorwedges

* Adding comment about the polydeg-choice

* Adapt constructor

* Add an example for TensorWedges

* Update example

* Update examples/dgmulti_3d/elixir_euler_tensorWedge.jl

Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com>

* Update examples/dgmulti_3d/elixir_euler_tensorWedge.jl

Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com>

* Update examples/dgmulti_3d/elixir_euler_tensorWedge.jl

Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com>

* Readd Gauss

Accidentaly deleted during Merge-Conflict resolvement online

* Remove explicit load of StartUpDG in elixir

* Update max_dt

* Apply suggestions from code review

Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com>

* Update comments

* Add tensorWedge-test

* Update src/solvers/dgmulti/dg.jl

Co-authored-by: Jesse Chan <1156048+jlchan@users.noreply.github.com>

* Renamed tensorwedge elixir

* Update testset title

* Remove typos

* temporal push

* fix max_dt dispatch

* Update TensorProductWedges-test

* Update test

* Address CI-Warning

* Update dt_polydeg_scaling

* Correct dt_polydeg_scaling

A previous commit changed it from (polydeg + 1) to polydeg

* Update DGMulti-call

* Rename file

* Change Path in test-file

* Apply suggestions from code review

Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com>

* Update examples/dgmulti_3d/elixir_advection_tensor_wedge.jl

Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com>

* Reformat src/Trixi.jl

* Format solvers/dgmulti/dg.jl

* Formatting

---------

Co-authored-by: Knapp <knap_da@sc-030182l.intra.dlr.de>
Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com>
Co-authored-by: Jesse Chan <1156048+jlchan@users.noreply.github.com>
Co-authored-by: Jesse Chan <jesse.chan@rice.edu>
---
 .../elixir_advection_tensor_wedge.jl          | 56 +++++++++++++++++++
 src/Trixi.jl                                  |  3 +-
 src/solvers/dgmulti/dg.jl                     | 11 ++--
 src/solvers/dgmulti/types.jl                  | 15 +++++
 test/test_dgmulti_3d.jl                       |  6 ++
 5 files changed, 86 insertions(+), 5 deletions(-)
 create mode 100644 examples/dgmulti_3d/elixir_advection_tensor_wedge.jl

diff --git a/examples/dgmulti_3d/elixir_advection_tensor_wedge.jl b/examples/dgmulti_3d/elixir_advection_tensor_wedge.jl
new file mode 100644
index 00000000000..4f43f2571a3
--- /dev/null
+++ b/examples/dgmulti_3d/elixir_advection_tensor_wedge.jl
@@ -0,0 +1,56 @@
+using OrdinaryDiffEq
+using Trixi
+using LinearAlgebra
+
+###############################################################################
+equations = LinearScalarAdvectionEquation3D(1.0, 1.0, 1.0)
+
+initial_condition = initial_condition_convergence_test
+
+# Define the polynomial degrees for the polynoms of the triangular base and the line
+# of the tensor-prism
+tensor_polydeg = (3, 4)
+
+dg = DGMulti(element_type = Wedge(),
+             approximation_type = Polynomial(),
+             surface_flux = flux_lax_friedrichs,
+             polydeg = tensor_polydeg)
+
+
+cells_per_dimension = (8, 8, 8)
+mesh = DGMultiMesh(dg, 
+                   cells_per_dimension,
+                   coordinates_min = (-1.0, -1.0, -1.0), 
+                   coordinates_max = (1.0, 1.0, 1.0),
+                   periodicity = true)
+
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, dg,
+                boundary_conditions=boundary_condition_periodic)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 5.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval=analysis_interval, uEltype=real(dg))
+
+alive_callback = AliveCallback(analysis_interval=analysis_interval)
+
+# The StepsizeCallback handles the re-calculation of the maximum Δt after each time step
+stepsize_callback = StepsizeCallback(cfl=1.0)
+
+callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, stepsize_callback)
+
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), dt = 1.0,
+            save_everystep=false, callback=callbacks);
+
+summary_callback() # print the timer summary
\ No newline at end of file
diff --git a/src/Trixi.jl b/src/Trixi.jl
index 86e349c7dad..66878f4b459 100644
--- a/src/Trixi.jl
+++ b/src/Trixi.jl
@@ -83,7 +83,8 @@ import SummationByPartsOperators: integrate, semidiscretize,
                                            upwind_operators
 
 # DGMulti solvers
-@reexport using StartUpDG: StartUpDG, Polynomial, Gauss, SBP, Line, Tri, Quad, Hex, Tet
+@reexport using StartUpDG: StartUpDG, Polynomial, Gauss, TensorProductWedge, SBP, Line, Tri,
+                           Quad, Hex, Tet, Wedge
 using StartUpDG: RefElemData, MeshData, AbstractElemShape
 
 # TODO: include_optimized
diff --git a/src/solvers/dgmulti/dg.jl b/src/solvers/dgmulti/dg.jl
index d51c7cabf9d..bc76aa1a9d2 100644
--- a/src/solvers/dgmulti/dg.jl
+++ b/src/solvers/dgmulti/dg.jl
@@ -226,6 +226,11 @@ function estimate_dt(mesh::DGMultiMesh, dg::DGMulti)
     return StartUpDG.estimate_h(rd, mesh.md) / StartUpDG.inverse_trace_constant(rd)
 end
 
+dt_polydeg_scaling(dg::DGMulti) = inv(dg.basis.N + 1)
+function dt_polydeg_scaling(dg::DGMulti{3, <:Wedge, <:TensorProductWedge})
+    inv(maximum(dg.basis.N) + 1)
+end
+
 # for the stepsize callback
 function max_dt(u, t, mesh::DGMultiMesh,
                 constant_speed::False, equations, dg::DGMulti{NDIMS},
@@ -247,8 +252,7 @@ function max_dt(u, t, mesh::DGMultiMesh,
     # `polydeg+1`. This is because `nnodes(dg)` returns the total number of
     # multi-dimensional nodes for DGMulti solver types, while `nnodes(dg)` returns
     # the number of 1D nodes for `DGSEM` solvers.
-    polydeg = rd.N
-    return 2 * dt_min / (polydeg + 1)
+    return 2 * dt_min * dt_polydeg_scaling(dg)
 end
 
 function max_dt(u, t, mesh::DGMultiMesh,
@@ -270,8 +274,7 @@ function max_dt(u, t, mesh::DGMultiMesh,
     # `polydeg+1`. This is because `nnodes(dg)` returns the total number of
     # multi-dimensional nodes for DGMulti solver types, while `nnodes(dg)` returns
     # the number of 1D nodes for `DGSEM` solvers.
-    polydeg = rd.N
-    return 2 * dt_min / (polydeg + 1)
+    return 2 * dt_min * dt_polydeg_scaling(dg)
 end
 
 # interpolates from solution coefficients to face quadrature points
diff --git a/src/solvers/dgmulti/types.jl b/src/solvers/dgmulti/types.jl
index c225e334e8e..f1f7b158dec 100644
--- a/src/solvers/dgmulti/types.jl
+++ b/src/solvers/dgmulti/types.jl
@@ -96,6 +96,21 @@ function DGMulti(; polydeg = nothing,
             polydeg = polydeg, kwargs...)
 end
 
+# dispatchable constructor for DGMulti using a TensorProductWedge
+function DGMulti(element_type::Wedge,
+                 approximation_type,
+                 volume_integral,
+                 surface_integral;
+                 polydeg::Tuple,
+                 kwargs...)
+    factor_a = RefElemData(Tri(), approximation_type, polydeg[1]; kwargs...)
+    factor_b = RefElemData(Line(), approximation_type, polydeg[2]; kwargs...)
+
+    tensor = TensorProductWedge(factor_a, factor_b)
+    rd = RefElemData(element_type, tensor; kwargs...)
+    return DG(rd, nothing, surface_integral, volume_integral)
+end
+
 # dispatchable constructor for DGMulti to allow for specialization
 function DGMulti(element_type::AbstractElemShape,
                  approximation_type,
diff --git a/test/test_dgmulti_3d.jl b/test/test_dgmulti_3d.jl
index 22c0a0fd3ba..68fa1d13304 100644
--- a/test/test_dgmulti_3d.jl
+++ b/test/test_dgmulti_3d.jl
@@ -135,6 +135,12 @@ isdir(outdir) && rm(outdir, recursive=true)
     )
   end
 
+  @trixi_testset "elixir_advection_tensor_wedge.jl" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_tensor_wedge.jl"),
+    l2 = [2.30487910e-04] ,
+    linf = [6.31795281e-04] )
+  end
+
 end
 
 # Clean up afterwards: delete Trixi.jl output directory

From a837dd7726f25a88c7bab8bed2dc6955f20534bf Mon Sep 17 00:00:00 2001
From: Andrew Winters <andrew.ross.winters@liu.se>
Date: Mon, 19 Jun 2023 14:50:27 +0200
Subject: [PATCH 3/3] Update test values for 1D discontinuous examples (#1511)

* set true discontinuous initial conditions for 1D

* update Burgers tests

* update MHD tests in 1D

* update structured_1d_dgsem/elixir_advection_shockcapturing elixir and test

* remove workaround in elixir_shallowwater_ec.jl and update test

* accidentally removed adjusted tolerance in structured mesh test

* remove workaround in elixir_shallowwater_shockcapturing.jl and update test

* remove workaround in elixir_shallowwater_well_balanced.jl and update tests

* bug fix in energy_total computation for ShallowWaterTwoLayerEquations1D

* remove workaround in elixir_shallowwater_twolayer_dam_break.jl and update test

* update docs test for adding a nonconservative equation

* update Nonconservative terms in 1D (linear advection) test

* update all necessary IdealGlmMhdMulticomponentEquations1D tests

* update all necessary CompressibleEulerMulticomponentEquations1D tests

* update necessary CompressibleEuler1D tests except for neural network shock capturing

* Update test values for NN-based SC

* add short description to the NEWS.md

* add extra test to avoid coverage drop

* Update src/solvers/dg.jl

Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com>

---------

Co-authored-by: Michael Schlottke-Lakemper <michael@sloede.com>
Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com>
---
 NEWS.md                                       |  1 +
 .../files/adding_nonconservative_equation.jl  |  2 +-
 .../elixir_advection_shockcapturing.jl        |  6 +-
 .../tree_1d_dgsem/elixir_burgers_shock.jl     |  6 +-
 ..._eulermulti_two_interacting_blast_waves.jl |  2 +-
 .../elixir_mhd_briowu_shock_tube.jl           |  2 +-
 examples/tree_1d_dgsem/elixir_mhdmulti_ec.jl  |  4 +-
 examples/tree_1d_dgsem/elixir_mhdmulti_es.jl  |  4 +-
 .../tree_1d_dgsem/elixir_shallowwater_ec.jl   | 71 ++++++---------
 .../elixir_shallowwater_shock_capturing.jl    | 81 +++++++----------
 .../elixir_shallowwater_twolayer_dam_break.jl | 77 +++++++----------
 .../elixir_shallowwater_well_balanced.jl      | 56 +++---------
 src/equations/shallow_water_two_layer_1d.jl   | 86 +++++++++----------
 src/solvers/dg.jl                             |  9 ++
 test/test_structured_1d.jl                    |  4 +-
 test/test_tree_1d.jl                          |  2 +-
 test/test_tree_1d_burgers.jl                  |  8 +-
 test/test_tree_1d_euler.jl                    | 44 +++++-----
 test/test_tree_1d_eulermulti.jl               | 32 ++++---
 test/test_tree_1d_mhd.jl                      | 12 +--
 test/test_tree_1d_mhdmulti.jl                 | 48 +++++------
 test/test_tree_1d_shallowwater.jl             | 28 +++---
 test/test_tree_1d_shallowwater_twolayer.jl    |  6 +-
 23 files changed, 267 insertions(+), 324 deletions(-)

diff --git a/NEWS.md b/NEWS.md
index 9b46ba565fe..35c7039b2ef 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -9,6 +9,7 @@ for human readability.
 #### Added
 
 - Experimental support for 3D parabolic diffusion terms has been added.
+- Capability to set truly discontinuous initial conditions in 1D.
 
 #### Changed
 
diff --git a/docs/literate/src/files/adding_nonconservative_equation.jl b/docs/literate/src/files/adding_nonconservative_equation.jl
index 08dd631058e..110fa486070 100644
--- a/docs/literate/src/files/adding_nonconservative_equation.jl
+++ b/docs/literate/src/files/adding_nonconservative_equation.jl
@@ -147,7 +147,7 @@ plot(sol)
 # above.
 
 error_1 = analysis_callback(sol).l2 |> first
-@test isapprox(error_1, 0.0002961027497) #src
+@test isapprox(error_1, 0.00029609575838969394) #src
 # Next, we increase the grid resolution by one refinement level and run the
 # simulation again.
 
diff --git a/examples/structured_1d_dgsem/elixir_advection_shockcapturing.jl b/examples/structured_1d_dgsem/elixir_advection_shockcapturing.jl
index 9a81acfe51c..313812fe08d 100644
--- a/examples/structured_1d_dgsem/elixir_advection_shockcapturing.jl
+++ b/examples/structured_1d_dgsem/elixir_advection_shockcapturing.jl
@@ -9,7 +9,9 @@ advection_velocity = 1.0
 """
     initial_condition_composite(x, t, equations::LinearScalarAdvectionEquation1D)
 
-Slight simplification of
+Wave form that is a combination of a Gaussian pulse, a square wave, a triangle wave,
+and half an ellipse with periodic boundary conditions.
+Slight simplification from
 - Jiang, Shu (1996)
   Efficient Implementation of Weighted ENO Schemes
   [DOI: 10.1006/jcph.1996.0130](https://doi.org/10.1006/jcph.1996.0130)
@@ -60,7 +62,7 @@ volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
 solver = DGSEM(basis, surface_flux, volume_integral)
 
 # Create curved mesh
-cells_per_dimension = (125,)
+cells_per_dimension = (120,)
 coordinates_min = (-1.0,) # minimum coordinate
 coordinates_max = (1.0,)  # maximum coordinate
 mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max,
diff --git a/examples/tree_1d_dgsem/elixir_burgers_shock.jl b/examples/tree_1d_dgsem/elixir_burgers_shock.jl
index 987fb320ad6..00b5314e19f 100644
--- a/examples/tree_1d_dgsem/elixir_burgers_shock.jl
+++ b/examples/tree_1d_dgsem/elixir_burgers_shock.jl
@@ -21,7 +21,7 @@ surface_flux = flux_lax_friedrichs
 volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
                                                  volume_flux_dg=surface_flux,
                                                  volume_flux_fv=surface_flux)
-                                                 
+
 solver = DGSEM(basis, surface_flux, volume_integral)
 
 coordinate_min = 0.0
@@ -59,7 +59,7 @@ end
 
 boundary_conditions = (x_neg=boundary_condition_inflow,
                        x_pos=boundary_condition_outflow)
-                       
+
 initial_condition = initial_condition_shock
 
 semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
@@ -79,7 +79,7 @@ analysis_callback = AnalysisCallback(semi, interval=analysis_interval)
 
 alive_callback = AliveCallback(analysis_interval=analysis_interval)
 
-stepsize_callback = StepsizeCallback(cfl=0.9)
+stepsize_callback = StepsizeCallback(cfl=0.85)
 
 
 callbacks = CallbackSet(summary_callback,
diff --git a/examples/tree_1d_dgsem/elixir_eulermulti_two_interacting_blast_waves.jl b/examples/tree_1d_dgsem/elixir_eulermulti_two_interacting_blast_waves.jl
index 353093e5f70..81966194180 100644
--- a/examples/tree_1d_dgsem/elixir_eulermulti_two_interacting_blast_waves.jl
+++ b/examples/tree_1d_dgsem/elixir_eulermulti_two_interacting_blast_waves.jl
@@ -88,7 +88,7 @@ ode = semidiscretize(semi, tspan)
 
 summary_callback = SummaryCallback()
 
-analysis_interval = 100
+analysis_interval = 1000
 
 analysis_callback = AnalysisCallback(semi, interval=analysis_interval)
 
diff --git a/examples/tree_1d_dgsem/elixir_mhd_briowu_shock_tube.jl b/examples/tree_1d_dgsem/elixir_mhd_briowu_shock_tube.jl
index 1c07fc4fdde..c5727109d92 100644
--- a/examples/tree_1d_dgsem/elixir_mhd_briowu_shock_tube.jl
+++ b/examples/tree_1d_dgsem/elixir_mhd_briowu_shock_tube.jl
@@ -94,7 +94,7 @@ amr_callback = AMRCallback(semi, amr_controller,
                            adapt_initial_condition=true,
                            adapt_initial_condition_only_refine=true)
 
-stepsize_callback = StepsizeCallback(cfl=0.8)
+stepsize_callback = StepsizeCallback(cfl=0.65)
 
 callbacks = CallbackSet(summary_callback,
                         analysis_callback, alive_callback,
diff --git a/examples/tree_1d_dgsem/elixir_mhdmulti_ec.jl b/examples/tree_1d_dgsem/elixir_mhdmulti_ec.jl
index 34fdce6634e..69ea0551bed 100644
--- a/examples/tree_1d_dgsem/elixir_mhdmulti_ec.jl
+++ b/examples/tree_1d_dgsem/elixir_mhdmulti_ec.jl
@@ -4,8 +4,8 @@ using Trixi
 
 ###############################################################################
 # semidiscretization of the ideal MHD equations
-equations = IdealGlmMhdMulticomponentEquations1D(gammas           = (2.0, 2.0, 2.0),
-                                                 gas_constants    = (2.0, 2.0, 2.0))
+equations = IdealGlmMhdMulticomponentEquations1D(gammas        = (2.0, 2.0, 2.0),
+                                                 gas_constants = (2.0, 2.0, 2.0))
 
 initial_condition = initial_condition_weak_blast_wave
 
diff --git a/examples/tree_1d_dgsem/elixir_mhdmulti_es.jl b/examples/tree_1d_dgsem/elixir_mhdmulti_es.jl
index 8ca32194b9e..93cf3e0fdb2 100644
--- a/examples/tree_1d_dgsem/elixir_mhdmulti_es.jl
+++ b/examples/tree_1d_dgsem/elixir_mhdmulti_es.jl
@@ -4,8 +4,8 @@ using Trixi
 
 ###############################################################################
 # semidiscretization of the ideal MHD equations
-equations = IdealGlmMhdMulticomponentEquations1D(gammas           = (2.0, 2.0, 2.0),
-                                                 gas_constants    = (2.0, 2.0, 2.0))
+equations = IdealGlmMhdMulticomponentEquations1D(gammas        = (2.0, 2.0, 2.0),
+                                                 gas_constants = (2.0, 2.0, 2.0))
 
 initial_condition = initial_condition_weak_blast_wave
 
diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_ec.jl b/examples/tree_1d_dgsem/elixir_shallowwater_ec.jl
index be6a2cb166c..1469afec1ca 100644
--- a/examples/tree_1d_dgsem/elixir_shallowwater_ec.jl
+++ b/examples/tree_1d_dgsem/elixir_shallowwater_ec.jl
@@ -8,9 +8,34 @@ using Trixi
 
 equations = ShallowWaterEquations1D(gravity_constant=9.81)
 
-# Note, this initial condition is used to compute errors in the analysis callback but the initialization is
-# overwritten by `initial_condition_ec_discontinuous_bottom` below.
-initial_condition = initial_condition_weak_blast_wave
+# Initial condition with a truly discontinuous water height, velocity, and bottom
+# topography function as an academic testcase for entropy conservation.
+# The errors from the analysis callback are not important but `∑∂S/∂U ⋅ Uₜ` should
+# be around machine roundoff.
+# Works as intended for TreeMesh1D with `initial_refinement_level=4`. If the mesh
+# refinement level is changed the initial condition below may need changed as well to
+# ensure that the discontinuities lie on an element interface.
+function initial_condition_ec_discontinuous_bottom(x, t, equations::ShallowWaterEquations1D)
+  # Set the background values
+  H = 4.25
+  v = 0.0
+  b = sin(x[1]) # arbitrary continuous function
+
+  # Setup the discontinuous water height and velocity
+  if x[1] >= 0.125 && x[1] <= 0.25
+    H = 5.0
+    v = 0.1882
+  end
+
+  # Setup a discontinuous bottom topography
+  if x[1] >= -0.25 && x[1] <= -0.125
+    b = 2.0 + 0.5 * sin(2.0 * pi * x[1])
+  end
+
+  return prim2cons(SVector(H, v, b), equations)
+end
+
+initial_condition = initial_condition_ec_discontinuous_bottom
 
 ###############################################################################
 # Get the DG approximation space
@@ -37,46 +62,6 @@ semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
 tspan = (0.0, 2.0)
 ode = semidiscretize(semi, tspan)
 
-###############################################################################
-# Workaround to set a discontinuous bottom topography and initial condition for debugging and testing.
-
-# alternative version of the initial conditinon used to setup a truly discontinuous
-# bottom topography function and initial condition for this academic testcase of entropy conservation.
-# The errors from the analysis callback are not important but `∑∂S/∂U ⋅ Uₜ` should be around machine roundoff
-# In contrast to the usual signature of initial conditions, this one get passed the
-# `element_id` explicitly. In particular, this initial conditions works as intended
-# only for the TreeMesh1D with `initial_refinement_level=4`.
-function initial_condition_ec_discontinuous_bottom(x, t, element_id, equations::ShallowWaterEquations1D)
-  # Set the background values
-  H = 4.25
-  v = 0.0
-  b = sin(x[1]) # arbitrary continuous function
-
-  # setup the discontinuous water height and velocity
-  if element_id == 10
-    H = 5.0
-    v = 0.1882
-  end
-
-  # Setup a discontinuous bottom topography using the element id number
-  if element_id == 7
-    b = 2.0 + 0.5 * sin(2.0 * pi * x[1])
-  end
-
-  return prim2cons(SVector(H, v, b), equations)
-end
-
-# point to the data we want to augment
-u = Trixi.wrap_array(ode.u0, semi)
-# reset the initial condition
-for element in eachelement(semi.solver, semi.cache)
-  for i in eachnode(semi.solver)
-    x_node = Trixi.get_node_coords(semi.cache.elements.node_coordinates, equations, semi.solver, i, element)
-    u_node = initial_condition_ec_discontinuous_bottom(x_node, first(tspan), element, equations)
-    Trixi.set_node_vars!(u, u_node, equations, semi.solver, i, element)
-  end
-end
-
 ###############################################################################
 # Callbacks
 
diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_shock_capturing.jl b/examples/tree_1d_dgsem/elixir_shallowwater_shock_capturing.jl
index 50241126a28..62346d7b5ab 100644
--- a/examples/tree_1d_dgsem/elixir_shallowwater_shock_capturing.jl
+++ b/examples/tree_1d_dgsem/elixir_shallowwater_shock_capturing.jl
@@ -7,24 +7,37 @@ using Trixi
 
 equations = ShallowWaterEquations1D(gravity_constant=9.812, H0=1.75)
 
-function initial_condition_stone_throw(x, t, equations::ShallowWaterEquations1D)
-    # Set up polar coordinates
-    inicenter = 0.15
-    x_norm = x[1] - inicenter[1]
-    r = abs(x_norm)
+# Initial condition with a truly discontinuous velocity and bottom topography.
+# Works as intended for TreeMesh1D with `initial_refinement_level=3`. If the mesh
+# refinement level is changed the initial condition below may need changed as well to
+# ensure that the discontinuities lie on an element interface.
+function initial_condition_stone_throw_discontinuous_bottom(x, t, equations::ShallowWaterEquations1D)
+
+  # Calculate primitive variables
+
+  # flat lake
+  H = equations.H0
+
+  # Discontinuous velocity
+  v = 0.0
+  if x[1] >= -0.75 && x[1] <= 0.0
+      v = -1.0
+  elseif x[1] >= 0.0 && x[1] <= 0.75
+      v = 1.0
+  end
 
-    # Calculate primitive variables
-    H = equations.H0
-    # v = 0.0 # for well-balanced test
-    v = r < 0.6 ? 1.75 : 0.0 # for stone throw
+  b = (  1.5 / exp( 0.5 * ((x[1] - 1.0)^2 ) )
+     + 0.75 / exp( 0.5 * ((x[1] + 1.0)^2 ) ) )
 
-    b = (  1.5 / exp( 0.5 * ((x[1] - 1.0)^2 ) )
-       + 0.75 / exp( 0.5 * ((x[1] + 1.0)^2 ) ) )
+  # Force a discontinuous bottom topography
+  if x[1] >= -1.5 && x[1] <= 0.0
+    b = 0.5
+  end
 
-    return prim2cons(SVector(H, v, b), equations)
+  return prim2cons(SVector(H, v, b), equations)
 end
 
-initial_condition = initial_condition_stone_throw
+initial_condition = initial_condition_stone_throw_discontinuous_bottom
 
 boundary_condition = boundary_condition_slip_wall
 
@@ -62,49 +75,13 @@ semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
                                     boundary_conditions = boundary_condition)
 
 ###############################################################################
-# ODE solvers, callbacks etc.
+# ODE solver
 
 tspan = (0.0, 3.0)
 ode = semidiscretize(semi, tspan)
 
-# Hack in a discontinuous bottom for a more interesting test
-function initial_condition_stone_throw_discontinuous_bottom(x, t, element_id, equations::ShallowWaterEquations1D)
-
-    inicenter = 0.15
-    x_norm = x[1] - inicenter[1]
-    r = abs(x_norm)
-
-    # Calculate primitive variables
-    H = equations.H0 # flat lake
-    # Discontinuous velocity set via element id number
-    v = 0.0
-    if element_id == 4
-        v = -1.0
-    elseif element_id == 5
-        v = 1.0
-    end
-
-    b = (  1.5 / exp( 0.5 * ((x[1] - 1.0)^2 ) )
-       + 0.75 / exp( 0.5 * ((x[1] + 1.0)^2 ) ) )
-
-    # Setup a discontinuous bottom topography using the element id number
-    if element_id == 3 || element_id == 4
-      b = 0.5
-    end
-
-    return prim2cons(SVector(H, v, b), equations)
-end
-
-# point to the data we want to augment
-u = Trixi.wrap_array(ode.u0, semi)
-# reset the initial condition
-for element in eachelement(semi.solver, semi.cache)
-  for i in eachnode(semi.solver)
-    x_node = Trixi.get_node_coords(semi.cache.elements.node_coordinates, equations, semi.solver, i, element)
-    u_node = initial_condition_stone_throw_discontinuous_bottom(x_node, first(tspan), element, equations)
-    Trixi.set_node_vars!(u, u_node, equations, semi.solver, i, element)
-  end
-end
+###############################################################################
+# Callbacks
 
 summary_callback = SummaryCallback()
 
diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_twolayer_dam_break.jl b/examples/tree_1d_dgsem/elixir_shallowwater_twolayer_dam_break.jl
index 67c1098b178..60770d158fa 100644
--- a/examples/tree_1d_dgsem/elixir_shallowwater_twolayer_dam_break.jl
+++ b/examples/tree_1d_dgsem/elixir_shallowwater_twolayer_dam_break.jl
@@ -3,20 +3,34 @@ using OrdinaryDiffEq
 using Trixi
 
 ###############################################################################
-# Semidiscretization of the two-layer shallow water equations for a dam break test with a 
-# discontinuous bottom topography function to test entropy conservation
+# Semidiscretization of the two-layer shallow water equations for a dam break
+# test with a discontinuous bottom topography function to test entropy conservation
 
 equations = ShallowWaterTwoLayerEquations1D(gravity_constant=9.81, H0=2.0, rho_upper=0.9, rho_lower=1.0)
 
-###############################################################################
-# Workaround to set a discontinuous bottom topography and initial condition.
+# Initial condition of a dam break with a discontinuous water heights and bottom topography.
+# Works as intended for TreeMesh1D with `initial_refinement_level=5`. If the mesh
+# refinement level is changed the initial condition below may need changed as well to
+# ensure that the discontinuities lie on an element interface.
+function initial_condition_dam_break(x, t, equations::ShallowWaterTwoLayerEquations1D)
+  v1_upper = 0.0
+  v1_lower = 0.0
+
+  # Set the discontinuity
+  if x[1] <= 10.0
+    H_lower = 2.0
+    H_upper = 4.0
+    b  = 0.0
+  else
+    H_lower = 1.5
+    H_upper = 3.0
+    b  = 0.5
+  end
+
+  return prim2cons(SVector(H_upper, v1_upper, H_lower, v1_lower, b), equations)
+end
 
-# This test case uses a special work around to setup a truly discontinuous bottom topography 
-# function and initial condition for this academic testcase of entropy conservation. First, a 
-# dummy initial condition is introduced to create the semidiscretization. Then the initial condition 
-# is reset with the true discontinuous values from initial_condition_dam_break. Note, that this
-# initial condition only works for TreeMesh1D with `initial_refinement_level=5`.
-initial_condition = initial_condition_convergence_test
+initial_condition = initial_condition_dam_break
 
 ###############################################################################
 # Get the DG approximation space
@@ -25,7 +39,6 @@ volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal)
 solver = DGSEM(polydeg=3, surface_flux=(flux_fjordholm_etal, flux_nonconservative_fjordholm_etal),
                volume_integral=VolumeIntegralFluxDifferencing(volume_flux))
 
-
 ###############################################################################
 # Get the TreeMesh and setup a non-periodic mesh
 
@@ -39,54 +52,22 @@ mesh = TreeMesh(coordinates_min, coordinates_max,
 boundary_condition = boundary_condition_slip_wall
 
 # create the semidiscretization object
-semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, 
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
                                     boundary_conditions=boundary_condition)
 
 ###############################################################################
-# ODE solvers, callbacks etc.
+# ODE solvers
 
 tspan = (0.0,0.4)
 ode = semidiscretize(semi, tspan)
 
-# Initial conditions dam break test case
-function initial_condition_dam_break(x, t, element_id, equations::ShallowWaterTwoLayerEquations1D)
-  v1_upper = 0.0
-  v1_lower = 0.0
-
-  # Set the discontinuity
-  if element_id <= 16
-    H_lower = 2.0
-    H_upper = 4.0
-    b  = 0.0
-  else
-    H_lower = 1.5
-    H_upper = 3.0
-    b  = 0.5
-  end
-
-  return prim2cons(SVector(H_upper, v1_upper, H_lower, v1_lower, b), equations)
-end
-
-
-# point to the data we want to augment
-u = Trixi.wrap_array(ode.u0, semi)
-# reset the initial condition
-for element in eachelement(semi.solver, semi.cache)
-  for i in eachnode(semi.solver)
-    x_node = Trixi.get_node_coords(semi.cache.elements.node_coordinates,
-      equations, semi.solver, i, element)
-    u_node = initial_condition_dam_break(x_node, first(tspan), element, equations)
-    Trixi.set_node_vars!(u, u_node, equations, semi.solver, i, element)
-  end
-end
-
-
-
+###############################################################################
+# Callbacks
 
 summary_callback = SummaryCallback()
 
 analysis_interval = 500
-analysis_callback = AnalysisCallback(semi, interval=analysis_interval, save_analysis=false, 
+analysis_callback = AnalysisCallback(semi, interval=analysis_interval, save_analysis=false,
     extra_analysis_integrals=(energy_total, energy_kinetic, energy_internal,))
 
 stepsize_callback = StepsizeCallback(cfl=1.0)
diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced.jl b/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced.jl
index 83835656839..e07bc04d76a 100644
--- a/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced.jl
+++ b/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced.jl
@@ -8,21 +8,28 @@ using Trixi
 
 equations = ShallowWaterEquations1D(gravity_constant=9.81, H0=3.25)
 
-# An initial condition with constant total water height and zero velocities to test well-balancedness.
-# Note, this routine is used to compute errors in the analysis callback but the initialization is
-# overwritten by `initial_condition_discontinuous_well_balancedness` below.
-function initial_condition_well_balancedness(x, t, equations::ShallowWaterEquations1D)
+# Setup a truly discontinuous bottom topography function for this academic
+# testcase of well-balancedness. The errors from the analysis callback are
+# not important but the error for this lake-at-rest test case
+# `∑|H0-(h+b)|` should be around machine roundoff.
+# Works as intended for TreeMesh1D with `initial_refinement_level=3`. If the mesh
+# refinement level is changed the initial condition below may need changed as well to
+# ensure that the discontinuities lie on an element interface.
+function initial_condition_discontinuous_well_balancedness(x, t, equations::ShallowWaterEquations1D)
   # Set the background values
   H = equations.H0
   v = 0.0
+  b = 0.0
 
-  # bottom topography inspired by from Pond.control in [HOHQMesh](https://github.com/trixi-framework/HOHQMesh)
-  b = (1.5 / exp( 0.5 * ((x[1] - 1.0)^2) )+ 0.75 / exp(0.5 * ((x[1] + 1.0)^2)))
+  # Setup a discontinuous bottom topography
+  if x[1] >= 0.5 && x[1] <= 0.75
+    b = 2.0 + 0.5 * sin(2.0 * pi * x[1])
+  end
 
   return prim2cons(SVector(H, v, b), equations)
 end
 
-initial_condition = initial_condition_well_balancedness
+initial_condition = initial_condition_discontinuous_well_balancedness
 
 ###############################################################################
 # Get the DG approximation space
@@ -50,41 +57,6 @@ semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
 tspan = (0.0, 100.0)
 ode = semidiscretize(semi, tspan)
 
-###############################################################################
-# Workaround to set a discontinuous bottom topography and initial condition for debugging and testing.
-
-# alternative version of the initial conditinon used to setup a truly discontinuous
-# bottom topography function for this academic testcase of well-balancedness.
-# The errors from the analysis callback are not important but the error for this lake at rest test case
-# `∑|H0-(h+b)|` should be around machine roundoff.
-# In contrast to the usual signature of initial conditions, this one get passed the
-# `element_id` explicitly. In particular, this initial conditions works as intended
-# only for the TreeMesh1D with `initial_refinement_level=3`.
-function initial_condition_discontinuous_well_balancedness(x, t, element_id, equations::ShallowWaterEquations1D)
-  # Set the background values
-  H = equations.H0
-  v = 0.0
-  b = 0.0
-
-  # Setup a discontinuous bottom topography using the element id number
-  if element_id == 7
-    b = 2.0 + 0.5 * sin(2.0 * pi * x[1])
-  end
-
-  return prim2cons(SVector(H, v, b), equations)
-end
-
-# point to the data we want to augment
-u = Trixi.wrap_array(ode.u0, semi)
-# reset the initial condition
-for element in eachelement(semi.solver, semi.cache)
-  for i in eachnode(semi.solver)
-    x_node = Trixi.get_node_coords(semi.cache.elements.node_coordinates, equations, semi.solver, i, element)
-    u_node = initial_condition_discontinuous_well_balancedness(x_node, first(tspan), element, equations)
-    Trixi.set_node_vars!(u, u_node, equations, semi.solver, i, element)
-  end
-end
-
 ###############################################################################
 # Callbacks
 
diff --git a/src/equations/shallow_water_two_layer_1d.jl b/src/equations/shallow_water_two_layer_1d.jl
index edf7d5e32ff..02899171509 100644
--- a/src/equations/shallow_water_two_layer_1d.jl
+++ b/src/equations/shallow_water_two_layer_1d.jl
@@ -11,28 +11,28 @@
 Two-Layer Shallow Water equations (2LSWE) in one space dimension. The equations are given by
 ```math
 \begin{alignat*}{4}
-&\frac{\partial}{\partial t}h_{upper} 
-&&+ \frac{\partial}{\partial x}\left(h_{upper} v_{1,upper}\right) 
+&\frac{\partial}{\partial t}h_{upper}
+&&+ \frac{\partial}{\partial x}\left(h_{upper} v_{1,upper}\right)
 &&= 0 \\
-&\frac{\partial}{\partial t}\left(h_{upper}v_{1,upper}\right) 
-&&+ \frac{\partial}{\partial x}\left(h_{upper}v_{1,upper}^2 + \dfrac{gh_{upper}^2}{2}\right) 
+&\frac{\partial}{\partial t}\left(h_{upper}v_{1,upper}\right)
+&&+ \frac{\partial}{\partial x}\left(h_{upper}v_{1,upper}^2 + \dfrac{gh_{upper}^2}{2}\right)
 &&= -gh_{upper}\frac{\partial}{\partial x}\left(b+h_{lower}\right)\\
-&\frac{\partial}{\partial t}h_{lower}  
-&&+ \frac{\partial}{\partial x}\left(h_{lower}v_{1,lower}\right) 
+&\frac{\partial}{\partial t}h_{lower}
+&&+ \frac{\partial}{\partial x}\left(h_{lower}v_{1,lower}\right)
 &&= 0 \\
-&\frac{\partial}{\partial t}\left(h_{lower}v_{1,lower}\right)  
-&&+ \frac{\partial}{\partial x}\left(h_{lower}v_{1,lower}^2 + \dfrac{gh_{lower}^2}{2}\right) 
+&\frac{\partial}{\partial t}\left(h_{lower}v_{1,lower}\right)
+&&+ \frac{\partial}{\partial x}\left(h_{lower}v_{1,lower}^2 + \dfrac{gh_{lower}^2}{2}\right)
 &&= -gh_{lower}\frac{\partial}{\partial x}\left(b+\dfrac{\rho_{upper}}{\rho_{lower}}h_{upper}\right).
 \end{alignat*}
 ```
-The unknown quantities of the 2LSWE are the water heights of the {lower} layer ``h_{lower}`` and the 
-{upper} layer ``h_{upper}`` with respective velocities ``v_{1,upper}`` and ``v_{1,lower}``. The gravitational constant is 
-denoted by `g`, the layer densitites by ``\rho_{upper}``and ``\rho_{lower}`` and the (possibly) variable 
-bottom topography function ``b(x)``. The conservative variable water height ``h_{lower}`` is measured 
-from the bottom topography ``b`` and ``h_{upper}`` relative to ``h_{lower}``, therefore one also defines the 
+The unknown quantities of the 2LSWE are the water heights of the {lower} layer ``h_{lower}`` and the
+{upper} layer ``h_{upper}`` with respective velocities ``v_{1,upper}`` and ``v_{1,lower}``. The gravitational constant is
+denoted by `g`, the layer densitites by ``\rho_{upper}``and ``\rho_{lower}`` and the (possibly) variable
+bottom topography function ``b(x)``. The conservative variable water height ``h_{lower}`` is measured
+from the bottom topography ``b`` and ``h_{upper}`` relative to ``h_{lower}``, therefore one also defines the
 total water heights as ``H_{upper} = h_{upper} + h_{upper} + b`` and ``H_{lower} = h_{lower} + b``.
 
-The densities must be chosen such that ``\rho_{upper} < \rho_{lower}``, to make sure that the heavier fluid 
+The densities must be chosen such that ``\rho_{upper} < \rho_{lower}``, to make sure that the heavier fluid
 ``\rho_{lower}`` is in the bottom layer and the lighter fluid ``\rho_{upper}`` in the {upper} layer.
 
 The additional quantity ``H_0`` is also available to store a reference value for the total water
@@ -41,13 +41,13 @@ height that is useful to set initial conditions or test the "lake-at-rest" well-
 The bottom topography function ``b(x)`` is set inside the initial condition routine
 for a particular problem setup.
 
-In addition to the unknowns, Trixi currently stores the bottom topography values at the 
-approximation points despite being fixed in time. This is done for convenience of computing the 
-bottom topography gradients on the fly during the approximation as well as computing auxiliary 
+In addition to the unknowns, Trixi currently stores the bottom topography values at the
+approximation points despite being fixed in time. This is done for convenience of computing the
+bottom topography gradients on the fly during the approximation as well as computing auxiliary
 quantities like the total water height ``H`` or the entropy variables.
 This affects the implementation and use of these equations in various ways:
 * The flux values corresponding to the bottom topography must be zero.
-* The bottom topography values must be included when defining initial conditions, boundary 
+* The bottom topography values must be included when defining initial conditions, boundary
   conditions or source terms.
 * [`AnalysisCallback`](@ref) analyzes this variable.
 * Trixi's visualization tools will visualize the bottom topography by default.
@@ -101,7 +101,7 @@ end
     initial_condition_convergence_test(x, t, equations::ShallowWaterTwoLayerEquations1D)
 
 A smooth initial condition used for convergence tests in combination with
-[`source_terms_convergence_test`](@ref) (and 
+[`source_terms_convergence_test`](@ref) (and
 [`BoundaryConditionDirichlet(initial_condition_convergence_test)`](@ref) in non-periodic domains).
 """
 function initial_condition_convergence_test(x, t,
@@ -121,9 +121,9 @@ end
 """
     source_terms_convergence_test(u, x, t, equations::ShallowWaterTwoLayerEquations1D)
 
-Source terms used for convergence tests in combination with 
-[`initial_condition_convergence_test`](@ref) 
-(and [`BoundaryConditionDirichlet(initial_condition_convergence_test)`](@ref) 
+Source terms used for convergence tests in combination with
+[`initial_condition_convergence_test`](@ref)
+(and [`BoundaryConditionDirichlet(initial_condition_convergence_test)`](@ref)
 in non-periodic domains).
 """
 @inline function source_terms_convergence_test(u, x, t,
@@ -167,8 +167,8 @@ the internal value.
 
 For details see Section 9.2.5 of the book:
 - Eleuterio F. Toro (2001)
-  Shock-Capturing Methods for Free-Surface Shallow Flows  
-  1st edition  
+  Shock-Capturing Methods for Free-Surface Shallow Flows
+  1st edition
   ISBN 0471987662
 """
 @inline function boundary_condition_slip_wall(u_inner, orientation_or_normal, direction,
@@ -219,7 +219,7 @@ end
 
 Non-symmetric two-point volume flux discretizing the nonconservative (source) term
 that contains the gradient of the bottom topography [`ShallowWaterTwoLayerEquations2D`](@ref) and an
-additional term that couples the momentum of both layers. This is a slightly modified version 
+additional term that couples the momentum of both layers. This is a slightly modified version
 to account for the additional source term compared to the standard SWE described in the paper.
 
 Further details are available in the paper:
@@ -238,7 +238,7 @@ Further details are available in the paper:
 
     z = zero(eltype(u_ll))
 
-    # Bottom gradient nonconservative term: (0, g*h_upper*(b+h_lower)_x, 
+    # Bottom gradient nonconservative term: (0, g*h_upper*(b+h_lower)_x,
     #                                        0, g*h_lower*(b+r*h_upper)_x, 0)
     f = SVector(z,
                 equations.gravity * h_upper_ll * (b_rr + h_lower_rr),
@@ -254,9 +254,9 @@ end
 
 !!! warning "Experimental code"
     This numerical flux is experimental and may change in any future release.
-                                    
-Non-symmetric two-point surface flux discretizing the nonconservative (source) term that contains 
-the gradients of the bottom topography and an additional term that couples the momentum of both 
+
+Non-symmetric two-point surface flux discretizing the nonconservative (source) term that contains
+the gradients of the bottom topography and an additional term that couples the momentum of both
 layers [`ShallowWaterTwoLayerEquations2D`](@ref).
 
 Further details are available in the paper:
@@ -286,7 +286,7 @@ formulation.
 
     z = zero(eltype(u_ll))
 
-    # Bottom gradient nonconservative term: (0, g*h_upper*(b+h_lower)_x, 
+    # Bottom gradient nonconservative term: (0, g*h_upper*(b+h_lower)_x,
     #                                        0, g*h_lower*(b+r*h_upper)_x, 0)
     f = SVector(z,
                 g * h_upper_ll * (b_ll + h_lower_ll) +
@@ -303,13 +303,13 @@ end
     flux_fjordholm_etal(u_ll, u_rr, orientation,
                         equations::ShallowWaterTwoLayerEquations1D)
 
-Total energy conservative (mathematical entropy for shallow water equations). When the bottom 
-topography is nonzero this should only be used as a surface flux otherwise the scheme will not be 
+Total energy conservative (mathematical entropy for shallow water equations). When the bottom
+topography is nonzero this should only be used as a surface flux otherwise the scheme will not be
 well-balanced. For well-balancedness in the volume flux use [`flux_wintermeyer_etal`](@ref).
 
 Details are available in Eq. (4.1) in the paper:
 - Ulrik S. Fjordholm, Siddhartha Mishra and Eitan Tadmor (2011)
-  Well-balanced and energy stable schemes for the shallow water equations with discontinuous 
+  Well-balanced and energy stable schemes for the shallow water equations with discontinuous
   topography [DOI: 10.1016/j.jcp.2011.03.042](https://doi.org/10.1016/j.jcp.2011.03.042)
 and the application to two layers is shown in the paper:
 - Ulrik Skre Fjordholm (2012)
@@ -348,11 +348,11 @@ end
 """
     flux_wintermeyer_etal(u_ll, u_rr, orientation,
                           equations::ShallowWaterTwoLayerEquations1D)
-                      
+
 Total energy conservative (mathematical entropy for two-layer 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). To obtain the flux for the
-two-layer shallow water equations the flux that is described in the paper for the normal shallow 
+two-layer shallow water equations the flux that is described in the paper for the normal shallow
 water equations is used within each layer.
 
 Further details are available in Theorem 1 of the paper:
@@ -391,8 +391,8 @@ end
     flux_es_fjordholm_etal(u_ll, u_rr, orientation,
                            equations::ShallowWaterTwoLayerEquations1D)
 
-Entropy stable surface flux for the two-layer shallow water equations. Uses the entropy 
-conservative flux_fjordholm_etal and adds a Lax-Friedrichs type dissipation dependent on the jump 
+Entropy stable surface flux for the two-layer shallow water equations. Uses the entropy
+conservative flux_fjordholm_etal and adds a Lax-Friedrichs type dissipation dependent on the jump
 of entropy variables.
 
 Further details are available in the paper:
@@ -460,10 +460,10 @@ formulation.
 end
 
 # Calculate approximation for maximum wave speed for local Lax-Friedrichs-type dissipation as the
-# maximum velocity magnitude plus the maximum speed of sound. This function uses approximate 
-# eigenvalues using the speed of the barotropic mode as there is no simple way to calculate them 
+# maximum velocity magnitude plus the maximum speed of sound. This function uses approximate
+# eigenvalues using the speed of the barotropic mode as there is no simple way to calculate them
 # analytically.
-# 
+#
 # A good overview of the derivation is given in:
 # -  Jonas Nycander, Andrew McC. Hogg, Leela M. Frankcombe (2008)
 #    Open boundary conditions for nonlinear channel Flows
@@ -488,7 +488,7 @@ end
     return (max(abs(v_m_ll) + c_ll, abs(v_m_rr) + c_rr))
 end
 
-# Specialized `DissipationLocalLaxFriedrichs` to avoid spurious dissipation in the bottom 
+# Specialized `DissipationLocalLaxFriedrichs` to avoid spurious dissipation in the bottom
 # topography
 @inline function (dissipation::DissipationLocalLaxFriedrichs)(u_ll, u_rr,
                                                               orientation_or_normal_direction,
@@ -530,7 +530,7 @@ end
 end
 
 # Convert conservative variables to entropy variables
-# Note, only the first four are the entropy variables, the fifth entry still just carries the 
+# Note, only the first four are the entropy variables, the fifth entry still just carries the
 # bottom topography values for convenience
 @inline function cons2entropy(u, equations::ShallowWaterTwoLayerEquations1D)
     h_upper, _, h_lower, _, b = u
@@ -567,7 +567,7 @@ end
 
 # Calculate total energy for a conservative state `cons`
 @inline function energy_total(cons, equations::ShallowWaterTwoLayerEquations1D)
-    h_upper, h_lower, h_v1_upper, h_v2_lower, b = cons
+    h_upper, h_v1_upper, h_lower, h_v2_lower, b = cons
     # Set new variables for better readability
     g = equations.gravity
     rho_upper = equations.rho_upper
diff --git a/src/solvers/dg.jl b/src/solvers/dg.jl
index 838fa2d5819..2536cfe0bf2 100644
--- a/src/solvers/dg.jl
+++ b/src/solvers/dg.jl
@@ -628,6 +628,15 @@ function compute_coefficients!(u, func, t, mesh::AbstractMesh{1}, equations, dg:
         for i in eachnode(dg)
             x_node = get_node_coords(cache.elements.node_coordinates, equations, dg, i,
                                      element)
+            # Changing the node positions passed to the initial condition by the minimum
+            # amount possible with the current type of floating point numbers allows setting
+            # discontinuous initial data in a simple way. In particular, a check like `if x < x_jump`
+            # works if the jump location `x_jump` is at the position of an interface.
+            if i == 1
+                x_node = SVector(nextfloat(x_node[1]))
+            elseif i == nnodes(dg)
+                x_node = SVector(prevfloat(x_node[1]))
+            end
             u_node = func(x_node, t, equations)
             set_node_vars!(u, u_node, equations, dg, i, element)
         end
diff --git a/test/test_structured_1d.jl b/test/test_structured_1d.jl
index a27d3c219e1..ec8c7a138d5 100644
--- a/test/test_structured_1d.jl
+++ b/test/test_structured_1d.jl
@@ -27,8 +27,8 @@ isdir(outdir) && rm(outdir, recursive=true)
 
   @trixi_testset "elixir_advection_shockcapturing.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_shockcapturing.jl"),
-      l2   = [0.08004076716881656],
-      linf = [0.6342577638501385],
+      l2   = [0.08015029105233593],
+      linf = [0.610709468736576],
       atol = 1.0e-5)
   end
 
diff --git a/test/test_tree_1d.jl b/test/test_tree_1d.jl
index e37e1efc3e6..7737a93a15a 100644
--- a/test/test_tree_1d.jl
+++ b/test/test_tree_1d.jl
@@ -271,7 +271,7 @@ end
   sol = solve(ode, Tsit5(), abstol=1.0e-6, reltol=1.0e-6,
               save_everystep=false, callback=callbacks);
 
-  @test analysis_callback(sol).l2 ≈ [0.00029610274971929974, 5.573684084938363e-6]
+  @test analysis_callback(sol).l2 ≈ [0.00029609575838969394, 5.5681704039507985e-6]
 end
 
 
diff --git a/test/test_tree_1d_burgers.jl b/test/test_tree_1d_burgers.jl
index 788c7ab4199..8c4cfaa406d 100644
--- a/test/test_tree_1d_burgers.jl
+++ b/test/test_tree_1d_burgers.jl
@@ -22,14 +22,14 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem")
 
   @trixi_testset "elixir_burgers_shock.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_burgers_shock.jl"),
-      l2   = [0.4407585104869119],
-      linf = [1.000000000000001])
+      l2   = [0.4422505602587537],
+      linf = [1.0000000000000009])
   end
 
   @trixi_testset "elixir_burgers_rarefaction.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_burgers_rarefaction.jl"),
-      l2   = [0.40287062735307044],
-      linf = [1.0042992585765542])
+      l2   = [0.4038224690923722],
+      linf = [1.0049201454652736])
   end
 end
 
diff --git a/test/test_tree_1d_euler.jl b/test/test_tree_1d_euler.jl
index 40f2a38b0e1..5fb74b80bce 100644
--- a/test/test_tree_1d_euler.jl
+++ b/test/test_tree_1d_euler.jl
@@ -22,8 +22,8 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem")
 
   @trixi_testset "elixir_euler_density_wave.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_density_wave.jl"),
-      l2   = [0.0011482554820185795, 0.00011482554830363504, 5.741277417754598e-6],
-      linf = [0.004090978306820037, 0.00040909783134346345, 2.0454891732413216e-5])
+      l2   = [0.0011482554820217855, 0.00011482554830323462, 5.741277429325267e-6],
+      linf = [0.004090978306812376, 0.0004090978313582294, 2.045489210189544e-5])
   end
 
   @trixi_testset "elixir_euler_density_wave.jl with initial_condition_constant" begin
@@ -41,14 +41,14 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem")
 
   @trixi_testset "elixir_euler_ec.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_ec.jl"),
-      l2   = [0.11915540925414216, 0.15489191247295198, 0.44543052524765375],
-      linf = [0.2751485868543495, 0.2712764982000735, 0.9951407418216425])
+      l2   = [0.11821957357197649, 0.15330089521538678, 0.4417674632047301],
+      linf = [0.24280567569982958, 0.29130548795961936, 0.8847009003152442])
   end
 
   @trixi_testset "elixir_euler_ec.jl with flux_kennedy_gruber" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_ec.jl"),
-      l2   = [0.07905582221868049, 0.10180958900546237, 0.29596551476711125],
-      linf = [0.23515297345769826, 0.2958208108392532, 0.8694224308790321],
+      l2   = [0.07803455838661963, 0.10032577312032283, 0.29228156303827935],
+      linf = [0.2549869853794955, 0.3376472164661263, 0.9650477546553962],
       maxiters = 10,
       surface_flux = flux_kennedy_gruber,
       volume_flux = flux_kennedy_gruber)
@@ -56,8 +56,8 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem")
 
   @trixi_testset "elixir_euler_ec.jl with flux_shima_etal" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_ec.jl"),
-      l2   = [0.07909267609417114, 0.1018246500951966, 0.2959649187481973],
-      linf = [0.23631829743146504, 0.2977756307879202, 0.8642794698697331],
+      l2   = [0.07800654460172655, 0.10030365573277883, 0.2921481199111959],
+      linf = [0.25408579350400395, 0.3388657679031271, 0.9776486386921928],
       maxiters = 10,
       surface_flux = flux_shima_etal,
       volume_flux = flux_shima_etal)
@@ -65,8 +65,8 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem")
 
   @trixi_testset "elixir_euler_ec.jl with flux_chandrashekar" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_ec.jl"),
-      l2   = [0.07905306555214126, 0.10181180378499956, 0.2959171937479504],
-      linf = [0.24057642004451651, 0.29691454643616433, 0.886425723870524],
+      l2   = [0.07801923089205756, 0.10039557434912669, 0.2922210399923278],
+      linf = [0.2576521982607225, 0.3409717926625057, 0.9772961936567048],
       maxiters = 10,
       surface_flux = flux_chandrashekar,
       volume_flux = flux_chandrashekar)
@@ -74,8 +74,8 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem")
 
   @trixi_testset "elixir_euler_ec.jl with flux_hll" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_ec.jl"),
-      l2   = [0.07959780803600519, 0.10342491934977621, 0.2978851659149904],
-      linf = [0.19228754121840885, 0.2524152253292552, 0.725604944702432],
+      l2   = [0.07852272782240548, 0.10209790867523805, 0.293873048809011],
+      linf = [0.19244768908604093, 0.2515941686151897, 0.7258000837553769],
       maxiters = 10,
       surface_flux = flux_hll,
       volume_flux = flux_ranocha)
@@ -83,8 +83,8 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem")
 
   @trixi_testset "elixir_euler_shockcapturing.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_shockcapturing.jl"),
-      l2   = [0.11665968950973675, 0.15105507394693413, 0.43503082674771115],
-      linf = [0.1867400345208743, 0.24621854448555328, 0.703826406555577])
+      l2   = [0.11606096465319675, 0.15028768943458806, 0.4328230323046703],
+      linf = [0.18031710091067965, 0.2351582421501841, 0.6776805692092567])
   end
 
   @trixi_testset "elixir_euler_sedov_blast_wave.jl" begin
@@ -96,8 +96,8 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem")
 
   @trixi_testset "elixir_euler_sedov_blast_wave_pure_fv.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_sedov_blast_wave_pure_fv.jl"),
-      l2   = [1.075075094036344, 0.06766902169711514, 0.9221426570128292],
-      linf = [3.3941512671408542, 0.16862631133303882, 2.6572394126490315],
+      l2   = [1.0735456065491455, 0.07131078703089379, 0.9205739468590453],
+      linf = [3.4296365168219216, 0.17635583964559245, 2.6574584326179505],
       # Let this test run longer to cover some lines in flux_hllc
       coverage_override = (maxiters=10^5, tspan=(0.0, 0.1)))
   end
@@ -129,22 +129,22 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem")
 
   @trixi_testset "elixir_euler_blast_wave.jl" begin
   @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_blast_wave.jl"),
-      l2   = [0.21651329948737183, 0.28091709900008616, 0.5580778880050432],
-      linf = [1.513525457073142, 1.5328754303137992, 2.0467706106669556],
+      l2   = [0.21934822867340323, 0.28131919126002686, 0.554361702716662],
+      linf = [1.5180897390290355, 1.3967085956620369, 2.0663825294019595],
       maxiters = 30)
   end
 
   @trixi_testset "elixir_euler_blast_wave_neuralnetwork_perssonperaire.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_blast_wave_neuralnetwork_perssonperaire.jl"),
-        l2   = [2.13605618e-01, 2.79953055e-01, 5.54424459e-01],
-        linf = [1.55151701e+00, 1.55696782e+00, 2.05525953e+00],
+        l2   = [0.21814833203212694, 0.2818328665444332, 0.5528379124720818],
+        linf = [1.5548653877320868, 1.4474018998129738, 2.071919577393772],
         maxiters = 30)
   end
 
   @trixi_testset "elixir_euler_blast_wave_neuralnetwork_rayhesthaven.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_blast_wave_neuralnetwork_rayhesthaven.jl"),
-        l2   = [2.18148857e-01, 2.83182959e-01, 5.59096194e-01],
-        linf = [1.62706876e+00, 1.61680275e+00, 2.05876517e+00],
+        l2   = [0.22054468879127423, 0.2828269190680846, 0.5542369885642424],
+        linf = [1.5623359741479623, 1.4290121654488288, 2.1040405133123072],
         maxiters = 30)
   end
 end
diff --git a/test/test_tree_1d_eulermulti.jl b/test/test_tree_1d_eulermulti.jl
index eac54a9372c..e880f98e2d0 100644
--- a/test/test_tree_1d_eulermulti.jl
+++ b/test/test_tree_1d_eulermulti.jl
@@ -11,39 +11,47 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem")
 
   @trixi_testset "elixir_eulermulti_ec.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_ec.jl"),
-      l2   = [1.54891912e-01, 4.45430525e-01, 1.70222013e-02, 3.40444026e-02, 6.80888053e-02],
-      linf = [2.71276498e-01, 9.95140742e-01, 3.93069410e-02, 7.86138820e-02, 1.57227764e-01])
+      l2   = [0.15330089521538684, 0.4417674632047301, 0.016888510510282385, 0.03377702102056477,
+              0.06755404204112954],
+      linf = [0.29130548795961864, 0.8847009003152357, 0.034686525099975274, 0.06937305019995055,
+              0.1387461003999011])
   end
 
   @trixi_testset "elixir_eulermulti_es.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_es.jl"),
-      l2   = [1.53387916e-01, 4.41585576e-01, 3.93605635e-02, 7.87211270e-02],
-      linf = [2.49632117e-01, 7.21088064e-01, 6.38328770e-02, 1.27665754e-01])
+      l2   = [0.1522380497572071, 0.43830846465313206, 0.03907262116499431, 0.07814524232998862],
+      linf = [0.24939193075537294, 0.7139395740052739, 0.06324208768391237, 0.12648417536782475])
   end
 
   @trixi_testset "elixir_eulermulti_convergence_ec.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_convergence_ec.jl"),
-      l2   = [8.57523604e-05, 1.63878043e-04, 1.94126993e-05, 3.88253986e-05],
-      linf = [3.05932773e-04, 6.24480393e-04, 7.25312144e-05, 1.45062429e-04])
+      l2   = [8.575236038539227e-5, 0.00016387804318585358, 1.9412699303977585e-5, 3.882539860795517e-5],
+      linf = [0.00030593277277124464, 0.0006244803933350696, 7.253121435135679e-5, 0.00014506242870271358])
   end
 
   @trixi_testset "elixir_eulermulti_convergence_es.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_convergence_es.jl"),
-      l2   = [1.8983933794407234e-5, 6.207744299844731e-5, 1.5466205761868047e-6, 3.0932411523736094e-6, 6.186482304747219e-6, 1.2372964609494437e-5],
-      linf = [0.00012014372605895218, 0.0003313207215800418, 6.50836791016296e-6, 1.301673582032592e-5, 2.603347164065184e-5, 5.206694328130368e-5])
+      l2   = [1.8983933794407234e-5, 6.207744299844731e-5, 1.5466205761868047e-6, 3.0932411523736094e-6,
+              6.186482304747219e-6, 1.2372964609494437e-5],
+      linf = [0.00012014372605895218, 0.0003313207215800418, 6.50836791016296e-6, 1.301673582032592e-5,
+              2.603347164065184e-5, 5.206694328130368e-5])
   end
 
   @trixi_testset "elixir_eulermulti_convergence_es.jl with flux_chandrashekar" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_convergence_es.jl"),
-      l2   = [1.88845048e-05, 5.49106005e-05, 9.42673716e-07, 1.88534743e-06, 3.77069486e-06, 7.54138973e-06],
-      linf = [1.16223512e-04, 3.07922197e-04, 3.21774233e-06, 6.43548465e-06, 1.28709693e-05, 2.57419386e-05],
+      l2   = [1.888450477353845e-5, 5.4910600482795386e-5, 9.426737161533622e-7, 1.8853474323067245e-6,
+              3.770694864613449e-6, 7.541389729226898e-6],
+      linf = [0.00011622351152063004, 0.0003079221967086099, 3.2177423254231563e-6, 6.435484650846313e-6,
+              1.2870969301692625e-5, 2.574193860338525e-5],
       volume_flux = flux_chandrashekar)
   end
 
   @trixi_testset "elixir_eulermulti_two_interacting_blast_waves.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulermulti_two_interacting_blast_waves.jl"),
-      l2   = [1.28886761e+00, 8.27133526e+01, 3.50680272e-03, 1.36987844e-02, 1.91795185e-02],
-      linf = [2.96413045e+01, 1.32258448e+03, 9.19191937e-02, 3.10929710e-01, 4.41798976e-01],
+      l2   = [1.288867611915533, 82.71335258388848, 0.00350680272313187, 0.013698784353152794,
+              0.019179518517518084],
+      linf = [29.6413044707026, 1322.5844802186496, 0.09191919374782143, 0.31092970966717925,
+              0.4417989757182038],
       tspan = (0.0, 0.0001))
   end
 
diff --git a/test/test_tree_1d_mhd.jl b/test/test_tree_1d_mhd.jl
index 938959831c1..e3a0cda3250 100644
--- a/test/test_tree_1d_mhd.jl
+++ b/test/test_tree_1d_mhd.jl
@@ -32,14 +32,14 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem")
 
   @trixi_testset "elixir_mhd_ec.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_ec.jl"),
-      l2   = [5.86009540e-02, 8.16048158e-02, 5.46791194e-02, 5.46791194e-02, 1.54509265e-01, 4.13046273e-17, 5.47637521e-02, 5.47637521e-02],
-      linf = [1.10014999e-01, 1.81982581e-01, 9.13611439e-02, 9.13611439e-02, 4.23831370e-01, 1.11022302e-16, 9.93731761e-02, 9.93731761e-02])
+      l2   = [0.05815183849746399, 0.08166807325621023, 0.054659228513541165, 0.054659228513541165, 0.15578125987042743, 4.130462730494e-17, 0.05465258887150046, 0.05465258887150046],
+      linf = [0.12165312668363826, 0.1901920742264952, 0.10059813883022554, 0.10059813883022554, 0.44079257431070706, 1.1102230246251565e-16, 0.10528911365809579, 0.10528911365809579])
   end
 
   @trixi_testset "elixir_mhd_briowu_shock_tube.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_briowu_shock_tube.jl"),
-      l2   = [0.17764301067932906, 0.19693621875378622, 0.3635136528288642, 0.0, 0.3757321708837591, 8.593007507325741e-16, 0.36473438378159656, 0.0],
-      linf = [0.5601530250396535, 0.43867368105486537, 1.0960903616351099, 0.0, 1.0551794137886303, 4.107825191113079e-15, 1.5374410890043144, 0.0],
+      l2   = [0.17477712356961989, 0.19489623595086944, 0.3596546157640463, 0.0, 0.3723215736814466, 1.2060075775846403e-15, 0.36276754492568164, 0.0],
+      linf = [0.5797109945880677, 0.4372991899547103, 1.0906536287185835, 0.0, 1.0526758874956808, 5.995204332975845e-15, 1.5122922036932964, 0.0],
       coverage_override = (maxiters=6,))
   end
 
@@ -51,8 +51,8 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem")
 
   @trixi_testset "elixir_mhd_ryujones_shock_tube.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_ryujones_shock_tube.jl"),
-      l2   = [2.34809441e-01, 3.92255943e-01, 8.23575546e-02, 1.75599624e-01, 9.61613519e-01, 6.60825891e-17, 2.15346454e-01, 1.07006529e-01],
-      linf = [6.40732148e-01, 9.44889516e-01, 3.54932707e-01, 8.54060243e-01, 2.07757711e+00, 1.11022302e-16, 4.92584725e-01, 2.49526561e-01],
+      l2   = [0.23469781891518154, 0.3916675299696121, 0.08245195301016353, 0.1745346945706147, 0.9606363432904367, 6.608258910237605e-17, 0.21542929107153735, 0.10705457908737925],
+      linf = [0.6447951791685409, 0.9461857095377463, 0.35074627554617605, 0.8515177411529542, 2.0770652030507053, 1.1102230246251565e-16, 0.49670855513788204, 0.24830199967863564],
       tspan = (0.0, 0.1))
   end
 
diff --git a/test/test_tree_1d_mhdmulti.jl b/test/test_tree_1d_mhdmulti.jl
index 2985e6d5663..5214ed26d38 100644
--- a/test/test_tree_1d_mhdmulti.jl
+++ b/test/test_tree_1d_mhdmulti.jl
@@ -11,33 +11,33 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem")
 
   @trixi_testset "elixir_mhdmulti_ec.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhdmulti_ec.jl"),
-      l2   = [0.08160481582829862, 0.05467911944326103, 0.05467911944326103, 0.15450926504459692,
-              4.130462730494e-17, 0.054763752050210085, 0.054763752050210085, 0.008371564857135208,
-              0.016743129714270416, 0.03348625942854083],
-      linf = [0.18198258075330706, 0.09136114386311774, 0.09136114386311774, 0.423831369951313,
-              1.1102230246251565e-16, 0.09937317613143604, 0.09937317613143604, 0.0157164284712992,
-              0.0314328569425984, 0.0628657138851968])
+      l2   = [0.08166807325620999, 0.054659228513541616, 0.054659228513541616, 0.15578125987042812,
+              4.130462730494e-17, 0.054652588871500665, 0.054652588871500665, 0.008307405499637766,
+              0.01661481099927553, 0.03322962199855106],
+      linf = [0.19019207422649645, 0.10059813883022888, 0.10059813883022888, 0.4407925743107146,
+              1.1102230246251565e-16, 0.10528911365809623, 0.10528911365809623, 0.01737901809766182,
+              0.03475803619532364, 0.06951607239064728])
   end
 
   @trixi_testset "elixir_mhdmulti_ec.jl with flux_derigs_etal" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhdmulti_ec.jl"),
-      l2   = [0.08153372259925547, 0.05464109003345891, 0.05464109003345891, 0.1540576724164453,
-              4.130462730494e-17, 0.054734930802131036, 0.054734930802131036, 0.008391254781284321,
-              0.016782509562568642, 0.033565019125137284],
-      linf = [0.17492544007323832, 0.09029632168248182, 0.09029632168248182, 0.40798609353896564,
-              1.1102230246251565e-16, 0.09872923637833075, 0.09872923637833075, 0.01609818847160674,
-              0.03219637694321348, 0.06439275388642696],
+      l2   = [0.08151404166186461, 0.054640238302693274, 0.054640238302693274, 0.15536125426328573,
+              4.130462730494e-17, 0.054665489963920275, 0.054665489963920275, 0.008308349501359825,
+              0.01661669900271965, 0.0332333980054393],
+      linf = [0.1824424257860952, 0.09734687137001484, 0.09734687137001484, 0.4243089502087325,
+              1.1102230246251565e-16, 0.09558639591092555, 0.09558639591092555, 0.017364773041550624,
+              0.03472954608310125, 0.0694590921662025],
       volume_flux = flux_derigs_etal)
   end
 
   @trixi_testset "elixir_mhdmulti_es.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhdmulti_es.jl"),
-      l2   = [0.07968782477167513, 0.05398115008116676, 0.05398115008116676, 0.15015281822439228,
-              4.130462730494e-17, 0.053629890024921495, 0.053629890024921495, 0.008279068245579706,
-              0.016558136491159413, 0.033116272982318826],
-      linf = [0.14118014632124837, 0.07820697032983395, 0.07820697032983395, 0.3390558674728652,
-              1.1102230246251565e-16, 0.06998787893467828, 0.06998787893467828, 0.014943825414763745,
-              0.02988765082952749, 0.05977530165905498])
+      l2   = [0.07994082660130175, 0.053940174914031976, 0.053940174914031976, 0.15165513559250643,
+              4.130462730494e-17, 0.05363207135290325, 0.05363207135290325, 0.008258265884659555,
+              0.01651653176931911, 0.03303306353863822],
+      linf = [0.14101014428198477, 0.07762441749521025, 0.07762441749521025, 0.3381334453289866,
+              1.1102230246251565e-16, 0.07003646400675223, 0.07003646400675223, 0.014962483760600165,
+              0.02992496752120033, 0.05984993504240066])
   end
 
   @trixi_testset "elixir_mhdmulti_convergence.jl" begin
@@ -52,12 +52,12 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem")
 
   @trixi_testset "elixir_mhdmulti_briowu_shock_tube.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhdmulti_briowu_shock_tube.jl"),
-      l2   = [0.1946577804333822, 0.3591196215528672, 0.0, 0.36875476066849383,
-              4.7644020131827105e-16, 0.36668249926193885, 0.0, 0.05775369214541893,
-              0.11550738429083786],
-      linf = [0.4345551123140612, 1.0874941615375844, 0.0, 1.0493729052116585,
-              3.219646771412954e-15, 1.5160434573973656, 0.0, 0.18616213071936066,
-              0.3723242614387213],
+      l2   = [0.1877830835572639, 0.3455841730726793, 0.0, 0.35413123388836687,
+              8.745556626531982e-16, 0.3629920109231055, 0.0, 0.05329005553971236,
+              0.10658011107942472],
+      linf = [0.4288187627971754, 1.0386547815614993, 0.0, 0.9541678878162702,
+              5.773159728050814e-15, 1.4595119339458051, 0.0, 0.18201910908829552,
+              0.36403821817659104],
       coverage_override = (maxiters=6,))
     end
 
diff --git a/test/test_tree_1d_shallowwater.jl b/test/test_tree_1d_shallowwater.jl
index f8901a3dcb6..1c3bac1fab6 100644
--- a/test/test_tree_1d_shallowwater.jl
+++ b/test/test_tree_1d_shallowwater.jl
@@ -10,22 +10,30 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem")
 @testset "Shallow Water" begin
   @trixi_testset "elixir_shallowwater_ec.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_ec.jl"),
-      l2   = [0.8122354510732459, 1.01586214815876, 0.43404255061704217],
-      linf = [1.4883285368551107, 3.8717508164234276, 1.7711213427919539],
+      l2   = [0.244729018751225, 0.8583565222389505, 0.07330427577586297],
+      linf = [2.1635021283528504, 3.8717508164234453, 1.7711213427919539],
+      tspan = (0.0, 0.25))
+  end
+
+  @trixi_testset "elixir_shallowwater_ec.jl with initial_condition_weak_blast_wave" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_ec.jl"),
+      l2   = [0.39464782107209717, 2.03880864210846, 4.1623084150546725e-10],
+      linf = [0.778905801278281, 3.2409883402608273, 7.419800190922032e-10],
+      initial_condition=initial_condition_weak_blast_wave,
       tspan = (0.0, 0.25))
   end
 
   @trixi_testset "elixir_shallowwater_well_balanced.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_well_balanced.jl"),
-      l2   = [1.2427984842961743, 1.0332499675061871e-14, 1.2427984842961741],
-      linf = [1.619041478244762, 1.266865149831811e-14, 1.6190414782447629],
+      l2   = [0.10416666834254829, 1.4352935256803184e-14, 0.10416666834254838],
+      linf = [1.9999999999999996, 3.248036646353028e-14, 2.0],
       tspan = (0.0, 0.25))
   end
 
   @trixi_testset "elixir_shallowwater_well_balanced.jl with FluxHydrostaticReconstruction" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_well_balanced.jl"),
-      l2   = [1.2427984842961743, 1.2663646513352053e-14, 1.2427984842961741],
-      linf = [1.619041478244762, 2.4566658711604395e-14, 1.6190414782447629],
+      l2   = [0.10416666834254835, 1.1891029971551825e-14, 0.10416666834254838],
+      linf = [2.0000000000000018, 2.4019608337954543e-14, 2.0],
       surface_flux=(FluxHydrostaticReconstruction(flux_lax_friedrichs, hydrostatic_reconstruction_audusse_etal), flux_nonconservative_audusse_etal),
       tspan = (0.0, 0.25))
   end
@@ -59,14 +67,14 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem")
       tspan = (0.0, 0.025))
   end
 
-  @trixi_testset "elixir_shallowwater_well_balanced_nonperiodic.jl with dirichlet boundary" begin
+  @trixi_testset "elixir_shallowwater_well_balanced_nonperiodic.jl with Dirichlet boundary" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_well_balanced_nonperiodic.jl"),
       l2   = [1.725964362045055e-8, 5.0427180314307505e-16, 1.7259643530442137e-8],
       linf = [3.844551077492042e-8, 3.469453422316143e-15, 3.844551077492042e-8],
       tspan = (0.0, 0.25))
   end
 
-  @trixi_testset "elixir_shallowwater_well_nonperiodic.jl with wall boundary" begin
+  @trixi_testset "elixir_shallowwater_well_balanced_nonperiodic.jl with wall boundary" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_well_balanced_nonperiodic.jl"),
       l2   = [1.7259643614361866e-8, 3.5519018243195145e-16, 1.7259643530442137e-8],
       linf = [3.844551010878661e-8, 9.846474508971374e-16, 3.844551077492042e-8],
@@ -76,8 +84,8 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem")
 
   @trixi_testset "elixir_shallowwater_shock_capturing.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_shock_capturing.jl"),
-      l2   = [0.2884024818919076, 0.5252262013521178, 0.2890348477852955],
-      linf = [0.7565706154863958, 2.076621603471687, 0.8646939843534258],
+      l2   = [0.07424140641160326, 0.2148642632748155, 0.0372579849000542],
+      linf = [1.1209754279344226, 1.3230788645853582, 0.8646939843534251],
       tspan = (0.0, 0.05))
   end
 end
diff --git a/test/test_tree_1d_shallowwater_twolayer.jl b/test/test_tree_1d_shallowwater_twolayer.jl
index 6c0ad2941cc..0d8a83806f9 100644
--- a/test/test_tree_1d_shallowwater_twolayer.jl
+++ b/test/test_tree_1d_shallowwater_twolayer.jl
@@ -38,9 +38,9 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem")
 
   @trixi_testset "elixir_shallowwater_twolayer_dam_break.jl with flux_lax_friedrichs" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_twolayer_dam_break.jl"),
-    l2    = [0.35490827242437256, 1.6715402155795918, 0.6960264969949427,
-             0.9351481433409805, 0.7938172946965545],
-    linf  = [0.6417127471419837, 1.9742107034120873, 1.135774587483082, 1.236125279347084, 1.1],
+    l2    = [0.10010269243463918, 0.5668733957648654, 0.08759617327649398,
+             0.4538443183566172, 0.013638618139749523],
+    linf  = [0.5854202777756559, 2.1278930820498934, 0.5193686074348809, 1.8071213168086229, 0.5],
     surface_flux = (flux_lax_friedrichs, flux_nonconservative_fjordholm_etal),
     tspan = (0.0, 0.25))
   end