From f3f397f520a362dbea6a7c67d1881881682515b8 Mon Sep 17 00:00:00 2001
From: KrisshChawla <krisshchawla@gmail.com>
Date: Tue, 29 Aug 2023 09:20:30 -0500
Subject: [PATCH 01/25] implementation of quasi shallow water equations 1d.

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

diff --git a/src/equations/shallow_water_quasi_1d.jl b/src/equations/shallow_water_quasi_1d.jl
new file mode 100644
index 00000000000..8b2f4358d84
--- /dev/null
+++ b/src/equations/shallow_water_quasi_1d.jl
@@ -0,0 +1,352 @@
+# 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"""
+    QuasiShallowWaterEquations1D(; gravity, H0 = 0, threshold_limiter = nothing threshold_wet = nothing)
+
+Quasi Shallow water equations (SWE) in one space dimension. The equations are given by
+```math
+\begin{aligned}
+  \frac{\partial}{\partial t}(a h) + \frac{\partial}{\partial x}(a h v) &= 0 \\
+    \frac{\partial}{\partial t}(a h v) + \frac{\partial}{\partial x}(a h v^2)
+    + g a h \frac{\partial}{\partial x}(h + b) &= 0
+\end{aligned}
+```
+The unknown quantities of the Quasi SWE are the water height ``h`` and the velocity ``v``.
+The gravitational constant is denoted by `g`, the (possibly) variable bottom topography function ``b(x)``, and channel width ''a(x)''. Conservative variable water height ``h`` is measured from the bottom topography ``b``, therefore one also defines the total water height as ``H = h + b``.
+
+The additional quantity ``H_0`` is also available to store a reference value for the total water height that
+is useful to set initial conditions or test the "lake-at-rest" well-balancedness.
+
+Also, there are two thresholds which prevent numerical problems as well as instabilities. Both of them do not
+have to be passed, as default values are defined within the struct. The first one, `threshold_limiter`, is
+used in [`PositivityPreservingLimiterShallowWater`](@ref) on the water height, as a (small) shift on the initial
+condition and cutoff before the next time step. The second one, `threshold_wet`, is applied on the water height to
+define when the flow is "wet" before calculating the numerical flux.
+
+The bottom topography function ``b(x)`` is set inside the initial condition routine
+for a particular problem setup. To test the conservative form of the SWE one can set the bottom topography
+variable `b` to zero and "a" to one. 
+
+In addition to the unknowns, Trixi.jl currently stores the bottom topography and channel width values at the approximation points 
+despite being fixed in time. This is done for convenience of computing the bottom topography gradients
+on the fly during the approximation as well as computing auxiliary quantities like the total water height ``H``
+or the entropy variables.
+This affects the implementation and use of these equations in various ways:
+* The flux values corresponding to the bottom topography and channel width must be zero.
+* The bottom topography and channel 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 bottom topography and channel width by default.
+"""
+
+struct QuasiShallowWaterEquations1D{RealT <: Real} <: Trixi.AbstractShallowWaterEquations{1, 4}
+    gravity::RealT # gravitational constant
+    H0::RealT      # constant "lake-at-rest" total water height
+    # `threshold_limiter` used in `PositivityPreservingLimiterShallowWater` on water height,
+    # as a (small) shift on the initial condition and cutoff before the next time step.
+    # Default is 500*eps() which in double precision is ≈1e-13.
+    threshold_limiter::RealT
+    # `threshold_wet` applied on water height to define when the flow is "wet"
+    # before calculating the numerical flux.
+    # Default is 5*eps() which in double precision is ≈1e-15.
+    threshold_wet::RealT
+end
+
+# Allow for flexibility to set the gravitational constant within an elixir depending on the
+# application where `gravity_constant=1.0` or `gravity_constant=9.81` are common values.
+# The reference total water height H0 defaults to 0.0 but is used for the "lake-at-rest"
+# well-balancedness test cases.
+# Strict default values for thresholds that performed well in many numerical experiments
+function QuasiShallowWaterEquations1D(; gravity_constant, H0 = zero(gravity_constant),
+                                 threshold_limiter = nothing, threshold_wet = nothing)
+    T = promote_type(typeof(gravity_constant), typeof(H0))
+    if threshold_limiter === nothing
+        threshold_limiter = 500 * eps(T)
+    end
+    if threshold_wet === nothing
+        threshold_wet = 5 * eps(T)
+    end
+    QuasiShallowWaterEquations1D(gravity_constant, H0, threshold_limiter, threshold_wet)
+end
+
+have_nonconservative_terms(::QuasiShallowWaterEquations1D) = Trixi.True()
+varnames(::typeof(cons2cons), ::QuasiShallowWaterEquations1D) = ("a_h", "a_h_v", "b", "a")
+# Note, we use the total water height, H = h + b, as the first primitive variable for easier
+# visualization and setting initial conditions
+varnames(::typeof(cons2prim), ::QuasiShallowWaterEquations1D) = ("H", "v", "b", "a")
+
+# Set initial conditions at physical location `x` for time `t`
+#need to revise
+"""
+    initial_condition_convergence_test(x, t, equations::QuasiShallowWaterEquations1D)
+
+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_convergence_test(x, t, equations::QuasiShallowWaterEquations1D)
+    #generate manufactured fucntion a(x)
+    # some constants are chosen such that the function is periodic on the domain [0,sqrt(2)]
+    Omega = sqrt(2) * pi 
+    H = 2.0 + 0.5 * sin(Omega * x[1] - t)
+    v = 0.25
+    b = 0.2 - 0.05 * sin(1 * sqrt(2) * pi *x[1])
+    a = 1 + 0.1 * cos(sqrt(2) * pi * x[1])   
+    return prim2cons(SVector(H, v, b, a), equations)
+end
+
+"""
+    source_terms_convergence_test(u, x, t, equations::QuasiShallowWaterEquations1D)
+
+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 bottom topography function
+`b(x) = 0.2 - 0.05 * sin(sqrt(2) * pi *x[1])` and channel width 'a(x)= 1 + 0.1 * cos(sqrt(2) * pi * x[1])'
+as defined in [`initial_condition_convergence_test`](@ref).
+"""
+
+@inline function source_terms_convergence_test(u, x, t,
+                                               equations::QuasiShallowWaterEquations1D)
+    # Same settings as in `initial_condition_convergence_test`. Some derivative simplify because
+    # this manufactured solution velocity is taken to be constant
+    Omega = sqrt(2) * pi 
+    H = 2.0 + 0.5 * sin(Omega * x[1] - t)
+    H_x = 0.5 * cos(Omega * x[1] - t) * Omega
+    H_t = -0.5 * cos(Omega * x[1] - t)
+
+    v = 0.25
+    
+    b =  0.2 - 0.05 * sin(sqrt(2) * pi *x[1])
+    b_x = -0.05 * cos(sqrt(2) * pi *x[1]) * sqrt(2) * pi
+    
+    a = 1 + 0.1 * cos(sqrt(2) * pi * x[1])
+    a_x = -0.1 * sin(sqrt(2) * pi * x[1]) * sqrt(2) * pi
+
+    du1 = a * H_t + v * (a_x * (H - b) + a * (H_x - b_x))
+    du2 = v * du1 + a * (equations.gravity * (H - b) * H_x)
+
+    return SVector(du1, du2, 0.0, 0.0)
+end  
+
+"""
+    boundary_condition_slip_wall(u_inner, orientation_or_normal, x, t, surface_flux_function,
+                                 equations::QuasiShallowWaterEquations1D)
+
+Create a boundary state by reflecting the normal velocity component and keep
+the tangential velocity component unchanged. The boundary water height is taken from
+the internal value.
+
+For details see Section 9.2.5 of the book:
+- Eleuterio F. Toro (2001)
+  Shock-Capturing Methods for Free-Surface Shallow Flows
+  1st edition
+  ISBN 0471987662
+"""
+
+function boundary_condition_slip_wall(u_inner, orientation_or_normal, direction,
+                                              x, t, surface_flux_function,
+                                              equations::QuasiShallowWaterEquations1D)
+    # create the "external" boundary solution state
+    u_boundary = SVector(u_inner[1],
+                         -u_inner[2],
+                         u_inner[3],
+                         u_inner[4])
+
+    # calculate the boundary flux
+    if iseven(direction) # u_inner is "left" of boundary, u_boundary is "right" of boundary
+        f = surface_flux_function(u_inner, u_boundary, orientation_or_normal, equations)
+    else # u_boundary is "left" of boundary, u_inner is "right" of boundary
+        f = surface_flux_function(u_boundary, u_inner, orientation_or_normal, equations)
+    end
+    return f
+end
+
+# Calculate 1D flux for a single point
+# Note, the bottom topography and channel width have no flux
+@inline function flux(u, orientation::Integer, equations::QuasiShallowWaterEquations1D)
+    a_h, a_h_v, _, a = u
+    h = waterheight(u, equations)
+    v = velocity(u, equations)
+
+    p = 0.5 * a * equations.gravity * h^2
+    
+    f1 = a_h_v
+    f2 = a_h_v * v + p
+
+    return SVector(f1, f2, zero(eltype(u)), zero(eltype(u)))
+end
+
+"""
+flux_nonconservative_chan_etal(u_ll, u_rr, orientation::Integer,
+                                                equations::QuasiShallowWaterEquations1D)
+
+Non-symmetric two-point volume flux discretizing the nonconservative (source) term
+that contains the gradient of the bottom topography [`QuasiShallowWaterEquations1D`](@ref) 
+and the channel 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::QuasiShallowWaterEquations1D)
+    a_h_ll, _, b_ll, a_ll = u_ll
+    a_h_rr, _, b_rr, a_rr = u_rr
+    
+    h_ll = waterheight(u_ll, equations)
+    h_rr = waterheight(u_rr, equations)
+    
+    z = zero(eltype(u_ll))
+    
+    return SVector(z, equations.gravity * a_ll * h_ll * (h_rr + b_rr), z, z)
+end
+
+"""
+    flux_chan_etal(u_ll, u_rr, orientation,
+                          equations::QuasiShallowWaterEquations1D)
+
+Total energy conservative (mathematical entropy for quasi shallow water equations) split form.
+When the bottom topography is nonzero this scheme will be well-balanced when used as a `volume_flux`.
+The `surface_flux` should still use, e.g., [`FluxPlusDissipation(flux_chan_etal, DissipationLocalLaxFriedrichs())`](@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
+  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::QuasiShallowWaterEquations1D)
+    a_h_ll, a_h_v_ll, _, _ = u_ll
+    a_h_rr, a_h_v_rr, _, _ = u_rr
+    
+    v_ll = velocity(u_ll, equations)
+    v_rr = velocity(u_rr, equations)
+
+    f1 = 0.5 * (a_h_v_ll + a_h_v_rr)
+    f2 = f1 * 0.5 * (v_ll + v_rr)
+    
+    return SVector(f1, f2, zero(eltype(u_ll)), zero(eltype(u_ll)))
+end
+
+# Calculate 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::QuasiShallowWaterEquations1D)
+    # Get the velocity quantities
+    v_ll = velocity(u_ll, equations)
+    v_rr = velocity(u_rr, equations)
+
+    # Calculate the wave celerity on the left and right
+    h_ll = waterheight(u_ll, equations)
+    h_rr = waterheight(u_rr, equations)
+    c_ll = sqrt(equations.gravity * h_ll)
+    c_rr = sqrt(equations.gravity * h_rr)
+
+    return max(abs(v_ll), abs(v_rr)) + max(c_ll, c_rr)
+end
+
+# Specialized `DissipationLocalLaxFriedrichs` to avoid spurious dissipation in the bottom topography
+@inline function (dissipation::DissipationLocalLaxFriedrichs)(u_ll, u_rr,
+                                                              orientation_or_normal_direction,
+                                                              equations::QuasiShallowWaterEquations1D)
+    λ = dissipation.max_abs_speed(u_ll, u_rr, orientation_or_normal_direction,
+                                  equations)
+    diss = -0.5 * λ * (u_rr - u_ll)
+    return SVector(diss[1], diss[2], zero(eltype(u_ll)), zero(eltype(u_ll)))
+end
+
+# Helper function to extract the velocity vector from the conservative variables
+@inline function velocity(u, equations::QuasiShallowWaterEquations1D)
+    a_h, a_h_v, _, _ = u
+    
+    v = a_h_v / a_h
+    
+    return v
+end
+
+# Convert conservative variables to primitive
+@inline function cons2prim(u, equations::QuasiShallowWaterEquations1D)
+    a_h, _, b, a = u
+    h = a_h / a
+    H = h + b
+    v = velocity(u, equations)
+    return SVector(H, v, b, a)
+end
+
+#Convert conservative variables to entropy variables
+# Note, only the first two are the entropy variables, the third and foruth entries still
+# just carry the bottom topography and channel width values for convenience
+@inline function Trixi.cons2entropy(u, equations::QuasiShallowWaterEquations1D)
+    a_h, a_h_v, b, a = u
+    h = waterheight(u, equations)
+    v = velocity(u, equations)
+    #entropy variables are the same as ones in standard equations
+    w1 = equations.gravity * (h + b) - 0.5 * v^2
+    w2 = v
+    
+    return SVector(w1, w2, b, a)
+end
+
+# Convert primitive to conservative variables
+@inline function prim2cons(prim, equations::QuasiShallowWaterEquations1D)
+    H, v, b, a = prim
+    
+    a_h = a * (H - b)
+    a_h_v = a_h * v
+    return SVector(a_h, a_h_v, b, a)
+end
+
+@inline function waterheight(u, equations::QuasiShallowWaterEquations1D)
+    return u[1] / u[4]
+end
+
+@inline function pressure(u, equations::QuasiShallowWaterEquations1D)
+    a = u[4]
+    h = waterheight(u, equations)
+    p = 0.5 * a * equations.gravity * h^2
+    return p
+end
+
+# Entropy function for the shallow water equations is the total energy
+@inline function entropy(cons, equations::QuasiShallowWaterEquations1D)
+    energy_total(cons, equations)
+end
+
+# Calculate total energy for a conservative state `cons`
+@inline function energy_total(cons, equations::QuasiShallowWaterEquations1D)
+    a_h, a_h_v, b, a = cons
+    e = (a_h_v^2) / (2 * a_h) + 0.5 * equations.gravity * (a_h^2 / a) + equations.gravity * a_h * b
+    return e
+end
+
+# Calculate kinetic energy for a conservative state `cons`
+@inline function energy_kinetic(u, equations::QuasiShallowWaterEquations1D)
+    a_h, a_h_v, _, _ = u
+    return (a_h_v^2) / (2 * a_h)
+end
+
+# Calculate potential energy for a conservative state `cons`
+@inline function energy_internal(cons, equations::QuasiShallowWaterEquations1D)
+    return energy_total(cons, equations) - energy_kinetic(cons, equations)
+end
+
+# Calculate the error for the "lake-at-rest" test case
+@inline function lake_at_rest_error(u, equations::QuasiShallowWaterEquations1D)
+    a_h, _, b, _ = u
+    h = waterheight(u, equations)
+    H0_wet_dry = max(equations.H0, b + equations.threshold_limiter)
+    return abs(H0_wet_dry - (h + b))
+end
+end # @muladd
\ No newline at end of file

