diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_multilayer_beach.jl b/examples/tree_1d_dgsem/elixir_shallowwater_multilayer_beach.jl new file mode 100644 index 0000000..ec3234f --- /dev/null +++ b/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 diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_multilayer_dam_break_es_dry.jl b/examples/tree_1d_dgsem/elixir_shallowwater_multilayer_dam_break_es_dry.jl new file mode 100644 index 0000000..e474494 --- /dev/null +++ b/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 diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_multilayer_parabolic_bowl.jl b/examples/tree_1d_dgsem/elixir_shallowwater_multilayer_parabolic_bowl.jl new file mode 100644 index 0000000..d48557e --- /dev/null +++ b/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 diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_multilayer_well_balanced_wet_dry.jl b/examples/tree_1d_dgsem/elixir_shallowwater_multilayer_well_balanced_wet_dry.jl new file mode 100644 index 0000000..b82ddba --- /dev/null +++ b/examples/tree_1d_dgsem/elixir_shallowwater_multilayer_well_balanced_wet_dry.jl @@ -0,0 +1,171 @@ + +using OrdinaryDiffEq +using Trixi +using TrixiShallowWater + +############################################################################### +# Semidiscretization of the two-layer shallow water equations to test well-balancedness + +equations = ShallowWaterMultiLayerEquations1D(gravity_constant = 1.0, H0 = 2.0, + rhos = (1.0, 3.0)) + +""" + initial_condition_twolayer_well_balanced_wet_dry(perturbation, equations:: ShallowWaterMultiLayerEquations1D) + initial_condition_twolayer_well_balanced_wet_dry(x, t, equations:: ShallowWaterMultiLayerEquations1D) + +Initial condition with a complex (discontinuous) bottom topography to test well-balancedness for a +two-layer shallow water system with dry states if `perturbation` is set to `false`. +Additionally, it is possible to set a perturbation in the lower layer to test perturbations from the +lake-at-rest condition by setting the `perturbation` variable to `true`. + +The initial condition is taken from the paper: + - S. Martínez-Aranda, A. Ramos-Pérez, P. García-Navarro (2020) + A 1D shallow-flow model for two-layer flows based on FORCE scheme with wet–dry treatment\n + [DOI:10.2166/hydro.2020.002](https://doi.org/10.2166/hydro.2020.002) +""" +function initial_condition_twolayer_well_balanced_wet_dry(perturbation::Bool, + equations::ShallowWaterMultiLayerEquations1D) + return function initial_condition_twolayer_well_balanced_wet_dry(x, t, equations) + # Set interface height for the upper layer + H_upper = 2.0 + + # Set interface height for the lower layer + H_lower = 1.0 + if perturbation == true + if x[1] <= 5.0 + nothing + elseif x[1] <= 15.0 + H_lower = 0.75 + elseif x[1] <= 35.0 + H_lower = 1.5 - (x[1] - 15.0) / 20.0 * 0.5 + elseif x[1] <= 40.0 + H_lower = 1.0 + (x[1] - 35.0) / 5.0 * 0.5 + elseif x[1] <= 55.0 + H_lower = 1.0 + elseif x[1] <= 70.0 + H_lower = 1.6 - (x[1] - 55.0) / 15.0 * 0.4 + elseif x[1] <= 75.0 + nothing + elseif x[1] <= 85.0 + H_lower = 1.0 - (x[1] - 75.0) / 10.0 * 0.3 + elseif x[1] <= 95.0 + H_lower = 1.1 + else + H_lower = 1.5 + end + end + + # Set the bottom topography + if x[1] <= 5.0 + b = 2.5 + elseif x[1] <= 20.0 + b = 0.5 + elseif x[1] <= 40.0 + b = 0.1 + (x[1] - 20.0) / 20.0 * 1.4 + elseif x[1] <= 60.0 + b = 1.0 - (x[1] - 40.0) / 20.0 * 1.3 + elseif x[1] <= 70.0 + b = 0.7 + elseif x[1] <= 75.0 + b = 2.2 + elseif x[1] <= 95.0 + b = 0.4 + else + b = 1.5 + end + + # Set zero initial velocity + v1_upper = zero(H_upper) + v1_lower = zero(H_upper) + + #= + 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_lower = max(H_lower, b + equations.threshold_limiter) + H_upper = max(H_upper, H_lower + equations.threshold_limiter) + + return prim2cons(SVector(H_upper, H_lower, v1_upper, v1_lower, b), equations) + end +end + +# Default value for the perturbation +perturbation = false + +initial_condition = initial_condition_twolayer_well_balanced_wet_dry(perturbation, + equations) + +############################################################################### +# 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 periodic mesh + +# The domain of interest [0.0, 100.0] is extended towards -28.0 to match the positions of the +# disctontinuities with TreeMesh. +coordinates_min = -28.0 +coordinates_max = 100.0 +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 7, + 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 = 0.5) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 1000, + save_initial_solution = true, + save_final_solution = true) + +stage_limiter! = PositivityPreservingLimiterShallowWater(variables = (Trixi.waterheight,)) + +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, + stepsize_callback) + +############################################################################### +# 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/TrixiShallowWater.jl b/src/TrixiShallowWater.jl index c3b0147..5522fc5 100644 --- a/src/TrixiShallowWater.jl +++ b/src/TrixiShallowWater.jl @@ -24,7 +24,7 @@ export ShallowWaterEquationsWetDry1D, ShallowWaterEquationsWetDry2D, export hydrostatic_reconstruction_chen_noelle, flux_nonconservative_chen_noelle, min_max_speed_chen_noelle, flux_hll_chen_noelle, - flux_ersing_etal, flux_es_ersing_etal + flux_ersing_etal, flux_es_ersing_etal, hydrostatic_reconstruction_ersing_etal export nlayers, eachlayer diff --git a/src/callbacks_stage/positivity_shallow_water.jl b/src/callbacks_stage/positivity_shallow_water.jl index c2c870c..8ed1630 100644 --- a/src/callbacks_stage/positivity_shallow_water.jl +++ b/src/callbacks_stage/positivity_shallow_water.jl @@ -10,10 +10,11 @@ The limiter is specifically designed for the shallow water equations. It is applied to all scalar `variables` in their given order -using the defined `threshold_limiter` from the [`ShallowWaterEquationsWetDry1D`](@ref) struct -or the [`ShallowWaterEquationsWetDry2D`](@ref) struct to determine the minimal acceptable values. +using the defined `threshold_limiter` from the equations struct +(e.g. in [`ShallowWaterEquationsWetDry1D`](@ref)) to determine the minimal acceptable values. The order of the `variables` is important and might have a strong influence -on the robustness. +on the robustness. The limiter is available for the [`ShallowWaterEquationsWetDry1D`](@ref), +[`ShallowWaterEquationsWetDry2D`](@ref), and [`ShallowWaterMultiLayerEquations1D`](@ref). As opposed to the standard version of the [`Trixi.PositivityPreservingLimiterZhangShu`](@extref), nodes with a water height below the `threshold_limiter` are treated in a special way. @@ -28,11 +29,28 @@ the `threshold_limiter` is applied again on all the DG nodes in order to avoid w In the case where the cell mean value is below the threshold before applying the limiter, there could still be dry nodes afterwards due to the logic of the limiter. +For the [`ShallowWaterMultiLayerEquations1D`](@ref) the implementation differs. In this case the +positivity limiter is applied layerwise and only the waterheight `h` is limited within each layer. +Furthermore, a velocity desingularization is applied after the limiting to avoid numerical problems +near dry states. Details about the desingularization strategy can be found in Section 2.2 of the paper +- A. Kurganov, G. Petrova (2007) + A second-order well-balanced positivity preserving central-upwind scheme for the Saint-Venant system + [doi: 10.4310/CMS.2007.v5.n1.a6](https://dx.doi.org/10.4310/CMS.2007.v5.n1.a6) + This fully-discrete positivity-preserving limiter is based on the work of - Zhang, Shu (2011) Maximum-principle-satisfying and positivity-preserving high-order schemes for conservation laws: survey and new developments [doi: 10.1098/rspa.2011.0153](https://doi.org/10.1098/rspa.2011.0153) +The specific implementation for the [`ShallowWaterMultiLayerEquations1D](@ref) is based on the work of +- Y. Xing, X. Zhang (2013) + Positivity-preserving well-balanced discontinuous Galerkin methods for the shallow water equations + on unstructured triangular meshes + [doi: 10.1007/s10915-012-9644-4](https://doi.org/10.1007/s10915-012-9644-4) + + +!!! warning "Experimental code" + This is an experimental feature and may change in future releases. """ struct PositivityPreservingLimiterShallowWater{N, Variables <: NTuple{N, Any}} variables::Variables @@ -62,7 +80,8 @@ end function limiter_shallow_water!(u, variables::NTuple{N, Any}, mesh, equations::Union{ShallowWaterEquationsWetDry1D, - ShallowWaterEquationsWetDry2D}, + ShallowWaterEquationsWetDry2D, + ShallowWaterMultiLayerEquations1D}, solver, cache) where {N} variable = first(variables) remaining_variables = Base.tail(variables) @@ -77,7 +96,8 @@ end function limiter_shallow_water!(u, variables::Tuple{}, mesh, equations::Union{ShallowWaterEquationsWetDry1D, - ShallowWaterEquationsWetDry2D}, + ShallowWaterEquationsWetDry2D, + ShallowWaterMultiLayerEquations1D}, 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 3494662..14b77e4 100644 --- a/src/callbacks_stage/positivity_shallow_water_dg1d.jl +++ b/src/callbacks_stage/positivity_shallow_water_dg1d.jl @@ -5,6 +5,8 @@ @muladd begin #! format: noindent +# !!! 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::ShallowWaterEquationsWetDry1D, @@ -84,4 +86,89 @@ function limiter_shallow_water!(u, threshold::Real, variable, return nothing end + +# Note that for the `ShallowWaterMultiLayerEquations1D` only the waterheight `h` is limited in +# each layer. Furthermore, a velocity desingularization is applied after the limiting to avoid +# numerical problems near dry states. +# +# !!! 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::ShallowWaterMultiLayerEquations1D, + dg::DGSEM, cache) + @unpack weights = dg.basis + + Trixi.@threaded for element in eachelement(dg, cache) + # Limit layerwise + for m in eachlayer(equations) + # 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)[m]) + 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. + # The waterheight `h` is limited independently in each layer. + value_mean = variable(u_mean, equations)[m] + 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)[m] + h_mean = waterheight(u_mean, equations)[m] + + u[m, i, element] = theta * h_node + (1 - theta) * h_mean + end + 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 = MVector(waterheight(u_node, equations)) + hv = MVector(momentum(u_node, equations)) + b = u_node[end] + + # Velocity desingularization + hv = h .* (2.0 * h .* hv) ./ + (h .^ 2 .+ max.(h .^ 2, equations.threshold_desingularization)) + + for i in eachlayer(equations) + # Ensure positivity and zero velocity at dry states + if h[i] <= threshold + h[i] = threshold + hv[i] = zero(eltype(u)) + end + end + + u_node = SVector(h..., hv..., b) + + set_node_vars!(u, u_node, equations, dg, i, element) + end + end + + return nothing +end end # @muladd diff --git a/src/callbacks_stage/positivity_shallow_water_dg2d.jl b/src/callbacks_stage/positivity_shallow_water_dg2d.jl index 812f453..69c47e1 100644 --- a/src/callbacks_stage/positivity_shallow_water_dg2d.jl +++ b/src/callbacks_stage/positivity_shallow_water_dg2d.jl @@ -5,6 +5,8 @@ @muladd begin #! format: noindent +# !!! 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{2}, equations::ShallowWaterEquationsWetDry2D, dg::DGSEM, diff --git a/src/equations/equations.jl b/src/equations/equations.jl index 89713b9..b2e69a6 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -42,4 +42,13 @@ Retrieve the number of layers from an equation instance of the `AbstractShallowW } NLAYERS end + +# TODO: Add suitable default thresholds for Float32 +# Provide default thresholds dependent on number format (Currently default thresholds are only provided +# for Float64) +default_threshold_partially_wet(::Type{Float64}) = 1e-4 +default_threshold_partially_wet(catchall) = throw(ArgumentError("threshold_partially_wet must be provided for non-Float64 types")) + +default_threshold_desingularization(::Type{Float64}) = 1e-10 +default_threshold_desingularization(catchall) = throw(ArgumentError("threshold_desingularization must be provided for non-Float64 types")) end # @muladd diff --git a/src/equations/shallow_water_multilayer_1d.jl b/src/equations/shallow_water_multilayer_1d.jl index 138768f..20067e6 100644 --- a/src/equations/shallow_water_multilayer_1d.jl +++ b/src/equations/shallow_water_multilayer_1d.jl @@ -30,6 +30,16 @@ nonconservative term, which has some benefits for the design of well-balanced sc The additional quantity ``H_0`` is also available to store a reference value for the total water height that is useful to set initial conditions or test the "lake-at-rest" well-balancedness. +Also, there are two thresholds which prevent numerical problems as well as instabilities. The limiters are +used in [`PositivityPreservingLimiterShallowWater`](@ref) on the water height. `threshold_limiter` +acts as a (small) shift on the initial condition and cutoff before the next time step, whereas +`threshold_desingularization` is used in the velocity desingularization. A third +`threshold_partially_wet` is applied on the water height to define "partially wet" elements in +[`IndicatorHennemannGassnerShallowWater`](@ref), that are then calculated with a pure FV method to +ensure well-balancedness. For `Float64` no threshold needs to be passed, as default values are +defined within the struct. For other number formats `threshold_desingularization` and `threshold_partially_wet` +must be provided. + The bottom topography function ``b(x)`` is set inside the initial condition routine for a particular problem setup. @@ -54,10 +64,16 @@ struct ShallowWaterMultiLayerEquations1D{NVARS, NLAYERS, RealT <: Real} <: AbstractShallowWaterMultiLayerEquations{1, NVARS, NLAYERS} gravity::RealT # gravitational constant H0::RealT # constant "lake-at-rest" total water height + 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 rhos::SVector{NLAYERS, RealT} # Vector of layer densities function ShallowWaterMultiLayerEquations1D{NVARS, NLAYERS, RealT}(gravity::RealT, H0::RealT, + threshold_limiter::RealT, + threshold_desingularization::RealT, + threshold_partially_wet::RealT, rhos::SVector{NLAYERS, RealT}) where { NVARS, @@ -70,7 +86,8 @@ struct ShallowWaterMultiLayerEquations1D{NVARS, NLAYERS, RealT <: Real} <: throw(ArgumentError("densities must be in increasing order (rhos[1] < rhos[2] < ... < rhos[NLAYERS])")) min(rhos...) > 0 || throw(ArgumentError("densities must be positive")) - new(gravity, H0, rhos) + new(gravity, H0, threshold_limiter, threshold_desingularization, + threshold_partially_wet, rhos) end end @@ -79,19 +96,38 @@ end # The reference total water height H0 defaults to 0.0 but is used for the "lake-at-rest" # well-balancedness test cases. function ShallowWaterMultiLayerEquations1D(; gravity_constant, - H0 = zero(gravity_constant), rhos) + H0 = zero(gravity_constant), + threshold_limiter = nothing, + threshold_desingularization = nothing, + threshold_partially_wet = nothing, + rhos) # Promote all variables to a common type _rhos = promote(rhos...) RealT = promote_type(eltype(_rhos), eltype(gravity_constant), eltype(H0)) __rhos = SVector(map(RealT, _rhos)) + # 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 + # Extract number of layers and variables NLAYERS = length(rhos) NVARS = 2 * NLAYERS + 1 return ShallowWaterMultiLayerEquations1D{NVARS, NLAYERS, RealT}(gravity_constant, - H0, __rhos) + H0, + threshold_limiter, + threshold_desingularization, + threshold_partially_wet, + __rhos) end @inline function Base.real(::ShallowWaterMultiLayerEquations1D{NVARS, NLAYERS, RealT}) where { @@ -113,13 +149,13 @@ function Trixi.varnames(::typeof(cons2cons), return (waterheight..., momentum..., "b") end -# We use the interface heights, H = ∑h + b as primitive variables for easier visualization and setting initial +# We use the total layer heights, H = ∑h + b as primitive variables for easier visualization and setting initial # conditions function Trixi.varnames(::typeof(cons2prim), equations::ShallowWaterMultiLayerEquations1D) - interface_height = ntuple(n -> "H" * string(n), Val(nlayers(equations))) + total_layer_height = ntuple(n -> "H" * string(n), Val(nlayers(equations))) velocity = ntuple(n -> "v" * string(n) * "_1", Val(nlayers(equations))) - return (interface_height..., velocity..., "b") + return (total_layer_height..., velocity..., "b") end # Set initial conditions at physical location `x` for time `t` @@ -358,6 +394,79 @@ In the two-layer setting this combination is equivalent to the fluxes in: return SVector(f) end +""" + hydrostatic_reconstruction_ersing_etal(u_ll, u_rr, equations::ShallowWaterMultiLayerEquations1D) + +A particular hydrostatic reconstruction of the water height and bottom topography to +guarantee well-balancedness in the presence of wet/dry transitions and entropy stability for the +[`ShallowWaterMultiLayerEquations1D`](@ref). +The reconstructed solution states `u_ll_star` and `u_rr_star` are used to evaluate the +surface numerical flux at the interface. The key idea is a piecewise linear reconstruction of the +bottom topography and water height interfaces using subcells, where the bottom topography is allowed +to be discontinuous. +Use in combination with the generic numerical flux routine [`Trixi.FluxHydrostaticReconstruction`](@extref). + +!!! warning "Experimental code" + This is an experimental feature and may change in future releases. +""" +@inline function hydrostatic_reconstruction_ersing_etal(u_ll, u_rr, + equations::ShallowWaterMultiLayerEquations1D) + # Unpack waterheight and bottom topographies + h_ll = waterheight(u_ll, equations) + h_rr = waterheight(u_rr, equations) + b_ll = u_ll[end] + b_rr = u_rr[end] + + # Get the velocities on either side + v_ll = MVector(velocity(u_ll, equations)) + v_rr = MVector(velocity(u_rr, equations)) + + threshold = equations.threshold_limiter + + # Ensure zero velocity at dry states + for i in eachlayer(equations) + if h_ll[i] <= threshold + v_ll[i] = zero(eltype(u_ll)) + end + if h_rr[i] <= threshold + v_rr[i] = zero(eltype(u_rr)) + end + end + + # 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 + b_ll_star = min(H_ll[1], max(b_rr, b_ll)) + b_rr_star = min(H_rr[1], max(b_rr, b_ll)) + + # Calculate reconstructed total layer heights + H_ll_star = max.(H_ll, b_ll_star) + H_rr_star = max.(H_rr, b_rr_star) + + # Initialize reconstructed waterheights + h_ll_star = zero(MVector{nlayers(equations), real(equations)}) + h_rr_star = zero(MVector{nlayers(equations), real(equations)}) + + # Reconstruct the waterheights + for i in eachlayer(equations) + if i == nlayers(equations) # The lowest layer is measured from the bottom topography + h_ll_star[i] = max(H_ll_star[i] - b_ll_star, threshold) + h_rr_star[i] = max(H_rr_star[i] - b_rr_star, threshold) + else + h_ll_star[i] = max(H_ll_star[i] - H_ll_star[i + 1], threshold) + h_rr_star[i] = max(H_rr_star[i] - H_rr_star[i + 1], threshold) + end + end + + # Reconstructed states + u_ll_star = SVector(h_ll_star..., (h_ll_star .* v_ll)..., b_ll_star) + u_rr_star = SVector(h_rr_star..., (h_rr_star .* v_rr)..., b_rr_star) + + return SVector(u_ll_star, u_rr_star) +end + # Specialized `DissipationLocalLaxFriedrichs` to avoid spurious dissipation in the bottom # topography @inline function (dissipation::DissipationLocalLaxFriedrichs)(u_ll, u_rr, @@ -401,7 +510,7 @@ end c_ll = sqrt(equations.gravity * sum(h_ll)) c_rr = sqrt(equations.gravity * sum(h_rr)) - return (max(abs(v_m_ll) + c_ll, abs(v_m_rr) + c_rr)) + return (max(abs(v_m_ll), abs(v_m_rr)) + max(c_ll, c_rr)) end # Convert conservative variables to primitive @@ -410,7 +519,7 @@ end h = waterheight(u, equations) b = u[end] - # Initialize interface height + # Initialize total layer height H = MVector{nlayers(equations), real(equations)}(undef) for i in reverse(eachlayer(equations)) if i == nlayers(equations) @@ -426,18 +535,20 @@ end # Convert primitive to conservative variables @inline function Trixi.prim2cons(prim, equations::ShallowWaterMultiLayerEquations1D) - H = prim[1:nlayers(equations)] - v = prim[(nlayers(equations) + 1):(2 * nlayers(equations))] + # To extract the total layer height and velocity we reuse the waterheight and momentum functions + # from the conservative variables. + H = waterheight(prim, equations) # For primitive variables this extracts the total layer height + v = momentum(prim, equations) # For primitive variables this extracts the velocity b = prim[end] # Calculate waterheight h = MVector{nlayers(equations), real(equations)}(undef) for i in eachlayer(equations) if i < nlayers(equations) - setindex!(h, H[i] - H[i + 1], i) + h[i] = H[i] - H[i + 1] else # The lowest layer is measured from the bottom topography - setindex!(h, H[i] - b, i) + h[i] = H[i] - b end end @@ -463,9 +574,9 @@ end # Calculate entropy variables in each layer for i in eachlayer(equations) - # Compute w1[i] = ρ[i] * g * (b + ∑h[k] + ∑σ[k] * h[k]), where σ[k] = ρ[k] / ρ[i] denotes - # the density ratio of different layers - w1 = equations.rhos[i] * (g * b - 0.5 * v[i]^2) + # Compute w1[i] = ρ[i] * g * (b + ∑h[k] + ∑σ[k] * h[k]) - 0.5 * ρ[i] * v[i]^2, + # where σ[k] = ρ[k] / ρ[i] denotes the density ratio of different layers + w1 = equations.rhos[i] * (g * b) - 0.5 * equations.rhos[i] * v[i]^2 for j in eachlayer(equations) if j < i w1 += equations.rhos[i] * g * @@ -539,14 +650,28 @@ end return energy_total(u, equations) - energy_kinetic(u, equations) end -# Calculate the error for the "lake-at-rest" test case where H = ∑h+b should -# be a constant value over time +@inline function Trixi.waterheight_pressure(u, + equations::ShallowWaterMultiLayerEquations1D) + h = waterheight(u, equations) + + return 0.5 * equations.gravity * sum(h .^ 3) +end + +# Calculate the error for the "lake-at-rest" test case where H = ∑h + b should +# be a constant value over time. +# Note, assumes there is a single reference water height `H0` with which to compare. @inline function Trixi.lake_at_rest_error(u, equations::ShallowWaterMultiLayerEquations1D) h = waterheight(u, equations) b = u[end] - return abs(equations.H0 - (sum(h) + b)) + # For well-balancedness testing with possible wet/dry regions the reference + # water height `H0` accounts for the possibility that the bottom topography + # can emerge out of the water as well as for the threshold offset to avoid + # division by a "hard" zero water heights as well. + H0_wet_dry = max(equations.H0, b + equations.threshold_limiter) + + return abs(H0_wet_dry - (sum(h) + b)) end # Helper function to set the layer values in the flux computation diff --git a/src/equations/shallow_water_wet_dry_1d.jl b/src/equations/shallow_water_wet_dry_1d.jl index bb5e726..8cca15f 100644 --- a/src/equations/shallow_water_wet_dry_1d.jl +++ b/src/equations/shallow_water_wet_dry_1d.jl @@ -28,7 +28,11 @@ Also, there are two thresholds which prevent numerical problems as well as insta have to be passed, as default values are defined within the struct. The first one, `threshold_limiter`, is used in [`PositivityPreservingLimiterShallowWater`](@ref) on the water height, as a (small) shift on the initial condition and cutoff before the next time step. The second one, `threshold_wet`, is applied on the water height to -define when the flow is "wet" before calculating the numerical flux. +define when the flow is "wet" before calculating the numerical flux. A third +`threshold_partially_wet` is applied on the water height to define "partially wet" elements in +[`IndicatorHennemannGassnerShallowWater`](@ref), that are then calculated with a pure FV method to +ensure well-balancedness. For `Float64` no threshold needs to be passed, as default values are +defined within the struct. For other number formats `threshold_partially_wet` must be provided. The bottom topography function ``b(x)`` is set inside the initial condition routine for a particular problem setup. To test the conservative form of the SWE one can set the bottom topography @@ -62,6 +66,10 @@ struct ShallowWaterEquationsWetDry1D{RealT <: Real} <: # before calculating the numerical flux. # Default is 5*eps() which in double precision is ≈1e-15. threshold_wet::RealT + # `threshold_partially_wet` used in `IndicatorHennemannGassnerShallowWater` on the water height + # to define "partially wet" elements. Those elements are calculated with a pure FV method to + # ensure well-balancedness. Default in double precision is 1e-4. + threshold_partially_wet::RealT # Standard shallow water equations for dispatch on Trixi.jl functions basic_swe::ShallowWaterEquations1D{RealT} end @@ -73,7 +81,8 @@ end # Strict default values for thresholds that performed well in many numerical experiments function ShallowWaterEquationsWetDry1D(; gravity_constant, H0 = zero(gravity_constant), threshold_limiter = nothing, - threshold_wet = nothing) + threshold_wet = nothing, + threshold_partially_wet = nothing) T = promote_type(typeof(gravity_constant), typeof(H0)) if threshold_limiter === nothing threshold_limiter = 500 * eps(T) @@ -81,13 +90,16 @@ function ShallowWaterEquationsWetDry1D(; gravity_constant, H0 = zero(gravity_con if threshold_wet === nothing threshold_wet = 5 * eps(T) end + if threshold_partially_wet === nothing + threshold_partially_wet = default_threshold_partially_wet(T) + end # Construct the standard SWE for dispatch. Even though the `basic_swe` already store the # gravity constant and the total water height, we store an extra copy in # `ShallowWaterEquationsWetDry1D` for convenience. basic_swe = ShallowWaterEquations1D(gravity_constant = gravity_constant, H0 = H0) ShallowWaterEquationsWetDry1D(gravity_constant, H0, threshold_limiter, - threshold_wet, basic_swe) + threshold_wet, threshold_partially_wet, basic_swe) end Trixi.have_nonconservative_terms(::ShallowWaterEquationsWetDry1D) = True() diff --git a/src/equations/shallow_water_wet_dry_2d.jl b/src/equations/shallow_water_wet_dry_2d.jl index 2d59294..363fceb 100644 --- a/src/equations/shallow_water_wet_dry_2d.jl +++ b/src/equations/shallow_water_wet_dry_2d.jl @@ -31,7 +31,11 @@ Also, there are two thresholds which prevent numerical problems as well as insta have to be passed, as default values are defined within the struct. The first one, `threshold_limiter`, is used in [`PositivityPreservingLimiterShallowWater`](@ref) on the water height, as a (small) shift on the initial condition and cutoff before the next time step. The second one, `threshold_wet`, is applied on the water height to -define when the flow is "wet" before calculating the numerical flux. +define when the flow is "wet" before calculating the numerical flux. A third +`threshold_partially_wet` is applied on the water height to define "partially wet" elements in +[`IndicatorHennemannGassnerShallowWater`](@ref), that are then calculated with a pure FV method to +ensure well-balancedness. For `Float64` no threshold needs to be passed, as default values are +defined within the struct. For other number formats `threshold_partially_wet` must be provided. The bottom topography function ``b(x,y)`` is set inside the initial condition routine for a particular problem setup. To test the conservative form of the SWE one can set the bottom topography @@ -65,6 +69,10 @@ struct ShallowWaterEquationsWetDry2D{RealT <: Real} <: # before calculating the numerical flux. # Default is 5*eps() which in double precision is ≈1e-15. threshold_wet::RealT + # `threshold_partially_wet` used in `IndicatorHennemannGassnerShallowWater` on the water height + # to define "partially wet" elements. Those elements are calculated with a pure FV method to + # ensure well-balancedness. Default in double precision is 1e-4. + threshold_partially_wet::RealT # Standard shallow water equations for dispatch on Trixi.jl functions basic_swe::ShallowWaterEquations2D{RealT} end @@ -76,7 +84,8 @@ end # Strict default values for thresholds that performed well in many numerical experiments function ShallowWaterEquationsWetDry2D(; gravity_constant, H0 = zero(gravity_constant), threshold_limiter = nothing, - threshold_wet = nothing) + threshold_wet = nothing, + threshold_partially_wet = nothing) T = promote_type(typeof(gravity_constant), typeof(H0)) if threshold_limiter === nothing threshold_limiter = 500 * eps(T) @@ -84,13 +93,16 @@ function ShallowWaterEquationsWetDry2D(; gravity_constant, H0 = zero(gravity_con if threshold_wet === nothing threshold_wet = 5 * eps(T) end + if threshold_partially_wet === nothing + threshold_partially_wet = default_threshold_partially_wet(T) + end # Construct the standard SWE for dispatch. Even though the `basic_swe` already store the # gravity constant and the total water height, we store an extra copy in # `ShallowWaterEquationsWetDry2D` for convenience. basic_swe = ShallowWaterEquations2D(gravity_constant = gravity_constant, H0 = H0) ShallowWaterEquationsWetDry2D(gravity_constant, H0, threshold_limiter, - threshold_wet, basic_swe) + threshold_wet, threshold_partially_wet, basic_swe) end Trixi.have_nonconservative_terms(::ShallowWaterEquationsWetDry2D) = True() diff --git a/src/solvers/indicators.jl b/src/solvers/indicators.jl index 99d90bd..884d19d 100644 --- a/src/solvers/indicators.jl +++ b/src/solvers/indicators.jl @@ -40,7 +40,8 @@ end # of the shallow water equations # 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::Trixi.AbstractShallowWaterEquations, +function IndicatorHennemannGassnerShallowWater(equations::Union{Trixi.AbstractShallowWaterEquations, + AbstractShallowWaterMultiLayerEquations}, basis; alpha_max = 0.5, alpha_min = 0.001, diff --git a/src/solvers/indicators_1d.jl b/src/solvers/indicators_1d.jl index 8a9f47d..5931e0e 100644 --- a/src/solvers/indicators_1d.jl +++ b/src/solvers/indicators_1d.jl @@ -5,13 +5,14 @@ @muladd begin #! format: noindent -# Modified indicator for ShallowWaterEquationsWetDry1D to apply full FV method on elements -# containing some "dry" LGL nodes. That is, if an element is partially "wet" then it becomes a -# full FV element. +# Modified indicator for ShallowWaterEquationsWetDry1D and ShallowWaterMultiLayerEquations1D to +# apply full FV method on elements containing some "dry" LGL nodes. That is, if an element is +# partially "wet" then it becomes a full FV element. function (indicator_hg::IndicatorHennemannGassnerShallowWater)(u::AbstractArray{<:Any, 3}, mesh, - equations::ShallowWaterEquationsWetDry1D, + equations::Union{ShallowWaterEquationsWetDry1D, + ShallowWaterMultiLayerEquations1D}, dg::DGSEM, cache; kwargs...) @unpack alpha_max, alpha_min, alpha_smooth, variable = indicator_hg @@ -28,21 +29,23 @@ function (indicator_hg::IndicatorHennemannGassnerShallowWater)(u::AbstractArray{ threshold = 0.5 * 10^(-1.8 * (nnodes(dg))^0.25) parameter_s = log((1 - 0.0001) / 0.0001) - # If the water height `h` at one LGL node is lower than `threshold_partially_wet` - # the indicator sets the element-wise blending factor alpha[element] = 1 - # via the local variable `indicator_wet`. In turn, this ensures that a pure - # FV method is used in partially wet elements and guarantees the well-balanced property. - # - # Hard-coded cut-off value of `threshold_partially_wet = 1e-4` was determined through many numerical experiments. - # Overall idea is to increase robustness when computing the velocity on (nearly) dry elements which - # could be "dangerous" due to division of conservative variables, e.g., v = hv / h. - # Here, the impact of the threshold on the number of elements being updated with FV is not that - # significant. However, its impact on the robustness is very significant. - # The value can be seen as a trade-off between accuracy and stability. - # Well-balancedness of the scheme on partially wet elements with hydrostatic reconstruction - # can only be proven for the FV method (see Chen and Noelle). - # Therefore we set alpha to one regardless of its given maximum value. - threshold_partially_wet = 1e-4 + #= + If the water height `h` at one LGL node is lower than `threshold_partially_wet` + the indicator sets the element-wise blending factor alpha[element] = 1 + via the local variable `indicator_wet`. In turn, this ensures that a pure + FV method is used in partially wet elements and guarantees the well-balanced property. + + Hard-coded cut-off value of `threshold_partially_wet = 1e-4` was determined through many numerical experiments. + Overall idea is to increase robustness when computing the velocity on (nearly) dry elements which + could be "dangerous" due to division of conservative variables, e.g., v = hv / h. + Here, the impact of the threshold on the number of elements being updated with FV is not that + significant. However, its impact on the robustness is very significant. + The value can be seen as a trade-off between accuracy and stability. + Well-balancedness of the scheme on partially wet elements with hydrostatic reconstruction + can only be proven for the FV method (see Chen and Noelle). + Therefore we set alpha to one regardless of its given maximum value. + =# + threshold_partially_wet = equations.threshold_partially_wet Trixi.@threaded for element in eachelement(dg, cache) indicator = indicator_threaded[Threads.threadid()] @@ -54,9 +57,10 @@ function (indicator_hg::IndicatorHennemannGassnerShallowWater)(u::AbstractArray{ # Calculate indicator variables at Gauss-Lobatto nodes for i in eachnode(dg) u_local = get_node_vars(u, equations, dg, i, element) - h, _, _ = u_local + h = waterheight(u_local, equations) - if h <= threshold_partially_wet + # Set indicator to FV if water height is below the threshold + if minimum(h) <= threshold_partially_wet indicator_wet = 0 end diff --git a/src/solvers/indicators_2d.jl b/src/solvers/indicators_2d.jl index 1e36414..db9295e 100644 --- a/src/solvers/indicators_2d.jl +++ b/src/solvers/indicators_2d.jl @@ -42,7 +42,7 @@ function (indicator_hg::IndicatorHennemannGassnerShallowWater)(u::AbstractArray{ # Well-balancedness of the scheme on partially wet elements with hydrostatic reconstruction # can only be proven for the FV method (see Chen and Noelle). # Therefore we set alpha to be one regardless of its given value from the modal indicator. - threshold_partially_wet = 1e-4 + threshold_partially_wet = equations.threshold_partially_wet Trixi.@threaded for element in eachelement(dg, cache) indicator = indicator_threaded[Threads.threadid()] diff --git a/test/test_tree_1d.jl b/test/test_tree_1d.jl index 34a2585..d7cde02 100644 --- a/test/test_tree_1d.jl +++ b/test/test_tree_1d.jl @@ -512,21 +512,21 @@ end # 2LSWE @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_multilayer_convergence.jl"), l2=[ - 0.0007413378395106811, - 0.0005405760373610748, - 8.772497321219398e-5, - 0.0006969915241438401, - 0.0004924779641660143, - 0.00011816455915431224, + 0.0007413423106723253, + 0.000540576469967569, + 8.772494905937258e-5, + 0.0006969877681271723, + 0.0004924792705658947, + 0.00011816401639126746, 0.000139126377287091, ], linf=[ - 0.002657085141947846, - 0.0028102251075707296, - 0.00022855706065982861, - 0.003310637567035979, - 0.002490976099376596, - 0.0005211086991627756, + 0.0026570823889522366, + 0.0028102271323960926, + 0.0002285569562520129, + 0.003310623577850169, + 0.0024909798619617285, + 0.0005211057995907487, 0.00021874455861881081, ], surface_flux=(flux_lax_friedrichs, @@ -542,6 +542,43 @@ end # 2LSWE end end + @trixi_testset "elixir_shallowwater_multilayer_convergence.jl with hydrostatic_reconstruction_ersing_etal" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_multilayer_convergence.jl"), + l2=[ + 0.0007413415478415575, + 0.0005405759470722189, + 8.772503316220759e-5, + 0.0006969870846204186, + 0.0004924774643042689, + 0.00011816402423681591, + 0.000139126377287091, + ], + linf=[ + 0.0026571096352370205, + 0.0028102212809490434, + 0.00022855702706781056, + 0.0033106637858648646, + 0.0024909780761511735, + 0.0005211084155695156, + 0.00021874455861881081, + ], + surface_flux=(FluxHydrostaticReconstruction(FluxPlusDissipation(flux_ersing_etal, + DissipationLocalLaxFriedrichs()), + hydrostatic_reconstruction_ersing_etal), + FluxHydrostaticReconstruction(flux_nonconservative_ersing_etal, + hydrostatic_reconstruction_ersing_etal)), + tspan=(0.0, 0.25)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end + @trixi_testset "elixir_shallowwater_multilayer_well_balanced.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_multilayer_well_balanced.jl"), @@ -574,6 +611,63 @@ end # 2LSWE end end + @trixi_testset "elixir_shallowwater_multilayer_well_balanced_wet_dry.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_multilayer_well_balanced_wet_dry.jl"), + l2=[ + 0.025515633605684252, + 0.020321228314358498, + 2.078418885421747e-16, + 1.5034457079135132e-17, + 0.04746160729122716, + ], + linf=[ + 0.9999999999999989, + 1.0, + 1.8096275952411747e-15, + 4.245490280714591e-16, + 2.0, + ], + tspan=(0.0, 0.25)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end + + @trixi_testset "elixir_shallowwater_multilayer_well_balanced_wet_dry.jl with perturbation" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_multilayer_well_balanced_wet_dry.jl"), + l2=[ + 0.031114402256402142, + 0.02924759823866292, + 0.0021533928763101534, + 0.01412031518684538, + 0.04746160729122716, + ], + linf=[ + 1.2499999999999991, + 0.998672198687742, + 0.04297148454930247, + 0.24675612377725484, + 2.0, + ], + tspan=(0.0, 0.25), + perturbation=true) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end + @trixi_testset "elixir_shallowwater_multilayer_dam_break_ec.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_multilayer_dam_break_ec.jl"), @@ -610,21 +704,21 @@ end # 2LSWE @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_multilayer_dam_break_es.jl"), l2=[ - 0.03715839782687889, - 0.03944992348741886, - 0.047363301719025856, - 0.28686684422935993, - 0.28693864344370795, - 0.21860509590943838, + 0.037157941241555185, + 0.03944955271506525, + 0.0473661084974011, + 0.28686995762998163, + 0.2869427613009073, + 0.21860712413448047, 0.013638618139749504, ], linf=[ - 0.18088706524197873, - 0.22218392524551045, - 0.4463223829540357, - 1.1462905426260082, - 1.14002097899466, - 0.9996443057617685, + 0.18060907942891546, + 0.2221666684315402, + 0.44639752069288585, + 1.146783259827299, + 1.140481663033156, + 0.9997998219394365, 0.4999999999999993, ], tspan=(0.0, 0.25)) @@ -637,6 +731,97 @@ end # 2LSWE @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 end end + + @trixi_testset "elixir_shallowwater_multilayer_dam_break_es_dry.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_multilayer_dam_break_es_dry.jl"), + l2=[ + 0.0539713484191238, + 0.05378401084239003, + 0.05343571027541134, + 0.05294334126567629, + 0.14409266517044308, + 0.2258690603183475, + 0.22395024283599815, + 0.22041584295217526, + 0.21549254848916716, + 0.6274449804372626, + 0.013638618139749523, + ], + linf=[ + 0.27418120077237185, + 0.2738885780879326, + 0.2733530772398547, + 0.2726148157344642, + 0.8243308057901404, + 0.6976533131655047, + 0.6927619903064818, + 0.6836813560099443, + 0.6708713913679563, + 3.0233275007119254, + 0.5, + ], + tspan=(0.0, 0.25)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end + + @trixi_testset "elixir_shallowwater_multilayer_beach.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_multilayer_beach.jl"), + l2=[ + 0.5837610281036327, + 3.7410178422784273, + 6.289710784508112e-8, + ], + linf=[ + 0.9917260088977286, + 7.0962810662023035, + 4.452604027704865e-7, + ], + tspan=(0.0, 0.25), + atol=1e-5) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end + + @trixi_testset "elixir_shallowwater_multilayer_parabolic_bowl.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_multilayer_parabolic_bowl.jl"), + l2=[ + 0.002063065256405533, + 0.00048255459387538827, + 4.0617503413126664e-17, + ], + linf=[ + 0.0036250458403655466, + 0.0008665322585535845, + 1.6653345369377348e-16, + ], + tspan=(0.0, 0.25), + atol=1e-9, + coverage_override=(maxiters = 15,)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end end # MLSWE end # TreeMesh1D diff --git a/test/test_unit.jl b/test/test_unit.jl index fa729f1..64002cf 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -269,6 +269,11 @@ end rhos = (0.1, 0.2, 0.3, 0.4)) @test_throws ArgumentError initial_condition_convergence_test(0.0, 0.0, equations) end + +@timed_testset "Exception check for default_threshold functions" begin + @test_throws ArgumentError TrixiShallowWater.default_threshold_partially_wet(Int64) + @test_throws ArgumentError TrixiShallowWater.default_threshold_desingularization(Int64) +end end # Unit tests end # module diff --git a/test/test_upstream.jl b/test/test_upstream.jl index 205fd16..1a5a713 100644 --- a/test/test_upstream.jl +++ b/test/test_upstream.jl @@ -66,6 +66,7 @@ isdir(outdir) && rm(outdir, recursive = true) @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 end end + # Shallow water wet/dry 2D # TreeMesh2D @trixi_testset "TreeMesh2D: elixir_shallowwater_conical_island.jl" begin