[WIP] Shallow-water exner equations with wetting/drying

examples/tree_1d_dgsem/elixir_shallowwater_exner_beach.jl

@@ -0,0 +1,123 @@
+
+using OrdinaryDiffEq
+using Trixi
+using TrixiShallowWater
+
+###############################################################################
+# Semidiscretization of the SWE-Exner equations
+
+# Equations with Grass model
+equations = ShallowWaterExnerEquations1D(gravity_constant = 9.812, H0 = 1.0,
+                                         rho_f = 0.5, rho_s = 1.0, porosity = 0.5,
+                                         friction = ManningFriction(n = 0.0),
+                                         sediment_model = GrassModel(A_g = 0.001),
+                                         threshold_desingularization = 1e-6)
+
+"""
+    initial_condition_beach(x, t, equations:: ShallowWaterExnerEquations1D)
+
+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á-Medvid’ová (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::ShallowWaterExnerEquations1D)
+    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
+    h_b = 0.01 + 99 / 409600 * 4^x[1]
+
+    if x[1] >= 6
+        H = h_b
+        v = 0.0
+    else
+        H = f
+        v = sqrt(equations.gravity / D) * H
+    end
+
+    H = max(H, h_b + equations.threshold_limiter)
+    return prim2cons(SVector(H, v, h_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 = (FluxPlusDissipation(flux_ersing_etal, DissipationLocalLaxFriedrichs()),
+                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 = 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)
+
+###############################################################################
+# 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 = 100
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
+                                     save_analysis = false)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(dt = 0.2,
+                                     save_initial_solution = true,
+                                     save_final_solution = true)
+
+stepsize_callback = StepsizeCallback(cfl = 0.3)
+
+callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution,
+                        stepsize_callback)
+
+stage_limiter! = PositivityPreservingLimiterShallowWater(variables = (Trixi.waterheight,))
+
+###############################################################################
+# run the simulation
+# use a Runge-Kutta method with CFL-based time step
+sol = solve(ode, SSPRK43(stage_limiter!);
+            ode_default_options()..., callback = callbacks, adaptive = false, dt = 1.0);
+
+summary_callback() # print the timer summary

examples/tree_1d_dgsem/elixir_shallowwater_exner_well_balanced_dry.jl

@@ -0,0 +1,107 @@
+
+using OrdinaryDiffEq
+using Trixi
+using TrixiShallowWater
+
+###############################################################################
+# Semidiscretization of the SWE-Exner equations with a discontinuous sediment bed 
+# to test well-balancedness
+
+# Equations with Grass model
+equations = ShallowWaterExnerEquations1D(gravity_constant = 9.81, H0 = 1.0,
+                                         rho_f = 0.5, rho_s = 1.0, porosity = 0.5,
+                                         friction = ManningFriction(n = 0.01),
+                                         sediment_model = GrassModel(A_g = 0.01))
+
+function initial_condition_steady_state(x, t,
+                                        equations::ShallowWaterExnerEquations1D)
+    hv = 0.0
+
+    # discontinuous sediment bed
+    if -2 < x[1] <= 0
+        h = 0.25
+        h_b = equations.threshold_limiter
+    else
+        h = equations.threshold_limiter
+        h_b = 0.5
+    end
+
+    return SVector(h, hv, h_b)
+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 = (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 = 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,
+                                    boundary_conditions = boundary_condition)
+
+###############################################################################
+# ODE solver
+tspan = (0.0, 2.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 = 5.0,
+                                     save_initial_solution = true,
+                                     save_final_solution = true)
+
+stepsize_callback = StepsizeCallback(cfl = 0.5)
+
+callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution,
+                        stepsize_callback)
+
+stage_limiter! = PositivityPreservingLimiterShallowWater(variables = (Trixi.waterheight,))
+
+###############################################################################
+# run the simulation
+# use a Runge-Kutta method with CFL-based time step
+sol = solve(ode, SSPRK43(stage_limiter!);
+            ode_default_options()..., callback = callbacks, adaptive = false, dt = 1.0);
+
+summary_callback() # print the timer summary

src/callbacks_stage/positivity_shallow_water.jl

@@ -82,7 +82,8 @@ function limiter_shallow_water!(u, variables::NTuple{N, Any},
                                 equations::Union{ShallowWaterEquationsWetDry1D,
                                                  ShallowWaterEquationsWetDry2D,
                                                  ShallowWaterMultiLayerEquations1D,
-                                                 ShallowWaterMultiLayerEquations2D},
+                                                 ShallowWaterMultiLayerEquations2D,
+                                                 ShallowWaterExnerEquations1D},
                                 solver, cache) where {N}
     variable = first(variables)
     remaining_variables = Base.tail(variables)
@@ -99,7 +100,8 @@ function limiter_shallow_water!(u, variables::Tuple{},
                                 equations::Union{ShallowWaterEquationsWetDry1D,
                                                  ShallowWaterEquationsWetDry2D,
                                                  ShallowWaterMultiLayerEquations1D,
-                                                 ShallowWaterMultiLayerEquations2D},
+                                                 ShallowWaterMultiLayerEquations2D,
+                                                 ShallowWaterExnerEquations1D},
                                 solver, cache)
     nothing
 end

src/callbacks_stage/positivity_shallow_water_dg1d.jl

@@ -171,4 +171,76 @@
 
     return nothing
 end
+
+# !!! warning "Experimental code"
+#     This is an experimental feature and may change in future releases.
+function limiter_shallow_water!(u, threshold::Real, variable,
+                                mesh::Trixi.AbstractMesh{1},
+                                equations::ShallowWaterExnerEquations1D,
+                                dg::DGSEM, cache)
+    @unpack weights = dg.basis
+
+    Trixi.@threaded for element in eachelement(dg, cache)
+        # determine minimum value
+        value_min = typemax(eltype(u))
+        for i in eachnode(dg)
+            u_node = get_node_vars(u, equations, dg, i, element)
+            value_min = min(value_min, variable(u_node, equations))
+        end
+
+        # detect if limiting is necessary
+        value_min < threshold - eps() || continue
+
+        # compute mean value
+        u_mean = zero(get_node_vars(u, equations, dg, 1, element))
+        for i in eachnode(dg)
+            u_node = get_node_vars(u, equations, dg, i, element)
+            u_mean += u_node * weights[i]
+        end
+
+        # note that the reference element is [-1,1]^ndims(dg), thus the weights sum to 2
+        u_mean = u_mean / 2^ndims(mesh)
+
+        # We compute the value directly with the mean values.
+        value_mean = variable(u_mean, equations)
+        theta = (value_mean - threshold) / (value_mean - value_min)
+
+        for i in eachnode(dg)
+            u_node = get_node_vars(u, equations, dg, i, element)
+            h_node = waterheight(u_node, equations)
+            h_mean = waterheight(u_mean, equations)
+
+            u[1, i, element] = theta * h_node + (1 - theta) * h_mean
+        end
+    end
+
+    # "Safety" application of the wet/dry thresholds over all the DG nodes
+    # on the current `element` after the limiting above in order to avoid dry nodes.
+    # If the value_mean < threshold before applying limiter, there
+    # could still be dry nodes afterwards due to logic of the limiting.
+    # Additionally, a velocity desingularization is applied to avoid numerical 
+    # problems near dry states.
+    Trixi.@threaded for element in eachelement(dg, cache)
+        for i in eachnode(dg)
+            u_node = get_node_vars(u, equations, dg, i, element)
+            h, hv, h_b = u_node
+
+            # Velocity desingularization
+            hv = h * (2.0 * h * hv) /
+                 (h^2 + max(h^2, equations.threshold_desingularization))
+
+            # Ensure positivity and zero velocity at dry states
+            if h <= threshold
+                h = threshold
+                hv = zero(eltype(u))
+            end
+
+            u_node = SVector(h, hv, h_b)
+
+            set_node_vars!(u, u_node, equations, dg, i, element)
+        end
+    end
+
+    return nothing
+end
 end # @muladd