From a97dabdd910eb0db8241320bd5307b47430e44f8 Mon Sep 17 00:00:00 2001
From: KrisshChawla <krisshchawla@gmail.com>
Date: Tue, 29 Aug 2023 09:28:59 -0500
Subject: [PATCH 02/25] added example elixer for shallow_water_quasi_1d

---
 ...xir_shallow_water_quasi_1D_manufactured.jl | 59 +++++++++++++++++++
 1 file changed, 59 insertions(+)
 create mode 100644 examples/tree_1d_dgsem/elixir_shallow_water_quasi_1D_manufactured.jl

diff --git a/examples/tree_1d_dgsem/elixir_shallow_water_quasi_1D_manufactured.jl b/examples/tree_1d_dgsem/elixir_shallow_water_quasi_1D_manufactured.jl
new file mode 100644
index 00000000000..7ef87e70fa8
--- /dev/null
+++ b/examples/tree_1d_dgsem/elixir_shallow_water_quasi_1D_manufactured.jl
@@ -0,0 +1,59 @@
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# Semidiscretization of the shallow water equations
+
+equations = QuasiShallowWaterEquations1D(gravity_constant=9.81)
+
+initial_condition = initial_condition_convergence_test
+
+
+###############################################################################
+# Get the DG approximation space
+
+volume_flux = (flux_chan_etal, flux_nonconservative_chan_etal) 
+surface_flux = (FluxPlusDissipation(flux_chan_etal, DissipationLocalLaxFriedrichs()) , flux_nonconservative_chan_etal)
+solver = DGSEM(polydeg=3, surface_flux=surface_flux,
+               volume_integral=VolumeIntegralFluxDifferencing(volume_flux))
+
+###############################################################################
+# Get the TreeMesh and setup a periodic mesh
+
+coordinates_min = 0.0
+coordinates_max = sqrt(2.0)
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level=3,
+                n_cells_max=10_000,
+                periodicity=true)
+
+# create the semi discretization object
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    source_terms=source_terms_convergence_test)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 1.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 500
+analysis_callback = AnalysisCallback(semi, interval=analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval=analysis_interval)
+
+save_solution = SaveSolutionCallback(interval=200,
+                                     save_initial_solution=true,
+                                     save_final_solution=true)
+
+callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution)
+
+###############################################################################
+# run the simulation
+
+# use a Runge-Kutta method with automatic (error based) time step size control
+sol = solve(ode, RDPK3SpFSAL49(); abstol=1.0e-8, reltol=1.0e-8,
+            ode_default_options()..., callback=callbacks);
+summary_callback() # print the timer summary
\ No newline at end of file

From b340a74872b00e8206806c1eced560dd849fb4b0 Mon Sep 17 00:00:00 2001
From: KrisshChawla <krisshchawla@gmail.com>
Date: Tue, 29 Aug 2023 17:03:16 -0500
Subject: [PATCH 03/25] changed the names of Quasi1d equations

---
 ...xir_shallow_water_quasi_1D_manufactured.jl |  2 +-
 src/equations/shallow_water_quasi_1d.jl       | 64 +++++++++----------
 2 files changed, 33 insertions(+), 33 deletions(-)

diff --git a/examples/tree_1d_dgsem/elixir_shallow_water_quasi_1D_manufactured.jl b/examples/tree_1d_dgsem/elixir_shallow_water_quasi_1D_manufactured.jl
index 7ef87e70fa8..4f906a79d27 100644
--- a/examples/tree_1d_dgsem/elixir_shallow_water_quasi_1D_manufactured.jl
+++ b/examples/tree_1d_dgsem/elixir_shallow_water_quasi_1D_manufactured.jl
@@ -4,7 +4,7 @@ using Trixi
 ###############################################################################
 # Semidiscretization of the shallow water equations
 
-equations = QuasiShallowWaterEquations1D(gravity_constant=9.81)
+equations = ShallowWaterEquationsQuasi1D(gravity_constant=9.81)
 
 initial_condition = initial_condition_convergence_test
 
diff --git a/src/equations/shallow_water_quasi_1d.jl b/src/equations/shallow_water_quasi_1d.jl
index 8b2f4358d84..f9cdd02bc4c 100644
--- a/src/equations/shallow_water_quasi_1d.jl
+++ b/src/equations/shallow_water_quasi_1d.jl
@@ -6,7 +6,7 @@
 #! format: noindent
     
 @doc raw"""
-    QuasiShallowWaterEquations1D(; gravity, H0 = 0, threshold_limiter = nothing threshold_wet = nothing)
+    ShallowWaterEquationsQuasi1D(; gravity, H0 = 0, threshold_limiter = nothing threshold_wet = nothing)
 
 Quasi Shallow water equations (SWE) in one space dimension. The equations are given by
 ```math
@@ -44,7 +44,7 @@ This affects the implementation and use of these equations in various ways:
 * Trixi.jl's visualization tools will visualize the bottom topography and channel width by default.
 """
 
-struct QuasiShallowWaterEquations1D{RealT <: Real} <: Trixi.AbstractShallowWaterEquations{1, 4}
+struct ShallowWaterEquationsQuasi1D{RealT <: Real} <: Trixi.AbstractShallowWaterEquations{1, 4}
     gravity::RealT # gravitational constant
     H0::RealT      # constant "lake-at-rest" total water height
     # `threshold_limiter` used in `PositivityPreservingLimiterShallowWater` on water height,
@@ -62,7 +62,7 @@ end
 # The reference total water height H0 defaults to 0.0 but is used for the "lake-at-rest"
 # well-balancedness test cases.
 # Strict default values for thresholds that performed well in many numerical experiments
