Extend CompressibleEulerQuasi1D and ShallowWaterQuasi1D to DGMulti

examples/dgmulti_1d/elixir_euler_quasi_1d.jl

@@ -0,0 +1,50 @@
+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_convergence_test
+
+surface_flux = (flux_chan_etal, flux_nonconservative_chan_etal)
+volume_flux = surface_flux
+dg = DGMulti(polydeg = 4, element_type = Line(), approximation_type = SBP(),
+             surface_integral = SurfaceIntegralWeakForm(surface_flux),
+             volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+
+cells_per_dimension = (8,)
+mesh = DGMultiMesh(dg, cells_per_dimension,
+                   coordinates_min = (-1.0,), coordinates_max = (1.0,), periodicity = true)
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, dg;
+                                    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, uEltype = real(dg))
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+stepsize_callback = StepsizeCallback(cfl = 0.8)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback,
+                        alive_callback,
+                        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

examples/dgmulti_1d/elixir_shallow_water_quasi_1d.jl

@@ -0,0 +1,50 @@
+using OrdinaryDiffEq
+using Trixi
+using ForwardDiff
+
+###############################################################################
+# 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)
+
+initial_condition = initial_condition_convergence_test
+
+volume_flux = (flux_chan_etal, flux_nonconservative_chan_etal)
+surface_flux = (FluxPlusDissipation(flux_chan_etal, DissipationLocalLaxFriedrichs()),
+                flux_nonconservative_chan_etal)
+
+dg = DGMulti(polydeg = 4, element_type = Line(), approximation_type = SBP(),
+             surface_integral = SurfaceIntegralWeakForm(surface_flux),
+             volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+
+cells_per_dimension = (8,)
+mesh = DGMultiMesh(dg, cells_per_dimension,
+                   coordinates_min = (0.0,), coordinates_max = (sqrt(2),),
+                   periodicity = true)
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, dg;
+                                    source_terms = source_terms_convergence_test)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 1.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval, uEltype = real(dg))
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback,
+                        alive_callback)
+
+###############################################################################
+# run the simulation
+
+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/compressible_euler_1d.jl

@@ -408,6 +408,11 @@ See also
     return SVector(f1, f2, f3)
 end
 
+# While `normal_direction` isn't strictly necessary in 1D, certain solvers assume that 
+# the normal component is incorporated into the numerical flux. 
+# 
+# See `flux(u, normal_direction::AbstractVector, equations::AbstractEquations{1})` for a 
+# similar implementation.
 @inline function flux_ranocha(u_ll, u_rr, normal_direction::AbstractVector,
                               equations::CompressibleEulerEquations1D)
     return normal_direction[1] * flux_ranocha(u_ll, u_rr, 1, equations)

src/equations/compressible_euler_quasi_1d.jl

@@ -161,8 +161,12 @@ end
 end
 
 """
-@inline function flux_nonconservative_chan_etal(u_ll, u_rr, orientation::Integer,
-                                                equations::CompressibleEulerEquationsQuasi1D)
+    flux_nonconservative_chan_etal(u_ll, u_rr, orientation::Integer,
+                                   equations::CompressibleEulerEquationsQuasi1D)
+    flux_nonconservative_chan_etal(u_ll, u_rr, normal_direction, 
+                                   equations::CompressibleEulerEquationsQuasi1D)
+    flux_nonconservative_chan_etal(u_ll, u_rr, normal_ll, normal_rr,
+                                   equations::CompressibleEulerEquationsQuasi1D)
 
 Non-symmetric two-point volume flux discretizing the nonconservative (source) term
 that contains the gradient of the pressure  [`CompressibleEulerEquationsQuasi1D`](@ref) 
@@ -190,6 +194,26 @@ Further details are available in the paper:
     return SVector(z, a_ll * p_avg, z, z)
 end
 
+# While `normal_direction` isn't strictly necessary in 1D, certain solvers assume that 
+# the normal component is incorporated into the numerical flux. 
+# 
+# See `flux(u, normal_direction::AbstractVector, equations::AbstractEquations{1})` for a 
+# similar implementation.
+@inline function flux_nonconservative_chan_etal(u_ll, u_rr,
+                                                normal_direction::AbstractVector,
+                                                equations::CompressibleEulerEquationsQuasi1D)
+    return normal_direction[1] *
+           flux_nonconservative_chan_etal(u_ll, u_rr, 1, equations)
+end
+
+@inline function flux_nonconservative_chan_etal(u_ll, u_rr,
+                                                normal_ll::AbstractVector,
+                                                normal_rr::AbstractVector,
+                                                equations::CompressibleEulerEquationsQuasi1D)
+    # normal_ll should be equal to normal_rr in 1D
+    return flux_nonconservative_chan_etal(u_ll, u_rr, normal_ll, equations)
+end
+
 """
 @inline function flux_chan_etal(u_ll, u_rr, orientation::Integer,
                                            equations::CompressibleEulerEquationsQuasi1D)
