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/examples/tree_2d_dgsem/elixir_shallowwater_wall.jl b/examples/tree_2d_dgsem/elixir_shallowwater_wall.jl new file mode 100644 index 00000000000..b8dbad50680 --- /dev/null +++ b/examples/tree_2d_dgsem/elixir_shallowwater_wall.jl @@ -0,0 +1,82 @@ +using OrdinaryDiffEq +using Trixi + +############################################################################### +# Semidiscretization of the shallow water equations + +equations = ShallowWaterEquations2D(gravity_constant = 9.81, H0 = 3.25) + +# An initial condition with a bottom topography and a perturbation in the waterheight to test +# boundary_condition_slip_wall +function initial_condition_perturbation(x, t, equations::ShallowWaterEquations2D) + # Set the background values + H = equations.H0 + v1 = 0.0 + v2 = 0.0 + + # Bottom topography + b = 1.5 * exp(-0.5 * ((x[1])^2 + (x[2])^2)) + # Waterheight perturbation + H = H + 0.5 * exp(-10.0 * ((x[1])^2 + (x[2])^2)) + + return prim2cons(SVector(H, v1, v2, b), equations) +end + +initial_condition = initial_condition_perturbation + +boundary_condition = boundary_condition_slip_wall + +############################################################################### +# Get the DG approximation space + +volume_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal) +surface_flux = (flux_lax_friedrichs, flux_nonconservative_ersing_etal) +solver = DGSEM(polydeg = 3, surface_flux = surface_flux, + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +############################################################################### +# Get the TreeMesh and setup a non-periodic mesh + +coordinates_min = (-1.0, -1.0) +coordinates_max = (1.0, 1.0) +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 4, + n_cells_max = 10_000, + periodicity = false) + +# create the semi discretization object +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_condition) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 0.25) +ode = semidiscretize(semi, tspan) + +############################################################################### +# Callbacks + +summary_callback = SummaryCallback() + +analysis_interval = 1000 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 1000, + save_initial_solution = true, + save_final_solution = true) + +stepsize_callback = StepsizeCallback(cfl = 1.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 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_2d.jl b/src/equations/shallow_water_2d.jl index e75c92a27d0..6728d7d5553 100644 --- a/src/equations/shallow_water_2d.jl +++ b/src/equations/shallow_water_2d.jl @@ -236,8 +236,12 @@ Should be used together with [`TreeMesh`](@ref). u_boundary = SVector(u_inner[1], u_inner[2], -u_inner[3], u_inner[4]) end - # compute and return the flux using `boundary_condition_slip_wall` routine above - flux = surface_flux_function(u_inner, u_boundary, orientation, equations) + # Calculate boundary flux + if iseven(direction) # u_inner is "left" of boundary, u_boundary is "right" of boundary + flux = surface_flux_function(u_inner, u_boundary, orientation, equations) + else # u_boundary is "left" of boundary, u_inner is "right" of boundary + flux = surface_flux_function(u_boundary, u_inner, orientation, equations) + end return flux end 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 diff --git a/test/test_tree_2d_shallowwater.jl b/test/test_tree_2d_shallowwater.jl index 58db7c5f35f..1f3dfbf5267 100644 --- a/test/test_tree_2d_shallowwater.jl +++ b/test/test_tree_2d_shallowwater.jl @@ -327,6 +327,31 @@ end @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 end end + +@trixi_testset "elixir_shallowwater_wall.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_wall.jl"), + l2=[ + 0.13517233723296504, + 0.20010876311162215, + 0.20010876311162223, + 2.719538414346464e-7, + ], + linf=[ + 0.5303607982988336, + 0.5080989745682338, + 0.5080989745682352, + 1.1301675764130437e-6, + ], + tspan=(0.0, 0.25)) + # 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