-function QuasiShallowWaterEquations1D(; gravity_constant, H0 = zero(gravity_constant),
+function ShallowWaterEquationsQuasi1D(; gravity_constant, H0 = zero(gravity_constant),
                                  threshold_limiter = nothing, threshold_wet = nothing)
     T = promote_type(typeof(gravity_constant), typeof(H0))
     if threshold_limiter === nothing
@@ -71,25 +71,25 @@ function QuasiShallowWaterEquations1D(; gravity_constant, H0 = zero(gravity_cons
     if threshold_wet === nothing
         threshold_wet = 5 * eps(T)
     end
-    QuasiShallowWaterEquations1D(gravity_constant, H0, threshold_limiter, threshold_wet)
+    ShallowWaterEquationsQuasi1D(gravity_constant, H0, threshold_limiter, threshold_wet)
 end
 
-have_nonconservative_terms(::QuasiShallowWaterEquations1D) = Trixi.True()
-varnames(::typeof(cons2cons), ::QuasiShallowWaterEquations1D) = ("a_h", "a_h_v", "b", "a")
+have_nonconservative_terms(::ShallowWaterEquationsQuasi1D) = Trixi.True()
+varnames(::typeof(cons2cons), ::ShallowWaterEquationsQuasi1D) = ("a_h", "a_h_v", "b", "a")
 # Note, we use the total water height, H = h + b, as the first primitive variable for easier
 # visualization and setting initial conditions
-varnames(::typeof(cons2prim), ::QuasiShallowWaterEquations1D) = ("H", "v", "b", "a")
+varnames(::typeof(cons2prim), ::ShallowWaterEquationsQuasi1D) = ("H", "v", "b", "a")
 
 # Set initial conditions at physical location `x` for time `t`
 #need to revise
 """
-    initial_condition_convergence_test(x, t, equations::QuasiShallowWaterEquations1D)
+    initial_condition_convergence_test(x, t, equations::ShallowWaterEquationsQuasi1D)
 
 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_convergence_test(x, t, equations::QuasiShallowWaterEquations1D)
+function initial_condition_convergence_test(x, t, equations::ShallowWaterEquationsQuasi1D)
     #generate manufactured fucntion a(x)
     # some constants are chosen such that the function is periodic on the domain [0,sqrt(2)]
     Omega = sqrt(2) * pi 
@@ -101,7 +101,7 @@ function initial_condition_convergence_test(x, t, equations::QuasiShallowWaterEq
 end
 
 """
-    source_terms_convergence_test(u, x, t, equations::QuasiShallowWaterEquations1D)
+    source_terms_convergence_test(u, x, t, equations::ShallowWaterEquationsQuasi1D)
 
 Source terms used for convergence tests in combination with
 [`initial_condition_convergence_test`](@ref)
@@ -113,7 +113,7 @@ as defined in [`initial_condition_convergence_test`](@ref).
 """
 
 @inline function source_terms_convergence_test(u, x, t,
-                                               equations::QuasiShallowWaterEquations1D)
+                                               equations::ShallowWaterEquationsQuasi1D)
     # Same settings as in `initial_condition_convergence_test`. Some derivative simplify because
     # this manufactured solution velocity is taken to be constant
     Omega = sqrt(2) * pi 
@@ -137,7 +137,7 @@ end
 
 """
     boundary_condition_slip_wall(u_inner, orientation_or_normal, x, t, surface_flux_function,
-                                 equations::QuasiShallowWaterEquations1D)
+                                 equations::ShallowWaterEquationsQuasi1D)
 
 Create a boundary state by reflecting the normal velocity component and keep
 the tangential velocity component unchanged. The boundary water height is taken from
@@ -152,7 +152,7 @@ For details see Section 9.2.5 of the book:
 
 function boundary_condition_slip_wall(u_inner, orientation_or_normal, direction,
                                               x, t, surface_flux_function,
-                                              equations::QuasiShallowWaterEquations1D)
+                                              equations::ShallowWaterEquationsQuasi1D)
     # create the "external" boundary solution state
     u_boundary = SVector(u_inner[1],
                          -u_inner[2],
@@ -170,7 +170,7 @@ end
 
 # Calculate 1D flux for a single point
 # Note, the bottom topography and channel width have no flux
-@inline function flux(u, orientation::Integer, equations::QuasiShallowWaterEquations1D)
+@inline function flux(u, orientation::Integer, equations::ShallowWaterEquationsQuasi1D)
     a_h, a_h_v, _, a = u
     h = waterheight(u, equations)
     v = velocity(u, equations)
@@ -185,10 +185,10 @@ end
 
 """
 flux_nonconservative_chan_etal(u_ll, u_rr, orientation::Integer,
-                                                equations::QuasiShallowWaterEquations1D)
+                                                equations::ShallowWaterEquationsQuasi1D)
 
 Non-symmetric two-point volume flux discretizing the nonconservative (source) term
-that contains the gradient of the bottom topography [`QuasiShallowWaterEquations1D`](@ref) 
+that contains the gradient of the bottom topography [`ShallowWaterEquationsQuasi1D`](@ref) 
 and the channel width.
 
 Further details are available in the paper:
@@ -199,7 +199,7 @@ Further details are available in the paper:
 """
 
 @inline function flux_nonconservative_chan_etal(u_ll, u_rr, orientation::Integer,
-                                                equations::QuasiShallowWaterEquations1D)
+                                                equations::ShallowWaterEquationsQuasi1D)
     a_h_ll, _, b_ll, a_ll = u_ll
     a_h_rr, _, b_rr, a_rr = u_rr
     
@@ -213,7 +213,7 @@ end
 
 """
     flux_chan_etal(u_ll, u_rr, orientation,
-                          equations::QuasiShallowWaterEquations1D)
+                          equations::ShallowWaterEquationsQuasi1D)
 
 Total energy conservative (mathematical entropy for quasi shallow water equations) split form.
 When the bottom topography is nonzero this scheme will be well-balanced when used as a `volume_flux`.
@@ -227,7 +227,7 @@ Further details are available in the paper:
 """
 
 @inline function flux_chan_etal(u_ll, u_rr, orientation::Integer, 
-                                equations::QuasiShallowWaterEquations1D)
+                                equations::ShallowWaterEquationsQuasi1D)
     a_h_ll, a_h_v_ll, _, _ = u_ll
     a_h_rr, a_h_v_rr, _, _ = u_rr
     
@@ -243,7 +243,7 @@ end
 # Calculate 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::QuasiShallowWaterEquations1D)
+                                     equations::ShallowWaterEquationsQuasi1D)
     # Get the velocity quantities
     v_ll = velocity(u_ll, equations)
     v_rr = velocity(u_rr, equations)
@@ -260,7 +260,7 @@ end
 # Specialized `DissipationLocalLaxFriedrichs` to avoid spurious dissipation in the bottom topography
 @inline function (dissipation::DissipationLocalLaxFriedrichs)(u_ll, u_rr,
                                                               orientation_or_normal_direction,
-                                                              equations::QuasiShallowWaterEquations1D)
+                                                              equations::ShallowWaterEquationsQuasi1D)
     λ = dissipation.max_abs_speed(u_ll, u_rr, orientation_or_normal_direction,
                                   equations)
     diss = -0.5 * λ * (u_rr - u_ll)
@@ -268,7 +268,7 @@ end
 end
 
 # Helper function to extract the velocity vector from the conservative variables
-@inline function velocity(u, equations::QuasiShallowWaterEquations1D)
+@inline function velocity(u, equations::ShallowWaterEquationsQuasi1D)
     a_h, a_h_v, _, _ = u
     
     v = a_h_v / a_h
@@ -277,7 +277,7 @@ end
 end
 
 # Convert conservative variables to primitive
-@inline function cons2prim(u, equations::QuasiShallowWaterEquations1D)
+@inline function cons2prim(u, equations::ShallowWaterEquationsQuasi1D)
     a_h, _, b, a = u
     h = a_h / a
     H = h + b
@@ -288,7 +288,7 @@ end
 #Convert conservative variables to entropy variables
 # Note, only the first two are the entropy variables, the third and foruth entries still
 # just carry the bottom topography and channel width values for convenience
-@inline function Trixi.cons2entropy(u, equations::QuasiShallowWaterEquations1D)
+@inline function Trixi.cons2entropy(u, equations::ShallowWaterEquationsQuasi1D)
     a_h, a_h_v, b, a = u
     h = waterheight(u, equations)
     v = velocity(u, equations)
@@ -300,7 +300,7 @@ end
 end
 
 # Convert primitive to conservative variables
-@inline function prim2cons(prim, equations::QuasiShallowWaterEquations1D)
+@inline function prim2cons(prim, equations::ShallowWaterEquationsQuasi1D)
     H, v, b, a = prim
     
     a_h = a * (H - b)
@@ -308,11 +308,11 @@ end
     return SVector(a_h, a_h_v, b, a)
 end
 
-@inline function waterheight(u, equations::QuasiShallowWaterEquations1D)
+@inline function waterheight(u, equations::ShallowWaterEquationsQuasi1D)
     return u[1] / u[4]
 end
 
-@inline function pressure(u, equations::QuasiShallowWaterEquations1D)
+@inline function pressure(u, equations::ShallowWaterEquationsQuasi1D)
     a = u[4]
     h = waterheight(u, equations)
     p = 0.5 * a * equations.gravity * h^2
@@ -320,30 +320,30 @@ end
 end
 
 # Entropy function for the shallow water equations is the total energy
-@inline function entropy(cons, equations::QuasiShallowWaterEquations1D)
+@inline function entropy(cons, equations::ShallowWaterEquationsQuasi1D)
     energy_total(cons, equations)
 end
 
 # Calculate total energy for a conservative state `cons`
-@inline function energy_total(cons, equations::QuasiShallowWaterEquations1D)
+@inline function energy_total(cons, equations::ShallowWaterEquationsQuasi1D)
     a_h, a_h_v, b, a = cons
     e = (a_h_v^2) / (2 * a_h) + 0.5 * equations.gravity * (a_h^2 / a) + equations.gravity * a_h * b
     return e
 end
 
 # Calculate kinetic energy for a conservative state `cons`
-@inline function energy_kinetic(u, equations::QuasiShallowWaterEquations1D)
+@inline function energy_kinetic(u, equations::ShallowWaterEquationsQuasi1D)
     a_h, a_h_v, _, _ = u
     return (a_h_v^2) / (2 * a_h)
 end
 
 # Calculate potential energy for a conservative state `cons`
-@inline function energy_internal(cons, equations::QuasiShallowWaterEquations1D)
+@inline function energy_internal(cons, equations::ShallowWaterEquationsQuasi1D)
     return energy_total(cons, equations) - energy_kinetic(cons, equations)
 end
 
 # Calculate the error for the "lake-at-rest" test case
-@inline function lake_at_rest_error(u, equations::QuasiShallowWaterEquations1D)
+@inline function lake_at_rest_error(u, equations::ShallowWaterEquationsQuasi1D)
     a_h, _, b, _ = u
     h = waterheight(u, equations)
     H0_wet_dry = max(equations.H0, b + equations.threshold_limiter)

From af25d90113cc84cc70b89e6485082956381baa14 Mon Sep 17 00:00:00 2001
From: KrisshChawla <krisshchawla@gmail.com>
Date: Tue, 29 Aug 2023 17:04:55 -0500
Subject: [PATCH 04/25] including and exported ShallowWaterEquationsQuasi1D

---
 src/Trixi.jl               | 3 ++-
 src/equations/equations.jl | 1 +
 2 files changed, 3 insertions(+), 1 deletion(-)

diff --git a/src/Trixi.jl b/src/Trixi.jl
index ec4d20558e5..63a773d5aa3 100644
--- a/src/Trixi.jl
+++ b/src/Trixi.jl
@@ -148,7 +148,8 @@ export AcousticPerturbationEquations2D,
        InviscidBurgersEquation1D,
        LatticeBoltzmannEquations2D, LatticeBoltzmannEquations3D,
        ShallowWaterEquations1D, ShallowWaterEquations2D,
-       ShallowWaterTwoLayerEquations1D, ShallowWaterTwoLayerEquations2D,
+       ShallowWaterTwoLayerEquations1D, ShallowWaterTwoLayerEquations2D, 
+       ShallowWaterEquationsQuasi1D,
        LinearizedEulerEquations2D
 
 export LaplaceDiffusion1D, LaplaceDiffusion2D,
diff --git a/src/equations/equations.jl b/src/equations/equations.jl
index 90b2cd62191..aea0ef6277f 100644
--- a/src/equations/equations.jl
+++ b/src/equations/equations.jl
@@ -350,6 +350,7 @@ include("shallow_water_1d.jl")
 include("shallow_water_2d.jl")
 include("shallow_water_two_layer_1d.jl")
 include("shallow_water_two_layer_2d.jl")
