From 823f19df49048956df553dc8cebce7a5213c7135 Mon Sep 17 00:00:00 2001 From: Patrick Ersing <114223904+patrickersing@users.noreply.github.com> Date: Fri, 19 Jan 2024 20:52:38 +0100 Subject: [PATCH 1/4] Add ShallowWaterEquationsWetDry1D (#10) * add equation SWEWetDry1D * add tests and elixirs * add OrdinaryDiffEq as a test dependency * comment functions related to HR_chen_noelle * fix tests * fix tests #2 * fix tests #3 * fix tests 4 * fix tests 5 * fix tests 6 * add unit tests and namespace conflict test * update unit tests, comment wave speed estimates * update unit test * add compat bounds * add some comments * clean some comments, remove wave speed estimates * comment unused lines * include examples folder for coverage test * add upstream tests * fix ci * Adjust comment * remove comments --- .JuliaFormatter.toml | 8 + .github/workflows/ci.yml | 2 + Project.toml | 6 + .../tree_1d_dgsem/elixir_shallowwater_ec.jl | 93 +++ .../elixir_shallowwater_source_terms.jl | 60 ++ ...xir_shallowwater_source_terms_dirichlet.jl | 63 ++ .../elixir_shallowwater_well_balanced.jl | 88 +++ ..._shallowwater_well_balanced_nonperiodic.jl | 77 +++ src/TrixiShallowWater.jl | 17 +- src/equations/equations.jl | 12 + src/equations/shallow_water_wet_dry_1d.jl | 637 ++++++++++++++++++ test/Project.toml | 13 + test/runtests.jl | 16 +- test/test_tree_1d.jl | 27 + test/test_tree_1d_shallowwater_wet_dry.jl | 392 +++++++++++ test/test_trixi.jl | 238 +++++++ test/test_unit.jl | 38 ++ test/test_upstream.jl | 51 ++ 18 files changed, 1828 insertions(+), 10 deletions(-) create mode 100644 .JuliaFormatter.toml create mode 100644 examples/tree_1d_dgsem/elixir_shallowwater_ec.jl create mode 100644 examples/tree_1d_dgsem/elixir_shallowwater_source_terms.jl create mode 100644 examples/tree_1d_dgsem/elixir_shallowwater_source_terms_dirichlet.jl create mode 100644 examples/tree_1d_dgsem/elixir_shallowwater_well_balanced.jl create mode 100644 examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_nonperiodic.jl create mode 100644 src/equations/equations.jl create mode 100644 src/equations/shallow_water_wet_dry_1d.jl create mode 100644 test/test_tree_1d.jl create mode 100644 test/test_tree_1d_shallowwater_wet_dry.jl create mode 100644 test/test_trixi.jl create mode 100644 test/test_unit.jl create mode 100644 test/test_upstream.jl diff --git a/.JuliaFormatter.toml b/.JuliaFormatter.toml new file mode 100644 index 0000000..8518d20 --- /dev/null +++ b/.JuliaFormatter.toml @@ -0,0 +1,8 @@ +# Use SciML style: https://github.com/SciML/SciMLStyle +style = "sciml" + +# Python style alignment. See https://github.com/domluna/JuliaFormatter.jl/pull/732. +yas_style_nesting = true + +# Align struct fields for better readability of large struct definitions +align_struct_field = true diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index a0295b0..de6f84c 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -62,6 +62,8 @@ jobs: with: coverage: true - uses: julia-actions/julia-processcoverage@v1 + with: + directories: src,examples - uses: codecov/codecov-action@v3 with: files: lcov.info diff --git a/Project.toml b/Project.toml index a2e8f3d..1ec3ece 100644 --- a/Project.toml +++ b/Project.toml @@ -4,8 +4,14 @@ authors = ["Andrew R. Winters ", "Michael Schlottke- version = "0.1.0-pre" [deps] +MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221" +Static = "aedffcd0-7271-4cad-89d0-dc628f76c6d3" +StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Trixi = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb" [compat] +MuladdMacro = "0.2.2" +Static = "0.3, 0.4, 0.5, 0.6, 0.7, 0.8" +StaticArrays = "1" Trixi = "0.5.17, 0.6" julia = "1.8" diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_ec.jl b/examples/tree_1d_dgsem/elixir_shallowwater_ec.jl new file mode 100644 index 0000000..8046736 --- /dev/null +++ b/examples/tree_1d_dgsem/elixir_shallowwater_ec.jl @@ -0,0 +1,93 @@ + +using OrdinaryDiffEq +using Trixi +using TrixiShallowWater + +############################################################################### +# Semidiscretization of the shallow water equations with a discontinuous +# bottom topography function for a fully wet configuration + +equations = ShallowWaterEquationsWetDry1D(gravity_constant = 9.81) + +# Initial condition with a truly discontinuous water height, velocity, and bottom +# topography function as an academic testcase for entropy conservation. +# The errors from the analysis callback are not important but `∑∂S/∂U ⋅ Uₜ` should +# be around machine roundoff. +# Works as intended for TreeMesh1D with `initial_refinement_level=4`. 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_ec_discontinuous_bottom(x, t, + equations::ShallowWaterEquationsWetDry1D) + # Set the background values + H = 4.25 + v = 0.0 + b = sin(x[1]) # arbitrary continuous function + + # Setup the discontinuous water height and velocity + if x[1] >= 0.125 && x[1] <= 0.25 + H = 5.0 + v = 0.1882 + end + + # Setup a discontinuous bottom topography + if x[1] >= -0.25 && x[1] <= -0.125 + b = 2.0 + 0.5 * sin(2.0 * pi * x[1]) + end + + return prim2cons(SVector(H, v, b), equations) +end + +initial_condition = initial_condition_ec_discontinuous_bottom + +############################################################################### +# Get the DG approximation space + +volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) +solver = DGSEM(polydeg = 4, + surface_flux = (flux_fjordholm_etal, flux_nonconservative_fjordholm_etal), + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +############################################################################### +# Get the TreeMesh and setup a periodic mesh + +coordinates_min = -1.0 +coordinates_max = 1.0 +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 4, + n_cells_max = 10_000) + +# Create the semi discretization object +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +############################################################################### +# ODE solver + +tspan = (0.0, 2.0) +ode = semidiscretize(semi, tspan) + +############################################################################### +# Callbacks + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 100, + save_initial_solution = true, + save_final_solution = true) + +stepsize_callback = StepsizeCallback(cfl = 3.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/tree_1d_dgsem/elixir_shallowwater_source_terms.jl b/examples/tree_1d_dgsem/elixir_shallowwater_source_terms.jl new file mode 100644 index 0000000..aded502 --- /dev/null +++ b/examples/tree_1d_dgsem/elixir_shallowwater_source_terms.jl @@ -0,0 +1,60 @@ + +using OrdinaryDiffEq +using Trixi +using TrixiShallowWater + +############################################################################### +# Semidiscretization of the shallow water equations for a fully wet configuration + +equations = ShallowWaterEquationsWetDry1D(gravity_constant = 9.81) + +initial_condition = initial_condition_convergence_test + +############################################################################### +# Get the DG approximation space + +volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) +solver = DGSEM(polydeg = 3, + surface_flux = (flux_lax_friedrichs, flux_nonconservative_fjordholm_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 = 200, + save_initial_solution = true, + save_final_solution = true) + +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution) + +############################################################################### +# run the simulation + +# use a Runge-Kutta method with automatic (error based) time step size control +sol = solve(ode, RDPK3SpFSAL49(); abstol = 1.0e-8, reltol = 1.0e-8, + ode_default_options()..., callback = callbacks); +summary_callback() # print the timer summary diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_source_terms_dirichlet.jl b/examples/tree_1d_dgsem/elixir_shallowwater_source_terms_dirichlet.jl new file mode 100644 index 0000000..e623077 --- /dev/null +++ b/examples/tree_1d_dgsem/elixir_shallowwater_source_terms_dirichlet.jl @@ -0,0 +1,63 @@ + +using OrdinaryDiffEq +using Trixi +using TrixiShallowWater + +############################################################################### +# Semidiscretization of the shallow water equations + +equations = ShallowWaterEquationsWetDry1D(gravity_constant = 9.81) + +initial_condition = initial_condition_convergence_test + +boundary_condition = BoundaryConditionDirichlet(initial_condition) + +############################################################################### +# Get the DG approximation space + +volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) +surface_flux = (flux_lax_friedrichs, flux_nonconservative_fjordholm_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 +coordinates_max = sqrt(2.0) +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 3, + n_cells_max = 10_000, + periodicity = false) + +# create the semi discretization object +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_condition, + 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 = 200, + save_initial_solution = true, + save_final_solution = true) + +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution) + +############################################################################### +# run the simulation + +# use a Runge-Kutta method with automatic (error based) time step size control +sol = solve(ode, RDPK3SpFSAL49(); abstol = 1.0e-8, reltol = 1.0e-8, + ode_default_options()..., callback = callbacks); +summary_callback() # print the timer summary diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced.jl b/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced.jl new file mode 100644 index 0000000..e78a867 --- /dev/null +++ b/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced.jl @@ -0,0 +1,88 @@ + +using OrdinaryDiffEq +using Trixi +using TrixiShallowWater + +############################################################################### +# semidiscretization of the shallow water equations with a discontinuous +# bottom topography function for a fully wet configuration + +equations = ShallowWaterEquationsWetDry1D(gravity_constant = 9.81, H0 = 3.25) + +# Setup a truly discontinuous bottom topography function for this academic +# testcase of well-balancedness. The errors from the analysis callback are +# not important but the error for this lake-at-rest test case +# `∑|H0-(h+b)|` should be around machine roundoff. +# Works as intended for TreeMesh1D with `initial_refinement_level=3`. 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_discontinuous_well_balancedness(x, t, + equations::ShallowWaterEquationsWetDry1D) + # Set the background values + H = equations.H0 + v = 0.0 + b = 0.0 + + # Setup a discontinuous bottom topography + if x[1] >= 0.5 && x[1] <= 0.75 + b = 2.0 + 0.5 * sin(2.0 * pi * x[1]) + end + + return prim2cons(SVector(H, v, b), equations) +end + +initial_condition = initial_condition_discontinuous_well_balancedness + +############################################################################### +# Get the DG approximation space + +volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) +surface_flux = (flux_fjordholm_etal, flux_nonconservative_fjordholm_etal) +solver = DGSEM(polydeg = 4, surface_flux = surface_flux, + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +############################################################################### +# Get the TreeMesh and setup a periodic mesh + +coordinates_min = -1.0 +coordinates_max = 1.0 +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 3, + n_cells_max = 10_000) + +# Create the semi discretization object +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +############################################################################### +# ODE solver + +tspan = (0.0, 100.0) +ode = semidiscretize(semi, tspan) + +############################################################################### +# Callbacks + +summary_callback = SummaryCallback() + +analysis_interval = 1000 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + extra_analysis_integrals = (lake_at_rest_error,)) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 1000, + save_initial_solution = true, + save_final_solution = true) + +stepsize_callback = StepsizeCallback(cfl = 3.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/tree_1d_dgsem/elixir_shallowwater_well_balanced_nonperiodic.jl b/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_nonperiodic.jl new file mode 100644 index 0000000..24c60c3 --- /dev/null +++ b/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_nonperiodic.jl @@ -0,0 +1,77 @@ + +using OrdinaryDiffEq +using Trixi +using TrixiShallowWater + +############################################################################### +# Semidiscretization of the shallow water equations for a fully wet configuration + +equations = ShallowWaterEquationsWetDry1D(gravity_constant = 1.0, H0 = 3.0) + +# An initial condition with constant total water height and zero velocities to test well-balancedness. +function initial_condition_well_balancedness(x, t, equations::ShallowWaterEquationsWetDry1D) + # Set the background values + H = equations.H0 + v = 0.0 + + b = (1.5 / exp(0.5 * ((x[1] - 1.0)^2)) + 0.75 / exp(0.5 * ((x[1] + 1.0)^2))) + + return prim2cons(SVector(H, v, b), equations) +end + +initial_condition = initial_condition_well_balancedness + +boundary_condition = BoundaryConditionDirichlet(initial_condition) + +############################################################################### +# Get the DG approximation space + +volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) +solver = DGSEM(polydeg = 4, surface_flux = (flux_hll, flux_nonconservative_fjordholm_etal), + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +############################################################################### +# Get the TreeMesh and setup a non-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 = false) + +# create the semi discretization object +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_condition) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 100.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 1000 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + save_analysis = true, + extra_analysis_integrals = (lake_at_rest_error,)) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 1000, + save_initial_solution = true, + save_final_solution = true) + +stepsize_callback = StepsizeCallback(cfl = 1.0) + +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, + stepsize_callback) + +############################################################################### +# run the simulation + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); +summary_callback() # print the timer summary diff --git a/src/TrixiShallowWater.jl b/src/TrixiShallowWater.jl index 87b673f..7cfe6da 100644 --- a/src/TrixiShallowWater.jl +++ b/src/TrixiShallowWater.jl @@ -1,10 +1,17 @@ module TrixiShallowWater - +# We decided to import only Trixi.jl and qualify symbols explicitly with e.g. `Trixi.function_name`. +# For more information, see +# https://github.com/trixi-framework/TrixiShallowWater.jl/pull/10#discussion_r1433720559 using Trixi: Trixi +using MuladdMacro: @muladd +using StaticArrays: SVector +using Static: True, False + +include("equations/equations.jl") -# Write your package code here -foo() = true -bar() = false -baz() = Trixi.examples_dir() +# export types/functions that define the public API of TrixiShallowWater.jl +export ShallowWaterEquationsWetDry1D +# TODO: These function are currently exported by Trixi.jl. Needs to be uncommented when removed from Trixi.jl +#export hydrostatic_reconstruction_chen_noelle, flux_nonconservative_chen_noelle, min_max_speed_chen_noelle end diff --git a/src/equations/equations.jl b/src/equations/equations.jl new file mode 100644 index 0000000..4d31451 --- /dev/null +++ b/src/equations/equations.jl @@ -0,0 +1,12 @@ +# 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 + +#################################################################################################### +# Include files with actual implementations for different systems of equations. + +include("shallow_water_wet_dry_1d.jl") +end # @muladd diff --git a/src/equations/shallow_water_wet_dry_1d.jl b/src/equations/shallow_water_wet_dry_1d.jl new file mode 100644 index 0000000..57313ff --- /dev/null +++ b/src/equations/shallow_water_wet_dry_1d.jl @@ -0,0 +1,637 @@ +# 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""" + ShallowWaterEquationsWetDry1D(; gravity, H0 = 0, threshold_limiter = nothing threshold_wet = nothing) + +Shallow water equations (SWE) in one space dimension. The equations are given by +```math +\begin{aligned} + \frac{\partial h}{\partial t} + \frac{\partial}{\partial x}(h v) &= 0 \\ + \frac{\partial}{\partial t}(h v) + \frac{\partial}{\partial x}\left(h v^2 + \frac{g}{2}h^2\right) + + g h \frac{\partial b}{\partial x} &= 0 +\end{aligned} +``` +The unknown quantities of the SWE are the water height ``h`` and the velocity ``v``. +The gravitational constant is denoted by `g` and the (possibly) variable bottom topography function ``b(x)``. +Conservative variable water height ``h`` is measured from the bottom topography ``b``, therefore one +also defines the total water height as ``H = h + b``. + +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. + +Also, there are two thresholds which prevent numerical problems as well as instabilities. Both of them do not +have to be passed, as default values are defined within the struct. The first one, `threshold_limiter`, is +used in [`PositivityPreservingLimiterShallowWater`](@ref) on the water height, as a (small) shift on the initial +condition and cutoff before the next time step. The second one, `threshold_wet`, is applied on the water height to +define when the flow is "wet" before calculating the numerical flux. + +The bottom topography function ``b(x)`` is set inside the initial condition routine +for a particular problem setup. To test the conservative form of the SWE one can set the bottom topography +variable `b` to zero. + +In addition to the unknowns, Trixi.jl 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.jl's visualization tools will visualize the bottom topography by default. + +References for the SWE are many but a good introduction is available in Chapter 13 of the book: +- Randall J. LeVeque (2002) + Finite Volume Methods for Hyperbolic Problems + [DOI: 10.1017/CBO9780511791253](https://doi.org/10.1017/CBO9780511791253) +""" +struct ShallowWaterEquationsWetDry1D{RealT <: Real} <: + Trixi.AbstractShallowWaterEquations{1, 3} + gravity::RealT # gravitational constant + H0::RealT # constant "lake-at-rest" total water height + # `threshold_limiter` used in `PositivityPreservingLimiterShallowWater` on water height, + # as a (small) shift on the initial condition and cutoff before the next time step. + # Default is 500*eps() which in double precision is ≈1e-13. + threshold_limiter::RealT + # `threshold_wet` applied on water height to define when the flow is "wet" + # before calculating the numerical flux. + # Default is 5*eps() which in double precision is ≈1e-15. + threshold_wet::RealT +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. +# Strict default values for thresholds that performed well in many numerical experiments +function ShallowWaterEquationsWetDry1D(; gravity_constant, H0 = zero(gravity_constant), + threshold_limiter = nothing, + threshold_wet = nothing) + T = promote_type(typeof(gravity_constant), typeof(H0)) + if threshold_limiter === nothing + threshold_limiter = 500 * eps(T) + end + if threshold_wet === nothing + threshold_wet = 5 * eps(T) + end + ShallowWaterEquationsWetDry1D(gravity_constant, H0, threshold_limiter, + threshold_wet) +end + +Trixi.have_nonconservative_terms(::ShallowWaterEquationsWetDry1D) = True() +function Trixi.varnames(::typeof(Trixi.cons2cons), ::ShallowWaterEquationsWetDry1D) + ("h", "h_v", "b") +end +# Note, we use the total water height, H = h + b, as the first primitive variable for easier +# visualization and setting initial conditions +function Trixi.varnames(::typeof(Trixi.cons2prim), ::ShallowWaterEquationsWetDry1D) + ("H", "v", "b") +end + +# TODO: Remove thresholds from ShallowWaterEquations1D after they have been moved + +# This equation set extends the basic ShallowWaterEquations1D from Trixi.jl with additional functionality +# for wet/dry transitions. Since many functions correspond to the fully wet case, we make use of the +# exisiting functionality and introduce a number of wrapper functions, that dispatch to the +# ShallowWaterEquations1D. + +# Set initial conditions at physical location `x` for time `t` +""" + initial_condition_convergence_test(x, t, equations::ShallowWaterEquationsWetDry1D) + +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::ShallowWaterEquationsWetDry1D) + return Trixi.initial_condition_convergence_test(x, t, + Trixi.ShallowWaterEquations1D(equations.gravity, + equations.H0, + eps(), + eps())) +end + +""" + source_terms_convergence_test(u, x, t, equations::ShallowWaterEquationsWetDry1D) + +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). + +This manufactured solution source term is specifically designed for the bottom topography function +`b(x) = 2.0 + 0.5 * sin(sqrt(2.0)*pi*x[1])` +as defined in [`initial_condition_convergence_test`](@ref). +""" + +@inline function Trixi.source_terms_convergence_test(u, x, t, + equations::ShallowWaterEquationsWetDry1D) + return Trixi.source_terms_convergence_test(u, x, t, + Trixi.ShallowWaterEquations1D(equations.gravity, + equations.H0, + eps(), + eps())) +end + +""" + initial_condition_weak_blast_wave(x, t, equations::ShallowWaterEquationsWetDry1D) + +A weak blast wave discontinuity useful for testing, e.g., total energy conservation. +Note for the shallow water equations to the total energy acts as a mathematical entropy function. +""" +function Trixi.initial_condition_weak_blast_wave(x, t, + equations::ShallowWaterEquationsWetDry1D) + return Trixi.initial_condition_weak_blast_wave(x, t, + Trixi.ShallowWaterEquations1D(equations.gravity, + equations.H0, + eps(), + eps())) +end + +""" + boundary_condition_slip_wall(u_inner, orientation_or_normal, x, t, surface_flux_function, + equations::ShallowWaterEquationsWetDry1D) + +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::ShallowWaterEquationsWetDry1D) + return Trixi.boundary_condition_slip_wall(u_inner, orientation_or_normal, direction, + x, t, + surface_flux_function, + Trixi.ShallowWaterEquations1D(equations.gravity, + equations.H0, + eps(), + eps())) +end + +# Calculate 1D flux for a single point +# Note, the bottom topography has no flux +@inline function Trixi.flux(u, orientation::Integer, + equations::ShallowWaterEquationsWetDry1D) + return Trixi.flux(u, orientation, + Trixi.ShallowWaterEquations1D(equations.gravity, equations.H0, + eps(), eps())) +end + +""" + flux_nonconservative_wintermeyer_etal(u_ll, u_rr, orientation::Integer, + equations::ShallowWaterEquationsWetDry1D) + +Non-symmetric two-point volume flux discretizing the nonconservative (source) term +that contains the gradient of the bottom topography [`ShallowWaterEquationsWetDry1D`](@ref). + +Further details are available in the paper:#include("numerical_fluxes.jl") +- 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_nonconservative_wintermeyer_etal(u_ll, u_rr, + orientation::Integer, + equations::ShallowWaterEquationsWetDry1D) + return Trixi.flux_nonconservative_wintermeyer_etal(u_ll, u_rr, orientation, + Trixi.ShallowWaterEquations1D(equations.gravity, + equations.H0, + eps(), + eps())) +end + +""" + flux_nonconservative_fjordholm_etal(u_ll, u_rr, orientation::Integer, + equations::ShallowWaterEquationsWetDry1D) + +Non-symmetric two-point surface flux discretizing the nonconservative (source) term of +that contains the gradient of the bottom topography [`ShallowWaterEquationsWetDry1D`](@ref). + +This contains additional terms compared to [`flux_nonconservative_wintermeyer_etal`](@ref) +that account for possible discontinuities in the bottom topography function. +Thus, this flux should be used in general at interfaces. For flux differencing volume terms, +[`flux_nonconservative_wintermeyer_etal`](@ref) is analytically equivalent but slightly +cheaper. + +Further details for the original finite volume formulation are available in +- Ulrik S. Fjordholm, Siddhartha Mishr and Eitan Tadmor (2011) + Well-balanced and energy stable schemes for the shallow water equations with discontinuous topography + [DOI: 10.1016/j.jcp.2011.03.042](https://doi.org/10.1016/j.jcp.2011.03.042) +and for curvilinear 2D case in 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_nonconservative_fjordholm_etal(u_ll, u_rr, + orientation::Integer, + equations::ShallowWaterEquationsWetDry1D) + return Trixi.flux_nonconservative_fjordholm_etal(u_ll, u_rr, orientation, + Trixi.ShallowWaterEquations1D(equations.gravity, + equations.H0, + eps(), + eps())) +end + +""" + flux_nonconservative_audusse_etal(u_ll, u_rr, orientation::Integer, + equations::ShallowWaterEquationsWetDry1D) + +Non-symmetric two-point surface flux that discretizes the nonconservative (source) term. +The discretization uses the `hydrostatic_reconstruction_audusse_etal` on the conservative +variables. + +This hydrostatic reconstruction ensures that the finite volume numerical fluxes remain +well-balanced for discontinuous bottom topographies [`ShallowWaterEquationsWetDry1D`](@ref). +Should be used together with [`FluxHydrostaticReconstruction`](@ref) and +[`hydrostatic_reconstruction_audusse_etal`](@ref) in the surface flux to ensure consistency. + +Further details on the hydrostatic reconstruction and its motivation can be found in +- Emmanuel Audusse, François Bouchut, Marie-Odile Bristeau, Rupert Klein, and Benoit Perthame (2004) + A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows + [DOI: 10.1137/S1064827503431090](https://doi.org/10.1137/S1064827503431090) +""" +@inline function Trixi.flux_nonconservative_audusse_etal(u_ll, u_rr, + orientation::Integer, + equations::ShallowWaterEquationsWetDry1D) + return Trixi.flux_nonconservative_audusse_etal(u_ll, u_rr, + orientation, + Trixi.ShallowWaterEquations1D(equations.gravity, + equations.H0, + eps(), + eps())) +end + +# TODO: This function is currently exported by Trixi.jl. Needs to be uncommented when removed from Trixi.jl +# """ +# flux_nonconservative_chen_noelle(u_ll, u_rr, +# orientation::Integer, +# equations::ShallowWaterEquationsWetDry1D) + +# Non-symmetric two-point surface flux that discretizes the nonconservative (source) term. +# The discretization uses the `hydrostatic_reconstruction_chen_noelle` on the conservative +# variables. + +# Should be used together with [`FluxHydrostaticReconstruction`](@ref) and +# [`hydrostatic_reconstruction_chen_noelle`](@ref) in the surface flux to ensure consistency. + +# Further details on the hydrostatic reconstruction and its motivation can be found in +# - Guoxian Chen and Sebastian Noelle (2017) +# A new hydrostatic reconstruction scheme based on subcell reconstructions +# [DOI:10.1137/15M1053074](https://dx.doi.org/10.1137/15M1053074) +# """ +# @inline function flux_nonconservative_chen_noelle(u_ll, u_rr, +# orientation::Integer, +# equations::ShallowWaterEquationsWetDry1D) + +# # Pull the water height and bottom topography on the left +# h_ll, _, b_ll = u_ll +# h_rr, _, b_rr = u_rr + +# H_ll = h_ll + b_ll +# H_rr = h_rr + b_rr + +# b_star = min(max(b_ll, b_rr), min(H_ll, H_rr)) + +# # Create the hydrostatic reconstruction for the left solution state +# u_ll_star, _ = hydrostatic_reconstruction_chen_noelle(u_ll, u_rr, equations) + +# # Copy the reconstructed water height for easier to read code +# h_ll_star = u_ll_star[1] + +# z = zero(eltype(u_ll)) +# # Includes two parts: +# # (i) Diagonal (consistent) term from the volume flux that uses `b_ll` to avoid +# # cross-averaging across a discontinuous bottom topography +# # (ii) True surface part that uses `h_ll` and `h_ll_star` to handle discontinuous bathymetry +# return SVector(z, +# equations.gravity * h_ll * b_ll - +# equations.gravity * (h_ll_star + h_ll) * (b_ll - b_star), +# z) +# end + +""" + flux_nonconservative_ersing_etal(u_ll, u_rr, orientation::Integer, + equations::ShallowWaterEquationsWetDry1D) + +Non-symmetric path-conservative two-point volume flux discretizing the nonconservative (source) term +that contains the gradient of the bottom topography [`ShallowWaterEquationsWetDry1D`](@ref). + +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::ShallowWaterEquationsWetDry1D) + return Trixi.flux_nonconservative_ersing_etal(u_ll, u_rr, orientation, + Trixi.ShallowWaterEquations1D(equations.gravity, + equations.H0, + eps(), + eps())) +end + +""" + flux_fjordholm_etal(u_ll, u_rr, orientation, + equations::ShallowWaterEquationsWetDry1D) + +Total energy conservative (mathematical entropy for shallow water equations). When the bottom topography +is nonzero this should only be used as a surface flux otherwise the scheme will not be well-balanced. +For well-balancedness in the volume flux use [`flux_wintermeyer_etal`](@ref). + +Details are available in Eq. (4.1) in the paper: +- Ulrik S. Fjordholm, Siddhartha Mishr and Eitan Tadmor (2011) + Well-balanced and energy stable schemes for the shallow water equations with discontinuous topography + [DOI: 10.1016/j.jcp.2011.03.042](https://doi.org/10.1016/j.jcp.2011.03.042) +""" +@inline function Trixi.flux_fjordholm_etal(u_ll, u_rr, orientation::Integer, + equations::ShallowWaterEquationsWetDry1D) + return Trixi.flux_fjordholm_etal(u_ll, u_rr, orientation, + Trixi.ShallowWaterEquations1D(equations.gravity, + equations.H0, eps(), + eps())) +end + +""" + flux_wintermeyer_etal(u_ll, u_rr, orientation, + equations::ShallowWaterEquationsWetDry1D) + +Total energy conservative (mathematical entropy for shallow water equations) split form. +When the bottom topography is nonzero this scheme will be well-balanced when used as a `volume_flux`. +The `surface_flux` should still use, e.g., [`flux_fjordholm_etal`](@ref). + +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::ShallowWaterEquationsWetDry1D) + return Trixi.flux_wintermeyer_etal(u_ll, u_rr, orientation, + Trixi.ShallowWaterEquations1D(equations.gravity, + equations.H0, + eps(), eps())) +end + +""" + hydrostatic_reconstruction_audusse_etal(u_ll, u_rr, orientation::Integer, + equations::ShallowWaterEquationsWetDry1D) + +A particular type of hydrostatic reconstruction on the water height to guarantee well-balancedness +for a general bottom topography [`ShallowWaterEquationsWetDry1D`](@ref). The reconstructed solution states +`u_ll_star` and `u_rr_star` variables are then used to evaluate the surface numerical flux at the interface. +Use in combination with the generic numerical flux routine [`FluxHydrostaticReconstruction`](@ref). + +Further details on this hydrostatic reconstruction and its motivation can be found in +- Emmanuel Audusse, François Bouchut, Marie-Odile Bristeau, Rupert Klein, and Benoit Perthame (2004) + A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows + [DOI: 10.1137/S1064827503431090](https://doi.org/10.1137/S1064827503431090) +""" +@inline function Trixi.hydrostatic_reconstruction_audusse_etal(u_ll, u_rr, + equations::ShallowWaterEquationsWetDry1D) + return Trixi.hydrostatic_reconstruction_audusse_etal(u_ll, u_rr, + Trixi.ShallowWaterEquations1D(equations.gravity, + equations.H0, + eps(), + eps())) +end + +# TODO: This function is currently exported by Trixi.jl. Needs to be uncommented when removed from Trixi.jl +# """ +# hydrostatic_reconstruction_chen_noelle(u_ll, u_rr, orientation::Integer, +# equations::ShallowWaterEquationsWetDry1D) + +# A particular type of hydrostatic reconstruction of the water height to guarantee well-balancedness +# for a general bottom topography of the [`ShallowWaterEquationsWetDry1D`](@ref). The reconstructed solution states +# `u_ll_star` and `u_rr_star` variables are used to evaluate the surface numerical flux at the interface. +# The key idea is a linear reconstruction of the bottom and water height at the interfaces using subcells. +# Use in combination with the generic numerical flux routine [`FluxHydrostaticReconstruction`](@ref). + +# Further details on this hydrostatic reconstruction and its motivation can be found in +# - Guoxian Chen and Sebastian Noelle (2017) +# A new hydrostatic reconstruction scheme based on subcell reconstructions +# [DOI:10.1137/15M1053074](https://dx.doi.org/10.1137/15M1053074) +# """ +# @inline function hydrostatic_reconstruction_chen_noelle(u_ll, u_rr, +# equations::ShallowWaterEquationsWetDry1D) +# # Unpack left and right water heights and bottom topographies +# h_ll, _, b_ll = u_ll +# h_rr, _, b_rr = u_rr + +# # Get the velocities on either side +# v_ll = velocity(u_ll, equations) +# v_rr = velocity(u_rr, equations) + +# H_ll = b_ll + h_ll +# H_rr = b_rr + h_rr + +# b_star = min(max(b_ll, b_rr), min(H_ll, H_rr)) + +# # Compute the reconstructed water heights +# h_ll_star = min(H_ll - b_star, h_ll) +# h_rr_star = min(H_rr - b_star, h_rr) + +# # Set the water height to be at least the value stored in the variable threshold after +# # the hydrostatic reconstruction is applied and before the numerical flux is calculated +# # to avoid numerical problem with arbitrary small values. Interfaces with a water height +# # lower or equal to the threshold can be declared as dry. +# # The default value for `threshold_wet` is ≈ 5*eps(), or 1e-15 in double precision, is set +# # in the `ShallowWaterEquationsWetDry1D` struct. This threshold value can be changed in the constructor +# # call of this equation struct in an elixir. +# threshold = equations.threshold_wet + +# if (h_ll_star <= threshold) +# h_ll_star = threshold +# v_ll = zero(v_ll) +# end + +# if (h_rr_star <= threshold) +# h_rr_star = threshold +# v_rr = zero(v_rr) +# end + +# # Create the conservative variables using the reconstruted water heights +# u_ll_star = SVector(h_ll_star, h_ll_star * v_ll, b_ll) +# u_rr_star = SVector(h_rr_star, h_rr_star * v_rr, b_rr) + +# return u_ll_star, u_rr_star +# end + +# Specialized `DissipationLocalLaxFriedrichs` to avoid spurious dissipation in the bottom topography +@inline function (dissipation::Trixi.DissipationLocalLaxFriedrichs)(u_ll, u_rr, + orientation_or_normal_direction, + equations::ShallowWaterEquationsWetDry1D) + return (dissipation::Trixi.DissipationLocalLaxFriedrichs)(u_ll, u_rr, + orientation_or_normal_direction, + Trixi.ShallowWaterEquations1D(equations.gravity, + equations.H0, + eps(), + eps())) +end + +# Specialized `FluxHLL` to avoid spurious dissipation in the bottom topography +@inline function (numflux::Trixi.FluxHLL)(u_ll, u_rr, orientation_or_normal_direction, + equations::ShallowWaterEquationsWetDry1D) + return (numflux::Trixi.FluxHLL)(u_ll, u_rr, orientation_or_normal_direction, + Trixi.ShallowWaterEquations1D(equations.gravity, + equations.H0, eps(), + eps())) +end + +# TODO: This function is currently exported by Trixi.jl. Needs to be uncommented when removed from Trixi.jl +# """ +# min_max_speed_chen_noelle(u_ll, u_rr, orientation::Integer, +# equations::ShallowWaterEquationsWetDry1D) + +# The approximated speeds for the HLL type numerical flux used by Chen and Noelle for their +# hydrostatic reconstruction. As they state in the paper, these speeds are chosen for the numerical +# flux to ensure positivity and to satisfy an entropy inequality. + +# Further details on this hydrostatic reconstruction and its motivation can be found in +# - Guoxian Chen and Sebastian Noelle (2017) +# A new hydrostatic reconstruction scheme based on subcell reconstructions +# [DOI:10.1137/15M1053074](https://dx.doi.org/10.1137/15M1053074) +# """ +# @inline function min_max_speed_chen_noelle(u_ll, u_rr, orientation::Integer, +# equations::ShallowWaterEquationsWetDry1D) +# # Get the velocity quantities +# v_ll = velocity(u_ll, equations) +# v_rr = velocity(u_rr, equations) + +# # Calculate the wave celerity on the left and right +# 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) + +# λ_min = min(v_ll - a_ll, v_rr - a_rr, zero(eltype(u_ll))) +# λ_max = max(v_ll + a_ll, v_rr + a_rr, zero(eltype(u_ll))) + +# return λ_min, λ_max +# end + +@inline function Trixi.max_abs_speeds(u, equations::ShallowWaterEquationsWetDry1D) + return Trixi.max_abs_speeds(u, + Trixi.ShallowWaterEquations1D(equations.gravity, + equations.H0, eps(), + eps())) +end + +# # Helper function to extract the velocity vector from the conservative variables +# @inline function Trixi.velocity(u, equations::ShallowWaterEquationsWetDry1D) +# return Trixi.velocity(u, +# Trixi.ShallowWaterEquations1D(equations.gravity, equations.H0, +# eps(), eps())) +# end + +# Convert conservative variables to primitive +@inline function Trixi.cons2prim(u, equations::ShallowWaterEquationsWetDry1D) + return Trixi.cons2prim(u, + Trixi.ShallowWaterEquations1D(equations.gravity, + equations.H0, eps(), eps())) +end + +# Convert conservative variables to entropy +# Note, only the first two are the entropy variables, the third entry still +# just carries the bottom topography values for convenience +@inline function Trixi.cons2entropy(u, equations::ShallowWaterEquationsWetDry1D) + return Trixi.cons2entropy(u, + Trixi.ShallowWaterEquations1D(equations.gravity, + equations.H0, eps(), eps())) +end + +# Convert entropy variables to conservative +@inline function Trixi.entropy2cons(w, equations::ShallowWaterEquationsWetDry1D) + return Trixi.entropy2cons(w, + Trixi.ShallowWaterEquations1D(equations.gravity, + equations.H0, eps(), eps())) +end + +# Convert primitive to conservative variables +@inline function Trixi.prim2cons(prim, equations::ShallowWaterEquationsWetDry1D) + return Trixi.prim2cons(prim, + Trixi.ShallowWaterEquations1D(equations.gravity, + equations.H0, eps(), eps())) +end + +# @inline function Trixi.waterheight(u, equations::ShallowWaterEquationsWetDry1D) +# return Trixi.waterheight(u, +# Trixi.ShallowWaterEquations1D(equations.gravity, +# equations.H0, eps(), eps())) +# end + +# @inline function Trixi.pressure(u, equations::ShallowWaterEquationsWetDry1D) +# return Trixi.pressure(u, +# Trixi.ShallowWaterEquations1D(equations.gravity, equations.H0, +# eps(), eps())) +# end + +# @inline function Trixi.waterheight_pressure(u, equations::ShallowWaterEquationsWetDry1D) +# return Trixi.waterheight_pressure(u, +# Trixi.ShallowWaterEquations1D(equations.gravity, +# equations.H0, eps(), +# eps())) +# end + +# Entropy function for the shallow water equations is the total energy +@inline function Trixi.entropy(cons, equations::ShallowWaterEquationsWetDry1D) + Trixi.energy_total(cons, equations) +end + +# Calculate total energy for a conservative state `cons` +@inline function Trixi.energy_total(cons, equations::ShallowWaterEquationsWetDry1D) + return Trixi.energy_total(cons, + Trixi.ShallowWaterEquations1D(equations.gravity, + equations.H0, eps(), eps())) +end + +# Calculate kinetic energy for a conservative state `cons` +@inline function Trixi.energy_kinetic(u, equations::ShallowWaterEquationsWetDry1D) + return Trixi.energy_kinetic(u, + Trixi.ShallowWaterEquations1D(equations.gravity, + equations.H0, eps(), + eps())) +end + +# Calculate potential energy for a conservative state `cons` +@inline function Trixi.energy_internal(cons, equations::ShallowWaterEquationsWetDry1D) + return Trixi.energy_total(cons, equations) - Trixi.energy_kinetic(cons, equations) +end + +# Calculate the error for the "lake-at-rest" test case where H = h+b should +# be a constant value over time. Note, assumes there is a single reference +# water height `H0` with which to compare. +@inline function Trixi.lake_at_rest_error(u, equations::ShallowWaterEquationsWetDry1D) + h, _, b = u + + # For well-balancedness testing with possible wet/dry regions the reference + # water height `H0` accounts for the possibility that the bottom topography + # can emerge out of the water as well as for the threshold offset to avoid + # division by a "hard" zero water heights as well. + H0_wet_dry = max(equations.H0, b + equations.threshold_limiter) + + return abs(H0_wet_dry - (h + b)) +end +end # @muladd diff --git a/test/Project.toml b/test/Project.toml index fbc64c5..cfac638 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -1,6 +1,19 @@ [deps] Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" Trixi = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb" +OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" [compat] +Test = "1" Trixi = "0.5, 0.6" +OrdinaryDiffEq = "6.49.1" + +[preferences.OrdinaryDiffEq] +PrecompileAutoSpecialize = false +PrecompileAutoSwitch = false +PrecompileDefaultSpecialize = false +PrecompileFunctionWrapperSpecialize = false +PrecompileLowStorage = false +PrecompileNoSpecialize = false +PrecompileNonStiff = false +PrecompileStiff = false diff --git a/test/runtests.jl b/test/runtests.jl index 909d0e2..8a0e590 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,3 +1,4 @@ +using Trixi using TrixiShallowWater using Test @@ -11,14 +12,19 @@ const TRIXI_NTHREADS = clamp(Sys.CPU_THREADS, 2, 3) @time @testset "TrixiShallowWater.jl tests" begin @time if TRIXI_TEST == "all" - @test TrixiShallowWater.foo() == true - @test TrixiShallowWater.bar() == false - @test TrixiShallowWater.baz() isa String + include("test_tree_1d_shallowwater_wet_dry.jl") + include("test_unit.jl") end @time if TRIXI_TEST == "all" || TRIXI_TEST == "upstream" - @testset "baz()" begin - @test TrixiShallowWater.baz() isa String + @testset "Namespace conflicts" begin + # Test for namespace conflicts between TrixiShallowWater.jl and Trixi.jl + for name in names(Trixi) + @test !(name in names(TrixiShallowWater)) + end end + + # Run upstream tests for each mesh and dimension to test compatibility with Trixi.jl + include("test_upstream.jl") end end diff --git a/test/test_tree_1d.jl b/test/test_tree_1d.jl new file mode 100644 index 0000000..7d96552 --- /dev/null +++ b/test/test_tree_1d.jl @@ -0,0 +1,27 @@ +module TestExamplesTree1D + +using Test +using TrixiShallowWater + +include("test_trixi.jl") + +EXAMPLES_DIR = pkgdir(TrixiShallowWater, "examples", "tree_1d_dgsem") + +# Start with a clean environment: remove Trixi.jl output directory if it exists +outdir = "out" +isdir(outdir) && rm(outdir, recursive = true) + +@testset "TreeMesh1D" begin +#! format: noindent + +# Run basic tests +@testset "Examples 1D" begin + # Shallow water + include("test_tree_1d_shallowwater_wet_dry.jl") +end + +# Clean up afterwards: delete Trixi.jl output directory +@test_nowarn rm(outdir, recursive = true) +end # TreeMesh1D + +end # module diff --git a/test/test_tree_1d_shallowwater_wet_dry.jl b/test/test_tree_1d_shallowwater_wet_dry.jl new file mode 100644 index 0000000..e2cec1a --- /dev/null +++ b/test/test_tree_1d_shallowwater_wet_dry.jl @@ -0,0 +1,392 @@ +module TestExamples1DShallowWaterWetDry + +using Test +using Trixi +using TrixiShallowWater + +include("test_trixi.jl") + +EXAMPLES_DIR = pkgdir(TrixiShallowWater, "examples", "tree_1d_dgsem") + +@testset "Shallow Water Wet Dry" 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 + 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 + 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 + 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 + 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 + end +end + +# TODO: Activate test when Wet_Dry functionality is moved +# @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 + +@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 + +@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=(flux_hll, 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 + +@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 + +@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 + +@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(flux_hll, + 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 + +@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)) + # 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_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), + 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 + +# TODO: Activate test when Wet_Dry functionality is moved +# @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 + +# TODO: Activate test when Wet_Dry functionality is moved +# @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 + +# TODO: Activate test when Wet_Dry functionality is moved +# @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 # module diff --git a/test/test_trixi.jl b/test/test_trixi.jl new file mode 100644 index 0000000..cebe216 --- /dev/null +++ b/test/test_trixi.jl @@ -0,0 +1,238 @@ +using Test: @test +import Trixi + +# Use a macro to avoid world age issues when defining new initial conditions etc. +# inside an elixir. +""" + @test_trixi_include(elixir; l2=nothing, linf=nothing, + atol=500*eps(), rtol=sqrt(eps()), + parameters...) + +Test Trixi by calling `trixi_include(elixir; parameters...)`. +By default, only the absence of error output is checked. +If `l2` or `linf` are specified, in addition the resulting L2/Linf errors +are compared approximately against these reference values, using `atol, rtol` +as absolute/relative tolerance. +""" +macro test_trixi_include(elixir, args...) + local l2 = get_kwarg(args, :l2, nothing) + local linf = get_kwarg(args, :linf, nothing) + local atol = get_kwarg(args, :atol, 500 * eps()) + local rtol = get_kwarg(args, :rtol, sqrt(eps())) + local skip_coverage = get_kwarg(args, :skip_coverage, false) + local coverage_override = expr_to_named_tuple(get_kwarg(args, :coverage_override, :())) + if !(:maxiters in keys(coverage_override)) + # maxiters in coverage_override defaults to 1 + coverage_override = (; coverage_override..., maxiters = 1) + end + + local cmd = string(Base.julia_cmd()) + local coverage = occursin("--code-coverage", cmd) && + !occursin("--code-coverage=none", cmd) + + local kwargs = Pair{Symbol, Any}[] + for arg in args + if (arg.head == :(=) && + !(arg.args[1] in (:l2, :linf, :atol, :rtol, :coverage_override, :skip_coverage)) + && !(coverage && arg.args[1] in keys(coverage_override))) + push!(kwargs, Pair(arg.args...)) + end + end + + if coverage + for key in keys(coverage_override) + push!(kwargs, Pair(key, coverage_override[key])) + end + end + + if coverage && skip_coverage + return quote + if Trixi.mpi_isroot() + println("═"^100) + println("Skipping coverage test of ", $elixir) + println("═"^100) + println("\n\n") + end + end + end + + quote + Trixi.mpi_isroot() && println("═"^100) + Trixi.mpi_isroot() && println($elixir) + + # if `maxiters` is set in tests, it is usually set to a small number to + # run only a few steps - ignore possible warnings coming from that + if any(==(:maxiters) ∘ first, $kwargs) + additional_ignore_content = [ + r"┌ Warning: Interrupted\. Larger maxiters is needed\..*\n└ @ SciMLBase .+\n", + r"┌ Warning: Interrupted\. Larger maxiters is needed\..*\n└ @ Trixi .+\n"] + else + additional_ignore_content = [] + end + + # evaluate examples in the scope of the module they're called from + @test_nowarn_mod trixi_include(@__MODULE__, $elixir; $kwargs...) additional_ignore_content + + # if present, compare l2 and linf errors against reference values + if !$coverage && (!isnothing($l2) || !isnothing($linf)) + l2_measured, linf_measured = analysis_callback(sol) + + if Trixi.mpi_isroot() && !isnothing($l2) + @test length($l2) == length(l2_measured) + for (l2_expected, l2_actual) in zip($l2, l2_measured) + @test isapprox(l2_expected, l2_actual, atol = $atol, rtol = $rtol) + end + end + + if Trixi.mpi_isroot() && !isnothing($linf) + @test length($linf) == length(linf_measured) + for (linf_expected, linf_actual) in zip($linf, linf_measured) + @test isapprox(linf_expected, linf_actual, atol = $atol, rtol = $rtol) + end + end + end + + Trixi.mpi_isroot() && println("═"^100) + Trixi.mpi_isroot() && println("\n\n") + end +end + +# Get the first value assigned to `keyword` in `args` and return `default_value` +# if there are no assignments to `keyword` in `args`. +function get_kwarg(args, keyword, default_value) + val = default_value + for arg in args + if arg.head == :(=) && arg.args[1] == keyword + val = arg.args[2] + break + end + end + return val +end + +function expr_to_named_tuple(expr) + result = (;) + + for arg in expr.args + if arg.head != :(=) + error("Invalid expression") + end + result = (; result..., arg.args[1] => arg.args[2]) + end + return result +end + +# Modified version of `@test_nowarn` that prints the content of `stderr` when +# it is not empty and ignores module replacements. +macro test_nowarn_mod(expr, additional_ignore_content = String[]) + quote + let fname = tempname() + try + ret = open(fname, "w") do f + redirect_stderr(f) do + $(esc(expr)) + end + end + stderr_content = read(fname, String) + if !isempty(stderr_content) + println("Content of `stderr`:\n", stderr_content) + end + + # Patterns matching the following ones will be ignored. Additional patterns + # passed as arguments can also be regular expressions, so we just use the + # type `Any` for `ignore_content`. + ignore_content = Any[ + # We need to ignore steady state information reported by our callbacks + r"┌ Info: Steady state tolerance reached\n│ steady_state_callback .+\n└ t = .+\n", + # We also ignore our own compilation messages + "[ Info: You just called `trixi_include`. Julia may now compile the code, please be patient.\n", + # TODO: Upstream (PlotUtils). This should be removed again once the + # deprecated stuff is fixed upstream. + "WARNING: importing deprecated binding Colors.RGB1 into Plots.\n", + "WARNING: importing deprecated binding Colors.RGB4 into Plots.\n", + r"┌ Warning: Keyword argument letter not supported with Plots.+\n└ @ Plots.+\n", + r"┌ Warning: `parse\(::Type, ::Coloarant\)` is deprecated.+\n│.+\n│.+\n└ @ Plots.+\n", + # TODO: Silence warning introduced by Flux v0.13.13. Should be properly fixed. + r"┌ Warning: Layer with Float32 parameters got Float64 input.+\n│.+\n│.+\n│.+\n└ @ Flux.+\n"] + append!(ignore_content, $additional_ignore_content) + for pattern in ignore_content + stderr_content = replace(stderr_content, pattern => "") + end + + # We also ignore simple module redefinitions for convenience. Thus, we + # check whether every line of `stderr_content` is of the form of a + # module replacement warning. + @test occursin(r"^(WARNING: replacing module .+\.\n)*$", stderr_content) + ret + finally + rm(fname, force = true) + end + end + end +end + +""" + @timed_testset "name of the testset" #= code to test #= + +Similar to `@testset`, but prints the name of the testset and its runtime +after execution. +""" +macro timed_testset(name, expr) + @assert name isa String + quote + local time_start = time_ns() + @testset $name $expr + local time_stop = time_ns() + if Trixi.mpi_isroot() + flush(stdout) + @info("Testset "*$name*" finished in " + *string(1.0e-9 * (time_stop - time_start))*" seconds.\n") + flush(stdout) + end + end +end + +""" + @trixi_testset "name of the testset" #= code to test #= + +Similar to `@testset`, but wraps the code inside a temporary module to avoid +namespace pollution. It also `include`s this file again to provide the +definition of `@test_trixi_include`. Moreover, it records the execution time +of the testset similarly to [`timed_testset`](@ref). +""" +macro trixi_testset(name, expr) + @assert name isa String + # TODO: `@eval` is evil + # We would like to use + # mod = gensym(name) + # ... + # module $mod + # to create new module names for every test set. However, this is not + # compatible with the dirty hack using `@eval` to get the mapping when + # loading structured, curvilinear meshes. Thus, we need to use a plain + # module name here. + quote + local time_start = time_ns() + @eval module TrixiTestModule + using Test + using Trixi + include(@__FILE__) + # We define `EXAMPLES_DIR` in (nearly) all test modules and use it to + # get the path to the elixirs to be tested. However, that's not required + # and we want to fail gracefully if it's not defined. + try + import ..EXAMPLES_DIR + catch + nothing + end + @testset $name $expr + end + local time_stop = time_ns() + if Trixi.mpi_isroot() + flush(stdout) + @info("Testset "*$name*" finished in " + *string(1.0e-9 * (time_stop - time_start))*" seconds.\n") + end + nothing + end +end diff --git a/test/test_unit.jl b/test/test_unit.jl new file mode 100644 index 0000000..9517443 --- /dev/null +++ b/test/test_unit.jl @@ -0,0 +1,38 @@ +module TestUnit + +using Test +using Trixi +using TrixiShallowWater + +include("test_trixi.jl") + +# Start with a clean environment: remove Trixi.jl output directory if it exists +outdir = "out" +isdir(outdir) && rm(outdir, recursive = true) + +# Run various unit (= non-elixir-triggered) tests +@testset "Unit tests" begin +#! format: noindent + +@timed_testset "Shallow water conversion between conservative/entropy variables" begin + H, v, b = 3.5, 0.25, 0.4 + + let equations = ShallowWaterEquationsWetDry1D(gravity_constant = 9.8) + cons_vars = Trixi.prim2cons(SVector(H, v, b), equations) + entropy_vars = Trixi.cons2entropy(cons_vars, equations) + @test cons_vars ≈ Trixi.entropy2cons(entropy_vars, equations) + + total_energy = Trixi.energy_total(cons_vars, equations) + @test total_energy ≈ Trixi.entropy(cons_vars, equations) + @test total_energy ≈ + Trixi.energy_internal(cons_vars, equations) + + energy_kinetic(cons_vars, equations) + # test tuple args + cons_vars = Trixi.prim2cons((H, v, b), equations) + entropy_vars = Trixi.cons2entropy(cons_vars, equations) + @test cons_vars ≈ Trixi.entropy2cons(entropy_vars, equations) + end +end +end + +end # module diff --git a/test/test_upstream.jl b/test/test_upstream.jl new file mode 100644 index 0000000..56c9850 --- /dev/null +++ b/test/test_upstream.jl @@ -0,0 +1,51 @@ +module TestExamplesUpstream + +using Test +using TrixiShallowWater + +include("test_trixi.jl") + +EXAMPLES_DIR = pkgdir(TrixiShallowWater, "examples") + +# Start with a clean environment: remove output directory if it exists +outdir = "out" +isdir(outdir) && rm(outdir, recursive = true) + +# Run upstream tests for each mesh and dimension to test compatibility with Trixi.jl +@testset "Upstream tests" begin +#! format: noindent + +# Run tests for TreeMesh +@testset "TreeMesh" begin + # Shallow water wet/dry 1D + @trixi_testset "1D-Test: elixir_shallowwater_well_balanced_nonperiodic.jl with wall boundary" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "tree_1d_dgsem", + "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), + 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 + +# Clean up afterwards: delete output directory +@test_nowarn rm(outdir, recursive = true) +end # Upstream tests + +end # module From 69818fff6803fd5346881cc7368632df680d9f9c Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Sun, 21 Jan 2024 20:03:15 +0100 Subject: [PATCH 2/4] coverallsapp at v2 (#9) Co-authored-by: Patrick Ersing <114223904+patrickersing@users.noreply.github.com> Co-authored-by: Andrew Winters --- .github/workflows/ci.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index de6f84c..6aa1eaa 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -67,7 +67,7 @@ jobs: - uses: codecov/codecov-action@v3 with: files: lcov.info - - uses: coverallsapp/github-action@master + - uses: coverallsapp/github-action@v2 with: github-token: ${{ secrets.GITHUB_TOKEN }} path-to-lcov: ./lcov.info From cd51eff4e9fa96e6fdd9c368421ce5b302b54c4c Mon Sep 17 00:00:00 2001 From: Patrick Ersing <114223904+patrickersing@users.noreply.github.com> Date: Wed, 24 Jan 2024 20:05:09 +0100 Subject: [PATCH 3/4] Fix downstream tests (#19) * remove thresholds from call signature * remove waterheight from call signature --- src/equations/shallow_water_wet_dry_1d.jl | 101 ++++++---------------- 1 file changed, 25 insertions(+), 76 deletions(-) diff --git a/src/equations/shallow_water_wet_dry_1d.jl b/src/equations/shallow_water_wet_dry_1d.jl index 57313ff..a19b6c9 100644 --- a/src/equations/shallow_water_wet_dry_1d.jl +++ b/src/equations/shallow_water_wet_dry_1d.jl @@ -112,10 +112,7 @@ A smooth initial condition used for convergence tests in combination with function Trixi.initial_condition_convergence_test(x, t, equations::ShallowWaterEquationsWetDry1D) return Trixi.initial_condition_convergence_test(x, t, - Trixi.ShallowWaterEquations1D(equations.gravity, - equations.H0, - eps(), - eps())) + Trixi.ShallowWaterEquations1D(gravity_constant = equations.gravity)) end """ @@ -133,10 +130,7 @@ as defined in [`initial_condition_convergence_test`](@ref). @inline function Trixi.source_terms_convergence_test(u, x, t, equations::ShallowWaterEquationsWetDry1D) return Trixi.source_terms_convergence_test(u, x, t, - Trixi.ShallowWaterEquations1D(equations.gravity, - equations.H0, - eps(), - eps())) + Trixi.ShallowWaterEquations1D(gravity_constant = equations.gravity)) end """ @@ -148,10 +142,7 @@ Note for the shallow water equations to the total energy acts as a mathematical function Trixi.initial_condition_weak_blast_wave(x, t, equations::ShallowWaterEquationsWetDry1D) return Trixi.initial_condition_weak_blast_wave(x, t, - Trixi.ShallowWaterEquations1D(equations.gravity, - equations.H0, - eps(), - eps())) + Trixi.ShallowWaterEquations1D(gravity_constant = equations.gravity)) end """ @@ -176,10 +167,7 @@ For details see Section 9.2.5 of the book: return Trixi.boundary_condition_slip_wall(u_inner, orientation_or_normal, direction, x, t, surface_flux_function, - Trixi.ShallowWaterEquations1D(equations.gravity, - equations.H0, - eps(), - eps())) + Trixi.ShallowWaterEquations1D(gravity_constant = equations.gravity)) end # Calculate 1D flux for a single point @@ -187,8 +175,7 @@ end @inline function Trixi.flux(u, orientation::Integer, equations::ShallowWaterEquationsWetDry1D) return Trixi.flux(u, orientation, - Trixi.ShallowWaterEquations1D(equations.gravity, equations.H0, - eps(), eps())) + Trixi.ShallowWaterEquations1D(gravity_constant = equations.gravity)) end """ @@ -208,10 +195,7 @@ Further details are available in the paper:#include("numerical_fluxes.jl") orientation::Integer, equations::ShallowWaterEquationsWetDry1D) return Trixi.flux_nonconservative_wintermeyer_etal(u_ll, u_rr, orientation, - Trixi.ShallowWaterEquations1D(equations.gravity, - equations.H0, - eps(), - eps())) + Trixi.ShallowWaterEquations1D(gravity_constant = equations.gravity)) end """ @@ -241,10 +225,7 @@ and for curvilinear 2D case in the paper: orientation::Integer, equations::ShallowWaterEquationsWetDry1D) return Trixi.flux_nonconservative_fjordholm_etal(u_ll, u_rr, orientation, - Trixi.ShallowWaterEquations1D(equations.gravity, - equations.H0, - eps(), - eps())) + Trixi.ShallowWaterEquations1D(gravity_constant = equations.gravity)) end """ @@ -270,10 +251,7 @@ Further details on the hydrostatic reconstruction and its motivation can be foun equations::ShallowWaterEquationsWetDry1D) return Trixi.flux_nonconservative_audusse_etal(u_ll, u_rr, orientation, - Trixi.ShallowWaterEquations1D(equations.gravity, - equations.H0, - eps(), - eps())) + Trixi.ShallowWaterEquations1D(gravity_constant = equations.gravity)) end # TODO: This function is currently exported by Trixi.jl. Needs to be uncommented when removed from Trixi.jl @@ -345,10 +323,7 @@ For further details see: orientation::Integer, equations::ShallowWaterEquationsWetDry1D) return Trixi.flux_nonconservative_ersing_etal(u_ll, u_rr, orientation, - Trixi.ShallowWaterEquations1D(equations.gravity, - equations.H0, - eps(), - eps())) + Trixi.ShallowWaterEquations1D(gravity_constant = equations.gravity)) end """ @@ -367,9 +342,7 @@ Details are available in Eq. (4.1) in the paper: @inline function Trixi.flux_fjordholm_etal(u_ll, u_rr, orientation::Integer, equations::ShallowWaterEquationsWetDry1D) return Trixi.flux_fjordholm_etal(u_ll, u_rr, orientation, - Trixi.ShallowWaterEquations1D(equations.gravity, - equations.H0, eps(), - eps())) + Trixi.ShallowWaterEquations1D(gravity_constant = equations.gravity)) end """ @@ -389,9 +362,7 @@ Further details are available in Theorem 1 of the paper: @inline function Trixi.flux_wintermeyer_etal(u_ll, u_rr, orientation::Integer, equations::ShallowWaterEquationsWetDry1D) return Trixi.flux_wintermeyer_etal(u_ll, u_rr, orientation, - Trixi.ShallowWaterEquations1D(equations.gravity, - equations.H0, - eps(), eps())) + Trixi.ShallowWaterEquations1D(gravity_constant = equations.gravity)) end """ @@ -411,10 +382,7 @@ Further details on this hydrostatic reconstruction and its motivation can be fou @inline function Trixi.hydrostatic_reconstruction_audusse_etal(u_ll, u_rr, equations::ShallowWaterEquationsWetDry1D) return Trixi.hydrostatic_reconstruction_audusse_etal(u_ll, u_rr, - Trixi.ShallowWaterEquations1D(equations.gravity, - equations.H0, - eps(), - eps())) + Trixi.ShallowWaterEquations1D(gravity_constant = equations.gravity)) end # TODO: This function is currently exported by Trixi.jl. Needs to be uncommented when removed from Trixi.jl @@ -484,19 +452,14 @@ end equations::ShallowWaterEquationsWetDry1D) return (dissipation::Trixi.DissipationLocalLaxFriedrichs)(u_ll, u_rr, orientation_or_normal_direction, - Trixi.ShallowWaterEquations1D(equations.gravity, - equations.H0, - eps(), - eps())) + Trixi.ShallowWaterEquations1D(gravity_constant = equations.gravity)) end # Specialized `FluxHLL` to avoid spurious dissipation in the bottom topography @inline function (numflux::Trixi.FluxHLL)(u_ll, u_rr, orientation_or_normal_direction, equations::ShallowWaterEquationsWetDry1D) return (numflux::Trixi.FluxHLL)(u_ll, u_rr, orientation_or_normal_direction, - Trixi.ShallowWaterEquations1D(equations.gravity, - equations.H0, eps(), - eps())) + Trixi.ShallowWaterEquations1D(gravity_constant = equations.gravity)) end # TODO: This function is currently exported by Trixi.jl. Needs to be uncommented when removed from Trixi.jl @@ -534,23 +497,19 @@ end @inline function Trixi.max_abs_speeds(u, equations::ShallowWaterEquationsWetDry1D) return Trixi.max_abs_speeds(u, - Trixi.ShallowWaterEquations1D(equations.gravity, - equations.H0, eps(), - eps())) + Trixi.ShallowWaterEquations1D(gravity_constant = equations.gravity)) end # # Helper function to extract the velocity vector from the conservative variables # @inline function Trixi.velocity(u, equations::ShallowWaterEquationsWetDry1D) # return Trixi.velocity(u, -# Trixi.ShallowWaterEquations1D(equations.gravity, equations.H0, -# eps(), eps())) +# Trixi.ShallowWaterEquations1D(gravity_constant = equations.gravity)) # end # Convert conservative variables to primitive @inline function Trixi.cons2prim(u, equations::ShallowWaterEquationsWetDry1D) return Trixi.cons2prim(u, - Trixi.ShallowWaterEquations1D(equations.gravity, - equations.H0, eps(), eps())) + Trixi.ShallowWaterEquations1D(gravity_constant = equations.gravity)) end # Convert conservative variables to entropy @@ -558,41 +517,34 @@ end # just carries the bottom topography values for convenience @inline function Trixi.cons2entropy(u, equations::ShallowWaterEquationsWetDry1D) return Trixi.cons2entropy(u, - Trixi.ShallowWaterEquations1D(equations.gravity, - equations.H0, eps(), eps())) + Trixi.ShallowWaterEquations1D(gravity_constant = equations.gravity)) end # Convert entropy variables to conservative @inline function Trixi.entropy2cons(w, equations::ShallowWaterEquationsWetDry1D) return Trixi.entropy2cons(w, - Trixi.ShallowWaterEquations1D(equations.gravity, - equations.H0, eps(), eps())) + Trixi.ShallowWaterEquations1D(gravity_constant = equations.gravity)) end # Convert primitive to conservative variables @inline function Trixi.prim2cons(prim, equations::ShallowWaterEquationsWetDry1D) return Trixi.prim2cons(prim, - Trixi.ShallowWaterEquations1D(equations.gravity, - equations.H0, eps(), eps())) + Trixi.ShallowWaterEquations1D(gravity_constant = equations.gravity)) end # @inline function Trixi.waterheight(u, equations::ShallowWaterEquationsWetDry1D) # return Trixi.waterheight(u, -# Trixi.ShallowWaterEquations1D(equations.gravity, -# equations.H0, eps(), eps())) +# Trixi.ShallowWaterEquations1D(gravity_constant = equations.gravity)) # end # @inline function Trixi.pressure(u, equations::ShallowWaterEquationsWetDry1D) # return Trixi.pressure(u, -# Trixi.ShallowWaterEquations1D(equations.gravity, equations.H0, -# eps(), eps())) +# Trixi.ShallowWaterEquations1D(gravity_constant = equations.gravity)) # end # @inline function Trixi.waterheight_pressure(u, equations::ShallowWaterEquationsWetDry1D) # return Trixi.waterheight_pressure(u, -# Trixi.ShallowWaterEquations1D(equations.gravity, -# equations.H0, eps(), -# eps())) +# Trixi.ShallowWaterEquations1D(gravity_constant = equations.gravity)) # end # Entropy function for the shallow water equations is the total energy @@ -603,16 +555,13 @@ end # Calculate total energy for a conservative state `cons` @inline function Trixi.energy_total(cons, equations::ShallowWaterEquationsWetDry1D) return Trixi.energy_total(cons, - Trixi.ShallowWaterEquations1D(equations.gravity, - equations.H0, eps(), eps())) + Trixi.ShallowWaterEquations1D(gravity_constant = equations.gravity)) end # Calculate kinetic energy for a conservative state `cons` @inline function Trixi.energy_kinetic(u, equations::ShallowWaterEquationsWetDry1D) return Trixi.energy_kinetic(u, - Trixi.ShallowWaterEquations1D(equations.gravity, - equations.H0, eps(), - eps())) + Trixi.ShallowWaterEquations1D(gravity_constant = equations.gravity)) end # Calculate potential energy for a conservative state `cons` From 6bb65dd8ea9f89aebd350810f23954c30ff144a2 Mon Sep 17 00:00:00 2001 From: Patrick Ersing <114223904+patrickersing@users.noreply.github.com> Date: Wed, 31 Jan 2024 14:25:58 +0100 Subject: [PATCH 4/4] Add spell check and format check workflows (#17) * add permissions for cache deletion * add format check workflow * add spell check * fix typo * apply formatter * Update .github/workflows/FormatCheck.yml * apply formatter --------- Co-authored-by: Andrew Winters Co-authored-by: Michael Schlottke-Lakemper --- .github/workflows/FormatCheck.yml | 41 +++++++++++++++++++++++ .github/workflows/SpellCheck.yml | 13 +++++++ .github/workflows/ci.yml | 5 +++ docs/make.jl | 34 +++++++++---------- src/equations/shallow_water_wet_dry_1d.jl | 2 +- test/runtests.jl | 4 +-- test/test_upstream.jl | 2 +- 7 files changed, 79 insertions(+), 22 deletions(-) create mode 100644 .github/workflows/FormatCheck.yml create mode 100644 .github/workflows/SpellCheck.yml diff --git a/.github/workflows/FormatCheck.yml b/.github/workflows/FormatCheck.yml new file mode 100644 index 0000000..1310d5d --- /dev/null +++ b/.github/workflows/FormatCheck.yml @@ -0,0 +1,41 @@ +name: format-check + +on: + push: + branches: + - 'main' + tags: '*' + pull_request: + +jobs: + check-format: + runs-on: ${{ matrix.os }} + strategy: + matrix: + julia-version: [1] + julia-arch: [x86] + os: [ubuntu-latest] + steps: + - uses: julia-actions/setup-julia@latest + with: + version: ${{ matrix.julia-version }} + + - uses: actions/checkout@v4 + - name: Install JuliaFormatter and format + # This will use the latest version by default but you can set the version like so: + # + # julia -e 'using Pkg; Pkg.add(PackageSpec(name = "JuliaFormatter", version = "0.13.0"))' + run: | + julia -e 'using Pkg; Pkg.add(PackageSpec(name = "JuliaFormatter"))' + julia -e 'using JuliaFormatter; format(".")' + - name: Format check + run: | + julia -e ' + out = Cmd(`git diff --name-only`) |> read |> String + if out == "" + exit(0) + else + @error "Some files have not been formatted !!!" + write(stdout, out) + exit(1) + end' \ No newline at end of file diff --git a/.github/workflows/SpellCheck.yml b/.github/workflows/SpellCheck.yml new file mode 100644 index 0000000..a780e97 --- /dev/null +++ b/.github/workflows/SpellCheck.yml @@ -0,0 +1,13 @@ +name: Spell Check + +on: [pull_request, workflow_dispatch] + +jobs: + typos-check: + name: Spell Check with Typos + runs-on: ubuntu-latest + steps: + - name: Checkout Actions Repository + uses: actions/checkout@v4 + - name: Check spelling + uses: crate-ci/typos@v1.16.26 diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 6aa1eaa..3ed7d99 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -28,6 +28,11 @@ concurrency: group: ${{ github.workflow }}-${{ github.ref }} cancel-in-progress: ${{ startsWith(github.ref, 'refs/pull/') }} +# needed to allow julia-actions/cache to delete old caches that it has created +permissions: + actions: write + contents: read + jobs: test: if: "!contains(github.event.head_commit.message, 'skip ci')" diff --git a/docs/make.jl b/docs/make.jl index 03bb874..bc59a9f 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,25 +1,23 @@ using TrixiShallowWater using Documenter -DocMeta.setdocmeta!(TrixiShallowWater, :DocTestSetup, :(using TrixiShallowWater); recursive=true) +DocMeta.setdocmeta!(TrixiShallowWater, :DocTestSetup, :(using TrixiShallowWater); + recursive = true) makedocs(; - modules=[TrixiShallowWater], - authors="Andrew R. Winters , Michael Schlottke-Lakemper ", - repo="https://github.com/trixi-framework/TrixiShallowWater.jl/blob/{commit}{path}#{line}", - sitename="TrixiShallowWater.jl", - format=Documenter.HTML(; - prettyurls=get(ENV, "CI", "false") == "true", - canonical="https://trixi-framework.github.io/TrixiShallowWater.jl", - edit_link="main", - assets=String[], - ), - pages=[ - "Home" => "index.md", - ], -) + modules = [TrixiShallowWater], + authors = "Andrew R. Winters , Michael Schlottke-Lakemper ", + repo = "https://github.com/trixi-framework/TrixiShallowWater.jl/blob/{commit}{path}#{line}", + sitename = "TrixiShallowWater.jl", + format = Documenter.HTML(; + prettyurls = get(ENV, "CI", "false") == "true", + canonical = "https://trixi-framework.github.io/TrixiShallowWater.jl", + edit_link = "main", + assets = String[],), + pages = [ + "Home" => "index.md", + ],) deploydocs(; - repo="github.com/trixi-framework/TrixiShallowWater.jl", - devbranch="main", -) + repo = "github.com/trixi-framework/TrixiShallowWater.jl", + devbranch = "main",) diff --git a/src/equations/shallow_water_wet_dry_1d.jl b/src/equations/shallow_water_wet_dry_1d.jl index a19b6c9..03cb118 100644 --- a/src/equations/shallow_water_wet_dry_1d.jl +++ b/src/equations/shallow_water_wet_dry_1d.jl @@ -97,7 +97,7 @@ end # This equation set extends the basic ShallowWaterEquations1D from Trixi.jl with additional functionality # for wet/dry transitions. Since many functions correspond to the fully wet case, we make use of the -# exisiting functionality and introduce a number of wrapper functions, that dispatch to the +# existing functionality and introduce a number of wrapper functions, that dispatch to the # ShallowWaterEquations1D. # Set initial conditions at physical location `x` for time `t` diff --git a/test/runtests.jl b/test/runtests.jl index 8a0e590..92a68d6 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -23,8 +23,8 @@ const TRIXI_NTHREADS = clamp(Sys.CPU_THREADS, 2, 3) @test !(name in names(TrixiShallowWater)) end end - + # Run upstream tests for each mesh and dimension to test compatibility with Trixi.jl - include("test_upstream.jl") + include("test_upstream.jl") end end diff --git a/test/test_upstream.jl b/test/test_upstream.jl index 56c9850..3beb128 100644 --- a/test/test_upstream.jl +++ b/test/test_upstream.jl @@ -20,7 +20,7 @@ isdir(outdir) && rm(outdir, recursive = true) # Shallow water wet/dry 1D @trixi_testset "1D-Test: elixir_shallowwater_well_balanced_nonperiodic.jl with wall boundary" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "tree_1d_dgsem", - "elixir_shallowwater_well_balanced_nonperiodic.jl"), + "elixir_shallowwater_well_balanced_nonperiodic.jl"), l2=[ 1.7259643614361866e-8, 3.5519018243195145e-16,