From 68148a7b3b7262fa20a8977866bf11e3f7ea108d Mon Sep 17 00:00:00 2001 From: patrickersing Date: Fri, 23 Aug 2024 10:49:47 +0200 Subject: [PATCH] Add wet/dry limiters, modified shock-capturing and HR --- .../elixir_shallowwater_exner_beach.jl | 123 ++++++++++++++++++ ...ir_shallowwater_exner_well_balanced_dry.jl | 107 +++++++++++++++ .../positivity_shallow_water.jl | 6 +- .../positivity_shallow_water_dg1d.jl | 72 ++++++++++ src/equations/shallow_water_exner_1d.jl | 82 +++++++++++- src/solvers/indicators.jl | 3 +- src/solvers/indicators_1d.jl | 3 +- 7 files changed, 390 insertions(+), 6 deletions(-) create mode 100644 examples/tree_1d_dgsem/elixir_shallowwater_exner_beach.jl create mode 100644 examples/tree_1d_dgsem/elixir_shallowwater_exner_well_balanced_dry.jl diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_exner_beach.jl b/examples/tree_1d_dgsem/elixir_shallowwater_exner_beach.jl new file mode 100644 index 0000000..7b5957f --- /dev/null +++ b/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 diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_exner_well_balanced_dry.jl b/examples/tree_1d_dgsem/elixir_shallowwater_exner_well_balanced_dry.jl new file mode 100644 index 0000000..6919249 --- /dev/null +++ b/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 diff --git a/src/callbacks_stage/positivity_shallow_water.jl b/src/callbacks_stage/positivity_shallow_water.jl index 80a5307..71e44fc 100644 --- a/src/callbacks_stage/positivity_shallow_water.jl +++ b/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 diff --git a/src/callbacks_stage/positivity_shallow_water_dg1d.jl b/src/callbacks_stage/positivity_shallow_water_dg1d.jl index 14b77e4..54479cf 100644 --- a/src/callbacks_stage/positivity_shallow_water_dg1d.jl +++ b/src/callbacks_stage/positivity_shallow_water_dg1d.jl @@ -171,4 +171,76 @@ function limiter_shallow_water!(u, threshold::Real, variable, 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 diff --git a/src/equations/shallow_water_exner_1d.jl b/src/equations/shallow_water_exner_1d.jl index cf77759..7cd917d 100644 --- a/src/equations/shallow_water_exner_1d.jl +++ b/src/equations/shallow_water_exner_1d.jl @@ -150,6 +150,9 @@ struct ShallowWaterExnerEquations1D{RealT <: Real, rho_f::RealT # density of fluid rho_s::RealT # density of sediment r::RealT # density ratio + threshold_limiter::RealT # threshold for the positivity-limiter + threshold_desingularization::RealT # threshold for velocity desingularization + threshold_partially_wet::RealT # threshold to define partially wet elements end # Allow for flexibility to set the gravitational constant within an elixir depending on the @@ -159,12 +162,32 @@ end function ShallowWaterExnerEquations1D(; gravity_constant, H0 = zero(gravity_constant), friction = ManningFriction(n = 0.0), sediment_model, - porosity, rho_f, rho_s) + porosity, rho_f, rho_s, + threshold_limiter = nothing, + threshold_desingularization = nothing, + threshold_partially_wet = nothing) + RealT = typeof(gravity_constant) + # Precompute common expressions for the porosity and density ratio porosity_inv = inv(1 - porosity) r = rho_f / rho_s + + # Set default values for thresholds + if threshold_limiter === nothing + threshold_limiter = 5 * eps(RealT) + end + if threshold_desingularization === nothing + threshold_desingularization = default_threshold_desingularization(RealT) + end + if threshold_partially_wet === nothing + threshold_partially_wet = default_threshold_partially_wet(RealT) + end + return ShallowWaterExnerEquations1D(gravity_constant, H0, friction, sediment_model, - porosity_inv, rho_f, rho_s, r) + porosity_inv, rho_f, rho_s, r, + threshold_limiter, + threshold_desingularization, + threshold_partially_wet) end Trixi.have_nonconservative_terms(::ShallowWaterExnerEquations1D) = True() @@ -407,6 +430,61 @@ To obtain an entropy stable formulation the `surface_flux` can be set as return SVector(f1, f2, f3) end +# TODO: This reconstruction is not well-balanced if used with LLF-Dissipation due to the discontinuous +# sediment layer. Instead this currently applies the Audusse reconstruction, but without accounting +# for the nonconservative terms. We should either find find a dissipation that is WB or use the Audusse +# reconstruction with a corresponding nonconservative term. +@inline function hydrostatic_reconstruction_ersing_etal(u_ll, u_rr, + equations::ShallowWaterExnerEquations1D) + # Unpack waterheight and bottom topographies + h_ll = waterheight(u_ll, equations) + h_rr = waterheight(u_rr, equations) + h_b_ll = u_ll[end] + h_b_rr = u_rr[end] + + # Get the velocities on either side + v_ll = velocity(u_ll, equations) + v_rr = velocity(u_rr, equations) + + threshold = equations.threshold_limiter + + # Calculate total layer heights + H_ll = waterheight(cons2prim(u_ll, equations), equations) + H_rr = waterheight(cons2prim(u_rr, equations), equations) + + # Discontinuous reconstruction of the bottom topography - This is not WB, because the discontinuous + # reconstruction introduces non-zero dissipation at steady states + #h_b_ll_star = min(H_ll, max(h_b_rr, h_b_ll)) + #h_b_rr_star = min(H_rr, max(h_b_rr, h_b_ll)) + + # Audusse style reconstruction is needed for WB - note that the NC terms from the reconstruction are + # not accounted for in nonconservative_flux_ersing_etal + h_b_ll_star = max(h_b_rr, h_b_ll) + h_b_rr_star = max(h_b_rr, h_b_ll) + + # Calculate reconstructed total layer heights + H_ll_star = max(H_ll, h_b_ll_star) + H_rr_star = max(H_rr, h_b_rr_star) + + # Reconstruct the waterheights + h_ll_star = max(H_ll_star - h_b_ll_star, threshold) + h_rr_star = max(H_rr_star - h_b_rr_star, threshold) + + # Ensure zero velocities at dry states + if h_ll_star <= threshold + v_ll = zero(eltype(u_ll)) + end + if h_rr_star <= threshold + v_rr = zero(eltype(u_rr)) + end + + # Reconstructed states + u_ll_star = SVector(h_ll_star, (h_ll_star * v_ll), h_b_ll_star) + u_rr_star = SVector(h_rr_star, (h_rr_star * v_rr), h_b_rr_star) + + return SVector(u_ll_star, u_rr_star) +end + """ dissipation_roe(u_ll, u_rr, orientation_or_normal_direction, equations::ShallowWaterExnerEquations1D) diff --git a/src/solvers/indicators.jl b/src/solvers/indicators.jl index 884d19d..fb77321 100644 --- a/src/solvers/indicators.jl +++ b/src/solvers/indicators.jl @@ -41,7 +41,8 @@ end # It modifies the shock-capturing indicator to use full FV method in dry elements # or partially dry elements containing a wet/dry transition. function IndicatorHennemannGassnerShallowWater(equations::Union{Trixi.AbstractShallowWaterEquations, - AbstractShallowWaterMultiLayerEquations}, + AbstractShallowWaterMultiLayerEquations, + ShallowWaterExnerEquations1D}, basis; alpha_max = 0.5, alpha_min = 0.001, diff --git a/src/solvers/indicators_1d.jl b/src/solvers/indicators_1d.jl index 5931e0e..06dd3db 100644 --- a/src/solvers/indicators_1d.jl +++ b/src/solvers/indicators_1d.jl @@ -12,7 +12,8 @@ function (indicator_hg::IndicatorHennemannGassnerShallowWater)(u::AbstractArray{ 3}, mesh, equations::Union{ShallowWaterEquationsWetDry1D, - ShallowWaterMultiLayerEquations1D}, + ShallowWaterMultiLayerEquations1D, + ShallowWaterExnerEquations1D}, dg::DGSEM, cache; kwargs...) @unpack alpha_max, alpha_min, alpha_smooth, variable = indicator_hg