+include("shallow_water_quasi_1d.jl")
 
 # CompressibleEulerEquations
 abstract type AbstractCompressibleEulerEquations{NDIMS, NVARS} <:

From 5defbf6cedd1345321e7aa8fdb7b916b706bd617 Mon Sep 17 00:00:00 2001
From: Jesse Chan <jesse.lee.chan@gmail.com>
Date: Tue, 29 Aug 2023 17:39:39 -0500
Subject: [PATCH 05/25] exporting flux_chan_etal and
 flux_chan_nonconservative_etal

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

diff --git a/src/Trixi.jl b/src/Trixi.jl
index 63a773d5aa3..63c0812cd48 100644
--- a/src/Trixi.jl
+++ b/src/Trixi.jl
@@ -165,6 +165,7 @@ export flux, flux_central, flux_lax_friedrichs, flux_hll, flux_hllc, flux_hlle,
        flux_kennedy_gruber, flux_shima_etal, flux_ec,
        flux_fjordholm_etal, flux_nonconservative_fjordholm_etal, flux_es_fjordholm_etal,
        flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal,
+       flux_chan_etal, flux_nonconservative_chan_etal,
        hydrostatic_reconstruction_audusse_etal, flux_nonconservative_audusse_etal,
 # TODO: TrixiShallowWater: move anything with "chen_noelle" to new file
        hydrostatic_reconstruction_chen_noelle, flux_nonconservative_chen_noelle,

From c0455707e4a741a5de2c397411ec15d52eb17a8f Mon Sep 17 00:00:00 2001
From: Jesse Chan <jesse.lee.chan@gmail.com>
Date: Tue, 29 Aug 2023 17:39:47 -0500
Subject: [PATCH 06/25] minor comment fix

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

diff --git a/src/equations/shallow_water_quasi_1d.jl b/src/equations/shallow_water_quasi_1d.jl
index f9cdd02bc4c..6924bab0766 100644
--- a/src/equations/shallow_water_quasi_1d.jl
+++ b/src/equations/shallow_water_quasi_1d.jl
@@ -90,7 +90,7 @@ A smooth initial condition used for convergence tests in combination with
 (and [`BoundaryConditionDirichlet(initial_condition_convergence_test)`](@ref) in non-periodic domains).
 """
 function initial_condition_convergence_test(x, t, equations::ShallowWaterEquationsQuasi1D)
-    #generate manufactured fucntion a(x)
+    # generates a manufactured solution. 
     # some constants are chosen such that the function is periodic on the domain [0,sqrt(2)]
     Omega = sqrt(2) * pi 
     H = 2.0 + 0.5 * sin(Omega * x[1] - t)

From 09b2e4ea0fe8d54c85b4a97b162448b3e4c75a4e Mon Sep 17 00:00:00 2001
From: Jesse Chan <jesse.lee.chan@gmail.com>
Date: Tue, 29 Aug 2023 17:39:51 -0500
Subject: [PATCH 07/25] adding tests

---
 test/test_tree_1d_shallowwater.jl | 7 +++++++
 1 file changed, 7 insertions(+)

diff --git a/test/test_tree_1d_shallowwater.jl b/test/test_tree_1d_shallowwater.jl
index cafa17edd4c..b6332181af9 100644
--- a/test/test_tree_1d_shallowwater.jl
+++ b/test/test_tree_1d_shallowwater.jl
@@ -111,6 +111,13 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem")
       linf = [0.00041080213807871235, 0.00014823261488938177, 2.220446049250313e-16],
       tspan = (0.0, 0.05))
   end
+
+  @trixi_testset "elixir_shallow_water_quasi_1D_manufactured.jl" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallow_water_quasi_1D_manufactured.jl"),
+      l2 = [6.37048760275098e-5, 0.0002745658116815704, 4.436491725647962e-6, 8.872983451152218e-6], 
+      linf = [0.00026747526881631956, 0.0012106730729152249, 9.098379777500165e-6, 1.8196759554278685e-5]
+      tspan = (0.0, 0.05))
+  end  
 end
 
 end # module

From 71fab3ee846ef1fa39a9f17d486de58e9d6e9a97 Mon Sep 17 00:00:00 2001
From: Jesse Chan <1156048+jlchan@users.noreply.github.com>
Date: Tue, 29 Aug 2023 17:41:10 -0500
Subject: [PATCH 08/25] Apply suggestions from code review

---
 src/equations/shallow_water_quasi_1d.jl | 29 +++++++++----------------
 1 file changed, 10 insertions(+), 19 deletions(-)

diff --git a/src/equations/shallow_water_quasi_1d.jl b/src/equations/shallow_water_quasi_1d.jl
index 6924bab0766..fa23f72dd9d 100644
--- a/src/equations/shallow_water_quasi_1d.jl
+++ b/src/equations/shallow_water_quasi_1d.jl
@@ -8,7 +8,7 @@
 @doc raw"""
     ShallowWaterEquationsQuasi1D(; gravity, H0 = 0, threshold_limiter = nothing threshold_wet = nothing)
 
-Quasi Shallow water equations (SWE) in one space dimension. The equations are given by
+The quasi-1D shallow water equations (SWE). The equations are given by
 ```math
 \begin{aligned}
   \frac{\partial}{\partial t}(a h) + \frac{\partial}{\partial x}(a h v) &= 0 \\
@@ -16,8 +16,8 @@ Quasi Shallow water equations (SWE) in one space dimension. The equations are gi
     + g a h \frac{\partial}{\partial x}(h + b) &= 0
 \end{aligned}
 ```
-The unknown quantities of the Quasi SWE are the water height ``h`` and the velocity ``v``.
-The gravitational constant is denoted by `g`, the (possibly) variable bottom topography function ``b(x)``, and channel width ''a(x)''. Conservative variable water height ``h`` is measured from the bottom topography ``b``, therefore one also defines the total water height as ``H = h + b``.
+The unknown quantities of the Quasi-1D SWE are the water height ``h`` and the scaled velocity ``v``.
+The gravitational constant is denoted by `g`, the (possibly) variable bottom topography function ``b(x)``, and (possibly) variable channel width ``a(x)``. The water height ``h`` is measured from the bottom topography ``b``, therefore one also defines the total water height as ``H = h + b``.
 
 The additional quantity ``H_0`` is also available to store a reference value for the total water height that
 is useful to set initial conditions or test the "lake-at-rest" well-balancedness.
@@ -28,9 +28,9 @@ used in [`PositivityPreservingLimiterShallowWater`](@ref) on the water height, a
 condition and cutoff before the next time step. The second one, `threshold_wet`, is applied on the water height to
 define when the flow is "wet" before calculating the numerical flux.
 
-The bottom topography function ``b(x)`` is set inside the initial condition routine
+The bottom topography function ``b(x)`` and channel width ``a(x)`` are set inside the initial condition routine
 for a particular problem setup. To test the conservative form of the SWE one can set the bottom topography
-variable `b` to zero and "a" to one. 
+variable `b` to zero and ``a`` to one. 
 
 In addition to the unknowns, Trixi.jl currently stores the bottom topography and channel width values at the approximation points 
 despite being fixed in time. This is done for convenience of computing the bottom topography gradients
@@ -312,35 +312,26 @@ end
     return u[1] / u[4]
 end
 
-@inline function pressure(u, equations::ShallowWaterEquationsQuasi1D)
-    a = u[4]
-    h = waterheight(u, equations)
-    p = 0.5 * a * equations.gravity * h^2
-    return p
-end
 
 # Entropy function for the shallow water equations is the total energy
 @inline function entropy(cons, equations::ShallowWaterEquationsQuasi1D)
-    energy_total(cons, equations)
+    a = cons[4]
+    return a * energy_total(cons, equations)
 end
 
 # Calculate total energy for a conservative state `cons`
 @inline function energy_total(cons, equations::ShallowWaterEquationsQuasi1D)
     a_h, a_h_v, b, a = cons
-    e = (a_h_v^2) / (2 * a_h) + 0.5 * equations.gravity * (a_h^2 / a) + equations.gravity * a_h * b
+    e = (a_h_v^2) / (2 * a * a_h) + 0.5 * equations.gravity * (a_h^2 / a) + equations.gravity * a_h * b
     return e
 end
 
 # Calculate kinetic energy for a conservative state `cons`
 @inline function energy_kinetic(u, equations::ShallowWaterEquationsQuasi1D)
-    a_h, a_h_v, _, _ = u
-    return (a_h_v^2) / (2 * a_h)
+    a_h, a_h_v, b, a = u
+    return (a_h_v^2) / (2 * a * a_h) # 0.5 * v^2
 end
 
-# Calculate potential energy for a conservative state `cons`
-@inline function energy_internal(cons, equations::ShallowWaterEquationsQuasi1D)
-    return energy_total(cons, equations) - energy_kinetic(cons, equations)
-end
 
 # Calculate the error for the "lake-at-rest" test case
 @inline function lake_at_rest_error(u, equations::ShallowWaterEquationsQuasi1D)

From 3362c0441bbd6d8289403053ef0d784adf2c31dd Mon Sep 17 00:00:00 2001
From: Jesse Chan <1156048+jlchan@users.noreply.github.com>
Date: Tue, 29 Aug 2023 17:45:45 -0500
Subject: [PATCH 09/25] Apply suggestions from code review

---
 src/equations/shallow_water_quasi_1d.jl | 20 ++++----------------
 1 file changed, 4 insertions(+), 16 deletions(-)

