From 0975ed6d8d381c2b276a3b37329e57277f832a6f Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Thu, 18 Jul 2024 16:15:21 +0200
Subject: [PATCH 01/30] Pick up where Gregor left

---
 src/Trixi.jl                               |   2 +-
 src/solvers/dg.jl                          |  58 ++++++++++++
 src/solvers/dgsem_tree/dg.jl               |   3 +
 src/solvers/dgsem_tree/dg_1d.jl            | 103 ++++++++++++++++++++-
 src/solvers/dgsem_tree/finite_volume_O2.jl |  76 +++++++++++++++
 5 files changed, 239 insertions(+), 3 deletions(-)
 create mode 100644 src/solvers/dgsem_tree/finite_volume_O2.jl

diff --git a/src/Trixi.jl b/src/Trixi.jl
index 1a509ed92d1..08bbfea1e16 100644
--- a/src/Trixi.jl
+++ b/src/Trixi.jl
@@ -234,7 +234,7 @@ export DG,
        FDSBP,
        VolumeIntegralWeakForm, VolumeIntegralStrongForm,
        VolumeIntegralFluxDifferencing,
-       VolumeIntegralPureLGLFiniteVolume,
+       VolumeIntegralPureLGLFiniteVolume, VolumeIntegralPureLGLFiniteVolumeO2,
        VolumeIntegralShockCapturingHG, IndicatorHennemannGassner,
        VolumeIntegralUpwind,
        SurfaceIntegralWeakForm, SurfaceIntegralStrongForm,
