Shallow water exner equations 1D (#51)

* add SWE-Exner1D

* fix docstrings

* add zero velocity workaround to roe flux

* fix typo

* add roots

* fix convergence tests

* fix convergence tests 2

* fix convergence tests 3

* update test_tree_1d reference values

* add tests to increase coverage

* Apply suggestions from code review

Co-authored-by: Andrew Winters <[email protected]>

* Modify docstrings, fix function q_s

---------

Co-authored-by: Andrew Winters <[email protected]>

Project.toml

@@ -6,16 +6,18 @@ version = "0.1.0-pre"
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221"
+Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
+Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
 Static = "aedffcd0-7271-4cad-89d0-dc628f76c6d3"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Trixi = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb"
-Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 
 [compat]
 LinearAlgebra = "1"
 MuladdMacro = "0.2.2"
+Printf = "1"
+Roots = "2.1.6"
 Static = "0.3, 0.4, 0.5, 0.6, 0.7, 0.8"
 StaticArrays = "1"
 Trixi = "0.7, 0.8"
 julia = "1.8"
-Printf = "1"

examples/tree_1d_dgsem/elixir_shallowwater_exner_channel.jl

@@ -0,0 +1,116 @@
+
+using OrdinaryDiffEq
+using Trixi
+using TrixiShallowWater
+using Roots
+
+###############################################################################
+# Semidiscretization of the shallow water exner equations for a channel flow problem 
+# with sediment transport
+
+equations = ShallowWaterExnerEquations1D(gravity_constant = 9.81, rho_f = 1.0, rho_s = 1.0,
+                                         porosity = 0.4,
+                                         sediment_model = GrassModel(A_g = 0.01))
+
+# Initial condition for for a channel flow problem over a sediment hump.
+# An asymptotic solution based on the method of characteristics was derived under a rigid-lid 
+# approximation in chapter 3.5.1 of the thesis:
+# - Justin Hudson (2001)
+#   "Numerical Techniques for Morphodynamic Modelling"
+#   [PhD Thesis, University of Reading](https://citeseerx.ist.psu.edu/document?repid=rep1&type=pdf&doi=f78dbae9cfbb12ae823975c6ce9d2585b40417ba)
+function initial_condition_channel(x, t,
+                                   equations::ShallowWaterExnerEquations1D)
+    (; sediment_model, porosity_inv) = equations
+
+    H_ref = 10.0    # Reference water level
+    hv = 10.0
+
+    # Use the method of characteristics to find the asymptotic solution for the bed height, see 
+    # Eq. 3.16 in the reference.
+    # First use an iterative method to predict x0 
+    f(x0) = x[1] - x0 -
+            sediment_model.A_g * porosity_inv * 3 * hv^3 * t *
+            (H_ref - sinpi((x0 - 300) / 200)^2)^(-4)
+    fx = Roots.ZeroProblem(f, 400.0)
+    x0 = Roots.solve(fx, atol = 0.0, rtol = 0.0)
+
+    # If the result is outside 300 < x < 500 the result is invalid and instead compute x0 from
+    if x0 > 500 || x0 < 300
+        x0 = x[1] -
+             sediment_model.A_g * porosity_inv * 3 * hv^3 * t * H_ref^(-4)
+    end
+
+    # Compute the sediment and water height
+    300 < x0 < 500 ? h_b = 0.1 + sinpi((x0 - 300) / 200)^2 : h_b = 0.1
+    h = H_ref - h_b
+
+    return SVector(h, hv, h_b)
+end
+
+initial_condition = initial_condition_channel
+
+###############################################################################
+# Get the DG approximation space
+
+volume_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal)
+surface_flux = (FluxPlusDissipation(flux_ersing_etal, dissipation_roe),
+                flux_nonconservative_ersing_etal)
+
+basis = LobattoLegendreBasis(6)
+
+indicator_sc = IndicatorHennemannGassner(equations, basis,
+                                         alpha_max = 0.5,
+                                         alpha_min = 0.001,
+                                         alpha_smooth = true,
+                                         variable = water_sediment_height)
+volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
+                                                 volume_flux_dg = volume_flux,
+                                                 volume_flux_fv = surface_flux)
+
+solver = DGSEM(basis, surface_flux, volume_integral)
+
+###############################################################################
+# Get the TreeMesh and setup a periodic mesh
+
+coordinates_min = 0.0
+coordinates_max = 1000.0
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 4,
+                n_cells_max = 10_000,
+                periodicity = true)
+
+# Create the semi discretization object
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+###############################################################################
+# ODE solver
+
+tspan = (0.0, 30_000.0)
+ode = semidiscretize(semi, tspan)
+
+###############################################################################
+# Callbacks
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 10000
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+stepsize_callback = StepsizeCallback(cfl = 0.8)
+
+save_solution = SaveSolutionCallback(dt = 10_000.0,
+                                     save_initial_solution = true,
+                                     save_final_solution = true)
+
+callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback,
+                        stepsize_callback, save_solution)
+
+###############################################################################
+# 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_exner_source_terms_grass.jl