diff --git a/src/equations/shallow_water_quasi_1d.jl b/src/equations/shallow_water_quasi_1d.jl
index fa23f72dd9d..c1338269323 100644
--- a/src/equations/shallow_water_quasi_1d.jl
+++ b/src/equations/shallow_water_quasi_1d.jl
@@ -44,7 +44,7 @@ This affects the implementation and use of these equations in various ways:
 * Trixi.jl's visualization tools will visualize the bottom topography and channel width by default.
 """
 
-struct ShallowWaterEquationsQuasi1D{RealT <: Real} <: Trixi.AbstractShallowWaterEquations{1, 4}
+struct ShallowWaterEquationsQuasi1D{RealT <: Real} <: AbstractShallowWaterEquations{1, 4}
     gravity::RealT # gravitational constant
     H0::RealT      # constant "lake-at-rest" total water height
     # `threshold_limiter` used in `PositivityPreservingLimiterShallowWater` on water height,
@@ -74,7 +74,7 @@ function ShallowWaterEquationsQuasi1D(; gravity_constant, H0 = zero(gravity_cons
     ShallowWaterEquationsQuasi1D(gravity_constant, H0, threshold_limiter, threshold_wet)
 end
 
-have_nonconservative_terms(::ShallowWaterEquationsQuasi1D) = Trixi.True()
+have_nonconservative_terms(::ShallowWaterEquationsQuasi1D) = True()
 varnames(::typeof(cons2cons), ::ShallowWaterEquationsQuasi1D) = ("a_h", "a_h_v", "b", "a")
 # Note, we use the total water height, H = h + b, as the first primitive variable for easier
 # visualization and setting initial conditions
@@ -288,11 +288,11 @@ end
 #Convert conservative variables to entropy variables
 # Note, only the first two are the entropy variables, the third and foruth entries still
 # just carry the bottom topography and channel width values for convenience
-@inline function Trixi.cons2entropy(u, equations::ShallowWaterEquationsQuasi1D)
+@inline function cons2entropy(u, equations::ShallowWaterEquationsQuasi1D)
     a_h, a_h_v, b, a = u
     h = waterheight(u, equations)
     v = velocity(u, equations)
-    #entropy variables are the same as ones in standard equations
+    #entropy variables are the same as ones in standard shallow water equations
     w1 = equations.gravity * (h + b) - 0.5 * v^2
     w2 = v
     
@@ -326,18 +326,6 @@ end
     return e
 end
 
-# Calculate kinetic energy for a conservative state `cons`
-@inline function energy_kinetic(u, equations::ShallowWaterEquationsQuasi1D)
-    a_h, a_h_v, b, a = u
-    return (a_h_v^2) / (2 * a * a_h) # 0.5 * v^2
-end
 
 
-# Calculate the error for the "lake-at-rest" test case
-@inline function lake_at_rest_error(u, equations::ShallowWaterEquationsQuasi1D)
-    a_h, _, b, _ = u
-    h = waterheight(u, equations)
-    H0_wet_dry = max(equations.H0, b + equations.threshold_limiter)
-    return abs(H0_wet_dry - (h + b))
-end
 end # @muladd
\ No newline at end of file

From 09ec4218ab40b2e0396d76b36ea9ea99dce7750c Mon Sep 17 00:00:00 2001
From: Jesse Chan <1156048+jlchan@users.noreply.github.com>
Date: Tue, 29 Aug 2023 17:47:58 -0500
Subject: [PATCH 10/25] Update src/equations/shallow_water_quasi_1d.jl

---
 src/equations/shallow_water_quasi_1d.jl | 32 -------------------------
 1 file changed, 32 deletions(-)

diff --git a/src/equations/shallow_water_quasi_1d.jl b/src/equations/shallow_water_quasi_1d.jl
index c1338269323..79d5a743ba2 100644
--- a/src/equations/shallow_water_quasi_1d.jl
+++ b/src/equations/shallow_water_quasi_1d.jl
@@ -135,38 +135,6 @@ as defined in [`initial_condition_convergence_test`](@ref).
     return SVector(du1, du2, 0.0, 0.0)
 end  
 
-"""
-    boundary_condition_slip_wall(u_inner, orientation_or_normal, x, t, surface_flux_function,
-                                 equations::ShallowWaterEquationsQuasi1D)
-
-Create a boundary state by reflecting the normal velocity component and keep
-the tangential velocity component unchanged. The boundary water height is taken from
-the internal value.
-
-For details see Section 9.2.5 of the book:
-- Eleuterio F. Toro (2001)
-  Shock-Capturing Methods for Free-Surface Shallow Flows
-  1st edition
-  ISBN 0471987662
-"""
-
-function boundary_condition_slip_wall(u_inner, orientation_or_normal, direction,
-                                              x, t, surface_flux_function,
-                                              equations::ShallowWaterEquationsQuasi1D)
-    # create the "external" boundary solution state
-    u_boundary = SVector(u_inner[1],
-                         -u_inner[2],
-                         u_inner[3],
-                         u_inner[4])
-
-    # calculate the boundary flux
-    if iseven(direction) # u_inner is "left" of boundary, u_boundary is "right" of boundary
-        f = surface_flux_function(u_inner, u_boundary, orientation_or_normal, equations)
-    else # u_boundary is "left" of boundary, u_inner is "right" of boundary
-        f = surface_flux_function(u_boundary, u_inner, orientation_or_normal, equations)
-    end
-    return f
-end
 
 # Calculate 1D flux for a single point
 # Note, the bottom topography and channel width have no flux

From 2d5c7047ebb7b0c6ee14930d0573fcad66d89171 Mon Sep 17 00:00:00 2001
From: Jesse Chan <jesse.lee.chan@gmail.com>
Date: Tue, 29 Aug 2023 17:49:57 -0500
Subject: [PATCH 11/25] formatting

---
 ...xir_shallow_water_quasi_1D_manufactured.jl | 36 +++++------
 src/Trixi.jl                                  |  2 +-
 src/equations/shallow_water_quasi_1d.jl       | 64 ++++++++++---------
 3 files changed, 54 insertions(+), 48 deletions(-)

diff --git a/examples/tree_1d_dgsem/elixir_shallow_water_quasi_1D_manufactured.jl b/examples/tree_1d_dgsem/elixir_shallow_water_quasi_1D_manufactured.jl
index 4f906a79d27..20ccc84c68b 100644
--- a/examples/tree_1d_dgsem/elixir_shallow_water_quasi_1D_manufactured.jl
+++ b/examples/tree_1d_dgsem/elixir_shallow_water_quasi_1D_manufactured.jl
@@ -4,18 +4,18 @@ using Trixi
 ###############################################################################
 # Semidiscretization of the shallow water equations
 
-equations = ShallowWaterEquationsQuasi1D(gravity_constant=9.81)
+equations = ShallowWaterEquationsQuasi1D(gravity_constant = 9.81)
 
 initial_condition = initial_condition_convergence_test
 
-
 ###############################################################################
 # Get the DG approximation space
 
-volume_flux = (flux_chan_etal, flux_nonconservative_chan_etal) 
-surface_flux = (FluxPlusDissipation(flux_chan_etal, DissipationLocalLaxFriedrichs()) , flux_nonconservative_chan_etal)
-solver = DGSEM(polydeg=3, surface_flux=surface_flux,
-               volume_integral=VolumeIntegralFluxDifferencing(volume_flux))
+volume_flux = (flux_chan_etal, flux_nonconservative_chan_etal)
+surface_flux = (FluxPlusDissipation(flux_chan_etal, DissipationLocalLaxFriedrichs()),
+                flux_nonconservative_chan_etal)
+solver = DGSEM(polydeg = 3, surface_flux = surface_flux,
+               volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
 
 ###############################################################################
 # Get the TreeMesh and setup a periodic mesh
@@ -23,13 +23,13 @@ solver = DGSEM(polydeg=3, surface_flux=surface_flux,
 coordinates_min = 0.0
 coordinates_max = sqrt(2.0)
 mesh = TreeMesh(coordinates_min, coordinates_max,
-                initial_refinement_level=3,
-                n_cells_max=10_000,
-                periodicity=true)
+                initial_refinement_level = 3,
+                n_cells_max = 10_000,
+                periodicity = true)
 
 # create the semi discretization object
 semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
-                                    source_terms=source_terms_convergence_test)
+                                    source_terms = source_terms_convergence_test)
 
 ###############################################################################
 # ODE solvers, callbacks etc.
@@ -40,13 +40,13 @@ ode = semidiscretize(semi, tspan)
 summary_callback = SummaryCallback()
 
 analysis_interval = 500
-analysis_callback = AnalysisCallback(semi, interval=analysis_interval)
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
 
-alive_callback = AliveCallback(analysis_interval=analysis_interval)
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
 
-save_solution = SaveSolutionCallback(interval=200,
-                                     save_initial_solution=true,
-                                     save_final_solution=true)
+save_solution = SaveSolutionCallback(interval = 200,
+                                     save_initial_solution = true,
+                                     save_final_solution = true)
 
 callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution)
 
@@ -54,6 +54,6 @@ callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, sav
 # run the simulation
 
 # use a Runge-Kutta method with automatic (error based) time step size control
-sol = solve(ode, RDPK3SpFSAL49(); abstol=1.0e-8, reltol=1.0e-8,
-            ode_default_options()..., callback=callbacks);
-summary_callback() # print the timer summary
\ No newline at end of file
+sol = solve(ode, RDPK3SpFSAL49(); abstol = 1.0e-8, reltol = 1.0e-8,
+            ode_default_options()..., callback = callbacks);
+summary_callback() # print the timer summary
diff --git a/src/Trixi.jl b/src/Trixi.jl
index 63c0812cd48..e171d90c924 100644
--- a/src/Trixi.jl
+++ b/src/Trixi.jl
@@ -148,7 +148,7 @@ export AcousticPerturbationEquations2D,
        InviscidBurgersEquation1D,
        LatticeBoltzmannEquations2D, LatticeBoltzmannEquations3D,
        ShallowWaterEquations1D, ShallowWaterEquations2D,
-       ShallowWaterTwoLayerEquations1D, ShallowWaterTwoLayerEquations2D, 
+       ShallowWaterTwoLayerEquations1D, ShallowWaterTwoLayerEquations2D,
        ShallowWaterEquationsQuasi1D,
        LinearizedEulerEquations2D
 
diff --git a/src/equations/shallow_water_quasi_1d.jl b/src/equations/shallow_water_quasi_1d.jl
index 6924bab0766..1ab86cfa449 100644
--- a/src/equations/shallow_water_quasi_1d.jl
+++ b/src/equations/shallow_water_quasi_1d.jl
@@ -4,7 +4,7 @@
 # See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
 @muladd begin
 #! format: noindent
-    
+
 @doc raw"""
     ShallowWaterEquationsQuasi1D(; gravity, H0 = 0, threshold_limiter = nothing threshold_wet = nothing)
 