@@ -230,6 +254,16 @@ Further details are available in the paper:
     return SVector(f1, f2, f3, zero(eltype(u_ll)))
 end
 
+# While `normal_direction` isn't strictly necessary in 1D, certain solvers assume that 
+# the normal component is incorporated into the numerical flux. 
+# 
+# See `flux(u, normal_direction::AbstractVector, equations::AbstractEquations{1})` for a 
+# similar implementation.
+@inline function flux_chan_etal(u_ll, u_rr, normal_direction::AbstractVector,
+                                equations::CompressibleEulerEquationsQuasi1D)
+    return normal_direction[1] * flux_chan_etal(u_ll, u_rr, 1, equations)
+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,

src/equations/shallow_water_quasi_1d.jl

@@ -152,6 +152,11 @@ end
 """
     flux_nonconservative_chan_etal(u_ll, u_rr, orientation::Integer,
                                    equations::ShallowWaterEquationsQuasi1D)
+    flux_nonconservative_chan_etal(u_ll, u_rr, normal_direction::AbstractVector,
+                                   equations::ShallowWaterEquationsQuasi1D)    
+    flux_nonconservative_chan_etal(u_ll, u_rr, 
+                                   normal_ll::AbstractVector, normal_rr::AbstractVector,
+                                   equations::ShallowWaterEquationsQuasi1D)    
 
 Non-symmetric two-point volume flux discretizing the nonconservative (source) term
 that contains the gradient of the bottom topography [`ShallowWaterEquationsQuasi1D`](@ref) 
@@ -176,6 +181,26 @@ Further details are available in the paper:
     return SVector(z, equations.gravity * a_ll * h_ll * (h_rr + b_rr), z, z)
 end
 
+# While `normal_direction` isn't strictly necessary in 1D, certain solvers assume that 
+# the normal component is incorporated into the numerical flux. 
+# 
+# See `flux(u, normal_direction::AbstractVector, equations::AbstractEquations{1})` for a 
+# similar implementation.
+@inline function flux_nonconservative_chan_etal(u_ll, u_rr,
+                                                normal_direction::AbstractVector,
+                                                equations::ShallowWaterEquationsQuasi1D)
+    return normal_direction[1] *
+           flux_nonconservative_chan_etal(u_ll, u_rr, 1, equations)
+end
+
+@inline function flux_nonconservative_chan_etal(u_ll, u_rr,
+                                                normal_ll::AbstractVector,
+                                                normal_rr::AbstractVector,
+                                                equations::ShallowWaterEquationsQuasi1D)
+    # normal_ll should be equal to normal_rr                                                
+    return flux_nonconservative_chan_etal(u_ll, u_rr, normal_ll, equations)
+end
+
 """
     flux_chan_etal(u_ll, u_rr, orientation,
                    equations::ShallowWaterEquationsQuasi1D)
@@ -204,6 +229,16 @@ Further details are available in the paper:
     return SVector(f1, f2, zero(eltype(u_ll)), zero(eltype(u_ll)))
 end
 
+# While `normal_direction` isn't strictly necessary in 1D, certain solvers assume that 
+# the normal component is incorporated into the numerical flux. 
+# 
+# See `flux(u, normal_direction::AbstractVector, equations::AbstractEquations{1})` for a 
+# similar implementation.
+@inline function flux_chan_etal(u_ll, u_rr, normal_direction::AbstractVector,
+                                equations::ShallowWaterEquationsQuasi1D)
+    return normal_direction[1] * flux_chan_etal(u_ll, u_rr, 1, equations)
+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,

test/test_dgmulti_1d.jl

@@ -162,6 +162,59 @@ end
     @test minimum(dg.basis.rst[1]) ≈ -1
     @test maximum(dg.basis.rst[1])≈1 atol=0.35
 end
+
+# test non-conservative systems
+@trixi_testset "elixir_shallow_water_quasi_1d.jl (SBP) " begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallow_water_quasi_1d.jl"),
+                        cells_per_dimension=(8,),
+                        approximation_type=SBP(),
+                        l2=[
+                            3.03001101100507e-6,
+                            1.692177335948727e-5,
+                            3.002634351734614e-16,
+                            1.1636653574178203e-15,
+                        ],
+                        linf=[
+                            1.2043401988570679e-5,
+                            5.346847010329059e-5,
+                            9.43689570931383e-16,
+                            2.220446049250313e-15,
+                        ])
+    # 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.jl (SBP) " begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_quasi_1d.jl"),
+                        cells_per_dimension=(8,),
+                        approximation_type=SBP(),
+                        l2=[
+                            1.633271343738687e-5,
+                            9.575385661756332e-6,
+                            1.2700331443128421e-5,
+                            0.0,
+                        ],
+                        linf=[
+                            7.304984704381567e-5,
+                            5.2365944135601694e-5,
+                            6.469559594934893e-5,
+                            0.0,
+                        ])
+    # 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
 
 # Clean up afterwards: delete Trixi.jl output directory