Fix flux for quasi-1D SWE (#1731)

* fix flux (remove non-conservative part)

* add test

* formatting

---------

Co-authored-by: Hendrik Ranocha <[email protected]>

examples/tree_1d_dgsem/elixir_shallowwater_quasi_1d_discontinuous.jl

@@ -0,0 +1,73 @@
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# 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)
+
+function initial_condition_discontinuity(x, t, equations::ShallowWaterEquationsQuasi1D)
+    H = 2 + 0.1 * exp(-25 * x[1]^2)
+    v = 0.0
+
+    if x[1] > 0
+        b = 0.1
+        a = 1.0
+    else
+        b = 0.0
+        a = 1.1
+    end
+
+    return prim2cons(SVector(H, v, b, a), equations)
+end
+
+initial_condition = initial_condition_discontinuity
+
+###############################################################################
+# Get the DG approximation space
+
+volume_flux = (flux_chan_etal, flux_nonconservative_chan_etal)
+surface_flux = (flux_lax_friedrichs, 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.5
+coordinates_max = 0.5
+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)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 0.25)
+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

src/equations/shallow_water_quasi_1d.jl

@@ -137,17 +137,14 @@ as defined in [`initial_condition_convergence_test`](@ref).
     return SVector(du1, du2, 0.0, 0.0)
 end
 
-# Calculate 1D flux for a single point
+# Calculate 1D conservative flux for a single point
 # Note, the bottom topography and channel width have no flux
 @inline function flux(u, orientation::Integer, equations::ShallowWaterEquationsQuasi1D)
-    a_h, a_h_v, _, a = u
-    h = waterheight(u, equations)
+    _, a_h_v, _, _ = u
     v = velocity(u, equations)
 
-    p = 0.5 * a * equations.gravity * h^2
-
     f1 = a_h_v
-    f2 = a_h_v * v + p
+    f2 = a_h_v * v
 
     return SVector(f1, f2, zero(eltype(u)), zero(eltype(u)))
 end

test/test_tree_1d_shallowwater.jl

@@ -387,6 +387,31 @@ end
         @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
     end
 end
+
+@trixi_testset "elixir_shallowwater_quasi_1d_discontinuous.jl" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR,
+                                 "elixir_shallowwater_quasi_1d_discontinuous.jl"),
+                        l2=[
+                            0.02843233740533314,
+                            0.14083324483705398,
+                            0.0054554472558998,
+                            0.005455447255899814,
+                        ],
+                        linf=[
+                            0.26095842440037487,
+                            0.45919004549253795,
+                            0.09999999999999983,
+                            0.10000000000000009,
+                        ],)
+    # 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