@@ -44,7 +44,8 @@ This affects the implementation and use of these equations in various ways:
 * Trixi.jl's visualization tools will visualize the bottom topography and channel width by default.
 """
 
-struct ShallowWaterEquationsQuasi1D{RealT <: Real} <: Trixi.AbstractShallowWaterEquations{1, 4}
+struct ShallowWaterEquationsQuasi1D{RealT <: Real} <:
+       Trixi.AbstractShallowWaterEquations{1, 4}
     gravity::RealT # gravitational constant
     H0::RealT      # constant "lake-at-rest" total water height
     # `threshold_limiter` used in `PositivityPreservingLimiterShallowWater` on water height,
@@ -63,7 +64,8 @@ end
 # well-balancedness test cases.
 # Strict default values for thresholds that performed well in many numerical experiments
 function ShallowWaterEquationsQuasi1D(; gravity_constant, H0 = zero(gravity_constant),
-                                 threshold_limiter = nothing, threshold_wet = nothing)
+                                      threshold_limiter = nothing,
+                                      threshold_wet = nothing)
     T = promote_type(typeof(gravity_constant), typeof(H0))
     if threshold_limiter === nothing
         threshold_limiter = 500 * eps(T)
@@ -75,7 +77,9 @@ function ShallowWaterEquationsQuasi1D(; gravity_constant, H0 = zero(gravity_cons
 end
 
 have_nonconservative_terms(::ShallowWaterEquationsQuasi1D) = Trixi.True()
-varnames(::typeof(cons2cons), ::ShallowWaterEquationsQuasi1D) = ("a_h", "a_h_v", "b", "a")
+function varnames(::typeof(cons2cons), ::ShallowWaterEquationsQuasi1D)
+    ("a_h", "a_h_v", "b", "a")
+end
 # Note, we use the total water height, H = h + b, as the first primitive variable for easier
 # visualization and setting initial conditions
 varnames(::typeof(cons2prim), ::ShallowWaterEquationsQuasi1D) = ("H", "v", "b", "a")
@@ -89,14 +93,15 @@ 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_convergence_test(x, t, equations::ShallowWaterEquationsQuasi1D)
+function initial_condition_convergence_test(x, t,
+                                            equations::ShallowWaterEquationsQuasi1D)
     # generates a manufactured solution. 
     # some constants are chosen such that the function is periodic on the domain [0,sqrt(2)]
-    Omega = sqrt(2) * pi 
+    Omega = sqrt(2) * pi
     H = 2.0 + 0.5 * sin(Omega * x[1] - t)
     v = 0.25
-    b = 0.2 - 0.05 * sin(1 * sqrt(2) * pi *x[1])
-    a = 1 + 0.1 * cos(sqrt(2) * pi * x[1])   
+    b = 0.2 - 0.05 * sin(1 * sqrt(2) * pi * x[1])
+    a = 1 + 0.1 * cos(sqrt(2) * pi * x[1])
     return prim2cons(SVector(H, v, b, a), equations)
 end
 
@@ -116,16 +121,16 @@ as defined in [`initial_condition_convergence_test`](@ref).
                                                equations::ShallowWaterEquationsQuasi1D)
     # Same settings as in `initial_condition_convergence_test`. Some derivative simplify because
     # this manufactured solution velocity is taken to be constant
-    Omega = sqrt(2) * pi 
+    Omega = sqrt(2) * pi
     H = 2.0 + 0.5 * sin(Omega * x[1] - t)
     H_x = 0.5 * cos(Omega * x[1] - t) * Omega
     H_t = -0.5 * cos(Omega * x[1] - t)
 
     v = 0.25
-    
-    b =  0.2 - 0.05 * sin(sqrt(2) * pi *x[1])
-    b_x = -0.05 * cos(sqrt(2) * pi *x[1]) * sqrt(2) * pi
-    
+
+    b = 0.2 - 0.05 * sin(sqrt(2) * pi * x[1])
+    b_x = -0.05 * cos(sqrt(2) * pi * x[1]) * sqrt(2) * pi
+
     a = 1 + 0.1 * cos(sqrt(2) * pi * x[1])
     a_x = -0.1 * sin(sqrt(2) * pi * x[1]) * sqrt(2) * pi
 
@@ -133,7 +138,7 @@ as defined in [`initial_condition_convergence_test`](@ref).
     du2 = v * du1 + a * (equations.gravity * (H - b) * H_x)
 
     return SVector(du1, du2, 0.0, 0.0)
-end  
+end
 
 """
     boundary_condition_slip_wall(u_inner, orientation_or_normal, x, t, surface_flux_function,
@@ -151,8 +156,8 @@ For details see Section 9.2.5 of the book:
 """
 
 function boundary_condition_slip_wall(u_inner, orientation_or_normal, direction,
-                                              x, t, surface_flux_function,
-                                              equations::ShallowWaterEquationsQuasi1D)
+                                      x, t, surface_flux_function,
+                                      equations::ShallowWaterEquationsQuasi1D)
     # create the "external" boundary solution state
     u_boundary = SVector(u_inner[1],
                          -u_inner[2],
@@ -176,7 +181,7 @@ end
     v = velocity(u, equations)
 
     p = 0.5 * a * equations.gravity * h^2
-    
+
     f1 = a_h_v
     f2 = a_h_v * v + p
 
@@ -202,12 +207,12 @@ Further details are available in the paper:
                                                 equations::ShallowWaterEquationsQuasi1D)
     a_h_ll, _, b_ll, a_ll = u_ll
     a_h_rr, _, b_rr, a_rr = u_rr
-    
+
     h_ll = waterheight(u_ll, equations)
     h_rr = waterheight(u_rr, equations)
-    
+
     z = zero(eltype(u_ll))
-    
+
     return SVector(z, equations.gravity * a_ll * h_ll * (h_rr + b_rr), z, z)
 end
 
@@ -226,17 +231,17 @@ 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, 
+@inline function flux_chan_etal(u_ll, u_rr, orientation::Integer,
                                 equations::ShallowWaterEquationsQuasi1D)
     a_h_ll, a_h_v_ll, _, _ = u_ll
     a_h_rr, a_h_v_rr, _, _ = u_rr
-    
+
     v_ll = velocity(u_ll, equations)
     v_rr = velocity(u_rr, equations)
 
     f1 = 0.5 * (a_h_v_ll + a_h_v_rr)
     f2 = f1 * 0.5 * (v_ll + v_rr)
-    
+
     return SVector(f1, f2, zero(eltype(u_ll)), zero(eltype(u_ll)))
 end
 
@@ -270,9 +275,9 @@ end
 # Helper function to extract the velocity vector from the conservative variables
 @inline function velocity(u, equations::ShallowWaterEquationsQuasi1D)
     a_h, a_h_v, _, _ = u
-    
+
     v = a_h_v / a_h
-    
+
     return v
 end
 
@@ -295,14 +300,14 @@ end
     #entropy variables are the same as ones in standard equations
     w1 = equations.gravity * (h + b) - 0.5 * v^2
     w2 = v
-    
+
     return SVector(w1, w2, b, a)
 end
 
 # Convert primitive to conservative variables
 @inline function prim2cons(prim, equations::ShallowWaterEquationsQuasi1D)
     H, v, b, a = prim
-    
+
     a_h = a * (H - b)
     a_h_v = a_h * v
     return SVector(a_h, a_h_v, b, a)
@@ -327,7 +332,8 @@ end
 # Calculate total energy for a conservative state `cons`
 @inline function energy_total(cons, equations::ShallowWaterEquationsQuasi1D)
     a_h, a_h_v, b, a = cons
-    e = (a_h_v^2) / (2 * a_h) + 0.5 * equations.gravity * (a_h^2 / a) + equations.gravity * a_h * b
+    e = (a_h_v^2) / (2 * a_h) + 0.5 * equations.gravity * (a_h^2 / a) +
+        equations.gravity * a_h * b
     return e
 end
 
@@ -349,4 +355,4 @@ end
     H0_wet_dry = max(equations.H0, b + equations.threshold_limiter)
     return abs(H0_wet_dry - (h + b))
 end
-end # @muladd
\ No newline at end of file
+end # @muladd

From d4f58b8897819edeeb8e1e06a60a8f3dd242be55 Mon Sep 17 00:00:00 2001
From: Jesse Chan <jesse.lee.chan@gmail.com>
Date: Tue, 29 Aug 2023 17:53:53 -0500
Subject: [PATCH 12/25] formatting

---
 src/equations/shallow_water_quasi_1d.jl | 15 +++++++--------
 1 file changed, 7 insertions(+), 8 deletions(-)

diff --git a/src/equations/shallow_water_quasi_1d.jl b/src/equations/shallow_water_quasi_1d.jl
index 0989e6e709d..4a07d5e96d1 100644
--- a/src/equations/shallow_water_quasi_1d.jl
+++ b/src/equations/shallow_water_quasi_1d.jl
@@ -44,7 +44,8 @@ This affects the implementation and use of these equations in various ways:
 * Trixi.jl's visualization tools will visualize the bottom topography and channel width by default.
 """
 
-struct ShallowWaterEquationsQuasi1D{RealT <: Real} <: AbstractShallowWaterEquations{1, 4}
+struct ShallowWaterEquationsQuasi1D{RealT <: Real} <:
+       AbstractShallowWaterEquations{1, 4}
     gravity::RealT # gravitational constant
     H0::RealT      # constant "lake-at-rest" total water height
     # `threshold_limiter` used in `PositivityPreservingLimiterShallowWater` on water height,
@@ -76,7 +77,9 @@ function ShallowWaterEquationsQuasi1D(; gravity_constant, H0 = zero(gravity_cons
 end
 
 have_nonconservative_terms(::ShallowWaterEquationsQuasi1D) = True()
-varnames(::typeof(cons2cons), ::ShallowWaterEquationsQuasi1D) = ("a_h", "a_h_v", "b", "a")
+function varnames(::typeof(cons2cons), ::ShallowWaterEquationsQuasi1D)
+    ("a_h", "a_h_v", "b", "a")
+end
 # Note, we use the total water height, H = h + b, as the first primitive variable for easier
 # visualization and setting initial conditions
 varnames(::typeof(cons2prim), ::ShallowWaterEquationsQuasi1D) = ("H", "v", "b", "a")
@@ -137,7 +140,6 @@ as defined in [`initial_condition_convergence_test`](@ref).
     return SVector(du1, du2, 0.0, 0.0)
 end
 
-
 # Calculate 1D flux for a single point
 # Note, the bottom topography and channel width have no flux
 @inline function flux(u, orientation::Integer, equations::ShallowWaterEquationsQuasi1D)
@@ -282,7 +284,6 @@ end
     return u[1] / u[4]
 end
 
-
 # Entropy function for the shallow water equations is the total energy
 @inline function entropy(cons, equations::ShallowWaterEquationsQuasi1D)
     a = cons[4]
@@ -292,10 +293,8 @@ end
 # Calculate total energy for a conservative state `cons`
 @inline function energy_total(cons, equations::ShallowWaterEquationsQuasi1D)
     a_h, a_h_v, b, a = cons
-    e = (a_h_v^2) / (2 * a * a_h) + 0.5 * equations.gravity * (a_h^2 / a) + equations.gravity * a_h * b
+    e = (a_h_v^2) / (2 * a * a_h) + 0.5 * equations.gravity * (a_h^2 / a) +
+        equations.gravity * a_h * b
     return e
 end
-
-
-
 end # @muladd

From d434e897bff362a75ac422639ab121b972fca355 Mon Sep 17 00:00:00 2001
From: Jesse Chan <jesse.lee.chan@gmail.com>
Date: Tue, 29 Aug 2023 19:18:32 -0500
Subject: [PATCH 13/25] forgot comma

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

diff --git a/test/test_tree_1d_shallowwater.jl b/test/test_tree_1d_shallowwater.jl
index b6332181af9..94212f0ba39 100644
--- a/test/test_tree_1d_shallowwater.jl
+++ b/test/test_tree_1d_shallowwater.jl
@@ -115,7 +115,7 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem")
   @trixi_testset "elixir_shallow_water_quasi_1D_manufactured.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallow_water_quasi_1D_manufactured.jl"),
       l2 = [6.37048760275098e-5, 0.0002745658116815704, 4.436491725647962e-6, 8.872983451152218e-6], 
-      linf = [0.00026747526881631956, 0.0012106730729152249, 9.098379777500165e-6, 1.8196759554278685e-5]
+      linf = [0.00026747526881631956, 0.0012106730729152249, 9.098379777500165e-6, 1.8196759554278685e-5],
       tspan = (0.0, 0.05))
   end  
 end

From 2a2cb4644eae65115dfa231af48c263f3a7eafb6 Mon Sep 17 00:00:00 2001
From: Krissh Chawla <127906314+KrisshChawla@users.noreply.github.com>
Date: Wed, 30 Aug 2023 08:31:41 -0500
Subject: [PATCH 14/25] Apply suggestions from code review

Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com>
---
 src/equations/shallow_water_quasi_1d.jl | 10 +++-------
 1 file changed, 3 insertions(+), 7 deletions(-)

diff --git a/src/equations/shallow_water_quasi_1d.jl b/src/equations/shallow_water_quasi_1d.jl
index 4a07d5e96d1..949bac15535 100644
--- a/src/equations/shallow_water_quasi_1d.jl
+++ b/src/equations/shallow_water_quasi_1d.jl
@@ -43,7 +43,6 @@ This affects the implementation and use of these equations in various ways:
 * [`AnalysisCallback`](@ref) analyzes this variable.
 * Trixi.jl's visualization tools will visualize the bottom topography and channel width by default.
 """
-
 struct ShallowWaterEquationsQuasi1D{RealT <: Real} <:
        AbstractShallowWaterEquations{1, 4}
     gravity::RealT # gravitational constant
@@ -116,7 +115,6 @@ This manufactured solution source term is specifically designed for the bottom t
 `b(x) = 0.2 - 0.05 * sin(sqrt(2) * pi *x[1])` and channel width 'a(x)= 1 + 0.1 * cos(sqrt(2) * pi * x[1])'
 as defined in [`initial_condition_convergence_test`](@ref).
 """
-
 @inline function source_terms_convergence_test(u, x, t,
                                                equations::ShallowWaterEquationsQuasi1D)
     # Same settings as in `initial_condition_convergence_test`. Some derivative simplify because
@@ -156,8 +154,8 @@ end
 end
 
 """
-flux_nonconservative_chan_etal(u_ll, u_rr, orientation::Integer,
-                                                equations::ShallowWaterEquationsQuasi1D)
+    flux_nonconservative_chan_etal(u_ll, u_rr, orientation::Integer,
+                                   equations::ShallowWaterEquationsQuasi1D)
 
 Non-symmetric two-point volume flux discretizing the nonconservative (source) term
 that contains the gradient of the bottom topography [`ShallowWaterEquationsQuasi1D`](@ref) 
@@ -169,7 +167,6 @@ Further details are available in the paper:
     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::ShallowWaterEquationsQuasi1D)
     a_h_ll, _, b_ll, a_ll = u_ll
@@ -185,7 +182,7 @@ end
 
 """
     flux_chan_etal(u_ll, u_rr, orientation,
-                          equations::ShallowWaterEquationsQuasi1D)
+                   equations::ShallowWaterEquationsQuasi1D)
 
 Total energy conservative (mathematical entropy for quasi shallow water equations) split form.
 When the bottom topography is nonzero this scheme will be well-balanced when used as a `volume_flux`.