@@ -0,0 +1,69 @@
+using OrdinaryDiffEq
+using Trixi
+using TrixiShallowWater
+
+###############################################################################
+# Semidiscretization of the SWE-Exner equations with source terms for convergence testing
+
+# Equations with Grass model
+equations = ShallowWaterExnerEquations1D(gravity_constant = 10.0, rho_f = 0.5,
+                                         rho_s = 1.0, porosity = 0.5,
+                                         friction = ManningFriction(n = 0.0),
+                                         sediment_model = GrassModel(A_g = 0.01))
+
+initial_condition = initial_condition_convergence_test
+
+###############################################################################
+# Get the DG approximation space
+
+volume_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal)
+surface_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal)
+
+solver = DGSEM(polydeg = 4,
+               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 = 2,
+                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 = 200
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+stepsize_callback = StepsizeCallback(cfl = 1.0)
+
+save_solution = SaveSolutionCallback(interval = 200,
+                                     save_initial_solution = true,
+                                     save_final_solution = true)
+
+callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback,
+                        stepsize_callback, save_solution)
+
+###############################################################################
+# 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_exner_source_terms_mpm.jl

@@ -0,0 +1,70 @@
+using OrdinaryDiffEq
+using Trixi
+using TrixiShallowWater
+
+###############################################################################
+# Semidiscretization of the SWE-Exner equations with source terms for convergence testing
+
+# Equations with Meyer-Peter-Mueller model
+equations = ShallowWaterExnerEquations1D(gravity_constant = 10.0, rho_f = 0.5,
+                                         rho_s = 1.0, porosity = 0.5,
+                                         friction = ManningFriction(n = 0.01),
+                                         sediment_model = MeyerPeterMueller(theta_c = 0.0,
+                                                                            d_s = 1e-3))
+
+initial_condition = initial_condition_convergence_test
+
+###############################################################################
+# Get the DG approximation space
+
+volume_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal)
+surface_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal)
+
+solver = DGSEM(polydeg = 4,
+               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 = 2,
+                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 = 200
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+stepsize_callback = StepsizeCallback(cfl = 1.0)
+
+save_solution = SaveSolutionCallback(interval = 200,
+                                     save_initial_solution = true,
+                                     save_final_solution = true)
+
+callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback,
+                        stepsize_callback, save_solution)
+
+###############################################################################
+# 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_exner_well_balanced.jl

@@ -0,0 +1,104 @@
+
+using OrdinaryDiffEq
+using Trixi
+using TrixiShallowWater
+
+###############################################################################
+# Semidiscretization of the SWE-Exner equations with a discontinuous sediment bed 
+# to test well-balancedness
+
+# Equations with Meyer-Peter-Mueller sedimentation model
+equations = ShallowWaterExnerEquations1D(gravity_constant = 10.0, H0 = 1.0,
+                                         rho_f = 0.5, rho_s = 1.0, porosity = 0.5,
+                                         friction = ManningFriction(n = 0.01),
+                                         sediment_model = MeyerPeterMueller(theta_c = 0.047,
+                                                                            d_s = 1e-3))
+
+function initial_condition_steady_state(x, t,
+                                        equations::ShallowWaterExnerEquations1D)
+    hv = 0.0
+
+    # discontinuous sediment bed
+    if -2.0 < x[1] < 0.0
+        h_s = 0.3
+    else
+        h_s = 0.1
+    end
+
+    h = 1.0 - h_s
+
+    return SVector(h, hv, h_s)
+end
+
+initial_condition = initial_condition_steady_state
+
+boundary_condition = boundary_condition_slip_wall
+
+###############################################################################
+# Get the DG approximation space
+
+volume_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal)
+surface_flux = (FluxPlusDissipation(flux_ersing_etal, DissipationLocalLaxFriedrichs()),
+                flux_nonconservative_ersing_etal)
+
+basis = LobattoLegendreBasis(3)
+
+indicator_sc = IndicatorHennemannGassner(equations, basis,
+                                         alpha_max = 0.5,
+                                         alpha_min = 0.001,
+                                         alpha_smooth = true,
+                                         variable = water_sediment_height)
+volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
+                                                 volume_flux_dg = volume_flux,
+                                                 volume_flux_fv = surface_flux)
+
+solver = DGSEM(basis, surface_flux, volume_integral)
+
+###############################################################################
+# Get the TreeMesh
+
+coordinates_min = -2.0
+coordinates_max = 2.0
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 4,
+                n_cells_max = 10_000,
+                periodicity = false)
+
+# Create the semi discretization object
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    source_terms = source_term_bottom_friction,
+                                    boundary_conditions = boundary_condition)
+
+###############################################################################
+# ODE solver
+
+tspan = (0.0, 200.0)
+ode = semidiscretize(semi, tspan)
+
+###############################################################################
+# Callbacks
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 10000
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
+                                     extra_analysis_integrals = (lake_at_rest_error,))
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(dt = 100,
+                                     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