[WIP] Add 1D multilayer shallow water equations (#30)

* add first implementation of  the ml-swe in 1D

* fix typo and formatting

* add well-balanced test and lake_at_rest_error

* add additional analysis functions and unit tests

* add wall_bc, specialized dissipation and tests

* switch to three-layer convergence test

* update reference values for tests

* Adjust comments

* apply formatter

* apply formatter

* add changes from code review

* fix comment

* add comment for energy_total

* add unit test for initial_condition_convergence

examples/tree_1d_dgsem/elixir_shallowwater_multilayer_convergence.jl

@@ -0,0 +1,61 @@
+
+using OrdinaryDiffEq
+using Trixi
+using TrixiShallowWater
+
+###############################################################################
+# Semidiscretization of the multilayer shallow water equations with three layers
+
+equations = ShallowWaterMultiLayerEquations1D(gravity_constant = 10.0,
+                                              rhos = (0.9, 1.0, 1.1))
+
+initial_condition = initial_condition_convergence_test
+
+###############################################################################
+# Get the DG approximation space
+
+volume_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal)
+solver = DGSEM(polydeg = 3,
+               surface_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal),
+               volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+
+###############################################################################
+# Get the TreeMesh and setup a periodic mesh
+
+coordinates_min = 0.0
+coordinates_max = sqrt(2.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,
+                                    source_terms = source_terms_convergence_test)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 1.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 500
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(interval = 500,
+                                     save_initial_solution = true,
+                                     save_final_solution = true)
+
+callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution)
+
+###############################################################################
+# run the simulation
+
+# use a Runge-Kutta method with automatic (error based) time step size control
+sol = solve(ode, RDPK3SpFSAL49(), abstol = 1.0e-8, reltol = 1.0e-8,
+            save_everystep = false, callback = callbacks);
+summary_callback() # print the timer summary

examples/tree_1d_dgsem/elixir_shallowwater_multilayer_dam_break_ec.jl

@@ -0,0 +1,92 @@
+
+using OrdinaryDiffEq
+using Trixi
+using TrixiShallowWater
+
+###############################################################################
+# Semidiscretization of the multilayer shallow water equations for a dam break
+# test with a discontinuous bottom topography function to test entropy conservation
+
+equations = ShallowWaterMultiLayerEquations1D(gravity_constant = 9.81,
+                                              rhos = (0.85, 0.9, 1.0))
+
+# Initial condition of a dam break with a discontinuous water heights and bottom topography.
+# Works as intended for TreeMesh1D with `initial_refinement_level=5`. If the mesh
+# refinement level is changed the initial condition below may need changed as well to
+# ensure that the discontinuities lie on an element interface.
+function initial_condition_dam_break(x, t, equations::ShallowWaterMultiLayerEquations1D)
+    v = [0.0, 0.0, 0.0]
+
+    # Set the discontinuity
+    if x[1] <= 10.0
+        H = [6.0, 4.0, 2.0]
+        b = 0.0
+    else
+        H = [5.5, 3.5, 1.5]
+        b = 0.5
+    end
+
+    return prim2cons(SVector(H..., v..., b), equations)
+end
+
+initial_condition = initial_condition_dam_break
+
+###############################################################################
+# Get the DG approximation space
+
+volume_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal)
+solver = DGSEM(polydeg = 3,
+               surface_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal),
+               volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+
+###############################################################################
+# Get the TreeMesh and setup a non-periodic mesh
+
+coordinates_min = 0.0
+coordinates_max = 20.0
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 5,
+                n_cells_max = 10000,
+                periodicity = false)
+
+boundary_condition = boundary_condition_slip_wall
+
+# create the semidiscretization object
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    boundary_conditions = boundary_condition)
+
+###############################################################################
+# ODE solvers
+
+tspan = (0.0, 0.4)
+ode = semidiscretize(semi, tspan)
+
+###############################################################################
+# Callbacks
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 500
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
+                                     save_analysis = false,
+                                     extra_analysis_integrals = (energy_total,
+                                                                 energy_kinetic,
+                                                                 energy_internal))
+
+stepsize_callback = StepsizeCallback(cfl = 1.0)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(interval = 500,
+                                     save_initial_solution = true,
+                                     save_final_solution = true)
+
+callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution)
+
+###############################################################################
+# run the simulation
+
+# use a Runge-Kutta method with automatic (error based) time step size control
+sol = solve(ode, RDPK3SpFSAL49(), abstol = 1.0e-8, reltol = 1.0e-8,
+            save_everystep = false, callback = callbacks);
+summary_callback() # print the timer summary

examples/tree_1d_dgsem/elixir_shallowwater_multilayer_dam_break_es.jl

