From ed0e8dcdd2aefc2b9715ac936f843fc47ca9009b Mon Sep 17 00:00:00 2001
From: Benedict Geihe <bgeihe@uni-koeln.de>
Date: Mon, 18 Dec 2023 13:37:13 +0100
Subject: [PATCH 01/30] add LMARS flux in 2D

---
 src/equations/compressible_euler_2d.jl | 87 ++++++++++++++++++++++++++
 src/equations/compressible_euler_3d.jl |  5 --
 2 files changed, 87 insertions(+), 5 deletions(-)

diff --git a/src/equations/compressible_euler_2d.jl b/src/equations/compressible_euler_2d.jl
index a992f99eaf4..0e56f934bbf 100644
--- a/src/equations/compressible_euler_2d.jl
+++ b/src/equations/compressible_euler_2d.jl
@@ -809,6 +809,93 @@ end
     return SVector(f1m, f2m, f3m, f4m)
 end
 
+"""
+    FluxLMARS(c)(u_ll, u_rr, orientation_or_normal_direction,
+                 equations::CompressibleEulerEquations2D)
+
+Low Mach number approximate Riemann solver (LMARS) for atmospheric flows using
+an estimate `c` of the speed of sound.
+
+References:
+- Xi Chen et al. (2013)
+  A Control-Volume Model of the Compressible Euler Equations with a Vertical
+  Lagrangian Coordinate
+  [DOI: 10.1175/MWR-D-12-00129.1](https://doi.org/10.1175/mwr-d-12-00129.1)
+"""
+struct FluxLMARS{SpeedOfSound}
+    # Estimate for the speed of sound
+    speed_of_sound::SpeedOfSound
+end
+
+@inline function (flux_lmars::FluxLMARS)(u_ll, u_rr, orientation::Integer,
+                                         equations::CompressibleEulerEquations2D)
+    c = flux_lmars.speed_of_sound
+
+    # Unpack left and right state
+    rho_ll, v1_ll, v2_ll, p_ll = cons2prim(u_ll, equations)
+    rho_rr, v1_rr, v2_rr, p_rr = cons2prim(u_rr, equations)
+
+    if orientation == 1
+        v_ll = v1_ll
+        v_rr = v1_rr
+    else # orientation == 2
+        v_ll = v2_ll
+        v_rr = v2_rr
+    end
+
+    rho = 0.5 * (rho_ll + rho_rr)
+    p = 0.5 * (p_ll + p_rr) - 0.5 * c * rho * (v_rr - v_ll)
+    v = 0.5 * (v_ll + v_rr) - 1 / (2 * c * rho) * (p_rr - p_ll)
+
+    if v >= 0
+        f1, f2, f3, f4 = v * u_ll
+    else
+        f1, f2, f3, f4 = v * u_rr
+    end
+
+    if orientation == 1
+        f2 += p
+    else # orientation == 2
+        f3 += p
+    end
+    f4 += p * v
+
+    return SVector(f1, f2, f3, f4)
+end
+
+@inline function (flux_lmars::FluxLMARS)(u_ll, u_rr, normal_direction::AbstractVector,
+                                         equations::CompressibleEulerEquations2D)
+    c = flux_lmars.speed_of_sound
+
+    # Unpack left and right state
+    rho_ll, v1_ll, v2_ll, p_ll = cons2prim(u_ll, equations)
+    rho_rr, v1_rr, v2_rr, p_rr = cons2prim(u_rr, equations)
+
+    v_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2]
+    v_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2]
+
+    # Note that this is the same as computing v_ll and v_rr with a normalized normal vector
+    # and then multiplying v by `norm_` again, but this version is slightly faster.
+    norm_ = norm(normal_direction)
+
+    rho = 0.5 * (rho_ll + rho_rr)
+    p = 0.5 * (p_ll + p_rr) - 0.5 * c * rho * (v_rr - v_ll) / norm_
+    v = 0.5 * (v_ll + v_rr) - 1 / (2 * c * rho) * (p_rr - p_ll) * norm_
+
+    if v >= 0
+        f1, f2, f3, f4 = u_ll * v
+        # f4 += p_ll * v ???
+    else
+        f1, f2, f3, f4 = u_rr * v
+        # f4 += p_rr * v ???
+    end
+
+    return SVector(f1,
+                   f2 + p * normal_direction[1],
+                   f3 + p * normal_direction[2],
+                   f4 + p * v)
+end
+
 """
     splitting_vanleer_haenel(u, orientation::Integer,
                              equations::CompressibleEulerEquations2D)
diff --git a/src/equations/compressible_euler_3d.jl b/src/equations/compressible_euler_3d.jl
index fc56f58025b..66d35edef43 100644
--- a/src/equations/compressible_euler_3d.jl
+++ b/src/equations/compressible_euler_3d.jl
@@ -944,11 +944,6 @@ References:
   Lagrangian Coordinate
   [DOI: 10.1175/MWR-D-12-00129.1](https://doi.org/10.1175/mwr-d-12-00129.1)
 """
-struct FluxLMARS{SpeedOfSound}
-    # Estimate for the speed of sound
-    speed_of_sound::SpeedOfSound
-end
-
 @inline function (flux_lmars::FluxLMARS)(u_ll, u_rr, orientation::Integer,
                                          equations::CompressibleEulerEquations3D)
     c = flux_lmars.speed_of_sound

From 77f0d461bbac221b864f0bffd913ccce37bdcfb2 Mon Sep 17 00:00:00 2001
From: Benedict Geihe <bgeihe@uni-koeln.de>
Date: Mon, 18 Dec 2023 14:00:47 +0100
Subject: [PATCH 02/30] add dry air warm bubble test case

---
 .../elixir_euler_warm_bubble.jl               | 124 ++++++++++++++++++
 1 file changed, 124 insertions(+)
 create mode 100644 examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl

diff --git a/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl b/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl
new file mode 100644
index 00000000000..abe96ecd70a
--- /dev/null
+++ b/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl
@@ -0,0 +1,124 @@
+using OrdinaryDiffEq
+using Trixi
+
+# Physical constants
+g::Float64   = 9.81           # gravity of earth
+p_0::Float64 = 100_000.0      # reference pressure
+c_p::Float64 = 1004.0         # heat capacity for constant pressure (dry air)
+c_v::Float64 = 717.0          # heat capacity for constant volume (dry air)
+R            = c_p - c_v      # gas constant (dry air)
+gamma        = c_p / c_v      # heat capacity ratio (dry air)
+
+# Warm bubble test from
+# Wicker, L. J., and Skamarock, W. C., 1998: A time-splitting scheme
+# for the elastic equations incorporating second-order Runge–Kutta
+# time differencing. Mon. Wea. Rev., 126, 1992–1999.
+function initial_condition_warm_bubble(x, t, equations::CompressibleEulerEquations2D)
+    # center of perturbation
+    xc = 10000.0
+    zc = 2000.0
+    # radius of perturbation
+    rc = 2000.0
+    # distance of current x to center of perturbation
+    r = sqrt((x[1] - xc)^2 + (x[2] - zc)^2)
+
+    # perturbation in potential temperature
+    θ_ref = 300.0
+    Δθ = 0.0
+    if r <= rc
+        Δθ = 2 * cospi(0.5*r/rc)^2
+    end
+    θ = θ_ref + Δθ # potential temperature
+
+    # Exner pressure, solves hydrostatic equation for x[2]
+    exner = 1 - g / (c_p * θ) * x[2]
+
+    # pressure
+    p = p_0 * exner^(c_p/R)
+
+    # temperature
+    T = θ * exner
+    
+    # density
+    rho = p / (R*T)
+    
+    v1 = 20.0
+    v2 = 0.0
+    E = c_v * T + 0.5 * (v1^2 + v2^2)  
+    return SVector(rho, rho * v1, rho * v2, rho * E)
+end
+
+@inline function source_terms_gravity(u, x, t,
+                                      equations::CompressibleEulerEquations2D)
+    rho, _, rho_v2, _ = u
+    return SVector(zero(eltype(u)), zero(eltype(u)), -g * rho, -g * rho_v2)
+end
+
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+
+equations = CompressibleEulerEquations2D(gamma)
+
+boundary_conditions = (x_neg=boundary_condition_periodic,
+                       x_pos=boundary_condition_periodic,
+                       y_neg=boundary_condition_slip_wall,
+                       y_pos=boundary_condition_slip_wall)
+
+polydeg = 4
+basis = LobattoLegendreBasis(polydeg)
+
+surface_flux = FluxLMARS(340.0)
+volume_flux = flux_ranocha
+volume_integral = VolumeIntegralFluxDifferencing(volume_flux)
+
+solver = DGSEM(basis, surface_flux, volume_integral)
+
+coordinates_min = (     0.0,      0.0)
+coordinates_max = (20_000.0, 10_000.0)
+
+cells_per_dimension = (64, 32)
+mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_warm_bubble, solver,
+                                    source_terms = source_terms_gravity,
+                                    boundary_conditions = boundary_conditions)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 1000.0)
+
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 1000
+
+analysis_callback = AnalysisCallback(semi, interval=analysis_interval,
+                                     extra_analysis_errors=(:entropy_conservation_error,))
+
+alive_callback = AliveCallback(analysis_interval=analysis_interval)
+
+save_solution = SaveSolutionCallback(interval=analysis_interval,
+                                     save_initial_solution=true,
+                                     save_final_solution=true,
+                                     output_directory="out.struct_lmars_ra",
+                                     solution_variables=cons2prim)
+
+stepsize_callback = StepsizeCallback(cfl=0.2)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback,
+                        alive_callback,
+                        save_solution,
+                        stepsize_callback)
+
+###############################################################################
+# run the simulation
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false), 
+            maxiters=1.0e7,
+            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep=false, callback=callbacks);
+
+summary_callback()

