From deb1fb6a5388302d8a27ed80a31e5981abadf162 Mon Sep 17 00:00:00 2001
From: KrisshChawla <krisshchawla@gmail.com>
Date: Fri, 24 Nov 2023 00:29:59 -0600
Subject: [PATCH 01/39] implementation of quasi 1d compressible Euler Equation

---
 src/equations/compressible_euler_quasi_1d.jl | 364 +++++++++++++++++++
 1 file changed, 364 insertions(+)
 create mode 100644 src/equations/compressible_euler_quasi_1d.jl

diff --git a/src/equations/compressible_euler_quasi_1d.jl b/src/equations/compressible_euler_quasi_1d.jl
new file mode 100644
index 00000000000..b8ca29fa0e3
--- /dev/null
+++ b/src/equations/compressible_euler_quasi_1d.jl
@@ -0,0 +1,364 @@
+# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
+# Since these FMAs can increase the performance of many numerical algorithms,
+# we need to opt-in explicitly.
+# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
+@muladd begin
+#! format: noindent
+
+@doc raw"""
+    CompressibleEulerEquations1D(gamma)
+
+The compressible Euler equations
+```math
+\frac{\partial}{\partial t}
+\begin{pmatrix}
+a \rho \\ a \rho v_1 \\ a e
+\end{pmatrix}
++
+\frac{\partial}{\partial x}
+\begin{pmatrix}
+a \rho v_1 \\ a \rho v_1^2 \\ a v_1 (e +p)
+\end{pmatrix}
++ 
+a \frac{\partial}{\partial x}
+\begin{pmatrix}
+0 \\ p \\ 0    
+\end{pmatrix}
+=
+\begin{pmatrix}
+0 \\ 0 \\ 0
+\end{pmatrix}
+```
+for an ideal gas with ratio of specific heats `gamma` in one space dimension.
+Here, ``\rho`` is the density, ``v_1`` the velocity, ``e`` the specific total energy **rather than** specific internal energy, ``a`` the (possibly) variable nozzle width, and
+```math
+p = (\gamma - 1) \left( e - \frac{1}{2} \rho v_1^2 \right)
+```
+the pressure.
+
+The nozzle width function ``a(x)`` is set inside the initial condition routine
+for a particular problem setup. To test the conservative form of the compressible Euler equations one can set the 
+nozzle width variable ``a`` to one. 
+
+In addition to the unknowns, Trixi.jl currently stores the nozzle width values at the approximation points 
+despite being fixed in time.
+This affects the implementation and use of these equations in various ways:
+* The flux values corresponding to the nozzle width must be zero.
+* The nozzle width values must be included when defining initial conditions, boundary conditions or
+  source terms.
+* [`AnalysisCallback`](@ref) analyzes this variable.
+* Trixi.jl's visualization tools will visualize the nozzle width by default.
+"""
+struct CompressibleEulerEquationsQuasi1D{RealT <: Real} <:
+       Trixi.AbstractCompressibleEulerEquations{1, 4}
+    gamma::RealT               # ratio of specific heats
+    inv_gamma_minus_one::RealT # = inv(gamma - 1); can be used to write slow divisions as fast multiplications
+
+    function CompressibleEulerEquationsQuasi1D(gamma)
+        γ, inv_gamma_minus_one = promote(gamma, inv(gamma - 1))
+        new{typeof(γ)}(γ, inv_gamma_minus_one)
+    end
+end
+
+have_nonconservative_terms(::CompressibleEulerEquationsQuasi1D) = Trixi.True()
+function varnames(::typeof(cons2cons), ::CompressibleEulerEquationsQuasi1D)
+    ("a_rho", "a_rho_v1", "a_e", "a")
+end
+varnames(::typeof(cons2prim), ::CompressibleEulerEquationsQuasi1D) = ("rho", "v1", "p", "a")
+
+"""
+    initial_condition_convergance_test(x, t, equations::CompressibleEulerEquationsQuasi1D)
+
+A smooth initial condition used for convergence tests in combination with
+[`source_terms_convergence_test`](@ref)
+(and [`BoundaryConditionDirichlet(initial_condition_convergence_test)`](@ref) in non-periodic domains).
+"""
+function initial_condition_convergance_test(x, t, equations::CompressibleEulerEquationsQuasi1D)
+    c = 2
+    A = 0.1
+    L = 2
+    f = 1 / L
+    ω = 2 * pi * f
+    ini = c + A * sin(ω * (x[1] - t))
+
+    rho = ini
+    v1 = 1.0
+    e = ini^2 / rho 
+    p = (equations.gamma - 1) * (e - 0.5 * rho * v1^2)
+    a = 1.5 - 0.5 * cos(x[1] * pi)
+    
+    return prim2cons(SVector(rho, v1, p, a), equations)
+end
+
+"""
+    source_terms_convergence_test(u, x, t, equations::CompressibleEulerEquationsQuasi1D)
+
+Source terms used for convergence tests in combination with
+[`initial_condition_convergence_test`](@ref)
+(and [`BoundaryConditionDirichlet(initial_condition_convergence_test)`](@ref) in non-periodic domains).
+
+This manufactured solution source term is specifically designed for the mozzle width 'a(x) = 1.5 - 0.5 * cos(x[1] * pi)'
+as defined in [`initial_condition_convergance_test`](@ref).
+"""
+@inline function source_terms_convergence_test(u, x, t,
+                                               equations::CompressibleEulerEquationsQuasi1D)
+    # Same settings as in `initial_condition`. 
+    #Derivatives calculated with ForwardDiff.jl
+    c = 2
+    A = 0.1
+    L = 2
+    f = 1 / L
+    ω = 2 * pi * f
+    x1, = x
+    ini(x1,t) = c + A * sin(ω * (x1 - t))
+
+    rho(x1,t) = ini(x1,t)
+    v1(x1, t) = 1.0
+    e(x1,t) = ini(x1,t)^2 / rho(x1,t)
+    p1(x1,t) = (equations.gamma - 1) * (e(x1,t) - 0.5 * rho(x1,t) * v1(x1,t)^2)
+    a(x1,t) =  1.5 - 0.5 * cos(x1 * pi)
+    
+    arho(x1,t) = a(x1,t) * rho(x1,t)
+    arhou(x1,t) = arho(x1,t) * v1(x1,t)
+    aE(x1,t) = a(x1,t) * e(x1,t)
+    
+    darho_dt(x1,t) = ForwardDiff.derivative(t->arho(x1,t), t)
+    darhou_dx(x1,t) = ForwardDiff.derivative(x1->arhou(x1,t), x1)
+    
+    arhouu(x1,t) = arhou(x1,t) * v1(x1,t)
+    darhou_dt(x1,t) = ForwardDiff.derivative(t->arhou(x1,t), t)
+    darhouu_dx(x1,t) = ForwardDiff.derivative(x1->arhouu(x1,t), x1)
+    dp1_dx(x1,t) = ForwardDiff.derivative(x1->p1(x1,t), x1)
+    
+    auEp(x1,t) = a(x1,t) * v1(x1,t) * (e(x1,t) + p1(x1,t))
+    daE_dt(x1,t) = ForwardDiff.derivative(t->aE(x1, t), t)
+    dauEp_dx(x1,t) = ForwardDiff.derivative(x1->auEp(x1,t), x1)
+
+
+    du1 = darho_dt(x1,t) + darhou_dx(x1,t)
+    du2 = darhou_dt(x1,t) + darhouu_dx(x1,t) + a(x1,t) * dp1_dx(x1,t) 
+    du3 = daE_dt(x1,t) + dauEp_dx(x1,t)
+    
+    return SVector(du1, du2, du3, 0.0)
+end
+
+"""
+    boundary_condition_slip_wall(u_inner, orientation, direction, x, t,
+                                 surface_flux_function, equations::CompressibleEulerEquations1D)
+Determine the boundary numerical surface flux for a slip wall condition.
+Imposes a zero normal velocity at the wall.
+Density is taken from the internal solution state and pressure is computed as an
+exact solution of a 1D Riemann problem. Further details about this boundary state
+are available in the paper:
+- J. J. W. van der Vegt and H. van der Ven (2002)
+  Slip flow boundary conditions in discontinuous Galerkin discretizations of
+  the Euler equations of gas dynamics
+  [PDF](https://reports.nlr.nl/bitstream/handle/10921/692/TP-2002-300.pdf?sequence=1)
+
+  Should be used together with [`TreeMesh`](@ref).
+"""
+@inline function boundary_condition_slip_wall(u_inner, orientation,
+                                              direction, x, t,
+                                              surface_flux_function,
+                                              equations::CompressibleEulerEquationsQuasi1D)
+    # compute the primitive variables
+    rho_local, v_normal, p_local, a_local = cons2prim(u_inner, equations)
+
+    if isodd(direction) # flip sign of normal to make it outward pointing
+        v_normal *= -1
+    end
+
+    # Get the solution of the pressure Riemann problem
+    # See Section 6.3.3 of
+    # Eleuterio F. Toro (2009)
+    # Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction
+    # [DOI: 10.1007/b79761](https://doi.org/10.1007/b79761)
+    if v_normal <= 0.0
+        sound_speed = sqrt(equations.gamma * p_local / rho_local) # local sound speed
+        p_star = p_local *
+                 (1 + 0.5 * (equations.gamma - 1) * v_normal / sound_speed)^(2 *
+                                                                             equations.gamma *
+                                                                             equations.inv_gamma_minus_one)
+    else # v_normal > 0.0
+        A = 2 / ((equations.gamma + 1) * rho_local)
+        B = p_local * (equations.gamma - 1) / (equations.gamma + 1)
+        p_star = p_local +
+                 0.5 * v_normal / A *
+                 (v_normal + sqrt(v_normal^2 + 4 * A * (p_local + B)))
+    end
+
+    # For the slip wall we directly set the flux as the normal velocity is zero
+    return SVector(zero(eltype(u_inner)),
+                   p_star,
+                   zero(eltype(u_inner)),
+                   zero(eltype(u_inner)))
+end
+
+# Calculate 1D flux for a single point
+@inline function Trixi.flux(u, orientation::Integer, equations::CompressibleEulerEquationsQuasi1D)
+    a_rho, a_rho_v1, a_e, a = u
+    rho, v1, p, a = cons2prim(u, equations)
+    e = a_e / a
+    
+    # Ignore orientation since it is always "1" in 1D
+    f1 = a_rho_v1
+    f2 = a_rho_v1 * v1 
+    f3 = a * v1 * (e + p)
+    
+    return SVector(f1, f2, f3, zero(eltype(u)))
+end
+
+"""
+@inline function flux_nonconservative_chan_etal(u_ll, u_rr, orientation::Integer,
+                                                equations::CompressibleEulerEquationsQuasi1D)
+
+Non-symmetric two-point volume flux discretizing the nonconservative (source) term
+that contains the gradient of the pressure  [`CompressibleEulerEquationsQuasi1D`](@ref) 
+and the nozzle width.
+
+Further details are available in the paper:
+- Jesse Chan, Khemraj Shukla, Xinhui Wu, Ruofeng Liu, Prani Nalluri (2023)
+    High order entropy stable schemes for the quasi-one-dimensional
+    shallow water and compressible Euler equations
+    [DOI: 10.48550/arXiv.2307.12089](https://doi.org/10.48550/arXiv.2307.12089)    
+"""
+@inline function flux_nonconservative_chan_etal(u_ll, u_rr, orientation::Integer,
+                                                equations::CompressibleEulerEquationsQuasi1D)
+    #Variables
+    _, _, p_ll, a_ll = cons2prim(u_ll, equations)
+    _, _, p_rr, _ = cons2prim(u_rr, equations)
+    
+    p_avg = 0.5 * (p_ll + p_rr)
+
+    z = zero(eltype(u_ll))
+    
+    return SVector(z, 2 * a_ll * p_avg, z, z)
+end
+
+"""
+@inline function flux_chan_etal(u_ll, u_rr, orientation::Integer,
+                                           equations::CompressibleEulerEquationsQuasi1D)
+
+Total energy conservative (mathematical entropy for quasi 1D compressible Euler equations) split form.
+
+Further details are available in the paper:
+- Jesse Chan, Khemraj Shukla, Xinhui Wu, Ruofeng Liu, Prani Nalluri (2023) 
+  High order entropy stable schemes for the quasi-one-dimensional
+  shallow water and compressible Euler equations
+  [DOI: 10.48550/arXiv.2307.12089](https://doi.org/10.48550/arXiv.2307.12089)     
+"""
+@inline function flux_chan_etal(u_ll, u_rr, orientation::Integer,
+                                           equations::CompressibleEulerEquationsQuasi1D)
+    # Unpack left and right state
+    rho_ll, v1_ll, p_ll, a_ll = cons2prim(u_ll, equations)
+    rho_rr, v1_rr, p_rr, a_rr = cons2prim(u_rr, equations)
+
+    # Compute the necessary mean values
+    rho_mean = Trixi.ln_mean(rho_ll, rho_rr)
+    # Algebraically equivalent to `inv_ln_mean(rho_ll / p_ll, rho_rr / p_rr)`
+    # in exact arithmetic since
+    #     log((ϱₗ/pₗ) / (ϱᵣ/pᵣ)) / (ϱₗ/pₗ - ϱᵣ/pᵣ)
+    #   = pₗ pᵣ log((ϱₗ pᵣ) / (ϱᵣ pₗ)) / (ϱₗ pᵣ - ϱᵣ pₗ)
+    inv_rho_p_mean = p_ll * p_rr * Trixi.inv_ln_mean(rho_ll * p_rr, rho_rr * p_ll)
+    v1_avg = 0.5 * (v1_ll + v1_rr)
+    a_v1_avg = 0.5 * (a_ll * v1_ll + a_rr * v1_rr) 
+    p_avg = 0.5 * (p_ll + p_rr)
+    velocity_square_avg = 0.5 * (v1_ll * v1_rr)
+
+    # Calculate fluxes
+    # Ignore orientation since it is always "1" in 1D
+    f1 = rho_mean * a_v1_avg
+    f2 = rho_mean * a_v1_avg * v1_avg
+    f3 = f1 * (velocity_square_avg + inv_rho_p_mean * equations.inv_gamma_minus_one) +
+         0.5 * (p_ll * a_rr * v1_rr + p_rr * a_ll * v1_ll)
+    
+    return SVector(f1, f2, f3, zero(eltype(u_ll)))
+end
+
+# Calculate estimates for maximum wave speed for local Lax-Friedrichs-type dissipation as the
+# maximum velocity magnitude plus the maximum speed of sound
+@inline function max_abs_speed_naive(u_ll, u_rr, orientation::Integer,
+                                     equations::CompressibleEulerEquationsQuasi1D)
+    a_rho_ll, a_rho_v1_ll, a_e_ll, a_ll = u_ll
+    a_rho_rr, a_rho_v1_rr, a_e_rr, a_rr = u_rr
+
+    # Calculate primitive variables and speed of sound
+    rho_ll = a_rho_ll / a_ll
+    e_ll = a_e_ll / a_ll
+    v1_ll = a_rho_v1_ll / a_rho_ll
+    v_mag_ll = abs(v1_ll)
+    p_ll = (equations.gamma - 1) * (e_ll - 0.5 * rho_ll * v_mag_ll^2)
+    c_ll = sqrt(equations.gamma * p_ll / rho_ll)
+    rho_rr = a_rho_rr / a_rr
+    e_rr = a_e_rr / a_rr
+    v1_rr = a_rho_v1_rr / a_rho_rr
+    v_mag_rr = abs(v1_rr)
+    p_rr = (equations.gamma - 1) * (e_rr - 0.5 * rho_rr * v_mag_rr^2)
+    c_rr = sqrt(equations.gamma * p_rr / rho_rr)
+
+    λ_max = max(v_mag_ll, v_mag_rr) + max(c_ll, c_rr)
+end
+
+@inline function Trixi.max_abs_speeds(u, equations::CompressibleEulerEquationsQuasi1D)
+    a_rho, a_rho_v1, a_e, a = u
+    rho = a_rho / a
+    v1 = a_rho_v1 / a_rho
+    e = a_e / a 
+    p = (equations.gamma - 1) * (e - 0.5 * rho * v1^2)
+    c = sqrt(equations.gamma * p / rho)
+
+    return (abs(v1) + c,)
+end
+
+# Convert conservative variables to primitive
+@inline function cons2prim(u, equations::CompressibleEulerEquationsQuasi1D)
+    a_rho, a_rho_v1, a_e, a = u    
+    q = cons2prim(u, CompressibleEulerEquations1D(equations.gamma))
+    
+    return SVector(q[1] / a, q[2], q[3] / a, a)
+end
+
+# Convert conservative variables to entropy
+@inline function Trixi.cons2entropy(u, equations::CompressibleEulerEquationsQuasi1D)
+    a_rho, a_rho_v1, a_e, a = u
+    q = cons2entropy(u, CompressibleEulerEquations1D(equations.gamma))
+        
+    w1 = q[1] -log(a)
+    w2 = q[2]
+    w3 = q[3]
+    
+    return SVector(w1, w2, w3, a)
+end
+
+# Convert primitive to conservative variables
+@inline function prim2cons(u, equations::CompressibleEulerEquationsQuasi1D)
+    rho, v1, p, a = u
+    q = prim2cons(u, CompressibleEulerEquations1D(equations.gamma))
+
+    return SVector(a * q[1], a * q[2], a * q[3] , a)
+end
+
+@inline function Trixi.density(u, equations::CompressibleEulerEquationsQuasi1D)
+    rho = u[1] / u[4]
+    return rho
+end
+
+@inline function Trixi.pressure(u, equations::CompressibleEulerEquationsQuasi1D)
+    a_rho, a_rho_v1, a_e, a = u
+    rho = a_rho / a
+    v1 = a_rho_v1 / a_rho
+    e = a_e / a 
+    p = (equations.gamma - 1) * (e - 0.5 * rho * v1^2)
+    return p
+end
+
+@inline function Trixi.density_pressure(u, equations::CompressibleEulerEquationsQuasi1D)
+    a_rho, a_rho_v1, a_e, a = u
+    rho = a_rho / a
+    v1 = a_rho_v1 / a_rho
+    rho_v1 = a_rho_v1 / a
+    e = a_e / a 
+    rho_times_p = (equations.gamma - 1) * (e - 0.5 * rho * v1^2) * rho
+    return rho_times_p
+end
+end # @muladd
\ No newline at end of file

