From 8e89df78daf0c98360fdbaa4be0055b500c2af2d Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Fri, 10 Nov 2023 21:42:31 -0600 Subject: [PATCH 1/3] fix flux (remove non-conservative part) --- src/equations/shallow_water_quasi_1d.jl | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/src/equations/shallow_water_quasi_1d.jl b/src/equations/shallow_water_quasi_1d.jl index 217a764e173..d3935f0e75f 100644 --- a/src/equations/shallow_water_quasi_1d.jl +++ b/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 From 819d28587af463e3f625cf0bd423717ce14406d7 Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Fri, 10 Nov 2023 21:42:34 -0600 Subject: [PATCH 2/3] add test --- ...xir_shallowwater_quasi_1d_discontinuous.jl | 75 +++++++++++++++++++ test/test_tree_1d_shallowwater.jl | 25 +++++++ 2 files changed, 100 insertions(+) create mode 100644 examples/tree_1d_dgsem/elixir_shallowwater_quasi_1d_discontinuous.jl diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_quasi_1d_discontinuous.jl b/examples/tree_1d_dgsem/elixir_shallowwater_quasi_1d_discontinuous.jl new file mode 100644 index 00000000000..191e027360b --- /dev/null +++ b/examples/tree_1d_dgsem/elixir_shallowwater_quasi_1d_discontinuous.jl @@ -0,0 +1,75 @@ +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, .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 + diff --git a/test/test_tree_1d_shallowwater.jl b/test/test_tree_1d_shallowwater.jl index 7ec3089d33a..fd69a0c1d0e 100644 --- a/test/test_tree_1d_shallowwater.jl +++ b/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 From 37a7f1811f33255504590dc52d9694ecaf085bdb Mon Sep 17 00:00:00 2001 From: Jesse Chan Date: Fri, 10 Nov 2023 21:44:45 -0600 Subject: [PATCH 3/3] formatting --- .../elixir_shallowwater_quasi_1d_discontinuous.jl | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_quasi_1d_discontinuous.jl b/examples/tree_1d_dgsem/elixir_shallowwater_quasi_1d_discontinuous.jl index 191e027360b..76c04759389 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_quasi_1d_discontinuous.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_quasi_1d_discontinuous.jl @@ -8,7 +8,6 @@ using Trixi 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 @@ -37,7 +36,7 @@ solver = DGSEM(polydeg = 3, surface_flux = surface_flux, # Get the TreeMesh and setup a periodic mesh coordinates_min = -0.5 -coordinates_max = 0.5 +coordinates_max = 0.5 mesh = TreeMesh(coordinates_min, coordinates_max, initial_refinement_level = 3, n_cells_max = 10_000, @@ -49,7 +48,7 @@ semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) ############################################################################### # ODE solvers, callbacks etc. -tspan = (0.0, .25) +tspan = (0.0, 0.25) ode = semidiscretize(semi, tspan) summary_callback = SummaryCallback() @@ -72,4 +71,3 @@ callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, sav sol = solve(ode, RDPK3SpFSAL49(); abstol = 1.0e-8, reltol = 1.0e-8, ode_default_options()..., callback = callbacks); summary_callback() # print the timer summary -