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

@@ -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

.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

Project.toml

@@ -4,8 +4,14 @@ authors = ["Andrew R. Winters <[email protected]>", "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"

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

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

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

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