From 62734c0469c2aeb888db3bec2742f81e9019af4a Mon Sep 17 00:00:00 2001
From: KrisshChawla <krisshchawla@gmail.com>
Date: Fri, 24 Nov 2023 00:47:09 -0600
Subject: [PATCH 02/39] added example elixir for quasi 1d compressible Euler
 equations

---
 .../elixir_euler_quasi_1d_source_terms.jl     | 60 +++++++++++++++++++
 1 file changed, 60 insertions(+)
 create mode 100644 examples/tree_1d_dgsem/elixir_euler_quasi_1d_source_terms.jl

diff --git a/examples/tree_1d_dgsem/elixir_euler_quasi_1d_source_terms.jl b/examples/tree_1d_dgsem/elixir_euler_quasi_1d_source_terms.jl
new file mode 100644
index 00000000000..964519aba9a
--- /dev/null
+++ b/examples/tree_1d_dgsem/elixir_euler_quasi_1d_source_terms.jl
@@ -0,0 +1,60 @@
+using OrdinaryDiffEq
+using Trixi
+using ForwardDiff
+
+###############################################################################
+# Semidiscretization of the quasi 1d compressible Euler equations
+# See Chan et al.  https://doi.org/10.48550/arXiv.2307.12089 for details
+
+equations = CompressibleEulerEquationsQuasi1D(1.4)
+
+initial_condition = initial_condition_convergance_test
+
+volume_flux = (flux_chan_etal, flux_nonconservative_chan_etal)
+surface_flux = volume_flux
+solver = DGSEM(polydeg=4, surface_flux=surface_flux,
+               volume_integral=VolumeIntegralFluxDifferencing(volume_flux))
+
+coordinates_min = -1.0
+coordinates_max = 1.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
\ No newline at end of file

From 34af78221c56bb4df017f153603cd62eec5a4d91 Mon Sep 17 00:00:00 2001
From: KrisshChawla <krisshchawla@gmail.com>
Date: Fri, 24 Nov 2023 00:58:03 -0600
Subject: [PATCH 03/39] added example elixir with a discontinuous initial
 condition

---
 .../elixir_euler_quasi_1d_discontinuous.jl    | 85 +++++++++++++++++++
 1 file changed, 85 insertions(+)
 create mode 100644 examples/tree_1d_dgsem/elixir_euler_quasi_1d_discontinuous.jl

diff --git a/examples/tree_1d_dgsem/elixir_euler_quasi_1d_discontinuous.jl b/examples/tree_1d_dgsem/elixir_euler_quasi_1d_discontinuous.jl
new file mode 100644
index 00000000000..c4e3c9cf0bf
--- /dev/null
+++ b/examples/tree_1d_dgsem/elixir_euler_quasi_1d_discontinuous.jl
@@ -0,0 +1,85 @@
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# Semidiscretization of the quasi 1d compressible Euler equations
+# See Chan et al.  https://doi.org/10.48550/arXiv.2307.12089 for details
+
+equations = CompressibleEulerEquationsQuasi1D(1.4)
+
+"""
+    initial_condition_discontinuity(x, t, equations::CompressibleEulerEquations1D)
+
+An entropy conservation verification initial condition taken from
+- Jesse Chan, Khemraj Shukla, Xinhui Wu, Ruofeng Liu, Prani Nalluri (2023)
+    High order entropy stable schemes for the quasi-one-dimensional
+    shallow water and compressible Euler equations
+    [DOI: 10.48550/arXiv.2307.12089](https://doi.org/10.48550/arXiv.2307.12089)   
+"""
+function initial_condition_discontinuity(x, t,
+                                            equations::CompressibleEulerEquationsQuasi1D)
+    rho = (x[1] < 0) ? 3.4718 : 2.0
+    v1 = (x[1] < 0) ? -2.5923 : -3.0
+    p = (x[1] < 0) ? 5.7118 : 2.639
+    a = (x[1] < 0) ? 1.0 : 1.5
+    
+    return prim2cons(SVector(rho, v1, p, a), equations)
+end
+
+initial_condition = initial_condition_discontinuity
+
+volume_flux = (flux_chan_etal, flux_nonconservative_chan_etal)
+surface_flux = (flux_lax_friedrichs, flux_nonconservative_chan_etal)
+
+basis = LobattoLegendreBasis(3)
+indicator_sc = IndicatorHennemannGassner(equations, basis,
+                                         alpha_max = 0.5,
+                                         alpha_min = 0.001,
+                                         alpha_smooth = true,
+                                         variable = density_pressure)
+volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
+                                                 volume_flux_dg = volume_flux,
+                                                 volume_flux_fv = surface_flux)
+solver = DGSEM(basis, surface_flux, volume_integral)
+
+coordinates_min = (-1.0,)
+coordinates_max = (1.0,)
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 6,
+                n_cells_max = 10_000)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+###############################################################################
+# 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)
+
+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.5)
+
+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
\ No newline at end of file

From bb3815dcf69f89619d0856d40532be72897b701a Mon Sep 17 00:00:00 2001
From: KrisshChawla <krisshchawla@gmail.com>
Date: Fri, 24 Nov 2023 10:12:00 -0600
Subject: [PATCH 04/39] including and exported
 CompressibleEulerEquationsQuasi1D

---
 src/Trixi.jl               |  1 +
 src/equations/equations.jl |  1 +
 test/test_tree_1d_euler.jl | 49 ++++++++++++++++++++++++++++++++++++++
 3 files changed, 51 insertions(+)

diff --git a/src/Trixi.jl b/src/Trixi.jl
index b8110cf5bdd..e7b849e2642 100644
--- a/src/Trixi.jl
+++ b/src/Trixi.jl
@@ -139,6 +139,7 @@ export AcousticPerturbationEquations2D,
        CompressibleEulerEquations3D,
        CompressibleEulerMulticomponentEquations1D,
        CompressibleEulerMulticomponentEquations2D,
+       CompressibleEulerEquationsQuasi1D,
        IdealGlmMhdEquations1D, IdealGlmMhdEquations2D, IdealGlmMhdEquations3D,
        IdealGlmMhdMulticomponentEquations1D, IdealGlmMhdMulticomponentEquations2D,
        HyperbolicDiffusionEquations1D, HyperbolicDiffusionEquations2D,
diff --git a/src/equations/equations.jl b/src/equations/equations.jl
index ba2ad1cd1cd..efbcd40d434 100644
--- a/src/equations/equations.jl
+++ b/src/equations/equations.jl
@@ -402,6 +402,7 @@ abstract type AbstractCompressibleEulerEquations{NDIMS, NVARS} <:
 include("compressible_euler_1d.jl")
 include("compressible_euler_2d.jl")
 include("compressible_euler_3d.jl")
+include("compressible_euler_quasi_1d.jl")
 
 # CompressibleEulerMulticomponentEquations
 abstract type AbstractCompressibleEulerMulticomponentEquations{NDIMS, NVARS, NCOMP} <:
diff --git a/test/test_tree_1d_euler.jl b/test/test_tree_1d_euler.jl
index 6cd7998ab02..5dc81da0261 100644
--- a/test/test_tree_1d_euler.jl
+++ b/test/test_tree_1d_euler.jl
@@ -393,6 +393,55 @@ end
         @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
     end
 end
