From 4644c7e7b9e6bd15c73f1dd131b39d03c6e48d4f Mon Sep 17 00:00:00 2001 From: Patrick Ersing <114223904+patrickersing@users.noreply.github.com> Date: Wed, 13 Mar 2024 10:45:33 +0100 Subject: [PATCH] Move two-layer shallow water equations from Trixi.jl (#28) * add two-layer swe * add elixirs * add tests * add unit test for input arguments of the 2LSWE --- ...lixir_shallowwater_twolayer_convergence.jl | 61 ++ .../elixir_shallowwater_twolayer_dam_break.jl | 95 +++ ...xir_shallowwater_twolayer_well_balanced.jl | 87 ++ ...lixir_shallowwater_twolayer_convergence.jl | 61 ++ ...xir_shallowwater_twolayer_well_balanced.jl | 82 ++ ...lixir_shallowwater_twolayer_convergence.jl | 64 ++ .../elixir_shallowwater_twolayer_dam_break.jl | 148 ++++ ...xir_shallowwater_twolayer_well_balanced.jl | 82 ++ src/TrixiShallowWater.jl | 9 +- src/equations/equations.jl | 3 + src/equations/shallow_water_two_layer_1d.jl | 507 +++++++++++ src/equations/shallow_water_two_layer_2d.jl | 804 ++++++++++++++++++ src/equations/shallow_water_wet_dry_1d.jl | 8 +- src/equations/shallow_water_wet_dry_2d.jl | 12 +- test/test_structured_2d.jl | 240 +++--- test/test_tree_1d.jl | 775 +++++++++-------- test/test_tree_2d.jl | 718 +++++++++------- test/test_unit.jl | 37 + test/test_unstructured_2d.jl | 754 ++++++++-------- 19 files changed, 3382 insertions(+), 1165 deletions(-) create mode 100644 examples/tree_1d_dgsem/elixir_shallowwater_twolayer_convergence.jl create mode 100644 examples/tree_1d_dgsem/elixir_shallowwater_twolayer_dam_break.jl create mode 100644 examples/tree_1d_dgsem/elixir_shallowwater_twolayer_well_balanced.jl create mode 100644 examples/tree_2d_dgsem/elixir_shallowwater_twolayer_convergence.jl create mode 100644 examples/tree_2d_dgsem/elixir_shallowwater_twolayer_well_balanced.jl create mode 100644 examples/unstructured_2d_dgsem/elixir_shallowwater_twolayer_convergence.jl create mode 100644 examples/unstructured_2d_dgsem/elixir_shallowwater_twolayer_dam_break.jl create mode 100644 examples/unstructured_2d_dgsem/elixir_shallowwater_twolayer_well_balanced.jl create mode 100644 src/equations/shallow_water_two_layer_1d.jl create mode 100644 src/equations/shallow_water_two_layer_2d.jl diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_twolayer_convergence.jl b/examples/tree_1d_dgsem/elixir_shallowwater_twolayer_convergence.jl new file mode 100644 index 0000000..0a8e285 --- /dev/null +++ b/examples/tree_1d_dgsem/elixir_shallowwater_twolayer_convergence.jl @@ -0,0 +1,61 @@ + +using OrdinaryDiffEq +using Trixi +using TrixiShallowWater + +############################################################################### +# Semidiscretization of the two-layer shallow water equations + +equations = ShallowWaterTwoLayerEquations1D(gravity_constant = 10.0, rho_upper = 0.9, + rho_lower = 1.0) + +initial_condition = initial_condition_convergence_test + +############################################################################### +# Get the DG approximation space + +volume_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal) +solver = DGSEM(polydeg = 3, + surface_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal), + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +############################################################################### +# Get the TreeMesh and setup a periodic mesh + +coordinates_min = 0.0 +coordinates_max = sqrt(2.0) +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, + source_terms = source_terms_convergence_test) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 1.0) +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 = 500, + 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, + save_everystep = false, callback = callbacks); +summary_callback() # print the timer summary diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_twolayer_dam_break.jl b/examples/tree_1d_dgsem/elixir_shallowwater_twolayer_dam_break.jl new file mode 100644 index 0000000..0554ef6 --- /dev/null +++ b/examples/tree_1d_dgsem/elixir_shallowwater_twolayer_dam_break.jl @@ -0,0 +1,95 @@ + +using OrdinaryDiffEq +using Trixi +using TrixiShallowWater + +############################################################################### +# Semidiscretization of the two-layer shallow water equations for a dam break +# test with a discontinuous bottom topography function to test entropy conservation + +equations = ShallowWaterTwoLayerEquations1D(gravity_constant = 9.81, H0 = 2.0, + rho_upper = 0.9, rho_lower = 1.0) + +# Initial condition of a dam break with a discontinuous water heights and bottom topography. +# Works as intended for TreeMesh1D with `initial_refinement_level=5`. If the mesh +# refinement level is changed the initial condition below may need changed as well to +# ensure that the discontinuities lie on an element interface. +function initial_condition_dam_break(x, t, equations::ShallowWaterTwoLayerEquations1D) + v1_upper = 0.0 + v1_lower = 0.0 + + # Set the discontinuity + if x[1] <= 10.0 + H_lower = 2.0 + H_upper = 4.0 + b = 0.0 + else + H_lower = 1.5 + H_upper = 3.0 + b = 0.5 + end + + return prim2cons(SVector(H_upper, v1_upper, H_lower, v1_lower, b), equations) +end + +initial_condition = initial_condition_dam_break + +############################################################################### +# Get the DG approximation space + +volume_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal) +solver = DGSEM(polydeg = 3, + surface_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal), + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +############################################################################### +# Get the TreeMesh and setup a non-periodic mesh + +coordinates_min = 0.0 +coordinates_max = 20.0 +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 5, + n_cells_max = 10000, + periodicity = false) + +boundary_condition = boundary_condition_slip_wall + +# create the semidiscretization object +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_condition) + +############################################################################### +# ODE solvers + +tspan = (0.0, 0.4) +ode = semidiscretize(semi, tspan) + +############################################################################### +# Callbacks + +summary_callback = SummaryCallback() + +analysis_interval = 500 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + save_analysis = false, + extra_analysis_integrals = (energy_total, + energy_kinetic, + energy_internal)) + +stepsize_callback = StepsizeCallback(cfl = 1.0) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 500, + 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, + save_everystep = false, callback = callbacks); +summary_callback() # print the timer summary diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_twolayer_well_balanced.jl b/examples/tree_1d_dgsem/elixir_shallowwater_twolayer_well_balanced.jl new file mode 100644 index 0000000..c206c63 --- /dev/null +++ b/examples/tree_1d_dgsem/elixir_shallowwater_twolayer_well_balanced.jl @@ -0,0 +1,87 @@ + +using OrdinaryDiffEq +using Trixi +using TrixiShallowWater + +############################################################################### +# Semidiscretization of the two-layer shallow water equations to test well-balancedness + +equations = ShallowWaterTwoLayerEquations1D(gravity_constant = 1.0, H0 = 0.6, + rho_upper = 0.9, rho_lower = 1.0) + +""" + initial_condition_fjordholm_well_balanced(x, t, equations::ShallowWaterTwoLayerEquations1D) + +Initial condition to test well balanced with a bottom topography from Fjordholm +- Ulrik Skre Fjordholm (2012) + Energy conservative and stable schemes for the two-layer shallow water equations. + [DOI: 10.1142/9789814417099_0039](https://doi.org/10.1142/9789814417099_0039) +""" +function initial_condition_fjordholm_well_balanced(x, t, + equations::ShallowWaterTwoLayerEquations1D) + inicenter = 0.5 + x_norm = x[1] - inicenter + r = abs(x_norm) + + H_lower = 0.5 + H_upper = 0.6 + v1_upper = 0.0 + v1_lower = 0.0 + b = r <= 0.1 ? 0.2 * (cos(10 * pi * (x[1] - 0.5)) + 1) : 0.0 + return prim2cons(SVector(H_upper, v1_upper, H_lower, v1_lower, b), equations) +end + +initial_condition = initial_condition_fjordholm_well_balanced + +############################################################################### +# Get the DG approximation space + +volume_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal) +solver = DGSEM(polydeg = 3, + surface_flux = (flux_es_ersing_etal, flux_nonconservative_ersing_etal), + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +############################################################################### +# Get the TreeMesh and setup a periodic mesh + +coordinates_min = 0.0 +coordinates_max = 1.0 +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 4, + 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, 10.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 1000 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + save_analysis = false, + extra_analysis_integrals = (lake_at_rest_error,)) + +stepsize_callback = StepsizeCallback(cfl = 1.0) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 1000, + save_initial_solution = true, + save_final_solution = true) + +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/examples/tree_2d_dgsem/elixir_shallowwater_twolayer_convergence.jl b/examples/tree_2d_dgsem/elixir_shallowwater_twolayer_convergence.jl new file mode 100644 index 0000000..cc63a61 --- /dev/null +++ b/examples/tree_2d_dgsem/elixir_shallowwater_twolayer_convergence.jl @@ -0,0 +1,61 @@ + +using OrdinaryDiffEq +using Trixi +using TrixiShallowWater + +############################################################################### +# Semidiscretization of the two-layer shallow water equations + +equations = ShallowWaterTwoLayerEquations2D(gravity_constant = 10.0, rho_upper = 0.9, + rho_lower = 1.0) + +initial_condition = initial_condition_convergence_test + +############################################################################### +# Get the DG approximation space + +volume_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal) +solver = DGSEM(polydeg = 3, + surface_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal), + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +############################################################################### +# Get the TreeMesh and setup a periodic mesh + +coordinates_min = (0.0, 0.0) +coordinates_max = (sqrt(2.0), sqrt(2.0)) +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 3, + n_cells_max = 20_000, + periodicity = true) + +# Create the semi discretization object +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + source_terms = source_terms_convergence_test) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 1.0) +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 = 500, + 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, + save_everystep = false, callback = callbacks); +summary_callback() # print the timer summary diff --git a/examples/tree_2d_dgsem/elixir_shallowwater_twolayer_well_balanced.jl b/examples/tree_2d_dgsem/elixir_shallowwater_twolayer_well_balanced.jl new file mode 100644 index 0000000..21bbfcf --- /dev/null +++ b/examples/tree_2d_dgsem/elixir_shallowwater_twolayer_well_balanced.jl @@ -0,0 +1,82 @@ + +using OrdinaryDiffEq +using Trixi +using TrixiShallowWater + +############################################################################### +# Semidiscretization of the two-layer shallow water equations with a bottom topography function +# to test well-balancedness + +equations = ShallowWaterTwoLayerEquations2D(gravity_constant = 9.81, H0 = 0.6, + rho_upper = 0.9, rho_lower = 1.0) + +# An initial condition with constant total water height, zero velocities and a bottom topography to +# test well-balancedness +function initial_condition_well_balanced(x, t, equations::ShallowWaterTwoLayerEquations2D) + H_lower = 0.5 + H_upper = 0.6 + v1_upper = 0.0 + v2_upper = 0.0 + v1_lower = 0.0 + v2_lower = 0.0 + b = (((x[1] - 0.5)^2 + (x[2] - 0.5)^2) < 0.04 ? + 0.2 * (cos(4 * pi * sqrt((x[1] - 0.5)^2 + (x[2] + + -0.5)^2)) + 1) : 0.0) + + return prim2cons(SVector(H_upper, v1_upper, v2_upper, H_lower, v1_lower, v2_lower, b), + equations) +end + +initial_condition = initial_condition_well_balanced + +############################################################################### +# Get the DG approximation space + +volume_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal) +surface_flux = (flux_es_ersing_etal, flux_nonconservative_ersing_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.0, 0.0) +coordinates_max = (1.0, 1.0) +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 solver + +tspan = (0.0, 10.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 1000 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + extra_analysis_integrals = (lake_at_rest_error,)) + +stepsize_callback = StepsizeCallback(cfl = 1.0) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 1000, + save_initial_solution = true, + save_final_solution = true) + +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/examples/unstructured_2d_dgsem/elixir_shallowwater_twolayer_convergence.jl b/examples/unstructured_2d_dgsem/elixir_shallowwater_twolayer_convergence.jl new file mode 100644 index 0000000..3783710 --- /dev/null +++ b/examples/unstructured_2d_dgsem/elixir_shallowwater_twolayer_convergence.jl @@ -0,0 +1,64 @@ + +using OrdinaryDiffEq +using Trixi +using TrixiShallowWater + +############################################################################### +# Semidiscretization of the two-layer shallow water equations with a periodic +# bottom topography function (set in the initial conditions) + +equations = ShallowWaterTwoLayerEquations2D(gravity_constant = 10.0, rho_upper = 0.9, + rho_lower = 1.0) + +initial_condition = initial_condition_convergence_test + +############################################################################### +# Get the DG approximation space + +volume_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal) +surface_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal) +solver = DGSEM(polydeg = 6, surface_flux = surface_flux, + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +############################################################################### +# This setup is for the curved, split form convergence test on a periodic domain + +# Get the unstructured quad mesh from a file (downloads the file if not available locally) +mesh_file = Trixi.download("https://gist.githubusercontent.com/andrewwinters5000/8f8cd23df27fcd494553f2a89f3c1ba4/raw/85e3c8d976bbe57ca3d559d653087b0889535295/mesh_alfven_wave_with_twist_and_flip.mesh", + joinpath(@__DIR__, "mesh_alfven_wave_with_twist_and_flip.mesh")) + +mesh = UnstructuredMesh2D(mesh_file, periodicity = true) + +# Create the semidiscretization object +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + source_terms = source_terms_convergence_test) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 1.0) +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 = 500, + 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/examples/unstructured_2d_dgsem/elixir_shallowwater_twolayer_dam_break.jl b/examples/unstructured_2d_dgsem/elixir_shallowwater_twolayer_dam_break.jl new file mode 100644 index 0000000..f01e839 --- /dev/null +++ b/examples/unstructured_2d_dgsem/elixir_shallowwater_twolayer_dam_break.jl @@ -0,0 +1,148 @@ + +using OrdinaryDiffEq +using Trixi +using TrixiShallowWater + +############################################################################### +# Semidiscretization of the two-layer shallow water equations for a dam break test with a +# discontinuous bottom topography function to test energy conservation + +equations = ShallowWaterTwoLayerEquations2D(gravity_constant = 1.0, rho_upper = 0.9, + rho_lower = 1.0) + +# This test case uses a special work around to setup a truly discontinuous bottom topography +# function and initial condition for this academic testcase of entropy conservation. First, a +# dummy initial_condition_dam_break is introduced to create the semidiscretization. Then the initial +# condition is reset with the true discontinuous values from initial_condition_discontinuous_dam_break. + +function initial_condition_dam_break(x, t, equations::ShallowWaterTwoLayerEquations2D) + if x[1] < sqrt(2) / 2 + H_upper = 1.0 + H_lower = 0.6 + b = 0.1 + else + H_upper = 0.9 + H_lower = 0.5 + b = 0.0 + end + + v1_upper = 0.0 + v2_upper = 0.0 + v1_lower = 0.0 + v2_lower = 0.0 + return prim2cons(SVector(H_upper, v1_upper, v2_upper, H_lower, v1_lower, v2_lower, b), + equations) +end + +initial_condition = initial_condition_dam_break + +boundary_condition_constant = BoundaryConditionDirichlet(initial_condition_dam_break) + +############################################################################### +# Get the DG approximation space + +volume_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal) +surface_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal) +solver = DGSEM(polydeg = 6, surface_flux = surface_flux, + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +############################################################################### +# Get the unstructured quad mesh from a file (downloads the file if not available locally) +mesh_file = Trixi.download("https://gist.githubusercontent.com/andrewwinters5000/8f8cd23df27fcd494553f2a89f3c1ba4/raw/85e3c8d976bbe57ca3d559d653087b0889535295/mesh_alfven_wave_with_twist_and_flip.mesh", + joinpath(@__DIR__, "mesh_alfven_wave_with_twist_and_flip.mesh")) + +mesh = UnstructuredMesh2D(mesh_file, periodicity = false) + +# Boundary conditions +boundary_condition = Dict(:Top => boundary_condition_slip_wall, + :Left => boundary_condition_slip_wall, + :Right => boundary_condition_slip_wall, + :Bottom => boundary_condition_slip_wall) + +# Create the semi discretization object +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, + solver, boundary_conditions = boundary_condition) + +############################################################################### +# ODE solver + +tspan = (0.0, 0.5) +ode = semidiscretize(semi, tspan) + +############################################################################### +# Workaround to set a discontinuous bottom topography and initial condition for debugging and testing. + +# alternative version of the initial conditinon used to setup a truly discontinuous +# test case and initial condition. +# In contrast to the usual signature of initial conditions, this one get passed the +# `element_id` explicitly. In particular, this initial conditions works as intended +# only for the specific mesh loaded above! + +function initial_condition_discontinuous_dam_break(x, t, element_id, + equations::ShallowWaterTwoLayerEquations2D) + # Constant values + v1_upper = 0.0 + v2_upper = 0.0 + v1_lower = 0.0 + v2_lower = 0.0 + + # Left side of discontinuity + IDs = [1, 2, 5, 6, 9, 10, 13, 14] + if element_id in IDs + H_upper = 1.0 + H_lower = 0.6 + b = 0.0 + # Right side of discontinuity + else + H_upper = 0.9 + H_lower = 0.5 + b = 0.1 + end + + return prim2cons(SVector(H_upper, v1_upper, v2_upper, H_lower, v1_lower, v2_lower, b), + equations) +end + +# point to the data we want to augment +u = Trixi.wrap_array(ode.u0, semi) +# reset the initial condition +for element in eachelement(semi.solver, semi.cache) + for j in eachnode(semi.solver), i in eachnode(semi.solver) + x_node = Trixi.get_node_coords(semi.cache.elements.node_coordinates, equations, + semi.solver, i, j, element) + u_node = initial_condition_discontinuous_dam_break(x_node, first(tspan), element, + equations) + Trixi.set_node_vars!(u, u_node, equations, semi.solver, i, j, element) + end +end + +############################################################################### +# Callbacks + +summary_callback = SummaryCallback() + +analysis_interval = 500 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + save_analysis = false, + extra_analysis_integrals = (energy_total, + energy_kinetic, + energy_internal)) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 500, + 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/examples/unstructured_2d_dgsem/elixir_shallowwater_twolayer_well_balanced.jl b/examples/unstructured_2d_dgsem/elixir_shallowwater_twolayer_well_balanced.jl new file mode 100644 index 0000000..fe5b168 --- /dev/null +++ b/examples/unstructured_2d_dgsem/elixir_shallowwater_twolayer_well_balanced.jl @@ -0,0 +1,82 @@ + +using OrdinaryDiffEq +using Trixi +using TrixiShallowWater + +############################################################################### +# Semidiscretization of the two-layer shallow water equations with a discontinuous bottom +# topography to test well-balancedness + +equations = ShallowWaterTwoLayerEquations2D(gravity_constant = 1.0, H0 = 0.6, + rho_upper = 0.9, rho_lower = 1.0) + +# An initial condition with constant total water height, zero velocities and a bottom topography to +# test well-balancedness +function initial_condition_well_balanced(x, t, equations::ShallowWaterTwoLayerEquations2D) + H_lower = 0.5 + H_upper = 0.6 + v1_upper = 0.0 + v2_upper = 0.0 + v1_lower = 0.0 + v2_lower = 0.0 + + # Bottom Topography + b = (((x[1] - 0.5)^2 + (x[2] - 0.5)^2) < 0.04 ? + 0.2 * (cos(4 * pi * sqrt((x[1] - 0.5)^2 + (x[2] + + -0.5)^2)) + 1) : 0.0) + return prim2cons(SVector(H_upper, v1_upper, v2_upper, H_lower, v1_lower, v2_lower, b), + equations) +end + +initial_condition = initial_condition_well_balanced + +############################################################################### +# Get the DG approximation space + +volume_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal) +surface_flux = (flux_es_ersing_etal, flux_nonconservative_ersing_etal) +solver = DGSEM(polydeg = 6, surface_flux = surface_flux, + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +############################################################################### +# This setup is for the curved, split form well-balancedness testing + +# Get the unstructured quad mesh from a file (downloads the file if not available locally) +mesh_file = Trixi.download("https://gist.githubusercontent.com/andrewwinters5000/8f8cd23df27fcd494553f2a89f3c1ba4/raw/85e3c8d976bbe57ca3d559d653087b0889535295/mesh_alfven_wave_with_twist_and_flip.mesh", + joinpath(@__DIR__, "mesh_alfven_wave_with_twist_and_flip.mesh")) + +mesh = UnstructuredMesh2D(mesh_file, periodicity = true) + +# Create the semi discretization object +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +############################################################################### +# ODE solver + +tspan = (0.0, 10.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 1000 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + extra_analysis_integrals = (lake_at_rest_error,)) + +stepsize_callback = StepsizeCallback(cfl = 1.0) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 1000, + save_initial_solution = true, + save_final_solution = true) + +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/TrixiShallowWater.jl b/src/TrixiShallowWater.jl index c065f79..3932671 100644 --- a/src/TrixiShallowWater.jl +++ b/src/TrixiShallowWater.jl @@ -6,9 +6,9 @@ module TrixiShallowWater # https://github.com/trixi-framework/TrixiShallowWater.jl/pull/10#discussion_r1433720559 using Trixi # Import additional symbols that are not exported by Trixi.jl -using Trixi: get_node_vars, set_node_vars! +using Trixi: get_node_vars, set_node_vars!, waterheight using MuladdMacro: @muladd -using StaticArrays: SVector +using StaticArrays: SVector, @SMatrix using Static: True, False using LinearAlgebra: norm @@ -18,12 +18,15 @@ include("callbacks_stage/callbacks_stage.jl") include("solvers/indicators.jl") # Export types/functions that define the public API of TrixiShallowWater.jl -export ShallowWaterEquationsWetDry1D, ShallowWaterEquationsWetDry2D +export ShallowWaterEquationsWetDry1D, ShallowWaterEquationsWetDry2D, + ShallowWaterTwoLayerEquations1D, ShallowWaterTwoLayerEquations2D export hydrostatic_reconstruction_chen_noelle, flux_nonconservative_chen_noelle, min_max_speed_chen_noelle, flux_hll_chen_noelle +export flux_es_ersing_etal + export PositivityPreservingLimiterShallowWater export IndicatorHennemannGassnerShallowWater diff --git a/src/equations/equations.jl b/src/equations/equations.jl index 9fb8148..076a994 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -10,4 +10,7 @@ include("shallow_water_wet_dry_1d.jl") include("shallow_water_wet_dry_2d.jl") + +include("shallow_water_two_layer_1d.jl") +include("shallow_water_two_layer_2d.jl") end # @muladd diff --git a/src/equations/shallow_water_two_layer_1d.jl b/src/equations/shallow_water_two_layer_1d.jl new file mode 100644 index 0000000..cd114b4 --- /dev/null +++ b/src/equations/shallow_water_two_layer_1d.jl @@ -0,0 +1,507 @@ +# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). +# Since these FMAs can increase the performance of many numerical algorithms, +# we need to opt-in explicitly. +# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. +@muladd begin +#! format: noindent + +@doc raw""" + ShallowWaterTwoLayerEquations1D(gravity, H0, rho_upper, rho_lower) + +Two-Layer Shallow Water equations (2LSWE) in one space dimension. The equations are given by +```math +\begin{alignat*}{4} +&\frac{\partial}{\partial t}h_{upper} +&&+ \frac{\partial}{\partial x}\left(h_{upper} v_{1,upper}\right) +&&= 0 \\ +&\frac{\partial}{\partial t}\left(h_{upper}v_{1,upper}\right) +&&+ \frac{\partial}{\partial x}\left(h_{upper}v_{1,upper}^2 + \dfrac{gh_{upper}^2}{2}\right) +&&= -gh_{upper}\frac{\partial}{\partial x}\left(b+h_{lower}\right)\\ +&\frac{\partial}{\partial t}h_{lower} +&&+ \frac{\partial}{\partial x}\left(h_{lower}v_{1,lower}\right) +&&= 0 \\ +&\frac{\partial}{\partial t}\left(h_{lower}v_{1,lower}\right) +&&+ \frac{\partial}{\partial x}\left(h_{lower}v_{1,lower}^2 + \dfrac{gh_{lower}^2}{2}\right) +&&= -gh_{lower}\frac{\partial}{\partial x}\left(b+\dfrac{\rho_{upper}}{\rho_{lower}}h_{upper}\right). +\end{alignat*} +``` +The unknown quantities of the 2LSWE are the water heights of the {lower} layer ``h_{lower}`` and the +{upper} layer ``h_{upper}`` with respective velocities ``v_{1,upper}`` and ``v_{1,lower}``. The gravitational constant is +denoted by `g`, the layer densitites by ``\rho_{upper}``and ``\rho_{lower}`` and the (possibly) variable +bottom topography function ``b(x)``. The conservative variable water height ``h_{lower}`` is measured +from the bottom topography ``b`` and ``h_{upper}`` relative to ``h_{lower}``, therefore one also defines the +total water heights as ``H_{upper} = h_{upper} + h_{upper} + b`` and ``H_{lower} = h_{lower} + b``. + +The densities must be chosen such that ``\rho_{upper} < \rho_{lower}``, to make sure that the heavier fluid +``\rho_{lower}`` is in the bottom layer and the lighter fluid ``\rho_{upper}`` in the {upper} layer. + +The additional quantity ``H_0`` is also available to store a reference value for the total water +height that is useful to set initial conditions or test the "lake-at-rest" well-balancedness. + +The bottom topography function ``b(x)`` is set inside the initial condition routine +for a particular problem setup. + +In addition to the unknowns, Trixi currently stores the bottom topography values at the +approximation points despite being fixed in time. This is done for convenience of computing the +bottom topography gradients on the fly during the approximation as well as computing auxiliary +quantities like the total water height ``H`` or the entropy variables. +This affects the implementation and use of these equations in various ways: +* The flux values corresponding to the bottom topography must be zero. +* The bottom topography values must be included when defining initial conditions, boundary + conditions or source terms. +* [`AnalysisCallback`](@ref) analyzes this variable. +* Trixi's visualization tools will visualize the bottom topography by default. + +A good introduction for the 2LSWE is available in Chapter 12 of the book: +- Benoit Cushman-Roisin (2011)\ + Introduction to geophyiscal fluid dynamics: physical and numerical aspects\ + \ + ISBN: 978-0-12-088759-0 +""" +struct ShallowWaterTwoLayerEquations1D{RealT <: Real} <: + Trixi.AbstractShallowWaterEquations{1, 5} + gravity::RealT # gravitational constant + H0::RealT # constant "lake-at-rest" total water height + rho_upper::RealT # lower layer density + rho_lower::RealT # upper layer density + r::RealT # ratio of rho_upper / rho_lower +end + +# Allow for flexibility to set the gravitational constant within an elixir depending on the +# application where `gravity_constant=1.0` or `gravity_constant=9.81` are common values. +# The reference total water height H0 defaults to 0.0 but is used for the "lake-at-rest" +# well-balancedness test cases. Densities must be specified such that rho_upper <= rho_lower. +function ShallowWaterTwoLayerEquations1D(; gravity_constant, + H0 = zero(gravity_constant), rho_upper, + rho_lower) + # Assign density ratio if rho_upper <= rho_lower + if rho_upper > rho_lower + error("Invalid input: Densities must be chosen such that rho_upper < rho_lower") + else + r = rho_upper / rho_lower + end + ShallowWaterTwoLayerEquations1D(gravity_constant, H0, rho_upper, rho_lower, r) +end + +Trixi.have_nonconservative_terms(::ShallowWaterTwoLayerEquations1D) = True() +function Trixi.varnames(::typeof(cons2cons), ::ShallowWaterTwoLayerEquations1D) + ("h_upper", "h_v_upper", + "h_lower", "h_v_lower", "b") +end +# Note, we use the total water height, H_lower = h_upper + h_lower + b, and first layer total height +# H_upper = h_upper + b as the first primitive variable for easier visualization and setting initial +# conditions +function Trixi.varnames(::typeof(cons2prim), ::ShallowWaterTwoLayerEquations1D) + ("H_upper", "v_upper", + "H_lower", "v_lower", "b") +end + +# Set initial conditions at physical location `x` for time `t` +""" + initial_condition_convergence_test(x, t, equations::ShallowWaterTwoLayerEquations1D) + +A smooth initial condition used for convergence tests in combination with +[`source_terms_convergence_test`](@ref) (and +[`BoundaryConditionDirichlet(initial_condition_convergence_test)`](@ref) in non-periodic domains). +""" +function Trixi.initial_condition_convergence_test(x, t, + equations::ShallowWaterTwoLayerEquations1D) + # some constants are chosen such that the function is periodic on the domain [0,sqrt(2)] + ω = 2.0 * pi * sqrt(2.0) + + H_lower = 2.0 + 0.1 * sin(ω * x[1] + t) + H_upper = 4.0 + 0.1 * cos(ω * x[1] + t) + v_lower = 1.0 + v_upper = 0.9 + b = 1.0 + 0.1 * cos(2.0 * ω * x[1]) + + return prim2cons(SVector(H_upper, v_upper, H_lower, v_lower, b), equations) +end + +""" + source_terms_convergence_test(u, x, t, equations::ShallowWaterTwoLayerEquations1D) + +Source terms used for convergence tests in combination with +[`initial_condition_convergence_test`](@ref) +(and [`BoundaryConditionDirichlet(initial_condition_convergence_test)`](@ref) +in non-periodic domains). +""" +@inline function Trixi.source_terms_convergence_test(u, x, t, + equations::ShallowWaterTwoLayerEquations1D) + # Same settings as in `initial_condition_convergence_test`. Some derivative simplify because + # this manufactured solution velocity is taken to be constant + ω = 2 * pi * sqrt(2.0) + + du1 = (-0.1 * cos(t + ω * x[1]) - 0.1 * sin(t + ω * x[1]) - + 0.09 * ω * cos(t + ω * x[1]) + + -0.09 * ω * sin(t + ω * x[1])) + du2 = (5.0 * (-0.1 * ω * cos(t + ω * x[1]) - 0.1 * ω * sin(t + ω * x[1])) * + (4.0 + 0.2 * cos(t + ω * x[1]) + + -0.2 * sin(t + ω * x[1])) + + 0.1 * ω * (20.0 + cos(t + ω * x[1]) - sin(t + ω * x[1])) * + cos(t + + ω * x[1]) - 0.09 * cos(t + ω * x[1]) - 0.09 * sin(t + ω * x[1]) - + 0.081 * ω * cos(t + ω * x[1]) + + -0.081 * ω * sin(t + ω * x[1])) + du3 = 0.1 * cos(t + ω * x[1]) + 0.1 * ω * cos(t + ω * x[1]) + + 0.2 * ω * sin(2.0 * ω * x[1]) + du4 = ((10.0 + sin(t + ω * x[1]) - cos(2ω * x[1])) * + (-0.09 * ω * cos(t + ω * x[1]) - 0.09 * ω * sin(t + + ω * x[1]) - + 0.2 * ω * sin(2 * ω * x[1])) + 0.1 * cos(t + ω * x[1]) + + 0.1 * ω * cos(t + ω * x[1]) + + 5.0 * (0.1 * ω * cos(t + ω * x[1]) + 0.2 * ω * sin(2.0 * ω * x[1])) * + (2.0 + 0.2 * sin(t + ω * x[1]) + + -0.2 * cos(2.0 * ω * x[1])) + 0.2 * ω * sin(2.0 * ω * x[1])) + + return SVector(du1, du2, du3, du4, zero(eltype(u))) +end + +""" + boundary_condition_slip_wall(u_inner, orientation_or_normal, x, t, surface_flux_function, + equations::ShallowWaterTwoLayerEquations1D) + +Create a boundary state by reflecting the normal velocity component and keep +the tangential velocity component unchanged. The boundary water height is taken from +the internal value. + +For details see Section 9.2.5 of the book: +- Eleuterio F. Toro (2001) + Shock-Capturing Methods for Free-Surface Shallow Flows + 1st edition + ISBN 0471987662 +""" +@inline function Trixi.boundary_condition_slip_wall(u_inner, orientation_or_normal, + direction, + x, t, surface_flux_function, + equations::ShallowWaterTwoLayerEquations1D) + # create the "external" boundary solution state + u_boundary = SVector(u_inner[1], + -u_inner[2], + u_inner[3], + -u_inner[4], + u_inner[5]) + + # calculate the boundary flux + if iseven(direction) # u_inner is "left" of boundary, u_boundary is "right" of boundary + f = surface_flux_function(u_inner, u_boundary, orientation_or_normal, equations) + else # u_boundary is "left" of boundary, u_inner is "right" of boundary + f = surface_flux_function(u_boundary, u_inner, orientation_or_normal, equations) + end + return f +end + +# Calculate 1D flux for a single point +# Note, the bottom topography has no flux +@inline function Trixi.flux(u, orientation::Integer, + equations::ShallowWaterTwoLayerEquations1D) + h_upper, h_v_upper, h_lower, h_v_lower, _ = u + + # Calculate velocities + v_upper, v_lower = velocity(u, equations) + # Calculate pressure + p_upper = 0.5 * equations.gravity * h_upper^2 + p_lower = 0.5 * equations.gravity * h_lower^2 + + f1 = h_v_upper + f2 = h_v_upper * v_upper + p_upper + f3 = h_v_lower + f4 = h_v_lower * v_lower + p_lower + + return SVector(f1, f2, f3, f4, zero(eltype(u))) +end + +""" + flux_nonconservative_ersing_etal(u_ll, u_rr, orientation::Integer, + equations::ShallowWaterTwoLayerEquations1D) + +Non-symmetric path-conservative two-point volume flux discretizing the nonconservative (source) term +that contains the gradient of the bottom topography [`ShallowWaterTwoLayerEquations1D`](@ref) and an +additional term that couples the momentum of both layers. + +This is a modified version of [`flux_nonconservative_wintermeyer_etal`](@ref) that gives entropy +conservation and well-balancedness in both the volume and surface when combined with +[`flux_wintermeyer_etal`](@ref). + +For further details see: +- Patrick Ersing, Andrew R. Winters (2023) + An entropy stable discontinuous Galerkin method for the two-layer shallow water equations on + curvilinear meshes + [DOI: 10.48550/arXiv.2306.12699](https://doi.org/10.48550/arXiv.2306.12699) +""" +@inline function Trixi.flux_nonconservative_ersing_etal(u_ll, u_rr, + orientation::Integer, + equations::ShallowWaterTwoLayerEquations1D) + # Pull the necessary left and right state information + h_upper_ll, h_lower_ll = waterheight(u_ll, equations) + h_upper_rr, h_lower_rr = waterheight(u_rr, equations) + b_rr = u_rr[5] + b_ll = u_ll[5] + + # Calculate jumps + h_upper_jump = (h_upper_rr - h_upper_ll) + h_lower_jump = (h_lower_rr - h_lower_ll) + b_jump = (b_rr - b_ll) + + z = zero(eltype(u_ll)) + + # Bottom gradient nonconservative term: (0, g*h_upper*(b+h_lower)_x, + # 0, g*h_lower*(b+r*h_upper)_x, 0) + f = SVector(z, + equations.gravity * h_upper_ll * (b_jump + h_lower_jump), + z, + equations.gravity * h_lower_ll * (b_jump + equations.r * h_upper_jump), + z) + return f +end + +""" + flux_wintermeyer_etal(u_ll, u_rr, orientation, + equations::ShallowWaterTwoLayerEquations1D) + +Total energy conservative (mathematical entropy for two-layer shallow water equations) split form. +When the bottom topography is nonzero this scheme will be well-balanced when used with the +nonconservative [`flux_nonconservative_ersing_etal`](@ref). To obtain the flux for the +two-layer shallow water equations the flux that is described in the paper for the normal shallow +water equations is used within each layer. + +Further details are available in Theorem 1 of the paper: +- Niklas Wintermeyer, Andrew R. Winters, Gregor J. Gassner and David A. Kopriva (2017) + An entropy stable nodal discontinuous Galerkin method for the two dimensional + shallow water equations on unstructured curvilinear meshes with discontinuous bathymetry + [DOI: 10.1016/j.jcp.2017.03.036](https://doi.org/10.1016/j.jcp.2017.03.036) +""" +@inline function Trixi.flux_wintermeyer_etal(u_ll, u_rr, + orientation::Integer, + equations::ShallowWaterTwoLayerEquations1D) + # Unpack left and right state + h_upper_ll, h_v_upper_ll, h_lower_ll, h_v_lower_ll, _ = u_ll + h_upper_rr, h_v_upper_rr, h_lower_rr, h_v_lower_rr, _ = u_rr + + # Get the velocities on either side + v_upper_ll, v_lower_ll = velocity(u_ll, equations) + v_upper_rr, v_lower_rr = velocity(u_rr, equations) + + # Average each factor of products in flux + v_upper_avg = 0.5 * (v_upper_ll + v_upper_rr) + v_lower_avg = 0.5 * (v_lower_ll + v_lower_rr) + p_upper_avg = 0.5 * equations.gravity * h_upper_ll * h_upper_rr + p_lower_avg = 0.5 * equations.gravity * h_lower_ll * h_lower_rr + + # Calculate fluxes + f1 = 0.5 * (h_v_upper_ll + h_v_upper_rr) + f2 = f1 * v_upper_avg + p_upper_avg + f3 = 0.5 * (h_v_lower_ll + h_v_lower_rr) + f4 = f3 * v_lower_avg + p_lower_avg + + return SVector(f1, f2, f3, f4, zero(eltype(u_ll))) +end + +""" + flux_es_ersing_etal(u_ll, u_rr, orientation_or_normal_direction, + equations::ShallowWaterTwoLayerEquations1D) +Entropy stable surface flux for the two-layer shallow water equations. Uses the entropy conservative +[`flux_wintermeyer_etal`](@ref) and adds a Lax-Friedrichs type dissipation dependent on the jump of +entropy variables. + +For further details see: +- Patrick Ersing, Andrew R. Winters (2023) + An entropy stable discontinuous Galerkin method for the two-layer shallow water equations on + curvilinear meshes + [DOI: 10.48550/arXiv.2306.12699](https://doi.org/10.48550/arXiv.2306.12699) +""" +@inline function flux_es_ersing_etal(u_ll, u_rr, + orientation::Integer, + equations::ShallowWaterTwoLayerEquations1D) + # Compute entropy conservative flux but without the bottom topography + f_ec = flux_wintermeyer_etal(u_ll, u_rr, + orientation, + equations) + + # Get maximum signal velocity + λ = max_abs_speed_naive(u_ll, u_rr, orientation, equations) + # Get entropy variables but without the bottom topography + q_rr = cons2entropy(u_rr, equations) + q_ll = cons2entropy(u_ll, equations) + + # Average values from left and right + u_avg = (u_ll + u_rr) / 2 + + # Introduce variables for better readability + rho_upper = equations.rho_upper + rho_lower = equations.rho_lower + g = equations.gravity + drho = rho_upper - rho_lower + + # Compute entropy Jacobian coefficients + h11 = -rho_lower / (g * rho_upper * drho) + h12 = -rho_lower * u_avg[2] / (g * rho_upper * u_avg[1] * drho) + h13 = 1.0 / (g * drho) + h14 = u_avg[4] / (g * u_avg[3] * drho) + h21 = -rho_lower * u_avg[2] / (g * rho_upper * u_avg[1] * drho) + h22 = ((g * rho_upper * u_avg[1]^3 - g * rho_lower * u_avg[1]^3 + + -rho_lower * u_avg[2]^2) / (g * rho_upper * u_avg[1]^2 * drho)) + h23 = u_avg[2] / (g * u_avg[1] * drho) + h24 = u_avg[2] * u_avg[4] / (g * u_avg[1] * u_avg[3] * drho) + h31 = 1.0 / (g * drho) + h32 = u_avg[2] / (g * u_avg[1] * drho) + h33 = -1.0 / (g * drho) + h34 = -u_avg[4] / (g * u_avg[3] * drho) + h41 = u_avg[4] / (g * u_avg[3] * drho) + h42 = u_avg[2] * u_avg[4] / (g * u_avg[1] * u_avg[3] * drho) + h43 = -u_avg[4] / (g * u_avg[3] * drho) + h44 = ((g * rho_upper * u_avg[3]^3 - g * rho_lower * u_avg[3]^3 + + -rho_lower * u_avg[4]^2) / (g * rho_lower * u_avg[3]^2 * drho)) + + # Entropy Jacobian matrix + H = @SMatrix [[h11;; h12;; h13;; h14;; 0]; + [h21;; h22;; h23;; h24;; 0]; + [h31;; h32;; h33;; h34;; 0]; + [h41;; h42;; h43;; h44;; 0]; + [0;; 0;; 0;; 0;; 0]] + + # Add dissipation to entropy conservative flux to obtain entropy stable flux + f_es = f_ec - 0.5 * λ * H * (q_rr - q_ll) + + return SVector(f_es[1], f_es[2], f_es[3], f_es[4], zero(eltype(u_ll))) +end + +# Calculate approximation for maximum wave speed for local Lax-Friedrichs-type dissipation as the +# maximum velocity magnitude plus the maximum speed of sound. This function uses approximate +# eigenvalues using the speed of the barotropic mode as there is no simple way to calculate them +# analytically. +# +# A good overview of the derivation is given in: +# - Jonas Nycander, Andrew McC. Hogg, Leela M. Frankcombe (2008) +# Open boundary conditions for nonlinear channel Flows +# [DOI: 10.1016/j.ocemod.2008.06.003](https://doi.org/10.1016/j.ocemod.2008.06.003) +@inline function Trixi.max_abs_speed_naive(u_ll, u_rr, + orientation::Integer, + equations::ShallowWaterTwoLayerEquations1D) + # Unpack left and right state + h_upper_ll, h_v_upper_ll, h_lower_ll, h_v_lower_ll, _ = u_ll + h_upper_rr, h_v_upper_rr, h_lower_rr, h_v_lower_rr, _ = u_rr + + # Get the averaged velocity + v_m_ll = (h_v_upper_ll + h_v_lower_ll) / (h_upper_ll + h_lower_ll) + v_m_rr = (h_v_upper_rr + h_v_lower_rr) / (h_upper_rr + h_lower_rr) + + # Calculate the wave celerity on the left and right + h_upper_ll, h_lower_ll = waterheight(u_ll, equations) + h_upper_rr, h_lower_rr = waterheight(u_rr, equations) + c_ll = sqrt(equations.gravity * (h_upper_ll + h_lower_ll)) + c_rr = sqrt(equations.gravity * (h_upper_rr + h_lower_rr)) + + return (max(abs(v_m_ll) + c_ll, abs(v_m_rr) + c_rr)) +end + +# Specialized `DissipationLocalLaxFriedrichs` to avoid spurious dissipation in the bottom +# topography +@inline function (dissipation::DissipationLocalLaxFriedrichs)(u_ll, u_rr, + orientation_or_normal_direction, + equations::ShallowWaterTwoLayerEquations1D) + λ = dissipation.max_abs_speed(u_ll, u_rr, orientation_or_normal_direction, + equations) + diss = -0.5 * λ * (u_rr - u_ll) + return SVector(diss[1], diss[2], diss[3], diss[4], zero(eltype(u_ll))) +end + +# Absolute speed of the barotropic mode +@inline function Trixi.max_abs_speeds(u, equations::ShallowWaterTwoLayerEquations1D) + h_upper, h_v_upper, h_lower, h_v_lower, _ = u + + # Calculate averaged velocity of both layers + v_m = (h_v_upper + h_v_lower) / (h_upper + h_lower) + c = sqrt(equations.gravity * (h_upper + h_lower)) + + return (abs(v_m) + c) +end + +# Helper function to extract the velocity vector from the conservative variables +@inline function Trixi.velocity(u, equations::ShallowWaterTwoLayerEquations1D) + h_upper, h_v_upper, h_lower, h_v_lower, _ = u + + v_upper = h_v_upper / h_upper + v_lower = h_v_lower / h_lower + return SVector(v_upper, v_lower) +end + +# Convert conservative variables to primitive +@inline function Trixi.cons2prim(u, equations::ShallowWaterTwoLayerEquations1D) + h_upper, _, h_lower, _, b = u + + H_lower = h_lower + b + H_upper = h_lower + h_upper + b + v_upper, v_lower = velocity(u, equations) + return SVector(H_upper, v_upper, H_lower, v_lower, b) +end + +# Convert conservative variables to entropy variables +# Note, only the first four are the entropy variables, the fifth entry still just carries the +# bottom topography values for convenience +@inline function Trixi.cons2entropy(u, equations::ShallowWaterTwoLayerEquations1D) + h_upper, _, h_lower, _, b = u + v_upper, v_lower = velocity(u, equations) + + w1 = (equations.rho_upper * + (equations.gravity * (h_upper + h_lower + b) - 0.5 * v_upper^2)) + w2 = equations.rho_upper * v_upper + w3 = (equations.rho_lower * + (equations.gravity * (equations.r * h_upper + h_lower + b) - 0.5 * v_lower^2)) + w4 = equations.rho_lower * v_lower + return SVector(w1, w2, w3, w4, b) +end + +# Convert primitive to conservative variables +@inline function Trixi.prim2cons(prim, equations::ShallowWaterTwoLayerEquations1D) + H_upper, v_upper, H_lower, v_lower, b = prim + + h_lower = H_lower - b + h_upper = H_upper - h_lower - b + h_v_upper = h_upper * v_upper + h_v_lower = h_lower * v_lower + return SVector(h_upper, h_v_upper, h_lower, h_v_lower, b) +end + +@inline function Trixi.waterheight(u, equations::ShallowWaterTwoLayerEquations1D) + return SVector(u[1], u[3]) +end + +# Entropy function for the shallow water equations is the total energy +@inline function Trixi.entropy(cons, equations::ShallowWaterTwoLayerEquations1D) + energy_total(cons, equations) +end + +# Calculate total energy for a conservative state `cons` +@inline function Trixi.energy_total(cons, equations::ShallowWaterTwoLayerEquations1D) + h_upper, h_v_upper, h_lower, h_v_lower, b = cons + # Set new variables for better readability + g = equations.gravity + rho_upper = equations.rho_upper + rho_lower = equations.rho_lower + + e = (0.5 * rho_upper * (h_v_upper^2 / h_upper + g * h_upper^2) + + 0.5 * rho_lower * (h_v_lower^2 / h_lower + g * h_lower^2) + + g * rho_lower * h_lower * b + g * rho_upper * h_upper * (h_lower + b)) + return e +end + +# Calculate kinetic energy for a conservative state `cons` +@inline function Trixi.energy_kinetic(u, equations::ShallowWaterTwoLayerEquations1D) + h_upper, h_v_upper, h_lower, h_v_lower, _ = u + return (0.5 * equations.rho_upper * h_v_upper^2 / h_upper + + 0.5 * equations.rho_lower * h_v_lower^2 / h_lower) +end + +# Calculate potential energy for a conservative state `cons` +@inline function Trixi.energy_internal(cons, equations::ShallowWaterTwoLayerEquations1D) + return energy_total(cons, equations) - energy_kinetic(cons, equations) +end + +# Calculate the error for the "lake-at-rest" test case where H = h_upper+h_lower+b should +# be a constant value over time +@inline function Trixi.lake_at_rest_error(u, equations::ShallowWaterTwoLayerEquations1D) + h_upper, _, h_lower, _, b = u + return abs(equations.H0 - (h_upper + h_lower + b)) +end +end # @muladd diff --git a/src/equations/shallow_water_two_layer_2d.jl b/src/equations/shallow_water_two_layer_2d.jl new file mode 100644 index 0000000..ac0df7d --- /dev/null +++ b/src/equations/shallow_water_two_layer_2d.jl @@ -0,0 +1,804 @@ +# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). +# Since these FMAs can increase the performance of many numerical algorithms, +# we need to opt-in explicitly. +# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. +@muladd begin +#! format: noindent + +@doc raw""" + ShallowWaterTwoLayerEquations2D(gravity, H0, rho_upper, rho_lower) + +Two-Layer Shallow water equations (2LSWE) in two space dimension. The equations are given by +```math +\begin{alignat*}{8} +&\frac{\partial}{\partial t}h_{upper} +&&+ \frac{\partial}{\partial x}\left(h_{upper} v_{1,upper}\right) +&&+ \frac{\partial}{\partial y}\left(h_{upper} v_{2,upper}\right) \quad +&&= \quad 0 \\ +&\frac{\partial}{\partial t}\left(h_{upper} v_{1,upper}\right) +&&+ \frac{\partial}{\partial x}\left(h_{upper} v_{1,upper}^2 + \frac{gh_{upper}^2}{2}\right) +&&+ \frac{\partial}{\partial y}\left(h_{upper} v_{1,upper} v_{2,upper}\right) \quad +&&= -gh_{upper}\frac{\partial}{\partial x}\left(b+h_{lower}\right) \\ +&\frac{\partial}{\partial t}\left(h_{upper} v_{2,upper}\right) +&&+ \frac{\partial}{\partial x}\left(h_{upper} v_{1,upper} v_{2,upper}\right) +&&+ \frac{\partial}{\partial y}\left(h_{upper} v_{2,upper}^2 + \frac{gh_{upper}^2}{2}\right) +&&= -gh_{upper}\frac{\partial}{\partial y}\left(b+h_{lower}\right)\\ +&\frac{\partial}{\partial t}h_{lower} +&&+ \frac{\partial}{\partial x}\left(h_{lower} v_{1,lower}\right) +&&+ \frac{\partial}{\partial y}\left(h_{lower} v_{2,lower}\right) +&&= \quad 0 \\ +&\frac{\partial}{\partial t}\left(h_{lower} v_{1,lower}\right) +&&+ \frac{\partial}{\partial x}\left(h_{lower} v_{1,lower}^2 + \frac{gh_{lower}^2}{2}\right) +&&+ \frac{\partial}{\partial y}\left(h_{lower} v_{1,lower} v_{2,lower}\right) +&&= -gh_{lower}\frac{\partial}{\partial x}\left(b+\frac{\rho_{upper}}{\rho_{lower}} h_{upper}\right)\\ +&\frac{\partial}{\partial t}\left(h_{lower} v_{2,lower}\right) +&&+ \frac{\partial}{\partial x}\left(h_{lower} v_{1,lower} v_{2,lower}\right) +&&+ \frac{\partial}{\partial y}\left(h_{lower} v_{2,lower}^2 + \frac{gh_{lower}^2}{2}\right) +&&= -gh_{lower}\frac{\partial}{\partial y}\left(b+\frac{\rho_{upper}}{\rho_{lower}} h_{upper}\right) +\end{alignat*} +``` +The unknown quantities of the 2LSWE are the water heights of the lower layer ``h_{lower}`` and the +upper +layer ``h_{upper}`` and the respective velocities in x-direction ``v_{1,lower}`` and ``v_{1,upper}`` and in y-direction +``v_{2,lower}`` and ``v_{2,upper}``. The gravitational constant is denoted by `g`, the layer densitites by +``\rho_{upper}``and ``\rho_{lower}`` and the (possibly) variable bottom topography function by ``b(x)``. +Conservative variable water height ``h_{lower}`` is measured from the bottom topography ``b`` and ``h_{upper}`` +relative to ``h_{lower}``, therefore one also defines the total water heights as ``H_{lower} = h_{lower} + b`` and +``H_{upper} = h_{upper} + h_{lower} + b``. + +The densities must be chosen such that ``\rho_{upper} < \rho_{lower}``, to make sure that the heavier fluid +``\rho_{lower}`` is in the bottom layer and the lighter fluid ``\rho_{upper}`` in the upper layer. + +The additional quantity ``H_0`` is also available to store a reference value for the total water +height that is useful to set initial conditions or test the "lake-at-rest" well-balancedness. + +The bottom topography function ``b(x)`` is set inside the initial condition routine +for a particular problem setup. + +In addition to the unknowns, Trixi currently stores the bottom topography values at the +approximation points despite being fixed in time. This is done for convenience of computing the +bottom topography gradients on the fly during the approximation as well as computing auxiliary +quantities like the total water height ``H`` or the entropy variables. +This affects the implementation and use of these equations in various ways: +* The flux values corresponding to the bottom topography must be zero. +* The bottom topography values must be included when defining initial conditions, boundary + conditions or source terms. +* [`AnalysisCallback`](@ref) analyzes this variable. +* Trixi's visualization tools will visualize the bottom topography by default. + +A good introduction for the 2LSWE is available in Chapter 12 of the book: + - Benoit Cushman-Roisin (2011)\ + Introduction to geophyiscal fluid dynamics: physical and numerical aspects\ + \ + ISBN: 978-0-12-088759-0 +""" +struct ShallowWaterTwoLayerEquations2D{RealT <: Real} <: + Trixi.AbstractShallowWaterEquations{2, 7} + gravity::RealT # gravitational constant + H0::RealT # constant "lake-at-rest" total water height + rho_upper::RealT # lower layer density + rho_lower::RealT # upper layer density + r::RealT # ratio of rho_upper / rho_lower +end + +# Allow for flexibility to set the gravitational constant within an elixir depending on the +# application where `gravity_constant=1.0` or `gravity_constant=9.81` are common values. +# The reference total water height H0 defaults to 0.0 but is used for the "lake-at-rest" +# well-balancedness test cases. Densities must be specified such that rho_upper < rho_lower. +function ShallowWaterTwoLayerEquations2D(; gravity_constant, + H0 = zero(gravity_constant), rho_upper, + rho_lower) + # Assign density ratio if rho_upper <= rho_lower + if rho_upper > rho_lower + error("Invalid input: Densities must be chosen such that rho_upper < rho_lower") + else + r = rho_upper / rho_lower + end + ShallowWaterTwoLayerEquations2D(gravity_constant, H0, rho_upper, rho_lower, r) +end + +Trixi.have_nonconservative_terms(::ShallowWaterTwoLayerEquations2D) = True() +function Trixi.varnames(::typeof(cons2cons), ::ShallowWaterTwoLayerEquations2D) + ("h_upper", "h_v1_upper", "h_v2_upper", "h_lower", "h_v1_lower", "h_v2_lower", "b") +end +# Note, we use the total water height, H_upper = h_upper + h_lower + b, and first layer total height +# H_lower = h_lower + b as the first primitive variable for easier visualization and setting initial +# conditions +function Trixi.varnames(::typeof(cons2prim), ::ShallowWaterTwoLayerEquations2D) + ("H_upper", "v1_upper", "v2_upper", "H_lower", "v1_lower", "v2_lower", "b") +end + +# Set initial conditions at physical location `x` for time `t` +""" + initial_condition_convergence_test(x, t, equations::ShallowWaterTwoLayerEquations2D) + +A smooth initial condition used for convergence tests in combination with +[`source_terms_convergence_test`](@ref). Constants must be set to ``rho_{upper} = 0.9``, +``rho_{lower} = 1.0``, ``g = 10.0``. +""" +function Trixi.initial_condition_convergence_test(x, t, + equations::ShallowWaterTwoLayerEquations2D) + # some constants are chosen such that the function is periodic on the domain [0,sqrt(2)]^2] + ω = 2.0 * pi * sqrt(2.0) + + H_lower = 2.0 + 0.1 * sin(ω * x[1] + t) * cos(ω * x[2] + t) + H_upper = 4.0 + 0.1 * cos(ω * x[1] + t) * sin(ω * x[2] + t) + v1_lower = 1.0 + v1_upper = 0.9 + v2_lower = 0.9 + v2_upper = 1.0 + b = 1.0 + 0.1 * cos(0.5 * ω * x[1]) * sin(0.5 * ω * x[2]) + + return prim2cons(SVector(H_upper, v1_upper, v2_upper, H_lower, v1_lower, v2_lower, + b), equations) +end + +""" + source_terms_convergence_test(u, x, t, equations::ShallowWaterTwoLayerEquations2D) + +Source terms used for convergence tests in combination with +[`initial_condition_convergence_test`](@ref). +""" +@inline function Trixi.source_terms_convergence_test(u, x, t, + equations::ShallowWaterTwoLayerEquations2D) + # Same settings as in `initial_condition_convergence_test`. + # some constants are chosen such that the function is periodic on the domain [0,sqrt(2)]^2] + ω = 2.0 * pi * sqrt(2.0) + + # Source terms obtained with SymPy + du1 = 0.01 * ω * cos(t + ω * x[1]) * cos(t + ω * x[2]) + + 0.01 * ω * sin(t + ω * x[1]) * sin(t + ω * x[2]) + du2 = (5.0 * + (-0.1 * ω * cos(t + ω * x[1]) * cos(t + ω * x[2]) - + 0.1 * ω * sin(t + ω * x[1]) * sin(t + + ω * x[2])) * + (4.0 + 0.2cos(t + ω * x[1]) * sin(t + ω * x[2]) - + 0.2 * sin(t + ω * x[1]) * cos(t + + ω * x[2])) + + 0.009 * ω * cos(t + ω * x[1]) * cos(t + ω * x[2]) + + 0.009 * ω * sin(t + ω * x[1]) * sin(t + + ω * x[2]) + + 0.1 * ω * + (20.0 + cos(t + ω * x[1]) * sin(t + ω * x[2]) - + sin(t + ω * x[1]) * cos(t + + ω * x[2])) * cos(t + ω * x[1]) * cos(t + ω * x[2])) + du3 = (5.0 * + (0.1 * ω * cos(t + ω * x[1]) * cos(t + ω * x[2]) + + 0.1 * ω * sin(t + ω * x[1]) * sin(t + + ω * x[2])) * + (4.0 + 0.2 * cos(t + ω * x[1]) * sin(t + ω * x[2]) - + 0.2 * sin(t + ω * x[1]) * cos(t + + ω * x[2])) + + 0.01ω * cos(t + ω * x[1]) * cos(t + ω * x[2]) + + 0.01 * ω * sin(t + ω * x[1]) * sin(t + ω * x[2]) + + -0.1 * ω * + (20.0 + cos(t + ω * x[1]) * sin(t + ω * x[2]) - + sin(t + ω * x[1]) * cos(t + ω * x[2])) * sin(t + + ω * x[1]) * sin(t + ω * x[2])) + du4 = (0.1 * cos(t + ω * x[1]) * cos(t + ω * x[2]) + + 0.1 * ω * cos(t + ω * x[1]) * cos(t + ω * x[2]) + + 0.05 * ω * sin(0.5 * ω * x[1]) * sin(0.5 * ω * x[2]) - + 0.1 * sin(t + ω * x[1]) * sin(t + ω * x[2]) + + -0.045 * ω * cos(0.5 * ω * x[1]) * cos(0.5 * ω * x[2]) - + 0.09 * ω * sin(t + ω * x[1]) * sin(t + ω * x[2])) + du5 = ((10.0 + sin(t + ω * x[1]) * cos(t + ω * x[2]) - + cos(0.5 * ω * x[1]) * sin(0.5 * ω * x[2])) * (-0.09 * ω * cos(t + + ω * x[1]) * cos(t + ω * x[2]) - + 0.09 * ω * sin(t + ω * x[1]) * sin(t + ω * x[2]) + + -0.05 * ω * sin(0.5 * ω * x[1]) * sin(0.5 * ω * x[2])) + + 5.0 * + (0.1 * ω * cos(t + ω * x[1]) * cos(t + ω * x[2]) + + 0.05 * ω * sin(0.5 * ω * x[1]) * sin(0.5 * ω * x[2])) * + (2.0 + 0.2 * sin(t + ω * x[1]) * cos(t + ω * x[2]) + + -0.2 * cos(0.5 * ω * x[1]) * sin(0.5 * ω * x[2])) + + 0.1 * cos(t + ω * x[1]) * cos(t + ω * x[2]) + + 0.1 * ω * cos(t + + ω * x[1]) * cos(t + ω * x[2]) + + 0.05 * ω * sin(0.5 * ω * x[1]) * sin(0.5 * ω * x[2]) - + 0.1 * sin(t + + ω * x[1]) * sin(t + ω * x[2]) - + 0.045 * ω * cos(0.5 * ω * x[1]) * cos(0.5 * ω * x[2]) - + 0.09 * ω * sin(t + + ω * x[1]) * sin(t + ω * x[2])) + du6 = ((10.0 + sin(t + ω * x[1]) * cos(t + ω * x[2]) + + -cos(0.5 * ω * x[1]) * sin(0.5 * ω * x[2])) * + (0.05 * ω * cos(0.5 * ω * x[1]) * cos(0.5 * ω * x[2]) + + 0.09 * ω * cos(t + ω * x[1]) * cos(t + ω * x[2]) + + 0.09 * ω * sin(t + ω * x[1]) * sin(t + ω * x[2])) + + 5.0 * + (-0.05 * ω * cos(0.5 * ω * x[1]) * cos(0.5 * ω * x[2]) - + 0.1 * ω * sin(t + ω * x[1]) * sin(t + + ω * x[2])) * + (2.0 + 0.2 * sin(t + ω * x[1]) * cos(t + ω * x[2]) + + -0.2 * cos(0.5 * ω * x[1]) * sin(0.5 * ω * x[2])) + + 0.09cos(t + ω * x[1]) * cos(t + ω * x[2]) + + 0.09 * ω * cos(t + ω * x[1]) * cos(t + ω * x[2]) + + 0.045 * ω * sin(0.5 * ω * x[1]) * sin(0.5 * ω * x[2]) + + -0.09 * sin(t + ω * x[1]) * sin(t + ω * x[2]) - + 0.0405 * ω * cos(0.5 * ω * x[1]) * cos(0.5 * ω * x[2]) + + -0.081 * ω * sin(t + ω * x[1]) * sin(t + ω * x[2])) + + return SVector(du1, du2, du3, du4, du5, du6, zero(eltype(u))) +end + +""" + boundary_condition_slip_wall(u_inner, normal_direction, x, t, surface_flux_function, + equations::ShallowWaterTwoLayerEquations2D) + +Create a boundary state by reflecting the normal velocity component and keep +the tangential velocity component unchanged. The boundary water height is taken from +the internal value. + +For details see Section 9.2.5 of the book: +- Eleuterio F. Toro (2001) + Shock-Capturing Methods for Free-Surface Shallow Flows + 1st edition + ISBN 0471987662 +""" +@inline function Trixi.boundary_condition_slip_wall(u_inner, + normal_direction::AbstractVector, + x, t, surface_flux_function, + equations::ShallowWaterTwoLayerEquations2D) + # normalize the outward pointing direction + normal = normal_direction / norm(normal_direction) + + # compute the normal velocity + v_normal_upper = normal[1] * u_inner[2] + normal[2] * u_inner[3] + v_normal_lower = normal[1] * u_inner[5] + normal[2] * u_inner[6] + + # create the "external" boundary solution state + u_boundary = SVector(u_inner[1], + u_inner[2] - 2.0 * v_normal_upper * normal[1], + u_inner[3] - 2.0 * v_normal_upper * normal[2], + u_inner[4], + u_inner[5] - 2.0 * v_normal_lower * normal[1], + u_inner[6] - 2.0 * v_normal_lower * normal[2], + u_inner[7]) + + # calculate the boundary flux + flux = surface_flux_function(u_inner, u_boundary, normal_direction, equations) + return flux +end + +# Calculate 1D flux for a single point +# Note, the bottom topography has no flux +@inline function Trixi.flux(u, orientation::Integer, + equations::ShallowWaterTwoLayerEquations2D) + h_upper, h_v1_upper, h_v2_upper, h_lower, h_v1_lower, h_v2_lower, _ = u + + # Calculate velocities + v1_upper, v2_upper, v1_lower, v2_lower = velocity(u, equations) + + # Calculate pressure + p_upper = 0.5 * equations.gravity * h_upper^2 + p_lower = 0.5 * equations.gravity * h_lower^2 + + # Calculate fluxes depending on orientation + if orientation == 1 + f1 = h_v1_upper + f2 = h_v1_upper * v1_upper + p_upper + f3 = h_v1_upper * v2_upper + f4 = h_v1_lower + f5 = h_v1_lower * v1_lower + p_lower + f6 = h_v1_lower * v2_lower + else + f1 = h_v2_upper + f2 = h_v2_upper * v1_upper + f3 = h_v2_upper * v2_upper + p_upper + f4 = h_v2_lower + f5 = h_v2_lower * v1_lower + f6 = h_v2_lower * v2_lower + p_lower + end + return SVector(f1, f2, f3, f4, f5, f6, zero(eltype(u))) +end + +# Calculate 1D flux for a single point in the normal direction +# Note, this directional vector is not normalized and the bottom topography has no flux +@inline function Trixi.flux(u, normal_direction::AbstractVector, + equations::ShallowWaterTwoLayerEquations2D) + h_upper, h_lower = waterheight(u, equations) + v1_upper, v2_upper, v1_lower, v2_lower = velocity(u, equations) + + v_normal_upper = v1_upper * normal_direction[1] + v2_upper * normal_direction[2] + v_normal_lower = v1_lower * normal_direction[1] + v2_lower * normal_direction[2] + h_v_upper_normal = h_upper * v_normal_upper + h_v_lower_normal = h_lower * v_normal_lower + + p_upper = 0.5 * equations.gravity * h_upper^2 + p_lower = 0.5 * equations.gravity * h_lower^2 + + f1 = h_v_upper_normal + f2 = h_v_upper_normal * v1_upper + p_upper * normal_direction[1] + f3 = h_v_upper_normal * v2_upper + p_upper * normal_direction[2] + f4 = h_v_lower_normal + f5 = h_v_lower_normal * v1_lower + p_lower * normal_direction[1] + f6 = h_v_lower_normal * v2_lower + p_lower * normal_direction[2] + + return SVector(f1, f2, f3, f4, f5, f6, zero(eltype(u))) +end + +""" + flux_nonconservative_ersing_etal(u_ll, u_rr, orientation::Integer, + equations::ShallowWaterTwoLayerEquations2D) + flux_nonconservative_ersing_etal(u_ll, u_rr, + normal_direction_ll::AbstractVector, + normal_direction_average::AbstractVector, + equations::ShallowWaterTwoLayerEquations2D) + +!!! warning "Experimental code" + This numerical flux is experimental and may change in any future release. + +Non-symmetric path-conservative two-point volume flux discretizing the nonconservative (source) term +that contains the gradient of the bottom topography [`ShallowWaterTwoLayerEquations2D`](@ref) and an +additional term that couples the momentum of both layers. + +This is a modified version of [`flux_nonconservative_wintermeyer_etal`](@ref) that gives entropy +conservation and well-balancedness in both the volume and surface when combined with +[`flux_wintermeyer_etal`](@ref). + +For further details see: +- Patrick Ersing, Andrew R. Winters (2023) + An entropy stable discontinuous Galerkin method for the two-layer shallow water equations on + curvilinear meshes + [DOI: 10.48550/arXiv.2306.12699](https://doi.org/10.48550/arXiv.2306.12699) +""" +@inline function Trixi.flux_nonconservative_ersing_etal(u_ll, u_rr, + orientation::Integer, + equations::ShallowWaterTwoLayerEquations2D) + # Pull the necessary left and right state information + h_upper_ll, h_lower_ll = waterheight(u_ll, equations) + h_upper_rr, h_lower_rr = waterheight(u_rr, equations) + b_rr = u_rr[7] + b_ll = u_ll[7] + + # Calculate jumps + h_upper_jump = (h_upper_rr - h_upper_ll) + h_lower_jump = (h_lower_rr - h_lower_ll) + b_jump = (b_rr - b_ll) + + z = zero(eltype(u_ll)) + + # Bottom gradient nonconservative term: (0, g*h_upper*(b + h_lower)_x, g*h_upper*(b + h_lower)_y , + # 0, g*h_lower*(b + r*h_upper)_x, + # g*h_lower*(b + r*h_upper)_y, 0) + if orientation == 1 + f = SVector(z, + equations.gravity * h_upper_ll * (b_jump + h_lower_jump), + z, z, + equations.gravity * h_lower_ll * + (b_jump + equations.r * h_upper_jump), + z, z) + else # orientation == 2 + f = SVector(z, z, + equations.gravity * h_upper_ll * (b_jump + h_lower_jump), + z, z, + equations.gravity * h_lower_ll * + (b_jump + equations.r * h_upper_jump), + z) + end + + return f +end + +@inline function Trixi.flux_nonconservative_ersing_etal(u_ll, u_rr, + normal_direction_ll::AbstractVector, + normal_direction_average::AbstractVector, + equations::ShallowWaterTwoLayerEquations2D) + # Pull the necessary left and right state information + h_upper_ll, h_lower_ll = waterheight(u_ll, equations) + h_upper_rr, h_lower_rr = waterheight(u_rr, equations) + b_rr = u_rr[7] + b_ll = u_ll[7] + + # Calculate jumps + h_upper_jump = (h_upper_rr - h_upper_ll) + h_lower_jump = (h_lower_rr - h_lower_ll) + b_jump = (b_rr - b_ll) + + # Note this routine only uses the `normal_direction_average` and the average of the + # bottom topography to get a quadratic split form DG gradient on curved elements + return SVector(zero(eltype(u_ll)), + normal_direction_average[1] * equations.gravity * h_upper_ll * + (b_jump + h_lower_jump), + normal_direction_average[2] * equations.gravity * h_upper_ll * + (b_jump + h_lower_jump), + zero(eltype(u_ll)), + normal_direction_average[1] * equations.gravity * h_lower_ll * + (b_jump + equations.r * h_upper_jump), + normal_direction_average[2] * equations.gravity * h_lower_ll * + (b_jump + equations.r * h_upper_jump), + zero(eltype(u_ll))) +end + +""" + flux_wintermeyer_etal(u_ll, u_rr, orientation, + equations::ShallowWaterTwoLayerEquations2D) + flux_wintermeyer_etal(u_ll, u_rr, + normal_direction::AbstractVector, + equations::ShallowWaterTwoLayerEquations2D) + +Total energy conservative (mathematical entropy for two-layer shallow water equations) split form. +When the bottom topography is nonzero this scheme will be well-balanced when used with the +nonconservative [`flux_nonconservative_ersing_etal`](@ref). To obtain the flux for the +two-layer shallow water equations the flux that is described in the paper for the normal shallow +water equations is used within each layer. + +Further details are available in Theorem 1 of the paper: +- Niklas Wintermeyer, Andrew R. Winters, Gregor J. Gassner and David A. Kopriva (2017) + An entropy stable nodal discontinuous Galerkin method for the two dimensional + shallow water equations on unstructured curvilinear meshes with discontinuous bathymetry + [DOI: 10.1016/j.jcp.2017.03.036](https://doi.org/10.1016/j.jcp.2017.03.036) +""" +@inline function Trixi.flux_wintermeyer_etal(u_ll, u_rr, + orientation::Integer, + equations::ShallowWaterTwoLayerEquations2D) + # Unpack left and right state + h_upper_ll, h_v1_upper_ll, h_v2_upper_ll, h_lower_ll, h_v1_lower_ll, h_v2_lower_ll, _ = u_ll + h_upper_rr, h_v1_upper_rr, h_v2_upper_rr, h_lower_rr, h_v1_lower_rr, h_v2_lower_rr, _ = u_rr + + # Get the velocities on either side + v1_upper_ll, v2_upper_ll, v1_lower_ll, v2_lower_ll = velocity(u_ll, equations) + v1_upper_rr, v2_upper_rr, v1_lower_rr, v2_lower_rr = velocity(u_rr, equations) + + # Average each factor of products in flux + v1_upper_avg = 0.5 * (v1_upper_ll + v1_upper_rr) + v1_lower_avg = 0.5 * (v1_lower_ll + v1_lower_rr) + v2_upper_avg = 0.5 * (v2_upper_ll + v2_upper_rr) + v2_lower_avg = 0.5 * (v2_lower_ll + v2_lower_rr) + p_upper_avg = 0.5 * equations.gravity * h_upper_ll * h_upper_rr + p_lower_avg = 0.5 * equations.gravity * h_lower_ll * h_lower_rr + + # Calculate fluxes depending on orientation + if orientation == 1 + f1 = 0.5 * (h_v1_upper_ll + h_v1_upper_rr) + f2 = f1 * v1_upper_avg + p_upper_avg + f3 = f1 * v2_upper_avg + f4 = 0.5 * (h_v1_lower_ll + h_v1_lower_rr) + f5 = f4 * v1_lower_avg + p_lower_avg + f6 = f4 * v2_lower_avg + else + f1 = 0.5 * (h_v2_upper_ll + h_v2_upper_rr) + f2 = f1 * v1_upper_avg + f3 = f1 * v2_upper_avg + p_upper_avg + f4 = 0.5 * (h_v2_lower_ll + h_v2_lower_rr) + f5 = f4 * v1_lower_avg + f6 = f4 * v2_lower_avg + p_lower_avg + end + + return SVector(f1, f2, f3, f4, f5, f6, zero(eltype(u_ll))) +end + +@inline function Trixi.flux_wintermeyer_etal(u_ll, u_rr, + normal_direction::AbstractVector, + equations::ShallowWaterTwoLayerEquations2D) + # Unpack left and right state + h_upper_ll, h_v1_upper_ll, h_v2_upper_ll, h_lower_ll, h_v1_lower_ll, h_v2_lower_ll, _ = u_ll + h_upper_rr, h_v1_upper_rr, h_v2_upper_rr, h_lower_rr, h_v1_lower_rr, h_v2_lower_rr, _ = u_rr + + # Get the velocities on either side + v1_upper_ll, v2_upper_ll, v1_lower_ll, v2_lower_ll = velocity(u_ll, equations) + v1_upper_rr, v2_upper_rr, v1_lower_rr, v2_lower_rr = velocity(u_rr, equations) + + # Average each factor of products in flux + v1_upper_avg = 0.5 * (v1_upper_ll + v1_upper_rr) + v1_lower_avg = 0.5 * (v1_lower_ll + v1_lower_rr) + v2_upper_avg = 0.5 * (v2_upper_ll + v2_upper_rr) + v2_lower_avg = 0.5 * (v2_lower_ll + v2_lower_rr) + p_upper_avg = 0.5 * equations.gravity * h_upper_ll * h_upper_rr + p_lower_avg = 0.5 * equations.gravity * h_lower_ll * h_lower_rr + h_v1_upper_avg = 0.5 * (h_v1_upper_ll + h_v1_upper_rr) + h_v2_upper_avg = 0.5 * (h_v2_upper_ll + h_v2_upper_rr) + h_v1_lower_avg = 0.5 * (h_v1_lower_ll + h_v1_lower_rr) + h_v2_lower_avg = 0.5 * (h_v2_lower_ll + h_v2_lower_rr) + + # Calculate fluxes depending on normal_direction + f1 = h_v1_upper_avg * normal_direction[1] + h_v2_upper_avg * normal_direction[2] + f2 = f1 * v1_upper_avg + p_upper_avg * normal_direction[1] + f3 = f1 * v2_upper_avg + p_upper_avg * normal_direction[2] + f4 = h_v1_lower_avg * normal_direction[1] + h_v2_lower_avg * normal_direction[2] + f5 = f4 * v1_lower_avg + p_lower_avg * normal_direction[1] + f6 = f4 * v2_lower_avg + p_lower_avg * normal_direction[2] + + return SVector(f1, f2, f3, f4, f5, f6, zero(eltype(u_ll))) +end + +""" + flux_es_ersing_etal(u_ll, u_rr, orientation_or_normal_direction, + equations::ShallowWaterTwoLayerEquations2D) + +Entropy stable surface flux for the two-layer shallow water equations. Uses the entropy conservative +[`flux_wintermeyer_etal`](@ref) and adds a Lax-Friedrichs type dissipation dependent on the jump of +entropy variables. + +For further details see: +- Patrick Ersing, Andrew R. Winters (2023) + An entropy stable discontinuous Galerkin method for the two-layer shallow water equations on + curvilinear meshes + [DOI: 10.48550/arXiv.2306.12699](https://doi.org/10.48550/arXiv.2306.12699) +""" +@inline function flux_es_ersing_etal(u_ll, u_rr, + orientation_or_normal_direction, + equations::ShallowWaterTwoLayerEquations2D) + # Compute entropy conservative flux but without the bottom topography + f_ec = flux_wintermeyer_etal(u_ll, u_rr, + orientation_or_normal_direction, + equations) + + # Get maximum signal velocity + λ = max_abs_speed_naive(u_ll, u_rr, orientation_or_normal_direction, equations) + + # Get entropy variables but without the bottom topography + q_rr = cons2entropy(u_rr, equations) + q_ll = cons2entropy(u_ll, equations) + + # Average values from left and right + u_avg = (u_ll + u_rr) / 2 + + # Introduce variables for better readability + rho_upper = equations.rho_upper + rho_lower = equations.rho_lower + g = equations.gravity + drho = rho_upper - rho_lower + + # Compute entropy Jacobian coefficients + h11 = -rho_lower / (g * rho_upper * drho) + h12 = -rho_lower * u_avg[2] / (g * rho_upper * u_avg[1] * drho) + h13 = -rho_lower * u_avg[3] / (g * rho_upper * u_avg[1] * drho) + h14 = 1.0 / (g * drho) + h15 = u_avg[5] / (g * u_avg[4] * drho) + h16 = u_avg[6] / (g * u_avg[4] * drho) + h21 = -rho_lower * u_avg[2] / (g * rho_upper * u_avg[1] * drho) + h22 = ((g * rho_upper * u_avg[1]^3 - g * rho_lower * u_avg[1]^3 + + -rho_lower * u_avg[2]^2) / (g * rho_upper * u_avg[1]^2 * drho)) + h23 = -rho_lower * u_avg[2] * u_avg[3] / (g * rho_upper * u_avg[1]^2 * drho) + h24 = u_avg[2] / (g * u_avg[1] * drho) + h25 = u_avg[2] * u_avg[5] / (g * u_avg[1] * u_avg[4] * drho) + h26 = u_avg[2] * u_avg[6] / (g * u_avg[1] * u_avg[4] * drho) + h31 = -rho_lower * u_avg[3] / (g * rho_upper * u_avg[1] * drho) + h32 = -rho_lower * u_avg[2] * u_avg[3] / (g * rho_upper * u_avg[1]^2 * drho) + h33 = ((g * rho_upper * u_avg[1]^3 - g * rho_lower * u_avg[1]^3 + + -rho_lower * u_avg[3]^2) / (g * rho_upper * u_avg[1]^2 * drho)) + h34 = u_avg[3] / (g * u_avg[1] * drho) + h35 = u_avg[3] * u_avg[5] / (g * u_avg[1] * u_avg[4] * drho) + h36 = u_avg[3] * u_avg[6] / (g * u_avg[1] * u_avg[4] * drho) + h41 = 1.0 / (g * drho) + h42 = u_avg[2] / (g * u_avg[1] * drho) + h43 = u_avg[3] / (g * u_avg[1] * drho) + h44 = -1.0 / (g * drho) + h45 = -u_avg[5] / (g * u_avg[4] * drho) + h46 = -u_avg[6] / (g * u_avg[4] * drho) + h51 = u_avg[5] / (g * u_avg[4] * drho) + h52 = u_avg[2] * u_avg[5] / (g * u_avg[1] * u_avg[4] * drho) + h53 = u_avg[3] * u_avg[5] / (g * u_avg[1] * u_avg[4] * drho) + h54 = -u_avg[5] / (g * u_avg[4] * drho) + h55 = ((g * rho_upper * u_avg[4]^3 - g * rho_lower * u_avg[4]^3 + + -rho_lower * u_avg[5]^2) / (g * rho_lower * u_avg[4]^2 * drho)) + h56 = -u_avg[5] * u_avg[6] / (g * u_avg[4]^2 * drho) + h61 = u_avg[6] / (g * u_avg[4] * drho) + h62 = u_avg[2] * u_avg[6] / (g * u_avg[1] * u_avg[4] * drho) + h63 = u_avg[3] * u_avg[6] / (g * u_avg[1] * u_avg[4] * drho) + h64 = -u_avg[6] / (g * u_avg[4] * drho) + h65 = -u_avg[5] * u_avg[6] / (g * u_avg[4]^2 * drho) + h66 = ((g * rho_upper * u_avg[4]^3 - g * rho_lower * u_avg[4]^3 + + -rho_lower * u_avg[6]^2) / (g * rho_lower * u_avg[4]^2 * drho)) + + # Entropy Jacobian matrix + H = @SMatrix [[h11;; h12;; h13;; h14;; h15;; h16;; 0]; + [h21;; h22;; h23;; h24;; h25;; h26;; 0]; + [h31;; h32;; h33;; h34;; h35;; h36;; 0]; + [h41;; h42;; h43;; h44;; h45;; h46;; 0]; + [h51;; h52;; h53;; h54;; h55;; h56;; 0]; + [h61;; h62;; h63;; h64;; h65;; h66;; 0]; + [0;; 0;; 0;; 0;; 0;; 0;; 0]] + + # Add dissipation to entropy conservative flux to obtain entropy stable flux + f_es = f_ec - 0.5 * λ * H * (q_rr - q_ll) + + return SVector(f_es[1], f_es[2], f_es[3], f_es[4], f_es[5], f_es[6], + zero(eltype(u_ll))) +end + +# Calculate approximation for maximum wave speed for local Lax-Friedrichs-type dissipation as the +# maximum velocity magnitude plus the maximum speed of sound. This function uses approximate +# eigenvalues using the speed of the barotropic mode as there is no simple way to calculate them +# analytically. +# +# A good overview of the derivation is given in: +# - Jonas Nycander, Andrew McC. Hogg, Leela M. Frankcombe (2008) +# Open boundary conditions for nonlinear channel Flows +# [DOI: 10.1016/j.ocemod.2008.06.003](https://doi.org/10.1016/j.ocemod.2008.06.003) +@inline function Trixi.max_abs_speed_naive(u_ll, u_rr, + orientation::Integer, + equations::ShallowWaterTwoLayerEquations2D) + # Unpack left and right state + h_upper_ll, h_v1_upper_ll, h_v2_upper_ll, h_lower_ll, h_v1_lower_ll, h_v2_lower_ll, _ = u_ll + h_upper_rr, h_v1_upper_rr, h_v2_upper_rr, h_lower_rr, h_v1_lower_rr, h_v2_lower_rr, _ = u_rr + + # Calculate averaged velocity of both layers + if orientation == 1 + v_m_ll = (h_v1_upper_ll + h_v1_lower_ll) / (h_upper_ll + h_lower_ll) + v_m_rr = (h_v1_upper_rr + h_v1_lower_rr) / (h_upper_rr + h_lower_rr) + else + v_m_ll = (h_v2_upper_ll + h_v2_lower_ll) / (h_upper_ll + h_lower_ll) + v_m_rr = (h_v2_upper_rr + h_v2_lower_rr) / (h_upper_rr + h_lower_rr) + end + + # Calculate the wave celerity on the left and right + h_upper_ll, h_lower_ll = waterheight(u_ll, equations) + h_upper_rr, h_lower_rr = waterheight(u_rr, equations) + + c_ll = sqrt(equations.gravity * (h_upper_ll + h_lower_ll)) + c_rr = sqrt(equations.gravity * (h_upper_rr + h_lower_rr)) + + return (max(abs(v_m_ll), abs(v_m_rr)) + max(c_ll, c_rr)) +end + +@inline function Trixi.max_abs_speed_naive(u_ll, u_rr, + normal_direction::AbstractVector, + equations::ShallowWaterTwoLayerEquations2D) + # Unpack left and right state + h_upper_ll, _, _, h_lower_ll, _, _, _ = u_ll + h_upper_rr, _, _, h_lower_rr, _, _, _ = u_rr + + # Extract and compute the velocities in the normal direction + v1_upper_ll, v2_upper_ll, v1_lower_ll, v2_lower_ll = velocity(u_ll, equations) + v1_upper_rr, v2_upper_rr, v1_lower_rr, v2_lower_rr = velocity(u_rr, equations) + + v_upper_dot_n_ll = v1_upper_ll * normal_direction[1] + + v2_upper_ll * normal_direction[2] + v_upper_dot_n_rr = v1_upper_rr * normal_direction[1] + + v2_upper_rr * normal_direction[2] + v_lower_dot_n_ll = v1_lower_ll * normal_direction[1] + + v2_lower_ll * normal_direction[2] + v_lower_dot_n_rr = v1_lower_rr * normal_direction[1] + + v2_lower_rr * normal_direction[2] + + # Calculate averaged velocity of both layers + v_m_ll = (v_upper_dot_n_ll * h_upper_ll + v_lower_dot_n_ll * h_lower_ll) / + (h_upper_ll + h_lower_ll) + v_m_rr = (v_upper_dot_n_rr * h_upper_rr + v_lower_dot_n_rr * h_lower_rr) / + (h_upper_rr + h_lower_rr) + + # Compute the wave celerity on the left and right + h_upper_ll, h_lower_ll = waterheight(u_ll, equations) + h_upper_rr, h_lower_rr = waterheight(u_rr, equations) + + c_ll = sqrt(equations.gravity * (h_upper_ll + h_lower_ll)) + c_rr = sqrt(equations.gravity * (h_upper_rr + h_lower_rr)) + + # The normal velocities are already scaled by the norm + return max(abs(v_m_ll), abs(v_m_rr)) + max(c_ll, c_rr) * norm(normal_direction) +end + +# Specialized `DissipationLocalLaxFriedrichs` to avoid spurious dissipation in the bottom topography +@inline function (dissipation::DissipationLocalLaxFriedrichs)(u_ll, u_rr, + orientation_or_normal_direction, + equations::ShallowWaterTwoLayerEquations2D) + λ = dissipation.max_abs_speed(u_ll, u_rr, orientation_or_normal_direction, + equations) + diss = -0.5 * λ * (u_rr - u_ll) + return SVector(diss[1], diss[2], diss[3], diss[4], diss[5], diss[6], + zero(eltype(u_ll))) +end + +# Absolute speed of the barotropic mode +@inline function Trixi.max_abs_speeds(u, equations::ShallowWaterTwoLayerEquations2D) + h_upper, h_v1_upper, h_v2_upper, h_lower, h_v1_lower, h_v2_lower, _ = u + + # Calculate averaged velocity of both layers + v1_m = (h_v1_upper + h_v1_lower) / (h_upper + h_lower) + v2_m = (h_v2_upper + h_v2_lower) / (h_upper + h_lower) + + h_upper, h_lower = waterheight(u, equations) + v1_upper, v2_upper, v1_lower, v2_lower = velocity(u, equations) + + c = sqrt(equations.gravity * (h_upper + h_lower)) + return (max(abs(v1_m) + c, abs(v1_upper), abs(v1_lower)), + max(abs(v2_m) + c, abs(v2_upper), abs(v2_lower))) +end + +# Helper function to extract the velocity vector from the conservative variables +@inline function Trixi.velocity(u, equations::ShallowWaterTwoLayerEquations2D) + h_upper, h_v1_upper, h_v2_upper, h_lower, h_v1_lower, h_v2_lower, _ = u + + v1_upper = h_v1_upper / h_upper + v2_upper = h_v2_upper / h_upper + v1_lower = h_v1_lower / h_lower + v2_lower = h_v2_lower / h_lower + + return SVector(v1_upper, v2_upper, v1_lower, v2_lower) +end + +# Convert conservative variables to primitive +@inline function Trixi.cons2prim(u, equations::ShallowWaterTwoLayerEquations2D) + h_upper, _, _, h_lower, _, _, b = u + + H_lower = h_lower + b + H_upper = h_lower + h_upper + b + v1_upper, v2_upper, v1_lower, v2_lower = velocity(u, equations) + + return SVector(H_upper, v1_upper, v2_upper, H_lower, v1_lower, v2_lower, b) +end + +# Convert conservative variables to entropy variables +# Note, only the first four are the entropy variables, the fifth entry still just carries the bottom +# topography values for convenience. +# In contrast to general usage the entropy variables are denoted with q instead of w, because w is +# already used for velocity in y-Direction +@inline function Trixi.cons2entropy(u, equations::ShallowWaterTwoLayerEquations2D) + h_upper, _, _, h_lower, _, _, b = u + # Assign new variables for better readability + rho_upper = equations.rho_upper + rho_lower = equations.rho_lower + v1_upper, v2_upper, v1_lower, v2_lower = velocity(u, equations) + + w1 = (rho_upper * (equations.gravity * (h_upper + h_lower + b) + + -0.5 * (v1_upper^2 + v2_upper^2))) + w2 = rho_upper * v1_upper + w3 = rho_upper * v2_upper + w4 = (rho_lower * (equations.gravity * (equations.r * h_upper + h_lower + b) + + -0.5 * (v1_lower^2 + v2_lower^2))) + w5 = rho_lower * v1_lower + w6 = rho_lower * v2_lower + return SVector(w1, w2, w3, w4, w5, w6, b) +end + +# Convert primitive to conservative variables +@inline function Trixi.prim2cons(prim, equations::ShallowWaterTwoLayerEquations2D) + H_upper, v1_upper, v2_upper, H_lower, v1_lower, v2_lower, b = prim + + h_lower = H_lower - b + h_upper = H_upper - h_lower - b + h_v1_upper = h_upper * v1_upper + h_v2_upper = h_upper * v2_upper + h_v1_lower = h_lower * v1_lower + h_v2_lower = h_lower * v2_lower + return SVector(h_upper, h_v1_upper, h_v2_upper, h_lower, h_v1_lower, h_v2_lower, b) +end + +@inline function Trixi.waterheight(u, equations::ShallowWaterTwoLayerEquations2D) + return SVector(u[1], u[4]) +end + +# Entropy function for the shallow water equations is the total energy +@inline function Trixi.entropy(cons, equations::ShallowWaterTwoLayerEquations2D) + energy_total(cons, equations) +end + +# Calculate total energy for a conservative state `cons` +@inline function Trixi.energy_total(cons, equations::ShallowWaterTwoLayerEquations2D) + h_upper, h_v1_upper, h_v2_upper, h_lower, h_v2_lower, h_v2_lower, b = cons + g = equations.gravity + rho_upper = equations.rho_upper + rho_lower = equations.rho_lower + + e = (0.5 * rho_upper * + (h_v1_upper^2 / h_upper + h_v2_upper^2 / h_upper + g * h_upper^2) + + 0.5 * rho_lower * + (h_v2_lower^2 / h_lower + h_v2_lower^2 / h_lower + g * h_lower^2) + + g * rho_lower * h_lower * b + g * rho_upper * h_upper * (h_lower + b)) + return e +end + +# Calculate kinetic energy for a conservative state `cons` +@inline function Trixi.energy_kinetic(u, equations::ShallowWaterTwoLayerEquations2D) + h_upper, h_v1_upper, h_v2_upper, h_lower, h_v2_lower, h_v2_lower, _ = u + + return (0.5 * equations.rho_upper * h_v1_upper^2 / h_upper + + 0.5 * equations.rho_upper * h_v2_upper^2 / h_upper + + 0.5 * equations.rho_lower * h_v2_lower^2 / h_lower + + 0.5 * equations.rho_lower * h_v2_lower^2 / h_lower) +end + +# Calculate potential energy for a conservative state `cons` +@inline function Trixi.energy_internal(cons, equations::ShallowWaterTwoLayerEquations2D) + return energy_total(cons, equations) - energy_kinetic(cons, equations) +end + +# Calculate the error for the "lake-at-rest" test case where H = h_upper+h_lower+b should +# be a constant value over time +@inline function Trixi.lake_at_rest_error(u, equations::ShallowWaterTwoLayerEquations2D) + h_upper, _, _, h_lower, _, _, b = u + return abs(equations.H0 - (h_upper + h_lower + b)) +end +end # @muladd diff --git a/src/equations/shallow_water_wet_dry_1d.jl b/src/equations/shallow_water_wet_dry_1d.jl index 095041d..f3eaa69 100644 --- a/src/equations/shallow_water_wet_dry_1d.jl +++ b/src/equations/shallow_water_wet_dry_1d.jl @@ -513,8 +513,8 @@ Further details on this hydrostatic reconstruction and its motivation can be fou v_rr = velocity(u_rr, equations) # Calculate the wave celerity on the left and right - h_ll = Trixi.waterheight(u_ll, equations) - h_rr = Trixi.waterheight(u_rr, equations) + h_ll = waterheight(u_ll, equations) + h_rr = waterheight(u_rr, equations) a_ll = sqrt(equations.gravity * h_ll) a_rr = sqrt(equations.gravity * h_rr) @@ -583,8 +583,8 @@ end end @inline function Trixi.waterheight(u, equations::ShallowWaterEquationsWetDry1D) - return Trixi.waterheight(u, - equations.basic_swe) + return waterheight(u, + equations.basic_swe) end @inline function Trixi.pressure(u, equations::ShallowWaterEquationsWetDry1D) diff --git a/src/equations/shallow_water_wet_dry_2d.jl b/src/equations/shallow_water_wet_dry_2d.jl index a9edf53..7762168 100644 --- a/src/equations/shallow_water_wet_dry_2d.jl +++ b/src/equations/shallow_water_wet_dry_2d.jl @@ -682,9 +682,9 @@ the reference below. The definition of the wave speeds are given in Equation (2. """ @inline function min_max_speed_chen_noelle(u_ll, u_rr, orientation::Integer, equations::ShallowWaterEquationsWetDry2D) - h_ll = Trixi.waterheight(u_ll, equations) + h_ll = waterheight(u_ll, equations) v1_ll, v2_ll = velocity(u_ll, equations) - h_rr = Trixi.waterheight(u_rr, equations) + h_rr = waterheight(u_rr, equations) v1_rr, v2_rr = velocity(u_rr, equations) a_ll = sqrt(equations.gravity * h_ll) @@ -703,9 +703,9 @@ end @inline function min_max_speed_chen_noelle(u_ll, u_rr, normal_direction::AbstractVector, equations::ShallowWaterEquationsWetDry2D) - h_ll = Trixi.waterheight(u_ll, equations) + h_ll = waterheight(u_ll, equations) v1_ll, v2_ll = velocity(u_ll, equations) - h_rr = Trixi.waterheight(u_rr, equations) + h_rr = waterheight(u_rr, equations) v1_rr, v2_rr = velocity(u_rr, equations) v_normal_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2] @@ -797,8 +797,8 @@ end end @inline function Trixi.waterheight(u, equations::ShallowWaterEquationsWetDry2D) - return Trixi.waterheight(u, - equations.basic_swe) + return waterheight(u, + equations.basic_swe) end @inline function Trixi.pressure(u, equations::ShallowWaterEquationsWetDry2D) diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl index 6312471..506d903 100644 --- a/test/test_structured_2d.jl +++ b/test/test_structured_2d.jl @@ -15,132 +15,138 @@ isdir(outdir) && rm(outdir, recursive = true) @testset "StructuredMesh2D" begin #! format: noindent -@trixi_testset "elixir_shallowwater_source_terms.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"), - l2=[ - 0.0017285599436729316, - 0.025584610912606776, - 0.028373834961180594, - 6.274146767730866e-5, - ], - linf=[ - 0.012972309788264802, - 0.108283714215621, - 0.15831585777928936, - 0.00018196759554722775, - ], - tspan=(0.0, 0.05)) - # 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 +@testset "Shallow Water Wet/Dry" begin + @trixi_testset "elixir_shallowwater_source_terms.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_source_terms.jl"), + l2=[ + 0.0017285599436729316, + 0.025584610912606776, + 0.028373834961180594, + 6.274146767730866e-5, + ], + linf=[ + 0.012972309788264802, + 0.108283714215621, + 0.15831585777928936, + 0.00018196759554722775, + ], + tspan=(0.0, 0.05)) + # 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 -@trixi_testset "elixir_shallowwater_well_balanced.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_well_balanced.jl"), - l2=[ - 0.7920927046419308, - 9.92129670988898e-15, - 1.0118635033124588e-14, - 0.7920927046419308, - ], - linf=[ - 2.408429868800133, - 5.5835419986809516e-14, - 5.448874313931364e-14, - 2.4084298688001335, - ], - 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 + @trixi_testset "elixir_shallowwater_well_balanced.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_well_balanced.jl"), + l2=[ + 0.7920927046419308, + 9.92129670988898e-15, + 1.0118635033124588e-14, + 0.7920927046419308, + ], + linf=[ + 2.408429868800133, + 5.5835419986809516e-14, + 5.448874313931364e-14, + 2.4084298688001335, + ], + 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 -@trixi_testset "elixir_shallowwater_well_balanced_wet_dry.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, - "elixir_shallowwater_well_balanced_wet_dry.jl"), - l2=[ - 0.019731646454942086, - 1.0694532773278277e-14, - 1.1969913383405568e-14, - 0.0771517260037954, - ], - linf=[ - 0.4999999999998892, - 6.067153702623552e-14, - 4.4849667259339357e-14, - 1.9999999999999993, - ], - 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 + @trixi_testset "elixir_shallowwater_well_balanced_wet_dry.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_well_balanced_wet_dry.jl"), + l2=[ + 0.019731646454942086, + 1.0694532773278277e-14, + 1.1969913383405568e-14, + 0.0771517260037954, + ], + linf=[ + 0.4999999999998892, + 6.067153702623552e-14, + 4.4849667259339357e-14, + 1.9999999999999993, + ], + 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 -@trixi_testset "elixir_shallowwater_conical_island.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_conical_island.jl"), - l2=[ - 0.04593154164306353, - 0.1644534881916908, - 0.16445348819169076, - 0.0011537702354532122, - ], - linf=[ - 0.21100717610846442, - 0.9501592344310412, - 0.950159234431041, - 0.021790250683516296, - ], - tspan=(0.0, 0.025)) - # 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 + @trixi_testset "elixir_shallowwater_conical_island.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_conical_island.jl"), + l2=[ + 0.04593154164306353, + 0.1644534881916908, + 0.16445348819169076, + 0.0011537702354532122, + ], + linf=[ + 0.21100717610846442, + 0.9501592344310412, + 0.950159234431041, + 0.021790250683516296, + ], + tspan=(0.0, 0.025)) + # 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 -@trixi_testset "elixir_shallowwater_parabolic_bowl.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_parabolic_bowl.jl"), - l2=[ - 0.00015285369980313484, - 1.9536806395943226e-5, - 9.936906607758672e-5, - 5.0686313334616055e-15, - ], - linf=[ - 0.003316119030459211, - 0.0005075409427972817, - 0.001986721761060583, - 4.701794509287538e-14, - ], - tspan=(0.0, 0.025), cells_per_dimension=(40, 40)) - # 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 + @trixi_testset "elixir_shallowwater_parabolic_bowl.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_parabolic_bowl.jl"), + l2=[ + 0.00015285369980313484, + 1.9536806395943226e-5, + 9.936906607758672e-5, + 5.0686313334616055e-15, + ], + linf=[ + 0.003316119030459211, + 0.0005075409427972817, + 0.001986721761060583, + 4.701794509287538e-14, + ], + tspan=(0.0, 0.025), cells_per_dimension=(40, 40)) + # 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 +end # SWE +end # StructuredMesh2D # Clean up afterwards: delete TrixiShallowWater.jl output directory @test_nowarn rm(outdir, recursive = true) diff --git a/test/test_tree_1d.jl b/test/test_tree_1d.jl index e5a1dee..b46b4f9 100644 --- a/test/test_tree_1d.jl +++ b/test/test_tree_1d.jl @@ -15,384 +15,465 @@ isdir(outdir) && rm(outdir, recursive = true) @testset "TreeMesh1D" begin #! format: noindent -@trixi_testset "elixir_shallowwater_ec.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_ec.jl"), - l2=[0.244729018751225, 0.8583565222389505, 0.07330427577586297], - linf=[ - 2.1635021283528504, - 3.8717508164234453, - 1.7711213427919539, - ], - 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 +@testset "Shallow Water Wet/Dry" begin + @trixi_testset "elixir_shallowwater_ec.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_ec.jl"), + l2=[ + 0.244729018751225, + 0.8583565222389505, + 0.07330427577586297, + ], + linf=[ + 2.1635021283528504, + 3.8717508164234453, + 1.7711213427919539, + ], + 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 -@trixi_testset "elixir_shallowwater_ec.jl with initial_condition_weak_blast_wave" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_ec.jl"), - l2=[ - 0.39464782107209717, - 2.03880864210846, - 4.1623084150546725e-10, - ], - linf=[ - 0.778905801278281, - 3.2409883402608273, - 7.419800190922032e-10, - ], - initial_condition=initial_condition_weak_blast_wave, - 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 + @trixi_testset "elixir_shallowwater_ec.jl with initial_condition_weak_blast_wave" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_ec.jl"), + l2=[ + 0.39464782107209717, + 2.03880864210846, + 4.1623084150546725e-10, + ], + linf=[ + 0.778905801278281, + 3.2409883402608273, + 7.419800190922032e-10, + ], + initial_condition=initial_condition_weak_blast_wave, + 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 -@trixi_testset "elixir_shallowwater_well_balanced.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_well_balanced.jl"), - l2=[ - 0.10416666834254829, - 1.4352935256803184e-14, - 0.10416666834254838, - ], - linf=[1.9999999999999996, 3.248036646353028e-14, 2.0], - 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 + @trixi_testset "elixir_shallowwater_well_balanced.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_well_balanced.jl"), + l2=[ + 0.10416666834254829, + 1.4352935256803184e-14, + 0.10416666834254838, + ], + linf=[1.9999999999999996, 3.248036646353028e-14, 2.0], + 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 -@trixi_testset "elixir_shallowwater_well_balanced.jl with FluxHydrostaticReconstruction" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_well_balanced.jl"), - l2=[ - 0.10416666834254835, - 1.1891029971551825e-14, - 0.10416666834254838, - ], - linf=[2.0000000000000018, 2.4019608337954543e-14, 2.0], - surface_flux=(FluxHydrostaticReconstruction(flux_lax_friedrichs, - hydrostatic_reconstruction_audusse_etal), - flux_nonconservative_audusse_etal), - 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 + @trixi_testset "elixir_shallowwater_well_balanced.jl with FluxHydrostaticReconstruction" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_well_balanced.jl"), + l2=[ + 0.10416666834254835, + 1.1891029971551825e-14, + 0.10416666834254838, + ], + linf=[2.0000000000000018, 2.4019608337954543e-14, 2.0], + surface_flux=(FluxHydrostaticReconstruction(flux_lax_friedrichs, + hydrostatic_reconstruction_audusse_etal), + flux_nonconservative_audusse_etal), + 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 -@trixi_testset "elixir_shallowwater_well_balanced.jl with flux_nonconservative_ersing_etal" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_well_balanced.jl"), - l2=[ - 0.10416666834254838, - 1.6657566141935285e-14, - 0.10416666834254838, - ], - linf=[2.0000000000000004, 3.0610625110157164e-14, 2.0], - surface_flux=(flux_wintermeyer_etal, - flux_nonconservative_ersing_etal), - volume_flux=(flux_wintermeyer_etal, - flux_nonconservative_ersing_etal), - 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 + @trixi_testset "elixir_shallowwater_well_balanced.jl with flux_nonconservative_ersing_etal" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_well_balanced.jl"), + l2=[ + 0.10416666834254838, + 1.6657566141935285e-14, + 0.10416666834254838, + ], + linf=[2.0000000000000004, 3.0610625110157164e-14, 2.0], + surface_flux=(flux_wintermeyer_etal, + flux_nonconservative_ersing_etal), + volume_flux=(flux_wintermeyer_etal, + flux_nonconservative_ersing_etal), + 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 -@trixi_testset "elixir_shallowwater_well_balanced_wet_dry.jl with FluxHydrostaticReconstruction" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, - "elixir_shallowwater_well_balanced_wet_dry.jl"), - l2=[ - 0.00965787167169024, - 5.345454081916856e-14, - 0.03857583749209928, - ], - linf=[ - 0.4999999999998892, - 2.2447689894899726e-13, - 1.9999999999999714, - ], - tspan=(0.0, 0.25), - # Soften the tolerance as test results vary between different CPUs - atol=1000 * eps()) - # 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 + @trixi_testset "elixir_shallowwater_well_balanced_wet_dry.jl with FluxHydrostaticReconstruction" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_well_balanced_wet_dry.jl"), + l2=[ + 0.00965787167169024, + 5.345454081916856e-14, + 0.03857583749209928, + ], + linf=[ + 0.4999999999998892, + 2.2447689894899726e-13, + 1.9999999999999714, + ], + tspan=(0.0, 0.25), + # Soften the tolerance as test results vary between different CPUs + atol=1000 * eps()) + # 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 -@trixi_testset "elixir_shallowwater_source_terms.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"), - l2=[ - 0.0022363707373868713, - 0.01576799981934617, - 4.436491725585346e-5, - ], - linf=[ - 0.00893601803417754, - 0.05939797350246456, - 9.098379777405796e-5, - ], - tspan=(0.0, 0.025)) - # 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 + @trixi_testset "elixir_shallowwater_source_terms.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_source_terms.jl"), + l2=[ + 0.0022363707373868713, + 0.01576799981934617, + 4.436491725585346e-5, + ], + linf=[ + 0.00893601803417754, + 0.05939797350246456, + 9.098379777405796e-5, + ], + tspan=(0.0, 0.025)) + # 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 -@trixi_testset "elixir_shallowwater_source_terms.jl with flux_hll" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"), - l2=[ - 0.0022758146627220154, - 0.015864082886204556, - 4.436491725585346e-5, - ], - linf=[ - 0.008457195427364006, - 0.057201667446161064, - 9.098379777405796e-5, - ], - tspan=(0.0, 0.025), - surface_flux=(FluxHLL(min_max_speed_naive), - flux_nonconservative_fjordholm_etal)) - # 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 + @trixi_testset "elixir_shallowwater_source_terms.jl with flux_hll" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_source_terms.jl"), + l2=[ + 0.0022758146627220154, + 0.015864082886204556, + 4.436491725585346e-5, + ], + linf=[ + 0.008457195427364006, + 0.057201667446161064, + 9.098379777405796e-5, + ], + tspan=(0.0, 0.025), + surface_flux=(FluxHLL(min_max_speed_naive), + flux_nonconservative_fjordholm_etal)) + # 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 -@trixi_testset "elixir_shallowwater_source_terms.jl with flux_nonconservative_ersing_etal" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"), - l2=[ - 0.005774284062933275, - 0.017408601639513584, - 4.43649172561843e-5, - ], - linf=[ - 0.01639116193303547, - 0.05102877460799604, - 9.098379777450205e-5, - ], - surface_flux=(flux_wintermeyer_etal, - flux_nonconservative_ersing_etal), - volume_flux=(flux_wintermeyer_etal, - flux_nonconservative_ersing_etal), - tspan=(0.0, 0.025)) - # 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 + @trixi_testset "elixir_shallowwater_source_terms.jl with flux_nonconservative_ersing_etal" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_source_terms.jl"), + l2=[ + 0.005774284062933275, + 0.017408601639513584, + 4.43649172561843e-5, + ], + linf=[ + 0.01639116193303547, + 0.05102877460799604, + 9.098379777450205e-5, + ], + surface_flux=(flux_wintermeyer_etal, + flux_nonconservative_ersing_etal), + volume_flux=(flux_wintermeyer_etal, + flux_nonconservative_ersing_etal), + tspan=(0.0, 0.025)) + # 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 -@trixi_testset "elixir_shallowwater_source_terms_dirichlet.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, - "elixir_shallowwater_source_terms_dirichlet.jl"), - l2=[ - 0.0022851099219788917, - 0.01560453773635554, - 4.43649172558535e-5, - ], - linf=[ - 0.008934615705174398, - 0.059403169140869405, - 9.098379777405796e-5, - ], - tspan=(0.0, 0.025)) - # 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 + @trixi_testset "elixir_shallowwater_source_terms_dirichlet.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_source_terms_dirichlet.jl"), + l2=[ + 0.0022851099219788917, + 0.01560453773635554, + 4.43649172558535e-5, + ], + linf=[ + 0.008934615705174398, + 0.059403169140869405, + 9.098379777405796e-5, + ], + tspan=(0.0, 0.025)) + # 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 -@trixi_testset "elixir_shallowwater_source_terms_dirichlet.jl with FluxHydrostaticReconstruction" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, - "elixir_shallowwater_source_terms_dirichlet.jl"), - l2=[ - 0.0022956052733432287, - 0.015540053559855601, - 4.43649172558535e-5, - ], - linf=[ - 0.008460440313118323, - 0.05720939349382359, - 9.098379777405796e-5, - ], - surface_flux=(FluxHydrostaticReconstruction(FluxHLL(min_max_speed_naive), - hydrostatic_reconstruction_audusse_etal), - flux_nonconservative_audusse_etal), - tspan=(0.0, 0.025)) - # 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 + @trixi_testset "elixir_shallowwater_source_terms_dirichlet.jl with FluxHydrostaticReconstruction" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_source_terms_dirichlet.jl"), + l2=[ + 0.0022956052733432287, + 0.015540053559855601, + 4.43649172558535e-5, + ], + linf=[ + 0.008460440313118323, + 0.05720939349382359, + 9.098379777405796e-5, + ], + surface_flux=(FluxHydrostaticReconstruction(FluxHLL(min_max_speed_naive), + hydrostatic_reconstruction_audusse_etal), + flux_nonconservative_audusse_etal), + tspan=(0.0, 0.025)) + # 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 -@trixi_testset "elixir_shallowwater_well_balanced_nonperiodic.jl with Dirichlet boundary" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, - "elixir_shallowwater_well_balanced_nonperiodic.jl"), - l2=[ - 1.725964362045055e-8, - 5.0427180314307505e-16, - 1.7259643530442137e-8, - ], - linf=[ - 3.844551077492042e-8, - 3.469453422316143e-15, - 3.844551077492042e-8, - ], - tspan=(0.0, 0.25), - surface_flux=(FluxHLL(min_max_speed_naive), - flux_nonconservative_fjordholm_etal),) - # 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 + @trixi_testset "elixir_shallowwater_well_balanced_nonperiodic.jl with Dirichlet boundary" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_well_balanced_nonperiodic.jl"), + l2=[ + 1.725964362045055e-8, + 5.0427180314307505e-16, + 1.7259643530442137e-8, + ], + linf=[ + 3.844551077492042e-8, + 3.469453422316143e-15, + 3.844551077492042e-8, + ], + tspan=(0.0, 0.25), + surface_flux=(FluxHLL(min_max_speed_naive), + flux_nonconservative_fjordholm_etal),) + # 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 -@trixi_testset "elixir_shallowwater_well_balanced_nonperiodic.jl with wall boundary" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, - "elixir_shallowwater_well_balanced_nonperiodic.jl"), - l2=[ - 1.7259643614361866e-8, - 3.5519018243195145e-16, - 1.7259643530442137e-8, - ], - linf=[ - 3.844551010878661e-8, - 9.846474508971374e-16, - 3.844551077492042e-8, - ], - tspan=(0.0, 0.25), - surface_flux=(FluxHLL(min_max_speed_naive), - flux_nonconservative_fjordholm_etal), - boundary_condition=boundary_condition_slip_wall) - # 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 + @trixi_testset "elixir_shallowwater_well_balanced_nonperiodic.jl with wall boundary" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_well_balanced_nonperiodic.jl"), + l2=[ + 1.7259643614361866e-8, + 3.5519018243195145e-16, + 1.7259643530442137e-8, + ], + linf=[ + 3.844551010878661e-8, + 9.846474508971374e-16, + 3.844551077492042e-8, + ], + tspan=(0.0, 0.25), + surface_flux=(FluxHLL(min_max_speed_naive), + flux_nonconservative_fjordholm_etal), + boundary_condition=boundary_condition_slip_wall) + # 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 -@trixi_testset "elixir_shallowwater_shock_capturing.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, - "elixir_shallowwater_shock_capturing.jl"), - l2=[0.07424140641160326, 0.2148642632748155, 0.0372579849000542], - linf=[ - 1.1209754279344226, - 1.3230788645853582, - 0.8646939843534251, - ], - tspan=(0.0, 0.05)) - # 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 + @trixi_testset "elixir_shallowwater_shock_capturing.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_shock_capturing.jl"), + l2=[ + 0.07424140641160326, + 0.2148642632748155, + 0.0372579849000542, + ], + linf=[ + 1.1209754279344226, + 1.3230788645853582, + 0.8646939843534251, + ], + tspan=(0.0, 0.05)) + # 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 -@trixi_testset "elixir_shallowwater_beach.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_beach.jl"), - l2=[ - 0.17979210479598923, - 1.2377495706611434, - 6.289818963361573e-8, - ], - linf=[ - 0.845938394800688, - 3.3740800777086575, - 4.4541473087633676e-7, - ], - tspan=(0.0, 0.05), - atol=3e-10) # see https://github.com/trixi-framework/Trixi.jl/issues/1617 - # 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 + @trixi_testset "elixir_shallowwater_beach.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_beach.jl"), + l2=[ + 0.17979210479598923, + 1.2377495706611434, + 6.289818963361573e-8, + ], + linf=[ + 0.845938394800688, + 3.3740800777086575, + 4.4541473087633676e-7, + ], + tspan=(0.0, 0.05), + atol=3e-10) # see https://github.com/trixi-framework/Trixi.jl/issues/1617 + # 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 -@trixi_testset "elixir_shallowwater_parabolic_bowl.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_parabolic_bowl.jl"), - l2=[ - 8.965981683033589e-5, - 1.8565707397810857e-5, - 4.1043039226164336e-17, - ], - linf=[ - 0.00041080213807871235, - 0.00014823261488938177, - 2.220446049250313e-16, - ], - tspan=(0.0, 0.05)) - # 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 + @trixi_testset "elixir_shallowwater_parabolic_bowl.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_parabolic_bowl.jl"), + l2=[ + 8.965981683033589e-5, + 1.8565707397810857e-5, + 4.1043039226164336e-17, + ], + linf=[ + 0.00041080213807871235, + 0.00014823261488938177, + 2.220446049250313e-16, + ], + tspan=(0.0, 0.05)) + # 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 # SWE + +@testset "Two-Layer Shallow Water" begin + @trixi_testset "elixir_shallowwater_twolayer_convergence.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_twolayer_convergence.jl"), + l2=[0.005012009872109003, 0.002091035326731071, + 0.005049271397924551, + 0.0024633066562966574, 0.0004744186597732739], + linf=[0.0213772149343594, 0.005385752427290447, + 0.02175023787351349, + 0.008212004668840978, 0.0008992474511784199], + 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 + + @trixi_testset "elixir_shallowwater_twolayer_well_balanced.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_twolayer_well_balanced.jl"), + l2=[8.949288784402005e-16, 4.0636427176237915e-17, + 0.001002881985401548, + 2.133351105037203e-16, 0.0010028819854016578], + linf=[2.6229018956769323e-15, 1.878451903240623e-16, + 0.005119880996670156, + 8.003199803957679e-16, 0.005119880996670666], + 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 + + @trixi_testset "elixir_shallowwater_twolayer_dam_break.jl with flux_lax_friedrichs" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_twolayer_dam_break.jl"), + l2=[0.1000774903431289, 0.5670692949571057, + 0.08764242501014498, + 0.45412307886094555, 0.013638618139749523], + linf=[0.586718937495144, 2.1215606128311584, + 0.5185911311186155, + 1.820382495072612, 0.5], + surface_flux=(flux_lax_friedrichs, + flux_nonconservative_ersing_etal), + 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 # 2LSWE end # TreeMesh1D # Clean up afterwards: delete TrixiShallowWater.jl output directory diff --git a/test/test_tree_2d.jl b/test/test_tree_2d.jl index e0c324c..214e84f 100644 --- a/test/test_tree_2d.jl +++ b/test/test_tree_2d.jl @@ -15,347 +15,433 @@ isdir(outdir) && rm(outdir, recursive = true) @testset "TreeMesh2D" begin #! format: noindent -@trixi_testset "elixir_shallowwater_ec.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_ec.jl"), - l2=[ - 0.991181203601035, - 0.734130029040644, - 0.7447696147162621, - 0.5875351036989047, - ], - linf=[ - 2.0117744577945413, - 2.9962317608172127, - 2.6554999727293653, - 3.0, - ], - 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 +@testset "Shallow Water Wet/Dry" begin + @trixi_testset "elixir_shallowwater_ec.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_ec.jl"), + l2=[ + 0.991181203601035, + 0.734130029040644, + 0.7447696147162621, + 0.5875351036989047, + ], + linf=[ + 2.0117744577945413, + 2.9962317608172127, + 2.6554999727293653, + 3.0, + ], + 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 -@trixi_testset "elixir_shallowwater_well_balanced.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_well_balanced.jl"), - l2=[ - 0.9130579602987144, - 1.0602847041965408e-14, - 1.082225645390032e-14, - 0.9130579602987147, - ], - linf=[ - 2.113062037615659, - 4.6613606802974e-14, - 5.4225772771633196e-14, - 2.1130620376156584, - ], - 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 + @trixi_testset "elixir_shallowwater_well_balanced.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_well_balanced.jl"), + l2=[ + 0.9130579602987144, + 1.0602847041965408e-14, + 1.082225645390032e-14, + 0.9130579602987147, + ], + linf=[ + 2.113062037615659, + 4.6613606802974e-14, + 5.4225772771633196e-14, + 2.1130620376156584, + ], + 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 -@trixi_testset "elixir_shallowwater_well_balanced_wall.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, - "elixir_shallowwater_well_balanced_wall.jl"), - l2=[ - 0.9130579602987144, - 1.0602847041965408e-14, - 1.082225645390032e-14, - 0.9130579602987147, - ], - linf=[ - 2.113062037615659, - 4.6613606802974e-14, - 5.4225772771633196e-14, - 2.1130620376156584, - ], - 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 + @trixi_testset "elixir_shallowwater_well_balanced_wall.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_well_balanced_wall.jl"), + l2=[ + 0.9130579602987144, + 1.0602847041965408e-14, + 1.082225645390032e-14, + 0.9130579602987147, + ], + linf=[ + 2.113062037615659, + 4.6613606802974e-14, + 5.4225772771633196e-14, + 2.1130620376156584, + ], + 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 -@trixi_testset "elixir_shallowwater_well_balanced.jl with FluxHydrostaticReconstruction" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_well_balanced.jl"), - l2=[ - 0.9130579602987147, - 9.68729463970494e-15, - 9.694538537436981e-15, - 0.9130579602987147, - ], - linf=[ - 2.1130620376156584, - 2.3875905654916432e-14, - 2.2492839032269154e-14, - 2.1130620376156584, - ], - surface_flux=(FluxHydrostaticReconstruction(flux_lax_friedrichs, - hydrostatic_reconstruction_audusse_etal), - flux_nonconservative_audusse_etal), - 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 + @trixi_testset "elixir_shallowwater_well_balanced.jl with FluxHydrostaticReconstruction" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_well_balanced.jl"), + l2=[ + 0.9130579602987147, + 9.68729463970494e-15, + 9.694538537436981e-15, + 0.9130579602987147, + ], + linf=[ + 2.1130620376156584, + 2.3875905654916432e-14, + 2.2492839032269154e-14, + 2.1130620376156584, + ], + surface_flux=(FluxHydrostaticReconstruction(flux_lax_friedrichs, + hydrostatic_reconstruction_audusse_etal), + flux_nonconservative_audusse_etal), + 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 -@trixi_testset "elixir_shallowwater_well_balanced.jl with flux_nonconservative_ersing_etal" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_well_balanced.jl"), - l2=[ - 0.9130579602987146, - 1.0323158914614244e-14, - 1.0276096319430528e-14, - 0.9130579602987147, - ], - linf=[ - 2.11306203761566, - 4.063916419044386e-14, - 3.694484044448245e-14, - 2.1130620376156584, - ], - surface_flux=(flux_wintermeyer_etal, - flux_nonconservative_ersing_etal), - volume_flux=(flux_wintermeyer_etal, - flux_nonconservative_ersing_etal), - 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 + @trixi_testset "elixir_shallowwater_well_balanced.jl with flux_nonconservative_ersing_etal" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_well_balanced.jl"), + l2=[ + 0.9130579602987146, + 1.0323158914614244e-14, + 1.0276096319430528e-14, + 0.9130579602987147, + ], + linf=[ + 2.11306203761566, + 4.063916419044386e-14, + 3.694484044448245e-14, + 2.1130620376156584, + ], + surface_flux=(flux_wintermeyer_etal, + flux_nonconservative_ersing_etal), + volume_flux=(flux_wintermeyer_etal, + flux_nonconservative_ersing_etal), + 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 -@trixi_testset "elixir_shallowwater_well_balanced_wet_dry.jl with FluxHydrostaticReconstruction" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, - "elixir_shallowwater_well_balanced_wet_dry.jl"), - l2=[ - 0.030186039395610056, - 2.513287752536758e-14, - 1.3631397744897607e-16, - 0.10911781485920438, - ], - linf=[ - 0.49999999999993505, - 5.5278950497971455e-14, - 7.462550826772548e-16, - 2.0, - ], - 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 + @trixi_testset "elixir_shallowwater_well_balanced_wet_dry.jl with FluxHydrostaticReconstruction" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_well_balanced_wet_dry.jl"), + l2=[ + 0.030186039395610056, + 2.513287752536758e-14, + 1.3631397744897607e-16, + 0.10911781485920438, + ], + linf=[ + 0.49999999999993505, + 5.5278950497971455e-14, + 7.462550826772548e-16, + 2.0, + ], + 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 -@trixi_testset "elixir_shallowwater_source_terms.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"), - l2=[ - 0.001868474306068482, - 0.01731687445878443, - 0.017649083171490863, - 6.274146767717023e-5, - ], - linf=[ - 0.016962486402209986, - 0.08768628853889782, - 0.09038488750767648, - 0.0001819675955490041, - ], - tspan=(0.0, 0.025)) - # 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 + @trixi_testset "elixir_shallowwater_source_terms.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_source_terms.jl"), + l2=[ + 0.001868474306068482, + 0.01731687445878443, + 0.017649083171490863, + 6.274146767717023e-5, + ], + linf=[ + 0.016962486402209986, + 0.08768628853889782, + 0.09038488750767648, + 0.0001819675955490041, + ], + tspan=(0.0, 0.025)) + # 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 -@trixi_testset "elixir_shallowwater_source_terms_dirichlet.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, - "elixir_shallowwater_source_terms_dirichlet.jl"), - l2=[ - 0.0018746929418489125, - 0.017332321628469628, - 0.01634953679145536, - 6.274146767717023e-5, - ], - linf=[ - 0.016262353691956388, - 0.08726160620859424, - 0.09043621801418844, - 0.0001819675955490041, - ], - tspan=(0.0, 0.025)) - # 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 + @trixi_testset "elixir_shallowwater_source_terms_dirichlet.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_source_terms_dirichlet.jl"), + l2=[ + 0.0018746929418489125, + 0.017332321628469628, + 0.01634953679145536, + 6.274146767717023e-5, + ], + linf=[ + 0.016262353691956388, + 0.08726160620859424, + 0.09043621801418844, + 0.0001819675955490041, + ], + tspan=(0.0, 0.025)) + # 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 -@trixi_testset "elixir_shallowwater_source_terms.jl with flux_hll" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"), - l2=[ - 0.0018957692481057034, - 0.016943229710439864, - 0.01755623297390675, - 6.274146767717414e-5, - ], - linf=[ - 0.015156105797771602, - 0.07964811135780492, - 0.0839787097210376, - 0.0001819675955490041, - ], - tspan=(0.0, 0.025), - surface_flux=(FluxHLL(min_max_speed_naive), - flux_nonconservative_fjordholm_etal)) - # 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 + @trixi_testset "elixir_shallowwater_source_terms.jl with flux_hll" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_source_terms.jl"), + l2=[ + 0.0018957692481057034, + 0.016943229710439864, + 0.01755623297390675, + 6.274146767717414e-5, + ], + linf=[ + 0.015156105797771602, + 0.07964811135780492, + 0.0839787097210376, + 0.0001819675955490041, + ], + tspan=(0.0, 0.025), + surface_flux=(FluxHLL(min_max_speed_naive), + flux_nonconservative_fjordholm_etal)) + # 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 -@trixi_testset "elixir_shallowwater_source_terms.jl with flux_nonconservative_ersing_etal" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"), - l2=[ - 0.002471853426064005, - 0.05619168608950033, - 0.11844727575152562, - 6.274146767730281e-5, - ], - linf=[ - 0.014332922987500218, - 0.2141204806174546, - 0.5392313755637872, - 0.0001819675955490041, - ], - surface_flux=(flux_wintermeyer_etal, - flux_nonconservative_ersing_etal), - volume_flux=(flux_wintermeyer_etal, - flux_nonconservative_ersing_etal), - 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 + @trixi_testset "elixir_shallowwater_source_terms.jl with flux_nonconservative_ersing_etal" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_source_terms.jl"), + l2=[ + 0.002471853426064005, + 0.05619168608950033, + 0.11844727575152562, + 6.274146767730281e-5, + ], + linf=[ + 0.014332922987500218, + 0.2141204806174546, + 0.5392313755637872, + 0.0001819675955490041, + ], + surface_flux=(flux_wintermeyer_etal, + flux_nonconservative_ersing_etal), + volume_flux=(flux_wintermeyer_etal, + flux_nonconservative_ersing_etal), + 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 -@trixi_testset "elixir_shallowwater_conical_island.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_conical_island.jl"), - l2=[ - 0.0459315416430658, - 0.1644534881916991, - 0.16445348819169914, - 0.0011537702354532694, - ], - linf=[ - 0.21100717610846464, - 0.9501592344310412, - 0.9501592344310417, - 0.021790250683516282, - ], - tspan=(0.0, 0.025)) - # 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 + @trixi_testset "elixir_shallowwater_conical_island.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_conical_island.jl"), + l2=[ + 0.0459315416430658, + 0.1644534881916991, + 0.16445348819169914, + 0.0011537702354532694, + ], + linf=[ + 0.21100717610846464, + 0.9501592344310412, + 0.9501592344310417, + 0.021790250683516282, + ], + tspan=(0.0, 0.025)) + # 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 -@trixi_testset "elixir_shallowwater_parabolic_bowl.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_parabolic_bowl.jl"), - l2=[ - 0.00025345501281482687, - 4.4525120338817177e-5, - 0.00015991819160294247, - 7.750412064917294e-15, - ], - linf=[ - 0.004664246019836723, - 0.0004972780116736669, - 0.0028735707270457628, - 6.866729407306593e-14, - ], - tspan=(0.0, 0.025), - basis=LobattoLegendreBasis(3)) - # 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 + @trixi_testset "elixir_shallowwater_parabolic_bowl.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_parabolic_bowl.jl"), + l2=[ + 0.00025345501281482687, + 4.4525120338817177e-5, + 0.00015991819160294247, + 7.750412064917294e-15, + ], + linf=[ + 0.004664246019836723, + 0.0004972780116736669, + 0.0028735707270457628, + 6.866729407306593e-14, + ], + tspan=(0.0, 0.025), + basis=LobattoLegendreBasis(3)) + # 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 -@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 + @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 # SWE + +@testset "Two-Layer Shallow Water" begin + @trixi_testset "elixir_shallowwater_twolayer_convergence.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_twolayer_convergence.jl"), + l2=[0.0004016779699408397, 0.005466339651545468, + 0.006148841330156112, + 0.0002882339012602492, 0.0030120142442780313, + 0.002680752838455618, + 8.873630921431545e-6], + linf=[0.002788654460984752, 0.01484602033450666, + 0.017572229756493973, + 0.0016010835493927011, 0.009369847995372549, + 0.008407961775489636, + 3.361991620143279e-5], + 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 + + @trixi_testset "elixir_shallowwater_twolayer_well_balanced.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_twolayer_well_balanced.jl"), + l2=[3.2935164267930016e-16, 4.6800825611195103e-17, + 4.843057532147818e-17, + 0.0030769233188015013, 1.4809161150389857e-16, + 1.509071695038043e-16, + 0.0030769233188014935], + linf=[2.248201624865942e-15, 2.346382070278936e-16, + 2.208565017494899e-16, + 0.026474051138910493, 9.237568031609006e-16, + 7.520758026187046e-16, + 0.026474051138910267], + 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 + + @trixi_testset "elixir_shallowwater_twolayer_well_balanced with flux_lax_friedrichs.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_twolayer_well_balanced.jl"), + l2=[2.0525741072929735e-16, 6.000589392730905e-17, + 6.102759428478984e-17, + 0.0030769233188014905, 1.8421386173122792e-16, + 1.8473184927121752e-16, + 0.0030769233188014935], + linf=[7.355227538141662e-16, 2.960836949170518e-16, + 4.2726562436938764e-16, + 0.02647405113891016, 1.038795478061861e-15, + 1.0401789378532516e-15, + 0.026474051138910267], + surface_flux=(flux_lax_friedrichs, + flux_nonconservative_ersing_etal), + 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 # 2LSWE end # TreeMesh2D # Clean up afterwards: delete TrixiShallowWater.jl output directory diff --git a/test/test_unit.jl b/test/test_unit.jl index 9d77370..6741761 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -57,6 +57,33 @@ end end end +@timed_testset "Two-layer shallow water conversion between conservative/entropy variables" begin + H_upper, v1_upper, v2_upper, H_lower, v1_lower, v2_lower, b = 3.5, 0.25, 0.13, 2.5, + 0.1, 0.37, 0.4 + + let equations = ShallowWaterTwoLayerEquations1D(gravity_constant = 9.8, + rho_upper = 0.9, rho_lower = 1.0) + cons_vars = prim2cons(SVector(H_upper, v1_upper, H_lower, v1_lower, b), + equations) + total_energy = energy_total(cons_vars, equations) + @test total_energy ≈ entropy(cons_vars, equations) + @test total_energy ≈ + energy_internal(cons_vars, equations) + + energy_kinetic(cons_vars, equations) + end + + let equations = ShallowWaterTwoLayerEquations2D(gravity_constant = 9.8, + rho_upper = 0.9, rho_lower = 1.0) + cons_vars = prim2cons(SVector(H_upper, v1_upper, v2_upper, H_lower, v1_lower, + v2_lower, b), equations) + total_energy = energy_total(cons_vars, equations) + @test total_energy ≈ entropy(cons_vars, equations) + @test total_energy ≈ + energy_internal(cons_vars, equations) + + energy_kinetic(cons_vars, equations) + end +end + @timed_testset "Consistency check for waterheight_pressure" begin H, v1, v2, b = 3.5, 0.25, 0.1, 0.4 @@ -146,6 +173,16 @@ end end end end + +@timed_testset "Exception check for the 2LSWE" begin + error_message = "Invalid input: Densities must be chosen such that rho_upper < rho_lower" + @test_throws error_message ShallowWaterTwoLayerEquations1D(gravity_constant = 9.81, + rho_upper = 1.0, + rho_lower = 0.9) + @test_throws error_message ShallowWaterTwoLayerEquations2D(gravity_constant = 9.81, + rho_upper = 1.0, + rho_lower = 0.9) +end end # Unit tests end # module diff --git a/test/test_unstructured_2d.jl b/test/test_unstructured_2d.jl index c2ca180..181e5aa 100644 --- a/test/test_unstructured_2d.jl +++ b/test/test_unstructured_2d.jl @@ -15,404 +15,414 @@ isdir(outdir) && rm(outdir, recursive = true) @testset "UnstructuredMesh2D" begin #! format: noindent -@trixi_testset "elixir_shallowwater_ec.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_ec.jl"), - l2=[ - 0.6106939484178353, - 0.48586236867426724, - 0.48234490854514356, - 0.29467422718511727, - ], - linf=[ - 2.775979948281604, - 3.1721242154451548, - 3.5713448319601393, - 2.052861364219655, - ], - 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 +@testset "Shallow Water Wet/Dry" begin + @trixi_testset "elixir_shallowwater_ec.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_ec.jl"), + l2=[ + 0.6106939484178353, + 0.48586236867426724, + 0.48234490854514356, + 0.29467422718511727, + ], + linf=[ + 2.775979948281604, + 3.1721242154451548, + 3.5713448319601393, + 2.052861364219655, + ], + 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 -@trixi_testset "elixir_shallowwater_well_balanced.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_well_balanced.jl"), - l2=[ - 1.2164292510839076, - 2.6118925543469468e-12, - 1.1636046671473883e-12, - 1.2164292510839079, - ], - linf=[ - 1.5138512282315846, - 4.998482888288039e-11, - 2.0246214978154587e-11, - 1.513851228231574, - ], - 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 + @trixi_testset "elixir_shallowwater_well_balanced.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_well_balanced.jl"), + l2=[ + 1.2164292510839076, + 2.6118925543469468e-12, + 1.1636046671473883e-12, + 1.2164292510839079, + ], + linf=[ + 1.5138512282315846, + 4.998482888288039e-11, + 2.0246214978154587e-11, + 1.513851228231574, + ], + 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 -@trixi_testset "elixir_shallowwater_well_balanced.jl with FluxHydrostaticReconstruction" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_well_balanced.jl"), - l2=[ - 1.2164292510839085, - 1.2643106818778908e-12, - 7.46884905098358e-13, - 1.2164292510839079, - ], - linf=[ - 1.513851228231562, - 1.6287765844373185e-11, - 6.8766999132716964e-12, - 1.513851228231574, - ], - surface_flux=(FluxHydrostaticReconstruction(flux_lax_friedrichs, - hydrostatic_reconstruction_audusse_etal), - flux_nonconservative_audusse_etal), - tspan=(0.0, 0.2)) - # 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 + @trixi_testset "elixir_shallowwater_well_balanced.jl with FluxHydrostaticReconstruction" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_well_balanced.jl"), + l2=[ + 1.2164292510839085, + 1.2643106818778908e-12, + 7.46884905098358e-13, + 1.2164292510839079, + ], + linf=[ + 1.513851228231562, + 1.6287765844373185e-11, + 6.8766999132716964e-12, + 1.513851228231574, + ], + surface_flux=(FluxHydrostaticReconstruction(flux_lax_friedrichs, + hydrostatic_reconstruction_audusse_etal), + flux_nonconservative_audusse_etal), + tspan=(0.0, 0.2)) + # 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 -@trixi_testset "elixir_shallowwater_well_balanced.jl with flux_nonconservative_ersing_etal" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_well_balanced.jl"), - l2=[ - 1.2164292510839083, - 2.590643638636187e-12, - 1.0945471514840143e-12, - 1.2164292510839079, - ], - linf=[ - 1.5138512282315792, - 5.0276441977281156e-11, - 1.9816934589292803e-11, - 1.513851228231574, - ], - surface_flux=(flux_wintermeyer_etal, - flux_nonconservative_ersing_etal), - volume_flux=(flux_wintermeyer_etal, - flux_nonconservative_ersing_etal), - 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 + @trixi_testset "elixir_shallowwater_well_balanced.jl with flux_nonconservative_ersing_etal" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_well_balanced.jl"), + l2=[ + 1.2164292510839083, + 2.590643638636187e-12, + 1.0945471514840143e-12, + 1.2164292510839079, + ], + linf=[ + 1.5138512282315792, + 5.0276441977281156e-11, + 1.9816934589292803e-11, + 1.513851228231574, + ], + surface_flux=(flux_wintermeyer_etal, + flux_nonconservative_ersing_etal), + volume_flux=(flux_wintermeyer_etal, + flux_nonconservative_ersing_etal), + 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 -@trixi_testset "elixir_shallowwater_source_terms.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"), - l2=[ - 0.0011197623982310795, - 0.04456344888447023, - 0.014317376629669337, - 5.089218476758975e-6, - ], - linf=[ - 0.007835284004819698, - 0.3486891284278597, - 0.11242778979399048, - 2.6407324614119432e-5, - ], - tspan=(0.0, 0.025)) - # 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 + @trixi_testset "elixir_shallowwater_source_terms.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_source_terms.jl"), + l2=[ + 0.0011197623982310795, + 0.04456344888447023, + 0.014317376629669337, + 5.089218476758975e-6, + ], + linf=[ + 0.007835284004819698, + 0.3486891284278597, + 0.11242778979399048, + 2.6407324614119432e-5, + ], + tspan=(0.0, 0.025)) + # 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 -@trixi_testset "elixir_shallowwater_source_terms.jl with FluxHydrostaticReconstruction" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"), - l2=[ - 0.0011197139793938152, - 0.015430259691310781, - 0.017081031802719724, - 5.089218476758271e-6, - ], - linf=[ - 0.014300809338967824, - 0.12783372461225184, - 0.17625472321992852, - 2.6407324614341476e-5, - ], - surface_flux=(FluxHydrostaticReconstruction(FluxHLL(min_max_speed_naive), - hydrostatic_reconstruction_audusse_etal), - flux_nonconservative_audusse_etal), - tspan=(0.0, 0.025)) - # 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 + @trixi_testset "elixir_shallowwater_source_terms.jl with FluxHydrostaticReconstruction" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_source_terms.jl"), + l2=[ + 0.0011197139793938152, + 0.015430259691310781, + 0.017081031802719724, + 5.089218476758271e-6, + ], + linf=[ + 0.014300809338967824, + 0.12783372461225184, + 0.17625472321992852, + 2.6407324614341476e-5, + ], + surface_flux=(FluxHydrostaticReconstruction(FluxHLL(min_max_speed_naive), + hydrostatic_reconstruction_audusse_etal), + flux_nonconservative_audusse_etal), + tspan=(0.0, 0.025)) + # 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 -@trixi_testset "elixir_shallowwater_source_terms.jl with flux_nonconservative_ersing_etal" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"), - l2=[ - 0.0011196687776346434, - 0.044562672453443995, - 0.014306265289763618, - 5.089218476759981e-6, - ], - linf=[ - 0.007825021762002393, - 0.348550815397918, - 0.1115517935018282, - 2.6407324614341476e-5, - ], - surface_flux=(flux_wintermeyer_etal, - flux_nonconservative_ersing_etal), - volume_flux=(flux_wintermeyer_etal, - flux_nonconservative_ersing_etal), - tspan=(0.0, 0.025)) - # 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 + @trixi_testset "elixir_shallowwater_source_terms.jl with flux_nonconservative_ersing_etal" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_source_terms.jl"), + l2=[ + 0.0011196687776346434, + 0.044562672453443995, + 0.014306265289763618, + 5.089218476759981e-6, + ], + linf=[ + 0.007825021762002393, + 0.348550815397918, + 0.1115517935018282, + 2.6407324614341476e-5, + ], + surface_flux=(flux_wintermeyer_etal, + flux_nonconservative_ersing_etal), + volume_flux=(flux_wintermeyer_etal, + flux_nonconservative_ersing_etal), + tspan=(0.0, 0.025)) + # 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 -@trixi_testset "elixir_shallowwater_source_terms.jl with flux_hll" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"), - l2=[ - 0.0011197139793938727, - 0.015430259691311309, - 0.017081031802719554, - 5.089218476759981e-6, - ], - linf=[ - 0.014300809338967824, - 0.12783372461224918, - 0.17625472321993918, - 2.6407324614341476e-5, - ], - surface_flux=(FluxHLL(min_max_speed_naive), - flux_nonconservative_fjordholm_etal), - tspan=(0.0, 0.025)) - # 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 + @trixi_testset "elixir_shallowwater_source_terms.jl with flux_hll" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_source_terms.jl"), + l2=[ + 0.0011197139793938727, + 0.015430259691311309, + 0.017081031802719554, + 5.089218476759981e-6, + ], + linf=[ + 0.014300809338967824, + 0.12783372461224918, + 0.17625472321993918, + 2.6407324614341476e-5, + ], + surface_flux=(FluxHLL(min_max_speed_naive), + flux_nonconservative_fjordholm_etal), + tspan=(0.0, 0.025)) + # 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 -@trixi_testset "elixir_shallowwater_dirichlet.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_dirichlet.jl"), - l2=[ - 1.1577518608940115e-5, - 4.867189932537344e-13, - 4.647273240470541e-13, - 1.1577518608933468e-5, - ], - linf=[ - 8.394063878602864e-5, - 1.1469760027632646e-10, - 1.1146619484429974e-10, - 8.394063879602065e-5, - ], - tspan=(0.0, 2.0), - surface_flux=(FluxHLL(min_max_speed_naive), - flux_nonconservative_fjordholm_etal)) - # 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 + @trixi_testset "elixir_shallowwater_dirichlet.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_dirichlet.jl"), + l2=[ + 1.1577518608940115e-5, + 4.867189932537344e-13, + 4.647273240470541e-13, + 1.1577518608933468e-5, + ], + linf=[ + 8.394063878602864e-5, + 1.1469760027632646e-10, + 1.1146619484429974e-10, + 8.394063879602065e-5, + ], + tspan=(0.0, 2.0), + surface_flux=(FluxHLL(min_max_speed_naive), + flux_nonconservative_fjordholm_etal)) + # 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 -@trixi_testset "elixir_shallowwater_wall_bc_shockcapturing.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, - "elixir_shallowwater_wall_bc_shockcapturing.jl"), - l2=[ - 0.04444388691670699, - 0.1527771788033111, - 0.1593763537203512, - 6.225080476986749e-8, - ], - linf=[ - 0.6526506870169639, - 1.980765893182952, - 2.4807635459119757, - 3.982097158683473e-7, - ], - tspan=(0.0, 0.05), - surface_flux=(FluxHydrostaticReconstruction(FluxHLL(min_max_speed_naive), - hydrostatic_reconstruction_audusse_etal), - flux_nonconservative_audusse_etal)) - # 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 + @trixi_testset "elixir_shallowwater_wall_bc_shockcapturing.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_wall_bc_shockcapturing.jl"), + l2=[ + 0.04444388691670699, + 0.1527771788033111, + 0.1593763537203512, + 6.225080476986749e-8, + ], + linf=[ + 0.6526506870169639, + 1.980765893182952, + 2.4807635459119757, + 3.982097158683473e-7, + ], + tspan=(0.0, 0.05), + surface_flux=(FluxHydrostaticReconstruction(FluxHLL(min_max_speed_naive), + hydrostatic_reconstruction_audusse_etal), + flux_nonconservative_audusse_etal)) + # 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 -@trixi_testset "elixir_shallowwater_ec_shockcapturing.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, - "elixir_shallowwater_ec_shockcapturing.jl"), - l2=[ - 0.6124656312639043, - 0.504371951785709, - 0.49180896200746366, - 0.29467422718511727, - ], - linf=[ - 2.7639232436274392, - 3.3985508653311767, - 3.3330308209196224, - 2.052861364219655, - ], - 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 + @trixi_testset "elixir_shallowwater_ec_shockcapturing.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_ec_shockcapturing.jl"), + l2=[ + 0.6124656312639043, + 0.504371951785709, + 0.49180896200746366, + 0.29467422718511727, + ], + linf=[ + 2.7639232436274392, + 3.3985508653311767, + 3.3330308209196224, + 2.052861364219655, + ], + 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 -@trixi_testset "elixir_shallowwater_three_mound_dam_break.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, - "elixir_shallowwater_three_mound_dam_break.jl"), - l2=[ - 0.0892957892027502, - 0.30648836484407915, - 2.28712547616214e-15, - 0.0008778654298684622, - ], - linf=[ - 0.850329472915091, - 2.330631694956507, - 5.783660020252348e-14, - 0.04326237921249021, - ], - basis=LobattoLegendreBasis(3), - 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 + @trixi_testset "elixir_shallowwater_three_mound_dam_break.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_three_mound_dam_break.jl"), + l2=[ + 0.0892957892027502, + 0.30648836484407915, + 2.28712547616214e-15, + 0.0008778654298684622, + ], + linf=[ + 0.850329472915091, + 2.330631694956507, + 5.783660020252348e-14, + 0.04326237921249021, + ], + basis=LobattoLegendreBasis(3), + 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 # SWE -# TODO: activate when twolayer equations are available -# @trixi_testset "elixir_shallowwater_twolayer_convergence.jl" begin -# @test_trixi_include(joinpath(EXAMPLES_DIR, -# "elixir_shallowwater_twolayer_convergence.jl"), -# l2=[0.0007935561625451243, 0.008825315509943844, -# 0.002429969315645897, -# 0.0007580145888686304, 0.004495741879625235, -# 0.0015758146898767814, -# 6.849532064729749e-6], -# linf=[0.0059205195991136605, 0.08072126590166251, -# 0.03463806075399023, -# 0.005884818649227186, 0.042658506561995546, -# 0.014125956138838602, 2.5829318284764646e-5], -# 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 +@testset "Two-Layer Shallow Water" begin + @trixi_testset "elixir_shallowwater_twolayer_convergence.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_twolayer_convergence.jl"), + l2=[0.0007935561625451243, 0.008825315509943844, + 0.002429969315645897, + 0.0007580145888686304, 0.004495741879625235, + 0.0015758146898767814, + 6.849532064729749e-6], + linf=[0.0059205195991136605, 0.08072126590166251, + 0.03463806075399023, + 0.005884818649227186, 0.042658506561995546, + 0.014125956138838602, 2.5829318284764646e-5], + 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 -# @trixi_testset "elixir_shallowwater_twolayer_well_balanced.jl" begin -# @test_trixi_include(joinpath(EXAMPLES_DIR, -# "elixir_shallowwater_twolayer_well_balanced.jl"), -# l2=[4.706532184998499e-16, 1.1215950712872183e-15, -# 6.7822712922421565e-16, -# 0.002192812926266047, 5.506855295923691e-15, -# 3.3105180099689275e-15, -# 0.0021928129262660085], -# linf=[4.468647674116255e-15, 1.3607872120431166e-14, -# 9.557155049520056e-15, -# 0.024280130945632084, 6.68910907640583e-14, -# 4.7000983997100496e-14, -# 0.024280130945632732], -# 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 + @trixi_testset "elixir_shallowwater_twolayer_well_balanced.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_twolayer_well_balanced.jl"), + l2=[4.706532184998499e-16, 1.1215950712872183e-15, + 6.7822712922421565e-16, + 0.002192812926266047, 5.506855295923691e-15, + 3.3105180099689275e-15, + 0.0021928129262660085], + linf=[4.468647674116255e-15, 1.3607872120431166e-14, + 9.557155049520056e-15, + 0.024280130945632084, 6.68910907640583e-14, + 4.7000983997100496e-14, + 0.024280130945632732], + 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 -# @trixi_testset "elixir_shallowwater_twolayer_dam_break.jl with flux_lax_friedrichs" begin -# @test_trixi_include(joinpath(EXAMPLES_DIR, -# "elixir_shallowwater_twolayer_dam_break.jl"), -# l2=[0.012447632879122346, 0.012361250464676683, -# 0.0009551519536340908, -# 0.09119400061322577, 0.015276216721920347, -# 0.0012126995108983853, 0.09991983966647647], -# linf=[0.044305765721807444, 0.03279620980615845, -# 0.010754320388190101, -# 0.111309922939555, 0.03663360204931427, -# 0.014332822306649284, -# 0.10000000000000003], -# surface_flux=(flux_lax_friedrichs, -# flux_nonconservative_ersing_etal), -# 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 + @trixi_testset "elixir_shallowwater_twolayer_dam_break.jl with flux_lax_friedrichs" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_twolayer_dam_break.jl"), + l2=[0.012447632879122346, 0.012361250464676683, + 0.0009551519536340908, + 0.09119400061322577, 0.015276216721920347, + 0.0012126995108983853, 0.09991983966647647], + linf=[0.044305765721807444, 0.03279620980615845, + 0.010754320388190101, + 0.111309922939555, 0.03663360204931427, + 0.014332822306649284, + 0.10000000000000003], + surface_flux=(flux_lax_friedrichs, + flux_nonconservative_ersing_etal), + 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 # 2LSWE +end # UnstructuredMesh2D # Clean up afterwards: delete TrixiShallowWater.jl output directory @test_nowarn rm(outdir, recursive = true)