@@ -197,7 +194,6 @@ Further details are available in the paper:
   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::ShallowWaterEquationsQuasi1D)
     a_h_ll, a_h_v_ll, _, _ = u_ll

From 32d028966d5e9b50b227c4fff76328d61809ec0f Mon Sep 17 00:00:00 2001
From: Krissh Chawla <127906314+KrisshChawla@users.noreply.github.com>
Date: Wed, 30 Aug 2023 18:21:58 -0500
Subject: [PATCH 15/25] renamed example elixir to
 elixir_shallow_water_quasi_1d_source_terms.jl

---
 ...ufactured.jl => elixir_shallow_water_quasi_1d_source_terms.jl} | 0
 1 file changed, 0 insertions(+), 0 deletions(-)
 rename examples/tree_1d_dgsem/{elixir_shallow_water_quasi_1D_manufactured.jl => elixir_shallow_water_quasi_1d_source_terms.jl} (100%)

diff --git a/examples/tree_1d_dgsem/elixir_shallow_water_quasi_1D_manufactured.jl b/examples/tree_1d_dgsem/elixir_shallow_water_quasi_1d_source_terms.jl
similarity index 100%
rename from examples/tree_1d_dgsem/elixir_shallow_water_quasi_1D_manufactured.jl
rename to examples/tree_1d_dgsem/elixir_shallow_water_quasi_1d_source_terms.jl

From 0bc2e2927cb889ffa7067097123ccd698b7c06aa Mon Sep 17 00:00:00 2001
From: Krissh Chawla <127906314+KrisshChawla@users.noreply.github.com>
Date: Wed, 30 Aug 2023 18:43:12 -0500
Subject: [PATCH 16/25] Apply suggestions from code review

Co-authored-by: Andrew Winters <andrew.ross.winters@liu.se>
---
 src/equations/shallow_water_quasi_1d.jl | 19 ++++++++++---------
 1 file changed, 10 insertions(+), 9 deletions(-)

diff --git a/src/equations/shallow_water_quasi_1d.jl b/src/equations/shallow_water_quasi_1d.jl
index 949bac15535..e17c350a6db 100644
--- a/src/equations/shallow_water_quasi_1d.jl
+++ b/src/equations/shallow_water_quasi_1d.jl
@@ -99,8 +99,8 @@ function initial_condition_convergence_test(x, t,
     Omega = sqrt(2) * pi
     H = 2.0 + 0.5 * sin(Omega * x[1] - t)
     v = 0.25
-    b = 0.2 - 0.05 * sin(1 * sqrt(2) * pi * x[1])
-    a = 1 + 0.1 * cos(sqrt(2) * pi * x[1])
+    b = 0.2 - 0.05 * sin(Omega * x[1])
+    a = 1 + 0.1 * cos(Omega * x[1])
     return prim2cons(SVector(H, v, b, a), equations)
 end
 
@@ -126,11 +126,11 @@ as defined in [`initial_condition_convergence_test`](@ref).
 
     v = 0.25
 
-    b = 0.2 - 0.05 * sin(sqrt(2) * pi * x[1])
-    b_x = -0.05 * cos(sqrt(2) * pi * x[1]) * sqrt(2) * pi
+    b = 0.2 - 0.05 * sin(Omega * x[1])
+    b_x = -0.05 * cos(Omega * x[1]) * Omega
 
-    a = 1 + 0.1 * cos(sqrt(2) * pi * x[1])
-    a_x = -0.1 * sin(sqrt(2) * pi * x[1]) * sqrt(2) * pi
+    a = 1 + 0.1 * cos(Omega * x[1])
+    a_x = -0.1 * sin(Omega * x[1]) * Omega
 
     du1 = a * H_t + v * (a_x * (H - b) + a * (H_x - b_x))
     du2 = v * du1 + a * (equations.gravity * (H - b) * H_x)
@@ -184,7 +184,7 @@ end
     flux_chan_etal(u_ll, u_rr, orientation,
                    equations::ShallowWaterEquationsQuasi1D)
 
-Total energy conservative (mathematical entropy for quasi shallow water equations) split form.
+Total energy conservative (mathematical entropy for quasi 1D shallow water equations) split form.
 When the bottom topography is nonzero this scheme will be well-balanced when used as a `volume_flux`.
 The `surface_flux` should still use, e.g., [`FluxPlusDissipation(flux_chan_etal, DissipationLocalLaxFriedrichs())`](@ref).
 
@@ -226,6 +226,7 @@ end
 end
 
 # Specialized `DissipationLocalLaxFriedrichs` to avoid spurious dissipation in the bottom topography
+# and channel width
 @inline function (dissipation::DissipationLocalLaxFriedrichs)(u_ll, u_rr,
                                                               orientation_or_normal_direction,
                                                               equations::ShallowWaterEquationsQuasi1D)
@@ -253,8 +254,8 @@ end
     return SVector(H, v, b, a)
 end
 
-#Convert conservative variables to entropy variables
-# Note, only the first two are the entropy variables, the third and foruth entries still
+# Convert conservative variables to entropy variables
+# Note, only the first two are the entropy variables, the third and fourth entries still
 # just carry the bottom topography and channel width values for convenience
 @inline function cons2entropy(u, equations::ShallowWaterEquationsQuasi1D)
     a_h, a_h_v, b, a = u

From 5330cc540952cecdc9b5350777a135a8236d7416 Mon Sep 17 00:00:00 2001
From: Krissh Chawla <127906314+KrisshChawla@users.noreply.github.com>
Date: Wed, 30 Aug 2023 18:45:15 -0500
Subject: [PATCH 17/25] Update test_tree_1d_shallowwater.jl with renamed
 example elixir

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

diff --git a/test/test_tree_1d_shallowwater.jl b/test/test_tree_1d_shallowwater.jl
index 94212f0ba39..3d4e648a7f9 100644
--- a/test/test_tree_1d_shallowwater.jl
+++ b/test/test_tree_1d_shallowwater.jl
@@ -112,7 +112,7 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem")
       tspan = (0.0, 0.05))
   end
 
-  @trixi_testset "elixir_shallow_water_quasi_1D_manufactured.jl" begin
+  @trixi_testset "elixir_shallow_water_quasi_1d_source_terms.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallow_water_quasi_1D_manufactured.jl"),
       l2 = [6.37048760275098e-5, 0.0002745658116815704, 4.436491725647962e-6, 8.872983451152218e-6], 
       linf = [0.00026747526881631956, 0.0012106730729152249, 9.098379777500165e-6, 1.8196759554278685e-5],

From 875b835e3e8d8437914492db1f1dbdc5a666f1af Mon Sep 17 00:00:00 2001
From: Krissh Chawla <127906314+KrisshChawla@users.noreply.github.com>
Date: Thu, 31 Aug 2023 08:29:23 -0500
Subject: [PATCH 18/25] comment fix

---
 src/equations/shallow_water_quasi_1d.jl | 1 -
 1 file changed, 1 deletion(-)

diff --git a/src/equations/shallow_water_quasi_1d.jl b/src/equations/shallow_water_quasi_1d.jl
index e17c350a6db..4893d78764e 100644
--- a/src/equations/shallow_water_quasi_1d.jl
+++ b/src/equations/shallow_water_quasi_1d.jl
@@ -84,7 +84,6 @@ end
 varnames(::typeof(cons2prim), ::ShallowWaterEquationsQuasi1D) = ("H", "v", "b", "a")
 
 # Set initial conditions at physical location `x` for time `t`
