Add wet/dry functionality for the Multilayer SWE in 1D (#38)

* add wet/dry functionality for the 1d mlswe

* add new elixirs to test wet/dry functionality

* add tests

* apply formatter

* update test values for new max_abs_speeds_naive

* additional wet/dry test with multiple layers

* fix typo

* merge branch 'main' into es_reconstruction_1d

* add reference for desingularization

* address changes from code review

* soften tolerances to fix macOS testing

* rework thresholds

* fix typo

* add unit test

* increase maxiters for single coverage test

* adjust comment

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

* update convergence test values

---------

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

examples/tree_1d_dgsem/elixir_shallowwater_multilayer_beach.jl

@@ -0,0 +1,126 @@
+
+using OrdinaryDiffEq
+using Trixi
+using TrixiShallowWater
+
+###############################################################################
+# Semidiscretization of the shallow water equations
+
+# By passing only a single value for rhos, the system recovers the standard shallow water equations
+equations = ShallowWaterMultiLayerEquations1D(gravity_constant = 9.812, rhos = 1.0)
+
+"""
+    initial_condition_beach(x, t, equations:: ShallowWaterMultiLayerEquationsWetDry1D)
+
+Initial condition to simulate a wave running towards a beach and crashing. Difficult test
+including both wetting and drying in the domain using slip wall boundary conditions.
+The bottom topography is altered to be differentiable on the domain [0,8] and
+differs from the reference below.
+
+The water height and speed functions used here, are adapted from the initial condition
+found in section 5.2 of the paper:
+  - Andreas Bollermann, Sebastian Noelle, Maria Lukáčová-Medvidová (2011)
+    Finite volume evolution Galerkin methods for the shallow water equations with dry beds\n
+    [DOI: 10.4208/cicp.220210.020710a](https://dx.doi.org/10.4208/cicp.220210.020710a)
+"""
+function initial_condition_beach(x, t, equations::ShallowWaterMultiLayerEquations1D)
+    D = 1
+    delta = 0.02
+    gamma = sqrt((3 * delta) / (4 * D))
+    x_a = sqrt((4 * D) / (3 * delta)) * acosh(sqrt(20))
+
+    f = D + 40 * delta * sech(gamma * (8 * x[1] - x_a))^2
+
+    # steep curved beach
+    b = 0.01 + 99 / 409600 * 4^x[1]
+
+    if x[1] >= 6
+        H = b
+        v = 0.0
+    else
+        H = f
+        v = sqrt(equations.gravity / D) * H
+    end
+
+    # It is mandatory to shift the water level at dry areas to make sure the water height h
+    # stays positive. The system would not be stable for h set to a hard 0 due to division by h in
+    # the computation of velocity, e.g., (h v) / h. Therefore, a small dry state threshold
+    # with a default value of 5*eps() ≈ 1e-15 in double precision, is set in the constructor above
+    # for the ShallowWaterMultiLayerEquations1D and added to the initial condition if h = 0.
+    # This default value can be changed within the constructor call depending on the simulation setup.
+    H = max(H, b + equations.threshold_limiter)
+    return prim2cons(SVector(H, v, b), equations)
+end
+
+initial_condition = initial_condition_beach
+boundary_condition = boundary_condition_slip_wall
+
+###############################################################################
+# Get the DG approximation space
+
+volume_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal)
+surface_flux = (FluxHydrostaticReconstruction(FluxPlusDissipation(flux_ersing_etal,
+                                                                  DissipationLocalLaxFriedrichs()),
+                                              hydrostatic_reconstruction_ersing_etal),
+                FluxHydrostaticReconstruction(flux_nonconservative_ersing_etal,
+                                              hydrostatic_reconstruction_ersing_etal))
+
+basis = LobattoLegendreBasis(3)
+
+indicator_sc = IndicatorHennemannGassnerShallowWater(equations, basis,
+                                                     alpha_max = 0.5,
+                                                     alpha_min = 0.001,
+                                                     alpha_smooth = true,
+                                                     variable = waterheight_pressure)
+volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
+                                                 volume_flux_dg = volume_flux,
+                                                 volume_flux_fv = surface_flux)
+
+solver = DGSEM(basis, surface_flux, volume_integral)
+
+###############################################################################
+# Create the TreeMesh for the domain [0, 8]
+
+coordinates_min = 0.0
+coordinates_max = 8.0
+
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 7,
+                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, 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 = (energy_kinetic,
+                                                                 energy_internal))
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(dt = 0.5,
+                                     save_initial_solution = true,
+                                     save_final_solution = true)
+
+callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution)
+
+stage_limiter! = PositivityPreservingLimiterShallowWater(variables = (Trixi.waterheight,))
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, SSPRK43(stage_limiter!);
+            ode_default_options()..., callback = callbacks);
+
+summary_callback() # print the timer summary