From c8b53606e7c94acb420eba2bbedd697545033ae8 Mon Sep 17 00:00:00 2001
From: Benedict Geihe <bgeihe@uni-koeln.de>
Date: Mon, 18 Dec 2023 15:43:24 +0100
Subject: [PATCH 03/30] get formatting right

---
 .../elixir_euler_warm_bubble.jl               | 61 +++++++++----------
 1 file changed, 30 insertions(+), 31 deletions(-)

diff --git a/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl b/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl
index abe96ecd70a..63991d45b88 100644
--- a/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl
+++ b/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl
@@ -2,12 +2,12 @@ using OrdinaryDiffEq
 using Trixi
 
 # Physical constants
-g::Float64   = 9.81           # gravity of earth
-p_0::Float64 = 100_000.0      # reference pressure
-c_p::Float64 = 1004.0         # heat capacity for constant pressure (dry air)
-c_v::Float64 = 717.0          # heat capacity for constant volume (dry air)
-R            = c_p - c_v      # gas constant (dry air)
-gamma        = c_p / c_v      # heat capacity ratio (dry air)
+g::Float64 = 9.81          # gravity of earth
+p_0::Float64 = 100_000.0   # reference pressure
+c_p::Float64 = 1004.0      # heat capacity for constant pressure (dry air)
+c_v::Float64 = 717.0       # heat capacity for constant volume (dry air)
+R = c_p - c_v              # gas constant (dry air)
+gamma = c_p / c_v          # heat capacity ratio (dry air)
 
 # Warm bubble test from
 # Wicker, L. J., and Skamarock, W. C., 1998: A time-splitting scheme
@@ -26,7 +26,7 @@ function initial_condition_warm_bubble(x, t, equations::CompressibleEulerEquatio
     θ_ref = 300.0
     Δθ = 0.0
     if r <= rc
-        Δθ = 2 * cospi(0.5*r/rc)^2
+        Δθ = 2 * cospi(0.5 * r / rc)^2
     end
     θ = θ_ref + Δθ # potential temperature
 
@@ -34,17 +34,17 @@ function initial_condition_warm_bubble(x, t, equations::CompressibleEulerEquatio
     exner = 1 - g / (c_p * θ) * x[2]
 
     # pressure
-    p = p_0 * exner^(c_p/R)
+    p = p_0 * exner^(c_p / R)
 
     # temperature
     T = θ * exner
-    
+
     # density
-    rho = p / (R*T)
-    
+    rho = p / (R * T)
+
     v1 = 20.0
     v2 = 0.0
-    E = c_v * T + 0.5 * (v1^2 + v2^2)  
+    E = c_v * T + 0.5 * (v1^2 + v2^2)
     return SVector(rho, rho * v1, rho * v2, rho * E)
 end
 
@@ -54,16 +54,15 @@ end
     return SVector(zero(eltype(u)), zero(eltype(u)), -g * rho, -g * rho_v2)
 end
 
-
 ###############################################################################
 # semidiscretization of the compressible Euler equations
 
 equations = CompressibleEulerEquations2D(gamma)
 
-boundary_conditions = (x_neg=boundary_condition_periodic,
-                       x_pos=boundary_condition_periodic,
-                       y_neg=boundary_condition_slip_wall,
-                       y_pos=boundary_condition_slip_wall)
+boundary_conditions = (x_neg = boundary_condition_periodic,
+                       x_pos = boundary_condition_periodic,
+                       y_neg = boundary_condition_slip_wall,
+                       y_pos = boundary_condition_slip_wall)
 
 polydeg = 4
 basis = LobattoLegendreBasis(polydeg)
@@ -74,7 +73,7 @@ volume_integral = VolumeIntegralFluxDifferencing(volume_flux)
 
 solver = DGSEM(basis, surface_flux, volume_integral)
 
-coordinates_min = (     0.0,      0.0)
+coordinates_min = (0.0, 0.0)
 coordinates_max = (20_000.0, 10_000.0)
 
 cells_per_dimension = (64, 32)
@@ -95,18 +94,18 @@ summary_callback = SummaryCallback()
 
 analysis_interval = 1000
 
-analysis_callback = AnalysisCallback(semi, interval=analysis_interval,
-                                     extra_analysis_errors=(:entropy_conservation_error,))
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
+                                     extra_analysis_errors = (:entropy_conservation_error,))
 
-alive_callback = AliveCallback(analysis_interval=analysis_interval)
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
 
-save_solution = SaveSolutionCallback(interval=analysis_interval,
-                                     save_initial_solution=true,
-                                     save_final_solution=true,
-                                     output_directory="out.struct_lmars_ra",
-                                     solution_variables=cons2prim)
+save_solution = SaveSolutionCallback(interval = analysis_interval,
+                                     save_initial_solution = true,
+                                     save_final_solution = true,
+                                     output_directory = "out.struct_lmars_ra",
+                                     solution_variables = cons2prim)
 
-stepsize_callback = StepsizeCallback(cfl=0.2)
+stepsize_callback = StepsizeCallback(cfl = 0.2)
 
 callbacks = CallbackSet(summary_callback,
                         analysis_callback,
@@ -116,9 +115,9 @@ callbacks = CallbackSet(summary_callback,
 
 ###############################################################################
 # run the simulation
-sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false), 
-            maxiters=1.0e7,
-            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
-            save_everystep=false, callback=callbacks);
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
+            maxiters = 1.0e7,
+            dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep = false, callback = callbacks);
 
 summary_callback()

From 786ccd990780d28a7171775a89864029b0f2410b Mon Sep 17 00:00:00 2001
From: Benedict Geihe <bgeihe@uni-koeln.de>
Date: Tue, 19 Dec 2023 13:02:11 +0100
Subject: [PATCH 04/30] remove +=

---
 src/equations/compressible_euler_2d.jl | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/equations/compressible_euler_2d.jl b/src/equations/compressible_euler_2d.jl
index 0e56f934bbf..8abf365840b 100644
--- a/src/equations/compressible_euler_2d.jl
+++ b/src/equations/compressible_euler_2d.jl
@@ -854,11 +854,11 @@ end
     end
 
     if orientation == 1
-        f2 += p
+        f2 = f2 + p
     else # orientation == 2
-        f3 += p
+        f3 = f3 + p
     end
-    f4 += p * v
+    f4 = f4 + p * v
 
     return SVector(f1, f2, f3, f4)
 end

From 62512dc395171ceca6cce9dbdc3db3b99a4ba141 Mon Sep 17 00:00:00 2001
From: Benedict Geihe <bgeihe@uni-koeln.de>
Date: Tue, 19 Dec 2023 13:07:08 +0100
Subject: [PATCH 05/30] add DOI

---
 examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl | 9 +++++----
 1 file changed, 5 insertions(+), 4 deletions(-)

diff --git a/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl b/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl
index 63991d45b88..5e3bba2e77b 100644
--- a/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl
+++ b/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl
@@ -9,10 +9,11 @@ c_v::Float64 = 717.0       # heat capacity for constant volume (dry air)
 R = c_p - c_v              # gas constant (dry air)
 gamma = c_p / c_v          # heat capacity ratio (dry air)
 
-# Warm bubble test from
-# Wicker, L. J., and Skamarock, W. C., 1998: A time-splitting scheme
-# for the elastic equations incorporating second-order Runge–Kutta
-# time differencing. Mon. Wea. Rev., 126, 1992–1999.
+# Warm bubble test case from
+# Wicker, L. J., and Skamarock, W. C.
+# A time-splitting scheme for the elastic equations incorporating second-order Runge–Kutta
+# time differencing
+# [DOI: 10.1175/1520-0493(1998)126%3C1992:ATSSFT%3E2.0.CO;2](https://doi.org/10.1175/1520-0493(1998)126%3C1992:ATSSFT%3E2.0.CO;2)
 function initial_condition_warm_bubble(x, t, equations::CompressibleEulerEquations2D)
     # center of perturbation
     xc = 10000.0

From 8043aafcae960eed263ebfb95c21fa80a52d58fa Mon Sep 17 00:00:00 2001
From: Benedict Geihe <bgeihe@uni-koeln.de>
Date: Tue, 19 Dec 2023 13:18:07 +0100
Subject: [PATCH 06/30] cleaned up variables (naming, scope)

---
 .../elixir_euler_warm_bubble.jl               | 34 +++++++++----------
 1 file changed, 17 insertions(+), 17 deletions(-)

diff --git a/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl b/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl
index 5e3bba2e77b..a962aade6e1 100644
--- a/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl
+++ b/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl
@@ -2,12 +2,10 @@ using OrdinaryDiffEq
 using Trixi
 
 # Physical constants
-g::Float64 = 9.81          # gravity of earth
-p_0::Float64 = 100_000.0   # reference pressure
-c_p::Float64 = 1004.0      # heat capacity for constant pressure (dry air)
-c_v::Float64 = 717.0       # heat capacity for constant volume (dry air)
-R = c_p - c_v              # gas constant (dry air)
-gamma = c_p / c_v          # heat capacity ratio (dry air)
+g::Float64 = 9.81            # gravity of earth
+c_p::Float64 = 1004.0        # heat capacity for constant pressure (dry air)
+c_v::Float64 = 717.0         # heat capacity for constant volume (dry air)
+gamma::Float64 = c_p / c_v   # heat capacity ratio (dry air)
 
 # Warm bubble test case from
 # Wicker, L. J., and Skamarock, W. C.
@@ -16,29 +14,31 @@ gamma = c_p / c_v          # heat capacity ratio (dry air)
 # [DOI: 10.1175/1520-0493(1998)126%3C1992:ATSSFT%3E2.0.CO;2](https://doi.org/10.1175/1520-0493(1998)126%3C1992:ATSSFT%3E2.0.CO;2)
 function initial_condition_warm_bubble(x, t, equations::CompressibleEulerEquations2D)
     # center of perturbation