-#need to revise
 """
     initial_condition_convergence_test(x, t, equations::ShallowWaterEquationsQuasi1D)
 

From c2f6c36ca6cb7fa927fcdf089a6619a6c2895445 Mon Sep 17 00:00:00 2001
From: Krissh Chawla <127906314+KrisshChawla@users.noreply.github.com>
Date: Thu, 31 Aug 2023 08:50:18 -0500
Subject: [PATCH 19/25] comment fix for
 elixir_shallow_water_quasi_1d_source_terms.jl

---
 .../elixir_shallow_water_quasi_1d_source_terms.jl              | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/examples/tree_1d_dgsem/elixir_shallow_water_quasi_1d_source_terms.jl b/examples/tree_1d_dgsem/elixir_shallow_water_quasi_1d_source_terms.jl
index 20ccc84c68b..72747c669e2 100644
--- a/examples/tree_1d_dgsem/elixir_shallow_water_quasi_1d_source_terms.jl
+++ b/examples/tree_1d_dgsem/elixir_shallow_water_quasi_1d_source_terms.jl
@@ -2,7 +2,8 @@ using OrdinaryDiffEq
 using Trixi
 
 ###############################################################################
-# Semidiscretization of the shallow water equations
+# Semidiscretization of the quasi 1d shallow water equations
+# See Chan et al.  https://doi.org/10.48550/arXiv.2307.12089 for details
 
 equations = ShallowWaterEquationsQuasi1D(gravity_constant = 9.81)
 

From 3df925513dbad9b00cb5484a8886ea16f5298f83 Mon Sep 17 00:00:00 2001
From: KrisshChawla <krisshchawla@gmail.com>
Date: Thu, 7 Sep 2023 19:34:05 -0500
Subject: [PATCH 20/25] Added well-balancedness test for shallow_water_quasi_1d

The initial condition in the elixir is intended to test a discontinuous channel width 'a(x)' and bottom topography 'b(x)' on a periodic mesh.
---
 ...ixir_shallowwater_quasi1d_well_balanced.jl | 81 +++++++++++++++++++
 1 file changed, 81 insertions(+)
 create mode 100644 examples/tree_1d_dgsem/elixir_shallowwater_quasi1d_well_balanced.jl

diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_quasi1d_well_balanced.jl b/examples/tree_1d_dgsem/elixir_shallowwater_quasi1d_well_balanced.jl
new file mode 100644
index 00000000000..94ddfa73b26
--- /dev/null
+++ b/examples/tree_1d_dgsem/elixir_shallowwater_quasi1d_well_balanced.jl
@@ -0,0 +1,81 @@
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the shallow water equations with a discontinuous
+# bottom topography function and channel width function
+
+equations = ShallowWaterEquationsQuasi1D(gravity_constant = 9.81, H0 = 2.0)
+
+# Setup a truly discontinuous bottom topography function and channel width for
+# this academic testcase of well-balancedness. The errors from the analysis
+# callback are not important but the error for this lake-at-rest test case
+# `∑|H0-(h+b)|` should be around machine roundoff.
+# Works as intended for TreeMesh1D with `initial_refinement_level=3`. If the mesh
+# refinement level is changed the initial condition below may need changed as well to
+# ensure that the discontinuities lie on an element interface.
+function initial_condition_discontinuous_well_balancedness(x, t,
+                                                           equations::ShallowWaterEquationsQuasi1D)
+    H = equations.H0
+    v = 0.0
+    b = 0.5 * (x[1] + 1)
+    a = 2 + x[1]
+
+    return prim2cons(SVector(H, v, b, a), equations)
+end
+
+initial_condition = initial_condition_discontinuous_well_balancedness
+
+###############################################################################
+# Get the DG approximation space
+
+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))
+
+###############################################################################
+# Get the TreeMesh and setup a periodic mesh
+
+coordinates_min = -1.0
+coordinates_max = 1.0
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 3,
+                n_cells_max = 10_000)
+
+# Create the semi discretization object
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+###############################################################################
+# ODE solver
+
+tspan = (0.0, 100.0)
+ode = semidiscretize(semi, tspan)
+
+###############################################################################
+# Callbacks
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 1000
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
+                                     extra_analysis_integrals = (lake_at_rest_error,))
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(interval = 1000,
+                                     save_initial_solution = true,
+                                     save_final_solution = true)
+
+stepsize_callback = StepsizeCallback(cfl = 3.0)
+
+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

From baef12da47760082b08f35104338368a886c7f2a Mon Sep 17 00:00:00 2001
From: KrisshChawla <krisshchawla@gmail.com>
Date: Thu, 7 Sep 2023 19:39:59 -0500
Subject: [PATCH 21/25] Added 'max_abs_speeds' function and
 'lake_at_rest_error'

---
 src/equations/shallow_water_quasi_1d.jl | 27 +++++++++++++++++++++++++
 1 file changed, 27 insertions(+)

diff --git a/src/equations/shallow_water_quasi_1d.jl b/src/equations/shallow_water_quasi_1d.jl
index 4893d78764e..217a764e173 100644
--- a/src/equations/shallow_water_quasi_1d.jl
+++ b/src/equations/shallow_water_quasi_1d.jl
@@ -235,6 +235,14 @@ end
     return SVector(diss[1], diss[2], zero(eltype(u_ll)), zero(eltype(u_ll)))
 end
 
+@inline function max_abs_speeds(u, equations::ShallowWaterEquationsQuasi1D)
+    h = waterheight(u, equations)
+    v = velocity(u, equations)
+
+    c = equations.gravity * sqrt(h)
+    return (abs(v) + c,)
+end
+
 # Helper function to extract the velocity vector from the conservative variables
 @inline function velocity(u, equations::ShallowWaterEquationsQuasi1D)
     a_h, a_h_v, _, _ = u
@@ -293,4 +301,23 @@ end
         equations.gravity * a_h * b
     return e
 end
+
+# Calculate the error for the "lake-at-rest" test case where H = h+b should
+# be a constant value over time. Note, assumes there is a single reference
+# water height `H0` with which to compare.
+#
+# TODO: TrixiShallowWater: where should `threshold_limiter` live? May need
+# to modify or have different versions of the `lake_at_rest_error` function
+@inline function lake_at_rest_error(u, equations::ShallowWaterEquationsQuasi1D)
+    _, _, b, _ = u
+    h = waterheight(u, equations)
+
+    # For well-balancedness testing with possible wet/dry regions the reference
+    # water height `H0` accounts for the possibility that the bottom topography
+    # can emerge out of the water as well as for the threshold offset to avoid
+    # division by a "hard" zero water heights as well.
+    H0_wet_dry = max(equations.H0, b + equations.threshold_limiter)
+
+    return abs(H0_wet_dry - (h + b))
+end
 end # @muladd

From c70ca5ce5e097cf32b891012bbbc93737ef76c91 Mon Sep 17 00:00:00 2001
From: KrisshChawla <krisshchawla@gmail.com>
Date: Thu, 7 Sep 2023 20:13:42 -0500
Subject: [PATCH 22/25] Updated test_tree_1d_shallowwater with quasi
 well-balancedness test

---
 test/test_tree_1d_shallowwater.jl | 9 ++++++++-
 1 file changed, 8 insertions(+), 1 deletion(-)

diff --git a/test/test_tree_1d_shallowwater.jl b/test/test_tree_1d_shallowwater.jl
index 3d4e648a7f9..69d207642f3 100644
--- a/test/test_tree_1d_shallowwater.jl
+++ b/test/test_tree_1d_shallowwater.jl
@@ -117,7 +117,14 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem")
       l2 = [6.37048760275098e-5, 0.0002745658116815704, 4.436491725647962e-6, 8.872983451152218e-6], 
       linf = [0.00026747526881631956, 0.0012106730729152249, 9.098379777500165e-6, 1.8196759554278685e-5],
       tspan = (0.0, 0.05))
-  end  
+  end
+  
+  @trixi_testset "elixir_shallowwater_quasi1d_well_balanced.jl" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_quasi1d_well_balanced.jl"),
+      l2 = [1.4250229186905198e-14, 2.495109919406496e-12, 7.408599286788738e-17, 2.7205812409138776e-16], 
+      linf = [5.284661597215745e-14, 2.74056233065078e-12, 2.220446049250313e-16, 8.881784197001252e-16],
+      tspan = (0.0, 100.0))
+  end
 end
 
 end # module

From edf45e04ef0a7fd710d964b8bcb9899864c14237 Mon Sep 17 00:00:00 2001
From: KrisshChawla <krisshchawla@gmail.com>
Date: Thu, 7 Sep 2023 20:24:54 -0500
Subject: [PATCH 23/25] File name fix in test_tree_1d_shallowwater

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

diff --git a/test/test_tree_1d_shallowwater.jl b/test/test_tree_1d_shallowwater.jl
index 69d207642f3..9bd0265f38e 100644
--- a/test/test_tree_1d_shallowwater.jl
+++ b/test/test_tree_1d_shallowwater.jl
@@ -113,7 +113,7 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem")
   end
 
   @trixi_testset "elixir_shallow_water_quasi_1d_source_terms.jl" begin
-    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallow_water_quasi_1D_manufactured.jl"),
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallow_water_quasi_1d_source_terms.jl"),
       l2 = [6.37048760275098e-5, 0.0002745658116815704, 4.436491725647962e-6, 8.872983451152218e-6], 
       linf = [0.00026747526881631956, 0.0012106730729152249, 9.098379777500165e-6, 1.8196759554278685e-5],
       tspan = (0.0, 0.05))

From cc4c50b562abf83ce5ca42bbd3b92ad8fdd2a240 Mon Sep 17 00:00:00 2001
From: Krissh Chawla <127906314+KrisshChawla@users.noreply.github.com>
Date: Thu, 7 Sep 2023 20:48:09 -0500
Subject: [PATCH 24/25] Update
 examples/tree_1d_dgsem/elixir_shallowwater_quasi1d_well_balanced.jl

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

diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_quasi1d_well_balanced.jl b/examples/tree_1d_dgsem/elixir_shallowwater_quasi1d_well_balanced.jl
index 94ddfa73b26..d9f1a52b500 100644
--- a/examples/tree_1d_dgsem/elixir_shallowwater_quasi1d_well_balanced.jl
+++ b/examples/tree_1d_dgsem/elixir_shallowwater_quasi1d_well_balanced.jl
@@ -18,6 +18,9 @@ function initial_condition_discontinuous_well_balancedness(x, t,
                                                            equations::ShallowWaterEquationsQuasi1D)
     H = equations.H0
     v = 0.0
+
+    # for a periodic domain, this choice of `b` and `a` mimic 
+    # discontinuity across the periodic boundary.
     b = 0.5 * (x[1] + 1)
     a = 2 + x[1]
 

From 7508ebc27f2c706c77549f43667838324846682d Mon Sep 17 00:00:00 2001
From: KrisshChawla <krisshchawla@gmail.com>
Date: Thu, 7 Sep 2023 20:51:32 -0500
Subject: [PATCH 25/25] Renamed to
 "elixir_shallowwater_quasi_1d_well_balanced.jl"

---
 ...anced.jl => elixir_shallowwater_quasi_1d_well_balanced.jl} | 0
 test/test_tree_1d_shallowwater.jl                             | 4 ++--
 2 files changed, 2 insertions(+), 2 deletions(-)
 rename examples/tree_1d_dgsem/{elixir_shallowwater_quasi1d_well_balanced.jl => elixir_shallowwater_quasi_1d_well_balanced.jl} (100%)

diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_quasi1d_well_balanced.jl b/examples/tree_1d_dgsem/elixir_shallowwater_quasi_1d_well_balanced.jl
similarity index 100%
rename from examples/tree_1d_dgsem/elixir_shallowwater_quasi1d_well_balanced.jl
rename to examples/tree_1d_dgsem/elixir_shallowwater_quasi_1d_well_balanced.jl
diff --git a/test/test_tree_1d_shallowwater.jl b/test/test_tree_1d_shallowwater.jl
index 6a99ae780e5..09fb2d9e432 100644
--- a/test/test_tree_1d_shallowwater.jl
+++ b/test/test_tree_1d_shallowwater.jl
@@ -120,8 +120,8 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem")
       tspan = (0.0, 0.05))
   end
   
-  @trixi_testset "elixir_shallowwater_quasi1d_well_balanced.jl" begin
-    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_quasi1d_well_balanced.jl"),
+  @trixi_testset "elixir_shallowwater_quasi_1d_well_balanced.jl" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_quasi_1d_well_balanced.jl"),
       l2 = [1.4250229186905198e-14, 2.495109919406496e-12, 7.408599286788738e-17, 2.7205812409138776e-16], 
       linf = [5.284661597215745e-14, 2.74056233065078e-12, 2.220446049250313e-16, 8.881784197001252e-16],
       tspan = (0.0, 100.0))