From 9f4cfbf8775499d4becc85973222bf03be5b6dd4 Mon Sep 17 00:00:00 2001 From: svengoldberg Date: Fri, 20 Jan 2023 15:00:26 +0100 Subject: [PATCH 01/44] HR of Chen and Noelle (1D) and edit SWE struct --- src/Trixi.jl | 4 +- src/equations/numerical_fluxes.jl | 21 +++++++- src/equations/shallow_water_1d.jl | 84 +++++++++++++++++++++++++++++-- 3 files changed, 102 insertions(+), 7 deletions(-) diff --git a/src/Trixi.jl b/src/Trixi.jl index 81bdee58d5..493b9d8108 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -140,16 +140,18 @@ export LaplaceDiffusion2D, export GradientVariablesPrimitive, GradientVariablesEntropy -export flux, flux_central, flux_lax_friedrichs, flux_hll, flux_hllc, flux_hlle, flux_godunov, +export flux, flux_central, flux_lax_friedrichs, flux_hll, flux_hll_cn, flux_hllc, flux_hlle, flux_godunov, flux_chandrashekar, flux_ranocha, flux_derigs_etal, flux_hindenlang_gassner, flux_nonconservative_powell, flux_kennedy_gruber, flux_shima_etal, flux_ec, flux_fjordholm_etal, flux_nonconservative_fjordholm_etal, flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal, hydrostatic_reconstruction_audusse_etal, flux_nonconservative_audusse_etal, + hydrostatic_reconstruction_chen_noelle, flux_nonconservative_chen_noelle, FluxPlusDissipation, DissipationGlobalLaxFriedrichs, DissipationLocalLaxFriedrichs, FluxLaxFriedrichs, max_abs_speed_naive, FluxHLL, min_max_speed_naive, + FluxHLLChenNoelle, min_max_speed_chen_noelle, FluxLMARS, FluxRotated, flux_shima_etal_turbo, flux_ranocha_turbo, diff --git a/src/equations/numerical_fluxes.jl b/src/equations/numerical_fluxes.jl index ff9596848b..d90eebcc8c 100644 --- a/src/equations/numerical_fluxes.jl +++ b/src/equations/numerical_fluxes.jl @@ -203,6 +203,7 @@ struct FluxHLL{MinMaxSpeed} end FluxHLL() = FluxHLL(min_max_speed_naive) +FluxHLLChenNoelle() = FluxHLL(min_max_speed_chen_noelle) """ min_max_speed_naive(u_ll, u_rr, orientation::Integer, equations) @@ -216,6 +217,7 @@ left and right states `u_ll, u_rr`, usually based only on the local wave speeds [DOI: 10.1137/1025002](https://doi.org/10.1137/1025002) """ function min_max_speed_naive end +function min_max_speed_chen_noelle end @inline function (numflux::FluxHLL)(u_ll, u_rr, orientation_or_normal_direction, equations) λ_min, λ_max = numflux.min_max_speed(u_ll, u_rr, orientation_or_normal_direction, equations) @@ -243,7 +245,7 @@ Base.show(io::IO, numflux::FluxHLL) = print(io, "FluxHLL(", numflux.min_max_spee See [`FluxHLL`](@ref). """ const flux_hll = FluxHLL() - +const flux_hll_cn = FluxHLLChenNoelle() """ @@ -305,12 +307,27 @@ end @inline function (numflux::FluxHydrostaticReconstruction)(u_ll, u_rr, orientation_or_normal_direction, - equations::AbstractEquations) + equations::AbstractEquations{1}) @unpack numerical_flux, hydrostatic_reconstruction = numflux # Create the reconstructed left/right solution states in conservative form u_ll_star, u_rr_star = hydrostatic_reconstruction(u_ll, u_rr, equations) + threshold = equations.threshold_wet + + # Apply threshold to cut off the height and the velocity at dry interfaces + if threshold > 0 + h_ll, _, b_ll = u_ll_star + h_rr, _, b_rr = u_rr_star + + if h_ll <= threshold + u_ll_star = SVector(threshold, 0, b_ll) + end + + if h_rr <= threshold + u_rr_star = SVector(threshold, 0, b_rr) + end + end # Use the reconstructed states to compute the numerical surface flux return numerical_flux(u_ll_star, u_rr_star, orientation_or_normal_direction, equations) end diff --git a/src/equations/shallow_water_1d.jl b/src/equations/shallow_water_1d.jl index 7ccd021dd9..2ac53892cf 100644 --- a/src/equations/shallow_water_1d.jl +++ b/src/equations/shallow_water_1d.jl @@ -47,14 +47,18 @@ References for the SWE are many but a good introduction is available in Chapter struct ShallowWaterEquations1D{RealT<:Real} <: AbstractShallowWaterEquations{1, 3} gravity::RealT # gravitational constant H0::RealT # constant "lake-at-rest" total water height + threshold_limiter::RealT # threshold to use in PositivityPreservingLimiterZhangShu on waterheight, + # as shift on the initial condition and cutoff before the next timestep + threshold_wet::RealT # threshold to be applied on waterheight before calculating numflux end # Allow for flexibility to set the gravitational constant within an elixir depending on the # application where `gravity_constant=1.0` or `gravity_constant=9.81` are common values. -# The reference total water height H0 defaults to 0.0 but is used for the "lake-at-rest" -# well-balancedness test cases -function ShallowWaterEquations1D(; gravity_constant, H0=0.0) - ShallowWaterEquations1D(gravity_constant, H0) +# The reference total water height H0 is an artefact from the old calculation of the lake_at_rest_error +# Strict default values for thresholds that performed great in several numerical experiments +function ShallowWaterEquations1D(; gravity_constant, H0=0.0, + threshold_limiter=1e-13, threshold_wet=1e-15) + ShallowWaterEquations1D(gravity_constant, H0, threshold_limiter, threshold_wet) end @@ -306,6 +310,35 @@ Further details on the hydrostatic reconstruction and its motivation can be foun z) end +@inline function flux_nonconservative_chen_noelle(u_ll, u_rr, + orientation::Integer, + equations::ShallowWaterEquations1D) + + # Pull the water height and bottom topography on the left + h_ll, _, b_ll = u_ll + h_rr, _, b_rr = u_rr + + H_ll = h_ll+b_ll + H_rr = h_rr+b_rr + + b_star = min( max( b_ll, b_rr ), min( H_ll, H_rr ) ) + + # Create the hydrostatic reconstruction for the left solution state + u_ll_star, _ = hydrostatic_reconstruction_chen_noelle(u_ll, u_rr, equations) + + # Copy the reconstructed water height for easier to read code + h_ll_star = u_ll_star[1] + + z = zero(eltype(u_ll)) + # Includes two parts: + # (i) Diagonal (consistent) term from the volume flux that uses `b_ll` to avoid + # cross-averaging across a discontinuous bottom topography + # (ii) True surface part that uses `h_ll` and `h_ll_star` to handle discontinuous bathymetry + return SVector(z, + equations.gravity * h_ll * b_ll - equations.gravity * (h_ll_star+h_ll) * (b_ll-b_star), + z) + +end """ flux_fjordholm_etal(u_ll, u_rr, orientation, @@ -408,6 +441,30 @@ Further details on this hydrostatic reconstruction and its motivation can be fou return u_ll_star, u_rr_star end +@inline function hydrostatic_reconstruction_chen_noelle(u_ll, u_rr, equations::ShallowWaterEquations1D) + # Unpack left and right water heights and bottom topographies + h_ll, _, b_ll = u_ll + h_rr, _, b_rr = u_rr + + # Get the velocities on either side + v_ll = velocity(u_ll, equations) + v_rr = velocity(u_rr, equations) + + H_ll = b_ll+h_ll + H_rr = b_rr+h_rr + + b_star = min( max( b_ll, b_rr ), min( H_ll, H_rr ) ) + + # Compute the reconstructed water heights + h_ll_star = min( H_ll-b_star, h_ll ) + h_rr_star = min( H_rr-b_star, h_rr ) + + # Create the conservative variables using the reconstruted water heights + u_ll_star = SVector( h_ll_star, h_ll_star * v_ll, b_ll ) + u_rr_star = SVector( h_rr_star, h_rr_star * v_rr, b_rr ) + + return u_ll_star, u_rr_star +end # Calculate maximum wave speed for local Lax-Friedrichs-type dissipation as the # maximum velocity magnitude plus the maximum speed of sound @@ -472,6 +529,25 @@ end end +@inline function min_max_speed_chen_noelle(u_ll, u_rr, orientation::Integer, equations::ShallowWaterEquations1D) + # Get the velocity quantities + v_ll = velocity(u_ll, equations) + v_rr = velocity(u_rr, equations) + + # Calculate the wave celerity on the left and right + h_ll = waterheight(u_ll, equations) + h_rr = waterheight(u_rr, equations) + + a_ll = sqrt(equations.gravity * h_ll) + a_rr = sqrt(equations.gravity * h_rr) + + λ_min = min( v_ll-a_ll, v_rr-a_rr, 0 ) + λ_max = max( v_ll+a_ll, v_rr+a_rr, 0 ) + + return λ_min, λ_max +end + + @inline function max_abs_speeds(u, equations::ShallowWaterEquations1D) h = waterheight(u, equations) v = velocity(u, equations) From 529fe9e8772b723fe3c0e5588a398930b03b91fb Mon Sep 17 00:00:00 2001 From: svengoldberg Date: Fri, 20 Jan 2023 15:02:39 +0100 Subject: [PATCH 02/44] Overload limiter (SWE 1D) to cut off waterheight --- .../positivity_zhang_shu_dg1d.jl | 58 +++++++++++++++++++ src/solvers/dg.jl | 7 +++ 2 files changed, 65 insertions(+) diff --git a/src/callbacks_stage/positivity_zhang_shu_dg1d.jl b/src/callbacks_stage/positivity_zhang_shu_dg1d.jl index 50d6b3f2c3..e116ad8139 100644 --- a/src/callbacks_stage/positivity_zhang_shu_dg1d.jl +++ b/src/callbacks_stage/positivity_zhang_shu_dg1d.jl @@ -43,5 +43,63 @@ function limiter_zhang_shu!(u, threshold::Real, variable, return nothing end +# Overload limiter for ShallowWaterEquations1D and cut off waterheight as well as velocity +# to avoid numerical problems/instabilities caused by velocities close to zero +function limiter_zhang_shu!(u, threshold::Real, variable, + mesh::AbstractMesh{1}, equations::ShallowWaterEquations1D, + dg::DGSEM, cache) + @unpack weights = dg.basis + + @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 || 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, as we assume that + # Jensen's inequality holds (e.g. pressure for compressible Euler equations). + 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) + + # Cut velocity off in case of waterheight beeing smaller than threshold + + h_node, h_v_node, b_node = u_node + h_mean, h_v_mean, _ = u_mean #b_mean is not used as b_node must not be overwritten + + # Set them both to zero to apply linear combination correctly + h_v_node = h_v_node * Int32(h_node > threshold) + + h_v_mean = h_v_mean * Int32(h_node > threshold) + + u_node = SVector(h_node, h_v_node, b_node) + u_mean = SVector(h_mean, h_v_mean, 0) + + # When velocity is cut off, the only value beeing averaged is the waterheight, + # because the velocity is set to zero and this value is passed. + # Otherwise, the velocity is averaged, as well. + set_node_vars!(u, theta * u_node + (1-theta) * u_mean, + (1,2), dg, i, element) + end + end + + return nothing +end + end # @muladd diff --git a/src/solvers/dg.jl b/src/solvers/dg.jl index f18fbca178..bf06f486dd 100644 --- a/src/solvers/dg.jl +++ b/src/solvers/dg.jl @@ -490,6 +490,13 @@ end return nothing end +@inline function set_node_vars!(u, u_node, variables::NTuple{N, Integer}, solver::DG, indices...) where N + for v in variables + u[v, indices...] = u_node[v] + end + return nothing +end + @inline function add_to_node_vars!(u, u_node, equations, solver::DG, indices...) for v in eachvariable(equations) u[v, indices...] += u_node[v] From e87dbf291f52d75397a4c0992a46c2039ff16752 Mon Sep 17 00:00:00 2001 From: svengoldberg Date: Fri, 20 Jan 2023 15:04:19 +0100 Subject: [PATCH 03/44] New indicatorHG (SWE 1D) to apply FV on dry cells --- src/solvers/dgsem_tree/indicators_1d.jl | 91 +++++++++++++++++++++++++ 1 file changed, 91 insertions(+) diff --git a/src/solvers/dgsem_tree/indicators_1d.jl b/src/solvers/dgsem_tree/indicators_1d.jl index f501c44fcb..555921443b 100644 --- a/src/solvers/dgsem_tree/indicators_1d.jl +++ b/src/solvers/dgsem_tree/indicators_1d.jl @@ -23,6 +23,97 @@ function create_cache(typ::Type{IndicatorHennemannGassner}, mesh, equations::Abs create_cache(typ, equations, dg.basis) end +# Overload indicator when ShallowWaterEquations1D is used to apply full FV method on cells +# containing dry LGL nodes +function (indicator_hg::IndicatorHennemannGassner)(u, mesh::Union{TreeMesh{1}, StructuredMesh{1}}, + equations::ShallowWaterEquations1D, dg::DGSEM, + cache; kwargs...) + @unpack alpha_max, alpha_min, alpha_smooth, variable = indicator_hg + @unpack alpha, alpha_tmp, indicator_threaded, modal_threaded = indicator_hg.cache + # TODO: Taal refactor, when to `resize!` stuff changed possibly by AMR? + # Shall we implement `resize!(semi::AbstractSemidiscretization, new_size)` + # or just `resize!` whenever we call the relevant methods as we do now? + resize!(alpha, nelements(dg, cache)) + if alpha_smooth + resize!(alpha_tmp, nelements(dg, cache)) + end + + # magic parameters + threshold = 0.5 * 10^(-1.8 * (nnodes(dg))^0.25) + parameter_s = log((1 - 0.0001)/0.0001) + + # If h at one LGL node is lower than threshold_wet, set alpha[element]=1 via indicator_wet + # to apply a full FV method and guarantee well balance property. + #= + Determination of hard coded treshold: Showed good results in numerical experiments. Idea is to + gain more stability of velocity (and therefore the whole PDE) on (nearly) dry cells which + could be counteracted by division of conservative variables (hv/v). + Here, the impact of the threshold on the number of cells beeing updated with FV is not that + huge. Rather huge is the impact on tha stability. + The value can be seen as a trade-off between accuracy and stability. + Well balancedness of scheme using hydrostatic reconstruction can only be proved for FV method. + Therefore we set alpha to one regardless its given maximum value. + =# + threshold_wet = 1e-4 + + @threaded for element in eachelement(dg, cache) + indicator = indicator_threaded[Threads.threadid()] + modal = modal_threaded[Threads.threadid()] + + # (Re-)set dummy variable for alpha_dry + indicator_wet = 1.0 + + # 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 + indicator_wet = min(indicator_wet, Int32(h > threshold_wet)) + indicator[i] = indicator_hg.variable(u_local, equations) + end + + # Convert to modal representation + multiply_scalar_dimensionwise!(modal, dg.basis.inverse_vandermonde_legendre, indicator) + + # Calculate total energies for all modes, without highest, without two highest + total_energy = zero(eltype(modal)) + for i in 1:nnodes(dg) + total_energy += modal[i]^2 + end + total_energy_clip1 = zero(eltype(modal)) + for i in 1:(nnodes(dg)-1) + total_energy_clip1 += modal[i]^2 + end + total_energy_clip2 = zero(eltype(modal)) + for i in 1:(nnodes(dg)-2) + total_energy_clip2 += modal[i]^2 + end + + # Calculate energy in higher modes + energy = max((total_energy - total_energy_clip1) / total_energy, + (total_energy_clip1 - total_energy_clip2) / total_energy_clip1) + + alpha_element = 1 / (1 + exp(-parameter_s / threshold * (energy - threshold))) + + # Take care of the case close to pure DG + if alpha_element < alpha_min + alpha_element = zero(alpha_element) + end + + # Take care of the case close to pure FV + if alpha_element > 1 - alpha_min + alpha_element = one(alpha_element) + end + + # Clip the maximum amount of FV allowed or set to one depending on indicator_wet + alpha[element] = indicator_wet * min(alpha_max, alpha_element) - (indicator_wet-1) + end + + if alpha_smooth + apply_smoothing!(mesh, alpha, alpha_tmp, dg, cache) + end + + return alpha +end function (indicator_hg::IndicatorHennemannGassner)(u, mesh::Union{TreeMesh{1}, StructuredMesh{1}}, equations, dg::DGSEM, cache; From fac597232bcc0d7b7e9c087b94a302ba8fcbb332 Mon Sep 17 00:00:00 2001 From: svengoldberg Date: Fri, 20 Jan 2023 15:07:00 +0100 Subject: [PATCH 04/44] Threshold in rhs! before calculation (SWE 1D) --- src/solvers/dgsem_tree/dg_1d.jl | 35 +++++++++++++++++++++++++++++++++ 1 file changed, 35 insertions(+) diff --git a/src/solvers/dgsem_tree/dg_1d.jl b/src/solvers/dgsem_tree/dg_1d.jl index 6a646b38c5..f731204001 100644 --- a/src/solvers/dgsem_tree/dg_1d.jl +++ b/src/solvers/dgsem_tree/dg_1d.jl @@ -76,6 +76,9 @@ function rhs!(du, u, t, # Reset du @trixi_timeit timer() "reset ∂u/∂t" reset_du!(du, dg, cache) + # Apply waterheight and velocity cut off for ShallowWaterEquations1D + apply_thresholds!(u, equations, dg, cache) + # Calculate volume integral @trixi_timeit timer() "volume integral" calc_volume_integral!( du, u, mesh, @@ -116,6 +119,38 @@ function rhs!(du, u, t, return nothing end +# Special case for SWE 1D: Apply the threshold to the waterheight and cut off velocity in dry cells +function apply_thresholds!(u, equations::ShallowWaterEquations1D, dg::DGSEM, cache) + + threshold = equations.threshold_limiter + + # If no threshold is set to zero, there is nothing to do + if threshold == 0 + return nothing + end + + @threaded for element in eachelement(dg, cache) + + for i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, element) + + h, v, b = u_node + + h = h * Int32(h > threshold) + threshold * Int32(h <= threshold) + v = v * Int32(h > threshold) + + u_node = SVector(h, v, b) + + set_node_vars!(u, u_node, equations, dg, i, element) + end + end + + return nothing +end + +function apply_thresholds!(u, equations::AbstractEquations{1}, dg::DGSEM, cache) + return nothing +end function calc_volume_integral!(du, u, mesh::Union{TreeMesh{1}, StructuredMesh{1}}, From 2443b8d2684d1918994b3da0d61c0daa18be4e12 Mon Sep 17 00:00:00 2001 From: svengoldberg Date: Fri, 20 Jan 2023 15:08:50 +0100 Subject: [PATCH 05/44] New lake_at_rest_error for SWE 1D --- src/callbacks_step/analysis.jl | 11 +++++++++++ src/callbacks_step/analysis_dg1d.jl | 14 ++++++++++++++ src/equations/shallow_water_1d.jl | 9 +++++++-- 3 files changed, 32 insertions(+), 2 deletions(-) diff --git a/src/callbacks_step/analysis.jl b/src/callbacks_step/analysis.jl index db27bc384e..cbdd3a8097 100644 --- a/src/callbacks_step/analysis.jl +++ b/src/callbacks_step/analysis.jl @@ -516,6 +516,17 @@ end pretty_form_utf(quantity) = get_name(quantity) pretty_form_ascii(quantity) = get_name(quantity) +# Overload analyze and pass the initial_condition as parameter for integrate +# Special case only used for lake_at_rest_error for AbstractShallowWaterEquations +# With the new wet-dry mechanics, the old calculation using a constant H0 as a reference value fails +# Imagine e.g. a dry hill beeing higher than the waterlevel +# New approach: Compare waterheight h at time T with the initial waterheight +function analyze(::typeof(lake_at_rest_error), du, u, t, semi::AbstractSemidiscretization) + mesh, equations, solver, cache = mesh_equations_solver_cache(semi) + @unpack initial_condition = semi + + integrate(lake_at_rest_error, u, mesh, equations, solver, cache, initial_condition, normalize=true) +end function entropy_timederivative end pretty_form_utf(::typeof(entropy_timederivative)) = "∑∂S/∂U ⋅ Uₜ" diff --git a/src/callbacks_step/analysis_dg1d.jl b/src/callbacks_step/analysis_dg1d.jl index e92701dc1f..3e3f08a107 100644 --- a/src/callbacks_step/analysis_dg1d.jl +++ b/src/callbacks_step/analysis_dg1d.jl @@ -175,6 +175,20 @@ function integrate(func::Func, u, end end +# Overload to pass the initial_condition as parameter for calculating lake_at_rest_error +function integrate(::typeof(lake_at_rest_error), u, + mesh::Union{TreeMesh{1}, StructuredMesh{1}}, + equations::ShallowWaterEquations1D, dg::DGSEM, cache, + initial_condition; normalize=true) + node_coordinates = cache.elements.node_coordinates + integrate_via_indices(u, mesh, equations, dg, cache; normalize=normalize) do u, i, element, equations, dg + x = get_node_coords(node_coordinates, equations, dg, i, element) + + u_exact = initial_condition(x, 0.0, equations) + u_local = get_node_vars(u, equations, dg, i, element) + return lake_at_rest_error(u_local, u_exact, equations) + end +end function analyze(::typeof(entropy_timederivative), du, u, t, mesh::Union{TreeMesh{1},StructuredMesh{1}}, equations, dg::DG, cache) diff --git a/src/equations/shallow_water_1d.jl b/src/equations/shallow_water_1d.jl index 2ac53892cf..644b203c21 100644 --- a/src/equations/shallow_water_1d.jl +++ b/src/equations/shallow_water_1d.jl @@ -658,9 +658,14 @@ end # Calculate the error for the "lake-at-rest" test case where H = h+b should # be a constant value over time -@inline function lake_at_rest_error(u, equations::ShallowWaterEquations1D) +@inline function lake_at_rest_error(u, u_exact, equations::ShallowWaterEquations1D) h, _, b = u - return abs(equations.H0 - (h + b)) + h_exact, _, b_exact= u_exact + + H = h + b + H_exact = h_exact + b_exact + + return abs(H - H_exact) end end # @muladd From 0decdaf1a7dc78b2817d9343515d4af7f08299ba Mon Sep 17 00:00:00 2001 From: svengoldberg Date: Fri, 20 Jan 2023 15:35:18 +0100 Subject: [PATCH 06/44] New wet/dry elixirs for testing scheme for SWE 1D --- .../elixir_shallowwater_beach.jl | 107 +++++++++++++++++ .../elixir_shallowwater_parabolic_bowl.jl | 109 ++++++++++++++++++ ...ixir_shallowwater_well_balanced_wet_dry.jl | 100 ++++++++++++++++ 3 files changed, 316 insertions(+) create mode 100644 examples/tree_1d_dgsem/elixir_shallowwater_beach.jl create mode 100644 examples/tree_1d_dgsem/elixir_shallowwater_parabolic_bowl.jl create mode 100644 examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_beach.jl b/examples/tree_1d_dgsem/elixir_shallowwater_beach.jl new file mode 100644 index 0000000000..273e8f26e6 --- /dev/null +++ b/examples/tree_1d_dgsem/elixir_shallowwater_beach.jl @@ -0,0 +1,107 @@ + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# Semidiscretization of the shallow water equations + +equations = ShallowWaterEquations1D(gravity_constant=9.812) +cfl = 0.2 + +function initial_condition_beach(x, t, equations:: ShallowWaterEquations1D) + D = 1 + δ = 0.02 + γ = sqrt((3*δ)/(4*D)) + xₐ = sqrt((4*D)/(3*δ)) * acosh(sqrt(20)) + + f = D + 40*δ * sech(γ*(8*x[1]-xₐ))^2 + + # steep wall + b = 0.01 + 99/409600 * 4^x[1] + + if x[1] >= 6 + H = b + v = 0 + else + H = f + v = sqrt(equations.gravity/D) * H + end + + 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_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) +surface_flux = (FluxHydrostaticReconstruction(flux_hll_cn, hydrostatic_reconstruction_chen_noelle), + flux_nonconservative_chen_noelle) + +basis = LobattoLegendreBasis(7) + +indicator_sc = IndicatorHennemannGassner(equations, basis, + alpha_max=.5, + alpha_min=.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=8, + 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, + lake_at_rest_error)) + +alive_callback = AliveCallback(analysis_interval=analysis_interval) + +save_solution = SaveSolutionCallback(interval=1000, + save_initial_solution=true, + save_final_solution=true) + +stepsize_callback = StepsizeCallback(cfl=cfl) + +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, + stepsize_callback) + +stage_limiter! = PositivityPreservingLimiterZhangShu(thresholds=(equations.threshold_limiter,), + variables=(Trixi.waterheight,)) + +############################################################################### +# run the simulation + +sol = solve(ode, SSPRK43(stage_limiter!), dt=1.0, + save_everystep=false, callback=callbacks, adaptive=false, maxiters=1e+9); + +summary_callback() # print the timer summary \ No newline at end of file diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_parabolic_bowl.jl b/examples/tree_1d_dgsem/elixir_shallowwater_parabolic_bowl.jl new file mode 100644 index 0000000000..724e5d96a6 --- /dev/null +++ b/examples/tree_1d_dgsem/elixir_shallowwater_parabolic_bowl.jl @@ -0,0 +1,109 @@ + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# Semidiscretization of the shallow water equations + +equations = ShallowWaterEquations1D(gravity_constant=9.81) +cfl = 0.3 + + +function parabolic_bowl_analytic_1D_H(gravity, x,t) + a = 1 + h_0 = 0.1 + σ = 0.5 + ω = sqrt(2*gravity*h_0)/a + + H = σ * h_0/a^2 * (2*x[1]*cos(ω*t) - σ) + h_0 + return H +end + +# Implemented based on Wintermeyer (here seen in 2D) +function initial_condition_parabolic_bowl(x, t, equations:: ShallowWaterEquations1D) + a = 1 + h_0 = 0.1 + σ = 0.5 + ω = sqrt(2*equations.gravity*h_0)/a + + v = -σ*ω*sin(ω*t) + + b = h_0 * x[1]^2/a^2 + + H = max(b, parabolic_bowl_analytic_1D_H(equations.gravity, x, t)) + + 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_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) +surface_flux = (FluxHydrostaticReconstruction(flux_hll_cn, hydrostatic_reconstruction_chen_noelle), + flux_nonconservative_chen_noelle) + +basis = LobattoLegendreBasis(5) + +indicator_sc = IndicatorHennemannGassner(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=10, + 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, + lake_at_rest_error)) + +alive_callback = AliveCallback(analysis_interval=analysis_interval) + +save_solution = SaveSolutionCallback(interval=1000, + save_initial_solution=true, + save_final_solution=true) + +stepsize_callback = StepsizeCallback(cfl=cfl) + +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, + stepsize_callback) + +stage_limiter! = PositivityPreservingLimiterZhangShu(thresholds=(equations.threshold_limiter,), + variables=(Trixi.waterheight,)) + +############################################################################### +# run the simulation + +sol = solve(ode, SSPRK43(stage_limiter!), dt=1.0, + save_everystep=false, callback=callbacks, adaptive=false, maxiters=1e+9); + +summary_callback() # print the timer summary \ No newline at end of file diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl b/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl new file mode 100644 index 0000000000..b63d13433b --- /dev/null +++ b/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl @@ -0,0 +1,100 @@ + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# Semidiscretization of the shallow water equations + +equations = ShallowWaterEquations1D(gravity_constant=9.812) +cfl = .3 + +function initial_condition_well_balanced_CN(x, t, equations:: ShallowWaterEquations1D) + v = 0 + b = sin(4*pi*x[1]) + 3 + + if x[1] >= 0.5 + b = sin(4*pi*x[1]) + 1 + end + + H = max(b, 2.5) + + if x[1] >= 0.5 + H = max(b, 1.5) + end + + H = max(H, b + equations.threshold_limiter) + return prim2cons(SVector(H, v, b), equations) +end + +initial_condition = initial_condition_well_balanced_CN + +############################################################################### +# Get the DG approximation space + +volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) +surface_flux = (FluxHydrostaticReconstruction(flux_hll_cn, hydrostatic_reconstruction_chen_noelle), + flux_nonconservative_chen_noelle) + +basis = LobattoLegendreBasis(5) + +indicator_sc = IndicatorHennemannGassner(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, 1] + +coordinates_min = 0. +coordinates_max = 1. + +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level=4, + n_cells_max=10_000, + ) + +# create the semi discretization object +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 25.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, + lake_at_rest_error)) + +alive_callback = AliveCallback(analysis_interval=analysis_interval) + +save_solution = SaveSolutionCallback(interval=1000, + save_initial_solution=true, + save_final_solution=true) + +stepsize_callback = StepsizeCallback(cfl=cfl) + +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, + stepsize_callback) + +stage_limiter! = PositivityPreservingLimiterZhangShu(thresholds=(equations.threshold_limiter,), + 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 \ No newline at end of file From 62c7265a8138336b4417a98f09c3dadca4ebec43 Mon Sep 17 00:00:00 2001 From: svengoldberg Date: Sat, 21 Jan 2023 18:09:35 +0100 Subject: [PATCH 07/44] Fixed MethodError; apply_thresholds! too strict --- src/solvers/dgsem_tree/dg_1d.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/solvers/dgsem_tree/dg_1d.jl b/src/solvers/dgsem_tree/dg_1d.jl index f731204001..bedc3c7ee1 100644 --- a/src/solvers/dgsem_tree/dg_1d.jl +++ b/src/solvers/dgsem_tree/dg_1d.jl @@ -148,7 +148,7 @@ function apply_thresholds!(u, equations::ShallowWaterEquations1D, dg::DGSEM, cac return nothing end -function apply_thresholds!(u, equations::AbstractEquations{1}, dg::DGSEM, cache) +function apply_thresholds!(u, equations::AbstractEquations{1}, dg, cache) return nothing end From 36517d234aea7e317f65cf3332debf70baf76f09 Mon Sep 17 00:00:00 2001 From: svengoldberg Date: Thu, 9 Feb 2023 14:36:30 +0100 Subject: [PATCH 08/44] Move threshold on volume integral in stage_limiter --- .../positivity_zhang_shu_dg1d.jl | 31 +++++++++++++--- src/solvers/dgsem_tree/dg_1d.jl | 35 ------------------- 2 files changed, 26 insertions(+), 40 deletions(-) diff --git a/src/callbacks_stage/positivity_zhang_shu_dg1d.jl b/src/callbacks_stage/positivity_zhang_shu_dg1d.jl index e116ad8139..094e04b03f 100644 --- a/src/callbacks_stage/positivity_zhang_shu_dg1d.jl +++ b/src/callbacks_stage/positivity_zhang_shu_dg1d.jl @@ -44,7 +44,7 @@ function limiter_zhang_shu!(u, threshold::Real, variable, end # Overload limiter for ShallowWaterEquations1D and cut off waterheight as well as velocity -# to avoid numerical problems/instabilities caused by velocities close to zero +# to avoid numerical problems/instabilities caused by velocities close to zero function limiter_zhang_shu!(u, threshold::Real, variable, mesh::AbstractMesh{1}, equations::ShallowWaterEquations1D, dg::DGSEM, cache) @@ -77,10 +77,10 @@ function limiter_zhang_shu!(u, threshold::Real, variable, for i in eachnode(dg) u_node = get_node_vars(u, equations, dg, i, element) - # Cut velocity off in case of waterheight beeing smaller than threshold + # Cut off velocity in case that the waterheight is smaller than the threshold h_node, h_v_node, b_node = u_node - h_mean, h_v_mean, _ = u_mean #b_mean is not used as b_node must not be overwritten + h_mean, h_v_mean, _ = u_mean # b_mean is not used as b_node must not be overwritten # Set them both to zero to apply linear combination correctly h_v_node = h_v_node * Int32(h_node > threshold) @@ -88,16 +88,37 @@ function limiter_zhang_shu!(u, threshold::Real, variable, h_v_mean = h_v_mean * Int32(h_node > threshold) u_node = SVector(h_node, h_v_node, b_node) - u_mean = SVector(h_mean, h_v_mean, 0) + u_mean = SVector(h_mean, h_v_mean, zero(eltype(u))) # When velocity is cut off, the only value beeing averaged is the waterheight, - # because the velocity is set to zero and this value is passed. + # because the velocity is set to zero and this value is passed. # Otherwise, the velocity is averaged, as well. + # Note that the auxiliary bottom topography variable `b` is never limited. set_node_vars!(u, theta * u_node + (1-theta) * u_mean, (1,2), dg, i, element) end end + # An extra "safety" check is done over all the degrees of + # freedom after the limiting 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 limiter. + @threaded for element in eachelement(dg, cache) + + for i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, element) + + h, v, b = u_node + + h = h * Int32(h > threshold) + threshold * Int32(h <= threshold) + v = v * Int32(h > threshold) + + u_node = SVector(h, v, b) + + set_node_vars!(u, u_node, equations, dg, i, element) + end + end + return nothing end diff --git a/src/solvers/dgsem_tree/dg_1d.jl b/src/solvers/dgsem_tree/dg_1d.jl index bedc3c7ee1..6a646b38c5 100644 --- a/src/solvers/dgsem_tree/dg_1d.jl +++ b/src/solvers/dgsem_tree/dg_1d.jl @@ -76,9 +76,6 @@ function rhs!(du, u, t, # Reset du @trixi_timeit timer() "reset ∂u/∂t" reset_du!(du, dg, cache) - # Apply waterheight and velocity cut off for ShallowWaterEquations1D - apply_thresholds!(u, equations, dg, cache) - # Calculate volume integral @trixi_timeit timer() "volume integral" calc_volume_integral!( du, u, mesh, @@ -119,38 +116,6 @@ function rhs!(du, u, t, return nothing end -# Special case for SWE 1D: Apply the threshold to the waterheight and cut off velocity in dry cells -function apply_thresholds!(u, equations::ShallowWaterEquations1D, dg::DGSEM, cache) - - threshold = equations.threshold_limiter - - # If no threshold is set to zero, there is nothing to do - if threshold == 0 - return nothing - end - - @threaded for element in eachelement(dg, cache) - - for i in eachnode(dg) - u_node = get_node_vars(u, equations, dg, i, element) - - h, v, b = u_node - - h = h * Int32(h > threshold) + threshold * Int32(h <= threshold) - v = v * Int32(h > threshold) - - u_node = SVector(h, v, b) - - set_node_vars!(u, u_node, equations, dg, i, element) - end - end - - return nothing -end - -function apply_thresholds!(u, equations::AbstractEquations{1}, dg, cache) - return nothing -end function calc_volume_integral!(du, u, mesh::Union{TreeMesh{1}, StructuredMesh{1}}, From 71501c382fa52f5c2ca9aa293b44bc7a08a00eb9 Mon Sep 17 00:00:00 2001 From: svengoldberg Date: Thu, 9 Feb 2023 14:41:30 +0100 Subject: [PATCH 09/44] Indentation, spacing and comments adjustment --- src/callbacks_step/analysis.jl | 8 ++++---- src/equations/numerical_fluxes.jl | 9 +++++++-- src/equations/shallow_water_1d.jl | 20 +++++++++++--------- src/solvers/dgsem_tree/indicators_1d.jl | 23 ++++++++++++----------- 4 files changed, 34 insertions(+), 26 deletions(-) diff --git a/src/callbacks_step/analysis.jl b/src/callbacks_step/analysis.jl index cbdd3a8097..d0a7ee7393 100644 --- a/src/callbacks_step/analysis.jl +++ b/src/callbacks_step/analysis.jl @@ -516,10 +516,10 @@ end pretty_form_utf(quantity) = get_name(quantity) pretty_form_ascii(quantity) = get_name(quantity) -# Overload analyze and pass the initial_condition as parameter for integrate -# Special case only used for lake_at_rest_error for AbstractShallowWaterEquations -# With the new wet-dry mechanics, the old calculation using a constant H0 as a reference value fails -# Imagine e.g. a dry hill beeing higher than the waterlevel +# Overload analyze and pass the initial_condition as parameter for integrate. +# Special case only used for lake_at_rest_error for AbstractShallowWaterEquations. +# With the new wet-dry mechanics, the old calculation using a constant H0 as a reference value fails. +# For example, when a hill extends to be higher than the waterlevel. # New approach: Compare waterheight h at time T with the initial waterheight function analyze(::typeof(lake_at_rest_error), du, u, t, semi::AbstractSemidiscretization) mesh, equations, solver, cache = mesh_equations_solver_cache(semi) diff --git a/src/equations/numerical_fluxes.jl b/src/equations/numerical_fluxes.jl index d90eebcc8c..7667f0a31c 100644 --- a/src/equations/numerical_fluxes.jl +++ b/src/equations/numerical_fluxes.jl @@ -313,6 +313,11 @@ end # Create the reconstructed left/right solution states in conservative form u_ll_star, u_rr_star = hydrostatic_reconstruction(u_ll, u_rr, equations) + # Set the water height to be at least the value stored in the variable threshold after + # the hydrostatic reconstruction is applied and before the numerical flux is calculated + # to avoid numerical problem with arbitrary small values. Interfaces with a water height + # lower or equal to the threshold can be declared as dry. + # The default value is set to 1e-15 and can be changed within the constructor call in an elixir. threshold = equations.threshold_wet # Apply threshold to cut off the height and the velocity at dry interfaces @@ -321,11 +326,11 @@ end h_rr, _, b_rr = u_rr_star if h_ll <= threshold - u_ll_star = SVector(threshold, 0, b_ll) + u_ll_star = SVector(threshold, zero(eltype(u_ll)), b_ll) end if h_rr <= threshold - u_rr_star = SVector(threshold, 0, b_rr) + u_rr_star = SVector(threshold, zero(eltype(u_rr)), b_rr) end end # Use the reconstructed states to compute the numerical surface flux diff --git a/src/equations/shallow_water_1d.jl b/src/equations/shallow_water_1d.jl index 644b203c21..9d010ef6c2 100644 --- a/src/equations/shallow_water_1d.jl +++ b/src/equations/shallow_water_1d.jl @@ -47,9 +47,11 @@ References for the SWE are many but a good introduction is available in Chapter struct ShallowWaterEquations1D{RealT<:Real} <: AbstractShallowWaterEquations{1, 3} gravity::RealT # gravitational constant H0::RealT # constant "lake-at-rest" total water height - threshold_limiter::RealT # threshold to use in PositivityPreservingLimiterZhangShu on waterheight, - # as shift on the initial condition and cutoff before the next timestep - threshold_wet::RealT # threshold to be applied on waterheight before calculating numflux + threshold_limiter::RealT # Threshold to use in PositivityPreservingLimiterZhangShu on waterheight, + # as a (small) shift on the initial condition and cutoff before the + # next time step. + threshold_wet::RealT # Threshold to be applied on water height to define when the flow is "wet" + # before calculating the numerical flux. end # Allow for flexibility to set the gravitational constant within an elixir depending on the @@ -318,8 +320,8 @@ end h_ll, _, b_ll = u_ll h_rr, _, b_rr = u_rr - H_ll = h_ll+b_ll - H_rr = h_rr+b_rr + H_ll = h_ll + b_ll + H_rr = h_rr + b_rr b_star = min( max( b_ll, b_rr ), min( H_ll, H_rr ) ) @@ -456,8 +458,8 @@ end b_star = min( max( b_ll, b_rr ), min( H_ll, H_rr ) ) # Compute the reconstructed water heights - h_ll_star = min( H_ll-b_star, h_ll ) - h_rr_star = min( H_rr-b_star, h_rr ) + h_ll_star = min( H_ll - b_star, h_ll ) + h_rr_star = min( H_rr - b_star, h_rr ) # Create the conservative variables using the reconstruted water heights u_ll_star = SVector( h_ll_star, h_ll_star * v_ll, b_ll ) @@ -541,8 +543,8 @@ end a_ll = sqrt(equations.gravity * h_ll) a_rr = sqrt(equations.gravity * h_rr) - λ_min = min( v_ll-a_ll, v_rr-a_rr, 0 ) - λ_max = max( v_ll+a_ll, v_rr+a_rr, 0 ) + λ_min = min( v_ll - a_ll, v_rr - a_rr, zero(eltype(u_ll)) ) + λ_max = max( v_ll + a_ll, v_rr + a_rr, zero(eltype(u_ll)) ) return λ_min, λ_max end diff --git a/src/solvers/dgsem_tree/indicators_1d.jl b/src/solvers/dgsem_tree/indicators_1d.jl index 555921443b..a5c2eef18f 100644 --- a/src/solvers/dgsem_tree/indicators_1d.jl +++ b/src/solvers/dgsem_tree/indicators_1d.jl @@ -43,17 +43,18 @@ function (indicator_hg::IndicatorHennemannGassner)(u, mesh::Union{TreeMesh{1}, S parameter_s = log((1 - 0.0001)/0.0001) # If h at one LGL node is lower than threshold_wet, set alpha[element]=1 via indicator_wet - # to apply a full FV method and guarantee well balance property. - #= - Determination of hard coded treshold: Showed good results in numerical experiments. Idea is to - gain more stability of velocity (and therefore the whole PDE) on (nearly) dry cells which - could be counteracted by division of conservative variables (hv/v). - Here, the impact of the threshold on the number of cells beeing updated with FV is not that - huge. Rather huge is the impact on tha stability. - The value can be seen as a trade-off between accuracy and stability. - Well balancedness of scheme using hydrostatic reconstruction can only be proved for FV method. - Therefore we set alpha to one regardless its given maximum value. - =# + # in order to apply a full FV method in the partially wet cell and guarantee the well-balanced property. + # + # Determination of hard coded threshold: Showed good results in many numerical experiments. + # The idea is to gain more stability when computing the velocity on (nearly) dry cells which + # could be counteracted by division of conservative variables, e.g., v = hv / h. + # Here, the impact of the threshold on the number of cells being updated with FV is not that + # significant. However, its impact on the stability is very significant. + # The value can be seen as a trade-off between accuracy and stability. + # Well-balancedness of the scheme on partially wet cells 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_wet = 1e-4 @threaded for element in eachelement(dg, cache) From bbb3f8c02360f799c2b687a08ebcd1901b1f756f Mon Sep 17 00:00:00 2001 From: svengoldberg Date: Thu, 9 Feb 2023 15:09:24 +0100 Subject: [PATCH 10/44] Renaming numerical HLL type flux (SWE 1D) --- examples/tree_1d_dgsem/elixir_shallowwater_beach.jl | 2 +- examples/tree_1d_dgsem/elixir_shallowwater_parabolic_bowl.jl | 2 +- .../elixir_shallowwater_well_balanced_wet_dry.jl | 2 +- src/Trixi.jl | 4 ++-- src/equations/numerical_fluxes.jl | 2 +- 5 files changed, 6 insertions(+), 6 deletions(-) diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_beach.jl b/examples/tree_1d_dgsem/elixir_shallowwater_beach.jl index 273e8f26e6..2825b215d1 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_beach.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_beach.jl @@ -38,7 +38,7 @@ boundary_condition = boundary_condition_slip_wall # Get the DG approximation space volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) -surface_flux = (FluxHydrostaticReconstruction(flux_hll_cn, hydrostatic_reconstruction_chen_noelle), +surface_flux = (FluxHydrostaticReconstruction(flux_hll_chen_noelle, hydrostatic_reconstruction_chen_noelle), flux_nonconservative_chen_noelle) basis = LobattoLegendreBasis(7) diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_parabolic_bowl.jl b/examples/tree_1d_dgsem/elixir_shallowwater_parabolic_bowl.jl index 724e5d96a6..6d58aa3e73 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_parabolic_bowl.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_parabolic_bowl.jl @@ -42,7 +42,7 @@ initial_condition = initial_condition_parabolic_bowl # Get the DG approximation space volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) -surface_flux = (FluxHydrostaticReconstruction(flux_hll_cn, hydrostatic_reconstruction_chen_noelle), +surface_flux = (FluxHydrostaticReconstruction(flux_hll_chen_noelle, hydrostatic_reconstruction_chen_noelle), flux_nonconservative_chen_noelle) basis = LobattoLegendreBasis(5) diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl b/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl index b63d13433b..280e85de9d 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl @@ -32,7 +32,7 @@ initial_condition = initial_condition_well_balanced_CN # Get the DG approximation space volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) -surface_flux = (FluxHydrostaticReconstruction(flux_hll_cn, hydrostatic_reconstruction_chen_noelle), +surface_flux = (FluxHydrostaticReconstruction(flux_hll_chen_noelle, hydrostatic_reconstruction_chen_noelle), flux_nonconservative_chen_noelle) basis = LobattoLegendreBasis(5) diff --git a/src/Trixi.jl b/src/Trixi.jl index 493b9d8108..2dfdd1f6cd 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -140,14 +140,14 @@ export LaplaceDiffusion2D, export GradientVariablesPrimitive, GradientVariablesEntropy -export flux, flux_central, flux_lax_friedrichs, flux_hll, flux_hll_cn, flux_hllc, flux_hlle, flux_godunov, +export flux, flux_central, flux_lax_friedrichs, flux_hll, flux_hllc, flux_hlle, flux_godunov, flux_chandrashekar, flux_ranocha, flux_derigs_etal, flux_hindenlang_gassner, flux_nonconservative_powell, flux_kennedy_gruber, flux_shima_etal, flux_ec, flux_fjordholm_etal, flux_nonconservative_fjordholm_etal, flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal, hydrostatic_reconstruction_audusse_etal, flux_nonconservative_audusse_etal, - hydrostatic_reconstruction_chen_noelle, flux_nonconservative_chen_noelle, + hydrostatic_reconstruction_chen_noelle, flux_nonconservative_chen_noelle, flux_hll_chen_noelle, FluxPlusDissipation, DissipationGlobalLaxFriedrichs, DissipationLocalLaxFriedrichs, FluxLaxFriedrichs, max_abs_speed_naive, FluxHLL, min_max_speed_naive, diff --git a/src/equations/numerical_fluxes.jl b/src/equations/numerical_fluxes.jl index 7667f0a31c..66ff71bbef 100644 --- a/src/equations/numerical_fluxes.jl +++ b/src/equations/numerical_fluxes.jl @@ -245,7 +245,7 @@ Base.show(io::IO, numflux::FluxHLL) = print(io, "FluxHLL(", numflux.min_max_spee See [`FluxHLL`](@ref). """ const flux_hll = FluxHLL() -const flux_hll_cn = FluxHLLChenNoelle() +const flux_hll_chen_noelle = FluxHLLChenNoelle() """ From a231148bf7265995c776ac0e537226fda7d93788 Mon Sep 17 00:00:00 2001 From: svengoldberg Date: Sat, 11 Feb 2023 08:13:53 +0100 Subject: [PATCH 11/44] Describing docs for Chen and Noelle HR (SWE 1D) --- src/equations/shallow_water_1d.jl | 45 +++++++++++++++++++++++++++++++ 1 file changed, 45 insertions(+) diff --git a/src/equations/shallow_water_1d.jl b/src/equations/shallow_water_1d.jl index 9d010ef6c2..fe0c5c4d8a 100644 --- a/src/equations/shallow_water_1d.jl +++ b/src/equations/shallow_water_1d.jl @@ -312,6 +312,23 @@ Further details on the hydrostatic reconstruction and its motivation can be foun z) end +""" + flux_nonconservative_chen_noelle(u_ll, u_rr, + orientation::Integer, + equations::ShallowWaterEquations1D) + +Non-symmetric two-point surface flux that discretizes the nonconservative (source) term. +The discretization uses the `hydrostatic_reconstruction_chen_noelle` on the conservative +variables. + +Should be used together with [`FluxHydrostaticReconstruction`](@ref) and +[`hydrostatic_reconstruction_chen_noelle`](@ref) in the surface flux to ensure consistency. + +Further details on the hydrostatic reconstruction and its motivation can be found in +- Guoxian Chen and Sebastian Noelle (2017) + A new hydrostatic reconstruction scheme based on subcell reconstructions + [DOI:10.1137/15M1053074](https://dx.doi.org/10.1137/15M1053074) +""" @inline function flux_nonconservative_chen_noelle(u_ll, u_rr, orientation::Integer, equations::ShallowWaterEquations1D) @@ -443,6 +460,21 @@ Further details on this hydrostatic reconstruction and its motivation can be fou return u_ll_star, u_rr_star end +""" + hydrostatic_reconstruction_chen_noelle(u_ll, u_rr, orientation::Integer, + equations::ShallowWaterEquations1D) + +A particular type of hydrostatic reconstruction on the water height to guarantee well-balancedness +for a general bottom topography [`ShallowWaterEquations1D`](@ref). The reconstructed solution states +`u_ll_star` and `u_rr_star` variables are used to evaluate the surface numerical flux at the interface. +The key idea is a linear reconstruction of the bottom and water height at the interfaces using subcells. +Use in combination with the generic numerical flux routine [`FluxHydrostaticReconstruction`](@ref). + +Further details on this hydrostatic reconstruction and its motivation can be found in +- Guoxian Chen and Sebastian Noelle (2017) + A new hydrostatic reconstruction scheme based on subcell reconstructions + [DOI:10.1137/15M1053074](https://dx.doi.org/10.1137/15M1053074) +""" @inline function hydrostatic_reconstruction_chen_noelle(u_ll, u_rr, equations::ShallowWaterEquations1D) # Unpack left and right water heights and bottom topographies h_ll, _, b_ll = u_ll @@ -531,6 +563,19 @@ end end +""" + min_max_speed_chen_noelle(u_ll, u_rr, orientation::Integer, + equations::ShallowWaterEquations1D) + +The approximated speeds for the HLL type numerical flux used by Chen and Noelle for their +hydrostatic reconstruction. As they state in the paper, those speeds are chosen for the numerical +flux to show positivity and an entropy inequality. + +Further details on this hydrostatic reconstruction and its motivation can be found in +- Guoxian Chen and Sebastian Noelle (2017) + A new hydrostatic reconstruction scheme based on subcell reconstructions + [DOI:10.1137/15M1053074](https://dx.doi.org/10.1137/15M1053074) +""" @inline function min_max_speed_chen_noelle(u_ll, u_rr, orientation::Integer, equations::ShallowWaterEquations1D) # Get the velocity quantities v_ll = velocity(u_ll, equations) From 3bb0447345d29f6063cb3ada855254e2f17fa8cf Mon Sep 17 00:00:00 2001 From: svengoldberg Date: Sat, 11 Feb 2023 08:20:26 +0100 Subject: [PATCH 12/44] Edit SWE 1D elixirs, error-based solver and docs --- .../elixir_shallowwater_beach.jl | 47 ++++++++++------- .../elixir_shallowwater_parabolic_bowl.jl | 51 ++++++++++--------- ...ixir_shallowwater_well_balanced_wet_dry.jl | 43 +++++++++++----- 3 files changed, 84 insertions(+), 57 deletions(-) diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_beach.jl b/examples/tree_1d_dgsem/elixir_shallowwater_beach.jl index 2825b215d1..62f688c364 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_beach.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_beach.jl @@ -6,27 +6,42 @@ using Trixi # Semidiscretization of the shallow water equations equations = ShallowWaterEquations1D(gravity_constant=9.812) -cfl = 0.2 +""" + initial_condition_beach(x, t, equations:: ShallowWaterEquations1D) +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 water height and speed function used here, are an adaption of 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) +The used bottom topography differs from the one out of the paper to be differentiable +""" function initial_condition_beach(x, t, equations:: ShallowWaterEquations1D) D = 1 - δ = 0.02 - γ = sqrt((3*δ)/(4*D)) - xₐ = sqrt((4*D)/(3*δ)) * acosh(sqrt(20)) + delta = 0.02 + gamma = sqrt((3 * delta) / (4 * D)) + x_a = sqrt((4 * D) / (3 * delta)) * acosh(sqrt(20)) - f = D + 40*δ * sech(γ*(8*x[1]-xₐ))^2 + f = D + 40 * delta * sech(gamma * (8 * x[1] - x_a))^2 # steep wall b = 0.01 + 99/409600 * 4^x[1] if x[1] >= 6 H = b - v = 0 + 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 + # (1e-13 per default, set in the constructor for the ShallowWaterEquations) is added 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 @@ -41,11 +56,11 @@ volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) surface_flux = (FluxHydrostaticReconstruction(flux_hll_chen_noelle, hydrostatic_reconstruction_chen_noelle), flux_nonconservative_chen_noelle) -basis = LobattoLegendreBasis(7) +basis = LobattoLegendreBasis(3) indicator_sc = IndicatorHennemannGassner(equations, basis, - alpha_max=.5, - alpha_min=.001, + alpha_max=0.5, + alpha_min=0.001, alpha_smooth=true, variable=waterheight_pressure) volume_integral = VolumeIntegralShockCapturingHG(indicator_sc; @@ -61,10 +76,9 @@ coordinates_min = 0.0 coordinates_max = 8.0 mesh = TreeMesh(coordinates_min, coordinates_max, - initial_refinement_level=8, + initial_refinement_level=7, n_cells_max=10_000, - periodicity=false - ) + periodicity=false) # create the semi discretization object semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, @@ -86,14 +100,11 @@ analysis_callback = AnalysisCallback(semi, interval=analysis_interval, save_anal alive_callback = AliveCallback(analysis_interval=analysis_interval) -save_solution = SaveSolutionCallback(interval=1000, +save_solution = SaveSolutionCallback(interval=250, save_initial_solution=true, save_final_solution=true) -stepsize_callback = StepsizeCallback(cfl=cfl) - -callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, - stepsize_callback) +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution) stage_limiter! = PositivityPreservingLimiterZhangShu(thresholds=(equations.threshold_limiter,), variables=(Trixi.waterheight,)) @@ -102,6 +113,6 @@ stage_limiter! = PositivityPreservingLimiterZhangShu(thresholds=(equations.thres # run the simulation sol = solve(ode, SSPRK43(stage_limiter!), dt=1.0, - save_everystep=false, callback=callbacks, adaptive=false, maxiters=1e+9); + save_everystep=false, callback=callbacks); summary_callback() # print the timer summary \ No newline at end of file diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_parabolic_bowl.jl b/examples/tree_1d_dgsem/elixir_shallowwater_parabolic_bowl.jl index 6d58aa3e73..d1e1c4a655 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_parabolic_bowl.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_parabolic_bowl.jl @@ -6,32 +6,37 @@ using Trixi # Semidiscretization of the shallow water equations equations = ShallowWaterEquations1D(gravity_constant=9.81) -cfl = 0.3 - -function parabolic_bowl_analytic_1D_H(gravity, x,t) - a = 1 - h_0 = 0.1 - σ = 0.5 - ω = sqrt(2*gravity*h_0)/a - - H = σ * h_0/a^2 * (2*x[1]*cos(ω*t) - σ) + h_0 - return H -end - -# Implemented based on Wintermeyer (here seen in 2D) +""" + initial_condition_parabolic_bowl(x, t, equations:: ShallowWaterEquations1D) + +Well-known initial condition to test the [`hydrostatic_reconstruction_chen_noelle`](@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 are taken out of section 6.2 from the paper: + - 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\n + [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:: ShallowWaterEquations1D) a = 1 h_0 = 0.1 - σ = 0.5 - ω = sqrt(2*equations.gravity*h_0)/a + sigma = 0.5 + ω = sqrt(2 * equations.gravity * h_0) / a - v = -σ*ω*sin(ω*t) + v = -sigma * ω * sin(ω * t) - b = h_0 * x[1]^2/a^2 + b = h_0 * x[1]^2 / a^2 - H = max(b, parabolic_bowl_analytic_1D_H(equations.gravity, x, t)) + 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 + # (1e-13 per default, set in the constructor for the ShallowWaterEquations) is added 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 @@ -66,8 +71,7 @@ coordinates_max = 2.0 mesh = TreeMesh(coordinates_min, coordinates_max, initial_refinement_level=10, - n_cells_max=10_000 - ) + n_cells_max=10_000) # create the semi discretization object semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) @@ -92,10 +96,7 @@ save_solution = SaveSolutionCallback(interval=1000, save_initial_solution=true, save_final_solution=true) -stepsize_callback = StepsizeCallback(cfl=cfl) - -callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, - stepsize_callback) +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution) stage_limiter! = PositivityPreservingLimiterZhangShu(thresholds=(equations.threshold_limiter,), variables=(Trixi.waterheight,)) @@ -104,6 +105,6 @@ stage_limiter! = PositivityPreservingLimiterZhangShu(thresholds=(equations.thres # run the simulation sol = solve(ode, SSPRK43(stage_limiter!), dt=1.0, - save_everystep=false, callback=callbacks, adaptive=false, maxiters=1e+9); + save_everystep=false, callback=callbacks); summary_callback() # print the timer summary \ No newline at end of file diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl b/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl index 280e85de9d..c44ff715fd 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl @@ -6,14 +6,25 @@ using Trixi # Semidiscretization of the shallow water equations equations = ShallowWaterEquations1D(gravity_constant=9.812) -cfl = .3 -function initial_condition_well_balanced_CN(x, t, equations:: ShallowWaterEquations1D) - v = 0 - b = sin(4*pi*x[1]) + 3 +""" + initial_condition_well_balanced_chen_noelle(x, t, equations:: ShallowWaterEquations1D) + +Initial condition with a complex (discontinuous) bottom topography to test the well-balanced +property for the [`hydrostatic_reconstruction_chen_noelle`](@ref) including dry areas within the +domain. The errors from the analysis callback are not important but the error for this +lake at rest test case `∑|H0-(h+b)|` should be around machine roundoff. +The initial condition was found in the section 5.2 of the paper: +- Guoxian Chen and Sebastian Noelle (2017) + A new hydrostatic reconstruction scheme based on subcell reconstructions + [DOI:10.1137/15M1053074](https://dx.doi.org/10.1137/15M1053074) +""" +function initial_condition_well_balanced_chen_noelle(x, t, equations:: ShallowWaterEquations1D) + v = 0.0 + b = sin(4 * pi * x[1]) + 3 if x[1] >= 0.5 - b = sin(4*pi*x[1]) + 1 + b = sin(4 * pi * x[1]) + 1 end H = max(b, 2.5) @@ -22,11 +33,16 @@ function initial_condition_well_balanced_CN(x, t, equations:: ShallowWaterEquati H = max(b, 1.5) 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 + # (1e-13 per default, set in the constructor for the ShallowWaterEquations) is added 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_well_balanced_CN +initial_condition = initial_condition_well_balanced_chen_noelle ############################################################################### # Get the DG approximation space @@ -35,7 +51,7 @@ volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) surface_flux = (FluxHydrostaticReconstruction(flux_hll_chen_noelle, hydrostatic_reconstruction_chen_noelle), flux_nonconservative_chen_noelle) -basis = LobattoLegendreBasis(5) +basis = LobattoLegendreBasis(3) indicator_sc = IndicatorHennemannGassner(equations, basis, alpha_max=0.5, @@ -51,13 +67,12 @@ solver = DGSEM(basis, surface_flux, volume_integral) ############################################################################### # Create the TreeMesh for the domain [0, 1] -coordinates_min = 0. -coordinates_max = 1. +coordinates_min = 0.0 +coordinates_max = 1.0 mesh = TreeMesh(coordinates_min, coordinates_max, initial_refinement_level=4, - n_cells_max=10_000, - ) + n_cells_max=10_000) # create the semi discretization object semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) @@ -65,7 +80,7 @@ semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) ############################################################################### # ODE solvers, callbacks etc. -tspan = (0.0, 25.0) +tspan = (0.0, 100.0) ode = semidiscretize(semi, tspan) @@ -79,11 +94,11 @@ analysis_callback = AnalysisCallback(semi, interval=analysis_interval, save_anal alive_callback = AliveCallback(analysis_interval=analysis_interval) -save_solution = SaveSolutionCallback(interval=1000, +save_solution = SaveSolutionCallback(interval=5000, save_initial_solution=true, save_final_solution=true) -stepsize_callback = StepsizeCallback(cfl=cfl) +stepsize_callback = StepsizeCallback(cfl=2.0) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, stepsize_callback) From 449af840c397dd0b5d2d625ec3c2cbac558558c3 Mon Sep 17 00:00:00 2001 From: svengoldberg Date: Sat, 11 Feb 2023 08:23:35 +0100 Subject: [PATCH 13/44] Including tests on new SWE 1D elixirs --- test/test_tree_1d_shallowwater.jl | 21 +++++++++++++++++++++ 1 file changed, 21 insertions(+) diff --git a/test/test_tree_1d_shallowwater.jl b/test/test_tree_1d_shallowwater.jl index d62d19c2a8..20f086bd02 100644 --- a/test/test_tree_1d_shallowwater.jl +++ b/test/test_tree_1d_shallowwater.jl @@ -30,6 +30,13 @@ EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "tree_1 tspan = (0.0, 0.25)) end + @trixi_testset "elixir_shallowwater_well_balanced_wet_dry.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_well_balanced_wet_dry.jl"), + l2 = [0.027869290160397523, 4.848702566024514e-14, 0.10911645796834736], + linf = [0.28423245872819947, 2.305820093827649e-13, 1.119394067096447], + tspan = (0.0, 0.25)) + end + @trixi_testset "elixir_shallowwater_source_terms.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"), l2 = [0.0022363707373868713, 0.01576799981934617, 4.436491725585346e-5], @@ -80,6 +87,20 @@ EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "tree_1 linf = [0.7565706154863958, 2.076621603471687, 0.8646939843534258], tspan = (0.0, 0.05)) end + + @trixi_testset "elixir_shallowwater_beach.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_beach.jl"), + l2 = [0.5837838749990294, 3.740552325994961, 6.289710785696625e-8], + linf = [0.9917232563137562, 7.096291692934598, 4.452604027704865e-7], + tspan = (0.0, 0.25)) + end + + @trixi_testset "elixir_shallowwater_parabolic_bowl.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_parabolic_bowl.jl"), + l2 = [0.012868413617725416, 0.007819724355957678, 3.8569975327636153e-17], + linf = [0.01910483682725715, 0.011878853710648146, 2.7755575615628914e-16], + tspan = (0.0, 0.75)) + end end end # module From 15930b19232ff7d59d0484222ee7a321ef3fb948 Mon Sep 17 00:00:00 2001 From: svengoldberg Date: Thu, 23 Feb 2023 14:43:48 +0100 Subject: [PATCH 14/44] New/reorganize positivity limiter (SWE 1D) --- .../elixir_shallowwater_beach.jl | 2 +- .../elixir_shallowwater_parabolic_bowl.jl | 2 +- ...ixir_shallowwater_well_balanced_wet_dry.jl | 2 +- src/Trixi.jl | 2 +- src/callbacks_stage/callbacks_stage.jl | 1 + .../positivity_shallow_water.jl | 82 ++++++++++++++++++ .../positivity_shallow_water_dg1d.jl | 86 +++++++++++++++++++ .../positivity_zhang_shu_dg1d.jl | 79 ----------------- 8 files changed, 173 insertions(+), 83 deletions(-) create mode 100644 src/callbacks_stage/positivity_shallow_water.jl create mode 100644 src/callbacks_stage/positivity_shallow_water_dg1d.jl diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_beach.jl b/examples/tree_1d_dgsem/elixir_shallowwater_beach.jl index 62f688c364..0110b2817f 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_beach.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_beach.jl @@ -106,7 +106,7 @@ save_solution = SaveSolutionCallback(interval=250, callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution) -stage_limiter! = PositivityPreservingLimiterZhangShu(thresholds=(equations.threshold_limiter,), +stage_limiter! = PositivityPreservingLimiterShallowWater(thresholds=(equations.threshold_limiter,), variables=(Trixi.waterheight,)) ############################################################################### diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_parabolic_bowl.jl b/examples/tree_1d_dgsem/elixir_shallowwater_parabolic_bowl.jl index d1e1c4a655..dedfad1586 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_parabolic_bowl.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_parabolic_bowl.jl @@ -98,7 +98,7 @@ save_solution = SaveSolutionCallback(interval=1000, callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution) -stage_limiter! = PositivityPreservingLimiterZhangShu(thresholds=(equations.threshold_limiter,), +stage_limiter! = PositivityPreservingLimiterShallowWater(thresholds=(equations.threshold_limiter,), variables=(Trixi.waterheight,)) ############################################################################### diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl b/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl index c44ff715fd..05b2e2d617 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl @@ -103,7 +103,7 @@ stepsize_callback = StepsizeCallback(cfl=2.0) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, stepsize_callback) -stage_limiter! = PositivityPreservingLimiterZhangShu(thresholds=(equations.threshold_limiter,), +stage_limiter! = PositivityPreservingLimiterShallowWater(thresholds=(equations.threshold_limiter,), variables=(Trixi.waterheight,)) ############################################################################### diff --git a/src/Trixi.jl b/src/Trixi.jl index 0e83245e22..1fb8c1246c 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -227,7 +227,7 @@ export ControllerThreeLevel, ControllerThreeLevelCombined, IndicatorLöhner, IndicatorLoehner, IndicatorMax, IndicatorNeuralNetwork, NeuralNetworkPerssonPeraire, NeuralNetworkRayHesthaven, NeuralNetworkCNN -export PositivityPreservingLimiterZhangShu +export PositivityPreservingLimiterZhangShu, PositivityPreservingLimiterShallowWater export trixi_include, examples_dir, get_examples, default_example, default_example_unstructured diff --git a/src/callbacks_stage/callbacks_stage.jl b/src/callbacks_stage/callbacks_stage.jl index f23f96eccf..9c99fc8ffd 100644 --- a/src/callbacks_stage/callbacks_stage.jl +++ b/src/callbacks_stage/callbacks_stage.jl @@ -6,6 +6,7 @@ include("positivity_zhang_shu.jl") +include("positivity_shallow_water.jl") end # @muladd diff --git a/src/callbacks_stage/positivity_shallow_water.jl b/src/callbacks_stage/positivity_shallow_water.jl new file mode 100644 index 0000000000..e0cda95913 --- /dev/null +++ b/src/callbacks_stage/positivity_shallow_water.jl @@ -0,0 +1,82 @@ +# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). + # Since these FMAs can increase the performance of many numerical algorithms, + # we need to opt-in explicitly. + # See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. + @muladd begin + + + """ + PositivityPreservingLimiterShallowWater(; threshold, variables) + + 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) + + It is specifically designed for `ShallowWaterEquations`. + The limiter is applied to all scalar `variables` in their given order + using the associated `thresholds` to determine the minimal acceptable values. + The order of the `variables` is important and might have a strong influence + on the robustness. + As opposed to the default version of the limiter, nodes with a waterheight below the threshold + are treated different. To avoid numerical problems caused by velocities close to zero in that case, + the velocity is cut off, as well, as the node can be identified as dry. The special feature of the + `ShallowWaterEquations`` used here is that the bottom topography is stored in `u`. However, this value + should not be changed, much less limited. That is why, the `set_node_vars!` function is not used + for the last conservation variable. + Afterwards for safety reasons, the wet/dry threshold is applied again on all the DG nodes of freedom + in order to avoid dry nodes. In the case where value_mean < threshold before applying + limiter, there could still be dry nodes afterwards due to the logic of the limiter. + """ + struct PositivityPreservingLimiterShallowWater{N, Thresholds<:NTuple{N,<:Real}, Variables<:NTuple{N,Any}} + thresholds::Thresholds + variables::Variables + end + + function PositivityPreservingLimiterShallowWater(; thresholds, variables) + PositivityPreservingLimiterShallowWater(thresholds, variables) + end + + + function (limiter!::PositivityPreservingLimiterShallowWater)( + u_ode, integrator, semi::AbstractSemidiscretization, t) + u = wrap_array(u_ode, semi) + @trixi_timeit timer() "positivity-preserving limiter" limiter_shallow_water!( + u, limiter!.thresholds, limiter!.variables, mesh_equations_solver_cache(semi)...) + end + + + # Iterate over tuples in a type-stable way using "lispy tuple programming", + # similar to https://stackoverflow.com/a/55849398: + # Iterating over tuples of different functions isn't type-stable in general + # but accessing the first element of a tuple is type-stable. Hence, it's good + # to process one element at a time and replace iteration by recursion here. + # Note that you shouldn't use this with too many elements per tuple since the + # compile times can increase otherwise - but a handful of elements per tuple + # is definitely fine. + function limiter_shallow_water!(u, thresholds::NTuple{N,<:Real}, variables::NTuple{N,Any}, + mesh, equations::Union{ShallowWaterEquations1D, ShallowWaterEquations2D}, + solver, cache) where {N} + threshold = first(thresholds) + remaining_thresholds = Base.tail(thresholds) + variable = first(variables) + remaining_variables = Base.tail(variables) + + limiter_shallow_water!(u, threshold, variable, mesh, equations, solver, cache) + limiter_shallow_water!(u, remaining_thresholds, remaining_variables, mesh, equations, solver, cache) + return nothing + end + + # terminate the type-stable iteration over tuples + function limiter_shallow_water!(u, thresholds::Tuple{}, variables::Tuple{}, + mesh, equations::Union{ShallowWaterEquations1D, ShallowWaterEquations2D}, + solver, cache) + nothing + end + + include("positivity_shallow_water_dg1d.jl") + + + + end # @muladd \ No newline at end of file diff --git a/src/callbacks_stage/positivity_shallow_water_dg1d.jl b/src/callbacks_stage/positivity_shallow_water_dg1d.jl new file mode 100644 index 0000000000..54bece7676 --- /dev/null +++ b/src/callbacks_stage/positivity_shallow_water_dg1d.jl @@ -0,0 +1,86 @@ +# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). +# Since these FMAs can increase the performance of many numerical algorithms, +# we need to opt-in explicitly. +# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. +@muladd begin + + +function limiter_shallow_water!(u, threshold::Real, variable, + mesh::AbstractMesh{1}, equations::ShallowWaterEquations1D, + dg::DGSEM, cache) + @unpack weights = dg.basis + + @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 || 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, as we assume that + # Jensen's inequality holds (e.g. pressure for compressible Euler equations). + 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) + + # Cut off velocity in case that the waterheight is smaller than the threshold + + h_node, h_v_node, b_node = u_node + h_mean, h_v_mean, _ = u_mean # b_mean is not used as b_node must not be overwritten + + # Set them both to zero to apply linear combination correctly + h_v_node = h_v_node * Int32(h_node > threshold) + + h_v_mean = h_v_mean * Int32(h_node > threshold) + + u_node = SVector(h_node, h_v_node, b_node) + u_mean = SVector(h_mean, h_v_mean, zero(eltype(u))) + + # When velocity is cut off, the only value beeing averaged is the waterheight, + # because the velocity is set to zero and this value is passed. + # Otherwise, the velocity is averaged, as well. + # Note that the auxiliary bottom topography variable `b` is never limited. + set_node_vars!(u, theta * u_node + (1-theta) * u_mean, + (1,2), dg, i, element) + end + end + + # An extra "safety" check is done over all the degrees of + # freedom after the limiting 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 limiter. + @threaded for element in eachelement(dg, cache) + + for i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, element) + + h, v, b = u_node + + h = h * Int32(h > threshold) + threshold * Int32(h <= threshold) + v = v * Int32(h > threshold) + + u_node = SVector(h, v, 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_zhang_shu_dg1d.jl b/src/callbacks_stage/positivity_zhang_shu_dg1d.jl index 094e04b03f..50d6b3f2c3 100644 --- a/src/callbacks_stage/positivity_zhang_shu_dg1d.jl +++ b/src/callbacks_stage/positivity_zhang_shu_dg1d.jl @@ -43,84 +43,5 @@ function limiter_zhang_shu!(u, threshold::Real, variable, return nothing end -# Overload limiter for ShallowWaterEquations1D and cut off waterheight as well as velocity -# to avoid numerical problems/instabilities caused by velocities close to zero -function limiter_zhang_shu!(u, threshold::Real, variable, - mesh::AbstractMesh{1}, equations::ShallowWaterEquations1D, - dg::DGSEM, cache) - @unpack weights = dg.basis - - @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 || 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, as we assume that - # Jensen's inequality holds (e.g. pressure for compressible Euler equations). - 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) - - # Cut off velocity in case that the waterheight is smaller than the threshold - - h_node, h_v_node, b_node = u_node - h_mean, h_v_mean, _ = u_mean # b_mean is not used as b_node must not be overwritten - - # Set them both to zero to apply linear combination correctly - h_v_node = h_v_node * Int32(h_node > threshold) - - h_v_mean = h_v_mean * Int32(h_node > threshold) - - u_node = SVector(h_node, h_v_node, b_node) - u_mean = SVector(h_mean, h_v_mean, zero(eltype(u))) - - # When velocity is cut off, the only value beeing averaged is the waterheight, - # because the velocity is set to zero and this value is passed. - # Otherwise, the velocity is averaged, as well. - # Note that the auxiliary bottom topography variable `b` is never limited. - set_node_vars!(u, theta * u_node + (1-theta) * u_mean, - (1,2), dg, i, element) - end - end - - # An extra "safety" check is done over all the degrees of - # freedom after the limiting 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 limiter. - @threaded for element in eachelement(dg, cache) - - for i in eachnode(dg) - u_node = get_node_vars(u, equations, dg, i, element) - - h, v, b = u_node - - h = h * Int32(h > threshold) + threshold * Int32(h <= threshold) - v = v * Int32(h > threshold) - - u_node = SVector(h, v, b) - - set_node_vars!(u, u_node, equations, dg, i, element) - end - end - - return nothing -end - end # @muladd From b0bb88d718ec22012b99434a73a5f5ac823d68c3 Mon Sep 17 00:00:00 2001 From: svengoldberg Date: Thu, 2 Mar 2023 22:21:05 +0100 Subject: [PATCH 15/44] Editing docs SWE 1D --- .../positivity_shallow_water.jl | 21 +++++++++++-------- src/equations/shallow_water_1d.jl | 10 +++++++-- 2 files changed, 20 insertions(+), 11 deletions(-) diff --git a/src/callbacks_stage/positivity_shallow_water.jl b/src/callbacks_stage/positivity_shallow_water.jl index e0cda95913..dd3171826c 100644 --- a/src/callbacks_stage/positivity_shallow_water.jl +++ b/src/callbacks_stage/positivity_shallow_water.jl @@ -19,15 +19,18 @@ using the associated `thresholds` to determine the minimal acceptable values. The order of the `variables` is important and might have a strong influence on the robustness. - As opposed to the default version of the limiter, nodes with a waterheight below the threshold - are treated different. To avoid numerical problems caused by velocities close to zero in that case, - the velocity is cut off, as well, as the node can be identified as dry. The special feature of the - `ShallowWaterEquations`` used here is that the bottom topography is stored in `u`. However, this value - should not be changed, much less limited. That is why, the `set_node_vars!` function is not used - for the last conservation variable. - Afterwards for safety reasons, the wet/dry threshold is applied again on all the DG nodes of freedom - in order to avoid dry nodes. In the case where value_mean < threshold before applying - limiter, there could still be dry nodes afterwards due to the logic of the limiter. + As opposed to the standard version of the [PositivityPreservingLimiterZhangShu](@ref), + nodes with a water height below the `threshold_limiter` are treated in a special way. + To avoid numerical problems caused by velocities close to zero, + the velocity is cut off, such that the node can be identified as "dry". The special feature of the + `ShallowWaterEquations` used here is that the bottom topography is stored as an additional + quantity in the solution vector `u`. However, the value of the bottom topography + should not be changed, much less limited. That is why, the `set_node_vars!` function is not applied + to the last conservation variable. + After the limiting process is applied to all degrees of freedom, for safety reasons, + the wet/dry threshold is applied again on all the DG nodes in order to avoid dry nodes. + In the case where value_mean < threshold before applying imiter, there could still be dry nodes + afterwards due to the logic of the limiter. """ struct PositivityPreservingLimiterShallowWater{N, Thresholds<:NTuple{N,<:Real}, Variables<:NTuple{N,Any}} thresholds::Thresholds diff --git a/src/equations/shallow_water_1d.jl b/src/equations/shallow_water_1d.jl index fe0c5c4d8a..13a02dd6ca 100644 --- a/src/equations/shallow_water_1d.jl +++ b/src/equations/shallow_water_1d.jl @@ -6,7 +6,7 @@ @doc raw""" - ShallowWaterEquations1D(gravity, H0) + ShallowWaterEquations1D(gravity, H0, threshold_limiter, threshold_wet) Shallow water equations (SWE) in one space dimension. The equations are given by ```math @@ -24,6 +24,12 @@ also defines the total water height as ``H = h + b``. 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. Both of them do not +have to be passed, as default values are defined within the struct. The first one, threshold_limiter, is +used in PositivityPreservingLimiterShallowWater 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. + 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 variable `b` to zero. @@ -47,7 +53,7 @@ References for the SWE are many but a good introduction is available in Chapter struct ShallowWaterEquations1D{RealT<:Real} <: AbstractShallowWaterEquations{1, 3} gravity::RealT # gravitational constant H0::RealT # constant "lake-at-rest" total water height - threshold_limiter::RealT # Threshold to use in PositivityPreservingLimiterZhangShu on waterheight, + threshold_limiter::RealT # Threshold to use in PositivityPreservingLimiterShallowWater on water height, # as a (small) shift on the initial condition and cutoff before the # next time step. threshold_wet::RealT # Threshold to be applied on water height to define when the flow is "wet" From ef0e6e684021d9a879e7b5d2b485d993d51ebab0 Mon Sep 17 00:00:00 2001 From: svengoldberg Date: Fri, 3 Mar 2023 18:02:06 +0100 Subject: [PATCH 16/44] Rearrange cut off at interfaces, edit tests SWE 1D --- src/equations/numerical_fluxes.jl | 22 +--------------------- src/equations/shallow_water_1d.jl | 17 +++++++++++++++-- test/test_tree_1d_shallowwater.jl | 8 ++++---- 3 files changed, 20 insertions(+), 27 deletions(-) diff --git a/src/equations/numerical_fluxes.jl b/src/equations/numerical_fluxes.jl index 66ff71bbef..8e9bbb4487 100644 --- a/src/equations/numerical_fluxes.jl +++ b/src/equations/numerical_fluxes.jl @@ -307,32 +307,12 @@ end @inline function (numflux::FluxHydrostaticReconstruction)(u_ll, u_rr, orientation_or_normal_direction, - equations::AbstractEquations{1}) + equations::AbstractEquations) @unpack numerical_flux, hydrostatic_reconstruction = numflux # Create the reconstructed left/right solution states in conservative form u_ll_star, u_rr_star = hydrostatic_reconstruction(u_ll, u_rr, equations) - # Set the water height to be at least the value stored in the variable threshold after - # the hydrostatic reconstruction is applied and before the numerical flux is calculated - # to avoid numerical problem with arbitrary small values. Interfaces with a water height - # lower or equal to the threshold can be declared as dry. - # The default value is set to 1e-15 and can be changed within the constructor call in an elixir. - threshold = equations.threshold_wet - - # Apply threshold to cut off the height and the velocity at dry interfaces - if threshold > 0 - h_ll, _, b_ll = u_ll_star - h_rr, _, b_rr = u_rr_star - - if h_ll <= threshold - u_ll_star = SVector(threshold, zero(eltype(u_ll)), b_ll) - end - - if h_rr <= threshold - u_rr_star = SVector(threshold, zero(eltype(u_rr)), b_rr) - end - end # Use the reconstructed states to compute the numerical surface flux return numerical_flux(u_ll_star, u_rr_star, orientation_or_normal_direction, equations) end diff --git a/src/equations/shallow_water_1d.jl b/src/equations/shallow_water_1d.jl index 13a02dd6ca..41f25f31bf 100644 --- a/src/equations/shallow_water_1d.jl +++ b/src/equations/shallow_water_1d.jl @@ -490,8 +490,8 @@ Further details on this hydrostatic reconstruction and its motivation can be fou v_ll = velocity(u_ll, equations) v_rr = velocity(u_rr, equations) - H_ll = b_ll+h_ll - H_rr = b_rr+h_rr + H_ll = b_ll + h_ll + H_rr = b_rr + h_rr b_star = min( max( b_ll, b_rr ), min( H_ll, H_rr ) ) @@ -499,6 +499,19 @@ Further details on this hydrostatic reconstruction and its motivation can be fou h_ll_star = min( H_ll - b_star, h_ll ) h_rr_star = min( H_rr - b_star, h_rr ) + # Set the water height to be at least the value stored in the variable threshold after + # the hydrostatic reconstruction is applied and before the numerical flux is calculated + # to avoid numerical problem with arbitrary small values. Interfaces with a water height + # lower or equal to the threshold can be declared as dry. + # The default value is set to 1e-15 and can be changed within the constructor call in an elixir. + threshold = equations.threshold_wet + + h_ll_star = h_ll_star * Int32(h_ll_star > threshold) + threshold * Int32(h_ll_star <= threshold) + h_rr_star = h_rr_star * Int32(h_rr_star > threshold) + threshold * Int32(h_rr_star <= threshold) + + v_ll = v_ll * Int32(h_ll_star > threshold) + v_rr = v_rr * Int32(h_rr_star > threshold) + # Create the conservative variables using the reconstruted water heights u_ll_star = SVector( h_ll_star, h_ll_star * v_ll, b_ll ) u_rr_star = SVector( h_rr_star, h_rr_star * v_rr, b_rr ) diff --git a/test/test_tree_1d_shallowwater.jl b/test/test_tree_1d_shallowwater.jl index 20f086bd02..2bb7ae469b 100644 --- a/test/test_tree_1d_shallowwater.jl +++ b/test/test_tree_1d_shallowwater.jl @@ -90,15 +90,15 @@ EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "tree_1 @trixi_testset "elixir_shallowwater_beach.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_beach.jl"), - l2 = [0.5837838749990294, 3.740552325994961, 6.289710785696625e-8], - linf = [0.9917232563137562, 7.096291692934598, 4.452604027704865e-7], + l2 = [0.5837839017922798, 3.7405541768121595, 6.289710785696625e-8], + linf = [0.9917266589187276, 7.096283536176454, 4.452604027704865e-7], tspan = (0.0, 0.25)) end @trixi_testset "elixir_shallowwater_parabolic_bowl.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_parabolic_bowl.jl"), - l2 = [0.012868413617725416, 0.007819724355957678, 3.8569975327636153e-17], - linf = [0.01910483682725715, 0.011878853710648146, 2.7755575615628914e-16], + l2 = [0.012868412255957612, 0.007819719632793341, 3.8569975327636153e-17], + linf = [0.01898697414997512, 0.011686441632120272, 2.7755575615628914e-16], tspan = (0.0, 0.75)) end end From 022e0ee59a35db766ad45422132805473fdb9159 Mon Sep 17 00:00:00 2001 From: svengoldberg Date: Fri, 3 Mar 2023 18:06:21 +0100 Subject: [PATCH 17/44] Edit docs, add Ref --- src/callbacks_stage/positivity_shallow_water.jl | 8 ++++---- src/equations/shallow_water_1d.jl | 4 ++-- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/callbacks_stage/positivity_shallow_water.jl b/src/callbacks_stage/positivity_shallow_water.jl index dd3171826c..67fb9ec7e5 100644 --- a/src/callbacks_stage/positivity_shallow_water.jl +++ b/src/callbacks_stage/positivity_shallow_water.jl @@ -18,9 +18,9 @@ The limiter is applied to all scalar `variables` in their given order using the associated `thresholds` to determine the minimal acceptable values. The order of the `variables` is important and might have a strong influence - on the robustness. - As opposed to the standard version of the [PositivityPreservingLimiterZhangShu](@ref), - nodes with a water height below the `threshold_limiter` are treated in a special way. + on the robustness. + As opposed to the standard version of the [`PositivityPreservingLimiterZhangShu`](@ref), + nodes with a water height below the `threshold_limiter` are treated in a special way. To avoid numerical problems caused by velocities close to zero, the velocity is cut off, such that the node can be identified as "dry". The special feature of the `ShallowWaterEquations` used here is that the bottom topography is stored as an additional @@ -28,7 +28,7 @@ should not be changed, much less limited. That is why, the `set_node_vars!` function is not applied to the last conservation variable. After the limiting process is applied to all degrees of freedom, for safety reasons, - the wet/dry threshold is applied again on all the DG nodes in order to avoid dry nodes. + the wet/dry threshold is applied again on all the DG nodes in order to avoid dry nodes. In the case where value_mean < threshold before applying imiter, there could still be dry nodes afterwards due to the logic of the limiter. """ diff --git a/src/equations/shallow_water_1d.jl b/src/equations/shallow_water_1d.jl index 41f25f31bf..ceb43b0429 100644 --- a/src/equations/shallow_water_1d.jl +++ b/src/equations/shallow_water_1d.jl @@ -26,8 +26,8 @@ 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. Both of them do not have to be passed, as default values are defined within the struct. The first one, threshold_limiter, is -used in PositivityPreservingLimiterShallowWater 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 +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. The bottom topography function ``b(x)`` is set inside the initial condition routine From 0bcf93faa21b2617ada00e184c090955ade4260b Mon Sep 17 00:00:00 2001 From: svengoldberg Date: Fri, 3 Mar 2023 23:29:27 +0100 Subject: [PATCH 18/44] Add 2D lake-at-rest-error logic to pass 2D tests --- src/callbacks_step/analysis_dg2d.jl | 15 ++++++++++++++- src/equations/shallow_water_2d.jl | 11 +++++++++++ 2 files changed, 25 insertions(+), 1 deletion(-) diff --git a/src/callbacks_step/analysis_dg2d.jl b/src/callbacks_step/analysis_dg2d.jl index 453474675f..45b90e6855 100644 --- a/src/callbacks_step/analysis_dg2d.jl +++ b/src/callbacks_step/analysis_dg2d.jl @@ -194,7 +194,20 @@ function integrate(func::Func, u, end end - +#Overload to pass the initial_condition as parameter for calculating lake_at_rest_error + function integrate(::typeof(lake_at_rest_error), u, + mesh::Union{TreeMesh{2}, StructuredMesh{2}, UnstructuredMesh2D, P4estMesh{2}}, + equations::ShallowWaterEquations2D, dg::DG, cache, initial_condition; normalize=true) + node_coordinates = cache.elements.node_coordinates + integrate_via_indices(u, mesh, equations, dg, cache; normalize=normalize) do u, i, j, element, equations, dg + x = get_node_coords(node_coordinates, equations, dg, i, j, element) + + u_exact = initial_condition(x, 0.0, equations) + u_local = get_node_vars(u, equations, dg, i, j, element) + return lake_at_rest_error(u_local, u_exact, equations) + end + end + function analyze(::typeof(entropy_timederivative), du, u, t, mesh::Union{TreeMesh{2}, StructuredMesh{2}, UnstructuredMesh2D, P4estMesh{2}}, equations, dg::DG, cache) diff --git a/src/equations/shallow_water_2d.jl b/src/equations/shallow_water_2d.jl index c78f79c5ab..00e15bbc11 100644 --- a/src/equations/shallow_water_2d.jl +++ b/src/equations/shallow_water_2d.jl @@ -880,4 +880,15 @@ end return abs(equations.H0 - (h + b)) end +@inline function lake_at_rest_error(u, u_exact, equations::ShallowWaterEquations2D) + h, _, _, b = u + + h_exact, _, _, b_exact= u_exact + + H = h + b + H_exact = h_exact + b_exact + + return abs(H - H_exact) +end + end # @muladd From 6906fc6b16b420e94ac4eacfc474c7408c6f9232 Mon Sep 17 00:00:00 2001 From: svengoldberg Date: Wed, 19 Apr 2023 15:29:51 +0200 Subject: [PATCH 19/44] Fixed typo. Confusing name, but correct math --- src/callbacks_stage/positivity_shallow_water_dg1d.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/callbacks_stage/positivity_shallow_water_dg1d.jl b/src/callbacks_stage/positivity_shallow_water_dg1d.jl index 54bece7676..a71f4d0ff0 100644 --- a/src/callbacks_stage/positivity_shallow_water_dg1d.jl +++ b/src/callbacks_stage/positivity_shallow_water_dg1d.jl @@ -68,12 +68,12 @@ function limiter_shallow_water!(u, threshold::Real, variable, for i in eachnode(dg) u_node = get_node_vars(u, equations, dg, i, element) - h, v, b = u_node + h, hv, b = u_node h = h * Int32(h > threshold) + threshold * Int32(h <= threshold) - v = v * Int32(h > threshold) + hv = hv * Int32(h > threshold) - u_node = SVector(h, v, b) + u_node = SVector(h, hv, b) set_node_vars!(u, u_node, equations, dg, i, element) end From 5bfd23ce963c47da314c05c61cdb4e7e986078ef Mon Sep 17 00:00:00 2001 From: svengoldberg Date: Tue, 25 Apr 2023 09:44:43 +0200 Subject: [PATCH 20/44] Correction of comments and docstrings --- src/callbacks_stage/positivity_shallow_water.jl | 2 +- src/equations/shallow_water_1d.jl | 10 +++++----- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/callbacks_stage/positivity_shallow_water.jl b/src/callbacks_stage/positivity_shallow_water.jl index 67fb9ec7e5..2a9a6a60e8 100644 --- a/src/callbacks_stage/positivity_shallow_water.jl +++ b/src/callbacks_stage/positivity_shallow_water.jl @@ -29,7 +29,7 @@ to the last conservation variable. After the limiting process is applied to all degrees of freedom, for safety reasons, the wet/dry threshold is applied again on all the DG nodes in order to avoid dry nodes. - In the case where value_mean < threshold before applying imiter, there could still be dry nodes + In the case where value_mean < threshold before limiting, there could still be dry nodes afterwards due to the logic of the limiter. """ struct PositivityPreservingLimiterShallowWater{N, Thresholds<:NTuple{N,<:Real}, Variables<:NTuple{N,Any}} diff --git a/src/equations/shallow_water_1d.jl b/src/equations/shallow_water_1d.jl index ceb43b0429..69a4408ce2 100644 --- a/src/equations/shallow_water_1d.jl +++ b/src/equations/shallow_water_1d.jl @@ -25,9 +25,9 @@ The additional quantity ``H_0`` is also available to store a reference value for 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. Both of them do not -have to be passed, as default values are defined within the struct. The first one, threshold_limiter, is +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 +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. The bottom topography function ``b(x)`` is set inside the initial condition routine @@ -360,7 +360,7 @@ Further details on the hydrostatic reconstruction and its motivation can be foun # cross-averaging across a discontinuous bottom topography # (ii) True surface part that uses `h_ll` and `h_ll_star` to handle discontinuous bathymetry return SVector(z, - equations.gravity * h_ll * b_ll - equations.gravity * (h_ll_star+h_ll) * (b_ll-b_star), + equations.gravity * h_ll * b_ll - equations.gravity * (h_ll_star + h_ll) * (b_ll - b_star), z) end @@ -587,8 +587,8 @@ end equations::ShallowWaterEquations1D) The approximated speeds for the HLL type numerical flux used by Chen and Noelle for their -hydrostatic reconstruction. As they state in the paper, those speeds are chosen for the numerical -flux to show positivity and an entropy inequality. +hydrostatic reconstruction. As they state in the paper, these speeds are chosen for the numerical +flux to ensure positivity and to satisfy an entropy inequality. Further details on this hydrostatic reconstruction and its motivation can be found in - Guoxian Chen and Sebastian Noelle (2017) From 4a500fbe2946a231232406cd1e116d24f667d1d3 Mon Sep 17 00:00:00 2001 From: svengoldberg Date: Thu, 4 May 2023 12:52:09 +0200 Subject: [PATCH 21/44] Reset lake-at-rest error computation strategy. New version of error only in wet-dry elixir (special case) Update test values as the bottom is now truly discontinuous --- ...ixir_shallowwater_well_balanced_wet_dry.jl | 96 ++++++++++++++++--- src/callbacks_step/analysis.jl | 11 --- src/callbacks_step/analysis_dg1d.jl | 14 --- src/callbacks_step/analysis_dg2d.jl | 14 --- src/equations/shallow_water_1d.jl | 27 +++--- src/equations/shallow_water_2d.jl | 19 ++-- test/test_tree_1d_shallowwater.jl | 6 +- 7 files changed, 109 insertions(+), 78 deletions(-) diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl b/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl index 05b2e2d617..750c092291 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl @@ -11,11 +11,11 @@ equations = ShallowWaterEquations1D(gravity_constant=9.812) initial_condition_well_balanced_chen_noelle(x, t, equations:: ShallowWaterEquations1D) Initial condition with a complex (discontinuous) bottom topography to test the well-balanced -property for the [`hydrostatic_reconstruction_chen_noelle`](@ref) including dry areas within the -domain. The errors from the analysis callback are not important but the error for this +property for the [`hydrostatic_reconstruction_chen_noelle`](@ref) including dry areas within the +domain. The errors from the analysis callback are not important but the error for this lake at rest test case `∑|H0-(h+b)|` should be around machine roundoff. The initial condition was found in the section 5.2 of the paper: -- Guoxian Chen and Sebastian Noelle (2017) +- Guoxian Chen and Sebastian Noelle (2017) A new hydrostatic reconstruction scheme based on subcell reconstructions [DOI:10.1137/15M1053074](https://dx.doi.org/10.1137/15M1053074) """ @@ -34,9 +34,9 @@ function initial_condition_well_balanced_chen_noelle(x, t, equations:: ShallowWa 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 + # 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 - # (1e-13 per default, set in the constructor for the ShallowWaterEquations) is added if h = 0. + # (1e-13 per default, set in the constructor for the ShallowWaterEquations) is added 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) @@ -71,7 +71,7 @@ coordinates_min = 0.0 coordinates_max = 1.0 mesh = TreeMesh(coordinates_min, coordinates_max, - initial_refinement_level=4, + initial_refinement_level=6, n_cells_max=10_000) # create the semi discretization object @@ -80,17 +80,41 @@ semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) ############################################################################### # ODE solvers, callbacks etc. -tspan = (0.0, 100.0) +tspan = (0.0, 25.0) ode = semidiscretize(semi, tspan) +function initial_condition_discontinuous_well_balancedness(x, t, element_id, equations:: ShallowWaterEquations1D) + v = 0.0 + b = sin(4 * pi * x[1]) + 1 + H = max(b, 1.5) + + # for refinement_level=6 + if (element_id >= 1 && element_id <= 32) + b = sin(4 * pi * x[1]) + 3 + H = max(b, 2.5) + end + + # Shift the water level by the amount `equations.threshold_limiter` (1e-13 per default) + # to avoid division by a hard 0 value in the water height `h` when computing velocities. + H = max(H, b + equations.threshold_limiter) + return prim2cons(SVector(H, v, b), equations) +end + +# point to the data we want to augment +u = Trixi.wrap_array(ode.u0, semi) +# reset the initial condition +for element in eachelement(semi.solver, semi.cache) + for i in eachnode(semi.solver) + x_node = Trixi.get_node_coords(semi.cache.elements.node_coordinates, equations, semi.solver, i, element) + u_node = initial_condition_discontinuous_well_balancedness(x_node, first(tspan), element, equations) + Trixi.set_node_vars!(u, u_node, equations, semi.solver, i, element) + end +end summary_callback = SummaryCallback() analysis_interval = 1000 -analysis_callback = AnalysisCallback(semi, interval=analysis_interval, save_analysis=false, - extra_analysis_integrals=(energy_kinetic, - energy_internal, - lake_at_rest_error)) +analysis_callback = AnalysisCallback(semi, interval=analysis_interval, save_analysis=false) alive_callback = AliveCallback(analysis_interval=analysis_interval) @@ -98,7 +122,7 @@ save_solution = SaveSolutionCallback(interval=5000, save_initial_solution=true, save_final_solution=true) -stepsize_callback = StepsizeCallback(cfl=2.0) +stepsize_callback = StepsizeCallback(cfl=1.5) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, stepsize_callback) @@ -112,4 +136,50 @@ stage_limiter! = PositivityPreservingLimiterShallowWater(thresholds=(equations.t sol = solve(ode, SSPRK43(stage_limiter!), dt=1.0, save_everystep=false, callback=callbacks, adaptive=false); -summary_callback() # print the timer summary \ No newline at end of file +summary_callback() # print the timer summary + +############################################################################### +# workaround to compute the well-balancedness error for this partiular problem +# that has two reference water heights. One for a lake to the left of the +# discontinuous bottom topography `H0_upper = 2.5` and another for a lake to the +# right of the discontinuous bottom topography `H0_lower = 1.5`. + +# Declare a special version of the function to compute the lake-at-rest error +# OBS! The reference water height values are hardcoded for convenience. +function lake_at_rest_error_two_level(u, element_id, equations::ShallowWaterEquations1D) + h, _, b = u + + # 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. + + # element_id check for refinement_level=6 + if (element_id >= 1 && element_id <= 32) + H0_wet_dry = max( 2.5 , b + equations.threshold_limiter ) + else + H0_wet_dry = max( 1.5 , b + equations.threshold_limiter ) + end + + return abs(H0_wet_dry - (h + b)) +end + +# point to the data we want to analyze +u = Trixi.wrap_array(sol[end], semi) +# Perform the actual integration of the well-balancedness error over the domain +l1_well_balance_error = Trixi.integrate_via_indices(u, mesh, equations, semi.solver, semi.cache; normalize=true) do u, i, element, equations, solver + u_local = Trixi.get_node_vars(u, equations, solver, i, element) + return lake_at_rest_error_two_level(u_local, element, equations) +end + +# report the well-balancedness lake-at-rest error to the screen +if Trixi.mpi_isroot() + Trixi.mpi_println("─"^100) + Trixi.mpi_println(" Lake-at-rest error for '", Trixi.get_name(equations), "' with ", Trixi.summary(solver), + " at final time " * Trixi.@sprintf("%10.8e", tspan[end])) + + Trixi.@printf(" %-12s:", Trixi.pretty_form_utf(lake_at_rest_error)) + Trixi.@printf(" % 10.8e", l1_well_balance_error) + Trixi.mpi_println() + Trixi.mpi_println("─"^100) +end \ No newline at end of file diff --git a/src/callbacks_step/analysis.jl b/src/callbacks_step/analysis.jl index e29dac278d..79a442f06c 100644 --- a/src/callbacks_step/analysis.jl +++ b/src/callbacks_step/analysis.jl @@ -548,17 +548,6 @@ end pretty_form_utf(quantity) = get_name(quantity) pretty_form_ascii(quantity) = get_name(quantity) -# Overload analyze and pass the initial_condition as parameter for integrate. -# Special case only used for lake_at_rest_error for AbstractShallowWaterEquations. -# With the new wet-dry mechanics, the old calculation using a constant H0 as a reference value fails. -# For example, when a hill extends to be higher than the waterlevel. -# New approach: Compare waterheight h at time T with the initial waterheight -function analyze(::typeof(lake_at_rest_error), du, u, t, semi::AbstractSemidiscretization) - mesh, equations, solver, cache = mesh_equations_solver_cache(semi) - @unpack initial_condition = semi - - integrate(lake_at_rest_error, u, mesh, equations, solver, cache, initial_condition, normalize=true) -end # Special analyze for `SemidiscretizationHyperbolicParabolic` such that # precomputed gradients are available. For now only implemented for the `enstrophy` diff --git a/src/callbacks_step/analysis_dg1d.jl b/src/callbacks_step/analysis_dg1d.jl index 3e3f08a107..e92701dc1f 100644 --- a/src/callbacks_step/analysis_dg1d.jl +++ b/src/callbacks_step/analysis_dg1d.jl @@ -175,20 +175,6 @@ function integrate(func::Func, u, end end -# Overload to pass the initial_condition as parameter for calculating lake_at_rest_error -function integrate(::typeof(lake_at_rest_error), u, - mesh::Union{TreeMesh{1}, StructuredMesh{1}}, - equations::ShallowWaterEquations1D, dg::DGSEM, cache, - initial_condition; normalize=true) - node_coordinates = cache.elements.node_coordinates - integrate_via_indices(u, mesh, equations, dg, cache; normalize=normalize) do u, i, element, equations, dg - x = get_node_coords(node_coordinates, equations, dg, i, element) - - u_exact = initial_condition(x, 0.0, equations) - u_local = get_node_vars(u, equations, dg, i, element) - return lake_at_rest_error(u_local, u_exact, equations) - end -end function analyze(::typeof(entropy_timederivative), du, u, t, mesh::Union{TreeMesh{1},StructuredMesh{1}}, equations, dg::DG, cache) diff --git a/src/callbacks_step/analysis_dg2d.jl b/src/callbacks_step/analysis_dg2d.jl index 45b90e6855..27d7057c29 100644 --- a/src/callbacks_step/analysis_dg2d.jl +++ b/src/callbacks_step/analysis_dg2d.jl @@ -194,20 +194,6 @@ function integrate(func::Func, u, end end -#Overload to pass the initial_condition as parameter for calculating lake_at_rest_error - function integrate(::typeof(lake_at_rest_error), u, - mesh::Union{TreeMesh{2}, StructuredMesh{2}, UnstructuredMesh2D, P4estMesh{2}}, - equations::ShallowWaterEquations2D, dg::DG, cache, initial_condition; normalize=true) - node_coordinates = cache.elements.node_coordinates - integrate_via_indices(u, mesh, equations, dg, cache; normalize=normalize) do u, i, j, element, equations, dg - x = get_node_coords(node_coordinates, equations, dg, i, j, element) - - u_exact = initial_condition(x, 0.0, equations) - u_local = get_node_vars(u, equations, dg, i, j, element) - return lake_at_rest_error(u_local, u_exact, equations) - end - end - function analyze(::typeof(entropy_timederivative), du, u, t, mesh::Union{TreeMesh{2}, StructuredMesh{2}, UnstructuredMesh2D, P4estMesh{2}}, equations, dg::DG, cache) diff --git a/src/equations/shallow_water_1d.jl b/src/equations/shallow_water_1d.jl index 05c593bd16..24757e9172 100644 --- a/src/equations/shallow_water_1d.jl +++ b/src/equations/shallow_water_1d.jl @@ -54,7 +54,7 @@ struct ShallowWaterEquations1D{RealT<:Real} <: AbstractShallowWaterEquations{1, gravity::RealT # gravitational constant H0::RealT # constant "lake-at-rest" total water height threshold_limiter::RealT # Threshold to use in PositivityPreservingLimiterShallowWater on water height, - # as a (small) shift on the initial condition and cutoff before the + # as a (small) shift on the initial condition and cutoff before the # next time step. threshold_wet::RealT # Threshold to be applied on water height to define when the flow is "wet" # before calculating the numerical flux. @@ -331,7 +331,7 @@ Should be used together with [`FluxHydrostaticReconstruction`](@ref) and [`hydrostatic_reconstruction_chen_noelle`](@ref) in the surface flux to ensure consistency. Further details on the hydrostatic reconstruction and its motivation can be found in -- Guoxian Chen and Sebastian Noelle (2017) +- Guoxian Chen and Sebastian Noelle (2017) A new hydrostatic reconstruction scheme based on subcell reconstructions [DOI:10.1137/15M1053074](https://dx.doi.org/10.1137/15M1053074) """ @@ -477,7 +477,7 @@ The key idea is a linear reconstruction of the bottom and water height at the in Use in combination with the generic numerical flux routine [`FluxHydrostaticReconstruction`](@ref). Further details on this hydrostatic reconstruction and its motivation can be found in -- Guoxian Chen and Sebastian Noelle (2017) +- Guoxian Chen and Sebastian Noelle (2017) A new hydrostatic reconstruction scheme based on subcell reconstructions [DOI:10.1137/15M1053074](https://dx.doi.org/10.1137/15M1053074) """ @@ -586,12 +586,12 @@ end min_max_speed_chen_noelle(u_ll, u_rr, orientation::Integer, equations::ShallowWaterEquations1D) -The approximated speeds for the HLL type numerical flux used by Chen and Noelle for their +The approximated speeds for the HLL type numerical flux used by Chen and Noelle for their hydrostatic reconstruction. As they state in the paper, these speeds are chosen for the numerical flux to ensure positivity and to satisfy an entropy inequality. Further details on this hydrostatic reconstruction and its motivation can be found in -- Guoxian Chen and Sebastian Noelle (2017) +- Guoxian Chen and Sebastian Noelle (2017) A new hydrostatic reconstruction scheme based on subcell reconstructions [DOI:10.1137/15M1053074](https://dx.doi.org/10.1137/15M1053074) """ @@ -607,7 +607,7 @@ Further details on this hydrostatic reconstruction and its motivation can be fou a_ll = sqrt(equations.gravity * h_ll) a_rr = sqrt(equations.gravity * h_rr) - λ_min = min( v_ll - a_ll, v_rr - a_rr, zero(eltype(u_ll)) ) + λ_min = min( v_ll - a_ll, v_rr - a_rr, zero(eltype(u_ll)) ) λ_max = max( v_ll + a_ll, v_rr + a_rr, zero(eltype(u_ll)) ) return λ_min, λ_max @@ -723,15 +723,18 @@ end # Calculate the error for the "lake-at-rest" test case where H = h+b should -# be a constant value over time -@inline function lake_at_rest_error(u, u_exact, equations::ShallowWaterEquations1D) +# be a constant value over time. Note, assumes there is a single reference +# water height `H0` with which to compare. +@inline function lake_at_rest_error(u, equations::ShallowWaterEquations1D) h, _, b = u - h_exact, _, b_exact= u_exact - H = h + b - H_exact = h_exact + b_exact + # 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(H - H_exact) + return abs(H0_wet_dry - (h + b)) end end # @muladd diff --git a/src/equations/shallow_water_2d.jl b/src/equations/shallow_water_2d.jl index b339fd37f3..c9c5b34af6 100644 --- a/src/equations/shallow_water_2d.jl +++ b/src/equations/shallow_water_2d.jl @@ -874,21 +874,18 @@ end # Calculate the error for the "lake-at-rest" test case where H = h+b should -# be a constant value over time +# be a constant value over time. Note, assumes there is a single reference +# water height `H0` with which to compare. @inline function lake_at_rest_error(u, equations::ShallowWaterEquations2D) h, _, _, b = u - return abs(equations.H0 - (h + b)) -end - -@inline function lake_at_rest_error(u, u_exact, equations::ShallowWaterEquations2D) - h, _, _, b = u - - h_exact, _, _, b_exact= u_exact - H = h + b - H_exact = h_exact + b_exact + # + # TODO: Compute this the proper way with + # H0_wet_dry = max( equations.H0 , b + equations.threshold_limiter ) + # once PR for wet/dry `ShallowWaterEquations2D` is merged into this PR + # - return abs(H - H_exact) + return abs(equations.H0 - (h + b)) end end # @muladd diff --git a/test/test_tree_1d_shallowwater.jl b/test/test_tree_1d_shallowwater.jl index 2bb7ae469b..32d48dd7b5 100644 --- a/test/test_tree_1d_shallowwater.jl +++ b/test/test_tree_1d_shallowwater.jl @@ -30,10 +30,10 @@ EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "tree_1 tspan = (0.0, 0.25)) end - @trixi_testset "elixir_shallowwater_well_balanced_wet_dry.jl" begin + @trixi_testset "elixir_shallowwater_well_balanced_wet_dry.jl with FluxHydrostaticReconstruction" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_well_balanced_wet_dry.jl"), - l2 = [0.027869290160397523, 4.848702566024514e-14, 0.10911645796834736], - linf = [0.28423245872819947, 2.305820093827649e-13, 1.119394067096447], + l2 = [0.009657871671690443, 5.403349017029091e-15, 0.03857583749209983], + linf = [0.49999999999990025, 3.344735163283149e-14, 1.9999999999999998], tspan = (0.0, 0.25)) end From 5c1ac3855d9fbd4bcd7b300e131d666cf0e65f65 Mon Sep 17 00:00:00 2001 From: svengoldberg Date: Thu, 4 May 2023 14:38:48 +0200 Subject: [PATCH 22/44] Fix typo --- src/callbacks_stage/positivity_shallow_water_dg1d.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/callbacks_stage/positivity_shallow_water_dg1d.jl b/src/callbacks_stage/positivity_shallow_water_dg1d.jl index a71f4d0ff0..d444df88f9 100644 --- a/src/callbacks_stage/positivity_shallow_water_dg1d.jl +++ b/src/callbacks_stage/positivity_shallow_water_dg1d.jl @@ -50,7 +50,7 @@ function limiter_shallow_water!(u, threshold::Real, variable, u_node = SVector(h_node, h_v_node, b_node) u_mean = SVector(h_mean, h_v_mean, zero(eltype(u))) - # When velocity is cut off, the only value beeing averaged is the waterheight, + # When velocity is cut off, the only averaged value is the waterheight, # because the velocity is set to zero and this value is passed. # Otherwise, the velocity is averaged, as well. # Note that the auxiliary bottom topography variable `b` is never limited. From cdc86c14d3f333823b89681dae82a4998c0cf574 Mon Sep 17 00:00:00 2001 From: svengoldberg Date: Thu, 4 May 2023 15:18:44 +0200 Subject: [PATCH 23/44] Shorten test run for parabolic bowl 1D --- test/test_tree_1d_shallowwater.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/test_tree_1d_shallowwater.jl b/test/test_tree_1d_shallowwater.jl index 32d48dd7b5..09c7a0b6d1 100644 --- a/test/test_tree_1d_shallowwater.jl +++ b/test/test_tree_1d_shallowwater.jl @@ -97,9 +97,9 @@ EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "tree_1 @trixi_testset "elixir_shallowwater_parabolic_bowl.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_parabolic_bowl.jl"), - l2 = [0.012868412255957612, 0.007819719632793341, 3.8569975327636153e-17], - linf = [0.01898697414997512, 0.011686441632120272, 2.7755575615628914e-16], - tspan = (0.0, 0.75)) + l2 = [0.0020698668618855105, 0.0004971872608360286, 3.8569975327636153e-17], + linf = [0.0030589243995025744, 0.000840091286020108, 2.7755575615628914e-16], + tspan = (0.0, 0.25)) end end From af9f8cb467005491abc518d2c94b5156cf7d702e Mon Sep 17 00:00:00 2001 From: svengoldberg Date: Thu, 4 May 2023 23:01:06 +0200 Subject: [PATCH 24/44] Choose lower resolution for parabolic bowl and update test values --- examples/tree_1d_dgsem/elixir_shallowwater_parabolic_bowl.jl | 2 +- test/test_tree_1d_shallowwater.jl | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_parabolic_bowl.jl b/examples/tree_1d_dgsem/elixir_shallowwater_parabolic_bowl.jl index dedfad1586..80855c00b0 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_parabolic_bowl.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_parabolic_bowl.jl @@ -70,7 +70,7 @@ coordinates_min = -2.0 coordinates_max = 2.0 mesh = TreeMesh(coordinates_min, coordinates_max, - initial_refinement_level=10, + initial_refinement_level=7, n_cells_max=10_000) # create the semi discretization object diff --git a/test/test_tree_1d_shallowwater.jl b/test/test_tree_1d_shallowwater.jl index 09c7a0b6d1..11d55ef4c4 100644 --- a/test/test_tree_1d_shallowwater.jl +++ b/test/test_tree_1d_shallowwater.jl @@ -97,8 +97,8 @@ EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "tree_1 @trixi_testset "elixir_shallowwater_parabolic_bowl.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_parabolic_bowl.jl"), - l2 = [0.0020698668618855105, 0.0004971872608360286, 3.8569975327636153e-17], - linf = [0.0030589243995025744, 0.000840091286020108, 2.7755575615628914e-16], + l2 = [0.0020644908545138043, 0.0004917516092096569, 3.9831502760522817e-17], + linf = [0.003726955238821753, 0.0007740912880794988, 2.220446049250313e-16], tspan = (0.0, 0.25)) end end From d8a5bdb749849cec327d7e2772c06c439d6db1b8 Mon Sep 17 00:00:00 2001 From: svengoldberg Date: Fri, 5 May 2023 09:47:54 +0200 Subject: [PATCH 25/44] Further reduce resolution for parabolic bowl and update test values --- examples/tree_1d_dgsem/elixir_shallowwater_parabolic_bowl.jl | 2 +- test/test_tree_1d_shallowwater.jl | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_parabolic_bowl.jl b/examples/tree_1d_dgsem/elixir_shallowwater_parabolic_bowl.jl index 80855c00b0..18b33f055f 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_parabolic_bowl.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_parabolic_bowl.jl @@ -70,7 +70,7 @@ coordinates_min = -2.0 coordinates_max = 2.0 mesh = TreeMesh(coordinates_min, coordinates_max, - initial_refinement_level=7, + initial_refinement_level=6, n_cells_max=10_000) # create the semi discretization object diff --git a/test/test_tree_1d_shallowwater.jl b/test/test_tree_1d_shallowwater.jl index 11d55ef4c4..86441d9668 100644 --- a/test/test_tree_1d_shallowwater.jl +++ b/test/test_tree_1d_shallowwater.jl @@ -97,8 +97,8 @@ EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "tree_1 @trixi_testset "elixir_shallowwater_parabolic_bowl.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_parabolic_bowl.jl"), - l2 = [0.0020644908545138043, 0.0004917516092096569, 3.9831502760522817e-17], - linf = [0.003726955238821753, 0.0007740912880794988, 2.220446049250313e-16], + l2 = [0.002062937061879629, 0.00048160584495176244, 3.828565374146696e-17], + linf = [0.003858884888644061, 0.000799010577904765, 1.6653345369377348e-16], tspan = (0.0, 0.25)) end end From 7595ac2506d3b533133f44e55936460ea230f3a4 Mon Sep 17 00:00:00 2001 From: svengoldberg Date: Wed, 10 May 2023 12:06:30 +0200 Subject: [PATCH 26/44] Remove MPI from well-balanced test --- ...lixir_shallowwater_well_balanced_wet_dry.jl | 18 ++++++++---------- 1 file changed, 8 insertions(+), 10 deletions(-) diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl b/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl index 750c092291..f5f2564f0a 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl @@ -173,13 +173,11 @@ l1_well_balance_error = Trixi.integrate_via_indices(u, mesh, equations, semi.sol end # report the well-balancedness lake-at-rest error to the screen -if Trixi.mpi_isroot() - Trixi.mpi_println("─"^100) - Trixi.mpi_println(" Lake-at-rest error for '", Trixi.get_name(equations), "' with ", Trixi.summary(solver), - " at final time " * Trixi.@sprintf("%10.8e", tspan[end])) - - Trixi.@printf(" %-12s:", Trixi.pretty_form_utf(lake_at_rest_error)) - Trixi.@printf(" % 10.8e", l1_well_balance_error) - Trixi.mpi_println() - Trixi.mpi_println("─"^100) -end \ No newline at end of file +println("─"^100) +println(" Lake-at-rest error for '", Trixi.get_name(equations), "' with ", Trixi.summary(solver), + " at final time " * Trixi.@sprintf("%10.8e", tspan[end])) + +Trixi.@printf(" %-12s:", Trixi.pretty_form_utf(lake_at_rest_error)) +Trixi.@printf(" % 10.8e", l1_well_balance_error) +println() +println("─"^100) \ No newline at end of file From ae55a5c20d7059ca36f4cd287a9d7f6c75aa2517 Mon Sep 17 00:00:00 2001 From: svengoldberg Date: Fri, 12 May 2023 20:28:01 +0200 Subject: [PATCH 27/44] Simplify workaround to set discontinuity --- ...ixir_shallowwater_well_balanced_wet_dry.jl | 83 ++++++++++--------- 1 file changed, 45 insertions(+), 38 deletions(-) diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl b/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl index f5f2564f0a..8de4f05722 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl @@ -8,7 +8,7 @@ using Trixi equations = ShallowWaterEquations1D(gravity_constant=9.812) """ - initial_condition_well_balanced_chen_noelle(x, t, equations:: ShallowWaterEquations1D) +initial_condition_complex_bottom_well_balanced(x, t, equations:: ShallowWaterEquations1D) Initial condition with a complex (discontinuous) bottom topography to test the well-balanced property for the [`hydrostatic_reconstruction_chen_noelle`](@ref) including dry areas within the @@ -19,7 +19,7 @@ The initial condition was found in the section 5.2 of the paper: A new hydrostatic reconstruction scheme based on subcell reconstructions [DOI:10.1137/15M1053074](https://dx.doi.org/10.1137/15M1053074) """ -function initial_condition_well_balanced_chen_noelle(x, t, equations:: ShallowWaterEquations1D) +function initial_condition_complex_bottom_well_balanced(x, t, equations:: ShallowWaterEquations1D) v = 0.0 b = sin(4 * pi * x[1]) + 3 @@ -42,7 +42,7 @@ function initial_condition_well_balanced_chen_noelle(x, t, equations:: ShallowWa return prim2cons(SVector(H, v, b), equations) end -initial_condition = initial_condition_well_balanced_chen_noelle +initial_condition = initial_condition_complex_bottom_well_balanced ############################################################################### # Get the DG approximation space @@ -83,22 +83,16 @@ semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) tspan = (0.0, 25.0) ode = semidiscretize(semi, tspan) -function initial_condition_discontinuous_well_balancedness(x, t, element_id, equations:: ShallowWaterEquations1D) - v = 0.0 - b = sin(4 * pi * x[1]) + 1 - H = max(b, 1.5) - - # for refinement_level=6 - if (element_id >= 1 && element_id <= 32) - b = sin(4 * pi * x[1]) + 3 - H = max(b, 2.5) - end - # Shift the water level by the amount `equations.threshold_limiter` (1e-13 per default) - # to avoid division by a hard 0 value in the water height `h` when computing velocities. - H = max(H, b + equations.threshold_limiter) - return prim2cons(SVector(H, v, b), equations) -end +############################################################################### +# Workaround to set a discontinuous water and bottom topography for +# debugging and testing. Essentially, this is a slight augmentation of the +# `compute_coefficients` where the `x` node value passed here is slightly +# perturbed to the left / right in order to set a true discontinuity that avoids +# the doubled value of the LGL nodes at a particular element interface. +# +# Note! The errors from the analysis callback are not important but the error +# for this lake at rest test case `∑|H0-(h+b)|` should be near machine roundoff. # point to the data we want to augment u = Trixi.wrap_array(ode.u0, semi) @@ -106,7 +100,14 @@ u = Trixi.wrap_array(ode.u0, semi) for element in eachelement(semi.solver, semi.cache) for i in eachnode(semi.solver) x_node = Trixi.get_node_coords(semi.cache.elements.node_coordinates, equations, semi.solver, i, element) - u_node = initial_condition_discontinuous_well_balancedness(x_node, first(tspan), element, equations) + # We know that the discontinuity is on an interface. Slightly augment the x value by a factor + # of unit roundoff to avoid the repeated value from the LGL nodes at the interface. + if i == 1 + x_node = SVector(nextfloat(x_node[1], 1)) + elseif i == nnodes(semi.solver) + x_node = SVector(prevfloat(x_node[1], 1)) + end + u_node = initial_condition_complex_bottom_well_balanced(x_node, first(tspan), equations) Trixi.set_node_vars!(u, u_node, equations, semi.solver, i, element) end end @@ -146,30 +147,36 @@ summary_callback() # print the timer summary # Declare a special version of the function to compute the lake-at-rest error # OBS! The reference water height values are hardcoded for convenience. -function lake_at_rest_error_two_level(u, element_id, equations::ShallowWaterEquations1D) - h, _, b = u - - # 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. - - # element_id check for refinement_level=6 - if (element_id >= 1 && element_id <= 32) - H0_wet_dry = max( 2.5 , b + equations.threshold_limiter ) - else - H0_wet_dry = max( 1.5 , b + equations.threshold_limiter ) - end - - return abs(H0_wet_dry - (h + b)) -end +function lake_at_rest_error_two_level(u, x, equations::ShallowWaterEquations1D) + h, _, b = u + + # 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. + if x[1] < 0.5 + H0_wet_dry = max( 2.5 , b + equations.threshold_limiter ) + else + H0_wet_dry = max( 1.5 , b + equations.threshold_limiter ) + end + + return abs(H0_wet_dry - (h + b)) + end # point to the data we want to analyze u = Trixi.wrap_array(sol[end], semi) # Perform the actual integration of the well-balancedness error over the domain l1_well_balance_error = Trixi.integrate_via_indices(u, mesh, equations, semi.solver, semi.cache; normalize=true) do u, i, element, equations, solver - u_local = Trixi.get_node_vars(u, equations, solver, i, element) - return lake_at_rest_error_two_level(u_local, element, equations) + x_node = Trixi.get_node_coords(semi.cache.elements.node_coordinates, equations, solver, i, element) + # We know that the discontinuity is a vertical line. Slightly augment the x value by a factor + # of unit roundoff to avoid the repeted value from the LGL nodes at at interface. + if i == 1 + x_node = SVector(nextfloat(x_node[1], 1)) + elseif i == nnodes(semi.solver) + x_node = SVector(prevfloat(x_node[1], 1)) + end + u_local = Trixi.get_node_vars(u, equations, solver, i, element) + return lake_at_rest_error_two_level(u_local, x_node, equations) end # report the well-balancedness lake-at-rest error to the screen From ca4ca18a6473330dbaef5be3d185ba1a25938309 Mon Sep 17 00:00:00 2001 From: svengoldberg Date: Thu, 18 May 2023 14:43:53 +0200 Subject: [PATCH 28/44] Change structure of Chen&Noelle flux --- src/Trixi.jl | 3 +-- src/equations/numerical_fluxes.jl | 3 +-- 2 files changed, 2 insertions(+), 4 deletions(-) diff --git a/src/Trixi.jl b/src/Trixi.jl index afce259cf5..9f20d1b942 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -153,8 +153,7 @@ export flux, flux_central, flux_lax_friedrichs, flux_hll, flux_hllc, flux_hlle, hydrostatic_reconstruction_chen_noelle, flux_nonconservative_chen_noelle, flux_hll_chen_noelle, FluxPlusDissipation, DissipationGlobalLaxFriedrichs, DissipationLocalLaxFriedrichs, FluxLaxFriedrichs, max_abs_speed_naive, - FluxHLL, min_max_speed_naive, - FluxHLLChenNoelle, min_max_speed_chen_noelle, + FluxHLL, min_max_speed_naive, min_max_speed_chen_noelle, FluxLMARS, FluxRotated, flux_shima_etal_turbo, flux_ranocha_turbo, diff --git a/src/equations/numerical_fluxes.jl b/src/equations/numerical_fluxes.jl index 8e9bbb4487..8a0a91bb9b 100644 --- a/src/equations/numerical_fluxes.jl +++ b/src/equations/numerical_fluxes.jl @@ -203,7 +203,6 @@ struct FluxHLL{MinMaxSpeed} end FluxHLL() = FluxHLL(min_max_speed_naive) -FluxHLLChenNoelle() = FluxHLL(min_max_speed_chen_noelle) """ min_max_speed_naive(u_ll, u_rr, orientation::Integer, equations) @@ -245,7 +244,7 @@ Base.show(io::IO, numflux::FluxHLL) = print(io, "FluxHLL(", numflux.min_max_spee See [`FluxHLL`](@ref). """ const flux_hll = FluxHLL() -const flux_hll_chen_noelle = FluxHLLChenNoelle() +const flux_hll_chen_noelle = FluxHLL(min_max_speed_chen_noelle) """ From 1a734afefd33957157002062e02004fe9f8d8206 Mon Sep 17 00:00:00 2001 From: svengoldberg Date: Thu, 18 May 2023 14:48:27 +0200 Subject: [PATCH 29/44] Fix typos and indenting --- examples/tree_1d_dgsem/elixir_shallowwater_beach.jl | 8 ++++---- .../elixir_shallowwater_parabolic_bowl.jl | 8 ++++---- .../elixir_shallowwater_well_balanced_wet_dry.jl | 10 +++++----- src/equations/shallow_water_1d.jl | 4 ++-- 4 files changed, 15 insertions(+), 15 deletions(-) diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_beach.jl b/examples/tree_1d_dgsem/elixir_shallowwater_beach.jl index 0110b2817f..d6bdd7da01 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_beach.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_beach.jl @@ -94,9 +94,9 @@ summary_callback = SummaryCallback() analysis_interval = 1000 analysis_callback = AnalysisCallback(semi, interval=analysis_interval, save_analysis=false, - extra_analysis_integrals=(energy_kinetic, - energy_internal, - lake_at_rest_error)) + extra_analysis_integrals=(energy_kinetic, + energy_internal, + lake_at_rest_error)) alive_callback = AliveCallback(analysis_interval=analysis_interval) @@ -107,7 +107,7 @@ save_solution = SaveSolutionCallback(interval=250, callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution) stage_limiter! = PositivityPreservingLimiterShallowWater(thresholds=(equations.threshold_limiter,), - variables=(Trixi.waterheight,)) + variables=(Trixi.waterheight,)) ############################################################################### # run the simulation diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_parabolic_bowl.jl b/examples/tree_1d_dgsem/elixir_shallowwater_parabolic_bowl.jl index 18b33f055f..1846a957b1 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_parabolic_bowl.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_parabolic_bowl.jl @@ -86,9 +86,9 @@ summary_callback = SummaryCallback() analysis_interval = 1000 analysis_callback = AnalysisCallback(semi, interval=analysis_interval, save_analysis=false, - extra_analysis_integrals=(energy_kinetic, - energy_internal, - lake_at_rest_error)) + extra_analysis_integrals=(energy_kinetic, + energy_internal, + lake_at_rest_error)) alive_callback = AliveCallback(analysis_interval=analysis_interval) @@ -99,7 +99,7 @@ save_solution = SaveSolutionCallback(interval=1000, callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution) stage_limiter! = PositivityPreservingLimiterShallowWater(thresholds=(equations.threshold_limiter,), - variables=(Trixi.waterheight,)) + variables=(Trixi.waterheight,)) ############################################################################### # run the simulation diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl b/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl index 8de4f05722..9ffa5f954f 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl @@ -103,9 +103,9 @@ for element in eachelement(semi.solver, semi.cache) # We know that the discontinuity is on an interface. Slightly augment the x value by a factor # of unit roundoff to avoid the repeated value from the LGL nodes at the interface. if i == 1 - x_node = SVector(nextfloat(x_node[1], 1)) + x_node = SVector(nextfloat(x_node[1])) elseif i == nnodes(semi.solver) - x_node = SVector(prevfloat(x_node[1], 1)) + x_node = SVector(prevfloat(x_node[1])) end u_node = initial_condition_complex_bottom_well_balanced(x_node, first(tspan), equations) Trixi.set_node_vars!(u, u_node, equations, semi.solver, i, element) @@ -129,7 +129,7 @@ callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, sav stepsize_callback) stage_limiter! = PositivityPreservingLimiterShallowWater(thresholds=(equations.threshold_limiter,), - variables=(Trixi.waterheight,)) + variables=(Trixi.waterheight,)) ############################################################################### # run the simulation @@ -171,9 +171,9 @@ l1_well_balance_error = Trixi.integrate_via_indices(u, mesh, equations, semi.sol # We know that the discontinuity is a vertical line. Slightly augment the x value by a factor # of unit roundoff to avoid the repeted value from the LGL nodes at at interface. if i == 1 - x_node = SVector(nextfloat(x_node[1], 1)) + x_node = SVector(nextfloat(x_node[1])) elseif i == nnodes(semi.solver) - x_node = SVector(prevfloat(x_node[1], 1)) + x_node = SVector(prevfloat(x_node[1])) end u_local = Trixi.get_node_vars(u, equations, solver, i, element) return lake_at_rest_error_two_level(u_local, x_node, equations) diff --git a/src/equations/shallow_water_1d.jl b/src/equations/shallow_water_1d.jl index 24757e9172..9f68ffb654 100644 --- a/src/equations/shallow_water_1d.jl +++ b/src/equations/shallow_water_1d.jl @@ -470,8 +470,8 @@ end hydrostatic_reconstruction_chen_noelle(u_ll, u_rr, orientation::Integer, equations::ShallowWaterEquations1D) -A particular type of hydrostatic reconstruction on the water height to guarantee well-balancedness -for a general bottom topography [`ShallowWaterEquations1D`](@ref). The reconstructed solution states +A particular type of hydrostatic reconstruction of the water height to guarantee well-balancedness +for a general bottom topography of the [`ShallowWaterEquations1D`](@ref). The reconstructed solution states `u_ll_star` and `u_rr_star` variables are used to evaluate the surface numerical flux at the interface. The key idea is a linear reconstruction of the bottom and water height at the interfaces using subcells. Use in combination with the generic numerical flux routine [`FluxHydrostaticReconstruction`](@ref). From b1b9fb16219bd4ec1f52ba406101f8887757c496 Mon Sep 17 00:00:00 2001 From: svengoldberg Date: Thu, 18 May 2023 14:54:58 +0200 Subject: [PATCH 30/44] Adjust call of solve and use ode_default_options --- examples/tree_1d_dgsem/elixir_shallowwater_beach.jl | 4 ++-- .../tree_1d_dgsem/elixir_shallowwater_parabolic_bowl.jl | 4 ++-- .../elixir_shallowwater_well_balanced_wet_dry.jl | 4 ++-- test/test_tree_1d_shallowwater.jl | 8 ++++---- 4 files changed, 10 insertions(+), 10 deletions(-) diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_beach.jl b/examples/tree_1d_dgsem/elixir_shallowwater_beach.jl index d6bdd7da01..e7c8a36c29 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_beach.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_beach.jl @@ -112,7 +112,7 @@ stage_limiter! = PositivityPreservingLimiterShallowWater(thresholds=(equations.t ############################################################################### # run the simulation -sol = solve(ode, SSPRK43(stage_limiter!), dt=1.0, - save_everystep=false, callback=callbacks); +sol = solve(ode, SSPRK43(stage_limiter!); + ode_default_options()..., callback=callbacks); summary_callback() # print the timer summary \ No newline at end of file diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_parabolic_bowl.jl b/examples/tree_1d_dgsem/elixir_shallowwater_parabolic_bowl.jl index 1846a957b1..7f43346e0a 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_parabolic_bowl.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_parabolic_bowl.jl @@ -104,7 +104,7 @@ stage_limiter! = PositivityPreservingLimiterShallowWater(thresholds=(equations.t ############################################################################### # run the simulation -sol = solve(ode, SSPRK43(stage_limiter!), dt=1.0, - save_everystep=false, callback=callbacks); +sol = solve(ode, SSPRK43(stage_limiter!); + ode_default_options()..., callback=callbacks); summary_callback() # print the timer summary \ No newline at end of file diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl b/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl index 9ffa5f954f..e561091ca4 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl @@ -134,8 +134,8 @@ stage_limiter! = PositivityPreservingLimiterShallowWater(thresholds=(equations.t ############################################################################### # run the simulation -sol = solve(ode, SSPRK43(stage_limiter!), dt=1.0, - save_everystep=false, callback=callbacks, adaptive=false); +sol = solve(ode, SSPRK43(stage_limiter!); dt=1.0, + ode_default_options()..., callback=callbacks, adaptive=false); summary_callback() # print the timer summary diff --git a/test/test_tree_1d_shallowwater.jl b/test/test_tree_1d_shallowwater.jl index 86441d9668..09f9b338dd 100644 --- a/test/test_tree_1d_shallowwater.jl +++ b/test/test_tree_1d_shallowwater.jl @@ -90,15 +90,15 @@ EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "tree_1 @trixi_testset "elixir_shallowwater_beach.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_beach.jl"), - l2 = [0.5837839017922798, 3.7405541768121595, 6.289710785696625e-8], - linf = [0.9917266589187276, 7.096283536176454, 4.452604027704865e-7], + l2 = [0.5837835144970676, 3.740558401311243, 6.289710785696625e-8], + linf = [0.9917277949257116, 7.096293347373189, 4.452604027704865e-7], tspan = (0.0, 0.25)) end @trixi_testset "elixir_shallowwater_parabolic_bowl.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_parabolic_bowl.jl"), - l2 = [0.002062937061879629, 0.00048160584495176244, 3.828565374146696e-17], - linf = [0.003858884888644061, 0.000799010577904765, 1.6653345369377348e-16], + l2 = [0.0020629352664007466, 0.0004816398224727313, 3.828565374146696e-17], + linf = [0.0038611903839202823, 0.0008591988088082134, 1.6653345369377348e-16], tspan = (0.0, 0.25)) end end From f579a72f4dfb5b3a9a422242433fa3c6157a2a4d Mon Sep 17 00:00:00 2001 From: svengoldberg Date: Thu, 18 May 2023 15:02:49 +0200 Subject: [PATCH 31/44] Edit docstring --- .../positivity_shallow_water.jl | 146 +++++++++--------- 1 file changed, 72 insertions(+), 74 deletions(-) diff --git a/src/callbacks_stage/positivity_shallow_water.jl b/src/callbacks_stage/positivity_shallow_water.jl index 2a9a6a60e8..d1e2792659 100644 --- a/src/callbacks_stage/positivity_shallow_water.jl +++ b/src/callbacks_stage/positivity_shallow_water.jl @@ -5,81 +5,79 @@ @muladd begin - """ - PositivityPreservingLimiterShallowWater(; threshold, variables) - - 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) - - It is specifically designed for `ShallowWaterEquations`. - The limiter is applied to all scalar `variables` in their given order - using the associated `thresholds` to determine the minimal acceptable values. - The order of the `variables` is important and might have a strong influence - on the robustness. - As opposed to the standard version of the [`PositivityPreservingLimiterZhangShu`](@ref), - nodes with a water height below the `threshold_limiter` are treated in a special way. - To avoid numerical problems caused by velocities close to zero, - the velocity is cut off, such that the node can be identified as "dry". The special feature of the - `ShallowWaterEquations` used here is that the bottom topography is stored as an additional - quantity in the solution vector `u`. However, the value of the bottom topography - should not be changed, much less limited. That is why, the `set_node_vars!` function is not applied - to the last conservation variable. - After the limiting process is applied to all degrees of freedom, for safety reasons, - the wet/dry threshold is applied again on all the DG nodes in order to avoid dry nodes. - In the case where value_mean < threshold before limiting, there could still be dry nodes - afterwards due to the logic of the limiter. - """ - struct PositivityPreservingLimiterShallowWater{N, Thresholds<:NTuple{N,<:Real}, Variables<:NTuple{N,Any}} - thresholds::Thresholds - variables::Variables - end - - function PositivityPreservingLimiterShallowWater(; thresholds, variables) - PositivityPreservingLimiterShallowWater(thresholds, variables) - end - - - function (limiter!::PositivityPreservingLimiterShallowWater)( - u_ode, integrator, semi::AbstractSemidiscretization, t) - u = wrap_array(u_ode, semi) - @trixi_timeit timer() "positivity-preserving limiter" limiter_shallow_water!( - u, limiter!.thresholds, limiter!.variables, mesh_equations_solver_cache(semi)...) - end - - - # Iterate over tuples in a type-stable way using "lispy tuple programming", - # similar to https://stackoverflow.com/a/55849398: - # Iterating over tuples of different functions isn't type-stable in general - # but accessing the first element of a tuple is type-stable. Hence, it's good - # to process one element at a time and replace iteration by recursion here. - # Note that you shouldn't use this with too many elements per tuple since the - # compile times can increase otherwise - but a handful of elements per tuple - # is definitely fine. - function limiter_shallow_water!(u, thresholds::NTuple{N,<:Real}, variables::NTuple{N,Any}, - mesh, equations::Union{ShallowWaterEquations1D, ShallowWaterEquations2D}, - solver, cache) where {N} - threshold = first(thresholds) - remaining_thresholds = Base.tail(thresholds) - variable = first(variables) - remaining_variables = Base.tail(variables) - - limiter_shallow_water!(u, threshold, variable, mesh, equations, solver, cache) - limiter_shallow_water!(u, remaining_thresholds, remaining_variables, mesh, equations, solver, cache) - return nothing - end - - # terminate the type-stable iteration over tuples - function limiter_shallow_water!(u, thresholds::Tuple{}, variables::Tuple{}, - mesh, equations::Union{ShallowWaterEquations1D, ShallowWaterEquations2D}, - solver, cache) - nothing - end +""" + PositivityPreservingLimiterShallowWater(; threshold, variables) - include("positivity_shallow_water_dg1d.jl") +The limiter is specifically designed for `ShallowWaterEquations`. +It is applied to all scalar `variables` in their given order +using the associated `thresholds` to determine the minimal acceptable values. +The order of the `variables` is important and might have a strong influence +on the robustness. +As opposed to the standard version of the [`PositivityPreservingLimiterZhangShu`](@ref), +nodes with a water height below the `threshold_limiter` are treated in a special way. +To avoid numerical problems caused by velocities close to zero, +the velocity is cut off, such that the node can be identified as "dry". The special feature of the +`ShallowWaterEquations` used here is that the bottom topography is stored as an additional +quantity in the solution vector `u`. However, the value of the bottom topography +should not be changed. That is why, it is not limited. +After the limiting process is applied to all degrees of freedom, for safety reasons, +the wet/dry threshold is applied again on all the DG nodes in order to avoid dry nodes. +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. +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) +""" +struct PositivityPreservingLimiterShallowWater{N, Thresholds<:NTuple{N,<:Real}, Variables<:NTuple{N,Any}} + thresholds::Thresholds + variables::Variables +end +function PositivityPreservingLimiterShallowWater(; thresholds, variables) + PositivityPreservingLimiterShallowWater(thresholds, variables) +end - end # @muladd \ No newline at end of file +function (limiter!::PositivityPreservingLimiterShallowWater)( + u_ode, integrator, semi::AbstractSemidiscretization, t) + u = wrap_array(u_ode, semi) + @trixi_timeit timer() "positivity-preserving limiter" limiter_shallow_water!( + u, limiter!.thresholds, limiter!.variables, mesh_equations_solver_cache(semi)...) +end + + +# Iterate over tuples in a type-stable way using "lispy tuple programming", +# similar to https://stackoverflow.com/a/55849398: +# Iterating over tuples of different functions isn't type-stable in general +# but accessing the first element of a tuple is type-stable. Hence, it's good +# to process one element at a time and replace iteration by recursion here. +# Note that you shouldn't use this with too many elements per tuple since the +# compile times can increase otherwise - but a handful of elements per tuple +# is definitely fine. +function limiter_shallow_water!(u, thresholds::NTuple{N,<:Real}, variables::NTuple{N,Any}, + mesh, equations::Union{ShallowWaterEquations1D, ShallowWaterEquations2D}, + solver, cache) where {N} + threshold = first(thresholds) + remaining_thresholds = Base.tail(thresholds) + variable = first(variables) + remaining_variables = Base.tail(variables) + + limiter_shallow_water!(u, threshold, variable, mesh, equations, solver, cache) + limiter_shallow_water!(u, remaining_thresholds, remaining_variables, mesh, equations, solver, cache) + return nothing +end + +# terminate the type-stable iteration over tuples +function limiter_shallow_water!(u, thresholds::Tuple{}, variables::Tuple{}, + mesh, equations::Union{ShallowWaterEquations1D, ShallowWaterEquations2D}, + solver, cache) + nothing +end + +include("positivity_shallow_water_dg1d.jl") + + + +end # @muladd \ No newline at end of file From 03daf4f84878b4f4e8c469edd87e94dcfa015d6f Mon Sep 17 00:00:00 2001 From: svengoldberg Date: Thu, 18 May 2023 15:35:03 +0200 Subject: [PATCH 32/44] Replace boolean with if, remove set_node_vars Shorten test runs on TreeMesh and UnstructuredMesh --- .../positivity_shallow_water_dg1d.jl | 17 ++++++++++------- src/equations/shallow_water_1d.jl | 12 ++++++++---- src/solvers/dg.jl | 7 ------- src/solvers/dgsem_tree/indicators_1d.jl | 8 ++++++-- test/test_tree_1d_shallowwater.jl | 12 ++++++------ 5 files changed, 30 insertions(+), 26 deletions(-) diff --git a/src/callbacks_stage/positivity_shallow_water_dg1d.jl b/src/callbacks_stage/positivity_shallow_water_dg1d.jl index d444df88f9..3107ea958b 100644 --- a/src/callbacks_stage/positivity_shallow_water_dg1d.jl +++ b/src/callbacks_stage/positivity_shallow_water_dg1d.jl @@ -43,19 +43,20 @@ function limiter_shallow_water!(u, threshold::Real, variable, h_mean, h_v_mean, _ = u_mean # b_mean is not used as b_node must not be overwritten # Set them both to zero to apply linear combination correctly - h_v_node = h_v_node * Int32(h_node > threshold) - - h_v_mean = h_v_mean * Int32(h_node > threshold) + if h_node <= threshold + h_v_node = 0 + h_v_mean = 0 + end u_node = SVector(h_node, h_v_node, b_node) - u_mean = SVector(h_mean, h_v_mean, zero(eltype(u))) + u_mean = SVector(h_mean, h_v_mean, b_node) # When velocity is cut off, the only averaged value is the waterheight, # because the velocity is set to zero and this value is passed. # Otherwise, the velocity is averaged, as well. # Note that the auxiliary bottom topography variable `b` is never limited. set_node_vars!(u, theta * u_node + (1-theta) * u_mean, - (1,2), dg, i, element) + equations, dg, i, element) end end @@ -70,8 +71,10 @@ function limiter_shallow_water!(u, threshold::Real, variable, h, hv, b = u_node - h = h * Int32(h > threshold) + threshold * Int32(h <= threshold) - hv = hv * Int32(h > threshold) + if h <= threshold + h = threshold + hv = 0 + end u_node = SVector(h, hv, b) diff --git a/src/equations/shallow_water_1d.jl b/src/equations/shallow_water_1d.jl index 9f68ffb654..d0f25bce2c 100644 --- a/src/equations/shallow_water_1d.jl +++ b/src/equations/shallow_water_1d.jl @@ -506,11 +506,15 @@ Further details on this hydrostatic reconstruction and its motivation can be fou # The default value is set to 1e-15 and can be changed within the constructor call in an elixir. threshold = equations.threshold_wet - h_ll_star = h_ll_star * Int32(h_ll_star > threshold) + threshold * Int32(h_ll_star <= threshold) - h_rr_star = h_rr_star * Int32(h_rr_star > threshold) + threshold * Int32(h_rr_star <= threshold) + if (h_ll_star <= threshold) + h_ll_star = threshold + v_ll = 0 + end - v_ll = v_ll * Int32(h_ll_star > threshold) - v_rr = v_rr * Int32(h_rr_star > threshold) + if (h_rr_star <= threshold) + h_rr_star = threshold + v_rr = 0 + end # Create the conservative variables using the reconstruted water heights u_ll_star = SVector( h_ll_star, h_ll_star * v_ll, b_ll ) diff --git a/src/solvers/dg.jl b/src/solvers/dg.jl index e6492cd5fe..29de151f14 100644 --- a/src/solvers/dg.jl +++ b/src/solvers/dg.jl @@ -491,13 +491,6 @@ end return nothing end -@inline function set_node_vars!(u, u_node, variables::NTuple{N, Integer}, solver::DG, indices...) where N - for v in variables - u[v, indices...] = u_node[v] - end - return nothing -end - @inline function add_to_node_vars!(u, u_node, equations, solver::DG, indices...) for v in eachvariable(equations) u[v, indices...] += u_node[v] diff --git a/src/solvers/dgsem_tree/indicators_1d.jl b/src/solvers/dgsem_tree/indicators_1d.jl index 81dfb132ec..19f702b805 100644 --- a/src/solvers/dgsem_tree/indicators_1d.jl +++ b/src/solvers/dgsem_tree/indicators_1d.jl @@ -62,13 +62,17 @@ function (indicator_hg::IndicatorHennemannGassner)(u, mesh::Union{TreeMesh{1}, S modal = modal_threaded[Threads.threadid()] # (Re-)set dummy variable for alpha_dry - indicator_wet = 1.0 + indicator_wet = 1 # 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 - indicator_wet = min(indicator_wet, Int32(h > threshold_wet)) + + if h <= threshold_wet + indicator_wet = 0 + end + indicator[i] = indicator_hg.variable(u_local, equations) end diff --git a/test/test_tree_1d_shallowwater.jl b/test/test_tree_1d_shallowwater.jl index 09f9b338dd..804de914e4 100644 --- a/test/test_tree_1d_shallowwater.jl +++ b/test/test_tree_1d_shallowwater.jl @@ -90,16 +90,16 @@ EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "tree_1 @trixi_testset "elixir_shallowwater_beach.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_beach.jl"), - l2 = [0.5837835144970676, 3.740558401311243, 6.289710785696625e-8], - linf = [0.9917277949257116, 7.096293347373189, 4.452604027704865e-7], - tspan = (0.0, 0.25)) + l2 = [0.1797896888850151, 1.2377413339448482, 6.290434242536467e-8], + linf = [0.8460686874664077, 3.3741901787491937, 4.4526743714357053e-7], + tspan = (0.0, 0.05)) end @trixi_testset "elixir_shallowwater_parabolic_bowl.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_parabolic_bowl.jl"), - l2 = [0.0020629352664007466, 0.0004816398224727313, 3.828565374146696e-17], - linf = [0.0038611903839202823, 0.0008591988088082134, 1.6653345369377348e-16], - tspan = (0.0, 0.25)) + l2 = [8.96598168303412e-5, 1.85657073978462e-5, 3.828565374146696e-17], + linf = [0.00041080213812394266, 0.0001482326148893917, 1.6653345369377348e-16], + tspan = (0.0, 0.05)) end end From bdc8c7243727592cd56761aa98eabaa4cceb3bd3 Mon Sep 17 00:00:00 2001 From: svengoldberg Date: Fri, 19 May 2023 15:30:23 +0200 Subject: [PATCH 33/44] Update comment regarding H0 for lake-at-rest error --- src/equations/shallow_water_1d.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/equations/shallow_water_1d.jl b/src/equations/shallow_water_1d.jl index d0f25bce2c..859f15ab81 100644 --- a/src/equations/shallow_water_1d.jl +++ b/src/equations/shallow_water_1d.jl @@ -62,7 +62,8 @@ end # Allow for flexibility to set the gravitational constant within an elixir depending on the # application where `gravity_constant=1.0` or `gravity_constant=9.81` are common values. -# The reference total water height H0 is an artefact from the old calculation of the lake_at_rest_error +# The reference total water height H0 defaults to 0.0 but is used for the "lake-at-rest" +# well-balancedness test cases. # Strict default values for thresholds that performed great in several numerical experiments function ShallowWaterEquations1D(; gravity_constant, H0=0.0, threshold_limiter=1e-13, threshold_wet=1e-15) From ef094efa878b480a8257827910828b0285e298a4 Mon Sep 17 00:00:00 2001 From: svengoldberg Date: Fri, 19 May 2023 15:38:36 +0200 Subject: [PATCH 34/44] Add the original source to the parabolic bowl test --- .../elixir_shallowwater_parabolic_bowl.jl | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_parabolic_bowl.jl b/examples/tree_1d_dgsem/elixir_shallowwater_parabolic_bowl.jl index 7f43346e0a..ffcdcc76c4 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_parabolic_bowl.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_parabolic_bowl.jl @@ -14,11 +14,16 @@ Well-known initial condition to test the [`hydrostatic_reconstruction_chen_noell 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 are taken out of section 6.2 from the paper: +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\n + [DOI: 10.1017/S0022112081001882](https://doi.org/10.1017/S0022112081001882). + +The particular setup is used as the one-dimensional version out of section 6.2 from the paper: - Niklas Wintermeyer, Andrew R. Winters, Gregor J. Gassner and Timothy Warburton (2018) - An entropy stable discontinuous Galerkin method for the shallow water equations on + An entropy stable discontinuous Galerkin method for the shallow water equations on curvilinear meshes with wet/dry fronts accelerated by GPUs\n - [DOI: 10.1016/j.jcp.2018.08.038](https://doi.org/10.1016/j.jcp.2018.08.038) + [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:: ShallowWaterEquations1D) a = 1 From f99b0167ee3df0edfd7ce2d5f6b294aa0de51f63 Mon Sep 17 00:00:00 2001 From: svengoldberg Date: Thu, 15 Jun 2023 15:01:09 +0200 Subject: [PATCH 35/44] Outsource HG shock capturing indicator for SWE Create different function to compute indicator Edit comments Change wet/dry clipping to if-else logic --- .../elixir_shallowwater_beach.jl | 10 +-- .../elixir_shallowwater_parabolic_bowl.jl | 10 +-- ...ixir_shallowwater_well_balanced_wet_dry.jl | 10 +-- src/Trixi.jl | 2 +- src/solvers/dgsem_tree/indicators.jl | 73 +++++++++++++++++++ src/solvers/dgsem_tree/indicators_1d.jl | 37 ++++++---- 6 files changed, 110 insertions(+), 32 deletions(-) diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_beach.jl b/examples/tree_1d_dgsem/elixir_shallowwater_beach.jl index e7c8a36c29..482d89c4fc 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_beach.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_beach.jl @@ -58,11 +58,11 @@ surface_flux = (FluxHydrostaticReconstruction(flux_hll_chen_noelle, hydrostatic_ basis = LobattoLegendreBasis(3) -indicator_sc = IndicatorHennemannGassner(equations, basis, - alpha_max=0.5, - alpha_min=0.001, - alpha_smooth=true, - variable=waterheight_pressure) +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) diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_parabolic_bowl.jl b/examples/tree_1d_dgsem/elixir_shallowwater_parabolic_bowl.jl index ffcdcc76c4..9023d01eb6 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_parabolic_bowl.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_parabolic_bowl.jl @@ -57,11 +57,11 @@ surface_flux = (FluxHydrostaticReconstruction(flux_hll_chen_noelle, hydrostatic_ basis = LobattoLegendreBasis(5) -indicator_sc = IndicatorHennemannGassner(equations, basis, - alpha_max=0.5, - alpha_min=0.001, - alpha_smooth=true, - variable=waterheight_pressure) +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) diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl b/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl index e561091ca4..f68d4887a4 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl @@ -53,11 +53,11 @@ surface_flux = (FluxHydrostaticReconstruction(flux_hll_chen_noelle, hydrostatic_ basis = LobattoLegendreBasis(3) -indicator_sc = IndicatorHennemannGassner(equations, basis, - alpha_max=0.5, - alpha_min=0.001, - alpha_smooth=true, - variable=waterheight_pressure) +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) diff --git a/src/Trixi.jl b/src/Trixi.jl index 762d436cc2..81fc872b4a 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -201,7 +201,7 @@ export DG, VolumeIntegralWeakForm, VolumeIntegralStrongForm, VolumeIntegralFluxDifferencing, VolumeIntegralPureLGLFiniteVolume, - VolumeIntegralShockCapturingHG, IndicatorHennemannGassner, + VolumeIntegralShockCapturingHG, IndicatorHennemannGassner, IndicatorHennemannGassnerShallowWater, VolumeIntegralUpwind, SurfaceIntegralWeakForm, SurfaceIntegralStrongForm, SurfaceIntegralUpwind, diff --git a/src/solvers/dgsem_tree/indicators.jl b/src/solvers/dgsem_tree/indicators.jl index 30d3b2c044..3ff14a5ff0 100644 --- a/src/solvers/dgsem_tree/indicators.jl +++ b/src/solvers/dgsem_tree/indicators.jl @@ -136,6 +136,79 @@ function (indicator_hg::IndicatorHennemannGassner)(u, mesh, equations, dg::DGSEM return alpha end +""" + IndicatorHennemannGassnerShallowWater(equations::AbstractEquations, basis; + alpha_max=0.5, + alpha_min=0.001, + alpha_smooth=true, + variable) + +Modified version of the [`IndicatorHennemannGassner`](@ref) +indicator used for shock-capturing for shallow water equations. After +the element-wise values for the blending factors are computed an additional check +is made to see if the element is partially wet. In this case, partially wet elements +are set to use the pure finite volume scheme that is guaranteed to be well-balanced +for this wet/dry transition state of the flow regime. + +See also [`VolumeIntegralShockCapturingHG`](@ref). + +## References + +- Hennemann, Gassner (2020) + "A provably entropy stable subcell shock capturing approach for high order split form DG" + [arXiv: 2008.12044](https://arxiv.org/abs/2008.12044) +""" +struct IndicatorHennemannGassnerShallowWater{RealT<:Real, Variable, Cache} <: AbstractIndicator + alpha_max::RealT + alpha_min::RealT + alpha_smooth::Bool + variable::Variable + cache::Cache +end + +# this method is used when the indicator is constructed as for shock-capturing volume integrals +# of the shallow water equations +# It modifies the shock-capturing indicator to use full FV method in dry cells +function IndicatorHennemannGassnerShallowWater(equations::AbstractShallowWaterEquations, basis; + alpha_max=0.5, + alpha_min=0.001, + alpha_smooth=true, + variable) + alpha_max, alpha_min = promote(alpha_max, alpha_min) + cache = create_cache(IndicatorHennemannGassner, equations, basis) + IndicatorHennemannGassnerShallowWater{typeof(alpha_max), typeof(variable), typeof(cache)}( + alpha_max, alpha_min, alpha_smooth, variable, cache) +end + + + +function Base.show(io::IO, indicator::IndicatorHennemannGassnerShallowWater) + @nospecialize indicator # reduce precompilation time + + print(io, "IndicatorHennemannGassnerShallowWater(") + print(io, indicator.variable) + print(io, ", alpha_max=", indicator.alpha_max) + print(io, ", alpha_min=", indicator.alpha_min) + print(io, ", alpha_smooth=", indicator.alpha_smooth) + print(io, ")") +end + +function Base.show(io::IO, ::MIME"text/plain", indicator::IndicatorHennemannGassnerShallowWater) + @nospecialize indicator # reduce precompilation time + + if get(io, :compact, false) + show(io, indicator) + else + setup = [ + "indicator variable" => indicator.variable, + "max. α" => indicator.alpha_max, + "min. α" => indicator.alpha_min, + "smooth α" => (indicator.alpha_smooth ? "yes" : "no"), + ] + summary_box(io, "IndicatorHennemannGassnerShallowWater", setup) + end +end + """ IndicatorLöhner (equivalent to IndicatorLoehner) diff --git a/src/solvers/dgsem_tree/indicators_1d.jl b/src/solvers/dgsem_tree/indicators_1d.jl index 59f7e173c0..e186019471 100644 --- a/src/solvers/dgsem_tree/indicators_1d.jl +++ b/src/solvers/dgsem_tree/indicators_1d.jl @@ -25,9 +25,9 @@ end # Overload indicator when ShallowWaterEquations1D is used to apply full FV method on cells # containing dry LGL nodes -function (indicator_hg::IndicatorHennemannGassner)(u, mesh::Union{TreeMesh{1}, StructuredMesh{1}}, - equations::ShallowWaterEquations1D, dg::DGSEM, - cache; kwargs...) +function (indicator_hg::IndicatorHennemannGassnerShallowWater)(u, mesh::Union{TreeMesh{1}, StructuredMesh{1}}, + equations::ShallowWaterEquations1D, dg::DGSEM, + cache; kwargs...) @unpack alpha_max, alpha_min, alpha_smooth, variable = indicator_hg @unpack alpha, alpha_tmp, indicator_threaded, modal_threaded = indicator_hg.cache # TODO: Taal refactor, when to `resize!` stuff changed possibly by AMR? @@ -42,19 +42,20 @@ function (indicator_hg::IndicatorHennemannGassner)(u, mesh::Union{TreeMesh{1}, S threshold = 0.5 * 10^(-1.8 * (nnodes(dg))^0.25) parameter_s = log((1 - 0.0001)/0.0001) - # If h at one LGL node is lower than threshold_wet, set alpha[element]=1 via indicator_wet - # in order to apply a full FV method in the partially wet cell and guarantee the well-balanced property. + # If the water height `h` at one LGL node is lower than `threshold_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 cells and guarantees the well-balanced property. # - # Determination of hard coded threshold: Showed good results in many numerical experiments. - # The idea is to gain more stability when computing the velocity on (nearly) dry cells which - # could be counteracted by division of conservative variables, e.g., v = hv / h. - # Here, the impact of the threshold on the number of cells being updated with FV is not that - # significant. However, its impact on the stability is very significant. - # The value can be seen as a trade-off between accuracy and stability. - # Well-balancedness of the scheme on partially wet cells 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. - # + # Hard-coded cut-off value of `threshold_wet = 1e-4` was determined through many numerical experiments. + # Overall idea is to increase robustness when computing the velocity on (nearly) dry cells 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 cells 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 cells 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_wet = 1e-4 @threaded for element in eachelement(dg, cache) @@ -110,7 +111,11 @@ function (indicator_hg::IndicatorHennemannGassner)(u, mesh::Union{TreeMesh{1}, S end # Clip the maximum amount of FV allowed or set to one depending on indicator_wet - alpha[element] = indicator_wet * min(alpha_max, alpha_element) - (indicator_wet-1) + if indicator_wet == 0 + alpha[element] = 1 + else # Element is not defined as dry but wet + alpha[element] = min(alpha_max, alpha_element) + end end if alpha_smooth From abade120c401f2665ac1fe5237b7529b22c43134 Mon Sep 17 00:00:00 2001 From: svengoldberg Date: Thu, 15 Jun 2023 15:06:56 +0200 Subject: [PATCH 36/44] Move limiterthreshold into function & edit docs Threshold was a passed variable in elixir before. Now, it is taken right from the SWE struct in the limiter Edit docs --- .../elixir_shallowwater_beach.jl | 3 +- .../elixir_shallowwater_parabolic_bowl.jl | 3 +- ...ixir_shallowwater_well_balanced_wet_dry.jl | 3 +- .../positivity_shallow_water.jl | 44 +++++++++---------- 4 files changed, 24 insertions(+), 29 deletions(-) diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_beach.jl b/examples/tree_1d_dgsem/elixir_shallowwater_beach.jl index 482d89c4fc..57c731aafb 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_beach.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_beach.jl @@ -106,8 +106,7 @@ save_solution = SaveSolutionCallback(interval=250, callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution) -stage_limiter! = PositivityPreservingLimiterShallowWater(thresholds=(equations.threshold_limiter,), - variables=(Trixi.waterheight,)) +stage_limiter! = PositivityPreservingLimiterShallowWater(variables=(Trixi.waterheight,)) ############################################################################### # run the simulation diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_parabolic_bowl.jl b/examples/tree_1d_dgsem/elixir_shallowwater_parabolic_bowl.jl index 9023d01eb6..4d4449d94b 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_parabolic_bowl.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_parabolic_bowl.jl @@ -103,8 +103,7 @@ save_solution = SaveSolutionCallback(interval=1000, callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution) -stage_limiter! = PositivityPreservingLimiterShallowWater(thresholds=(equations.threshold_limiter,), - variables=(Trixi.waterheight,)) +stage_limiter! = PositivityPreservingLimiterShallowWater(variables=(Trixi.waterheight,)) ############################################################################### # run the simulation diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl b/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl index f68d4887a4..6fbb8de29e 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl @@ -128,8 +128,7 @@ stepsize_callback = StepsizeCallback(cfl=1.5) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, stepsize_callback) -stage_limiter! = PositivityPreservingLimiterShallowWater(thresholds=(equations.threshold_limiter,), - variables=(Trixi.waterheight,)) +stage_limiter! = PositivityPreservingLimiterShallowWater(variables=(Trixi.waterheight,)) ############################################################################### # run the simulation diff --git a/src/callbacks_stage/positivity_shallow_water.jl b/src/callbacks_stage/positivity_shallow_water.jl index d1e2792659..62f0020d7e 100644 --- a/src/callbacks_stage/positivity_shallow_water.jl +++ b/src/callbacks_stage/positivity_shallow_water.jl @@ -1,42 +1,42 @@ # By default, Julia/LLVM does not use fused multiply-add operations (FMAs). - # Since these FMAs can increase the performance of many numerical algorithms, - # we need to opt-in explicitly. - # See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. - @muladd begin +# Since these FMAs can increase the performance of many numerical algorithms, +# we need to opt-in explicitly. +# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. +@muladd begin """ - PositivityPreservingLimiterShallowWater(; threshold, variables) + PositivityPreservingLimiterShallowWater(; variables) -The limiter is specifically designed for `ShallowWaterEquations`. +The limiter is specifically designed for the shallow water equations. It is applied to all scalar `variables` in their given order -using the associated `thresholds` to determine the minimal acceptable values. +using the defined `threshold_limiter` from the [`ShallowWaterEquations1D`](@ref) struct +or the [`ShallowWaterEquations2D`](@ref) struct to determine the minimal acceptable values. The order of the `variables` is important and might have a strong influence on the robustness. -As opposed to the standard version of the [`PositivityPreservingLimiterZhangShu`](@ref), -nodes with a water height below the `threshold_limiter` are treated in a special way. +As opposed to the standard version of the [`PositivityPreservingLimiterZhangShu`](@ref), +nodes with a water height below the `threshold_limiter` are treated in a special way. To avoid numerical problems caused by velocities close to zero, the velocity is cut off, such that the node can be identified as "dry". The special feature of the `ShallowWaterEquations` used here is that the bottom topography is stored as an additional quantity in the solution vector `u`. However, the value of the bottom topography should not be changed. That is why, it is not limited. After the limiting process is applied to all degrees of freedom, for safety reasons, -the wet/dry threshold is applied again on all the DG nodes in order to avoid dry nodes. +the `threshold_limiter` is applied again on all the DG nodes in order to avoid water height below. 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. -This fully-discrete positivity-preserving limiter is based on the work of +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) """ -struct PositivityPreservingLimiterShallowWater{N, Thresholds<:NTuple{N,<:Real}, Variables<:NTuple{N,Any}} - thresholds::Thresholds +struct PositivityPreservingLimiterShallowWater{N, Variables<:NTuple{N,Any}} variables::Variables end -function PositivityPreservingLimiterShallowWater(; thresholds, variables) - PositivityPreservingLimiterShallowWater(thresholds, variables) +function PositivityPreservingLimiterShallowWater(;variables) + PositivityPreservingLimiterShallowWater(variables) end @@ -44,7 +44,7 @@ function (limiter!::PositivityPreservingLimiterShallowWater)( u_ode, integrator, semi::AbstractSemidiscretization, t) u = wrap_array(u_ode, semi) @trixi_timeit timer() "positivity-preserving limiter" limiter_shallow_water!( - u, limiter!.thresholds, limiter!.variables, mesh_equations_solver_cache(semi)...) + u, limiter!.variables, mesh_equations_solver_cache(semi)...) end @@ -56,22 +56,20 @@ end # Note that you shouldn't use this with too many elements per tuple since the # compile times can increase otherwise - but a handful of elements per tuple # is definitely fine. -function limiter_shallow_water!(u, thresholds::NTuple{N,<:Real}, variables::NTuple{N,Any}, +function limiter_shallow_water!(u, variables::NTuple{N,Any}, mesh, equations::Union{ShallowWaterEquations1D, ShallowWaterEquations2D}, solver, cache) where {N} - threshold = first(thresholds) - remaining_thresholds = Base.tail(thresholds) variable = first(variables) remaining_variables = Base.tail(variables) - limiter_shallow_water!(u, threshold, variable, mesh, equations, solver, cache) - limiter_shallow_water!(u, remaining_thresholds, remaining_variables, mesh, equations, solver, cache) + limiter_shallow_water!(u, equations.threshold_limiter, variable, mesh, equations, solver, cache) + limiter_shallow_water!(u, remaining_variables, mesh, equations, solver, cache) return nothing end # terminate the type-stable iteration over tuples -function limiter_shallow_water!(u, thresholds::Tuple{}, variables::Tuple{}, - mesh, equations::Union{ShallowWaterEquations1D, ShallowWaterEquations2D}, +function limiter_shallow_water!(u, variables::Tuple{}, + mesh, equations::Union{ShallowWaterEquations1D, ShallowWaterEquations2D}, solver, cache) nothing end From 3c306b1d9788df577b883fbc3ec8e422d88210aa Mon Sep 17 00:00:00 2001 From: svengoldberg Date: Thu, 15 Jun 2023 15:12:39 +0200 Subject: [PATCH 37/44] Move new limiter safety check in same element loop --- src/callbacks_stage/positivity_shallow_water_dg1d.jl | 12 +++++------- 1 file changed, 5 insertions(+), 7 deletions(-) diff --git a/src/callbacks_stage/positivity_shallow_water_dg1d.jl b/src/callbacks_stage/positivity_shallow_water_dg1d.jl index 3107ea958b..eddfb723e1 100644 --- a/src/callbacks_stage/positivity_shallow_water_dg1d.jl +++ b/src/callbacks_stage/positivity_shallow_water_dg1d.jl @@ -55,16 +55,14 @@ function limiter_shallow_water!(u, threshold::Real, variable, # because the velocity is set to zero and this value is passed. # Otherwise, the velocity is averaged, as well. # Note that the auxiliary bottom topography variable `b` is never limited. - set_node_vars!(u, theta * u_node + (1-theta) * u_mean, + set_node_vars!(u, theta * u_node + (1 - theta) * u_mean, equations, dg, i, element) end - end - # An extra "safety" check is done over all the degrees of - # freedom after the limiting 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 limiter. - @threaded for element in eachelement(dg, cache) + # An extra "safety" check is done over all the degrees of + # freedom after the limiting 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 limiter. for i in eachnode(dg) u_node = get_node_vars(u, equations, dg, i, element) From 096d58e115a4255192d9f914ff7e1ec2d1fe215e Mon Sep 17 00:00:00 2001 From: svengoldberg Date: Thu, 15 Jun 2023 15:21:54 +0200 Subject: [PATCH 38/44] Adjust default threshold values --- .../elixir_shallowwater_beach.jl | 5 ++-- .../elixir_shallowwater_parabolic_bowl.jl | 5 ++-- ...ixir_shallowwater_well_balanced_wet_dry.jl | 3 +- src/equations/shallow_water_1d.jl | 29 +++++++++++++------ test/test_tree_1d_shallowwater.jl | 12 ++++---- 5 files changed, 34 insertions(+), 20 deletions(-) diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_beach.jl b/examples/tree_1d_dgsem/elixir_shallowwater_beach.jl index 57c731aafb..d025430b48 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_beach.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_beach.jl @@ -38,9 +38,10 @@ function initial_condition_beach(x, t, equations:: ShallowWaterEquations1D) 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 + # 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 - # (1e-13 per default, set in the constructor for the ShallowWaterEquations) is added if h = 0. + # with a default value of 500*eps() ≈ 1e-13 in double precision, is set in the constructor above + # for the ShallowWaterEquations 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) diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_parabolic_bowl.jl b/examples/tree_1d_dgsem/elixir_shallowwater_parabolic_bowl.jl index 4d4449d94b..e0acaeb24f 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_parabolic_bowl.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_parabolic_bowl.jl @@ -38,9 +38,10 @@ function initial_condition_parabolic_bowl(x, t, equations:: ShallowWaterEquation 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 + # 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 - # (1e-13 per default, set in the constructor for the ShallowWaterEquations) is added if h = 0. + # with a default value of 500*eps() ≈ 1e-13 in double precision, is set in the constructor above + # for the ShallowWaterEquations 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) diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl b/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl index 6fbb8de29e..ff1404663b 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl @@ -36,7 +36,8 @@ function initial_condition_complex_bottom_well_balanced(x, t, equations:: Shallo # 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 - # (1e-13 per default, set in the constructor for the ShallowWaterEquations) is added if h = 0. + # with a default value of 500*eps() ≈ 1e-13 in double precision, is set in the constructor above + # for the ShallowWaterEquations 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) diff --git a/src/equations/shallow_water_1d.jl b/src/equations/shallow_water_1d.jl index 859f15ab81..69ccad2b11 100644 --- a/src/equations/shallow_water_1d.jl +++ b/src/equations/shallow_water_1d.jl @@ -53,20 +53,29 @@ References for the SWE are many but a good introduction is available in Chapter struct ShallowWaterEquations1D{RealT<:Real} <: AbstractShallowWaterEquations{1, 3} gravity::RealT # gravitational constant H0::RealT # constant "lake-at-rest" total water height - threshold_limiter::RealT # Threshold to use in PositivityPreservingLimiterShallowWater on water height, - # as a (small) shift on the initial condition and cutoff before the - # next time step. - threshold_wet::RealT # Threshold to be applied on water height to define when the flow is "wet" - # before calculating the numerical flux. + threshold_limiter::RealT # Threshold to use in `PositivityPreservingLimiterShallowWater` on water height, + # as a (small) shift on the initial condition and cutoff before the + # next time step. + # Default in double precision is 500*eps() ≈ 1e-13. + threshold_wet::RealT # Threshold to be applied on water height to define when the flow is "wet" + # before calculating the numerical flux. + # Default in double precision is 5*eps() ≈ 1e-15. end # Allow for flexibility to set the gravitational constant within an elixir depending on the # application where `gravity_constant=1.0` or `gravity_constant=9.81` are common values. # The reference total water height H0 defaults to 0.0 but is used for the "lake-at-rest" # well-balancedness test cases. -# Strict default values for thresholds that performed great in several numerical experiments -function ShallowWaterEquations1D(; gravity_constant, H0=0.0, - threshold_limiter=1e-13, threshold_wet=1e-15) +# Strict default values for thresholds that performed well in many numerical experiments +function ShallowWaterEquations1D(; gravity_constant, H0=zero(gravity_constant), + threshold_limiter=nothing, threshold_wet=nothing) + T = promote_type(typeof(gravity_constant), typeof(H0)) + if threshold_limiter === nothing + threshold_limiter = 500 * eps(T) + end + if threshold_wet === nothing + threshold_wet = 5 * eps(T) + end ShallowWaterEquations1D(gravity_constant, H0, threshold_limiter, threshold_wet) end @@ -504,7 +513,9 @@ Further details on this hydrostatic reconstruction and its motivation can be fou # the hydrostatic reconstruction is applied and before the numerical flux is calculated # to avoid numerical problem with arbitrary small values. Interfaces with a water height # lower or equal to the threshold can be declared as dry. - # The default value is set to 1e-15 and can be changed within the constructor call in an elixir. + # The default value for `threshold_wet` is ≈ 5*eps(), or 1e-15 in double precision, is set + # in the `ShallowWaterEquations1D` struct. This threshold value can be changed in the constructor + # call of this eqaution struct in an elixir. threshold = equations.threshold_wet if (h_ll_star <= threshold) diff --git a/test/test_tree_1d_shallowwater.jl b/test/test_tree_1d_shallowwater.jl index d4fdaa29b9..b0914db986 100644 --- a/test/test_tree_1d_shallowwater.jl +++ b/test/test_tree_1d_shallowwater.jl @@ -32,8 +32,8 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem") @trixi_testset "elixir_shallowwater_well_balanced_wet_dry.jl with FluxHydrostaticReconstruction" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_well_balanced_wet_dry.jl"), - l2 = [0.009657871671690443, 5.403349017029091e-15, 0.03857583749209983], - linf = [0.49999999999990025, 3.344735163283149e-14, 1.9999999999999998], + l2 = [0.00965787167169024, 5.345454081916856e-14, 0.03857583749209928], + linf = [0.4999999999998892, 2.2447689894899726e-13, 1.9999999999999714], tspan = (0.0, 0.25)) end @@ -90,15 +90,15 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem") @trixi_testset "elixir_shallowwater_beach.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_beach.jl"), - l2 = [0.1797896888850151, 1.2377413339448482, 6.290434242536467e-8], - linf = [0.8460686874664077, 3.3741901787491937, 4.4526743714357053e-7], + l2 = [0.17978981087947196, 1.237741563904404, 6.289469352788183e-8], + linf = [0.8460686638571695, 3.374190220014994, 4.452604152049844e-7], tspan = (0.0, 0.05)) end @trixi_testset "elixir_shallowwater_parabolic_bowl.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_parabolic_bowl.jl"), - l2 = [8.96598168303412e-5, 1.85657073978462e-5, 3.828565374146696e-17], - linf = [0.00041080213812394266, 0.0001482326148893917, 1.6653345369377348e-16], + l2 = [8.965981683033589e-5, 1.8565707397810857e-5, 4.1043039226164336e-17], + linf = [0.00041080213807871235, 0.00014823261488938177, 2.220446049250313e-16], tspan = (0.0, 0.05)) end end From eeb4dc46509f2ca2a2e9990a8cec84e5f3bcca29 Mon Sep 17 00:00:00 2001 From: svengoldberg Date: Thu, 15 Jun 2023 15:24:13 +0200 Subject: [PATCH 39/44] Remove type instability --- src/callbacks_stage/positivity_shallow_water_dg1d.jl | 6 +++--- src/equations/shallow_water_1d.jl | 4 ++-- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/callbacks_stage/positivity_shallow_water_dg1d.jl b/src/callbacks_stage/positivity_shallow_water_dg1d.jl index eddfb723e1..37b801a54f 100644 --- a/src/callbacks_stage/positivity_shallow_water_dg1d.jl +++ b/src/callbacks_stage/positivity_shallow_water_dg1d.jl @@ -44,8 +44,8 @@ function limiter_shallow_water!(u, threshold::Real, variable, # Set them both to zero to apply linear combination correctly if h_node <= threshold - h_v_node = 0 - h_v_mean = 0 + h_v_node = zero(eltype(u)) + h_v_mean = zero(eltype(u)) end u_node = SVector(h_node, h_v_node, b_node) @@ -71,7 +71,7 @@ function limiter_shallow_water!(u, threshold::Real, variable, if h <= threshold h = threshold - hv = 0 + hv = zero(eltype(u)) end u_node = SVector(h, hv, b) diff --git a/src/equations/shallow_water_1d.jl b/src/equations/shallow_water_1d.jl index 69ccad2b11..8115aed376 100644 --- a/src/equations/shallow_water_1d.jl +++ b/src/equations/shallow_water_1d.jl @@ -520,12 +520,12 @@ Further details on this hydrostatic reconstruction and its motivation can be fou if (h_ll_star <= threshold) h_ll_star = threshold - v_ll = 0 + v_ll = zero(v_ll) end if (h_rr_star <= threshold) h_rr_star = threshold - v_rr = 0 + v_rr = zero(v_rr) end # Create the conservative variables using the reconstruted water heights From ac14ca3850e954e98d927994c33d8aa6b3e67a0b Mon Sep 17 00:00:00 2001 From: svengoldberg Date: Thu, 15 Jun 2023 15:26:57 +0200 Subject: [PATCH 40/44] Import Printf package for terminal output --- .../elixir_shallowwater_well_balanced_wet_dry.jl | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl b/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl index ff1404663b..62ac97a536 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl @@ -1,6 +1,7 @@ using OrdinaryDiffEq using Trixi +using Printf: @printf, @sprintf ############################################################################### # Semidiscretization of the shallow water equations @@ -181,10 +182,10 @@ end # report the well-balancedness lake-at-rest error to the screen println("─"^100) -println(" Lake-at-rest error for '", Trixi.get_name(equations), "' with ", Trixi.summary(solver), - " at final time " * Trixi.@sprintf("%10.8e", tspan[end])) +println(" Lake-at-rest error for '", Trixi.get_name(equations), "' with ", summary(solver), + " at final time " * @sprintf("%10.8e", tspan[end])) -Trixi.@printf(" %-12s:", Trixi.pretty_form_utf(lake_at_rest_error)) -Trixi.@printf(" % 10.8e", l1_well_balance_error) +@printf(" %-12s:", Trixi.pretty_form_utf(lake_at_rest_error)) +@printf(" % 10.8e", l1_well_balance_error) println() println("─"^100) \ No newline at end of file From a10e65d07471106d3863fd2f2ee90233771099c8 Mon Sep 17 00:00:00 2001 From: svengoldberg Date: Thu, 15 Jun 2023 15:27:28 +0200 Subject: [PATCH 41/44] Edit docs --- .../elixir_shallowwater_parabolic_bowl.jl | 20 +++++++++---------- ...ixir_shallowwater_well_balanced_wet_dry.jl | 6 +++--- src/equations/shallow_water_1d.jl | 4 ++-- 3 files changed, 15 insertions(+), 15 deletions(-) diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_parabolic_bowl.jl b/examples/tree_1d_dgsem/elixir_shallowwater_parabolic_bowl.jl index e0acaeb24f..b090e5cee2 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_parabolic_bowl.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_parabolic_bowl.jl @@ -11,19 +11,19 @@ equations = ShallowWaterEquations1D(gravity_constant=9.81) initial_condition_parabolic_bowl(x, t, equations:: ShallowWaterEquations1D) Well-known initial condition to test the [`hydrostatic_reconstruction_chen_noelle`](@ref) and its -wet-dry mechanics. This test has analytical solutions. The initial condition is defined by the +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\n - [DOI: 10.1017/S0022112081001882](https://doi.org/10.1017/S0022112081001882). - -The particular setup is used as the one-dimensional version out of section 6.2 from the paper: - - 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\n - [DOI: 10.1016/j.jcp.2018.08.038](https://doi.org/10.1016/j.jcp.2018.08.038). +- 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:: ShallowWaterEquations1D) a = 1 diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl b/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl index 62ac97a536..ac108a8c69 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl @@ -9,13 +9,13 @@ using Printf: @printf, @sprintf equations = ShallowWaterEquations1D(gravity_constant=9.812) """ -initial_condition_complex_bottom_well_balanced(x, t, equations:: ShallowWaterEquations1D) + initial_condition_complex_bottom_well_balanced(x, t, equations:: ShallowWaterEquations1D) Initial condition with a complex (discontinuous) bottom topography to test the well-balanced property for the [`hydrostatic_reconstruction_chen_noelle`](@ref) including dry areas within the domain. The errors from the analysis callback are not important but the error for this -lake at rest test case `∑|H0-(h+b)|` should be around machine roundoff. -The initial condition was found in the section 5.2 of the paper: +lake-at-rest test case `∑|H0-(h+b)|` should be around machine roundoff. +The initial condition is taken from Section 5.2 of the paper: - Guoxian Chen and Sebastian Noelle (2017) A new hydrostatic reconstruction scheme based on subcell reconstructions [DOI:10.1137/15M1053074](https://dx.doi.org/10.1137/15M1053074) diff --git a/src/equations/shallow_water_1d.jl b/src/equations/shallow_water_1d.jl index 8115aed376..9fe07e9268 100644 --- a/src/equations/shallow_water_1d.jl +++ b/src/equations/shallow_water_1d.jl @@ -26,7 +26,7 @@ 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. Both of them do not 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 +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. @@ -448,7 +448,7 @@ end A particular type of hydrostatic reconstruction on the water height to guarantee well-balancedness for a general bottom topography [`ShallowWaterEquations1D`](@ref). The reconstructed solution states -`u_ll_star` and `u_rr_star` variables are used to evaluate the surface numerical flux at the interface. +`u_ll_star` and `u_rr_star` variables are then used to evaluate the surface numerical flux at the interface. Use in combination with the generic numerical flux routine [`FluxHydrostaticReconstruction`](@ref). Further details on this hydrostatic reconstruction and its motivation can be found in From 10e5012463e2ff92bcac5ac2661997e41dd7fb8d Mon Sep 17 00:00:00 2001 From: svengoldberg Date: Tue, 20 Jun 2023 21:46:32 +0200 Subject: [PATCH 42/44] Add Printf package to the test/Project.toml Used for printing lake-at-rest error in well-balancedness test --- test/Project.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/test/Project.toml b/test/Project.toml index 7d38641522..cae1d4ff39 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -17,6 +17,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" +Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" From d8e45c8cfd4d06540a9507610a3cb41128a66f40 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Thu, 6 Jul 2023 07:34:52 +0200 Subject: [PATCH 43/44] apply JuliaFormatter --- src/equations/shallow_water_1d.jl | 192 ++++++++++++------------ src/solvers/dgsem_tree/indicators_1d.jl | 175 ++++++++++----------- 2 files changed, 187 insertions(+), 180 deletions(-) diff --git a/src/equations/shallow_water_1d.jl b/src/equations/shallow_water_1d.jl index 44c69625f9..938f654fa6 100644 --- a/src/equations/shallow_water_1d.jl +++ b/src/equations/shallow_water_1d.jl @@ -50,17 +50,17 @@ References for the SWE are many but a good introduction is available in Chapter Finite Volume Methods for Hyperbolic Problems [DOI: 10.1017/CBO9780511791253](https://doi.org/10.1017/CBO9780511791253) """ -struct ShallowWaterEquations1D{RealT<:Real} <: AbstractShallowWaterEquations{1, 3} - gravity::RealT # gravitational constant - H0::RealT # constant "lake-at-rest" total water height - # `threshold_limiter` used in `PositivityPreservingLimiterShallowWater` on water height, - # as a (small) shift on the initial condition and cutoff before the next time step. - # Default is 500*eps() which in double precision is ≈1e-13. - threshold_limiter::RealT - # `threshold_wet` applied on water height to define when the flow is "wet" - # before calculating the numerical flux. - # Default is 5*eps() which in double precision is ≈1e-15. - threshold_wet::RealT +struct ShallowWaterEquations1D{RealT <: Real} <: AbstractShallowWaterEquations{1, 3} + gravity::RealT # gravitational constant + H0::RealT # constant "lake-at-rest" total water height + # `threshold_limiter` used in `PositivityPreservingLimiterShallowWater` on water height, + # as a (small) shift on the initial condition and cutoff before the next time step. + # Default is 500*eps() which in double precision is ≈1e-13. + threshold_limiter::RealT + # `threshold_wet` applied on water height to define when the flow is "wet" + # before calculating the numerical flux. + # Default is 5*eps() which in double precision is ≈1e-15. + threshold_wet::RealT end # Allow for flexibility to set the gravitational constant within an elixir depending on the @@ -68,16 +68,16 @@ end # The reference total water height H0 defaults to 0.0 but is used for the "lake-at-rest" # well-balancedness test cases. # Strict default values for thresholds that performed well in many numerical experiments -function ShallowWaterEquations1D(; gravity_constant, H0=zero(gravity_constant), - threshold_limiter=nothing, threshold_wet=nothing) - T = promote_type(typeof(gravity_constant), typeof(H0)) - if threshold_limiter === nothing - threshold_limiter = 500 * eps(T) - end - if threshold_wet === nothing - threshold_wet = 5 * eps(T) - end - ShallowWaterEquations1D(gravity_constant, H0, threshold_limiter, threshold_wet) +function ShallowWaterEquations1D(; gravity_constant, H0 = zero(gravity_constant), + threshold_limiter = nothing, threshold_wet = nothing) + T = promote_type(typeof(gravity_constant), typeof(H0)) + if threshold_limiter === nothing + threshold_limiter = 500 * eps(T) + end + if threshold_wet === nothing + threshold_wet = 5 * eps(T) + end + ShallowWaterEquations1D(gravity_constant, H0, threshold_limiter, threshold_wet) end have_nonconservative_terms(::ShallowWaterEquations1D) = True() @@ -351,30 +351,30 @@ Further details on the hydrostatic reconstruction and its motivation can be foun orientation::Integer, equations::ShallowWaterEquations1D) - # Pull the water height and bottom topography on the left - h_ll, _, b_ll = u_ll - h_rr, _, b_rr = u_rr - - H_ll = h_ll + b_ll - H_rr = h_rr + b_rr + # Pull the water height and bottom topography on the left + h_ll, _, b_ll = u_ll + h_rr, _, b_rr = u_rr - b_star = min( max( b_ll, b_rr ), min( H_ll, H_rr ) ) + H_ll = h_ll + b_ll + H_rr = h_rr + b_rr - # Create the hydrostatic reconstruction for the left solution state - u_ll_star, _ = hydrostatic_reconstruction_chen_noelle(u_ll, u_rr, equations) + b_star = min(max(b_ll, b_rr), min(H_ll, H_rr)) - # Copy the reconstructed water height for easier to read code - h_ll_star = u_ll_star[1] + # Create the hydrostatic reconstruction for the left solution state + u_ll_star, _ = hydrostatic_reconstruction_chen_noelle(u_ll, u_rr, equations) - z = zero(eltype(u_ll)) - # Includes two parts: - # (i) Diagonal (consistent) term from the volume flux that uses `b_ll` to avoid - # cross-averaging across a discontinuous bottom topography - # (ii) True surface part that uses `h_ll` and `h_ll_star` to handle discontinuous bathymetry - return SVector(z, - equations.gravity * h_ll * b_ll - equations.gravity * (h_ll_star + h_ll) * (b_ll - b_star), - z) + # Copy the reconstructed water height for easier to read code + h_ll_star = u_ll_star[1] + z = zero(eltype(u_ll)) + # Includes two parts: + # (i) Diagonal (consistent) term from the volume flux that uses `b_ll` to avoid + # cross-averaging across a discontinuous bottom topography + # (ii) True surface part that uses `h_ll` and `h_ll_star` to handle discontinuous bathymetry + return SVector(z, + equations.gravity * h_ll * b_ll - + equations.gravity * (h_ll_star + h_ll) * (b_ll - b_star), + z) end """ @@ -495,48 +495,49 @@ Further details on this hydrostatic reconstruction and its motivation can be fou A new hydrostatic reconstruction scheme based on subcell reconstructions [DOI:10.1137/15M1053074](https://dx.doi.org/10.1137/15M1053074) """ -@inline function hydrostatic_reconstruction_chen_noelle(u_ll, u_rr, equations::ShallowWaterEquations1D) - # Unpack left and right water heights and bottom topographies - h_ll, _, b_ll = u_ll - h_rr, _, b_rr = u_rr - - # Get the velocities on either side - v_ll = velocity(u_ll, equations) - v_rr = velocity(u_rr, equations) - - H_ll = b_ll + h_ll - H_rr = b_rr + h_rr +@inline function hydrostatic_reconstruction_chen_noelle(u_ll, u_rr, + equations::ShallowWaterEquations1D) + # Unpack left and right water heights and bottom topographies + h_ll, _, b_ll = u_ll + h_rr, _, b_rr = u_rr - b_star = min( max( b_ll, b_rr ), min( H_ll, H_rr ) ) + # Get the velocities on either side + v_ll = velocity(u_ll, equations) + v_rr = velocity(u_rr, equations) - # Compute the reconstructed water heights - h_ll_star = min( H_ll - b_star, h_ll ) - h_rr_star = min( H_rr - b_star, h_rr ) + H_ll = b_ll + h_ll + H_rr = b_rr + h_rr - # Set the water height to be at least the value stored in the variable threshold after - # the hydrostatic reconstruction is applied and before the numerical flux is calculated - # to avoid numerical problem with arbitrary small values. Interfaces with a water height - # lower or equal to the threshold can be declared as dry. - # The default value for `threshold_wet` is ≈ 5*eps(), or 1e-15 in double precision, is set - # in the `ShallowWaterEquations1D` struct. This threshold value can be changed in the constructor - # call of this eqaution struct in an elixir. - threshold = equations.threshold_wet + b_star = min(max(b_ll, b_rr), min(H_ll, H_rr)) - if (h_ll_star <= threshold) - h_ll_star = threshold - v_ll = zero(v_ll) - end + # Compute the reconstructed water heights + h_ll_star = min(H_ll - b_star, h_ll) + h_rr_star = min(H_rr - b_star, h_rr) + + # Set the water height to be at least the value stored in the variable threshold after + # the hydrostatic reconstruction is applied and before the numerical flux is calculated + # to avoid numerical problem with arbitrary small values. Interfaces with a water height + # lower or equal to the threshold can be declared as dry. + # The default value for `threshold_wet` is ≈ 5*eps(), or 1e-15 in double precision, is set + # in the `ShallowWaterEquations1D` struct. This threshold value can be changed in the constructor + # call of this eqaution struct in an elixir. + threshold = equations.threshold_wet + + if (h_ll_star <= threshold) + h_ll_star = threshold + v_ll = zero(v_ll) + end - if (h_rr_star <= threshold) - h_rr_star = threshold - v_rr = zero(v_rr) - end + if (h_rr_star <= threshold) + h_rr_star = threshold + v_rr = zero(v_rr) + end - # Create the conservative variables using the reconstruted water heights - u_ll_star = SVector( h_ll_star, h_ll_star * v_ll, b_ll ) - u_rr_star = SVector( h_rr_star, h_rr_star * v_rr, b_rr ) + # Create the conservative variables using the reconstruted water heights + u_ll_star = SVector(h_ll_star, h_ll_star * v_ll, b_ll) + u_rr_star = SVector(h_rr_star, h_rr_star * v_rr, b_rr) - return u_ll_star, u_rr_star + return u_ll_star, u_rr_star end # Calculate maximum wave speed for local Lax-Friedrichs-type dissipation as the @@ -603,7 +604,6 @@ end return λ_min, λ_max end - """ min_max_speed_chen_noelle(u_ll, u_rr, orientation::Integer, equations::ShallowWaterEquations1D) @@ -617,25 +617,25 @@ Further details on this hydrostatic reconstruction and its motivation can be fou A new hydrostatic reconstruction scheme based on subcell reconstructions [DOI:10.1137/15M1053074](https://dx.doi.org/10.1137/15M1053074) """ -@inline function min_max_speed_chen_noelle(u_ll, u_rr, orientation::Integer, equations::ShallowWaterEquations1D) - # Get the velocity quantities - v_ll = velocity(u_ll, equations) - v_rr = velocity(u_rr, equations) +@inline function min_max_speed_chen_noelle(u_ll, u_rr, orientation::Integer, + equations::ShallowWaterEquations1D) + # Get the velocity quantities + v_ll = velocity(u_ll, equations) + v_rr = velocity(u_rr, equations) - # Calculate the wave celerity on the left and right - h_ll = waterheight(u_ll, equations) - h_rr = waterheight(u_rr, equations) + # Calculate the wave celerity on the left and right + h_ll = waterheight(u_ll, equations) + h_rr = waterheight(u_rr, equations) - a_ll = sqrt(equations.gravity * h_ll) - a_rr = sqrt(equations.gravity * h_rr) + a_ll = sqrt(equations.gravity * h_ll) + a_rr = sqrt(equations.gravity * h_rr) - λ_min = min( v_ll - a_ll, v_rr - a_rr, zero(eltype(u_ll)) ) - λ_max = max( v_ll + a_ll, v_rr + a_rr, zero(eltype(u_ll)) ) + λ_min = min(v_ll - a_ll, v_rr - a_rr, zero(eltype(u_ll))) + λ_max = max(v_ll + a_ll, v_rr + a_rr, zero(eltype(u_ll))) - return λ_min, λ_max + return λ_min, λ_max end - @inline function max_abs_speeds(u, equations::ShallowWaterEquations1D) h = waterheight(u, equations) v = velocity(u, equations) @@ -737,14 +737,14 @@ end # be a constant value over time. Note, assumes there is a single reference # water height `H0` with which to compare. @inline function lake_at_rest_error(u, equations::ShallowWaterEquations1D) - h, _, b = u + h, _, b = u - # 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 ) + # 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 - (h + b)) + return abs(H0_wet_dry - (h + b)) end end # @muladd diff --git a/src/solvers/dgsem_tree/indicators_1d.jl b/src/solvers/dgsem_tree/indicators_1d.jl index 12e856ee12..f6c6ea6275 100644 --- a/src/solvers/dgsem_tree/indicators_1d.jl +++ b/src/solvers/dgsem_tree/indicators_1d.jl @@ -27,104 +27,111 @@ end # Modified indicator for ShallowWaterEquations1D to apply full FV method on cells # 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, mesh::Union{TreeMesh{1}, StructuredMesh{1}}, - equations::ShallowWaterEquations1D, dg::DGSEM, +function (indicator_hg::IndicatorHennemannGassnerShallowWater)(u, + mesh::Union{TreeMesh{1}, + StructuredMesh{ + 1 + } + }, + equations::ShallowWaterEquations1D, + dg::DGSEM, cache; kwargs...) - @unpack alpha_max, alpha_min, alpha_smooth, variable = indicator_hg - @unpack alpha, alpha_tmp, indicator_threaded, modal_threaded = indicator_hg.cache - # TODO: Taal refactor, when to `resize!` stuff changed possibly by AMR? - # Shall we implement `resize!(semi::AbstractSemidiscretization, new_size)` - # or just `resize!` whenever we call the relevant methods as we do now? - resize!(alpha, nelements(dg, cache)) - if alpha_smooth - resize!(alpha_tmp, nelements(dg, cache)) - end - - # magic parameters - 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_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 cells and guarantees the well-balanced property. - # - # Hard-coded cut-off value of `threshold_wet = 1e-4` was determined through many numerical experiments. - # Overall idea is to increase robustness when computing the velocity on (nearly) dry cells 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 cells 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 cells 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_wet = 1e-4 - - @threaded for element in eachelement(dg, cache) - indicator = indicator_threaded[Threads.threadid()] - modal = modal_threaded[Threads.threadid()] - - # (Re-)set dummy variable for alpha_dry - indicator_wet = 1 + @unpack alpha_max, alpha_min, alpha_smooth, variable = indicator_hg + @unpack alpha, alpha_tmp, indicator_threaded, modal_threaded = indicator_hg.cache + # TODO: Taal refactor, when to `resize!` stuff changed possibly by AMR? + # Shall we implement `resize!(semi::AbstractSemidiscretization, new_size)` + # or just `resize!` whenever we call the relevant methods as we do now? + resize!(alpha, nelements(dg, cache)) + if alpha_smooth + resize!(alpha_tmp, nelements(dg, cache)) + end - # 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 + # magic parameters + 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_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 cells and guarantees the well-balanced property. + # + # Hard-coded cut-off value of `threshold_wet = 1e-4` was determined through many numerical experiments. + # Overall idea is to increase robustness when computing the velocity on (nearly) dry cells 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 cells 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 cells 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_wet = 1e-4 - if h <= threshold_wet - indicator_wet = 0 - end + @threaded for element in eachelement(dg, cache) + indicator = indicator_threaded[Threads.threadid()] + modal = modal_threaded[Threads.threadid()] - indicator[i] = indicator_hg.variable(u_local, equations) - end + # (Re-)set dummy variable for alpha_dry + indicator_wet = 1 - # Convert to modal representation - multiply_scalar_dimensionwise!(modal, dg.basis.inverse_vandermonde_legendre, indicator) + # 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 - # Calculate total energies for all modes, without highest, without two highest - total_energy = zero(eltype(modal)) - for i in 1:nnodes(dg) - total_energy += modal[i]^2 - end - total_energy_clip1 = zero(eltype(modal)) - for i in 1:(nnodes(dg)-1) - total_energy_clip1 += modal[i]^2 - end - total_energy_clip2 = zero(eltype(modal)) - for i in 1:(nnodes(dg)-2) - total_energy_clip2 += modal[i]^2 - end + if h <= threshold_wet + indicator_wet = 0 + end - # Calculate energy in higher modes - energy = max((total_energy - total_energy_clip1) / total_energy, - (total_energy_clip1 - total_energy_clip2) / total_energy_clip1) + indicator[i] = indicator_hg.variable(u_local, equations) + end - alpha_element = 1 / (1 + exp(-parameter_s / threshold * (energy - threshold))) + # Convert to modal representation + multiply_scalar_dimensionwise!(modal, dg.basis.inverse_vandermonde_legendre, + indicator) - # Take care of the case close to pure DG - if alpha_element < alpha_min - alpha_element = zero(alpha_element) - end + # Calculate total energies for all modes, without highest, without two highest + total_energy = zero(eltype(modal)) + for i in 1:nnodes(dg) + total_energy += modal[i]^2 + end + total_energy_clip1 = zero(eltype(modal)) + for i in 1:(nnodes(dg) - 1) + total_energy_clip1 += modal[i]^2 + end + total_energy_clip2 = zero(eltype(modal)) + for i in 1:(nnodes(dg) - 2) + total_energy_clip2 += modal[i]^2 + end - # Take care of the case close to pure FV - if alpha_element > 1 - alpha_min - alpha_element = one(alpha_element) - end + # Calculate energy in higher modes + energy = max((total_energy - total_energy_clip1) / total_energy, + (total_energy_clip1 - total_energy_clip2) / total_energy_clip1) + + alpha_element = 1 / (1 + exp(-parameter_s / threshold * (energy - threshold))) + + # Take care of the case close to pure DG + if alpha_element < alpha_min + alpha_element = zero(alpha_element) + end - # Clip the maximum amount of FV allowed or set to one depending on indicator_wet - if indicator_wet == 0 - alpha[element] = 1 - else # Element is not defined as dry but wet - alpha[element] = min(alpha_max, alpha_element) + # Take care of the case close to pure FV + if alpha_element > 1 - alpha_min + alpha_element = one(alpha_element) + end + + # Clip the maximum amount of FV allowed or set to one depending on indicator_wet + if indicator_wet == 0 + alpha[element] = 1 + else # Element is not defined as dry but wet + alpha[element] = min(alpha_max, alpha_element) + end end - end - if alpha_smooth - apply_smoothing!(mesh, alpha, alpha_tmp, dg, cache) - end + if alpha_smooth + apply_smoothing!(mesh, alpha, alpha_tmp, dg, cache) + end - return alpha + return alpha end # Use this function barrier and unpack inside to avoid passing closures to Polyester.jl From 0792bd32565204020afad118ac3dbbe357ae0113 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Thu, 6 Jul 2023 08:23:57 +0200 Subject: [PATCH 44/44] simplify elixir as we can set discontinuous ICs in 1D. Also update beach test values --- .../elixir_shallowwater_beach.jl | 8 ++-- ...ixir_shallowwater_well_balanced_wet_dry.jl | 37 ++----------------- src/solvers/dgsem_tree/indicators_1d.jl | 13 +++---- test/test_tree_1d_shallowwater.jl | 4 +- 4 files changed, 15 insertions(+), 47 deletions(-) diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_beach.jl b/examples/tree_1d_dgsem/elixir_shallowwater_beach.jl index d025430b48..c920333bef 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_beach.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_beach.jl @@ -11,7 +11,7 @@ equations = ShallowWaterEquations1D(gravity_constant=9.812) initial_condition_beach(x, t, equations:: ShallowWaterEquations1D) 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 water height and speed function used here, are an adaption of the initial condition +The water height and speed function used here, are an adaption of 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 @@ -34,7 +34,7 @@ function initial_condition_beach(x, t, equations:: ShallowWaterEquations1D) v = 0.0 else H = f - v = sqrt(equations.gravity/D) * H + 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 @@ -69,7 +69,7 @@ volume_integral = VolumeIntegralShockCapturingHG(indicator_sc; volume_flux_fv=surface_flux) solver = DGSEM(basis, surface_flux, volume_integral) - + ############################################################################### # Create the TreeMesh for the domain [0, 8] @@ -101,7 +101,7 @@ analysis_callback = AnalysisCallback(semi, interval=analysis_interval, save_anal alive_callback = AliveCallback(analysis_interval=analysis_interval) -save_solution = SaveSolutionCallback(interval=250, +save_solution = SaveSolutionCallback(dt=0.5, save_initial_solution=true, save_final_solution=true) diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl b/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl index ac108a8c69..c68c54e5ad 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl @@ -65,7 +65,7 @@ volume_integral = VolumeIntegralShockCapturingHG(indicator_sc; volume_flux_fv=surface_flux) solver = DGSEM(basis, surface_flux, volume_integral) - + ############################################################################### # Create the TreeMesh for the domain [0, 1] @@ -85,38 +85,9 @@ semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) tspan = (0.0, 25.0) ode = semidiscretize(semi, tspan) - -############################################################################### -# Workaround to set a discontinuous water and bottom topography for -# debugging and testing. Essentially, this is a slight augmentation of the -# `compute_coefficients` where the `x` node value passed here is slightly -# perturbed to the left / right in order to set a true discontinuity that avoids -# the doubled value of the LGL nodes at a particular element interface. -# -# Note! The errors from the analysis callback are not important but the error -# for this lake at rest test case `∑|H0-(h+b)|` should be near machine roundoff. - -# point to the data we want to augment -u = Trixi.wrap_array(ode.u0, semi) -# reset the initial condition -for element in eachelement(semi.solver, semi.cache) - for i in eachnode(semi.solver) - x_node = Trixi.get_node_coords(semi.cache.elements.node_coordinates, equations, semi.solver, i, element) - # We know that the discontinuity is on an interface. Slightly augment the x value by a factor - # of unit roundoff to avoid the repeated value from the LGL nodes at the interface. - if i == 1 - x_node = SVector(nextfloat(x_node[1])) - elseif i == nnodes(semi.solver) - x_node = SVector(prevfloat(x_node[1])) - end - u_node = initial_condition_complex_bottom_well_balanced(x_node, first(tspan), equations) - Trixi.set_node_vars!(u, u_node, equations, semi.solver, i, element) - end -end - summary_callback = SummaryCallback() -analysis_interval = 1000 +analysis_interval = 5000 analysis_callback = AnalysisCallback(semi, interval=analysis_interval, save_analysis=false) alive_callback = AliveCallback(analysis_interval=analysis_interval) @@ -125,7 +96,7 @@ save_solution = SaveSolutionCallback(interval=5000, save_initial_solution=true, save_final_solution=true) -stepsize_callback = StepsizeCallback(cfl=1.5) +stepsize_callback = StepsizeCallback(cfl=1.5) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, stepsize_callback) @@ -141,7 +112,7 @@ sol = solve(ode, SSPRK43(stage_limiter!); dt=1.0, summary_callback() # print the timer summary ############################################################################### -# workaround to compute the well-balancedness error for this partiular problem +# Workaround to compute the well-balancedness error for this particular problem # that has two reference water heights. One for a lake to the left of the # discontinuous bottom topography `H0_upper = 2.5` and another for a lake to the # right of the discontinuous bottom topography `H0_lower = 1.5`. diff --git a/src/solvers/dgsem_tree/indicators_1d.jl b/src/solvers/dgsem_tree/indicators_1d.jl index f6c6ea6275..04252b68bd 100644 --- a/src/solvers/dgsem_tree/indicators_1d.jl +++ b/src/solvers/dgsem_tree/indicators_1d.jl @@ -27,15 +27,12 @@ end # Modified indicator for ShallowWaterEquations1D to apply full FV method on cells # 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, - mesh::Union{TreeMesh{1}, - StructuredMesh{ - 1 - } - }, +function (indicator_hg::IndicatorHennemannGassnerShallowWater)(u::AbstractArray{<:Any, 3 + }, + mesh, equations::ShallowWaterEquations1D, - dg::DGSEM, - cache; kwargs...) + dg::DGSEM, cache; + kwargs...) @unpack alpha_max, alpha_min, alpha_smooth, variable = indicator_hg @unpack alpha, alpha_tmp, indicator_threaded, modal_threaded = indicator_hg.cache # TODO: Taal refactor, when to `resize!` stuff changed possibly by AMR? diff --git a/test/test_tree_1d_shallowwater.jl b/test/test_tree_1d_shallowwater.jl index e1201d4e2f..573d02f363 100644 --- a/test/test_tree_1d_shallowwater.jl +++ b/test/test_tree_1d_shallowwater.jl @@ -98,8 +98,8 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem") @trixi_testset "elixir_shallowwater_beach.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_beach.jl"), - l2 = [0.17978981087947196, 1.237741563904404, 6.289469352788183e-8], - linf = [0.8460686638571695, 3.374190220014994, 4.452604152049844e-7], + l2 = [0.17979210479598923, 1.2377495706611434, 6.289818963361573e-8], + linf = [0.845938394800688, 3.3740800777086575, 4.4541473087633676e-7], tspan = (0.0, 0.05)) end