-    xc = 10000.0
-    zc = 2000.0
+    center_x = 10000.0
+    center_z = 2000.0
     # radius of perturbation
-    rc = 2000.0
+    radius = 2000.0
     # distance of current x to center of perturbation
-    r = sqrt((x[1] - xc)^2 + (x[2] - zc)^2)
+    r = sqrt((x[1] - center_x)^2 + (x[2] - center_z)^2)
 
     # perturbation in potential temperature
-    θ_ref = 300.0
-    Δθ = 0.0
-    if r <= rc
-        Δθ = 2 * cospi(0.5 * r / rc)^2
+    potential_temperature_ref = 300.0
+    potential_temperature_perturbation = 0.0
+    if r <= radius
+        potential_temperature_perturbation = 2 * cospi(0.5 * r / radius)^2
     end
-    θ = θ_ref + Δθ # potential temperature
+    potential_temperature = potential_temperature_ref + potential_temperature_perturbation
 
     # Exner pressure, solves hydrostatic equation for x[2]
-    exner = 1 - g / (c_p * θ) * x[2]
+    exner = 1 - g / (c_p * potential_temperature) * x[2]
 
     # pressure
+    p_0 = 100_000.0  # reference pressure
+    R = c_p - c_v    # gas constant (dry air)
     p = p_0 * exner^(c_p / R)
 
     # temperature
-    T = θ * exner
+    T = potential_temperature * exner
 
     # density
     rho = p / (R * T)

From b6e527bd63767e0a124a274d57e986df9a88b766 Mon Sep 17 00:00:00 2001
From: Benedict Geihe <bgeihe@uni-koeln.de>
Date: Tue, 19 Dec 2023 13:18:34 +0100
Subject: [PATCH 07/30] reduce run time of test case

---
 examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl | 8 +++++---
 1 file changed, 5 insertions(+), 3 deletions(-)

diff --git a/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl b/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl
index a962aade6e1..d250136d5d3 100644
--- a/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl
+++ b/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl
@@ -77,7 +77,8 @@ solver = DGSEM(basis, surface_flux, volume_integral)
 coordinates_min = (0.0, 0.0)
 coordinates_max = (20_000.0, 10_000.0)
 
-cells_per_dimension = (64, 32)
+# increase resolution to (64, 32)
+cells_per_dimension = (32, 16)
 mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max)
 
 semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_warm_bubble, solver,
@@ -87,7 +88,8 @@ semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_warm_bubb
 ###############################################################################
 # ODE solvers, callbacks etc.
 
-tspan = (0.0, 1000.0)
+# change final time to 1000.0
+tspan = (0.0, 10.0)
 
 ode = semidiscretize(semi, tspan)
 
@@ -103,7 +105,7 @@ alive_callback = AliveCallback(analysis_interval = analysis_interval)
 save_solution = SaveSolutionCallback(interval = analysis_interval,
                                      save_initial_solution = true,
                                      save_final_solution = true,
-                                     output_directory = "out.struct_lmars_ra",
+                                     output_directory = "out",
                                      solution_variables = cons2prim)
 
 stepsize_callback = StepsizeCallback(cfl = 0.2)

From 7cfe952baa099842d8a359c8ff797c1bd06558a0 Mon Sep 17 00:00:00 2001
From: Benedict Geihe <bgeihe@uni-koeln.de>
Date: Thu, 21 Dec 2023 09:14:33 +0100
Subject: [PATCH 08/30] Revert "reduce run time of test case"

This reverts commit b6e527bd63767e0a124a274d57e986df9a88b766.
---
 examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl | 8 +++-----
 1 file changed, 3 insertions(+), 5 deletions(-)

diff --git a/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl b/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl
index d250136d5d3..a962aade6e1 100644
--- a/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl
+++ b/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl
@@ -77,8 +77,7 @@ solver = DGSEM(basis, surface_flux, volume_integral)
 coordinates_min = (0.0, 0.0)
 coordinates_max = (20_000.0, 10_000.0)
 
-# increase resolution to (64, 32)
-cells_per_dimension = (32, 16)
+cells_per_dimension = (64, 32)
 mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max)
 
 semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_warm_bubble, solver,
@@ -88,8 +87,7 @@ semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_warm_bubb
 ###############################################################################
 # ODE solvers, callbacks etc.
 
-# change final time to 1000.0
-tspan = (0.0, 10.0)
+tspan = (0.0, 1000.0)
 
 ode = semidiscretize(semi, tspan)
 
@@ -105,7 +103,7 @@ alive_callback = AliveCallback(analysis_interval = analysis_interval)
 save_solution = SaveSolutionCallback(interval = analysis_interval,
                                      save_initial_solution = true,
                                      save_final_solution = true,
-                                     output_directory = "out",
+                                     output_directory = "out.struct_lmars_ra",
                                      solution_variables = cons2prim)
 
 stepsize_callback = StepsizeCallback(cfl = 0.2)

From 03db8f134a7d6beb7668b3ef87840592e4424eec Mon Sep 17 00:00:00 2001
From: Benedict Geihe <bgeihe@uni-koeln.de>
Date: Thu, 21 Dec 2023 09:15:24 +0100
Subject: [PATCH 09/30] change output folder

---
 examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl b/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl
index a962aade6e1..24ba7d1aad0 100644
--- a/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl
+++ b/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl
@@ -103,7 +103,7 @@ alive_callback = AliveCallback(analysis_interval = analysis_interval)
 save_solution = SaveSolutionCallback(interval = analysis_interval,
                                      save_initial_solution = true,
                                      save_final_solution = true,
-                                     output_directory = "out.struct_lmars_ra",
+                                     output_directory = "out",
                                      solution_variables = cons2prim)
 
 stepsize_callback = StepsizeCallback(cfl = 0.2)

From 78f22c31874b3037462d20d8f7896cfaf44e8162 Mon Sep 17 00:00:00 2001
From: Benedict Geihe <bgeihe@uni-koeln.de>
Date: Thu, 21 Dec 2023 09:16:00 +0100
Subject: [PATCH 10/30] change energy term in LMARS solver to use p_l/r

---
 src/equations/compressible_euler_2d.jl | 13 +++++++++----
 src/equations/compressible_euler_3d.jl | 20 +++++++++++++-------
 2 files changed, 22 insertions(+), 11 deletions(-)

diff --git a/src/equations/compressible_euler_2d.jl b/src/equations/compressible_euler_2d.jl
index 8abf365840b..f7c3b5c7fa2 100644
--- a/src/equations/compressible_euler_2d.jl
+++ b/src/equations/compressible_euler_2d.jl
@@ -847,10 +847,14 @@ end
     p = 0.5 * (p_ll + p_rr) - 0.5 * c * rho * (v_rr - v_ll)
     v = 0.5 * (v_ll + v_rr) - 1 / (2 * c * rho) * (p_rr - p_ll)
 
+    # We treat the energy term in analogy to the potential temperature term in the paper by
+    # Chen et al., i.e. we use p_ll and p_rr, and not p
     if v >= 0
         f1, f2, f3, f4 = v * u_ll
+        f4 = f4 + p_ll * v
     else
         f1, f2, f3, f4 = v * u_rr
+        f4 = f4 + p_rr * v
     end
 
     if orientation == 1
@@ -858,7 +862,6 @@ end
     else # orientation == 2
         f3 = f3 + p
     end
-    f4 = f4 + p * v
 
     return SVector(f1, f2, f3, f4)
 end
@@ -882,18 +885,20 @@ end
     p = 0.5 * (p_ll + p_rr) - 0.5 * c * rho * (v_rr - v_ll) / norm_
     v = 0.5 * (v_ll + v_rr) - 1 / (2 * c * rho) * (p_rr - p_ll) * norm_
 
+    # We treat the energy term in analogy to the potential temperature term in the paper by
+    # Chen et al., i.e. we use p_ll and p_rr, and not p
     if v >= 0
         f1, f2, f3, f4 = u_ll * v
-        # f4 += p_ll * v ???
+        f4 = f4 + p_ll * v
     else
         f1, f2, f3, f4 = u_rr * v