examples/tree_1d_dgsem/elixir_shallowwater_multilayer_dam_break_es_dry.jl

@@ -0,0 +1,124 @@
+
+using OrdinaryDiffEq
+using Trixi
+using TrixiShallowWater
+
+###############################################################################
+# Semidiscretization of the multilayer shallow water equations for a dam break
+# test over a dry domain with a discontinuous bottom topography function
+
+equations = ShallowWaterMultiLayerEquations1D(gravity_constant = 9.81,
+                                              rhos = (0.8, 0.85, 0.9, 0.95, 1.0))
+
+# Initial condition of a dam break with a discontinuous water heights and bottom topography.
+# To test the wet/dry functionality this test case considers a dam break over a dry domain with 
+# multiple layers.
+# 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)
+    # Set the discontinuity
+    if x[1] <= 10.0
+        H = [3.5, 3.0, 2.5, 2.0, 1.5]
+        b = 0.0
+    else
+        # Right side of the domain is dry
+        b = 0.5
+        H = [b, b, b, b, b]
+    end
+
+    # It is mandatory to shift the water level at dry areas to make sure the water height h
+    # stays positive. The system would not be stable for h set to a hard 0 due to division by h in
+    # the computation of velocity, e.g., (h v) / h. Therefore, a small dry state threshold
+    # with a default value of 5*eps() ≈ 1e-15 in double precision, is set in the constructor above
+    # for the ShallowWaterMultiLayerEquations1D and added to the initial condition if h = 0.
+    # This default value can be changed within the constructor call depending on the simulation setup.
+    for i in reverse(eachlayer(equations))
+        if i == nlayers(equations)
+            H[i] = max(H[i], b + equations.threshold_limiter)
+        else
+            H[i] = max(H[i], H[i + 1] + equations.threshold_limiter)
+        end
+    end
+
+    # Initialize zero velocity
+    v = zero(H)
+
+    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)
+surface_flux = (FluxHydrostaticReconstruction(FluxPlusDissipation(flux_ersing_etal,
+                                                                  DissipationLocalLaxFriedrichs()),
+                                              hydrostatic_reconstruction_ersing_etal),
+                FluxHydrostaticReconstruction(flux_nonconservative_ersing_etal,
+                                              hydrostatic_reconstruction_ersing_etal))
+
+basis = LobattoLegendreBasis(3)
+
+indicator_sc = IndicatorHennemannGassnerShallowWater(equations, basis,
+                                                     alpha_max = 0.5,
+                                                     alpha_min = 0.001,
+                                                     alpha_smooth = true,
+                                                     variable = waterheight_pressure)
+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 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, 1.0)
+ode = semidiscretize(semi, tspan)
+
+###############################################################################
+# Callbacks
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 500
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
+                                     save_analysis = false)
+
+stepsize_callback = StepsizeCallback(cfl = 0.2)
+
+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,
+                        stepsize_callback, save_solution)
+
+stage_limiter! = PositivityPreservingLimiterShallowWater(variables = (Trixi.waterheight,))
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, SSPRK43(stage_limiter!), dt = 1.0,
+            save_everystep = false, callback = callbacks, adaptive = false);
+summary_callback() # print the timer summary

examples/tree_1d_dgsem/elixir_shallowwater_multilayer_parabolic_bowl.jl