+
+@trixi_testset "elixir_euler_quasi_1d_source_terms.jl" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_quasi_1d_source_terms.jl"),
+                        l2=[
+                            3.876288364589425e-7,
+                            2.224704317890825e-7,
+                            2.964004223514882e-7,
+                            5.2716983399807875e-8,
+                        ],
+                        linf=[
+                            2.392511842419509e-6,
+                            1.360369332736866e-6,
+                            1.8218888442333991e-6,
+                            1.1166012159335992e-7,
+                        ])
+    # 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_euler_quasi_1d_discontinuous.jl" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR,
+                                 "elixir_euler_quasi_1d_discontinuous.jl"),
+                        l2=[
+                            0.04551042115634054,
+                            0.03675058478890224,
+                            0.24689859591319177,
+                            0.0368449418082899,
+                        ],
+                        linf=[
+                            0.3313374853025661,
+                            0.11621933362147452,
+                            1.8274030135685742,
+                            0.28045939999015945,
+                        ])
+    # 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

From f5c15beed686afc3a0631e49f96289f7f2595fd7 Mon Sep 17 00:00:00 2001
From: KrisshChawla <krisshchawla@gmail.com>
Date: Fri, 24 Nov 2023 10:24:21 -0600
Subject: [PATCH 05/39] formatting

---
 .../elixir_euler_quasi_1d_discontinuous.jl    |   6 +-
 .../elixir_euler_quasi_1d_source_terms.jl     |   8 +-
 src/equations/compressible_euler_quasi_1d.jl  | 114 +++++++++---------
 3 files changed, 66 insertions(+), 62 deletions(-)