@@ -0,0 +1,94 @@
+
+using OrdinaryDiffEq
+using Trixi
+using TrixiShallowWater
+
+###############################################################################
+# Semidiscretization of the multilayer shallow water equations for a dam break
+# test with a discontinuous bottom topography function for an entropy stable flux
+
+equations = ShallowWaterMultiLayerEquations1D(gravity_constant = 9.81,
+                                              rhos = (0.85, 0.9, 1.0))
+
+# Initial condition of a dam break with a discontinuous water heights and bottom topography.
+# Works as intended for TreeMesh1D with `initial_refinement_level=5`. If the mesh
+# refinement level is changed the initial condition below may need changed as well to
+# ensure that the discontinuities lie on an element interface.
+function initial_condition_dam_break(x, t, equations::ShallowWaterMultiLayerEquations1D)
+    v = [0.0, 0.0, 0.0]
+
+    # Set the discontinuity
+    if x[1] <= 10.0
+        H = [6.0, 4.0, 2.0]
+        b = 0.0
+    else
+        H = [5.5, 3.5, 1.5]
+        b = 0.5
+    end
+
+    return prim2cons(SVector(H..., v..., b), equations)
+end
+
+initial_condition = initial_condition_dam_break
+
+###############################################################################
+# Get the DG approximation space
+
+volume_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal)
+solver = DGSEM(polydeg = 3,
+               surface_flux = (FluxPlusDissipation(flux_ersing_etal,
+                                                   DissipationLocalLaxFriedrichs()),
+                               flux_nonconservative_ersing_etal),
+               volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+
+###############################################################################
+# Get the TreeMesh and setup a non-periodic mesh
+
+coordinates_min = 0.0
+coordinates_max = 20.0
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 5,
+                n_cells_max = 10000,
+                periodicity = false)
+
+boundary_condition = boundary_condition_slip_wall
+
+# create the semidiscretization object
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    boundary_conditions = boundary_condition)
+
+###############################################################################
+# ODE solvers
+
+tspan = (0.0, 0.4)
+ode = semidiscretize(semi, tspan)
+
+###############################################################################
+# Callbacks
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 500
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
+                                     save_analysis = false,
+                                     extra_analysis_integrals = (energy_total,
+                                                                 energy_kinetic,
+                                                                 energy_internal))
+
+stepsize_callback = StepsizeCallback(cfl = 1.0)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(interval = 500,
+                                     save_initial_solution = true,
+                                     save_final_solution = true)
+
+callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution)
+
+###############################################################################
+# run the simulation
+
+# use a Runge-Kutta method with automatic (error based) time step size control
+sol = solve(ode, RDPK3SpFSAL49(), abstol = 1.0e-8, reltol = 1.0e-8,
+            save_everystep = false, callback = callbacks);
+summary_callback() # print the timer summary

examples/tree_1d_dgsem/elixir_shallowwater_multilayer_well_balanced.jl

@@ -0,0 +1,86 @@
+
+using OrdinaryDiffEq
+using Trixi
+using TrixiShallowWater
+
+###############################################################################
+# Semidiscretization of the multilayer shallow water equations to test well-balancedness
+
+equations = ShallowWaterMultiLayerEquations1D(gravity_constant = 1.0, H0 = 0.7,
+                                              rhos = (0.8, 0.9, 1.0))
+
+"""
+    initial_condition_fjordholm_well_balanced(x, t, equations::ShallowWaterMultiLayerEquations1D)
+
+Initial condition to test well balanced with a bottom topography adapted from Fjordholm
+- Ulrik Skre Fjordholm (2012)
+  Energy conservative and stable schemes for the two-layer shallow water equations.
+  [DOI: 10.1142/9789814417099_0039](https://doi.org/10.1142/9789814417099_0039)
+"""
+function initial_condition_fjordholm_well_balanced(x, t,
+                                                   equations::ShallowWaterMultiLayerEquations1D)
+    inicenter = 0.5
+    x_norm = x[1] - inicenter
+    r = abs(x_norm)
+
+    H = [0.7, 0.6, 0.5]
+    v = [0.0, 0.0, 0.0]
+    b = r <= 0.1 ? 0.2 * (cos(10 * pi * (x[1] - 0.5)) + 1) : 0.0
+
+    return prim2cons(SVector(H..., v..., b), equations)
+end
+
+initial_condition = initial_condition_fjordholm_well_balanced
+
+###############################################################################
+# Get the DG approximation space
+
+volume_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal)
+solver = DGSEM(polydeg = 3,
+               surface_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal),
+               volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+
+###############################################################################
+# Get the TreeMesh and setup a periodic mesh
+
+coordinates_min = 0.0
+coordinates_max = 1.0
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 4,
+                n_cells_max = 10_000,
+                periodicity = true)
+
+# create the semi discretization object
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 10.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 1000
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
+                                     save_analysis = false,
+                                     extra_analysis_integrals = (lake_at_rest_error,))
+
+stepsize_callback = StepsizeCallback(cfl = 1.0)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(interval = 1000,
+                                     save_initial_solution = true,
+                                     save_final_solution = true)
+
+callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution,
+                        stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
+            dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep = false, callback = callbacks);
+summary_callback() # print the timer summary

src/TrixiShallowWater.jl

@@ -8,7 +8,7 @@ using Trixi
 # Import additional symbols that are not exported by Trixi.jl
 using Trixi: get_node_vars, set_node_vars!, waterheight
 using MuladdMacro: @muladd
-using StaticArrays: SVector, @SMatrix
+using StaticArrays: SVector, @SMatrix, MVector
 using Static: True, False
 using LinearAlgebra: norm
 
@@ -19,15 +19,17 @@ include("solvers/indicators.jl")
 
 # Export types/functions that define the public API of TrixiShallowWater.jl
 export ShallowWaterEquationsWetDry1D, ShallowWaterEquationsWetDry2D,
-       ShallowWaterTwoLayerEquations1D, ShallowWaterTwoLayerEquations2D
+       ShallowWaterTwoLayerEquations1D, ShallowWaterTwoLayerEquations2D,
+       ShallowWaterMultiLayerEquations1D
 
 export hydrostatic_reconstruction_chen_noelle, flux_nonconservative_chen_noelle,
-       min_max_speed_chen_noelle,
-       flux_hll_chen_noelle
+       min_max_speed_chen_noelle, flux_hll_chen_noelle,
+       flux_ersing_etal, flux_es_ersing_etal
 
-export flux_es_ersing_etal
+export nlayers, eachlayer
 
 export PositivityPreservingLimiterShallowWater
+
 export IndicatorHennemannGassnerShallowWater
 
 end