@@ -0,0 +1,124 @@
+
+using OrdinaryDiffEq
+using Trixi
+using TrixiShallowWater
+
+###############################################################################
+# Semidiscretization of the shallow water equations
+
+# By passing only a single value for rhos, the system recovers the standard shallow water equations
+equations = ShallowWaterMultiLayerEquations1D(gravity_constant = 9.81, rhos = 1.0)
+
+"""
+    initial_condition_parabolic_bowl(x, t, equations:: ShallowWaterMultiLayerEquations1D)
+
+Well-known initial condition to test the [`hydrostatic_reconstruction_ersing_etal`](@ref) and its
+wet-dry mechanics. This test has analytical solutions. The initial condition is defined by the
+analytical solution at time t=0. The bottom topography defines a bowl and the water level is given
+by an oscillating lake.
+
+The original test and its analytical solution in two dimensions were first presented in
+- William C. Thacker (1981)
+  Some exact solutions to the nonlinear shallow-water wave equations
+  [DOI: 10.1017/S0022112081001882](https://doi.org/10.1017/S0022112081001882).
+
+The particular setup below is taken from Section 6.2 of
+- Niklas Wintermeyer, Andrew R. Winters, Gregor J. Gassner and Timothy Warburton (2018)
+  An entropy stable discontinuous Galerkin method for the shallow water equations on
+  curvilinear meshes with wet/dry fronts accelerated by GPUs
+  [DOI: 10.1016/j.jcp.2018.08.038](https://doi.org/10.1016/j.jcp.2018.08.038).
+"""
+function initial_condition_parabolic_bowl(x, t,
+                                          equations::ShallowWaterMultiLayerEquations1D)
+    a = 1
+    h_0 = 0.1
+    sigma = 0.5
+    ω = sqrt(2 * equations.gravity * h_0) / a
+
+    v = -sigma * ω * sin(ω * t)
+
+    b = h_0 * x[1]^2 / a^2
+
+    H = sigma * h_0 / a^2 * (2 * x[1] * cos(ω * t) - sigma) + h_0
+
+    #= 
+    It is mandatory to shift the water level at dry areas to make sure the water height h
+    stays positive. The system would not be stable for h set to a hard 0 due to division by h in
+    the computation of velocity, e.g., (h v) / h. Therefore, a small dry state threshold
+    with a default value of 5*eps() ≈ 1e-15 in double precision, is set in the constructor above
+    for the ShallowWaterMultiLayerEquations1D and added to the initial condition if h = 0.
+    This default value can be changed within the constructor call depending on the simulation setup.
+    =#
+    H = max(H, b + equations.threshold_limiter)
+    return prim2cons(SVector(H, v, b), equations)
+end
+
+initial_condition = initial_condition_parabolic_bowl
+
+###############################################################################
+# Get the DG approximation space
+
+volume_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal)
+surface_flux = (FluxHydrostaticReconstruction(FluxPlusDissipation(flux_ersing_etal,
+                                                                  DissipationLocalLaxFriedrichs()),
+                                              hydrostatic_reconstruction_ersing_etal),
+                FluxHydrostaticReconstruction(flux_nonconservative_ersing_etal,
+                                              hydrostatic_reconstruction_ersing_etal))
+
+basis = LobattoLegendreBasis(5)
+
+indicator_sc = IndicatorHennemannGassnerShallowWater(equations, basis,
+                                                     alpha_max = 0.5,
+                                                     alpha_min = 0.001,
+                                                     alpha_smooth = true,
+                                                     variable = waterheight_pressure)
+volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
+                                                 volume_flux_dg = volume_flux,
+                                                 volume_flux_fv = surface_flux)
+
+solver = DGSEM(basis, surface_flux, volume_integral)
+
+###############################################################################
+# Create the TreeMesh for the domain [-2, 2]
+
+coordinates_min = -2.0
+coordinates_max = 2.0
+
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 6,
+                n_cells_max = 10_000)
+
+# 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 = (energy_kinetic,
+                                                                 energy_internal))
+
+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)
+
+stage_limiter! = PositivityPreservingLimiterShallowWater(variables = (Trixi.waterheight,))
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, SSPRK43(stage_limiter!);
+            ode_default_options()..., callback = callbacks);
+
+summary_callback() # print the timer summary