diff --git a/examples/tree_1d_dgsem/elixir_euler_quasi_1d_discontinuous.jl b/examples/tree_1d_dgsem/elixir_euler_quasi_1d_discontinuous.jl
index c4e3c9cf0bf..8f0196b3d2e 100644
--- a/examples/tree_1d_dgsem/elixir_euler_quasi_1d_discontinuous.jl
+++ b/examples/tree_1d_dgsem/elixir_euler_quasi_1d_discontinuous.jl
@@ -17,12 +17,12 @@ An entropy conservation verification initial condition taken from
     [DOI: 10.48550/arXiv.2307.12089](https://doi.org/10.48550/arXiv.2307.12089)   
 """
 function initial_condition_discontinuity(x, t,
-                                            equations::CompressibleEulerEquationsQuasi1D)
+                                         equations::CompressibleEulerEquationsQuasi1D)
     rho = (x[1] < 0) ? 3.4718 : 2.0
     v1 = (x[1] < 0) ? -2.5923 : -3.0
     p = (x[1] < 0) ? 5.7118 : 2.639
     a = (x[1] < 0) ? 1.0 : 1.5
-    
+
     return prim2cons(SVector(rho, v1, p, a), equations)
 end
 
@@ -82,4 +82,4 @@ 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
\ No newline at end of file
+summary_callback() # print the timer summary
diff --git a/examples/tree_1d_dgsem/elixir_euler_quasi_1d_source_terms.jl b/examples/tree_1d_dgsem/elixir_euler_quasi_1d_source_terms.jl
index 964519aba9a..7b08a455892 100644
--- a/examples/tree_1d_dgsem/elixir_euler_quasi_1d_source_terms.jl
+++ b/examples/tree_1d_dgsem/elixir_euler_quasi_1d_source_terms.jl
@@ -12,8 +12,8 @@ initial_condition = initial_condition_convergance_test
 
 volume_flux = (flux_chan_etal, flux_nonconservative_chan_etal)
 surface_flux = volume_flux
-solver = DGSEM(polydeg=4, surface_flux=surface_flux,
-               volume_integral=VolumeIntegralFluxDifferencing(volume_flux))
+solver = DGSEM(polydeg = 4, surface_flux = surface_flux,
+               volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
 
 coordinates_min = -1.0
 coordinates_max = 1.0
@@ -22,7 +22,7 @@ mesh = TreeMesh(coordinates_min, coordinates_max,
                 n_cells_max = 10_000)
 
 semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
-                                    source_terms=source_terms_convergence_test)
+                                    source_terms = source_terms_convergence_test)
 
 ###############################################################################
 # ODE solvers, callbacks etc.
@@ -57,4 +57,4 @@ 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
\ No newline at end of file
+summary_callback() # print the timer summary
diff --git a/src/equations/compressible_euler_quasi_1d.jl b/src/equations/compressible_euler_quasi_1d.jl
index b8ca29fa0e3..e3361629b7a 100644
--- a/src/equations/compressible_euler_quasi_1d.jl
+++ b/src/equations/compressible_euler_quasi_1d.jl
@@ -6,9 +6,9 @@
 #! format: noindent
 
 @doc raw"""
-    CompressibleEulerEquations1D(gamma)
+    CompressibleEulerEquationsQuasi1D(gamma)
 
-The compressible Euler equations
+The quasi-1d compressible Euler equations
 ```math
 \frac{\partial}{\partial t}
 \begin{pmatrix}
@@ -30,7 +30,8 @@ a \frac{\partial}{\partial x}
 \end{pmatrix}
 ```
 for an ideal gas with ratio of specific heats `gamma` in one space dimension.
-Here, ``\rho`` is the density, ``v_1`` the velocity, ``e`` the specific total energy **rather than** specific internal energy, ``a`` the (possibly) variable nozzle width, and
+Here, ``\rho`` is the density, ``v_1`` the velocity, ``e`` the specific total energy **rather than** specific internal energy, 
+``a`` the (possibly) variable nozzle width, and
 ```math
 p = (\gamma - 1) \left( e - \frac{1}{2} \rho v_1^2 \right)
 ```
@@ -64,7 +65,9 @@ have_nonconservative_terms(::CompressibleEulerEquationsQuasi1D) = Trixi.True()
 function varnames(::typeof(cons2cons), ::CompressibleEulerEquationsQuasi1D)
     ("a_rho", "a_rho_v1", "a_e", "a")
 end
-varnames(::typeof(cons2prim), ::CompressibleEulerEquationsQuasi1D) = ("rho", "v1", "p", "a")
+function varnames(::typeof(cons2prim), ::CompressibleEulerEquationsQuasi1D)
+    ("rho", "v1", "p", "a")
+end
 
 """
     initial_condition_convergance_test(x, t, equations::CompressibleEulerEquationsQuasi1D)
@@ -73,7 +76,8 @@ A smooth initial condition used for convergence tests in combination with
 [`source_terms_convergence_test`](@ref)
 (and [`BoundaryConditionDirichlet(initial_condition_convergence_test)`](@ref) in non-periodic domains).
 """
-function initial_condition_convergance_test(x, t, equations::CompressibleEulerEquationsQuasi1D)
+function initial_condition_convergance_test(x, t,
+                                            equations::CompressibleEulerEquationsQuasi1D)
     c = 2
     A = 0.1
     L = 2
@@ -83,10 +87,10 @@ function initial_condition_convergance_test(x, t, equations::CompressibleEulerEq
 
     rho = ini
     v1 = 1.0
-    e = ini^2 / rho 
+    e = ini^2 / rho
     p = (equations.gamma - 1) * (e - 0.5 * rho * v1^2)
     a = 1.5 - 0.5 * cos(x[1] * pi)
-    
+
     return prim2cons(SVector(rho, v1, p, a), equations)
 end
 
@@ -102,7 +106,7 @@ as defined in [`initial_condition_convergance_test`](@ref).
 """
 @inline function source_terms_convergence_test(u, x, t,
                                                equations::CompressibleEulerEquationsQuasi1D)
-    # Same settings as in `initial_condition`. 
+    # Same settings as in `initial_condition_convergance_test`. 
     #Derivatives calculated with ForwardDiff.jl
     c = 2
     A = 0.1
@@ -110,41 +114,40 @@ as defined in [`initial_condition_convergance_test`](@ref).
     f = 1 / L
     ω = 2 * pi * f
     x1, = x
-    ini(x1,t) = c + A * sin(ω * (x1 - t))
+    ini(x1, t) = c + A * sin(ω * (x1 - t))
 
-    rho(x1,t) = ini(x1,t)
+    rho(x1, t) = ini(x1, t)
     v1(x1, t) = 1.0
-    e(x1,t) = ini(x1,t)^2 / rho(x1,t)
-    p1(x1,t) = (equations.gamma - 1) * (e(x1,t) - 0.5 * rho(x1,t) * v1(x1,t)^2)
-    a(x1,t) =  1.5 - 0.5 * cos(x1 * pi)
-    
-    arho(x1,t) = a(x1,t) * rho(x1,t)
-    arhou(x1,t) = arho(x1,t) * v1(x1,t)
-    aE(x1,t) = a(x1,t) * e(x1,t)
-    
-    darho_dt(x1,t) = ForwardDiff.derivative(t->arho(x1,t), t)
-    darhou_dx(x1,t) = ForwardDiff.derivative(x1->arhou(x1,t), x1)
-    
-    arhouu(x1,t) = arhou(x1,t) * v1(x1,t)
-    darhou_dt(x1,t) = ForwardDiff.derivative(t->arhou(x1,t), t)
-    darhouu_dx(x1,t) = ForwardDiff.derivative(x1->arhouu(x1,t), x1)
-    dp1_dx(x1,t) = ForwardDiff.derivative(x1->p1(x1,t), x1)
-    
-    auEp(x1,t) = a(x1,t) * v1(x1,t) * (e(x1,t) + p1(x1,t))
-    daE_dt(x1,t) = ForwardDiff.derivative(t->aE(x1, t), t)
-    dauEp_dx(x1,t) = ForwardDiff.derivative(x1->auEp(x1,t), x1)
-
-
-    du1 = darho_dt(x1,t) + darhou_dx(x1,t)
-    du2 = darhou_dt(x1,t) + darhouu_dx(x1,t) + a(x1,t) * dp1_dx(x1,t) 
-    du3 = daE_dt(x1,t) + dauEp_dx(x1,t)
-    
+    e(x1, t) = ini(x1, t)^2 / rho(x1, t)
+    p1(x1, t) = (equations.gamma - 1) * (e(x1, t) - 0.5 * rho(x1, t) * v1(x1, t)^2)
+    a(x1, t) = 1.5 - 0.5 * cos(x1 * pi)
+
+    arho(x1, t) = a(x1, t) * rho(x1, t)
+    arhou(x1, t) = arho(x1, t) * v1(x1, t)
+    aE(x1, t) = a(x1, t) * e(x1, t)
+
+    darho_dt(x1, t) = ForwardDiff.derivative(t -> arho(x1, t), t)
+    darhou_dx(x1, t) = ForwardDiff.derivative(x1 -> arhou(x1, t), x1)
+
+    arhouu(x1, t) = arhou(x1, t) * v1(x1, t)
+    darhou_dt(x1, t) = ForwardDiff.derivative(t -> arhou(x1, t), t)
+    darhouu_dx(x1, t) = ForwardDiff.derivative(x1 -> arhouu(x1, t), x1)
+    dp1_dx(x1, t) = ForwardDiff.derivative(x1 -> p1(x1, t), x1)
+
+    auEp(x1, t) = a(x1, t) * v1(x1, t) * (e(x1, t) + p1(x1, t))
+    daE_dt(x1, t) = ForwardDiff.derivative(t -> aE(x1, t), t)
+    dauEp_dx(x1, t) = ForwardDiff.derivative(x1 -> auEp(x1, t), x1)
+
+    du1 = darho_dt(x1, t) + darhou_dx(x1, t)
+    du2 = darhou_dt(x1, t) + darhouu_dx(x1, t) + a(x1, t) * dp1_dx(x1, t)
+    du3 = daE_dt(x1, t) + dauEp_dx(x1, t)
+
     return SVector(du1, du2, du3, 0.0)
 end
 
 """
     boundary_condition_slip_wall(u_inner, orientation, direction, x, t,
-                                 surface_flux_function, equations::CompressibleEulerEquations1D)
+                                 surface_flux_function, equations::CompressibleEulerEquationsQuasi1D)
 Determine the boundary numerical surface flux for a slip wall condition.
 Imposes a zero normal velocity at the wall.
 Density is taken from the internal solution state and pressure is computed as an
@@ -195,16 +198,17 @@ are available in the paper:
 end
 
 # Calculate 1D flux for a single point
-@inline function Trixi.flux(u, orientation::Integer, equations::CompressibleEulerEquationsQuasi1D)
+@inline function Trixi.flux(u, orientation::Integer,
+                            equations::CompressibleEulerEquationsQuasi1D)
     a_rho, a_rho_v1, a_e, a = u
     rho, v1, p, a = cons2prim(u, equations)
     e = a_e / a
-    
+
     # Ignore orientation since it is always "1" in 1D
     f1 = a_rho_v1
-    f2 = a_rho_v1 * v1 
+    f2 = a_rho_v1 * v1
     f3 = a * v1 * (e + p)
-    
+
     return SVector(f1, f2, f3, zero(eltype(u)))
 end
 
@@ -227,11 +231,11 @@ Further details are available in the paper:
     #Variables
     _, _, p_ll, a_ll = cons2prim(u_ll, equations)
     _, _, p_rr, _ = cons2prim(u_rr, equations)
-    
+
     p_avg = 0.5 * (p_ll + p_rr)
 
     z = zero(eltype(u_ll))
-    
+
     return SVector(z, 2 * a_ll * p_avg, z, z)
 end
 
@@ -248,7 +252,7 @@ Further details are available in the paper:
   [DOI: 10.48550/arXiv.2307.12089](https://doi.org/10.48550/arXiv.2307.12089)     
 """
 @inline function flux_chan_etal(u_ll, u_rr, orientation::Integer,
-                                           equations::CompressibleEulerEquationsQuasi1D)
+                                equations::CompressibleEulerEquationsQuasi1D)
     # Unpack left and right state
     rho_ll, v1_ll, p_ll, a_ll = cons2prim(u_ll, equations)
     rho_rr, v1_rr, p_rr, a_rr = cons2prim(u_rr, equations)
@@ -261,7 +265,7 @@ Further details are available in the paper:
     #   = pₗ pᵣ log((ϱₗ pᵣ) / (ϱᵣ pₗ)) / (ϱₗ pᵣ - ϱᵣ pₗ)
     inv_rho_p_mean = p_ll * p_rr * Trixi.inv_ln_mean(rho_ll * p_rr, rho_rr * p_ll)
     v1_avg = 0.5 * (v1_ll + v1_rr)
-    a_v1_avg = 0.5 * (a_ll * v1_ll + a_rr * v1_rr) 
+    a_v1_avg = 0.5 * (a_ll * v1_ll + a_rr * v1_rr)
     p_avg = 0.5 * (p_ll + p_rr)
     velocity_square_avg = 0.5 * (v1_ll * v1_rr)
 
@@ -271,7 +275,7 @@ Further details are available in the paper:
     f2 = rho_mean * a_v1_avg * v1_avg
     f3 = f1 * (velocity_square_avg + inv_rho_p_mean * equations.inv_gamma_minus_one) +
          0.5 * (p_ll * a_rr * v1_rr + p_rr * a_ll * v1_ll)
-    
+
     return SVector(f1, f2, f3, zero(eltype(u_ll)))
 end
 
@@ -303,7 +307,7 @@ end
     a_rho, a_rho_v1, a_e, a = u
     rho = a_rho / a
     v1 = a_rho_v1 / a_rho
-    e = a_e / a 
+    e = a_e / a
     p = (equations.gamma - 1) * (e - 0.5 * rho * v1^2)
     c = sqrt(equations.gamma * p / rho)
 
@@ -312,9 +316,9 @@ end
 
 # Convert conservative variables to primitive
 @inline function cons2prim(u, equations::CompressibleEulerEquationsQuasi1D)
-    a_rho, a_rho_v1, a_e, a = u    
+    a_rho, a_rho_v1, a_e, a = u
     q = cons2prim(u, CompressibleEulerEquations1D(equations.gamma))
-    
+
     return SVector(q[1] / a, q[2], q[3] / a, a)
 end
 
@@ -322,11 +326,11 @@ end
 @inline function Trixi.cons2entropy(u, equations::CompressibleEulerEquationsQuasi1D)
     a_rho, a_rho_v1, a_e, a = u
     q = cons2entropy(u, CompressibleEulerEquations1D(equations.gamma))
-        
-    w1 = q[1] -log(a)
+
+    w1 = q[1] - log(a)
     w2 = q[2]
     w3 = q[3]
-    
+
     return SVector(w1, w2, w3, a)
 end
 
@@ -335,7 +339,7 @@ end
     rho, v1, p, a = u
     q = prim2cons(u, CompressibleEulerEquations1D(equations.gamma))
 
-    return SVector(a * q[1], a * q[2], a * q[3] , a)
+    return SVector(a * q[1], a * q[2], a * q[3], a)
 end
 
 @inline function Trixi.density(u, equations::CompressibleEulerEquationsQuasi1D)
@@ -347,7 +351,7 @@ end
     a_rho, a_rho_v1, a_e, a = u
     rho = a_rho / a
     v1 = a_rho_v1 / a_rho
-    e = a_e / a 
+    e = a_e / a
     p = (equations.gamma - 1) * (e - 0.5 * rho * v1^2)
     return p
 end
@@ -357,8 +361,8 @@ end
     rho = a_rho / a
     v1 = a_rho_v1 / a_rho
     rho_v1 = a_rho_v1 / a
-    e = a_e / a 
+    e = a_e / a
     rho_times_p = (equations.gamma - 1) * (e - 0.5 * rho * v1^2) * rho
     return rho_times_p
 end
-end # @muladd
\ No newline at end of file
+end # @muladd

From 75c8dcab5f731396cbf23897d01ff1937f07c289 Mon Sep 17 00:00:00 2001
From: KrisshChawla <krisshchawla@gmail.com>
Date: Fri, 24 Nov 2023 17:21:44 -0600
Subject: [PATCH 06/39] added entropy conservative test

---
 .../tree_1d_dgsem/elixir_euler_quasi_1d_ec.jl | 68 +++++++++++++++++++
 1 file changed, 68 insertions(+)
 create mode 100644 examples/tree_1d_dgsem/elixir_euler_quasi_1d_ec.jl

diff --git a/examples/tree_1d_dgsem/elixir_euler_quasi_1d_ec.jl b/examples/tree_1d_dgsem/elixir_euler_quasi_1d_ec.jl
new file mode 100644
index 00000000000..81ca5c13023
--- /dev/null
+++ b/examples/tree_1d_dgsem/elixir_euler_quasi_1d_ec.jl
@@ -0,0 +1,68 @@
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# Semidiscretization of the quasi 1d compressible Euler equations
+# See Chan et al.  https://doi.org/10.48550/arXiv.2307.12089 for details
+
+equations = CompressibleEulerEquationsQuasi1D(1.4)
+
+function initial_condition_ec(x, t, equations::CompressibleEulerEquationsQuasi1D)
+    v1 = 0.1
+    rho = 2.0 + 0.1 * x[1]
+    p = 3.0
+    a = 2.0 + x[1]
+    
+    return prim2cons(SVector(rho, v1, p, a), equations)
+end
+
+initial_condition = initial_condition_ec
+
+volume_flux = (flux_chan_etal, flux_nonconservative_chan_etal)
+surface_flux = volume_flux
+
+basis = LobattoLegendreBasis(3)
+solver = DGSEM(polydeg = 4, surface_flux = surface_flux,
+               volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+
+coordinates_min = (-1.0,)
+coordinates_max = (1.0,)
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 6,
+                n_cells_max = 10_000)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+###############################################################################
+# 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)
+
+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.5)
+
+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
\ No newline at end of file

From 0ef0305de95b7c1e1d4e1a9779c143bcf42dd096 Mon Sep 17 00:00:00 2001
From: KrisshChawla <krisshchawla@gmail.com>
Date: Fri, 24 Nov 2023 17:27:42 -0600
Subject: [PATCH 07/39] fixing spelling

---
 .../tree_1d_dgsem/elixir_euler_quasi_1d_source_terms.jl   | 2 +-
 src/equations/compressible_euler_quasi_1d.jl              | 8 ++++----
 2 files changed, 5 insertions(+), 5 deletions(-)

diff --git a/examples/tree_1d_dgsem/elixir_euler_quasi_1d_source_terms.jl b/examples/tree_1d_dgsem/elixir_euler_quasi_1d_source_terms.jl
index 7b08a455892..3f49e2718f5 100644
--- a/examples/tree_1d_dgsem/elixir_euler_quasi_1d_source_terms.jl
+++ b/examples/tree_1d_dgsem/elixir_euler_quasi_1d_source_terms.jl
@@ -8,7 +8,7 @@ using ForwardDiff
 
 equations = CompressibleEulerEquationsQuasi1D(1.4)
 
-initial_condition = initial_condition_convergance_test
+initial_condition = initial_condition_convergence_test
 
 volume_flux = (flux_chan_etal, flux_nonconservative_chan_etal)
 surface_flux = volume_flux
diff --git a/src/equations/compressible_euler_quasi_1d.jl b/src/equations/compressible_euler_quasi_1d.jl
index e3361629b7a..34e6cafe2dc 100644
--- a/src/equations/compressible_euler_quasi_1d.jl
+++ b/src/equations/compressible_euler_quasi_1d.jl
@@ -70,13 +70,13 @@ function varnames(::typeof(cons2prim), ::CompressibleEulerEquationsQuasi1D)
 end
 
 """
-    initial_condition_convergance_test(x, t, equations::CompressibleEulerEquationsQuasi1D)
+    initial_condition_convergence_test(x, t, equations::CompressibleEulerEquationsQuasi1D)
 
 A smooth initial condition used for convergence tests in combination with
 [`source_terms_convergence_test`](@ref)
 (and [`BoundaryConditionDirichlet(initial_condition_convergence_test)`](@ref) in non-periodic domains).
 """
-function initial_condition_convergance_test(x, t,
+function initial_condition_convergence_test(x, t,
                                             equations::CompressibleEulerEquationsQuasi1D)
     c = 2
     A = 0.1
@@ -102,11 +102,11 @@ Source terms used for convergence tests in combination with
 (and [`BoundaryConditionDirichlet(initial_condition_convergence_test)`](@ref) in non-periodic domains).
 
 This manufactured solution source term is specifically designed for the mozzle width 'a(x) = 1.5 - 0.5 * cos(x[1] * pi)'
-as defined in [`initial_condition_convergance_test`](@ref).
+as defined in [`initial_condition_convergence_test`](@ref).
 """
 @inline function source_terms_convergence_test(u, x, t,
                                                equations::CompressibleEulerEquationsQuasi1D)
-    # Same settings as in `initial_condition_convergance_test`. 
+    # Same settings as in `initial_condition_convergence_test`. 
     #Derivatives calculated with ForwardDiff.jl
     c = 2
     A = 0.1

From f51e16153853c54a1b0e002d8f6e7b1783204ac4 Mon Sep 17 00:00:00 2001
From: KrisshChawla <krisshchawla@gmail.com>
Date: Fri, 24 Nov 2023 17:32:05 -0600
Subject: [PATCH 08/39] formatting

---
 examples/tree_1d_dgsem/elixir_euler_quasi_1d_ec.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/examples/tree_1d_dgsem/elixir_euler_quasi_1d_ec.jl b/examples/tree_1d_dgsem/elixir_euler_quasi_1d_ec.jl
index 81ca5c13023..04422dd6d29 100644
--- a/examples/tree_1d_dgsem/elixir_euler_quasi_1d_ec.jl
+++ b/examples/tree_1d_dgsem/elixir_euler_quasi_1d_ec.jl
@@ -12,7 +12,7 @@ function initial_condition_ec(x, t, equations::CompressibleEulerEquationsQuasi1D
     rho = 2.0 + 0.1 * x[1]
     p = 3.0
     a = 2.0 + x[1]
-    
+
     return prim2cons(SVector(rho, v1, p, a), equations)
 end
 
@@ -65,4 +65,4 @@ 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
\ No newline at end of file
+summary_callback() # print the timer summary

From 862f8dc11402e44d17e71f16e84e221875a5bbb2 Mon Sep 17 00:00:00 2001
From: Krissh Chawla <127906314+KrisshChawla@users.noreply.github.com>
Date: Fri, 24 Nov 2023 19:13:13 -0600
Subject: [PATCH 09/39] Update
 examples/tree_1d_dgsem/elixir_euler_quasi_1d_discontinuous.jl

Co-authored-by: Jesse Chan <1156048+jlchan@users.noreply.github.com>
---
 examples/tree_1d_dgsem/elixir_euler_quasi_1d_discontinuous.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/examples/tree_1d_dgsem/elixir_euler_quasi_1d_discontinuous.jl b/examples/tree_1d_dgsem/elixir_euler_quasi_1d_discontinuous.jl
index 8f0196b3d2e..949348afc13 100644
--- a/examples/tree_1d_dgsem/elixir_euler_quasi_1d_discontinuous.jl
+++ b/examples/tree_1d_dgsem/elixir_euler_quasi_1d_discontinuous.jl
@@ -10,7 +10,7 @@ equations = CompressibleEulerEquationsQuasi1D(1.4)
 """
     initial_condition_discontinuity(x, t, equations::CompressibleEulerEquations1D)
 
-An entropy conservation verification initial condition taken from
+A discontinuous initial condition taken from
 - Jesse Chan, Khemraj Shukla, Xinhui Wu, Ruofeng Liu, Prani Nalluri (2023)
     High order entropy stable schemes for the quasi-one-dimensional
     shallow water and compressible Euler equations

From 110d268e8b88c5d0155b094386fca301d79a06c5 Mon Sep 17 00:00:00 2001
From: Krissh Chawla <127906314+KrisshChawla@users.noreply.github.com>
Date: Fri, 24 Nov 2023 19:13:30 -0600
Subject: [PATCH 10/39] Update src/equations/compressible_euler_quasi_1d.jl

Co-authored-by: Jesse Chan <1156048+jlchan@users.noreply.github.com>
---
 src/equations/compressible_euler_quasi_1d.jl | 7 +++++--
 1 file changed, 5 insertions(+), 2 deletions(-)

diff --git a/src/equations/compressible_euler_quasi_1d.jl b/src/equations/compressible_euler_quasi_1d.jl
index 34e6cafe2dc..1b0be2bcd1f 100644
--- a/src/equations/compressible_euler_quasi_1d.jl
+++ b/src/equations/compressible_euler_quasi_1d.jl
@@ -232,11 +232,14 @@ Further details are available in the paper:
     _, _, p_ll, a_ll = cons2prim(u_ll, equations)
     _, _, p_rr, _ = cons2prim(u_rr, equations)
 
-    p_avg = 0.5 * (p_ll + p_rr)
+    # For flux differencing using non-conservative terms, we return the 
+    # non-conservative flux scaled by 2. This cancels with a factor of 0.5 
+    # in the arithmetic average of {p}.
+    p_avg = p_ll + p_rr 
 
     z = zero(eltype(u_ll))
 
-    return SVector(z, 2 * a_ll * p_avg, z, z)
+    return SVector(z, a_ll * p_avg, z, z)
 end
 
 """

From 59fdbbf0abc4bbe350950f21d2f0f3632634177a Mon Sep 17 00:00:00 2001
From: Krissh Chawla <127906314+KrisshChawla@users.noreply.github.com>
Date: Fri, 24 Nov 2023 19:13:56 -0600
Subject: [PATCH 11/39] Update src/equations/compressible_euler_quasi_1d.jl

Co-authored-by: Jesse Chan <1156048+jlchan@users.noreply.github.com>
---
 src/equations/compressible_euler_quasi_1d.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/equations/compressible_euler_quasi_1d.jl b/src/equations/compressible_euler_quasi_1d.jl
index 1b0be2bcd1f..db10cdaed1a 100644
--- a/src/equations/compressible_euler_quasi_1d.jl
+++ b/src/equations/compressible_euler_quasi_1d.jl
@@ -246,7 +246,7 @@ end
 @inline function flux_chan_etal(u_ll, u_rr, orientation::Integer,
                                            equations::CompressibleEulerEquationsQuasi1D)
 
-Total energy conservative (mathematical entropy for quasi 1D compressible Euler equations) split form.
+Conservative (symmetric) part of the entropy conservative flux for quasi 1D compressible Euler equations split form.
 
 Further details are available in the paper:
 - Jesse Chan, Khemraj Shukla, Xinhui Wu, Ruofeng Liu, Prani Nalluri (2023) 

From b13eb72121759863a56c823dba8fdb72cc1ff951 Mon Sep 17 00:00:00 2001
From: Krissh Chawla <127906314+KrisshChawla@users.noreply.github.com>
Date: Fri, 24 Nov 2023 19:14:19 -0600
Subject: [PATCH 12/39] Update src/equations/compressible_euler_quasi_1d.jl

Co-authored-by: Jesse Chan <1156048+jlchan@users.noreply.github.com>
---
 src/equations/compressible_euler_quasi_1d.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/equations/compressible_euler_quasi_1d.jl b/src/equations/compressible_euler_quasi_1d.jl
index db10cdaed1a..849706e73f0 100644
--- a/src/equations/compressible_euler_quasi_1d.jl
+++ b/src/equations/compressible_euler_quasi_1d.jl
@@ -247,7 +247,7 @@ end
                                            equations::CompressibleEulerEquationsQuasi1D)
 
 Conservative (symmetric) part of the entropy conservative flux for quasi 1D compressible Euler equations split form.
-
+This flux is a generalization of [`flux_ranocha`](@ref) for [`CompressibleEulerEquations1D`](@ref).
 Further details are available in the paper:
 - Jesse Chan, Khemraj Shukla, Xinhui Wu, Ruofeng Liu, Prani Nalluri (2023) 
   High order entropy stable schemes for the quasi-one-dimensional

From bd2a155a54048b6f232e805b2ac18bca0a0da920 Mon Sep 17 00:00:00 2001
From: KrisshChawla <krisshchawla@gmail.com>
Date: Fri, 24 Nov 2023 20:20:41 -0600
Subject: [PATCH 13/39] updated test_tree_1d_euler.jl and formatting

---
 .../tree_1d_dgsem/elixir_euler_quasi_1d_ec.jl | 15 ++++++++----
 test/test_tree_1d_euler.jl                    | 24 +++++++++++++++++++
 2 files changed, 34 insertions(+), 5 deletions(-)

diff --git a/examples/tree_1d_dgsem/elixir_euler_quasi_1d_ec.jl b/examples/tree_1d_dgsem/elixir_euler_quasi_1d_ec.jl
index 04422dd6d29..74efe9ce57c 100644
--- a/examples/tree_1d_dgsem/elixir_euler_quasi_1d_ec.jl
+++ b/examples/tree_1d_dgsem/elixir_euler_quasi_1d_ec.jl
@@ -2,11 +2,18 @@ using OrdinaryDiffEq
 using Trixi
 
 ###############################################################################
-# Semidiscretization of the quasi 1d compressible Euler equations
+# Semidiscretization of the quasi 1d compressible Euler equations with a discontinuous nozzle width function.
 # See Chan et al.  https://doi.org/10.48550/arXiv.2307.12089 for details
 
 equations = CompressibleEulerEquationsQuasi1D(1.4)
 
+# Setup a truly discontinuous density function and nozzle width for
+# this academic testcase of entropy conservation. The errors from the analysis
+# callback are not important but the entropy error for this test case
+# `∑∂S/∂U ⋅ Uₜ` should be around machine roundoff.
+# Works as intended for TreeMesh1D with `initial_refinement_level=6`. If the mesh
+# refinement level is changed the initial condition below may need changed as well to
+# ensure that the discontinuities lie on an element interface.
 function initial_condition_ec(x, t, equations::CompressibleEulerEquationsQuasi1D)
     v1 = 0.1
     rho = 2.0 + 0.1 * x[1]
@@ -20,8 +27,6 @@ initial_condition = initial_condition_ec
 
 volume_flux = (flux_chan_etal, flux_nonconservative_chan_etal)
 surface_flux = volume_flux
-
-basis = LobattoLegendreBasis(3)
 solver = DGSEM(polydeg = 4, surface_flux = surface_flux,
                volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
 
@@ -36,7 +41,7 @@ semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
 ###############################################################################
 # ODE solvers, callbacks etc.
 
-tspan = (0.0, 2.0)
+tspan = (0.0, 0.4)
 ode = semidiscretize(semi, tspan)
 
 summary_callback = SummaryCallback()
@@ -52,7 +57,7 @@ save_solution = SaveSolutionCallback(interval = 100,
                                      save_final_solution = true,
                                      solution_variables = cons2prim)
 
-stepsize_callback = StepsizeCallback(cfl = 0.5)
+stepsize_callback = StepsizeCallback(cfl = 0.8)
 
 callbacks = CallbackSet(summary_callback,
                         analysis_callback, alive_callback,
diff --git a/test/test_tree_1d_euler.jl b/test/test_tree_1d_euler.jl
index 5dc81da0261..629bea05e5a 100644
--- a/test/test_tree_1d_euler.jl
+++ b/test/test_tree_1d_euler.jl
@@ -442,6 +442,30 @@ end
         @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
     end
 end
+
+@trixi_testset "elixir_euler_quasi_1d_ec.jl" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_quasi_1d_ec.jl"),
+                        l2=[
+                            0.08889113985713934,
+                            0.16199235348889698,
+                            0.4031652436505426,
+                            2.9602775074723667e-16,
+                        ],
+                        linf=[
+                            0.28891355898287685,
+                            0.3752709888964392,
+                            0.8447710240240944,
+                            8.881784197001252e-16,
+                        ])
+    # 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

From b49a388854d2658aa6d78d156b19dcd4f0492497 Mon Sep 17 00:00:00 2001
From: Krissh Chawla <127906314+KrisshChawla@users.noreply.github.com>
Date: Sat, 25 Nov 2023 10:11:49 -0600
Subject: [PATCH 14/39] Update src/equations/compressible_euler_quasi_1d.jl

Co-authored-by: Daniel Doehring <doehringd2@gmail.com>
---
 src/equations/compressible_euler_quasi_1d.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/equations/compressible_euler_quasi_1d.jl b/src/equations/compressible_euler_quasi_1d.jl
index 849706e73f0..d1c4c6ff446 100644
--- a/src/equations/compressible_euler_quasi_1d.jl
+++ b/src/equations/compressible_euler_quasi_1d.jl
@@ -198,7 +198,7 @@ are available in the paper:
 end
 
 # Calculate 1D flux for a single point
-@inline function Trixi.flux(u, orientation::Integer,
+@inline function flux(u, orientation::Integer,
                             equations::CompressibleEulerEquationsQuasi1D)
     a_rho, a_rho_v1, a_e, a = u
     rho, v1, p, a = cons2prim(u, equations)

From 4522f4cc6de2bf10834bb6f2efe43a232cf17cd8 Mon Sep 17 00:00:00 2001
From: Krissh Chawla <127906314+KrisshChawla@users.noreply.github.com>
Date: Sat, 25 Nov 2023 10:12:06 -0600
Subject: [PATCH 15/39] Update src/equations/compressible_euler_quasi_1d.jl

Co-authored-by: Daniel Doehring <doehringd2@gmail.com>
---
 src/equations/compressible_euler_quasi_1d.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/equations/compressible_euler_quasi_1d.jl b/src/equations/compressible_euler_quasi_1d.jl
index d1c4c6ff446..f4721902a4e 100644
--- a/src/equations/compressible_euler_quasi_1d.jl
+++ b/src/equations/compressible_euler_quasi_1d.jl
@@ -306,7 +306,7 @@ end
     λ_max = max(v_mag_ll, v_mag_rr) + max(c_ll, c_rr)
 end
 
-@inline function Trixi.max_abs_speeds(u, equations::CompressibleEulerEquationsQuasi1D)
+@inline function max_abs_speeds(u, equations::CompressibleEulerEquationsQuasi1D)
     a_rho, a_rho_v1, a_e, a = u
     rho = a_rho / a
     v1 = a_rho_v1 / a_rho

From 20c0a101cebd3d4bf903e34a62fa5d3d8dd93c0a Mon Sep 17 00:00:00 2001
From: Krissh Chawla <127906314+KrisshChawla@users.noreply.github.com>
Date: Sat, 25 Nov 2023 10:12:20 -0600
Subject: [PATCH 16/39] Update src/equations/compressible_euler_quasi_1d.jl

Co-authored-by: Daniel Doehring <doehringd2@gmail.com>
---
 src/equations/compressible_euler_quasi_1d.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/equations/compressible_euler_quasi_1d.jl b/src/equations/compressible_euler_quasi_1d.jl
index f4721902a4e..9d7c3c80b30 100644
--- a/src/equations/compressible_euler_quasi_1d.jl
+++ b/src/equations/compressible_euler_quasi_1d.jl
@@ -345,7 +345,7 @@ end
     return SVector(a * q[1], a * q[2], a * q[3], a)
 end
 
-@inline function Trixi.density(u, equations::CompressibleEulerEquationsQuasi1D)
+@inline function density(u, equations::CompressibleEulerEquationsQuasi1D)
     rho = u[1] / u[4]
     return rho
 end

From ecb5345801e269bb04cf2c7ca29998649db274ea Mon Sep 17 00:00:00 2001
From: Krissh Chawla <127906314+KrisshChawla@users.noreply.github.com>
Date: Sat, 25 Nov 2023 10:12:31 -0600
Subject: [PATCH 17/39] Update src/equations/compressible_euler_quasi_1d.jl

Co-authored-by: Daniel Doehring <doehringd2@gmail.com>
---
 src/equations/compressible_euler_quasi_1d.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/equations/compressible_euler_quasi_1d.jl b/src/equations/compressible_euler_quasi_1d.jl
index 9d7c3c80b30..7d43acf446d 100644
--- a/src/equations/compressible_euler_quasi_1d.jl
+++ b/src/equations/compressible_euler_quasi_1d.jl
@@ -350,7 +350,7 @@ end
     return rho
 end
 
-@inline function Trixi.pressure(u, equations::CompressibleEulerEquationsQuasi1D)
+@inline function pressure(u, equations::CompressibleEulerEquationsQuasi1D)
     a_rho, a_rho_v1, a_e, a = u
     rho = a_rho / a
     v1 = a_rho_v1 / a_rho

From 6e1156518542e350a8b6e42d9195958313fa2ffa Mon Sep 17 00:00:00 2001
From: Krissh Chawla <127906314+KrisshChawla@users.noreply.github.com>
Date: Sat, 25 Nov 2023 10:12:41 -0600
Subject: [PATCH 18/39] Update src/equations/compressible_euler_quasi_1d.jl

Co-authored-by: Daniel Doehring <doehringd2@gmail.com>
---
 src/equations/compressible_euler_quasi_1d.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/equations/compressible_euler_quasi_1d.jl b/src/equations/compressible_euler_quasi_1d.jl
index 7d43acf446d..a4e2b2f1199 100644
--- a/src/equations/compressible_euler_quasi_1d.jl
+++ b/src/equations/compressible_euler_quasi_1d.jl
@@ -359,7 +359,7 @@ end
     return p
 end
 
-@inline function Trixi.density_pressure(u, equations::CompressibleEulerEquationsQuasi1D)
+@inline function density_pressure(u, equations::CompressibleEulerEquationsQuasi1D)
     a_rho, a_rho_v1, a_e, a = u
     rho = a_rho / a
     v1 = a_rho_v1 / a_rho

From 5f9c15ca565fed6b459fad79875fb4746623e6c5 Mon Sep 17 00:00:00 2001
From: Krissh Chawla <127906314+KrisshChawla@users.noreply.github.com>
Date: Sat, 25 Nov 2023 10:12:52 -0600
Subject: [PATCH 19/39] Update src/equations/compressible_euler_quasi_1d.jl

Co-authored-by: Daniel Doehring <doehringd2@gmail.com>
---
 src/equations/compressible_euler_quasi_1d.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/equations/compressible_euler_quasi_1d.jl b/src/equations/compressible_euler_quasi_1d.jl
index a4e2b2f1199..6312414ba3f 100644
--- a/src/equations/compressible_euler_quasi_1d.jl
+++ b/src/equations/compressible_euler_quasi_1d.jl
@@ -326,7 +326,7 @@ end
 end
 
 # Convert conservative variables to entropy
-@inline function Trixi.cons2entropy(u, equations::CompressibleEulerEquationsQuasi1D)
+@inline function cons2entropy(u, equations::CompressibleEulerEquationsQuasi1D)
     a_rho, a_rho_v1, a_e, a = u
     q = cons2entropy(u, CompressibleEulerEquations1D(equations.gamma))
 

From dd388409e9c8e177316cfb25a37a9d1e71006074 Mon Sep 17 00:00:00 2001
From: KrisshChawla <krisshchawla@gmail.com>
Date: Sat, 25 Nov 2023 11:27:43 -0600
Subject: [PATCH 20/39] adding consistency check for flux

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

diff --git a/test/test_unit.jl b/test/test_unit.jl
index d2e744da62f..8276ee87b53 100644
--- a/test/test_unit.jl
+++ b/test/test_unit.jl
@@ -654,6 +654,17 @@ end
     end
 end
 
+@timed_testset "Consistency check for flux_chan_etal: CEEQ" begin
+
+    # Set up equations and dummy conservative variables state
+    equations = CompressibleEulerEquationsQuasi1D(1.4)
+    u = SVector(1.1, 2.34, 5.5, 2.73)
+
+    orientations = [1]
+    for orientation in orientations
+        @test flux_chan_etal(u, u, orientation, equations) ≈ flux(u, orientation, equations)
+    end
+
 @timed_testset "Consistency check for HLL flux (naive): LEE" begin
     flux_hll = FluxHLL(min_max_speed_naive)
 

From 33bc713bf8071fa5088fc73e5f8c86971aeb2f25 Mon Sep 17 00:00:00 2001
From: KrisshChawla <krisshchawla@gmail.com>
Date: Sat, 25 Nov 2023 11:42:20 -0600
Subject: [PATCH 21/39] Update compressible_euler_quasi_1d.jl

---
 src/equations/compressible_euler_quasi_1d.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/equations/compressible_euler_quasi_1d.jl b/src/equations/compressible_euler_quasi_1d.jl
index 6312414ba3f..6656b4cbece 100644
--- a/src/equations/compressible_euler_quasi_1d.jl
+++ b/src/equations/compressible_euler_quasi_1d.jl
@@ -8,7 +8,7 @@
 @doc raw"""
     CompressibleEulerEquationsQuasi1D(gamma)
 
-The quasi-1d compressible Euler equations
+The quasi-1d compressible Euler equations (see Chan et al.  https://doi.org/10.48550/arXiv.2307.12089 for details)
 ```math
 \frac{\partial}{\partial t}
 \begin{pmatrix}

From 240690ebc4b34f70eebf4862a4894fe4eb436471 Mon Sep 17 00:00:00 2001
From: KrisshChawla <krisshchawla@gmail.com>
Date: Sat, 25 Nov 2023 11:53:01 -0600
Subject: [PATCH 22/39] formatting

---
 src/equations/compressible_euler_quasi_1d.jl | 4 ++--
 test/test_unit.jl                            | 4 +++-
 2 files changed, 5 insertions(+), 3 deletions(-)

diff --git a/src/equations/compressible_euler_quasi_1d.jl b/src/equations/compressible_euler_quasi_1d.jl
index 6656b4cbece..7e4b0228cff 100644
--- a/src/equations/compressible_euler_quasi_1d.jl
+++ b/src/equations/compressible_euler_quasi_1d.jl
@@ -199,7 +199,7 @@ end
 
 # Calculate 1D flux for a single point
 @inline function flux(u, orientation::Integer,
-                            equations::CompressibleEulerEquationsQuasi1D)
+                      equations::CompressibleEulerEquationsQuasi1D)
     a_rho, a_rho_v1, a_e, a = u
     rho, v1, p, a = cons2prim(u, equations)
     e = a_e / a
@@ -235,7 +235,7 @@ Further details are available in the paper:
     # For flux differencing using non-conservative terms, we return the 
     # non-conservative flux scaled by 2. This cancels with a factor of 0.5 
     # in the arithmetic average of {p}.
-    p_avg = p_ll + p_rr 
+    p_avg = p_ll + p_rr
 
     z = zero(eltype(u_ll))
 
diff --git a/test/test_unit.jl b/test/test_unit.jl
index 8276ee87b53..b3ed29d38e3 100644
--- a/test/test_unit.jl
+++ b/test/test_unit.jl
@@ -662,8 +662,10 @@ end
 
     orientations = [1]
     for orientation in orientations
-        @test flux_chan_etal(u, u, orientation, equations) ≈ flux(u, orientation, equations)
+        @test flux_chan_etal(u, u, orientation, equations) ≈
+              flux(u, orientation, equations)
     end
+end
 
 @timed_testset "Consistency check for HLL flux (naive): LEE" begin
     flux_hll = FluxHLL(min_max_speed_naive)

From 582078a991b55955c32d64977c4b85936538a847 Mon Sep 17 00:00:00 2001
From: Krissh Chawla <127906314+KrisshChawla@users.noreply.github.com>
Date: Sat, 25 Nov 2023 13:13:33 -0600
Subject: [PATCH 23/39] Update src/equations/compressible_euler_quasi_1d.jl

Co-authored-by: Jesse Chan <1156048+jlchan@users.noreply.github.com>
---
 src/equations/compressible_euler_quasi_1d.jl | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/equations/compressible_euler_quasi_1d.jl b/src/equations/compressible_euler_quasi_1d.jl
index 7e4b0228cff..e24d82cb5e2 100644
--- a/src/equations/compressible_euler_quasi_1d.jl
+++ b/src/equations/compressible_euler_quasi_1d.jl
@@ -317,12 +317,12 @@ end
     return (abs(v1) + c,)
 end
 
-# Convert conservative variables to primitive
+# Convert conservative variables to primitive. We use the convention that the primitive
+# variables for the quasi-1D equations are `(rho, v1, p)` (i.e., the same as the primitive
+# variables for `CompressibleEulerEquations1D`)
 @inline function cons2prim(u, equations::CompressibleEulerEquationsQuasi1D)
     a_rho, a_rho_v1, a_e, a = u
-    q = cons2prim(u, CompressibleEulerEquations1D(equations.gamma))
-
-    return SVector(q[1] / a, q[2], q[3] / a, a)
+    return cons2prim(SVector(a_rho, a_rho_v1, a_e) / a, CompressibleEulerEquations1D(equations.gamma))
 end
 
 # Convert conservative variables to entropy

From 819020831fa9e448950e06294a161a34167b1558 Mon Sep 17 00:00:00 2001
From: Krissh Chawla <127906314+KrisshChawla@users.noreply.github.com>
Date: Sat, 25 Nov 2023 13:14:19 -0600
Subject: [PATCH 24/39] Update src/equations/compressible_euler_quasi_1d.jl

Co-authored-by: Jesse Chan <1156048+jlchan@users.noreply.github.com>
---
 src/equations/compressible_euler_quasi_1d.jl | 24 +++++++++++++-------
 1 file changed, 16 insertions(+), 8 deletions(-)

diff --git a/src/equations/compressible_euler_quasi_1d.jl b/src/equations/compressible_euler_quasi_1d.jl
index e24d82cb5e2..a3685af11ce 100644
--- a/src/equations/compressible_euler_quasi_1d.jl
+++ b/src/equations/compressible_euler_quasi_1d.jl
@@ -325,16 +325,24 @@ end
     return cons2prim(SVector(a_rho, a_rho_v1, a_e) / a, CompressibleEulerEquations1D(equations.gamma))
 end
 
-# Convert conservative variables to entropy
-@inline function cons2entropy(u, equations::CompressibleEulerEquationsQuasi1D)
+# The entropy for the quasi-1D compressible Euler equations is the entropy for the
+# 1D compressible Euler equations scaled by the channel width `a`.
+@inline function entropy(u, equations::CompressibleEulerEquationsQuasi1D)
     a_rho, a_rho_v1, a_e, a = u
-    q = cons2entropy(u, CompressibleEulerEquations1D(equations.gamma))
-
-    w1 = q[1] - log(a)
-    w2 = q[2]
-    w3 = q[3]
+    return a * entropy(SVector(a_rho, a_rho_v1, a_e) / a, CompressibleEulerEquations1D(equations.gamma))
+end
 
-    return SVector(w1, w2, w3, a)
+# Convert conservative variables to entropy. The entropy variables for the 
+# quasi-1D compressible Euler equations are identical to the entropy variables
+# for the standard Euler equations for an appropriate definition of `entropy`.
+@inline function cons2entropy(u, equations::CompressibleEulerEquationsQuasi1D)
+    a_rho, a_rho_v1, a_e, a = u
+    w = cons2entropy(SVector(a_rho, a_rho_v1, a_e) / a, CompressibleEulerEquations1D(equations.gamma))
+    
+    # we follow the convention for other spatially-varying equations such as
+    # `ShallowWaterEquations1D` and return the spatially varying coefficient 
+    # `a` as the final entropy variable.
+    return SVector(w[1], w[2], w[3], a)
 end
 
 # Convert primitive to conservative variables

From f6c822ad971eb6495fddfccb30033bbc0d705345 Mon Sep 17 00:00:00 2001
From: Krissh Chawla <127906314+KrisshChawla@users.noreply.github.com>
Date: Sat, 25 Nov 2023 13:14:50 -0600
Subject: [PATCH 25/39] Update src/equations/compressible_euler_quasi_1d.jl

Co-authored-by: Jesse Chan <1156048+jlchan@users.noreply.github.com>
---
 src/equations/compressible_euler_quasi_1d.jl | 6 +-----
 1 file changed, 1 insertion(+), 5 deletions(-)

diff --git a/src/equations/compressible_euler_quasi_1d.jl b/src/equations/compressible_euler_quasi_1d.jl
index a3685af11ce..b99226f4abe 100644
--- a/src/equations/compressible_euler_quasi_1d.jl
+++ b/src/equations/compressible_euler_quasi_1d.jl
@@ -360,11 +360,7 @@ end
 
 @inline function pressure(u, equations::CompressibleEulerEquationsQuasi1D)
     a_rho, a_rho_v1, a_e, a = u
-    rho = a_rho / a
-    v1 = a_rho_v1 / a_rho
-    e = a_e / a
-    p = (equations.gamma - 1) * (e - 0.5 * rho * v1^2)
-    return p
+    return pressure(SVector(a_rho, a_rho_v1, a_e) / a, CompressibleEulerEquations1D(equations.gamma))
 end
 
 @inline function density_pressure(u, equations::CompressibleEulerEquationsQuasi1D)

From 5827a8730bb6fc0c1076160a659056689151afef Mon Sep 17 00:00:00 2001
From: Krissh Chawla <127906314+KrisshChawla@users.noreply.github.com>
Date: Sat, 25 Nov 2023 13:15:05 -0600
Subject: [PATCH 26/39] Update src/equations/compressible_euler_quasi_1d.jl

Co-authored-by: Jesse Chan <1156048+jlchan@users.noreply.github.com>
---
 src/equations/compressible_euler_quasi_1d.jl | 7 +------
 1 file changed, 1 insertion(+), 6 deletions(-)

diff --git a/src/equations/compressible_euler_quasi_1d.jl b/src/equations/compressible_euler_quasi_1d.jl
index b99226f4abe..f75c10538bf 100644
--- a/src/equations/compressible_euler_quasi_1d.jl
+++ b/src/equations/compressible_euler_quasi_1d.jl
@@ -365,11 +365,6 @@ end
 
 @inline function density_pressure(u, equations::CompressibleEulerEquationsQuasi1D)
     a_rho, a_rho_v1, a_e, a = u
-    rho = a_rho / a
-    v1 = a_rho_v1 / a_rho
-    rho_v1 = a_rho_v1 / a
-    e = a_e / a
-    rho_times_p = (equations.gamma - 1) * (e - 0.5 * rho * v1^2) * rho
-    return rho_times_p
+    return density_pressure(SVector(a_rho, a_rho_v1, a_e) / a, CompressibleEulerEquations(equations.gamma))
 end
 end # @muladd

From dc7ccf05dda3c52a3f7eda4e0ac376b830805651 Mon Sep 17 00:00:00 2001
From: KrisshChawla <krisshchawla@gmail.com>
Date: Sat, 25 Nov 2023 13:53:25 -0600
Subject: [PATCH 27/39] Update compressible_euler_quasi_1d and
 test_tree_1d_euler

---
 src/equations/compressible_euler_quasi_1d.jl | 10 +++--
 test/test_tree_1d_euler.jl                   | 40 ++++++++++----------
 2 files changed, 27 insertions(+), 23 deletions(-)

diff --git a/src/equations/compressible_euler_quasi_1d.jl b/src/equations/compressible_euler_quasi_1d.jl
index f75c10538bf..c17076fb641 100644
--- a/src/equations/compressible_euler_quasi_1d.jl
+++ b/src/equations/compressible_euler_quasi_1d.jl
@@ -322,14 +322,18 @@ end
 # variables for `CompressibleEulerEquations1D`)
 @inline function cons2prim(u, equations::CompressibleEulerEquationsQuasi1D)
     a_rho, a_rho_v1, a_e, a = u
-    return cons2prim(SVector(a_rho, a_rho_v1, a_e) / a, CompressibleEulerEquations1D(equations.gamma))
+    q = cons2prim(SVector(a_rho, a_rho_v1, a_e) / a, CompressibleEulerEquations1D(equations.gamma))
+
+    return SVector(q[1], q[2], q[3], a)
 end
 
 # The entropy for the quasi-1D compressible Euler equations is the entropy for the
 # 1D compressible Euler equations scaled by the channel width `a`.
 @inline function entropy(u, equations::CompressibleEulerEquationsQuasi1D)
     a_rho, a_rho_v1, a_e, a = u
-    return a * entropy(SVector(a_rho, a_rho_v1, a_e) / a, CompressibleEulerEquations1D(equations.gamma))
+    q = a * entropy(SVector(a_rho, a_rho_v1, a_e) / a, CompressibleEulerEquations1D(equations.gamma))
+
+    return SVector(q[1], q[2], q[3], a)
 end
 
 # Convert conservative variables to entropy. The entropy variables for the 
@@ -365,6 +369,6 @@ end
 
 @inline function density_pressure(u, equations::CompressibleEulerEquationsQuasi1D)
     a_rho, a_rho_v1, a_e, a = u
-    return density_pressure(SVector(a_rho, a_rho_v1, a_e) / a, CompressibleEulerEquations(equations.gamma))
+    return density_pressure(SVector(a_rho, a_rho_v1, a_e) / a, CompressibleEulerEquations1D(equations.gamma))
 end
 end # @muladd
diff --git a/test/test_tree_1d_euler.jl b/test/test_tree_1d_euler.jl
index 629bea05e5a..39a1f6e30ba 100644
--- a/test/test_tree_1d_euler.jl
+++ b/test/test_tree_1d_euler.jl
@@ -397,15 +397,15 @@ end
 @trixi_testset "elixir_euler_quasi_1d_source_terms.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_quasi_1d_source_terms.jl"),
                         l2=[
-                            3.876288364589425e-7,
-                            2.224704317890825e-7,
-                            2.964004223514882e-7,
+                            3.876288369618363e-7,
+                            2.2247043122302947e-7,
+                            2.964004224572679e-7,
                             5.2716983399807875e-8,
                         ],
                         linf=[
-                            2.392511842419509e-6,
-                            1.360369332736866e-6,
-                            1.8218888442333991e-6,
+                            2.3925118561862746e-6,
+                            1.3603693522767912e-6,
+                            1.821888865105592e-6,
                             1.1166012159335992e-7,
                         ])
     # Ensure that we do not have excessive memory allocations
@@ -422,16 +422,16 @@ end
     @test_trixi_include(joinpath(EXAMPLES_DIR,
                                  "elixir_euler_quasi_1d_discontinuous.jl"),
                         l2=[
-                            0.04551042115634054,
-                            0.03675058478890224,
-                            0.24689859591319177,
-                            0.0368449418082899,
+                            0.045510421156346015,
+                            0.036750584788912195,
+                            0.2468985959132176,
+                            0.03684494180829024,
                         ],
                         linf=[
-                            0.3313374853025661,
-                            0.11621933362147452,
-                            1.8274030135685742,
-                            0.28045939999015945,
+                            0.3313374853025697,
+                            0.11621933362158643,
+                            1.827403013568638,
+                            0.28045939999015723,
                         ])
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
@@ -446,15 +446,15 @@ end
 @trixi_testset "elixir_euler_quasi_1d_ec.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_quasi_1d_ec.jl"),
                         l2=[
-                            0.08889113985713934,
-                            0.16199235348889698,
-                            0.4031652436505426,
+                            0.08889113985713998,
+                            0.16199235348889673,
+                            0.40316524365054346,
                             2.9602775074723667e-16,
                         ],
                         linf=[
-                            0.28891355898287685,
-                            0.3752709888964392,
-                            0.8447710240240944,
+                            0.28891355898284043,
+                            0.3752709888964313,
+                            0.84477102402413,
                             8.881784197001252e-16,
                         ])
     # Ensure that we do not have excessive memory allocations

From d1ccb20b04f9fce4b5e978a0c2a8e59791403f34 Mon Sep 17 00:00:00 2001
From: KrisshChawla <krisshchawla@gmail.com>
Date: Sat, 25 Nov 2023 14:02:13 -0600
Subject: [PATCH 28/39] formatting

---
 src/equations/compressible_euler_quasi_1d.jl | 17 +++++++++++------
 1 file changed, 11 insertions(+), 6 deletions(-)

diff --git a/src/equations/compressible_euler_quasi_1d.jl b/src/equations/compressible_euler_quasi_1d.jl
index c17076fb641..029d67fb546 100644
--- a/src/equations/compressible_euler_quasi_1d.jl
+++ b/src/equations/compressible_euler_quasi_1d.jl
@@ -322,7 +322,8 @@ end
 # variables for `CompressibleEulerEquations1D`)
 @inline function cons2prim(u, equations::CompressibleEulerEquationsQuasi1D)
     a_rho, a_rho_v1, a_e, a = u
-    q = cons2prim(SVector(a_rho, a_rho_v1, a_e) / a, CompressibleEulerEquations1D(equations.gamma))
+    q = cons2prim(SVector(a_rho, a_rho_v1, a_e) / a,
+                  CompressibleEulerEquations1D(equations.gamma))
 
     return SVector(q[1], q[2], q[3], a)
 end
@@ -331,7 +332,8 @@ end
 # 1D compressible Euler equations scaled by the channel width `a`.
 @inline function entropy(u, equations::CompressibleEulerEquationsQuasi1D)
     a_rho, a_rho_v1, a_e, a = u
-    q = a * entropy(SVector(a_rho, a_rho_v1, a_e) / a, CompressibleEulerEquations1D(equations.gamma))
+    q = a * entropy(SVector(a_rho, a_rho_v1, a_e) / a,
+                CompressibleEulerEquations1D(equations.gamma))
 
     return SVector(q[1], q[2], q[3], a)
 end
@@ -341,8 +343,9 @@ end
 # for the standard Euler equations for an appropriate definition of `entropy`.
 @inline function cons2entropy(u, equations::CompressibleEulerEquationsQuasi1D)
     a_rho, a_rho_v1, a_e, a = u
-    w = cons2entropy(SVector(a_rho, a_rho_v1, a_e) / a, CompressibleEulerEquations1D(equations.gamma))
-    
+    w = cons2entropy(SVector(a_rho, a_rho_v1, a_e) / a,
+                     CompressibleEulerEquations1D(equations.gamma))
+
     # we follow the convention for other spatially-varying equations such as
     # `ShallowWaterEquations1D` and return the spatially varying coefficient 
     # `a` as the final entropy variable.
@@ -364,11 +367,13 @@ end
 
 @inline function pressure(u, equations::CompressibleEulerEquationsQuasi1D)
     a_rho, a_rho_v1, a_e, a = u
-    return pressure(SVector(a_rho, a_rho_v1, a_e) / a, CompressibleEulerEquations1D(equations.gamma))
+    return pressure(SVector(a_rho, a_rho_v1, a_e) / a,
+                    CompressibleEulerEquations1D(equations.gamma))
 end
 
 @inline function density_pressure(u, equations::CompressibleEulerEquationsQuasi1D)
     a_rho, a_rho_v1, a_e, a = u
-    return density_pressure(SVector(a_rho, a_rho_v1, a_e) / a, CompressibleEulerEquations1D(equations.gamma))
+    return density_pressure(SVector(a_rho, a_rho_v1, a_e) / a,
+                            CompressibleEulerEquations1D(equations.gamma))
 end
 end # @muladd

From 5854bd7fa86031957a48f0fe420f6a029bf691db Mon Sep 17 00:00:00 2001
From: Krissh Chawla <127906314+KrisshChawla@users.noreply.github.com>
Date: Sun, 26 Nov 2023 10:43:06 -0600
Subject: [PATCH 29/39] Update
 examples/tree_1d_dgsem/elixir_euler_quasi_1d_discontinuous.jl

Co-authored-by: Jesse Chan <1156048+jlchan@users.noreply.github.com>
---
 examples/tree_1d_dgsem/elixir_euler_quasi_1d_discontinuous.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/examples/tree_1d_dgsem/elixir_euler_quasi_1d_discontinuous.jl b/examples/tree_1d_dgsem/elixir_euler_quasi_1d_discontinuous.jl
index 949348afc13..cc4535be028 100644
--- a/examples/tree_1d_dgsem/elixir_euler_quasi_1d_discontinuous.jl
+++ b/examples/tree_1d_dgsem/elixir_euler_quasi_1d_discontinuous.jl
@@ -28,8 +28,8 @@ end
 
 initial_condition = initial_condition_discontinuity
 
-volume_flux = (flux_chan_etal, flux_nonconservative_chan_etal)
 surface_flux = (flux_lax_friedrichs, flux_nonconservative_chan_etal)
+volume_flux = (flux_chan_etal, flux_nonconservative_chan_etal)
 
 basis = LobattoLegendreBasis(3)
 indicator_sc = IndicatorHennemannGassner(equations, basis,

From 8be7c58a14a439c8a591602c9e04b704e0498a7b Mon Sep 17 00:00:00 2001
From: Krissh Chawla <127906314+KrisshChawla@users.noreply.github.com>
Date: Sun, 26 Nov 2023 10:43:27 -0600
Subject: [PATCH 30/39] Update
 examples/tree_1d_dgsem/elixir_euler_quasi_1d_ec.jl

Co-authored-by: Jesse Chan <1156048+jlchan@users.noreply.github.com>
---
 examples/tree_1d_dgsem/elixir_euler_quasi_1d_ec.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/examples/tree_1d_dgsem/elixir_euler_quasi_1d_ec.jl b/examples/tree_1d_dgsem/elixir_euler_quasi_1d_ec.jl
index 74efe9ce57c..ae1b2b24b62 100644
--- a/examples/tree_1d_dgsem/elixir_euler_quasi_1d_ec.jl
+++ b/examples/tree_1d_dgsem/elixir_euler_quasi_1d_ec.jl
@@ -25,8 +25,8 @@ end
 
 initial_condition = initial_condition_ec
 
-volume_flux = (flux_chan_etal, flux_nonconservative_chan_etal)
-surface_flux = volume_flux
+surface_flux = (flux_chan_etal, flux_nonconservative_chan_etal)
+volume_flux = surface_flux
 solver = DGSEM(polydeg = 4, surface_flux = surface_flux,
                volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
 

From e3f523a794b554b1beeb723601df8ed80c77b6bc Mon Sep 17 00:00:00 2001
From: Krissh Chawla <127906314+KrisshChawla@users.noreply.github.com>
Date: Sun, 26 Nov 2023 10:43:41 -0600
Subject: [PATCH 31/39] Update
 examples/tree_1d_dgsem/elixir_euler_quasi_1d_source_terms.jl

Co-authored-by: Jesse Chan <1156048+jlchan@users.noreply.github.com>
---
 examples/tree_1d_dgsem/elixir_euler_quasi_1d_source_terms.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/examples/tree_1d_dgsem/elixir_euler_quasi_1d_source_terms.jl b/examples/tree_1d_dgsem/elixir_euler_quasi_1d_source_terms.jl
index 3f49e2718f5..91bb1ba6e8c 100644
--- a/examples/tree_1d_dgsem/elixir_euler_quasi_1d_source_terms.jl
+++ b/examples/tree_1d_dgsem/elixir_euler_quasi_1d_source_terms.jl
@@ -10,8 +10,8 @@ equations = CompressibleEulerEquationsQuasi1D(1.4)
 
 initial_condition = initial_condition_convergence_test
 
-volume_flux = (flux_chan_etal, flux_nonconservative_chan_etal)
-surface_flux = volume_flux
+surface_flux = (flux_chan_etal, flux_nonconservative_chan_etal)
+volume_flux = surface_flux
 solver = DGSEM(polydeg = 4, surface_flux = surface_flux,
                volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
 

From f3db12c4d452195f369b9cdec704fc91f8893efb Mon Sep 17 00:00:00 2001
From: Jesse Chan <1156048+jlchan@users.noreply.github.com>
Date: Sun, 26 Nov 2023 11:20:55 -0600
Subject: [PATCH 32/39] Update src/equations/compressible_euler_quasi_1d.jl

Co-authored-by: Daniel Doehring <doehringd2@gmail.com>
---
 src/equations/compressible_euler_quasi_1d.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/equations/compressible_euler_quasi_1d.jl b/src/equations/compressible_euler_quasi_1d.jl
index 029d67fb546..1bab49e1ee5 100644
--- a/src/equations/compressible_euler_quasi_1d.jl
+++ b/src/equations/compressible_euler_quasi_1d.jl
@@ -261,7 +261,7 @@ Further details are available in the paper:
     rho_rr, v1_rr, p_rr, a_rr = cons2prim(u_rr, equations)
 
     # Compute the necessary mean values
-    rho_mean = Trixi.ln_mean(rho_ll, rho_rr)
+    rho_mean = ln_mean(rho_ll, rho_rr)
     # Algebraically equivalent to `inv_ln_mean(rho_ll / p_ll, rho_rr / p_rr)`
     # in exact arithmetic since
     #     log((ϱₗ/pₗ) / (ϱᵣ/pᵣ)) / (ϱₗ/pₗ - ϱᵣ/pᵣ)

From 2264693257522721cad3fdd1e59249c49f5d3c7b Mon Sep 17 00:00:00 2001
From: Jesse Chan <1156048+jlchan@users.noreply.github.com>
Date: Sun, 26 Nov 2023 11:21:03 -0600
Subject: [PATCH 33/39] Update src/equations/compressible_euler_quasi_1d.jl

Co-authored-by: Daniel Doehring <doehringd2@gmail.com>
---
 src/equations/compressible_euler_quasi_1d.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/equations/compressible_euler_quasi_1d.jl b/src/equations/compressible_euler_quasi_1d.jl
index 1bab49e1ee5..18e241f8826 100644
--- a/src/equations/compressible_euler_quasi_1d.jl
+++ b/src/equations/compressible_euler_quasi_1d.jl
@@ -266,7 +266,7 @@ Further details are available in the paper:
     # in exact arithmetic since
     #     log((ϱₗ/pₗ) / (ϱᵣ/pᵣ)) / (ϱₗ/pₗ - ϱᵣ/pᵣ)
     #   = pₗ pᵣ log((ϱₗ pᵣ) / (ϱᵣ pₗ)) / (ϱₗ pᵣ - ϱᵣ pₗ)
-    inv_rho_p_mean = p_ll * p_rr * Trixi.inv_ln_mean(rho_ll * p_rr, rho_rr * p_ll)
+    inv_rho_p_mean = p_ll * p_rr * inv_ln_mean(rho_ll * p_rr, rho_rr * p_ll)
     v1_avg = 0.5 * (v1_ll + v1_rr)
     a_v1_avg = 0.5 * (a_ll * v1_ll + a_rr * v1_rr)
     p_avg = 0.5 * (p_ll + p_rr)

From d2bddf42161fab6df20b0c321b2798f58b373197 Mon Sep 17 00:00:00 2001
From: Krissh Chawla <127906314+KrisshChawla@users.noreply.github.com>
Date: Sun, 26 Nov 2023 12:06:59 -0600
Subject: [PATCH 34/39] Update src/equations/compressible_euler_quasi_1d.jl

Co-authored-by: Daniel Doehring <doehringd2@gmail.com>
---
 src/equations/compressible_euler_quasi_1d.jl | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/equations/compressible_euler_quasi_1d.jl b/src/equations/compressible_euler_quasi_1d.jl
index 18e241f8826..be6a512b7a3 100644
--- a/src/equations/compressible_euler_quasi_1d.jl
+++ b/src/equations/compressible_euler_quasi_1d.jl
@@ -361,7 +361,8 @@ end
 end
 
 @inline function density(u, equations::CompressibleEulerEquationsQuasi1D)
-    rho = u[1] / u[4]
+    a_rho, _, _, a = u
+    rho = a_rho / a
     return rho
 end
 

From 2c4f3a45cf5581bd71fb911096dde083467ba1d6 Mon Sep 17 00:00:00 2001
From: KrisshChawla <krisshchawla@gmail.com>
Date: Sun, 26 Nov 2023 12:17:01 -0600
Subject: [PATCH 35/39] update boundary condition slip wall

---
 src/equations/compressible_euler_quasi_1d.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/equations/compressible_euler_quasi_1d.jl b/src/equations/compressible_euler_quasi_1d.jl
index be6a512b7a3..a8dcd16f536 100644
--- a/src/equations/compressible_euler_quasi_1d.jl
+++ b/src/equations/compressible_euler_quasi_1d.jl
@@ -167,7 +167,7 @@ are available in the paper:
     # compute the primitive variables
     rho_local, v_normal, p_local, a_local = cons2prim(u_inner, equations)
 
-    if isodd(direction) # flip sign of normal to make it outward pointing
+    if direction == 1 # flip sign of normal to make it outward pointing
         v_normal *= -1
     end
 

From b67db30e2b1149864a6ccdf1c8a32c7250da83d3 Mon Sep 17 00:00:00 2001
From: KrisshChawla <krisshchawla@gmail.com>
Date: Sun, 26 Nov 2023 12:21:35 -0600
Subject: [PATCH 36/39] update compressible_euler_quasi_1d.jl

---
 src/equations/compressible_euler_quasi_1d.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/equations/compressible_euler_quasi_1d.jl b/src/equations/compressible_euler_quasi_1d.jl
index a8dcd16f536..0d0dfd966f0 100644
--- a/src/equations/compressible_euler_quasi_1d.jl
+++ b/src/equations/compressible_euler_quasi_1d.jl
@@ -51,7 +51,7 @@ This affects the implementation and use of these equations in various ways:
 * Trixi.jl's visualization tools will visualize the nozzle width by default.
 """
 struct CompressibleEulerEquationsQuasi1D{RealT <: Real} <:
-       Trixi.AbstractCompressibleEulerEquations{1, 4}
+       AbstractCompressibleEulerEquations{1, 4}
     gamma::RealT               # ratio of specific heats
     inv_gamma_minus_one::RealT # = inv(gamma - 1); can be used to write slow divisions as fast multiplications
 
@@ -61,7 +61,7 @@ struct CompressibleEulerEquationsQuasi1D{RealT <: Real} <:
     end
 end
 
-have_nonconservative_terms(::CompressibleEulerEquationsQuasi1D) = Trixi.True()
+have_nonconservative_terms(::CompressibleEulerEquationsQuasi1D) = True()
 function varnames(::typeof(cons2cons), ::CompressibleEulerEquationsQuasi1D)
     ("a_rho", "a_rho_v1", "a_e", "a")
 end

From aaadbc83ccbf1c771618c5975510163446885380 Mon Sep 17 00:00:00 2001
From: Daniel Doehring <doehringd2@gmail.com>
Date: Sat, 2 Dec 2023 11:55:28 +0100
Subject: [PATCH 37/39] Update src/equations/compressible_euler_quasi_1d.jl

---
 src/equations/compressible_euler_quasi_1d.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/equations/compressible_euler_quasi_1d.jl b/src/equations/compressible_euler_quasi_1d.jl
index 0d0dfd966f0..4ba5dccc340 100644
--- a/src/equations/compressible_euler_quasi_1d.jl
+++ b/src/equations/compressible_euler_quasi_1d.jl
@@ -107,7 +107,7 @@ as defined in [`initial_condition_convergence_test`](@ref).
 @inline function source_terms_convergence_test(u, x, t,
                                                equations::CompressibleEulerEquationsQuasi1D)
     # Same settings as in `initial_condition_convergence_test`. 
-    #Derivatives calculated with ForwardDiff.jl
+    # Derivatives calculated with ForwardDiff.jl
     c = 2
     A = 0.1
     L = 2

From 4eb51bf5c4d887cf3295467401c256cfd5a8fc3d Mon Sep 17 00:00:00 2001
From: Daniel Doehring <doehringd2@gmail.com>
Date: Sat, 2 Dec 2023 11:55:41 +0100
Subject: [PATCH 38/39] Update src/equations/compressible_euler_quasi_1d.jl

---
 src/equations/compressible_euler_quasi_1d.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/equations/compressible_euler_quasi_1d.jl b/src/equations/compressible_euler_quasi_1d.jl
index 4ba5dccc340..ab661b8b35d 100644
--- a/src/equations/compressible_euler_quasi_1d.jl
+++ b/src/equations/compressible_euler_quasi_1d.jl
@@ -8,7 +8,7 @@
 @doc raw"""
     CompressibleEulerEquationsQuasi1D(gamma)
 
-The quasi-1d compressible Euler equations (see Chan et al.  https://doi.org/10.48550/arXiv.2307.12089 for details)
+The quasi-1d compressible Euler equations (see Chan et al.  [DOI: 10.48550/arXiv.2307.12089](https://doi.org/10.48550/arXiv.2307.12089)  for details)
 ```math
 \frac{\partial}{\partial t}
 \begin{pmatrix}

From db0e311b010ea2c692ca2e37191949205d0433c4 Mon Sep 17 00:00:00 2001
From: KrisshChawla <krisshchawla@gmail.com>
Date: Wed, 13 Dec 2023 16:59:26 -0600
Subject: [PATCH 39/39] remove boundary_condition_slip_wall

---
 src/equations/compressible_euler_quasi_1d.jl | 52 --------------------
 1 file changed, 52 deletions(-)

diff --git a/src/equations/compressible_euler_quasi_1d.jl b/src/equations/compressible_euler_quasi_1d.jl
index ab661b8b35d..0a543277ee4 100644
--- a/src/equations/compressible_euler_quasi_1d.jl
+++ b/src/equations/compressible_euler_quasi_1d.jl
@@ -145,58 +145,6 @@ as defined in [`initial_condition_convergence_test`](@ref).
     return SVector(du1, du2, du3, 0.0)
 end
 
-"""
-    boundary_condition_slip_wall(u_inner, orientation, direction, x, t,
-                                 surface_flux_function, equations::CompressibleEulerEquationsQuasi1D)
-Determine the boundary numerical surface flux for a slip wall condition.
-Imposes a zero normal velocity at the wall.
-Density is taken from the internal solution state and pressure is computed as an
-exact solution of a 1D Riemann problem. Further details about this boundary state
-are available in the paper:
-- J. J. W. van der Vegt and H. van der Ven (2002)
-  Slip flow boundary conditions in discontinuous Galerkin discretizations of
-  the Euler equations of gas dynamics
-  [PDF](https://reports.nlr.nl/bitstream/handle/10921/692/TP-2002-300.pdf?sequence=1)
-
-  Should be used together with [`TreeMesh`](@ref).
-"""
-@inline function boundary_condition_slip_wall(u_inner, orientation,
-                                              direction, x, t,
-                                              surface_flux_function,
-                                              equations::CompressibleEulerEquationsQuasi1D)
-    # compute the primitive variables
-    rho_local, v_normal, p_local, a_local = cons2prim(u_inner, equations)
-
-    if direction == 1 # flip sign of normal to make it outward pointing
-        v_normal *= -1
-    end
-
-    # Get the solution of the pressure Riemann problem
-    # See Section 6.3.3 of
-    # Eleuterio F. Toro (2009)
-    # Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction
-    # [DOI: 10.1007/b79761](https://doi.org/10.1007/b79761)
-    if v_normal <= 0.0
-        sound_speed = sqrt(equations.gamma * p_local / rho_local) # local sound speed
-        p_star = p_local *
-                 (1 + 0.5 * (equations.gamma - 1) * v_normal / sound_speed)^(2 *
-                                                                             equations.gamma *
-                                                                             equations.inv_gamma_minus_one)
-    else # v_normal > 0.0
-        A = 2 / ((equations.gamma + 1) * rho_local)
-        B = p_local * (equations.gamma - 1) / (equations.gamma + 1)
-        p_star = p_local +
-                 0.5 * v_normal / A *
-                 (v_normal + sqrt(v_normal^2 + 4 * A * (p_local + B)))
-    end
-
-    # For the slip wall we directly set the flux as the normal velocity is zero
-    return SVector(zero(eltype(u_inner)),
-                   p_star,
-                   zero(eltype(u_inner)),
-                   zero(eltype(u_inner)))
-end
-
 # Calculate 1D flux for a single point
 @inline function flux(u, orientation::Integer,
                       equations::CompressibleEulerEquationsQuasi1D)