-        # f4 += p_rr * v ???
+        f4 = f4 + p_rr * v
     end
 
     return SVector(f1,
                    f2 + p * normal_direction[1],
                    f3 + p * normal_direction[2],
-                   f4 + p * v)
+                   f4)
 end
 
 """
diff --git a/src/equations/compressible_euler_3d.jl b/src/equations/compressible_euler_3d.jl
index 66d35edef43..2a2b0119d03 100644
--- a/src/equations/compressible_euler_3d.jl
+++ b/src/equations/compressible_euler_3d.jl
@@ -967,10 +967,14 @@ References:
     p = 0.5 * (p_ll + p_rr) - 0.5 * c * rho * (v_rr - v_ll)
     v = 0.5 * (v_ll + v_rr) - 1 / (2 * c * rho) * (p_rr - p_ll)
 
+    # We treat the energy term in analogy to the potential temperature term in the paper by
+    # Chen et al., i.e. we use p_ll and p_rr, and not p
     if v >= 0
         f1, f2, f3, f4, f5 = v * u_ll
+        f5 = f5 + p_ll * v
     else
         f1, f2, f3, f4, f5 = v * u_rr
+        f5 = f5 + p_rr * v
     end
 
     if orientation == 1
@@ -980,7 +984,6 @@ References:
     else # orientation == 3
         f4 += p
     end
-    f5 += p * v
 
     return SVector(f1, f2, f3, f4, f5)
 end
@@ -1006,18 +1009,21 @@ end
     p = 0.5 * (p_ll + p_rr) - 0.5 * c * rho * (v_rr - v_ll) / norm_
     v = 0.5 * (v_ll + v_rr) - 1 / (2 * c * rho) * (p_rr - p_ll) * norm_
 
+    # We treat the energy term in analogy to the potential temperature term in the paper by
+    # Chen et al., i.e. we use p_ll and p_rr, and not p
     if v >= 0
         f1, f2, f3, f4, f5 = v * u_ll
+        f5 = f5 + p_ll * v
     else
         f1, f2, f3, f4, f5 = v * u_rr
+        f5 = f5 + p_rr * v
     end
 
-    f2 += p * normal_direction[1]
-    f3 += p * normal_direction[2]
-    f4 += p * normal_direction[3]
-    f5 += p * v
-
-    return SVector(f1, f2, f3, f4, f5)
+    return SVector(f1,
+                   f2 + p * normal_direction[1],
+                   f3 + p * normal_direction[2],
+                   f4 + p * normal_direction[3],
+                   f5)
 end
 
 # Calculate maximum wave speed for local Lax-Friedrichs-type dissipation as the

From 700dd30b9236b2b900fee4b83fd16f617f480730 Mon Sep 17 00:00:00 2001
From: Benedict Geihe <bgeihe@uni-koeln.de>
Date: Thu, 21 Dec 2023 09:29:31 +0100
Subject: [PATCH 11/30] add lmars consistency checks

---
 test/test_unit.jl | 48 ++++++++++++++++++++++++++++++++++++++++++++---
 1 file changed, 45 insertions(+), 3 deletions(-)

diff --git a/test/test_unit.jl b/test/test_unit.jl
index b3ed29d38e3..abb97c0b9ac 100644
--- a/test/test_unit.jl
+++ b/test/test_unit.jl
@@ -1247,6 +1247,49 @@ end
     end
 end
 
+@timed_testset "Consistency check for LMARS flux" begin
+    equations = CompressibleEulerEquations2D(1.4)
+    flux_lmars = FluxLMARS(340)
+
+    normal_directions = [SVector(1.0, 0.0),
+        SVector(0.0, 1.0),
+        SVector(0.5, -0.5),
+        SVector(-1.2, 0.3)]
+    orientations = [1, 2]
+    u_values = [SVector(1.0, 0.5, -0.7, 1.0),
+        SVector(1.5, -0.2, 0.1, 5.0)]
+
+    for u in u_values, orientation in orientations
+        @test flux_lmars(u, u, orientation, equations) ≈
+              flux_lmars(u, orientation, equations)
+    end
+
+    for u in u_values, normal_direction in normal_directions
+        @test flux_lmars(u, u, normal_direction, equations) ≈
+              flux_lmars(u, normal_direction, equations)
+    end
+
+    equations = CompressibleEulerEquations3D(1.4)
+    normal_directions = [SVector(1.0, 0.0, 0.0),
+        SVector(0.0, 1.0, 0.0),
+        SVector(0.0, 0.0, 1.0),
+        SVector(0.5, -0.5, 0.2),
+        SVector(-1.2, 0.3, 1.4)]
+    orientations = [1, 2, 3]
+    u_values = [SVector(1.0, 0.5, -0.7, 0.1, 1.0),
+        SVector(1.5, -0.2, 0.1, 0.2, 5.0)]
+
+    for u in u_values, orientation in orientations
+        @test flux_lmars(u, u, orientation, equations) ≈
+              flux_lmars(u, orientation, equations)
+    end
+
+    for u in u_values, normal_direction in normal_directions
+        @test flux_lmars(u, u, normal_direction, equations) ≈
+              flux_lmars(u, normal_direction, equations)
+    end
+end
+
 @testset "FluxRotated vs. direct implementation" begin
     @timed_testset "CompressibleEulerMulticomponentEquations2D" begin
         equations = CompressibleEulerMulticomponentEquations2D(gammas = (1.4, 1.4),
@@ -1280,7 +1323,7 @@ end
         u_values = [SVector(1.0, 0.5, -0.7, 1.0),
             SVector(1.5, -0.2, 0.1, 5.0)]
         fluxes = [flux_central, flux_ranocha, flux_shima_etal, flux_kennedy_gruber,
-            flux_hll, FluxHLL(min_max_speed_davis), flux_hlle]
+            FluxLMARS(340), flux_hll, FluxHLL(min_max_speed_davis), flux_hlle]
 
         for f_std in fluxes
             f_rot = FluxRotated(f_std)
@@ -1303,8 +1346,7 @@ end
         u_values = [SVector(1.0, 0.5, -0.7, 0.1, 1.0),
             SVector(1.5, -0.2, 0.1, 0.2, 5.0)]
         fluxes = [flux_central, flux_ranocha, flux_shima_etal, flux_kennedy_gruber,
-            FluxLMARS(340),
-            flux_hll, FluxHLL(min_max_speed_davis), flux_hlle]
+            FluxLMARS(340), flux_hll, FluxHLL(min_max_speed_davis), flux_hlle]
 
         for f_std in fluxes
             f_rot = FluxRotated(f_std)

From e25f2dc9a67823a84d160fa4f15b2e8b982b5bf8 Mon Sep 17 00:00:00 2001
From: Benedict Geihe <bgeihe@uni-koeln.de>
Date: Thu, 21 Dec 2023 12:25:37 +0100
Subject: [PATCH 12/30] switched to kennedy gruber flux

---
 examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl b/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl
index 24ba7d1aad0..f73e5681392 100644
--- a/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl
+++ b/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl
@@ -69,7 +69,7 @@ polydeg = 4
 basis = LobattoLegendreBasis(polydeg)
 
 surface_flux = FluxLMARS(340.0)
-volume_flux = flux_ranocha
+volume_flux = flux_kennedy_gruber
 volume_integral = VolumeIntegralFluxDifferencing(volume_flux)
 
 solver = DGSEM(basis, surface_flux, volume_integral)

From 115178960c0d78e645c151611cc84bed915e08cf Mon Sep 17 00:00:00 2001
From: Benedict Geihe <bgeihe@uni-koeln.de>
Date: Thu, 21 Dec 2023 12:26:04 +0100
Subject: [PATCH 13/30] add euler warm bubble elixir to tests

---
 test/test_structured_2d.jl | 27 +++++++++++++++++++++++++++
 1 file changed, 27 insertions(+)

diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl
index 1addc29e3e6..3a5b1c39f47 100644
--- a/test/test_structured_2d.jl
+++ b/test/test_structured_2d.jl
@@ -611,6 +611,33 @@ end
     end
 end
 
+@trixi_testset "elixir_euler_warm_bubble.jl" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR,
+                                 "elixir_euler_warm_bubble.jl"),
+                        l2=[
+                            0.00019389374978502348,
+                            0.030867307039472824,
+                            0.04541800001263111,
+                            43.89530131084914,
+                        ],
+                        linf=[
+                            0.0015906758532284737,
+                            0.17545022864421256,
+                            0.3729842228441566,
+                            307.76174927831744,
+                        ],
+                        cells_per_dimension=(32, 16),
+                        tspan=(0.0, 10.))
+    # 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)) < 100
+    end
+end
+
 @trixi_testset "elixir_eulerpolytropic_convergence.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulerpolytropic_convergence.jl"),
                         l2=[

From d91e0f61a28ea7f54d4bfe1132f4c84ac96f6185 Mon Sep 17 00:00:00 2001
From: Benedict Geihe <bgeihe@uni-koeln.de>
Date: Thu, 21 Dec 2023 13:16:58 +0100
Subject: [PATCH 14/30] adapt errors due to change flux

---
 test/test_p4est_3d.jl | 40 ++++++++++++++++++++--------------------
 1 file changed, 20 insertions(+), 20 deletions(-)

diff --git a/test/test_p4est_3d.jl b/test/test_p4est_3d.jl
index f2467f30204..2ec8defbcdc 100644
--- a/test/test_p4est_3d.jl
+++ b/test/test_p4est_3d.jl
@@ -352,18 +352,18 @@ end
     @test_trixi_include(joinpath(EXAMPLES_DIR,
                                  "elixir_euler_circular_wind_nonconforming.jl"),
                         l2=[
-                            1.573832094977477e-7,
-                            3.863090659429634e-5,
-                            3.867293305754584e-5,
-                            3.686550296950078e-5,
-                            0.05508968493733932,
+                            1.5737711609657832e-7,
+                            3.8630261900166194e-5,
+                            3.8672287531936816e-5,
+                            3.6865116098660796e-5,
+                            0.05508620970403884,
                         ],
                         linf=[
-                            2.2695202613887133e-6,
-                            0.0005314968179916946,
-                            0.0005314969614147458,
-                            0.0005130280733059617,
-                            0.7944959432352334,
+                            2.268845333053271e-6,
+                            0.000531462302113539,
+                            0.0005314624461298934,
+                            0.0005129931254772464,
+                            0.7942778058932163,
                         ],
                         tspan=(0.0, 2e2),
                         coverage_override=(trees_per_cube_face = (1, 1), polydeg = 3)) # Prevent long compile time in CI
@@ -381,18 +381,18 @@ end
     @test_trixi_include(joinpath(EXAMPLES_DIR,
                                  "elixir_euler_baroclinic_instability.jl"),
                         l2=[
-                            6.725065410642336e-7,
-                            0.00021710117340245454,
-                            0.000438679759422352,
-                            0.00020836356588024185,
-                            0.07602006689579247,
+                            6.725093801700048e-7,
+                            0.00021710076010951073,
+                            0.0004386796338203878,
+                            0.00020836270267103122,
+                            0.07601887903440395,
                         ],
                         linf=[
-                            1.9101671995258585e-5,
-                            0.029803626911022396,
-                            0.04847630924006063,
-                            0.022001371349740104,
-                            4.847761006938526,
+                            1.9107530539574924e-5,
+                            0.02980358831035801,
+                            0.048476331898047564,
+                            0.02200137344113612,
+                            4.848310144356219,
                         ],
                         tspan=(0.0, 1e2),
                         # Decrease tolerance of adaptive time stepping to get similar results across different systems

From b38d57e74518cd5e2e80b0954658ded020ce286d Mon Sep 17 00:00:00 2001
From: Benedict Geihe <bgeihe@uni-koeln.de>
Date: Thu, 21 Dec 2023 13:17:40 +0100
Subject: [PATCH 15/30] add warm bubble test with TreeMesh

---
 .../tree_2d_dgsem/elixir_euler_warm_bubble.jl | 126 ++++++++++++++++++
 test/test_tree_2d_euler.jl                    |  26 ++++
 2 files changed, 152 insertions(+)
 create mode 100644 examples/tree_2d_dgsem/elixir_euler_warm_bubble.jl

diff --git a/examples/tree_2d_dgsem/elixir_euler_warm_bubble.jl b/examples/tree_2d_dgsem/elixir_euler_warm_bubble.jl
new file mode 100644
index 00000000000..8395656872a
--- /dev/null
+++ b/examples/tree_2d_dgsem/elixir_euler_warm_bubble.jl
@@ -0,0 +1,126 @@
+using OrdinaryDiffEq
+using Trixi
+
+# Physical constants
+g::Float64 = 9.81            # gravity of earth
+c_p::Float64 = 1004.0        # heat capacity for constant pressure (dry air)
+c_v::Float64 = 717.0         # heat capacity for constant volume (dry air)
+gamma::Float64 = c_p / c_v   # heat capacity ratio (dry air)
+
+# Warm bubble test case from
+# Wicker, L. J., and Skamarock, W. C.
+# A time-splitting scheme for the elastic equations incorporating second-order Runge–Kutta
+# time differencing
+# [DOI: 10.1175/1520-0493(1998)126%3C1992:ATSSFT%3E2.0.CO;2](https://doi.org/10.1175/1520-0493(1998)126%3C1992:ATSSFT%3E2.0.CO;2)
+function initial_condition_warm_bubble(x, t, equations::CompressibleEulerEquations2D)
+    # center of perturbation
+    center_x = 10000.0
+    center_z = 2000.0
+    # radius of perturbation
+    radius = 2000.0
+    # distance of current x to center of perturbation
+    r = sqrt((x[1] - center_x)^2 + (x[2] - center_z)^2)
+
+    # perturbation in potential temperature
+    potential_temperature_ref = 300.0
+    potential_temperature_perturbation = 0.0
+    if r <= radius
+        potential_temperature_perturbation = 2 * cospi(0.5 * r / radius)^2
+    end
+    potential_temperature = potential_temperature_ref + potential_temperature_perturbation
+
+    # Exner pressure, solves hydrostatic equation for x[2]
+    exner = 1 - g / (c_p * potential_temperature) * x[2]
+
+    # pressure
+    p_0 = 100_000.0  # reference pressure
+    R = c_p - c_v    # gas constant (dry air)
+    p = p_0 * exner^(c_p / R)
+
+    # temperature
+    T = potential_temperature * exner
+
+    # density
+    rho = p / (R * T)
+
+    v1 = 20.0
+    v2 = 0.0
+    E = c_v * T + 0.5 * (v1^2 + v2^2)
+    return SVector(rho, rho * v1, rho * v2, rho * E)
+end
+
+@inline function source_terms_gravity(u, x, t,
+                                      equations::CompressibleEulerEquations2D)
+    rho, _, rho_v2, _ = u
+    return SVector(zero(eltype(u)), zero(eltype(u)), -g * rho, -g * rho_v2)
+end
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+
+equations = CompressibleEulerEquations2D(gamma)
+
+boundary_conditions = (x_neg = boundary_condition_periodic,
+                       x_pos = boundary_condition_periodic,
+                       y_neg = boundary_condition_slip_wall,
+                       y_pos = boundary_condition_slip_wall)
+
+polydeg = 4
+basis = LobattoLegendreBasis(polydeg)
+
+surface_flux = FluxLMARS(340.0)
+volume_flux = flux_kennedy_gruber
+volume_integral = VolumeIntegralFluxDifferencing(volume_flux)
+
+solver = DGSEM(basis, surface_flux, volume_integral)
+
+coordinates_min = (0.0, 0.0)
+coordinates_max = (20_000.0, 10_000.0)
+
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 6,
+                n_cells_max = 10_000,
+                periodicity = (true, false))
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_warm_bubble, solver,
+                                    source_terms = source_terms_gravity,
+                                    boundary_conditions = boundary_conditions)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 1000.0)
+
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 1000
+
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
+                                     extra_analysis_errors = (:entropy_conservation_error,))
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(interval = analysis_interval,
+                                     save_initial_solution = true,
+                                     save_final_solution = true,
+                                     output_directory = "out",
+                                     solution_variables = cons2prim)
+
+stepsize_callback = StepsizeCallback(cfl = 0.2)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback,
+                        alive_callback,
+                        save_solution,
+                        stepsize_callback)
+
+###############################################################################
+# run the simulation
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
+            maxiters = 1.0e7,
+            dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep = false, callback = callbacks);
+
+summary_callback()
diff --git a/test/test_tree_2d_euler.jl b/test/test_tree_2d_euler.jl
index 65899cd5263..39aa82320e1 100644
--- a/test/test_tree_2d_euler.jl
+++ b/test/test_tree_2d_euler.jl
@@ -834,6 +834,32 @@ end
     end
 end
 
+@trixi_testset "elixir_euler_warm_bubble.jl" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_warm_bubble.jl"),
+                        l2=[
+                            0.00013820060399128458,
+                            0.020833682516209616,
+                            0.03327372067821428,
+                            31.39884011181286,
+                        ],
+                        linf=[
+                            0.0016677536924397662,
+                            0.15617627577194781,
+                            0.3368120143100738,
+                            331.40871663423604,
+                        ],
+                        tspan=(0.0, 10.0),
+                        initial_refinement_level=4)
+    # 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)) < 100
+    end
+end
+
 # Coverage test for all initial conditions
 @testset "Compressible Euler: Tests for initial conditions" begin
     @trixi_testset "elixir_euler_vortex.jl one step with initial_condition_constant" begin

From e373c08591902f9e4d0bc60c9bce54c2cba9299c Mon Sep 17 00:00:00 2001
From: Benedict Geihe <bgeihe@uni-koeln.de>
Date: Thu, 21 Dec 2023 13:18:09 +0100
Subject: [PATCH 16/30] fix unit test

---
 test/test_unit.jl | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/test/test_unit.jl b/test/test_unit.jl
index abb97c0b9ac..45e625cb760 100644
--- a/test/test_unit.jl
+++ b/test/test_unit.jl
@@ -1261,12 +1261,12 @@ end
 
     for u in u_values, orientation in orientations
         @test flux_lmars(u, u, orientation, equations) ≈
-              flux_lmars(u, orientation, equations)
+              flux(u, orientation, equations)
     end
 
     for u in u_values, normal_direction in normal_directions
         @test flux_lmars(u, u, normal_direction, equations) ≈
-              flux_lmars(u, normal_direction, equations)
+              flux(u, normal_direction, equations)
     end
 
     equations = CompressibleEulerEquations3D(1.4)
@@ -1281,12 +1281,12 @@ end
 
     for u in u_values, orientation in orientations
         @test flux_lmars(u, u, orientation, equations) ≈
-              flux_lmars(u, orientation, equations)
+              flux(u, orientation, equations)
     end
 
     for u in u_values, normal_direction in normal_directions
         @test flux_lmars(u, u, normal_direction, equations) ≈
-              flux_lmars(u, normal_direction, equations)
+              flux(u, normal_direction, equations)
     end
 end
 

From 285f09b3008f5f212b00ab528cc11e57d4d1efef Mon Sep 17 00:00:00 2001
From: Benedict Geihe <bgeihe@uni-koeln.de>
Date: Thu, 21 Dec 2023 13:19:24 +0100
Subject: [PATCH 17/30] fix format

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

diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl
index 3a5b1c39f47..f6857a3df8a 100644
--- a/test/test_structured_2d.jl
+++ b/test/test_structured_2d.jl
@@ -627,7 +627,7 @@ end
                             307.76174927831744,
                         ],
                         cells_per_dimension=(32, 16),
-                        tspan=(0.0, 10.))
+                        tspan=(0.0, 10.0))
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
     let

From 63d66ce937b78b21b73cc865c7886e2d4f274f74 Mon Sep 17 00:00:00 2001
From: Benedict Geihe <bgeihe@uni-koeln.de>
Date: Tue, 9 Jan 2024 13:41:50 +0100
Subject: [PATCH 18/30] adapt polynomial degree and CFL number

---
 examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl | 4 ++--
 examples/tree_2d_dgsem/elixir_euler_warm_bubble.jl       | 4 ++--
 2 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl b/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl
index f73e5681392..be52f674a4e 100644
--- a/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl
+++ b/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl
@@ -65,7 +65,7 @@ boundary_conditions = (x_neg = boundary_condition_periodic,
                        y_neg = boundary_condition_slip_wall,
                        y_pos = boundary_condition_slip_wall)
 
-polydeg = 4
+polydeg = 3
 basis = LobattoLegendreBasis(polydeg)
 
 surface_flux = FluxLMARS(340.0)
@@ -106,7 +106,7 @@ save_solution = SaveSolutionCallback(interval = analysis_interval,
                                      output_directory = "out",
                                      solution_variables = cons2prim)
 
-stepsize_callback = StepsizeCallback(cfl = 0.2)
+stepsize_callback = StepsizeCallback(cfl = 1.0)
 
 callbacks = CallbackSet(summary_callback,
                         analysis_callback,
diff --git a/examples/tree_2d_dgsem/elixir_euler_warm_bubble.jl b/examples/tree_2d_dgsem/elixir_euler_warm_bubble.jl
index 8395656872a..15c6e6eb52b 100644
--- a/examples/tree_2d_dgsem/elixir_euler_warm_bubble.jl
+++ b/examples/tree_2d_dgsem/elixir_euler_warm_bubble.jl
@@ -65,7 +65,7 @@ boundary_conditions = (x_neg = boundary_condition_periodic,
                        y_neg = boundary_condition_slip_wall,
                        y_pos = boundary_condition_slip_wall)
 
-polydeg = 4
+polydeg = 3
 basis = LobattoLegendreBasis(polydeg)
 
 surface_flux = FluxLMARS(340.0)
@@ -108,7 +108,7 @@ save_solution = SaveSolutionCallback(interval = analysis_interval,
                                      output_directory = "out",
                                      solution_variables = cons2prim)
 
-stepsize_callback = StepsizeCallback(cfl = 0.2)
+stepsize_callback = StepsizeCallback(cfl = 1.0)
 
 callbacks = CallbackSet(summary_callback,
                         analysis_callback,

From 82d1d7617929ae817a873b4a38f33b67f8cdabc4 Mon Sep 17 00:00:00 2001
From: Benedict Geihe <bgeihe@uni-koeln.de>
Date: Tue, 9 Jan 2024 13:46:18 +0100
Subject: [PATCH 19/30] fix format

---
 test/test_unit.jl | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/test/test_unit.jl b/test/test_unit.jl
index fff6d915033..e8a8effbe29 100644
--- a/test/test_unit.jl
+++ b/test/test_unit.jl
@@ -1363,7 +1363,8 @@ end
         u_values = [SVector(1.0, 0.5, -0.7, 1.0),
             SVector(1.5, -0.2, 0.1, 5.0)]
         fluxes = [flux_central, flux_ranocha, flux_shima_etal, flux_kennedy_gruber,
-            FluxLMARS(340), flux_hll, FluxHLL(min_max_speed_davis), flux_hlle, flux_hllc]
+            FluxLMARS(340), flux_hll, FluxHLL(min_max_speed_davis), flux_hlle, flux_hllc,
+        ]
 
         for f_std in fluxes
             f_rot = FluxRotated(f_std)
@@ -1386,7 +1387,8 @@ end
         u_values = [SVector(1.0, 0.5, -0.7, 0.1, 1.0),
             SVector(1.5, -0.2, 0.1, 0.2, 5.0)]
         fluxes = [flux_central, flux_ranocha, flux_shima_etal, flux_kennedy_gruber,
-            FluxLMARS(340), flux_hll, FluxHLL(min_max_speed_davis), flux_hlle, flux_hllc]
+            FluxLMARS(340), flux_hll, FluxHLL(min_max_speed_davis), flux_hlle, flux_hllc,
+        ]
 
         for f_std in fluxes
             f_rot = FluxRotated(f_std)

From 17e11a21be992651731516e21e10c99774934a74 Mon Sep 17 00:00:00 2001
From: Benedict Geihe <bgeihe@uni-koeln.de>
Date: Tue, 9 Jan 2024 14:23:54 +0100
Subject: [PATCH 20/30] adapt tests due to changed parameters

---
 test/test_structured_2d.jl | 16 ++++++++--------
 test/test_tree_2d_euler.jl | 16 ++++++++--------
 2 files changed, 16 insertions(+), 16 deletions(-)

diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl
index f6857a3df8a..e5d45ebcc07 100644
--- a/test/test_structured_2d.jl
+++ b/test/test_structured_2d.jl
@@ -615,16 +615,16 @@ end
     @test_trixi_include(joinpath(EXAMPLES_DIR,
                                  "elixir_euler_warm_bubble.jl"),
                         l2=[
-                            0.00019389374978502348,
-                            0.030867307039472824,
-                            0.04541800001263111,
-                            43.89530131084914,
+                            0.00019387402388722496,
+                            0.03086514388623955,
+                            0.04541427917165,
+                            43.892826583444716,
                         ],
                         linf=[
-                            0.0015906758532284737,
-                            0.17545022864421256,
-                            0.3729842228441566,
-                            307.76174927831744,
+                            0.0015942305974430138,
+                            0.17449778969139373,
+                            0.3729704262394843,
+                            307.6706958565337,
                         ],
                         cells_per_dimension=(32, 16),
                         tspan=(0.0, 10.0))
diff --git a/test/test_tree_2d_euler.jl b/test/test_tree_2d_euler.jl
index 39aa82320e1..3e3f69d36c3 100644
--- a/test/test_tree_2d_euler.jl
+++ b/test/test_tree_2d_euler.jl
@@ -837,16 +837,16 @@ end
 @trixi_testset "elixir_euler_warm_bubble.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_warm_bubble.jl"),
                         l2=[
-                            0.00013820060399128458,
-                            0.020833682516209616,
-                            0.03327372067821428,
-                            31.39884011181286,
+                            0.0001379946769624388,
+                            0.02078779689715382,
+                            0.033237241571263176,
+                            31.36068872331705,
                         ],
                         linf=[
-                            0.0016677536924397662,
-                            0.15617627577194781,
-                            0.3368120143100738,
-                            331.40871663423604,
+                            0.0016286690573188434,
+                            0.15623770697198225,
+                            0.3341371832270615,
+                            331.5373488726036,
                         ],
                         tspan=(0.0, 10.0),
                         initial_refinement_level=4)

From 0c97c9f64118420d36c6df7b348039af14d68cb6 Mon Sep 17 00:00:00 2001
From: Benedict <135045760+bgeihe@users.noreply.github.com>
Date: Wed, 10 Jan 2024 09:05:53 +0100
Subject: [PATCH 21/30] Update src/equations/compressible_euler_2d.jl

Co-authored-by: Andrew Winters <andrew.ross.winters@liu.se>
---
 src/equations/compressible_euler_2d.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/equations/compressible_euler_2d.jl b/src/equations/compressible_euler_2d.jl
index fa2c9f5e94f..2d100b23bf0 100644
--- a/src/equations/compressible_euler_2d.jl
+++ b/src/equations/compressible_euler_2d.jl
@@ -847,7 +847,7 @@ end
     p = 0.5 * (p_ll + p_rr) - 0.5 * c * rho * (v_rr - v_ll)
     v = 0.5 * (v_ll + v_rr) - 1 / (2 * c * rho) * (p_rr - p_ll)
 
-    # We treat the energy term in analogy to the potential temperature term in the paper by
+    # We treat the energy term analogous to the potential temperature term in the paper by
     # Chen et al., i.e. we use p_ll and p_rr, and not p
     if v >= 0
         f1, f2, f3, f4 = v * u_ll

From fc4da198db06bc01da53db43e801c40f912fd57a Mon Sep 17 00:00:00 2001
From: Benedict <135045760+bgeihe@users.noreply.github.com>
Date: Wed, 10 Jan 2024 09:06:38 +0100
Subject: [PATCH 22/30] Update src/equations/compressible_euler_2d.jl

Co-authored-by: Andrew Winters <andrew.ross.winters@liu.se>
---
 src/equations/compressible_euler_2d.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/equations/compressible_euler_2d.jl b/src/equations/compressible_euler_2d.jl
index 2d100b23bf0..3c6f759db2b 100644
--- a/src/equations/compressible_euler_2d.jl
+++ b/src/equations/compressible_euler_2d.jl
@@ -885,7 +885,7 @@ end
     p = 0.5 * (p_ll + p_rr) - 0.5 * c * rho * (v_rr - v_ll) / norm_
     v = 0.5 * (v_ll + v_rr) - 1 / (2 * c * rho) * (p_rr - p_ll) * norm_
 
-    # We treat the energy term in analogy to the potential temperature term in the paper by
+    # We treat the energy term analogous to the potential temperature term in the paper by
     # Chen et al., i.e. we use p_ll and p_rr, and not p
     if v >= 0
         f1, f2, f3, f4 = u_ll * v

From 8b8f1e21827482c5e1a0301bb36525e8a79bd2f2 Mon Sep 17 00:00:00 2001
From: Benedict <135045760+bgeihe@users.noreply.github.com>
Date: Wed, 10 Jan 2024 09:08:04 +0100
Subject: [PATCH 23/30] Update src/equations/compressible_euler_3d.jl

Co-authored-by: Andrew Winters <andrew.ross.winters@liu.se>
---
 src/equations/compressible_euler_3d.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/equations/compressible_euler_3d.jl b/src/equations/compressible_euler_3d.jl
index db74f1f5485..54b36e0b832 100644
--- a/src/equations/compressible_euler_3d.jl
+++ b/src/equations/compressible_euler_3d.jl
@@ -967,7 +967,7 @@ References:
     p = 0.5 * (p_ll + p_rr) - 0.5 * c * rho * (v_rr - v_ll)
     v = 0.5 * (v_ll + v_rr) - 1 / (2 * c * rho) * (p_rr - p_ll)
 
-    # We treat the energy term in analogy to the potential temperature term in the paper by
+    # We treat the energy term analogous to the potential temperature term in the paper by
     # Chen et al., i.e. we use p_ll and p_rr, and not p
     if v >= 0
         f1, f2, f3, f4, f5 = v * u_ll

From 854eefda0bc2955938c6f7b79d77c610706b0fac Mon Sep 17 00:00:00 2001
From: Benedict <135045760+bgeihe@users.noreply.github.com>
Date: Wed, 10 Jan 2024 09:08:19 +0100
Subject: [PATCH 24/30] Update src/equations/compressible_euler_3d.jl

Co-authored-by: Andrew Winters <andrew.ross.winters@liu.se>
---
 src/equations/compressible_euler_3d.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/equations/compressible_euler_3d.jl b/src/equations/compressible_euler_3d.jl
index 54b36e0b832..292b912f009 100644
--- a/src/equations/compressible_euler_3d.jl
+++ b/src/equations/compressible_euler_3d.jl
@@ -1009,7 +1009,7 @@ end
     p = 0.5 * (p_ll + p_rr) - 0.5 * c * rho * (v_rr - v_ll) / norm_
     v = 0.5 * (v_ll + v_rr) - 1 / (2 * c * rho) * (p_rr - p_ll) * norm_
 
-    # We treat the energy term in analogy to the potential temperature term in the paper by
+    # We treat the energy term analogous to the potential temperature term in the paper by
     # Chen et al., i.e. we use p_ll and p_rr, and not p
     if v >= 0
         f1, f2, f3, f4, f5 = v * u_ll

From 2d55dc9c60033aad03784fa01fd3313418cbae04 Mon Sep 17 00:00:00 2001
From: Benedict Geihe <bgeihe@uni-koeln.de>
Date: Wed, 10 Jan 2024 12:49:48 +0100
Subject: [PATCH 25/30] correct test result

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

diff --git a/test/test_tree_2d_euler.jl b/test/test_tree_2d_euler.jl
index 3e3f69d36c3..61b5c54b5e9 100644
--- a/test/test_tree_2d_euler.jl
+++ b/test/test_tree_2d_euler.jl
@@ -846,7 +846,7 @@ end
                             0.0016286690573188434,
                             0.15623770697198225,
                             0.3341371832270615,
-                            331.5373488726036,
+                            334.5373488726036,
                         ],
                         tspan=(0.0, 10.0),
                         initial_refinement_level=4)

From 54d7de00213ec096fd0a910709d076ddcfd54737 Mon Sep 17 00:00:00 2001
From: Benedict Geihe <bgeihe@uni-koeln.de>
Date: Wed, 10 Jan 2024 14:53:23 +0100
Subject: [PATCH 26/30] use callable struct to hold parameters

Thanks sloede!
---
 .../elixir_euler_warm_bubble.jl               | 35 ++++++++++++------
 .../tree_2d_dgsem/elixir_euler_warm_bubble.jl | 37 +++++++++++++------
 2 files changed, 48 insertions(+), 24 deletions(-)

diff --git a/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl b/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl
index be52f674a4e..ba1f23c8171 100644
--- a/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl
+++ b/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl
@@ -1,18 +1,27 @@
 using OrdinaryDiffEq
 using Trixi
 
-# Physical constants
-g::Float64 = 9.81            # gravity of earth
-c_p::Float64 = 1004.0        # heat capacity for constant pressure (dry air)
-c_v::Float64 = 717.0         # heat capacity for constant volume (dry air)
-gamma::Float64 = c_p / c_v   # heat capacity ratio (dry air)
-
 # Warm bubble test case from
 # Wicker, L. J., and Skamarock, W. C.
 # A time-splitting scheme for the elastic equations incorporating second-order Runge–Kutta
 # time differencing
 # [DOI: 10.1175/1520-0493(1998)126%3C1992:ATSSFT%3E2.0.CO;2](https://doi.org/10.1175/1520-0493(1998)126%3C1992:ATSSFT%3E2.0.CO;2)
-function initial_condition_warm_bubble(x, t, equations::CompressibleEulerEquations2D)
+struct WarmBubbleSetup
+    # Physical constants
+    g::Float64       # gravity of earth
+    c_p::Float64     # heat capacity for constant pressure (dry air)
+    c_v::Float64     # heat capacity for constant volume (dry air)
+    gamma::Float64   # heat capacity ratio (dry air)
+
+    function WarmBubbleSetup(; g = 9.81, c_p = 1004.0, c_v = 717.0, gamma = c_p / c_v)
+        new(g, c_p, c_v, gamma)
+    end
+end
+
+# Initial condition
+function (setup::WarmBubbleSetup)(x, t, equations::CompressibleEulerEquations2D)
+    @unpack g, c_p, c_v = setup
+
     # center of perturbation
     center_x = 10000.0
     center_z = 2000.0
@@ -49,16 +58,18 @@ function initial_condition_warm_bubble(x, t, equations::CompressibleEulerEquatio
     return SVector(rho, rho * v1, rho * v2, rho * E)
 end
 
-@inline function source_terms_gravity(u, x, t,
-                                      equations::CompressibleEulerEquations2D)
+# Source terms
+@inline function (setup::WarmBubbleSetup)(u, x, t, equations::CompressibleEulerEquations2D)
+    @unpack g = setup
     rho, _, rho_v2, _ = u
     return SVector(zero(eltype(u)), zero(eltype(u)), -g * rho, -g * rho_v2)
 end
 
 ###############################################################################
 # semidiscretization of the compressible Euler equations
+warm_bubble_setup = WarmBubbleSetup()
 
-equations = CompressibleEulerEquations2D(gamma)
+equations = CompressibleEulerEquations2D(warm_bubble_setup.gamma)
 
 boundary_conditions = (x_neg = boundary_condition_periodic,
                        x_pos = boundary_condition_periodic,
@@ -80,8 +91,8 @@ coordinates_max = (20_000.0, 10_000.0)
 cells_per_dimension = (64, 32)
 mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max)
 
-semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_warm_bubble, solver,
-                                    source_terms = source_terms_gravity,
+semi = SemidiscretizationHyperbolic(mesh, equations, warm_bubble_setup, solver,
+                                    source_terms = warm_bubble_setup,
                                     boundary_conditions = boundary_conditions)
 
 ###############################################################################
diff --git a/examples/tree_2d_dgsem/elixir_euler_warm_bubble.jl b/examples/tree_2d_dgsem/elixir_euler_warm_bubble.jl
index 15c6e6eb52b..1064da50232 100644
--- a/examples/tree_2d_dgsem/elixir_euler_warm_bubble.jl
+++ b/examples/tree_2d_dgsem/elixir_euler_warm_bubble.jl
@@ -1,18 +1,27 @@
 using OrdinaryDiffEq
 using Trixi
 
-# Physical constants
-g::Float64 = 9.81            # gravity of earth
-c_p::Float64 = 1004.0        # heat capacity for constant pressure (dry air)
-c_v::Float64 = 717.0         # heat capacity for constant volume (dry air)
-gamma::Float64 = c_p / c_v   # heat capacity ratio (dry air)
-
 # Warm bubble test case from
 # Wicker, L. J., and Skamarock, W. C.
 # A time-splitting scheme for the elastic equations incorporating second-order Runge–Kutta
 # time differencing
 # [DOI: 10.1175/1520-0493(1998)126%3C1992:ATSSFT%3E2.0.CO;2](https://doi.org/10.1175/1520-0493(1998)126%3C1992:ATSSFT%3E2.0.CO;2)
-function initial_condition_warm_bubble(x, t, equations::CompressibleEulerEquations2D)
+struct WarmBubbleSetup
+    # Physical constants
+    g::Float64       # gravity of earth
+    c_p::Float64     # heat capacity for constant pressure (dry air)
+    c_v::Float64     # heat capacity for constant volume (dry air)
+    gamma::Float64   # heat capacity ratio (dry air)
+
+    function WarmBubbleSetup(; g = 9.81, c_p = 1004.0, c_v = 717.0, gamma = c_p / c_v)
+        new(g, c_p, c_v, gamma)
+    end
+end
+
+# Initial condition
+function (setup::WarmBubbleSetup)(x, t, equations::CompressibleEulerEquations2D)
+    @unpack g, c_p, c_v = setup
+
     # center of perturbation
     center_x = 10000.0
     center_z = 2000.0
@@ -49,16 +58,18 @@ function initial_condition_warm_bubble(x, t, equations::CompressibleEulerEquatio
     return SVector(rho, rho * v1, rho * v2, rho * E)
 end
 
-@inline function source_terms_gravity(u, x, t,
-                                      equations::CompressibleEulerEquations2D)
+# Source terms
+@inline function (setup::WarmBubbleSetup)(u, x, t, equations::CompressibleEulerEquations2D)
+    @unpack g = setup
     rho, _, rho_v2, _ = u
     return SVector(zero(eltype(u)), zero(eltype(u)), -g * rho, -g * rho_v2)
 end
 
 ###############################################################################
 # semidiscretization of the compressible Euler equations
+warm_bubble_setup = WarmBubbleSetup()
 
-equations = CompressibleEulerEquations2D(gamma)
+equations = CompressibleEulerEquations2D(warm_bubble_setup.gamma)
 
 boundary_conditions = (x_neg = boundary_condition_periodic,
                        x_pos = boundary_condition_periodic,
@@ -77,13 +88,15 @@ solver = DGSEM(basis, surface_flux, volume_integral)
 coordinates_min = (0.0, 0.0)
 coordinates_max = (20_000.0, 10_000.0)
 
+# Same coordinates as in examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl
+# However TreeMesh will generate a 20_000 x 20_000 square domain instead
 mesh = TreeMesh(coordinates_min, coordinates_max,
                 initial_refinement_level = 6,
                 n_cells_max = 10_000,
                 periodicity = (true, false))
 
-semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_warm_bubble, solver,
-                                    source_terms = source_terms_gravity,
+semi = SemidiscretizationHyperbolic(mesh, equations, warm_bubble_setup, solver,
+                                    source_terms = warm_bubble_setup,
                                     boundary_conditions = boundary_conditions)
 
 ###############################################################################

From 843173b344921d7388fe701112ce6a1d1ce19472 Mon Sep 17 00:00:00 2001
From: Benedict <135045760+bgeihe@users.noreply.github.com>
Date: Wed, 10 Jan 2024 16:50:25 +0100
Subject: [PATCH 27/30] Update
 examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl

Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com>
---
 .../elixir_euler_warm_bubble.jl                  | 16 ++++++++++++----
 1 file changed, 12 insertions(+), 4 deletions(-)

diff --git a/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl b/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl
index ba1f23c8171..06d187568ae 100644
--- a/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl
+++ b/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl
@@ -2,10 +2,18 @@ using OrdinaryDiffEq
 using Trixi
 
 # Warm bubble test case from
-# Wicker, L. J., and Skamarock, W. C.
-# A time-splitting scheme for the elastic equations incorporating second-order Runge–Kutta
-# time differencing
-# [DOI: 10.1175/1520-0493(1998)126%3C1992:ATSSFT%3E2.0.CO;2](https://doi.org/10.1175/1520-0493(1998)126%3C1992:ATSSFT%3E2.0.CO;2)
+# - Wicker, L. J., and Skamarock, W. C.
+#   A time-splitting scheme for the elastic equations incorporating 
+#   second-order Runge–Kutta time differencing
+#   [DOI: 10.1175/1520-0493(1998)126%3C1992:ATSSFT%3E2.0.CO;2](https://doi.org/10.1175/1520-0493(1998)126%3C1992:ATSSFT%3E2.0.CO;2)
+# See also
+# - Bryan and Fritsch (2002)
+#   A Benchmark Simulation for Moist Nonhydrostatic Numerical Models
+#   [DOI: 10.1175/1520-0493(2002)130<2917:ABSFMN>2.0.CO;2](https://doi.org/10.1175/1520-0493(2002)130<2917:ABSFMN>2.0.CO;2)
+# - Carpenter, Droegemeier, Woodward, Hane (1990)
+#   Application of the Piecewise Parabolic Method (PPM) to 
+#   Meteorological Modeling
+#   [DOI: 10.1175/1520-0493(1990)118<0586:AOTPPM>2.0.CO;2](https://doi.org/10.1175/1520-0493(1990)118<0586:AOTPPM>2.0.CO;2)
 struct WarmBubbleSetup
     # Physical constants
     g::Float64       # gravity of earth

From 623dae344196b62f8ad43bc46769429cd85b26db Mon Sep 17 00:00:00 2001
From: Benedict <135045760+bgeihe@users.noreply.github.com>
Date: Wed, 10 Jan 2024 16:53:32 +0100
Subject: [PATCH 28/30] Update
 examples/tree_2d_dgsem/elixir_euler_warm_bubble.jl

Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com>
---
 .../tree_2d_dgsem/elixir_euler_warm_bubble.jl    | 16 ++++++++++++----
 1 file changed, 12 insertions(+), 4 deletions(-)

diff --git a/examples/tree_2d_dgsem/elixir_euler_warm_bubble.jl b/examples/tree_2d_dgsem/elixir_euler_warm_bubble.jl
index 1064da50232..867c27ab944 100644
--- a/examples/tree_2d_dgsem/elixir_euler_warm_bubble.jl
+++ b/examples/tree_2d_dgsem/elixir_euler_warm_bubble.jl
@@ -2,10 +2,18 @@ using OrdinaryDiffEq
 using Trixi
 
 # Warm bubble test case from
-# Wicker, L. J., and Skamarock, W. C.
-# A time-splitting scheme for the elastic equations incorporating second-order Runge–Kutta
-# time differencing
-# [DOI: 10.1175/1520-0493(1998)126%3C1992:ATSSFT%3E2.0.CO;2](https://doi.org/10.1175/1520-0493(1998)126%3C1992:ATSSFT%3E2.0.CO;2)
+# - Wicker, L. J., and Skamarock, W. C.
+#   A time-splitting scheme for the elastic equations incorporating 
+#   second-order Runge–Kutta time differencing
+#   [DOI: 10.1175/1520-0493(1998)126%3C1992:ATSSFT%3E2.0.CO;2](https://doi.org/10.1175/1520-0493(1998)126%3C1992:ATSSFT%3E2.0.CO;2)
+# See also
+# - Bryan and Fritsch (2002)
+#   A Benchmark Simulation for Moist Nonhydrostatic Numerical Models
+#   [DOI: 10.1175/1520-0493(2002)130<2917:ABSFMN>2.0.CO;2](https://doi.org/10.1175/1520-0493(2002)130<2917:ABSFMN>2.0.CO;2)
+# - Carpenter, Droegemeier, Woodward, Hane (1990)
+#   Application of the Piecewise Parabolic Method (PPM) to 
+#   Meteorological Modeling
+#   [DOI: 10.1175/1520-0493(1990)118<0586:AOTPPM>2.0.CO;2](https://doi.org/10.1175/1520-0493(1990)118<0586:AOTPPM>2.0.CO;2)
 struct WarmBubbleSetup
     # Physical constants
     g::Float64       # gravity of earth

From bed6975aa40cf9a03f947569ea1ef2d63283c695 Mon Sep 17 00:00:00 2001
From: Benedict Geihe <bgeihe@uni-koeln.de>
Date: Thu, 11 Jan 2024 11:33:06 +0100
Subject: [PATCH 29/30] year of Wicker paper, comment on tspan

[no ci]
---
 examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl | 4 ++--
 examples/tree_2d_dgsem/elixir_euler_warm_bubble.jl       | 4 ++--
 2 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl b/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl
index 06d187568ae..b26b3249e9e 100644
--- a/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl
+++ b/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl
@@ -2,7 +2,7 @@ using OrdinaryDiffEq
 using Trixi
 
 # Warm bubble test case from
-# - Wicker, L. J., and Skamarock, W. C.
+# - Wicker, L. J., and Skamarock, W. C. (1998)
 #   A time-splitting scheme for the elastic equations incorporating 
 #   second-order Runge–Kutta time differencing
 #   [DOI: 10.1175/1520-0493(1998)126%3C1992:ATSSFT%3E2.0.CO;2](https://doi.org/10.1175/1520-0493(1998)126%3C1992:ATSSFT%3E2.0.CO;2)
@@ -106,7 +106,7 @@ semi = SemidiscretizationHyperbolic(mesh, equations, warm_bubble_setup, solver,
 ###############################################################################
 # ODE solvers, callbacks etc.
 
-tspan = (0.0, 1000.0)
+tspan = (0.0, 1000.0)  # 1000 seconds final time
 
 ode = semidiscretize(semi, tspan)
 
diff --git a/examples/tree_2d_dgsem/elixir_euler_warm_bubble.jl b/examples/tree_2d_dgsem/elixir_euler_warm_bubble.jl
index 867c27ab944..e3054090291 100644
--- a/examples/tree_2d_dgsem/elixir_euler_warm_bubble.jl
+++ b/examples/tree_2d_dgsem/elixir_euler_warm_bubble.jl
@@ -2,7 +2,7 @@ using OrdinaryDiffEq
 using Trixi
 
 # Warm bubble test case from
-# - Wicker, L. J., and Skamarock, W. C.
+# - Wicker, L. J., and Skamarock, W. C. (1998)
 #   A time-splitting scheme for the elastic equations incorporating 
 #   second-order Runge–Kutta time differencing
 #   [DOI: 10.1175/1520-0493(1998)126%3C1992:ATSSFT%3E2.0.CO;2](https://doi.org/10.1175/1520-0493(1998)126%3C1992:ATSSFT%3E2.0.CO;2)
@@ -110,7 +110,7 @@ semi = SemidiscretizationHyperbolic(mesh, equations, warm_bubble_setup, solver,
 ###############################################################################
 # ODE solvers, callbacks etc.
 
-tspan = (0.0, 1000.0)
+tspan = (0.0, 1000.0)  # 1000 seconds final time
 
 ode = semidiscretize(semi, tspan)
 

From cc43a5074cfb74fc9105be036f75315e8bea7c37 Mon Sep 17 00:00:00 2001
From: Benedict Geihe <bgeihe@uni-koeln.de>
Date: Thu, 11 Jan 2024 18:12:42 +0100
Subject: [PATCH 30/30] add comment on speed of sound

---
 examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl | 3 +++
 examples/tree_2d_dgsem/elixir_euler_warm_bubble.jl       | 3 +++
 2 files changed, 6 insertions(+)

diff --git a/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl b/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl
index b26b3249e9e..05c09d57530 100644
--- a/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl
+++ b/examples/structured_2d_dgsem/elixir_euler_warm_bubble.jl
@@ -87,7 +87,10 @@ boundary_conditions = (x_neg = boundary_condition_periodic,
 polydeg = 3
 basis = LobattoLegendreBasis(polydeg)
 
+# This is a good estimate for the speed of sound in this example.
+# Other values between 300 and 400 should work as well.
 surface_flux = FluxLMARS(340.0)
+
 volume_flux = flux_kennedy_gruber
 volume_integral = VolumeIntegralFluxDifferencing(volume_flux)
 
diff --git a/examples/tree_2d_dgsem/elixir_euler_warm_bubble.jl b/examples/tree_2d_dgsem/elixir_euler_warm_bubble.jl
index e3054090291..f2e14273ae7 100644
--- a/examples/tree_2d_dgsem/elixir_euler_warm_bubble.jl
+++ b/examples/tree_2d_dgsem/elixir_euler_warm_bubble.jl
@@ -87,7 +87,10 @@ boundary_conditions = (x_neg = boundary_condition_periodic,
 polydeg = 3
 basis = LobattoLegendreBasis(polydeg)
 
+# This is a good estimate for the speed of sound in this example.
+# Other values between 300 and 400 should work as well.
 surface_flux = FluxLMARS(340.0)
+
 volume_flux = flux_kennedy_gruber
 volume_integral = VolumeIntegralFluxDifferencing(volume_flux)