From e58215857c2b84f3418896840944c144bf832bb6 Mon Sep 17 00:00:00 2001 From: Jesse Chan <1156048+jlchan@users.noreply.github.com> Date: Thu, 4 Jan 2024 13:00:05 -0600 Subject: [PATCH] Extend `CompressibleEulerQuasi1D` and `ShallowWaterQuasi1D` to `DGMulti` (#1797) * adding DGMulti versions of fluxes * remove incorrect factor of 2 * add example and test * formatting * add comment * revert removing factor of 2 * formatting * add SWE quasi-1D test d * enable quasi1D SWE for DGMulti * add docstrings * formatting * Update src/equations/compressible_euler_quasi_1d.jl Co-authored-by: Hendrik Ranocha * adding comments explaining why `normal_direction` is included in 1D * Apply suggestions from code review Co-authored-by: Daniel Doehring --------- Co-authored-by: Hendrik Ranocha Co-authored-by: Daniel Doehring --- examples/dgmulti_1d/elixir_euler_quasi_1d.jl | 49 +++++++++++++++++ .../elixir_shallow_water_quasi_1d.jl | 49 +++++++++++++++++ src/equations/compressible_euler_1d.jl | 5 ++ src/equations/compressible_euler_quasi_1d.jl | 38 ++++++++++++- src/equations/shallow_water_quasi_1d.jl | 35 ++++++++++++ test/test_dgmulti_1d.jl | 53 +++++++++++++++++++ 6 files changed, 227 insertions(+), 2 deletions(-) create mode 100644 examples/dgmulti_1d/elixir_euler_quasi_1d.jl create mode 100644 examples/dgmulti_1d/elixir_shallow_water_quasi_1d.jl diff --git a/examples/dgmulti_1d/elixir_euler_quasi_1d.jl b/examples/dgmulti_1d/elixir_euler_quasi_1d.jl new file mode 100644 index 00000000000..19269fa925b --- /dev/null +++ b/examples/dgmulti_1d/elixir_euler_quasi_1d.jl @@ -0,0 +1,49 @@ +using OrdinaryDiffEq +using Trixi + +############################################################################### +# Semidiscretization of the quasi 1d compressible Euler equations +# See Chan et al. https://doi.org/10.48550/arXiv.2307.12089 for details + +equations = CompressibleEulerEquationsQuasi1D(1.4) + +initial_condition = 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 diff --git a/examples/dgmulti_1d/elixir_shallow_water_quasi_1d.jl b/examples/dgmulti_1d/elixir_shallow_water_quasi_1d.jl new file mode 100644 index 00000000000..85741f9dbd3 --- /dev/null +++ b/examples/dgmulti_1d/elixir_shallow_water_quasi_1d.jl @@ -0,0 +1,49 @@ +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) + +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 diff --git a/src/equations/compressible_euler_1d.jl b/src/equations/compressible_euler_1d.jl index 05c38ce791d..a50c896cd1c 100644 --- a/src/equations/compressible_euler_1d.jl +++ b/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) diff --git a/src/equations/compressible_euler_quasi_1d.jl b/src/equations/compressible_euler_quasi_1d.jl index 0a543277ee4..6844bf9bee5 100644 --- a/src/equations/compressible_euler_quasi_1d.jl +++ b/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, diff --git a/src/equations/shallow_water_quasi_1d.jl b/src/equations/shallow_water_quasi_1d.jl index d3935f0e75f..d52fbab841d 100644 --- a/src/equations/shallow_water_quasi_1d.jl +++ b/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, diff --git a/test/test_dgmulti_1d.jl b/test/test_dgmulti_1d.jl index 79ad64075b4..0363086341f 100644 --- a/test/test_dgmulti_1d.jl +++ b/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