diff --git a/src/solvers/dg.jl b/src/solvers/dg.jl
index fb4c8f182e0..ce759c8e298 100644
--- a/src/solvers/dg.jl
+++ b/src/solvers/dg.jl
@@ -184,6 +184,64 @@ function Base.show(io::IO, ::MIME"text/plain",
     end
 end
 
+"""
+    VolumeIntegralPureLGLFiniteVolumeO2(basis::Basis;
+                                        volume_flux_fv = flux_lax_friedrichs,
+                                        reconstruction_mode = reconstruction_small_stencil,
+                                        slope_limiter = minmod)
+
+This gives a formally O(2)-accurate finite volume scheme on an LGL-type subcell
+mesh (LGL = Legendre-Gauss-Lobatto).
+
+!!! warning "Experimental implementation"
+    This is an experimental feature and may change in future releases.
+
+## References
+
+- Hennemann, Gassner (2020)
+  "A provably entropy stable subcell shock capturing approach for high order split form DG"
+  [arXiv: 2008.12044](https://arxiv.org/abs/2008.12044)
+"""
+struct VolumeIntegralPureLGLFiniteVolumeO2{RealT, Basis, VolumeFluxFV, Reconstruction,
+                                           Limiter} <: AbstractVolumeIntegral
+    x_interfaces::Vector{RealT} # x-coordinates of the sub-cell element interfaces
+    volume_flux_fv::VolumeFluxFV # non-symmetric in general, e.g. entropy-dissipative
+    reconstruction_mode::Reconstruction # which type of FV reconstruction to use
+    slope_limiter::Limiter # which type of slope limiter function
+end
+
+function VolumeIntegralPureLGLFiniteVolumeO2(basis::Basis;
+                                             volume_flux_fv = flux_lax_friedrichs,
+                                             reconstruction_mode = reconstruction_small_stencil,
+                                             slope_limiter = minmod) where {Basis}
+    # Suffices to store only the intermediate boundaries of the sub-cell elements                                             
+    x_interfaces = cumsum(basis.weights)[1:(end - 1)] .- 1
+    VolumeIntegralPureLGLFiniteVolumeO2{eltype(basis.weights),
+                                        typeof(basis),
+                                        typeof(volume_flux_fv),
+                                        typeof(reconstruction_mode),
+                                        typeof(slope_limiter)}(x_interfaces,
+                                                               volume_flux_fv,
+                                                               reconstruction_mode,
+                                                               slope_limiter)
+end
+
+function Base.show(io::IO, ::MIME"text/plain",
+                   integral::VolumeIntegralPureLGLFiniteVolumeO2)
+    @nospecialize integral # reduce precompilation time
+
+    if get(io, :compact, false)
+        show(io, integral)
+    else
+        setup = [
+            "FV flux" => integral.volume_flux_fv,
+            "Reconstruction" => integral.reconstruction_mode,
+            "Slope limiter" => integral.slope_limiter,
+        ]
+        summary_box(io, "VolumeIntegralPureLGLFiniteVolumeO2", setup)
+    end
+end
+
 """
     VolumeIntegralSubcellLimiting(limiter;
                                   volume_flux_dg, volume_flux_fv)
diff --git a/src/solvers/dgsem_tree/dg.jl b/src/solvers/dgsem_tree/dg.jl
index 0993b3c9b85..36c93e28461 100644
--- a/src/solvers/dgsem_tree/dg.jl
+++ b/src/solvers/dgsem_tree/dg.jl
@@ -63,6 +63,9 @@ include("dg_parallel.jl")
 # Helper structs for parabolic AMR
 include("containers_viscous.jl")
 
+# Some standard functions for the canonical second-order FV (MUSCL) scheme
+include("finite_volume_O2.jl")
+
 # 1D DG implementation
 include("dg_1d.jl")
 include("dg_1d_parabolic.jl")
diff --git a/src/solvers/dgsem_tree/dg_1d.jl b/src/solvers/dgsem_tree/dg_1d.jl
index 4a0747d1c09..c8345149a43 100644
--- a/src/solvers/dgsem_tree/dg_1d.jl
+++ b/src/solvers/dgsem_tree/dg_1d.jl
@@ -61,8 +61,9 @@ end
 
 function create_cache(mesh::Union{TreeMesh{1}, StructuredMesh{1}, P4estMesh{1}},
                       equations,
-                      volume_integral::VolumeIntegralPureLGLFiniteVolume, dg::DG,
-                      uEltype)
+                      volume_integral::Union{VolumeIntegralPureLGLFiniteVolume,
+                                             VolumeIntegralPureLGLFiniteVolumeO2},
+                      dg::DG, uEltype)
     A2dp1_x = Array{uEltype, 2}
     fstar1_L_threaded = A2dp1_x[A2dp1_x(undef, nvariables(equations), nnodes(dg) + 1)
                                 for _ in 1:Threads.nthreads()]
@@ -309,6 +310,24 @@ function calc_volume_integral!(du, u,
     return nothing
 end
 
+function calc_volume_integral!(du, u,
+                               mesh::TreeMesh{1},
+                               nonconservative_terms, equations,
+                               volume_integral::VolumeIntegralPureLGLFiniteVolumeO2,
+                               dg::DGSEM, cache)
+    @unpack x_interfaces, volume_flux_fv, reconstruction_mode, slope_limiter = volume_integral
+
+    # Calculate LGL second-order FV volume integral
+    @threaded for element in eachelement(dg, cache)
+        fv_kernel!(du, u, mesh, nonconservative_terms, equations, volume_flux_fv,
+                   dg, cache, element,
+                   x_interfaces, reconstruction_mode, slope_limiter,
+                   true)
+    end
+
+    return nothing
+end
+
 @inline function fv_kernel!(du, u,
                             mesh::Union{TreeMesh{1}, StructuredMesh{1}},
                             nonconservative_terms, equations,
@@ -335,6 +354,35 @@ end
     return nothing
 end
 
+@inline function fv_kernel!(du, u,
+                            mesh::Union{TreeMesh{1}, StructuredMesh{1}},
+                            nonconservative_terms, equations,
+                            volume_flux_fv, dg::DGSEM, cache, element,
+                            x_interfaces, reconstruction_mode, slope_limiter,
+                            alpha = true)
+    @unpack fstar1_L_threaded, fstar1_R_threaded = cache
+    @unpack inverse_weights = dg.basis
+
+    # Calculate FV two-point fluxes
+    fstar1_L = fstar1_L_threaded[Threads.threadid()]
+    fstar1_R = fstar1_R_threaded[Threads.threadid()]
+    calcflux_fv!(fstar1_L, fstar1_R, u, mesh, nonconservative_terms, equations,
+                 volume_flux_fv,
+                 dg, element, cache,
+                 x_interfaces, reconstruction_mode, slope_limiter)
+
+    # Calculate FV volume integral contribution
+    for i in eachnode(dg)
+        for v in eachvariable(equations)
+            du[v, i, element] += (alpha *
+                                  (inverse_weights[i] *
+                                   (fstar1_L[v, i + 1] - fstar1_R[v, i])))
+        end
+    end
+
+    return nothing
+end
+
 @inline function calcflux_fv!(fstar1_L, fstar1_R, u::AbstractArray{<:Any, 3},
                               mesh::Union{TreeMesh{1}, StructuredMesh{1}},
                               nonconservative_terms::False,
@@ -388,6 +436,57 @@ end
     return nothing
 end
 
+@inline function calcflux_fv!(fstar1_L, fstar1_R, u::AbstractArray{<:Any, 3},
+                              mesh::Union{TreeMesh{1}, StructuredMesh{1}},
+                              nonconservative_terms::False,
+                              equations, volume_flux_fv, dg::DGSEM, element, cache,
+                              x_interfaces, reconstruction_mode, slope_limiter)
+    fstar1_L[:, 1] .= zero(eltype(fstar1_L))
+    fstar1_L[:, nnodes(dg) + 1] .= zero(eltype(fstar1_L))
+    fstar1_R[:, 1] .= zero(eltype(fstar1_R))
+    fstar1_R[:, nnodes(dg) + 1] .= zero(eltype(fstar1_R))
+
+    for i in 2:nnodes(dg)
+        #             Reference element:             
+        #  -1 -----------------0----------------- 1 -> x
+        # Gauss Lobatto Legendre nodes (schematic for k = 3):
+        #   .         .                 .         .
+        #   ^         ^                 ^         ^
+        # i - 2,    i - 1,               i,     i + 1
+        # mm        ll                   rr     pp
+        # Cell boundaries (schematic for k = 3): 
+        # (note that only the inner three boundaries are stored)
+        #  -1 -----------------0----------------- 1 -> x
+        #   |     |            |             |    |
+        # Cell index:
+        #      1         2              3       4
+
+        # Obtain unlimited values in primal variables
+        # TODO: If i - 2 = 0 go to neighbor element?
+        u_mm = cons2prim(get_node_vars(u, equations, dg, max(1, i - 2), element),
+                         equations)
+        u_ll = cons2prim(get_node_vars(u, equations, dg, i - 1, element), equations)
+        u_rr = cons2prim(get_node_vars(u, equations, dg, i, element), equations)
+        # TODO: If i + 1 > nnodes(dg) go to neighbor element?
+        u_pp = cons2prim(get_node_vars(u, equations, dg, min(nnodes(dg), i + 1),
+                                       element), equations)
+
+        # Reconstruct values with limiting
+        u_ll, u_rr = reconstruction_mode(u_mm, u_ll, u_rr, u_pp,
+                                         x_interfaces,
+                                         i, slope_limiter, dg)
+
+        # Convert primal variables back to conservative variables
+        flux = volume_flux_fv(prim2cons(u_ll, equations), prim2cons(u_rr, equations),
+                              1, equations) # orientation 1: x direction
+
+        set_node_vars!(fstar1_L, flux, equations, dg, i)
+        set_node_vars!(fstar1_R, flux, equations, dg, i)
+    end
+
+    return nothing
+end
+
 # We pass the `surface_integral` argument solely for dispatch
 function prolong2interfaces!(cache, u,
                              mesh::TreeMesh{1}, equations, surface_integral, dg::DG)
diff --git a/src/solvers/dgsem_tree/finite_volume_O2.jl b/src/solvers/dgsem_tree/finite_volume_O2.jl
new file mode 100644
index 00000000000..726979d3408
--- /dev/null
+++ b/src/solvers/dgsem_tree/finite_volume_O2.jl
@@ -0,0 +1,76 @@
+@inline function linear_reconstruction(u_ll, u_rr, ux_ll, ux_rr, x_ll, x_rr, x_interfaces, node_index)
+    # Linear reconstruction at the interface
+    u_ll = u_ll + ux_ll * (x_interfaces[node_index - 1] - x_ll)
+    u_rr = u_rr + ux_rr * (x_interfaces[node_index - 1] - x_rr)
+
+    return u_ll, u_rr
+end
+
+@inline function reconstruction_small_stencil(u_mm, u_ll, u_rr, u_pp,
+                                              x_interfaces,
+                                              node_index, limiter, dg)
+    #             Reference element:             
+    #  -1 -----------------0----------------- 1 -> x
+    # Gauss Lobatto Legendre nodes (schematic for k = 3):
+    #   .         .                 .         .
+    #   ^         ^                 ^         ^
+    # i - 2,    i - 1,               i,     i + 1
+    # mm       ll                   rr     pp
+    # Cell boundaries (schematic for k = 3): 
+    # (note that only the inner three boundaries are stored)
+    #  -1 -----------------0----------------- 1 -> x
+    #   |     |            |             |    |
+    # Cell index:
+    #      1         2              3       4
+
+    @unpack nodes = dg.basis
+    x_ll = nodes[node_index - 1]
+    x_rr = nodes[node_index]
+
+    # Middle element slope
+    ux_m = (u_rr - u_ll) / (x_rr - x_ll)
+    if node_index == 2 # Catch case mm == ll
+        ux_ll = ux_m
+    else
+        # Left element slope
+        ux_l = (u_ll - u_mm) / (x_ll - nodes[node_index - 2])
+        ux_ll = limiter.(ux_l, ux_m)
+    end
+
+    if node_index == nnodes(dg) # Catch case rr == pp
+        ux_rr = ux_m
+    else
+        # Right element slope
+        ux_r = (u_pp - u_rr) / (nodes[node_index + 1] - x_rr)
+        ux_rr = limiter.(ux_m, ux_r)
+    end
+
+    linear_reconstruction(u_ll, u_rr, ux_ll, ux_rr, x_ll, x_rr, x_interfaces, node_index)
+end
+
+@inline function reconstruction_constant(u_mm, u_ll, u_rr, u_pp,
+                                         x_interfaces,
+                                         node_index, limiter, dg)
+    return u_ll, u_rr
+end
+
+@inline function central_recon(sl, sr)
+    s = 0.5 * (sl + sr)
+    return s
+end
+
+@inline function minmod(sl, sr)
+    s = 0.0
+    if sign(sl) == sign(sr)
+        s = sign(sl) * min(abs(sl), abs(sr))
+    end
+    return s
+end
+
+@inline function monotonized_central(sl, sr)
+    s = 0.0
+    if sign(sl) == sign(sr)
+        s = sign(sl) * min(2 * abs(sl), 2 * abs(sr), 0.5 * abs(sl + sr))
+    end
+    return s
+end

From 5c8c32786f1e86807d099299c4e98b37ccc7968e Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Thu, 18 Jul 2024 16:32:15 +0200
Subject: [PATCH 02/30] preliminary example

---
 .../elixir_euler_convergence_pure_fvO2.jl     | 72 +++++++++++++++++++
 src/solvers/dgsem_tree/dg_1d.jl               |  2 +-
 2 files changed, 73 insertions(+), 1 deletion(-)
 create mode 100644 examples/structured_1d_dgsem/elixir_euler_convergence_pure_fvO2.jl

diff --git a/examples/structured_1d_dgsem/elixir_euler_convergence_pure_fvO2.jl b/examples/structured_1d_dgsem/elixir_euler_convergence_pure_fvO2.jl
new file mode 100644
index 00000000000..3d96a725098
--- /dev/null
+++ b/examples/structured_1d_dgsem/elixir_euler_convergence_pure_fvO2.jl
@@ -0,0 +1,72 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+
+equations = CompressibleEulerEquations1D(1.4)
+
+initial_condition = initial_condition_convergence_test
+
+polydeg = 3
+basis = LobattoLegendreBasis(polydeg)
+solver = DGSEM(polydeg = polydeg, surface_flux = flux_hllc,
+               volume_integral = VolumeIntegralPureLGLFiniteVolumeO2(basis, volume_flux_fv = flux_hllc))
+
+#=               
+coordinates_min = (0.0,) # minimum coordinate
+coordinates_max = (2.0,) # maximum coordinate
+cells_per_dimension = (16,)
+
+# Create curved mesh with 16 cells
+mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max)
+=#
+
+coordinates_min = 0.0
+coordinates_max = 2.0
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 4,
+                n_cells_max = 10_000)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    source_terms = source_terms_convergence_test)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 2.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
+                                     extra_analysis_errors = (:l2_error_primitive,
+                                                              :linf_error_primitive))
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(interval = 100,
+                                     save_initial_solution = true,
+                                     save_final_solution = true,
+                                     solution_variables = cons2prim)
+
+stepsize_callback = StepsizeCallback(cfl = 0.8)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback,
+                        save_solution,
+                        stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
+            dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep = false, callback = callbacks);
+summary_callback() # print the timer summary
+
+using Plots
+
+plot(sol)
\ No newline at end of file
diff --git a/src/solvers/dgsem_tree/dg_1d.jl b/src/solvers/dgsem_tree/dg_1d.jl
index c8345149a43..0d52d90e297 100644
--- a/src/solvers/dgsem_tree/dg_1d.jl
+++ b/src/solvers/dgsem_tree/dg_1d.jl
@@ -311,7 +311,7 @@ function calc_volume_integral!(du, u,
 end
 
 function calc_volume_integral!(du, u,
-                               mesh::TreeMesh{1},
+                               mesh::Union{TreeMesh{1}, StructuredMesh{1}},
                                nonconservative_terms, equations,
                                volume_integral::VolumeIntegralPureLGLFiniteVolumeO2,
                                dg::DGSEM, cache)

From 08784226dbcb365ca9ddd0b16921cf27b7ea21a1 Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Fri, 19 Jul 2024 12:18:38 +0200
Subject: [PATCH 03/30] more limiters

---
 .../elixir_euler_convergence_pure_fvO2.jl     |   3 +-
 src/Trixi.jl                                  |   2 +
 src/solvers/dgsem_tree/finite_volume_O2.jl    | 107 +++++++++++++++++-
 3 files changed, 108 insertions(+), 4 deletions(-)

diff --git a/examples/structured_1d_dgsem/elixir_euler_convergence_pure_fvO2.jl b/examples/structured_1d_dgsem/elixir_euler_convergence_pure_fvO2.jl
index 3d96a725098..5a083464048 100644
--- a/examples/structured_1d_dgsem/elixir_euler_convergence_pure_fvO2.jl
+++ b/examples/structured_1d_dgsem/elixir_euler_convergence_pure_fvO2.jl
@@ -12,7 +12,8 @@ initial_condition = initial_condition_convergence_test
 polydeg = 3
 basis = LobattoLegendreBasis(polydeg)
 solver = DGSEM(polydeg = polydeg, surface_flux = flux_hllc,
-               volume_integral = VolumeIntegralPureLGLFiniteVolumeO2(basis, volume_flux_fv = flux_hllc))
+               volume_integral = VolumeIntegralPureLGLFiniteVolumeO2(basis, volume_flux_fv = flux_hllc,
+               slope_limiter = superbee))
 
 #=               
 coordinates_min = (0.0,) # minimum coordinate
diff --git a/src/Trixi.jl b/src/Trixi.jl
index 08bbfea1e16..8dfd4dc1788 100644
--- a/src/Trixi.jl
+++ b/src/Trixi.jl
@@ -241,6 +241,8 @@ export DG,
        SurfaceIntegralUpwind,
        MortarL2
 
+export minmod, monotonized_central, superbee, vanLeer_limiter
+
 export VolumeIntegralSubcellLimiting, BoundsCheckCallback,
        SubcellLimiterIDP, SubcellLimiterIDPCorrection
 
diff --git a/src/solvers/dgsem_tree/finite_volume_O2.jl b/src/solvers/dgsem_tree/finite_volume_O2.jl
index 726979d3408..89468875379 100644
--- a/src/solvers/dgsem_tree/finite_volume_O2.jl
+++ b/src/solvers/dgsem_tree/finite_volume_O2.jl
@@ -48,29 +48,130 @@ end
     linear_reconstruction(u_ll, u_rr, ux_ll, ux_rr, x_ll, x_rr, x_interfaces, node_index)
 end
 
+"""
+    reconstruction_constant(u_mm, u_ll, u_rr, u_pp,
+                            x_interfaces,
+                            node_index, limiter, dg)
+
+Returns the constant "reconstructed" values at the interface `x_interfaces[node_index - 1]`
+obtained from constant polynomials.
+Formally O(1) accurate.
+"""
 @inline function reconstruction_constant(u_mm, u_ll, u_rr, u_pp,
                                          x_interfaces,
                                          node_index, limiter, dg)
     return u_ll, u_rr
 end
 
+"""
+    central_recon(sl, sr)
+
+Central, non-TVD reconstruction given left and right slopes `sl` and `sr`.
+Gives formally full order of accuracy at the expense of sacrificied stability.
+"""    
 @inline function central_recon(sl, sr)
     s = 0.5 * (sl + sr)
     return s
 end
 
+"""
+    minmod(sl, sr)
+
+Classic minmod limiter function for a TVD reconstruction given left and right slopes `sl` and `sr`.
+There are many different ways how the minmod limiter can be implemented.
+For reference, see for instance Eq. (6.27) in
+
+- Randall J. LeVeque (2002)
+  Finite Volume Methods for Hyperbolic Problems
+  [DOI: 10.1017/CBO9780511791253](https://doi.org/10.1017/CBO9780511791253)
+"""
 @inline function minmod(sl, sr)
-    s = 0.0
+    #=
     if sign(sl) == sign(sr)
         s = sign(sl) * min(abs(sl), abs(sr))
+    else
+        s = 0.0
     end
+
     return s
+    =#
+
+    return 0.5 * (sign(sl) + sign(sr)) * min(abs(sl), abs(sr))
 end
 
+"""
+    monotonized_central(sl, sr)
+
+Monotonized central limiter function for a TVD reconstruction given left and right slopes `sl` and `sr`.
+There are many different ways how the monotonized central limiter can be implemented.
+For reference, see for instance Eq. (6.29) in
+
+- Randall J. LeVeque (2002)
+  Finite Volume Methods for Hyperbolic Problems
+  [DOI: 10.1017/CBO9780511791253](https://doi.org/10.1017/CBO9780511791253)
+"""
 @inline function monotonized_central(sl, sr)
-    s = 0.0
+    # CARE: MC assumes equidistant grid in 0.5 * abs(sl + sr)!
+    #=
+    if sign(sl) == sign(sr)
+        s = sign(sl) * min(0.5 * abs(sl + sr), 2 * abs(sl), 2 * abs(sr))
+    else
+        s = 0.0
+    end
+    =#
+
+    # Use recursive property of minmod function
+    s = minmod(0.5 * (sl + sr), minmod(2 * sl, 2 * sr))
+
+    return s
+end
+
+@inline function maxmod(sl, sr)
+    #=
     if sign(sl) == sign(sr)
-        s = sign(sl) * min(2 * abs(sl), 2 * abs(sr), 0.5 * abs(sl + sr))
+        s = sign(sl) * max(abs(sl), abs(sr))
+    else
+        s = 0.0
     end
+
+    return s
+    =#
+
+    return 0.5 * (sign(sl) + sign(sr)) * max(abs(sl), abs(sr))
+end
+
+"""
+    superbee(sl, sr)
+
+Superbee limiter function for a TVD reconstruction given left and right slopes `sl` and `sr`.
+There are many different ways how the monotonized central limiter can be implemented.
+For reference, see for instance Eq. (6.28) in
+
+- Randall J. LeVeque (2002)
+  Finite Volume Methods for Hyperbolic Problems
+  [DOI: 10.1017/CBO9780511791253](https://doi.org/10.1017/CBO9780511791253)
+"""
+@inline function superbee(sl, sr)
+    s = maxmod(minmod(sl, 2 * sr), minmod(sl, 2 * sr))
     return s
 end
+
+"""
+    vanLeer_limiter(sl, sr)
+
+Symmetric Limiter by van Leer.
+See for reference page 70 in 
+
+- Siddhartha Mishra, Ulrik Skre Fjordholm and Rémi Abgrall
+  Numerical methods for conservation laws and related equations.
+  [Link](https://metaphor.ethz.ch/x/2019/hs/401-4671-00L/literature/mishra_hyperbolic_pdes.pdf)
+"""
+@inline function vanLeer_limiter(sl, sr)
+    if abs(sl) + abs(sr) > 0.0
+        s = (abs(sr) * sl + abs(sl) * sr) / (abs(sl) + abs(sr))
+    else
+        s = 0.0
+    end
+
+    return s
+end
\ No newline at end of file

From 80dd41bffdb3ff509917e334d87745f1c0f8f34e Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Fri, 19 Jul 2024 12:20:26 +0200
Subject: [PATCH 04/30] comments

---
 src/solvers/dgsem_tree/finite_volume_O2.jl | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/src/solvers/dgsem_tree/finite_volume_O2.jl b/src/solvers/dgsem_tree/finite_volume_O2.jl
index 89468875379..c10d9dffefd 100644
--- a/src/solvers/dgsem_tree/finite_volume_O2.jl
+++ b/src/solvers/dgsem_tree/finite_volume_O2.jl
@@ -48,6 +48,9 @@ end
     linear_reconstruction(u_ll, u_rr, ux_ll, ux_rr, x_ll, x_rr, x_interfaces, node_index)
 end
 
+# TODO: Different reconstructions, see
+# https://github.com/trixi-framework/Trixi.jl/pull/433/files
+
 """
     reconstruction_constant(u_mm, u_ll, u_rr, u_pp,
                             x_interfaces,
@@ -126,6 +129,7 @@ For reference, see for instance Eq. (6.29) in
     return s
 end
 
+# Note: This is NOT a limiter, just a helper for the `superbee` limiter below.
 @inline function maxmod(sl, sr)
     #=
     if sign(sl) == sign(sr)

From 059c66071c6480301b47504461d9cda1a87e0779 Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Fri, 19 Jul 2024 13:47:11 +0200
Subject: [PATCH 05/30] fmt

---
 src/solvers/dgsem_tree/finite_volume_O2.jl | 47 +++++++++++-----------
 1 file changed, 23 insertions(+), 24 deletions(-)

diff --git a/src/solvers/dgsem_tree/finite_volume_O2.jl b/src/solvers/dgsem_tree/finite_volume_O2.jl
index c10d9dffefd..92ebc97d5be 100644
--- a/src/solvers/dgsem_tree/finite_volume_O2.jl
+++ b/src/solvers/dgsem_tree/finite_volume_O2.jl
@@ -1,14 +1,31 @@
-@inline function linear_reconstruction(u_ll, u_rr, ux_ll, ux_rr, x_ll, x_rr, x_interfaces, node_index)
+@inline function linear_reconstruction(u_ll, u_rr, ux_ll, ux_rr,
+                                       x_ll, x_rr, x_interfaces, node_index)
     # Linear reconstruction at the interface
     u_ll = u_ll + ux_ll * (x_interfaces[node_index - 1] - x_ll)
     u_rr = u_rr + ux_rr * (x_interfaces[node_index - 1] - x_rr)
 
     return u_ll, u_rr
 end
+# TODO: Different reconstructions, see
+# https://github.com/trixi-framework/Trixi.jl/pull/433/files
+
+"""
+    reconstruction_constant(u_mm, u_ll, u_rr, u_pp,
+                            x_interfaces,
+                            node_index, limiter, dg)
+
+Returns the constant "reconstructed" values at the interface `x_interfaces[node_index - 1]`
+obtained from constant polynomials.
+Formally O(1) accurate.
+"""
+@inline function reconstruction_constant(u_mm, u_ll, u_rr, u_pp,
+                                         x_interfaces,
+                                         node_index, limiter, dg)
+    return u_ll, u_rr
+end
 
 @inline function reconstruction_small_stencil(u_mm, u_ll, u_rr, u_pp,
-                                              x_interfaces,
-                                              node_index, limiter, dg)
+                                              x_interfaces, node_index, limiter, dg)
     #             Reference element:             
     #  -1 -----------------0----------------- 1 -> x
     # Gauss Lobatto Legendre nodes (schematic for k = 3):
@@ -48,30 +65,12 @@ end
     linear_reconstruction(u_ll, u_rr, ux_ll, ux_rr, x_ll, x_rr, x_interfaces, node_index)
 end
 
-# TODO: Different reconstructions, see
-# https://github.com/trixi-framework/Trixi.jl/pull/433/files
-
-"""
-    reconstruction_constant(u_mm, u_ll, u_rr, u_pp,
-                            x_interfaces,
-                            node_index, limiter, dg)
-
-Returns the constant "reconstructed" values at the interface `x_interfaces[node_index - 1]`
-obtained from constant polynomials.
-Formally O(1) accurate.
-"""
-@inline function reconstruction_constant(u_mm, u_ll, u_rr, u_pp,
-                                         x_interfaces,
-                                         node_index, limiter, dg)
-    return u_ll, u_rr
-end
-
 """
     central_recon(sl, sr)
 
 Central, non-TVD reconstruction given left and right slopes `sl` and `sr`.
-Gives formally full order of accuracy at the expense of sacrificied stability.
-"""    
+Gives formally full order of accuracy at the expense of sacrificied nonlinear stability.
+"""
 @inline function central_recon(sl, sr)
     s = 0.5 * (sl + sr)
     return s
@@ -178,4 +177,4 @@ See for reference page 70 in
     end
 
     return s
-end
\ No newline at end of file
+end

From 7f4afa678463f74061121892760c5f3ced99cf89 Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Tue, 23 Jul 2024 07:33:24 +0200
Subject: [PATCH 06/30] continue

---
 src/Trixi.jl                               |  3 ++-
 src/solvers/dgsem_tree/finite_volume_O2.jl | 31 +++++++++++-----------
 2 files changed, 17 insertions(+), 17 deletions(-)

diff --git a/src/Trixi.jl b/src/Trixi.jl
index 8dfd4dc1788..f5c237a9044 100644
--- a/src/Trixi.jl
+++ b/src/Trixi.jl
@@ -241,7 +241,8 @@ export DG,
        SurfaceIntegralUpwind,
        MortarL2
 
-export minmod, monotonized_central, superbee, vanLeer_limiter
+export reconstruction_small_stencil,
+       minmod, monotonized_central, superbee, vanLeer_limiter
 
 export VolumeIntegralSubcellLimiting, BoundsCheckCallback,
        SubcellLimiterIDP, SubcellLimiterIDPCorrection
diff --git a/src/solvers/dgsem_tree/finite_volume_O2.jl b/src/solvers/dgsem_tree/finite_volume_O2.jl
index 92ebc97d5be..27b2ce8df74 100644
--- a/src/solvers/dgsem_tree/finite_volume_O2.jl
+++ b/src/solvers/dgsem_tree/finite_volume_O2.jl
@@ -6,8 +6,6 @@
 
     return u_ll, u_rr
 end
-# TODO: Different reconstructions, see
-# https://github.com/trixi-framework/Trixi.jl/pull/433/files
 
 """
     reconstruction_constant(u_mm, u_ll, u_rr, u_pp,
@@ -24,28 +22,29 @@ Formally O(1) accurate.
     return u_ll, u_rr
 end
 
+#             Reference element:             
+#  -1 -----------------0----------------- 1 -> x
+# Gauss Lobatto Legendre nodes (schematic for k = 3):
+#   .         .                 .         .
+#   ^         ^                 ^         ^
+# i - 2,    i - 1,               i,     i + 1
+# mm       ll                   rr     pp
+# Cell boundaries (schematic for k = 3): 
+# (note that only the inner three boundaries are stored)
+#  -1 -----------------0----------------- 1 -> x
+#   |     |            |             |    |
+# Cell index:
+#      1         2              3       4
+
 @inline function reconstruction_small_stencil(u_mm, u_ll, u_rr, u_pp,
                                               x_interfaces, node_index, limiter, dg)
-    #             Reference element:             
-    #  -1 -----------------0----------------- 1 -> x
-    # Gauss Lobatto Legendre nodes (schematic for k = 3):
-    #   .         .                 .         .
-    #   ^         ^                 ^         ^
-    # i - 2,    i - 1,               i,     i + 1
-    # mm       ll                   rr     pp
-    # Cell boundaries (schematic for k = 3): 
-    # (note that only the inner three boundaries are stored)
-    #  -1 -----------------0----------------- 1 -> x
-    #   |     |            |             |    |
-    # Cell index:
-    #      1         2              3       4
-
     @unpack nodes = dg.basis
     x_ll = nodes[node_index - 1]
     x_rr = nodes[node_index]
 
     # Middle element slope
     ux_m = (u_rr - u_ll) / (x_rr - x_ll)
+
     if node_index == 2 # Catch case mm == ll
         ux_ll = ux_m
     else

From 196854a6bce2e6aac10dee01f331ebc9bbcee1aa Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Tue, 23 Jul 2024 19:54:02 +0200
Subject: [PATCH 07/30] comments

---
 .../elixir_euler_convergence_pure_fvO2.jl     | 11 +------
 src/solvers/dgsem_tree/finite_volume_O2.jl    | 30 ++++++++++++-------
 2 files changed, 20 insertions(+), 21 deletions(-)

diff --git a/examples/structured_1d_dgsem/elixir_euler_convergence_pure_fvO2.jl b/examples/structured_1d_dgsem/elixir_euler_convergence_pure_fvO2.jl
index 5a083464048..6aece3cf278 100644
--- a/examples/structured_1d_dgsem/elixir_euler_convergence_pure_fvO2.jl
+++ b/examples/structured_1d_dgsem/elixir_euler_convergence_pure_fvO2.jl
@@ -48,16 +48,10 @@ analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
 
 alive_callback = AliveCallback(analysis_interval = analysis_interval)
 
-save_solution = SaveSolutionCallback(interval = 100,
-                                     save_initial_solution = true,
-                                     save_final_solution = true,
-                                     solution_variables = cons2prim)
-
 stepsize_callback = StepsizeCallback(cfl = 0.8)
 
 callbacks = CallbackSet(summary_callback,
                         analysis_callback, alive_callback,
-                        save_solution,
                         stepsize_callback)
 
 ###############################################################################
@@ -66,8 +60,5 @@ callbacks = CallbackSet(summary_callback,
 sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
             dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
             save_everystep = false, callback = callbacks);
-summary_callback() # print the timer summary
-
-using Plots
 
-plot(sol)
\ No newline at end of file
+summary_callback() # print the timer summary
\ No newline at end of file
diff --git a/src/solvers/dgsem_tree/finite_volume_O2.jl b/src/solvers/dgsem_tree/finite_volume_O2.jl
index 27b2ce8df74..eae95edd2eb 100644
--- a/src/solvers/dgsem_tree/finite_volume_O2.jl
+++ b/src/solvers/dgsem_tree/finite_volume_O2.jl
@@ -1,8 +1,8 @@
-@inline function linear_reconstruction(u_ll, u_rr, ux_ll, ux_rr,
+@inline function linear_reconstruction(u_ll, u_rr, s_l, s_r,
                                        x_ll, x_rr, x_interfaces, node_index)
     # Linear reconstruction at the interface
-    u_ll = u_ll + ux_ll * (x_interfaces[node_index - 1] - x_ll)
-    u_rr = u_rr + ux_rr * (x_interfaces[node_index - 1] - x_rr)
+    u_ll = u_ll + s_l * (x_interfaces[node_index - 1] - x_ll)
+    u_rr = u_rr + s_r * (x_interfaces[node_index - 1] - x_rr)
 
     return u_ll, u_rr
 end
@@ -36,6 +36,14 @@ end
 # Cell index:
 #      1         2              3       4
 
+"""
+    reconstruction_small_stencil(u_mm, u_ll, u_rr, u_pp,
+                                 x_interfaces, node_index, limiter, dg)
+
+Computes limited (linear) slopes on the subcells for a DGSEM element.
+The supplied `limiter` governs the choice of slopes given the nodal values
+`u_mm`, `u_ll`, `u_rr`, and `u_pp` at the (Gauss-Lobatto Legendre) nodes.
+"""
 @inline function reconstruction_small_stencil(u_mm, u_ll, u_rr, u_pp,
                                               x_interfaces, node_index, limiter, dg)
     @unpack nodes = dg.basis
@@ -43,25 +51,25 @@ end
     x_rr = nodes[node_index]
 
     # Middle element slope
-    ux_m = (u_rr - u_ll) / (x_rr - x_ll)
+    s_mm = (u_rr - u_ll) / (x_rr - x_ll)
 
     if node_index == 2 # Catch case mm == ll
-        ux_ll = ux_m
+        s_l = s_mm
     else
         # Left element slope
-        ux_l = (u_ll - u_mm) / (x_ll - nodes[node_index - 2])
-        ux_ll = limiter.(ux_l, ux_m)
+        s_ll = (u_ll - u_mm) / (x_ll - nodes[node_index - 2])
+        s_l = limiter.(s_ll, s_mm)
     end
 
     if node_index == nnodes(dg) # Catch case rr == pp
-        ux_rr = ux_m
+        s_r = s_mm
     else
         # Right element slope
-        ux_r = (u_pp - u_rr) / (nodes[node_index + 1] - x_rr)
-        ux_rr = limiter.(ux_m, ux_r)
+        s_rr = (u_pp - u_rr) / (nodes[node_index + 1] - x_rr)
+        s_r = limiter.(s_mm, s_rr)
     end
 
-    linear_reconstruction(u_ll, u_rr, ux_ll, ux_rr, x_ll, x_rr, x_interfaces, node_index)
+    linear_reconstruction(u_ll, u_rr, s_l, s_r, x_ll, x_rr, x_interfaces, node_index)
 end
 
 """

From 86aed4fd030241b4b3f0561bb7b66f695faf16bd Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Tue, 23 Jul 2024 20:01:28 +0200
Subject: [PATCH 08/30] print some more info

---
 src/solvers/dg.jl | 1 +
 1 file changed, 1 insertion(+)

diff --git a/src/solvers/dg.jl b/src/solvers/dg.jl
index ce759c8e298..1fbfde355be 100644
--- a/src/solvers/dg.jl
+++ b/src/solvers/dg.jl
@@ -237,6 +237,7 @@ function Base.show(io::IO, ::MIME"text/plain",
             "FV flux" => integral.volume_flux_fv,
             "Reconstruction" => integral.reconstruction_mode,
             "Slope limiter" => integral.slope_limiter,
+            "Subcell boundaries" => vcat([-1.0], integral.x_interfaces, [1.0]),
         ]
         summary_box(io, "VolumeIntegralPureLGLFiniteVolumeO2", setup)
     end

From 064fb10ccbbceb73784bf6290e0862968032d13d Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Tue, 23 Jul 2024 20:15:38 +0200
Subject: [PATCH 09/30] Add unit tests

---
 test/test_unit.jl | 31 +++++++++++++++++++++++++++++++
 1 file changed, 31 insertions(+)

diff --git a/test/test_unit.jl b/test/test_unit.jl
index 09cf506576c..f6dc84de86a 100644
--- a/test/test_unit.jl
+++ b/test/test_unit.jl
@@ -1731,6 +1731,37 @@ end
     @test isapprox(mu_control(u, equations_parabolic, T_ref, R_specific, C, mu_ref),
                    1.803e-5, atol = 5e-8)
 end
+
+@testset "Limiters" begin
+    sl = 1.0
+    sr = -1.0
+
+    @test minmod(sl, sr) == 0.0
+    @test monotonized_central(sl, sr) == 0.0
+    @test superbee(sl, sr) == 0.0
+    @test vanLeer_limiter(sl, sr) == 0.0
+
+    sr = 0.5
+    @test minmod(sl, sr) == 0.5
+    @test monotonized_central(sl, sr) == 0.75
+    @test superbee(sl, sr) == 1.0
+    @test isapprox(vanLeer_limiter(sl, sr), 2/3)
+
+    sl = -1.0
+    sr = 0.0
+
+    @test minmod(sl, sr) == 0.0
+    @test monotonized_central(sl, sr) == 0.0
+    @test superbee(sl, sr) == 0.0
+    @test vanLeer_limiter(sl, sr) == 0.0
+
+    sr = -0.8
+    @test minmod(sl, sr) == -0.8
+    @test monotonized_central(sl, sr) == -0.9
+    @test superbee(sl, sr) == -1.0
+    @test isapprox(vanLeer_limiter(sl, sr), -8/9)
+end
+
 end
 
 end #module

From 8b0b72922cd2dbfc90b5f3d58854a4da16a4b91d Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Tue, 23 Jul 2024 20:20:34 +0200
Subject: [PATCH 10/30] add comment

---
 src/solvers/dgsem_tree/dg_1d.jl | 9 ++++++---
 1 file changed, 6 insertions(+), 3 deletions(-)

diff --git a/src/solvers/dgsem_tree/dg_1d.jl b/src/solvers/dgsem_tree/dg_1d.jl
index 0d52d90e297..6fe94040f9b 100644
--- a/src/solvers/dgsem_tree/dg_1d.jl
+++ b/src/solvers/dgsem_tree/dg_1d.jl
@@ -462,16 +462,19 @@ end
         #      1         2              3       4
 
         # Obtain unlimited values in primal variables
-        # TODO: If i - 2 = 0 go to neighbor element?
+
+        # Note: If i - 2 = 0 we do not go to neighbor element, as one would do in a finite volume scheme.
+        # Here, we keep it purely cell-local, thus overshoots between elements are not ruled out.
         u_mm = cons2prim(get_node_vars(u, equations, dg, max(1, i - 2), element),
                          equations)
         u_ll = cons2prim(get_node_vars(u, equations, dg, i - 1, element), equations)
         u_rr = cons2prim(get_node_vars(u, equations, dg, i, element), equations)
-        # TODO: If i + 1 > nnodes(dg) go to neighbor element?
+        # Note: If i + 1 > nnodes(dg) we do not go to neighbor element, as one would do in a finite volume scheme.
+        # Here, we keep it purely cell-local, thus overshoots between elements are not ruled out.
         u_pp = cons2prim(get_node_vars(u, equations, dg, min(nnodes(dg), i + 1),
                                        element), equations)
 
-        # Reconstruct values with limiting
+        # Reconstruct values at interfaces with limiting
         u_ll, u_rr = reconstruction_mode(u_mm, u_ll, u_rr, u_pp,
                                          x_interfaces,
                                          i, slope_limiter, dg)

From 13138eed7a1da6bf514e86e38a415eff86517cca Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Tue, 23 Jul 2024 20:22:44 +0200
Subject: [PATCH 11/30] Remove some alternative limiter implementations.

---
 src/solvers/dgsem_tree/finite_volume_O2.jl | 30 +---------------------
 1 file changed, 1 insertion(+), 29 deletions(-)

diff --git a/src/solvers/dgsem_tree/finite_volume_O2.jl b/src/solvers/dgsem_tree/finite_volume_O2.jl
index eae95edd2eb..846212f6dc4 100644
--- a/src/solvers/dgsem_tree/finite_volume_O2.jl
+++ b/src/solvers/dgsem_tree/finite_volume_O2.jl
@@ -95,16 +95,6 @@ For reference, see for instance Eq. (6.27) in
   [DOI: 10.1017/CBO9780511791253](https://doi.org/10.1017/CBO9780511791253)
 """
 @inline function minmod(sl, sr)
-    #=
-    if sign(sl) == sign(sr)
-        s = sign(sl) * min(abs(sl), abs(sr))
-    else
-        s = 0.0
-    end
-
-    return s
-    =#
-
     return 0.5 * (sign(sl) + sign(sr)) * min(abs(sl), abs(sr))
 end
 
@@ -121,14 +111,6 @@ For reference, see for instance Eq. (6.29) in
 """
 @inline function monotonized_central(sl, sr)
     # CARE: MC assumes equidistant grid in 0.5 * abs(sl + sr)!
-    #=
-    if sign(sl) == sign(sr)
-        s = sign(sl) * min(0.5 * abs(sl + sr), 2 * abs(sl), 2 * abs(sr))
-    else
-        s = 0.0
-    end
-    =#
-
     # Use recursive property of minmod function
     s = minmod(0.5 * (sl + sr), minmod(2 * sl, 2 * sr))
 
@@ -137,16 +119,6 @@ end
 
 # Note: This is NOT a limiter, just a helper for the `superbee` limiter below.
 @inline function maxmod(sl, sr)
-    #=
-    if sign(sl) == sign(sr)
-        s = sign(sl) * max(abs(sl), abs(sr))
-    else
-        s = 0.0
-    end
-
-    return s
-    =#
-
     return 0.5 * (sign(sl) + sign(sr)) * max(abs(sl), abs(sr))
 end
 
@@ -169,7 +141,7 @@ end
 """
     vanLeer_limiter(sl, sr)
 
-Symmetric Limiter by van Leer.
+Symmetric limiter by van Leer.
 See for reference page 70 in 
 
 - Siddhartha Mishra, Ulrik Skre Fjordholm and Rémi Abgrall

From 2df840db1ea3639563cd1e1efab499010ac2d7f2 Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Wed, 24 Jul 2024 21:19:59 +0200
Subject: [PATCH 12/30] move, comments, fmt

---
 .../elixir_euler_convergence_pure_fvO2.jl          | 14 +++-----------
 src/solvers/dgsem_tree/dg_1d.jl                    |  4 ++--
 src/solvers/dgsem_tree/finite_volume_O2.jl         |  4 ++--
 3 files changed, 7 insertions(+), 15 deletions(-)
 rename examples/{structured_1d_dgsem => tree_1d_dgsem}/elixir_euler_convergence_pure_fvO2.jl (84%)

diff --git a/examples/structured_1d_dgsem/elixir_euler_convergence_pure_fvO2.jl b/examples/tree_1d_dgsem/elixir_euler_convergence_pure_fvO2.jl
similarity index 84%
rename from examples/structured_1d_dgsem/elixir_euler_convergence_pure_fvO2.jl
rename to examples/tree_1d_dgsem/elixir_euler_convergence_pure_fvO2.jl
index 6aece3cf278..37620a54548 100644
--- a/examples/structured_1d_dgsem/elixir_euler_convergence_pure_fvO2.jl
+++ b/examples/tree_1d_dgsem/elixir_euler_convergence_pure_fvO2.jl
@@ -11,19 +11,11 @@ initial_condition = initial_condition_convergence_test
 
 polydeg = 3
 basis = LobattoLegendreBasis(polydeg)
-solver = DGSEM(polydeg = polydeg, surface_flux = flux_hllc,
-               volume_integral = VolumeIntegralPureLGLFiniteVolumeO2(basis, volume_flux_fv = flux_hllc,
+surf_flux = flux_hllc
+solver = DGSEM(polydeg = polydeg, surface_flux = surf_flux,
+               volume_integral = VolumeIntegralPureLGLFiniteVolumeO2(basis, volume_flux_fv = surf_flux,
                slope_limiter = superbee))
 
-#=               
-coordinates_min = (0.0,) # minimum coordinate
-coordinates_max = (2.0,) # maximum coordinate
-cells_per_dimension = (16,)
-
-# Create curved mesh with 16 cells
-mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max)
-=#
-
 coordinates_min = 0.0
 coordinates_max = 2.0
 mesh = TreeMesh(coordinates_min, coordinates_max,
diff --git a/src/solvers/dgsem_tree/dg_1d.jl b/src/solvers/dgsem_tree/dg_1d.jl
index 6fe94040f9b..863ae0eda0b 100644
--- a/src/solvers/dgsem_tree/dg_1d.jl
+++ b/src/solvers/dgsem_tree/dg_1d.jl
@@ -452,8 +452,8 @@ end
         # Gauss Lobatto Legendre nodes (schematic for k = 3):
         #   .         .                 .         .
         #   ^         ^                 ^         ^
-        # i - 2,    i - 1,               i,     i + 1
-        # mm        ll                   rr     pp
+        # i - 2,    i - 1,              i,      i + 1
+        # mm        ll                  rr      pp
         # Cell boundaries (schematic for k = 3): 
         # (note that only the inner three boundaries are stored)
         #  -1 -----------------0----------------- 1 -> x
diff --git a/src/solvers/dgsem_tree/finite_volume_O2.jl b/src/solvers/dgsem_tree/finite_volume_O2.jl
index 846212f6dc4..328925173af 100644
--- a/src/solvers/dgsem_tree/finite_volume_O2.jl
+++ b/src/solvers/dgsem_tree/finite_volume_O2.jl
@@ -27,8 +27,8 @@ end
 # Gauss Lobatto Legendre nodes (schematic for k = 3):
 #   .         .                 .         .
 #   ^         ^                 ^         ^
-# i - 2,    i - 1,               i,     i + 1
-# mm       ll                   rr     pp
+# i - 2,    i - 1,              i,      i + 1
+# mm        ll                  rr      pp
 # Cell boundaries (schematic for k = 3): 
 # (note that only the inner three boundaries are stored)
 #  -1 -----------------0----------------- 1 -> x

From 4424ee5c52f60902f26cbda1a46b8ea3e79e034b Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Wed, 24 Jul 2024 21:29:44 +0200
Subject: [PATCH 13/30] Use second order timestepping

---
 .../elixir_euler_convergence_pure_fvO2.jl             | 11 ++++++-----
 1 file changed, 6 insertions(+), 5 deletions(-)

diff --git a/examples/tree_1d_dgsem/elixir_euler_convergence_pure_fvO2.jl b/examples/tree_1d_dgsem/elixir_euler_convergence_pure_fvO2.jl
index 37620a54548..f3af646b0ce 100644
--- a/examples/tree_1d_dgsem/elixir_euler_convergence_pure_fvO2.jl
+++ b/examples/tree_1d_dgsem/elixir_euler_convergence_pure_fvO2.jl
@@ -13,8 +13,9 @@ polydeg = 3
 basis = LobattoLegendreBasis(polydeg)
 surf_flux = flux_hllc
 solver = DGSEM(polydeg = polydeg, surface_flux = surf_flux,
-               volume_integral = VolumeIntegralPureLGLFiniteVolumeO2(basis, volume_flux_fv = surf_flux,
-               slope_limiter = superbee))
+               volume_integral = VolumeIntegralPureLGLFiniteVolumeO2(basis,
+                                                                     volume_flux_fv = surf_flux,
+                                                                     slope_limiter = superbee))
 
 coordinates_min = 0.0
 coordinates_max = 2.0
@@ -40,7 +41,7 @@ analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
 
 alive_callback = AliveCallback(analysis_interval = analysis_interval)
 
-stepsize_callback = StepsizeCallback(cfl = 0.8)
+stepsize_callback = StepsizeCallback(cfl = 0.5)
 
 callbacks = CallbackSet(summary_callback,
                         analysis_callback, alive_callback,
@@ -49,8 +50,8 @@ callbacks = CallbackSet(summary_callback,
 ###############################################################################
 # run the simulation
 
-sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
+sol = solve(ode, SSPRK22(),
             dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
             save_everystep = false, callback = callbacks);
 
-summary_callback() # print the timer summary
\ No newline at end of file
+summary_callback() # print the timer summary

From f42bb9f1905699d8bc2369e595b79c4548866066 Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Wed, 24 Jul 2024 22:56:38 +0200
Subject: [PATCH 14/30] debug superbee

---
 .../elixir_euler_convergence_pure_fvO2.jl          |  6 +++---
 src/solvers/dgsem_tree/finite_volume_O2.jl         |  5 +++--
 test/test_unit.jl                                  | 14 +++++++++++++-
 3 files changed, 19 insertions(+), 6 deletions(-)

diff --git a/examples/tree_1d_dgsem/elixir_euler_convergence_pure_fvO2.jl b/examples/tree_1d_dgsem/elixir_euler_convergence_pure_fvO2.jl
index f3af646b0ce..361819dcc20 100644
--- a/examples/tree_1d_dgsem/elixir_euler_convergence_pure_fvO2.jl
+++ b/examples/tree_1d_dgsem/elixir_euler_convergence_pure_fvO2.jl
@@ -15,7 +15,7 @@ surf_flux = flux_hllc
 solver = DGSEM(polydeg = polydeg, surface_flux = surf_flux,
                volume_integral = VolumeIntegralPureLGLFiniteVolumeO2(basis,
                                                                      volume_flux_fv = surf_flux,
-                                                                     slope_limiter = superbee))
+                                                                     slope_limiter = monotonized_central))
 
 coordinates_min = 0.0
 coordinates_max = 2.0
@@ -41,7 +41,7 @@ analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
 
 alive_callback = AliveCallback(analysis_interval = analysis_interval)
 
-stepsize_callback = StepsizeCallback(cfl = 0.5)
+stepsize_callback = StepsizeCallback(cfl = 0.4)
 
 callbacks = CallbackSet(summary_callback,
                         analysis_callback, alive_callback,
@@ -54,4 +54,4 @@ sol = solve(ode, SSPRK22(),
             dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
             save_everystep = false, callback = callbacks);
 
-summary_callback() # print the timer summary
+summary_callback() # print the timer summary
\ No newline at end of file
diff --git a/src/solvers/dgsem_tree/finite_volume_O2.jl b/src/solvers/dgsem_tree/finite_volume_O2.jl
index 328925173af..96cbc408cb3 100644
--- a/src/solvers/dgsem_tree/finite_volume_O2.jl
+++ b/src/solvers/dgsem_tree/finite_volume_O2.jl
@@ -43,6 +43,7 @@ end
 Computes limited (linear) slopes on the subcells for a DGSEM element.
 The supplied `limiter` governs the choice of slopes given the nodal values
 `u_mm`, `u_ll`, `u_rr`, and `u_pp` at the (Gauss-Lobatto Legendre) nodes.
+Formally O(2) accurate.
 """
 @inline function reconstruction_small_stencil(u_mm, u_ll, u_rr, u_pp,
                                               x_interfaces, node_index, limiter, dg)
@@ -126,7 +127,7 @@ end
     superbee(sl, sr)
 
 Superbee limiter function for a TVD reconstruction given left and right slopes `sl` and `sr`.
-There are many different ways how the monotonized central limiter can be implemented.
+There are many different ways how the superbee limiter can be implemented.
 For reference, see for instance Eq. (6.28) in
 
 - Randall J. LeVeque (2002)
@@ -134,7 +135,7 @@ For reference, see for instance Eq. (6.28) in
   [DOI: 10.1017/CBO9780511791253](https://doi.org/10.1017/CBO9780511791253)
 """
 @inline function superbee(sl, sr)
-    s = maxmod(minmod(sl, 2 * sr), minmod(sl, 2 * sr))
+    s = maxmod(minmod(sl, 2 * sr), minmod(2 * sl, sr))
     return s
 end
 
diff --git a/test/test_unit.jl b/test/test_unit.jl
index f6dc84de86a..b82eb07772a 100644
--- a/test/test_unit.jl
+++ b/test/test_unit.jl
@@ -1749,7 +1749,6 @@ end
 
     sl = -1.0
     sr = 0.0
-
     @test minmod(sl, sr) == 0.0
     @test monotonized_central(sl, sr) == 0.0
     @test superbee(sl, sr) == 0.0
@@ -1760,6 +1759,19 @@ end
     @test monotonized_central(sl, sr) == -0.9
     @test superbee(sl, sr) == -1.0
     @test isapprox(vanLeer_limiter(sl, sr), -8/9)
+
+    # Test symmetry
+    @test minmod(sr, sl) == -0.8
+    @test monotonized_central(sr, sl) == -0.9
+    @test superbee(sr, sl) == -1.0
+    @test isapprox(vanLeer_limiter(sr, sl), -8/9)
+
+    sl = 1.0
+    sr = 0.0
+    @test minmod(sl, sr) == 0.0
+    @test monotonized_central(sl, sr) == 0.0
+    @test superbee(sl, sr) == 0.0
+    @test vanLeer_limiter(sl, sr) == 0.0
 end
 
 end

From 1a444991114e7f792aa53baec4de0756f733830f Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Wed, 24 Jul 2024 22:59:03 +0200
Subject: [PATCH 15/30] prim2cons 1D Adv

---
 src/equations/linear_scalar_advection_1d.jl | 2 ++
 src/equations/linear_scalar_advection_2d.jl | 2 ++
 src/equations/linear_scalar_advection_3d.jl | 2 ++
 3 files changed, 6 insertions(+)

diff --git a/src/equations/linear_scalar_advection_1d.jl b/src/equations/linear_scalar_advection_1d.jl
index 743d2df870a..bd81336e544 100644
--- a/src/equations/linear_scalar_advection_1d.jl
+++ b/src/equations/linear_scalar_advection_1d.jl
@@ -216,6 +216,8 @@ end
 
 # Convert conservative variables to primitive
 @inline cons2prim(u, equation::LinearScalarAdvectionEquation1D) = u
+# Convert primitive variables to conservative
+@inline prim2cons(u, equation::LinearScalarAdvectionEquation1D) = u
 
 # Convert conservative variables to entropy variables
 @inline cons2entropy(u, equation::LinearScalarAdvectionEquation1D) = u
diff --git a/src/equations/linear_scalar_advection_2d.jl b/src/equations/linear_scalar_advection_2d.jl
index 5e4f8463f52..19429b5425f 100644
--- a/src/equations/linear_scalar_advection_2d.jl
+++ b/src/equations/linear_scalar_advection_2d.jl
@@ -277,6 +277,8 @@ end
 
 # Convert conservative variables to primitive
 @inline cons2prim(u, equation::LinearScalarAdvectionEquation2D) = u
+# Convert primitive variables to conservative
+@inline prim2cons(u, equation::LinearScalarAdvectionEquation1D) = u
 
 # Convert conservative variables to entropy variables
 @inline cons2entropy(u, equation::LinearScalarAdvectionEquation2D) = u
diff --git a/src/equations/linear_scalar_advection_3d.jl b/src/equations/linear_scalar_advection_3d.jl
index 088f934cc3e..2a530800959 100644
--- a/src/equations/linear_scalar_advection_3d.jl
+++ b/src/equations/linear_scalar_advection_3d.jl
@@ -196,6 +196,8 @@ end
 
 # Convert conservative variables to primitive
 @inline cons2prim(u, equation::LinearScalarAdvectionEquation3D) = u
+# Convert primitive variables to conservative
+@inline prim2cons(u, equation::LinearScalarAdvectionEquation1D) = u
 
 # Convert conservative variables to entropy variables
 @inline cons2entropy(u, equation::LinearScalarAdvectionEquation3D) = u

From 5bb0ab2d915627af9bbaec361a685d2e6f8d37f0 Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Wed, 24 Jul 2024 23:04:23 +0200
Subject: [PATCH 16/30] test

---
 test/test_tree_1d_euler.jl | 22 ++++++++++++++++++++++
 test/test_unit.jl          |  7 +++----
 2 files changed, 25 insertions(+), 4 deletions(-)

diff --git a/test/test_tree_1d_euler.jl b/test/test_tree_1d_euler.jl
index dc523586f89..31f285fa3a0 100644
--- a/test/test_tree_1d_euler.jl
+++ b/test/test_tree_1d_euler.jl
@@ -490,6 +490,28 @@ end
         @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
     end
 end
+
+@trixi_testset "elixir_euler_convergence_pure_fvO2.jl" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_convergence_pure_fvO2.jl"),
+                        l2=[
+                            0.00048190103584221393,
+                            0.0005306536705484938,
+                            0.0007902755183354382,
+                        ],
+                        linf=[
+                            0.0014512544822959939,
+                            0.0014178023249342697,
+                            0.002568866462516972,
+                        ])
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
+    let
+        t = sol.t[end]
+        u_ode = sol.u[end]
+        du_ode = similar(u_ode)
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+    end
+end
 end
 
 end # module
diff --git a/test/test_unit.jl b/test/test_unit.jl
index b82eb07772a..76a9acf89e3 100644
--- a/test/test_unit.jl
+++ b/test/test_unit.jl
@@ -1745,7 +1745,7 @@ end
     @test minmod(sl, sr) == 0.5
     @test monotonized_central(sl, sr) == 0.75
     @test superbee(sl, sr) == 1.0
-    @test isapprox(vanLeer_limiter(sl, sr), 2/3)
+    @test isapprox(vanLeer_limiter(sl, sr), 2 / 3)
 
     sl = -1.0
     sr = 0.0
@@ -1758,13 +1758,13 @@ end
     @test minmod(sl, sr) == -0.8
     @test monotonized_central(sl, sr) == -0.9
     @test superbee(sl, sr) == -1.0
-    @test isapprox(vanLeer_limiter(sl, sr), -8/9)
+    @test isapprox(vanLeer_limiter(sl, sr), -8 / 9)
 
     # Test symmetry
     @test minmod(sr, sl) == -0.8
     @test monotonized_central(sr, sl) == -0.9
     @test superbee(sr, sl) == -1.0
-    @test isapprox(vanLeer_limiter(sr, sl), -8/9)
+    @test isapprox(vanLeer_limiter(sr, sl), -8 / 9)
 
     sl = 1.0
     sr = 0.0
@@ -1773,7 +1773,6 @@ end
     @test superbee(sl, sr) == 0.0
     @test vanLeer_limiter(sl, sr) == 0.0
 end
-
 end
 
 end #module

From c2aab4ee4f965772dad7d1f51ebde46e39029aca Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Thu, 25 Jul 2024 00:00:28 +0200
Subject: [PATCH 17/30] fmt, typo

---
 examples/tree_1d_dgsem/elixir_euler_convergence_pure_fvO2.jl | 2 +-
 src/solvers/dgsem_tree/finite_volume_O2.jl                   | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/examples/tree_1d_dgsem/elixir_euler_convergence_pure_fvO2.jl b/examples/tree_1d_dgsem/elixir_euler_convergence_pure_fvO2.jl
index 361819dcc20..7961bf50769 100644
--- a/examples/tree_1d_dgsem/elixir_euler_convergence_pure_fvO2.jl
+++ b/examples/tree_1d_dgsem/elixir_euler_convergence_pure_fvO2.jl
@@ -54,4 +54,4 @@ sol = solve(ode, SSPRK22(),
             dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
             save_everystep = false, callback = callbacks);
 
-summary_callback() # print the timer summary
\ No newline at end of file
+summary_callback() # print the timer summary
diff --git a/src/solvers/dgsem_tree/finite_volume_O2.jl b/src/solvers/dgsem_tree/finite_volume_O2.jl
index 96cbc408cb3..fcf37480b31 100644
--- a/src/solvers/dgsem_tree/finite_volume_O2.jl
+++ b/src/solvers/dgsem_tree/finite_volume_O2.jl
@@ -77,7 +77,7 @@ end
     central_recon(sl, sr)
 
 Central, non-TVD reconstruction given left and right slopes `sl` and `sr`.
-Gives formally full order of accuracy at the expense of sacrificied nonlinear stability.
+Gives formally full order of accuracy at the expense of sacrificed nonlinear stability.
 """
 @inline function central_recon(sl, sr)
     s = 0.5 * (sl + sr)

From 3d591c535c978cbb19d2e22f7e79d5eb5df4ff17 Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Thu, 25 Jul 2024 00:11:05 +0200
Subject: [PATCH 18/30] typos

---
 src/equations/linear_scalar_advection_2d.jl | 2 +-
 src/equations/linear_scalar_advection_3d.jl | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/equations/linear_scalar_advection_2d.jl b/src/equations/linear_scalar_advection_2d.jl
index 19429b5425f..704b598779d 100644
--- a/src/equations/linear_scalar_advection_2d.jl
+++ b/src/equations/linear_scalar_advection_2d.jl
@@ -278,7 +278,7 @@ end
 # Convert conservative variables to primitive
 @inline cons2prim(u, equation::LinearScalarAdvectionEquation2D) = u
 # Convert primitive variables to conservative
-@inline prim2cons(u, equation::LinearScalarAdvectionEquation1D) = u
+@inline prim2cons(u, equation::LinearScalarAdvectionEquation2D) = u
 
 # Convert conservative variables to entropy variables
 @inline cons2entropy(u, equation::LinearScalarAdvectionEquation2D) = u
diff --git a/src/equations/linear_scalar_advection_3d.jl b/src/equations/linear_scalar_advection_3d.jl
index 2a530800959..defe1dd2d08 100644
--- a/src/equations/linear_scalar_advection_3d.jl
+++ b/src/equations/linear_scalar_advection_3d.jl
@@ -197,7 +197,7 @@ end
 # Convert conservative variables to primitive
 @inline cons2prim(u, equation::LinearScalarAdvectionEquation3D) = u
 # Convert primitive variables to conservative
-@inline prim2cons(u, equation::LinearScalarAdvectionEquation1D) = u
+@inline prim2cons(u, equation::LinearScalarAdvectionEquation3D) = u
 
 # Convert conservative variables to entropy variables
 @inline cons2entropy(u, equation::LinearScalarAdvectionEquation3D) = u

From fdc590ab4e715f798768048aa591e344d4bb00e8 Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Thu, 25 Jul 2024 06:36:05 +0200
Subject: [PATCH 19/30] some more tests

---
 src/Trixi.jl                               |  5 +++--
 src/solvers/dgsem_tree/finite_volume_O2.jl | 11 ++++++-----
 test/test_unit.jl                          |  2 ++
 3 files changed, 11 insertions(+), 7 deletions(-)

diff --git a/src/Trixi.jl b/src/Trixi.jl
index f5c237a9044..386ceec6312 100644
--- a/src/Trixi.jl
+++ b/src/Trixi.jl
@@ -241,8 +241,9 @@ export DG,
        SurfaceIntegralUpwind,
        MortarL2
 
-export reconstruction_small_stencil,
-       minmod, monotonized_central, superbee, vanLeer_limiter
+export reconstruction_small_stencil, reconstruction_constant,
+       minmod, monotonized_central, superbee, vanLeer_limiter,
+       central_slope
 
 export VolumeIntegralSubcellLimiting, BoundsCheckCallback,
        SubcellLimiterIDP, SubcellLimiterIDPCorrection
diff --git a/src/solvers/dgsem_tree/finite_volume_O2.jl b/src/solvers/dgsem_tree/finite_volume_O2.jl
index fcf37480b31..1cc97a7f721 100644
--- a/src/solvers/dgsem_tree/finite_volume_O2.jl
+++ b/src/solvers/dgsem_tree/finite_volume_O2.jl
@@ -17,8 +17,8 @@ obtained from constant polynomials.
 Formally O(1) accurate.
 """
 @inline function reconstruction_constant(u_mm, u_ll, u_rr, u_pp,
-                                         x_interfaces,
-                                         node_index, limiter, dg)
+                                         x_interfaces, node_index, 
+                                         limiter, dg)
     return u_ll, u_rr
 end
 
@@ -46,7 +46,8 @@ The supplied `limiter` governs the choice of slopes given the nodal values
 Formally O(2) accurate.
 """
 @inline function reconstruction_small_stencil(u_mm, u_ll, u_rr, u_pp,
-                                              x_interfaces, node_index, limiter, dg)
+                                              x_interfaces, node_index, 
+                                              limiter, dg)
     @unpack nodes = dg.basis
     x_ll = nodes[node_index - 1]
     x_rr = nodes[node_index]
@@ -74,12 +75,12 @@ Formally O(2) accurate.
 end
 
 """
-    central_recon(sl, sr)
+    central_slope(sl, sr)
 
 Central, non-TVD reconstruction given left and right slopes `sl` and `sr`.
 Gives formally full order of accuracy at the expense of sacrificed nonlinear stability.
 """
-@inline function central_recon(sl, sr)
+@inline function central_slope(sl, sr)
     s = 0.5 * (sl + sr)
     return s
 end
diff --git a/test/test_unit.jl b/test/test_unit.jl
index c38433987c2..3959ce7666b 100644
--- a/test/test_unit.jl
+++ b/test/test_unit.jl
@@ -1778,6 +1778,8 @@ end
     @test monotonized_central(sl, sr) == 0.0
     @test superbee(sl, sr) == 0.0
     @test vanLeer_limiter(sl, sr) == 0.0
+
+    @test central_slope(sl, sr) == 0.5
 end
 end
 

From 8aabfd50b183652413810f9768cdb349b71a7693 Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Thu, 25 Jul 2024 06:38:26 +0200
Subject: [PATCH 20/30] fmt

---
 src/solvers/dgsem_tree/finite_volume_O2.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/solvers/dgsem_tree/finite_volume_O2.jl b/src/solvers/dgsem_tree/finite_volume_O2.jl
index 1cc97a7f721..753735c2a39 100644
--- a/src/solvers/dgsem_tree/finite_volume_O2.jl
+++ b/src/solvers/dgsem_tree/finite_volume_O2.jl
@@ -17,7 +17,7 @@ obtained from constant polynomials.
 Formally O(1) accurate.
 """
 @inline function reconstruction_constant(u_mm, u_ll, u_rr, u_pp,
-                                         x_interfaces, node_index, 
+                                         x_interfaces, node_index,
                                          limiter, dg)
     return u_ll, u_rr
 end
@@ -46,7 +46,7 @@ The supplied `limiter` governs the choice of slopes given the nodal values
 Formally O(2) accurate.
 """
 @inline function reconstruction_small_stencil(u_mm, u_ll, u_rr, u_pp,
-                                              x_interfaces, node_index, 
+                                              x_interfaces, node_index,
                                               limiter, dg)
     @unpack nodes = dg.basis
     x_ll = nodes[node_index - 1]

From 072ca29a8982242c885bb3a01145f5971bfb0d31 Mon Sep 17 00:00:00 2001
From: Daniel Doehring <daniel.doehring@rwth-aachen.de>
Date: Sat, 10 Aug 2024 07:12:48 +0200
Subject: [PATCH 21/30] Update src/solvers/dgsem_tree/finite_volume_O2.jl

---
 src/solvers/dgsem_tree/finite_volume_O2.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/solvers/dgsem_tree/finite_volume_O2.jl b/src/solvers/dgsem_tree/finite_volume_O2.jl
index 753735c2a39..119f86c378d 100644
--- a/src/solvers/dgsem_tree/finite_volume_O2.jl
+++ b/src/solvers/dgsem_tree/finite_volume_O2.jl
@@ -24,7 +24,7 @@ end
 
 #             Reference element:             
 #  -1 -----------------0----------------- 1 -> x
-# Gauss Lobatto Legendre nodes (schematic for k = 3):
+# Gauss-Lobatto-Legendre nodes (schematic for k = 3):
 #   .         .                 .         .
 #   ^         ^                 ^         ^
 # i - 2,    i - 1,              i,      i + 1

From fa18a4741efe5e804a3daed4b3d96bc2452fbdcb Mon Sep 17 00:00:00 2001
From: Daniel Doehring <daniel.doehring@rwth-aachen.de>
Date: Sat, 10 Aug 2024 07:12:53 +0200
Subject: [PATCH 22/30] Update test/test_unit.jl

---
 test/test_unit.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/test_unit.jl b/test/test_unit.jl
index 3959ce7666b..ead32694cf2 100644
--- a/test/test_unit.jl
+++ b/test/test_unit.jl
@@ -1738,7 +1738,7 @@ end
                    1.803e-5, atol = 5e-8)
 end
 
-@testset "Limiters" begin
+@testset "Slope Limiters" begin
     sl = 1.0
     sr = -1.0
 

From a64bb60a970a478c341c705f0176cb69c0da3849 Mon Sep 17 00:00:00 2001
From: Daniel Doehring <daniel.doehring@rwth-aachen.de>
Date: Sat, 10 Aug 2024 07:13:01 +0200
Subject: [PATCH 23/30] Update src/solvers/dgsem_tree/dg_1d.jl

---
 src/solvers/dgsem_tree/dg_1d.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/solvers/dgsem_tree/dg_1d.jl b/src/solvers/dgsem_tree/dg_1d.jl
index 725dd46848a..2c55a8fc22f 100644
--- a/src/solvers/dgsem_tree/dg_1d.jl
+++ b/src/solvers/dgsem_tree/dg_1d.jl
@@ -449,7 +449,7 @@ end
     for i in 2:nnodes(dg)
         #             Reference element:             
         #  -1 -----------------0----------------- 1 -> x
-        # Gauss Lobatto Legendre nodes (schematic for k = 3):
+        # Gauss-Lobatto-Legendre nodes (schematic for k = 3):
         #   .         .                 .         .
         #   ^         ^                 ^         ^
         # i - 2,    i - 1,              i,      i + 1

From 610dc2487fb0bfb6000cfb14630837cd2d694fa7 Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Tue, 24 Sep 2024 17:20:26 +0200
Subject: [PATCH 24/30] fmt

---
 src/solvers/dg.jl          | 2 +-
 test/test_tree_1d_euler.jl | 4 ++--
 2 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/solvers/dg.jl b/src/solvers/dg.jl
index 11441ab3ac7..cf3dd08f251 100644
--- a/src/solvers/dg.jl
+++ b/src/solvers/dg.jl
@@ -237,7 +237,7 @@ function Base.show(io::IO, ::MIME"text/plain",
             "FV flux" => integral.volume_flux_fv,
             "Reconstruction" => integral.reconstruction_mode,
             "Slope limiter" => integral.slope_limiter,
-            "Subcell boundaries" => vcat([-1.0], integral.x_interfaces, [1.0]),
+            "Subcell boundaries" => vcat([-1.0], integral.x_interfaces, [1.0])
         ]
         summary_box(io, "VolumeIntegralPureLGLFiniteVolumeO2", setup)
     end
diff --git a/test/test_tree_1d_euler.jl b/test/test_tree_1d_euler.jl
index cf01545e0c7..68096329e2d 100644
--- a/test/test_tree_1d_euler.jl
+++ b/test/test_tree_1d_euler.jl
@@ -496,12 +496,12 @@ end
                         l2=[
                             0.00048190103584221393,
                             0.0005306536705484938,
-                            0.0007902755183354382,
+                            0.0007902755183354382
                         ],
                         linf=[
                             0.0014512544822959939,
                             0.0014178023249342697,
-                            0.002568866462516972,
+                            0.002568866462516972
                         ])
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)

From 699d7e5310f1ca781bffbfd7abca1b5e5e54b2fb Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Mon, 7 Oct 2024 12:23:39 +0200
Subject: [PATCH 25/30] add different recontruction mode

---
 ...xir_euler_source_terms_nonperiodic_fvO2.jl | 66 +++++++++++++++++
 .../elixir_euler_convergence_pure_fvO2.jl     |  6 +-
 src/Trixi.jl                                  |  3 +-
 src/solvers/dg.jl                             | 26 +++++--
 src/solvers/dgsem_tree/dg_1d.jl               |  6 +-
 src/solvers/dgsem_tree/finite_volume_O2.jl    | 72 +++++++++++++++----
 6 files changed, 153 insertions(+), 26 deletions(-)
 create mode 100644 examples/structured_1d_dgsem/elixir_euler_source_terms_nonperiodic_fvO2.jl

diff --git a/examples/structured_1d_dgsem/elixir_euler_source_terms_nonperiodic_fvO2.jl b/examples/structured_1d_dgsem/elixir_euler_source_terms_nonperiodic_fvO2.jl
new file mode 100644
index 00000000000..a1a9841eb0f
--- /dev/null
+++ b/examples/structured_1d_dgsem/elixir_euler_source_terms_nonperiodic_fvO2.jl
@@ -0,0 +1,66 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+
+equations = CompressibleEulerEquations1D(1.4)
+
+initial_condition = initial_condition_convergence_test
+
+source_terms = source_terms_convergence_test
+
+# you can either use a single function to impose the BCs weakly in all
+# 1*ndims == 2 directions or you can pass a tuple containing BCs for
+# each direction
+boundary_condition = BoundaryConditionDirichlet(initial_condition)
+boundary_conditions = (x_neg = boundary_condition,
+                       x_pos = boundary_condition)
+
+polydeg = 8 # Governs in this case only the number of subcells
+basis = LobattoLegendreBasis(polydeg)
+surf_flux = flux_hll
+solver = DGSEM(polydeg = polydeg, surface_flux = surf_flux,
+               volume_integral = VolumeIntegralPureLGLFiniteVolumeO2(basis,
+                                                                     volume_flux_fv = surf_flux,
+                                                                     reconstruction_mode = reconstruction_small_stencil_inner,
+                                                                     slope_limiter = monotonized_central))
+
+f1() = SVector(0.0)
+f2() = SVector(2.0)
+initial_refinement_level = 3 # required for convergence study
+cells_per_dimension = 2^initial_refinement_level
+mesh = StructuredMesh((cells_per_dimension,), (f1, f2), periodicity = false)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    source_terms = source_terms,
+                                    boundary_conditions = boundary_conditions)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 2.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
+                                     extra_analysis_errors = (:conservation_error,))
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+stepsize_callback = StepsizeCallback(cfl = 0.3)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback,
+                        stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, ORK256(),
+            dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep = false, callback = callbacks);
+summary_callback() # print the timer summary
diff --git a/examples/tree_1d_dgsem/elixir_euler_convergence_pure_fvO2.jl b/examples/tree_1d_dgsem/elixir_euler_convergence_pure_fvO2.jl
index 7961bf50769..df8bdae970d 100644
--- a/examples/tree_1d_dgsem/elixir_euler_convergence_pure_fvO2.jl
+++ b/examples/tree_1d_dgsem/elixir_euler_convergence_pure_fvO2.jl
@@ -9,12 +9,13 @@ equations = CompressibleEulerEquations1D(1.4)
 
 initial_condition = initial_condition_convergence_test
 
-polydeg = 3
+polydeg = 3 # Governs in this case only the number of subcells
 basis = LobattoLegendreBasis(polydeg)
 surf_flux = flux_hllc
 solver = DGSEM(polydeg = polydeg, surface_flux = surf_flux,
                volume_integral = VolumeIntegralPureLGLFiniteVolumeO2(basis,
                                                                      volume_flux_fv = surf_flux,
+                                                                     reconstruction_mode = reconstruction_small_stencil_full,
                                                                      slope_limiter = monotonized_central))
 
 coordinates_min = 0.0
@@ -37,7 +38,8 @@ summary_callback = SummaryCallback()
 analysis_interval = 100
 analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
                                      extra_analysis_errors = (:l2_error_primitive,
-                                                              :linf_error_primitive))
+                                                              :linf_error_primitive,
+                                                              :conservation_error))
 
 alive_callback = AliveCallback(analysis_interval = analysis_interval)
 
diff --git a/src/Trixi.jl b/src/Trixi.jl
index c480473bab4..5d53a58fd1f 100644
--- a/src/Trixi.jl
+++ b/src/Trixi.jl
@@ -242,7 +242,8 @@ export DG,
        SurfaceIntegralUpwind,
        MortarL2
 
-export reconstruction_small_stencil, reconstruction_constant,
+export reconstruction_small_stencil_inner, reconstruction_small_stencil_full, 
+       reconstruction_constant,
        minmod, monotonized_central, superbee, vanLeer_limiter,
        central_slope
 
diff --git a/src/solvers/dg.jl b/src/solvers/dg.jl
index cf3dd08f251..07765eeb9a0 100644
--- a/src/solvers/dg.jl
+++ b/src/solvers/dg.jl
@@ -187,20 +187,36 @@ end
 """
     VolumeIntegralPureLGLFiniteVolumeO2(basis::Basis;
                                         volume_flux_fv = flux_lax_friedrichs,
-                                        reconstruction_mode = reconstruction_small_stencil,
+                                        reconstruction_mode = reconstruction_small_stencil_inner,
                                         slope_limiter = minmod)
 
-This gives a formally O(2)-accurate finite volume scheme on an LGL-type subcell
+This gives an up to  O(2)-accurate finite volume scheme on an LGL-type subcell
 mesh (LGL = Legendre-Gauss-Lobatto).
+Depending on the `reconstruction_mode` and `slope_limiter`, experimental orders of convergence
+between 1 and 2 can be expected in practice.
+Currently, all reconstructions are purely cell-local, i.e., no neighboring elements are 
+queried at reconstruction stage.
+
+The non-boundary subcells are always reconstructed using the standard MUSCL-type reconstruction.
+For the subcells at the boundaries, two options are available:
+
+1) The unlimited slope is used on these cells. This gives full O(2) accuracy, but may lead to overshoots between cells.
+   The `reconstruction_mode` corresponding to this is `reconstruction_small_stencil_full`.
+2) On boundary subcells, the solution is represented using a constant value, thereby falling back to formally only O(1).
+   The `reconstruction_mode` corresponding to this is `reconstruction_small_stencil_inner`.
+   In the reference below, this is the recommended reconstruction mode and is thus used by default.
 
 !!! warning "Experimental implementation"
     This is an experimental feature and may change in future releases.
 
 ## References
 
-- Hennemann, Gassner (2020)
-  "A provably entropy stable subcell shock capturing approach for high order split form DG"
-  [arXiv: 2008.12044](https://arxiv.org/abs/2008.12044)
+See especially Sections 3.2 and 4 and Appendix D of the paper
+
+- Rueda-Ramírez, Hennemann, Hindenlang, Winters, & Gassner (2021).
+  "An entropy stable nodal discontinuous Galerkin method for the resistive MHD equations. 
+   Part II: Subcell finite volume shock capturing"
+  [JCP: 2021.110580](https://doi.org/10.1016/j.jcp.2021.110580)
 """
 struct VolumeIntegralPureLGLFiniteVolumeO2{RealT, Basis, VolumeFluxFV, Reconstruction,
                                            Limiter} <: AbstractVolumeIntegral
diff --git a/src/solvers/dgsem_tree/dg_1d.jl b/src/solvers/dgsem_tree/dg_1d.jl
index 681b396eedd..2b2e6cb5b11 100644
--- a/src/solvers/dgsem_tree/dg_1d.jl
+++ b/src/solvers/dgsem_tree/dg_1d.jl
@@ -461,7 +461,7 @@ end
         # Cell index:
         #      1         2              3       4
 
-        # Obtain unlimited values in primal variables
+        ## Obtain unlimited values in primal variables ##
 
         # Note: If i - 2 = 0 we do not go to neighbor element, as one would do in a finite volume scheme.
         # Here, we keep it purely cell-local, thus overshoots between elements are not ruled out.
@@ -474,12 +474,12 @@ end
         u_pp = cons2prim(get_node_vars(u, equations, dg, min(nnodes(dg), i + 1),
                                        element), equations)
 
-        # Reconstruct values at interfaces with limiting
+        ## Reconstruct values at interfaces with limiting ##
         u_ll, u_rr = reconstruction_mode(u_mm, u_ll, u_rr, u_pp,
                                          x_interfaces,
                                          i, slope_limiter, dg)
 
-        # Convert primal variables back to conservative variables
+        ## Convert primal variables back to conservative variables ##
         flux = volume_flux_fv(prim2cons(u_ll, equations), prim2cons(u_rr, equations),
                               1, equations) # orientation 1: x direction
 
diff --git a/src/solvers/dgsem_tree/finite_volume_O2.jl b/src/solvers/dgsem_tree/finite_volume_O2.jl
index 119f86c378d..0ee8a01b5ad 100644
--- a/src/solvers/dgsem_tree/finite_volume_O2.jl
+++ b/src/solvers/dgsem_tree/finite_volume_O2.jl
@@ -1,12 +1,3 @@
-@inline function linear_reconstruction(u_ll, u_rr, s_l, s_r,
-                                       x_ll, x_rr, x_interfaces, node_index)
-    # Linear reconstruction at the interface
-    u_ll = u_ll + s_l * (x_interfaces[node_index - 1] - x_ll)
-    u_rr = u_rr + s_r * (x_interfaces[node_index - 1] - x_rr)
-
-    return u_ll, u_rr
-end
-
 """
     reconstruction_constant(u_mm, u_ll, u_rr, u_pp,
                             x_interfaces,
@@ -22,6 +13,15 @@ Formally O(1) accurate.
     return u_ll, u_rr
 end
 
+@inline function linear_reconstruction(u_ll, u_rr, s_l, s_r,
+                                       x_ll, x_rr, x_interfaces, node_index)
+    # Linear reconstruction at the interface
+    u_ll = u_ll + s_l * (x_interfaces[node_index - 1] - x_ll)
+    u_rr = u_rr + s_r * (x_interfaces[node_index - 1] - x_rr)
+
+    return u_ll, u_rr
+end
+
 #             Reference element:             
 #  -1 -----------------0----------------- 1 -> x
 # Gauss-Lobatto-Legendre nodes (schematic for k = 3):
@@ -37,17 +37,18 @@ end
 #      1         2              3       4
 
 """
-    reconstruction_small_stencil(u_mm, u_ll, u_rr, u_pp,
-                                 x_interfaces, node_index, limiter, dg)
+    reconstruction_small_stencil_full(u_mm, u_ll, u_rr, u_pp,
+                                      x_interfaces, node_index, limiter, dg)
 
 Computes limited (linear) slopes on the subcells for a DGSEM element.
 The supplied `limiter` governs the choice of slopes given the nodal values
 `u_mm`, `u_ll`, `u_rr`, and `u_pp` at the (Gauss-Lobatto Legendre) nodes.
-Formally O(2) accurate.
+Formally O(2) accurate when used without a limiter, i.e., the `central_slope`.
+For the subcell boundaries the reconstructed slopes are not limited.
 """
-@inline function reconstruction_small_stencil(u_mm, u_ll, u_rr, u_pp,
-                                              x_interfaces, node_index,
-                                              limiter, dg)
+@inline function reconstruction_small_stencil_full(u_mm, u_ll, u_rr, u_pp,
+                                                   x_interfaces, node_index,
+                                                   limiter, dg)
     @unpack nodes = dg.basis
     x_ll = nodes[node_index - 1]
     x_rr = nodes[node_index]
@@ -74,6 +75,47 @@ Formally O(2) accurate.
     linear_reconstruction(u_ll, u_rr, s_l, s_r, x_ll, x_rr, x_interfaces, node_index)
 end
 
+"""
+    reconstruction_small_stencil_inner(u_mm, u_ll, u_rr, u_pp,
+                                       x_interfaces, node_index, limiter, dg)
+
+Computes limited (linear) slopes on the *inner* subcells for a DGSEM element.
+The supplied `limiter` governs the choice of slopes given the nodal values
+`u_mm`, `u_ll`, `u_rr`, and `u_pp` at the (Gauss-Lobatto Legendre) nodes.
+For the outer, i.e., boundary subcells, constant values are used, i.e, no reconstruction.
+This reduces the order of the scheme below 2.
+"""
+@inline function reconstruction_small_stencil_inner(u_mm, u_ll, u_rr, u_pp,
+                                                    x_interfaces, node_index,
+                                                    limiter, dg)
+    @unpack nodes = dg.basis
+    x_ll = nodes[node_index - 1]
+    x_rr = nodes[node_index]
+
+    # Middle element slope
+    s_mm = (u_rr - u_ll) / (x_rr - x_ll)
+
+    if node_index == 2 # Catch case mm == ll
+        # Do not reconstruct at the boundary
+        s_l = zero(s_mm)
+    else
+        # Left element slope
+        s_ll = (u_ll - u_mm) / (x_ll - nodes[node_index - 2])
+        s_l = limiter.(s_ll, s_mm)
+    end
+
+    if node_index == nnodes(dg) # Catch case rr == pp
+        # Do not reconstruct at the boundary
+        s_r = zero(s_mm)
+    else
+        # Right element slope
+        s_rr = (u_pp - u_rr) / (nodes[node_index + 1] - x_rr)
+        s_r = limiter.(s_mm, s_rr)
+    end
+
+    linear_reconstruction(u_ll, u_rr, s_l, s_r, x_ll, x_rr, x_interfaces, node_index)
+end
+
 """
     central_slope(sl, sr)
 

From 07733318bc956f9154679ea3d2d342c9d6ac92c5 Mon Sep 17 00:00:00 2001
From: Daniel Doehring <daniel.doehring@rwth-aachen.de>
Date: Mon, 7 Oct 2024 12:24:40 +0200
Subject: [PATCH 26/30] Update src/solvers/dgsem_tree/finite_volume_O2.jl
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Co-authored-by: Andrés Rueda-Ramírez <aruedara@uni-koeln.de>
---
 src/solvers/dgsem_tree/finite_volume_O2.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/solvers/dgsem_tree/finite_volume_O2.jl b/src/solvers/dgsem_tree/finite_volume_O2.jl
index 0ee8a01b5ad..efb27d4face 100644
--- a/src/solvers/dgsem_tree/finite_volume_O2.jl
+++ b/src/solvers/dgsem_tree/finite_volume_O2.jl
@@ -154,7 +154,7 @@ For reference, see for instance Eq. (6.29) in
   [DOI: 10.1017/CBO9780511791253](https://doi.org/10.1017/CBO9780511791253)
 """
 @inline function monotonized_central(sl, sr)
-    # CARE: MC assumes equidistant grid in 0.5 * abs(sl + sr)!
+    # CARE: MC assumes equidistant grid in 0.5 * (sl + sr)!
     # Use recursive property of minmod function
     s = minmod(0.5 * (sl + sr), minmod(2 * sl, 2 * sr))
 

From 9f71e16cf65162874b73984a47d806d0f76e7124 Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Mon, 7 Oct 2024 12:29:59 +0200
Subject: [PATCH 27/30] test + fmt

---
 src/Trixi.jl               |  2 +-
 test/test_structured_1d.jl | 23 +++++++++++++++++++++++
 2 files changed, 24 insertions(+), 1 deletion(-)

diff --git a/src/Trixi.jl b/src/Trixi.jl
index 5d53a58fd1f..8c0331c1892 100644
--- a/src/Trixi.jl
+++ b/src/Trixi.jl
@@ -242,7 +242,7 @@ export DG,
        SurfaceIntegralUpwind,
        MortarL2
 
-export reconstruction_small_stencil_inner, reconstruction_small_stencil_full, 
+export reconstruction_small_stencil_inner, reconstruction_small_stencil_full,
        reconstruction_constant,
        minmod, monotonized_central, superbee, vanLeer_limiter,
        central_slope
diff --git a/test/test_structured_1d.jl b/test/test_structured_1d.jl
index f64b8c9c065..b0ad7a5847f 100644
--- a/test/test_structured_1d.jl
+++ b/test/test_structured_1d.jl
@@ -139,6 +139,29 @@ end
     end
 end
 
+@trixi_testset "elixir_euler_source_terms_nonperiodic_fvO2.jl" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR,
+                                 "elixir_euler_source_terms_nonperiodic_fvO2.jl"),
+                        l2=[
+                            0.0003005552269336209,
+                            0.0004392132979796131,
+                            0.0008693370944987093
+                        ],
+                        linf=[
+                            0.0013999244590574556,
+                            0.001298120994041474,
+                            0.004533951299385386
+                        ])
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
+    let
+        t = sol.t[end]
+        u_ode = sol.u[end]
+        du_ode = similar(u_ode)
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+    end
+end
+
 @trixi_testset "elixir_linearizedeuler_characteristic_system.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR,
                                  "elixir_linearizedeuler_characteristic_system.jl"),

From 1f9bb57b8542187c8cee073396aa63f7381951f3 Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Mon, 7 Oct 2024 12:32:16 +0200
Subject: [PATCH 28/30] comments

---
 src/solvers/dgsem_tree/dg_1d.jl | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/src/solvers/dgsem_tree/dg_1d.jl b/src/solvers/dgsem_tree/dg_1d.jl
index 2b2e6cb5b11..850547d1c4f 100644
--- a/src/solvers/dgsem_tree/dg_1d.jl
+++ b/src/solvers/dgsem_tree/dg_1d.jl
@@ -464,13 +464,15 @@ end
         ## Obtain unlimited values in primal variables ##
 
         # Note: If i - 2 = 0 we do not go to neighbor element, as one would do in a finite volume scheme.
-        # Here, we keep it purely cell-local, thus overshoots between elements are not ruled out.
+        # Here, we keep it purely cell-local, thus overshoots between elements are not ruled out, 
+        # unless `reconstruction_small_stencil_inner` is used as the `reconstruction_mode`.
         u_mm = cons2prim(get_node_vars(u, equations, dg, max(1, i - 2), element),
                          equations)
         u_ll = cons2prim(get_node_vars(u, equations, dg, i - 1, element), equations)
         u_rr = cons2prim(get_node_vars(u, equations, dg, i, element), equations)
         # Note: If i + 1 > nnodes(dg) we do not go to neighbor element, as one would do in a finite volume scheme.
-        # Here, we keep it purely cell-local, thus overshoots between elements are not ruled out.
+        # Here, we keep it purely cell-local, thus overshoots between elements are not ruled out,
+        # unless `reconstruction_small_stencil_inner` is used as the `reconstruction_mode`.
         u_pp = cons2prim(get_node_vars(u, equations, dg, min(nnodes(dg), i + 1),
                                        element), equations)
 

From 20fe628af9f03684d948eead38edec3af4f7baa6 Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Mon, 7 Oct 2024 12:42:20 +0200
Subject: [PATCH 29/30] correct way cells dim

---
 .../elixir_euler_source_terms_nonperiodic_fvO2.jl            | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/examples/structured_1d_dgsem/elixir_euler_source_terms_nonperiodic_fvO2.jl b/examples/structured_1d_dgsem/elixir_euler_source_terms_nonperiodic_fvO2.jl
index a1a9841eb0f..3c1b9917a29 100644
--- a/examples/structured_1d_dgsem/elixir_euler_source_terms_nonperiodic_fvO2.jl
+++ b/examples/structured_1d_dgsem/elixir_euler_source_terms_nonperiodic_fvO2.jl
@@ -29,9 +29,8 @@ solver = DGSEM(polydeg = polydeg, surface_flux = surf_flux,
 
 f1() = SVector(0.0)
 f2() = SVector(2.0)
-initial_refinement_level = 3 # required for convergence study
-cells_per_dimension = 2^initial_refinement_level
-mesh = StructuredMesh((cells_per_dimension,), (f1, f2), periodicity = false)
+cells_per_dimension = (8,)
+mesh = StructuredMesh(cells_per_dimension, (f1, f2), periodicity = false)
 
 semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
                                     source_terms = source_terms,

From ba708b9c8013dca9a52d4268fed10ed45fc929d3 Mon Sep 17 00:00:00 2001
From: Daniel_Doehring <daniel.doehring@rwth-aachen.de>
Date: Mon, 7 Oct 2024 14:05:49 +0200
Subject: [PATCH 30/30] increase coverage

---
 test/test_unit.jl | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/test/test_unit.jl b/test/test_unit.jl
index 42f42693d89..7ba84262ce4 100644
--- a/test/test_unit.jl
+++ b/test/test_unit.jl
@@ -1742,6 +1742,11 @@ end
     sl = 1.0
     sr = -1.0
 
+    # Test for coverage
+    dummy = 42
+    @test reconstruction_constant(dummy, sl, sr, dummy, dummy, dummy, dummy, dummy) ==
+          (sl, sr)
+
     @test minmod(sl, sr) == 0.0
     @test monotonized_central(sl, sr) == 0.0
     @test superbee(sl, sr) == 0.0