From ee96ce99e11d1990a94d556969c321c76d8a3c00 Mon Sep 17 00:00:00 2001 From: patrickersing Date: Tue, 9 Apr 2024 11:21:18 +0200 Subject: [PATCH 01/17] add wet/dry functionality for the 1d mlswe --- src/TrixiShallowWater.jl | 2 +- .../positivity_shallow_water.jl | 26 +++- .../positivity_shallow_water_dg1d.jl | 87 ++++++++++++ .../positivity_shallow_water_dg2d.jl | 2 + src/equations/shallow_water_multilayer_1d.jl | 134 ++++++++++++++++-- src/solvers/indicators.jl | 3 +- src/solvers/indicators_1d.jl | 44 +++--- 7 files changed, 260 insertions(+), 38 deletions(-) diff --git a/src/TrixiShallowWater.jl b/src/TrixiShallowWater.jl index ac74b97..eef2ece 100644 --- a/src/TrixiShallowWater.jl +++ b/src/TrixiShallowWater.jl @@ -24,7 +24,7 @@ export ShallowWaterEquationsWetDry1D, ShallowWaterEquationsWetDry2D, export hydrostatic_reconstruction_chen_noelle, flux_nonconservative_chen_noelle, min_max_speed_chen_noelle, flux_hll_chen_noelle, - flux_ersing_etal, flux_es_ersing_etal + flux_ersing_etal, flux_es_ersing_etal, hydrostatic_reconstruction_ersing_etal export nlayers, eachlayer diff --git a/src/callbacks_stage/positivity_shallow_water.jl b/src/callbacks_stage/positivity_shallow_water.jl index 17b8388..62b03a4 100644 --- a/src/callbacks_stage/positivity_shallow_water.jl +++ b/src/callbacks_stage/positivity_shallow_water.jl @@ -8,12 +8,16 @@ """ PositivityPreservingLimiterShallowWater(; variables) +!!! warning "Experimental code" + This is an experimental feature and may change in future releases. + The limiter is specifically designed for the shallow water equations. It is applied to all scalar `variables` in their given order -using the defined `threshold_limiter` from the [`ShallowWaterEquationsWetDry1D`](@ref) struct -or the [`ShallowWaterEquationsWetDry2D`](@ref) struct to determine the minimal acceptable values. +using the defined `threshold_limiter` from the equations struct +(e.g. in [`ShallowWaterEquationsWetDry1D`](@ref)) to determine the minimal acceptable values. The order of the `variables` is important and might have a strong influence -on the robustness. +on the robustness. The limiter is availble for the [`ShallowWaterEquationsWetDry1D`](@ref), +[`ShallowWaterEquationsWetDry2D`](@ref), and [`ShallowWaterMultiLayerEquations1D`](@ref). 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. @@ -28,11 +32,21 @@ the `threshold_limiter` is applied again on all the DG nodes in order to avoid w In the case where the cell mean value is below the threshold before applying the limiter, there could still be dry nodes afterwards due to the logic of the limiter. +For the [`ShallowWaterMultiLayerEquations1D`](@ref) the implementation differs. In this case the +positivity limiter is applied layerwise and only the waterheight `h` is limited within each layer. +Furthermore, a velocity desingularization is applied after the limiting to avoid numerical problems +near dry states. + This fully-discrete positivity-preserving limiter is based on the work of - Zhang, Shu (2011) Maximum-principle-satisfying and positivity-preserving high-order schemes for conservation laws: survey and new developments [doi: 10.1098/rspa.2011.0153](https://doi.org/10.1098/rspa.2011.0153) +The specific implementation for the [`ShallowWaterMultiLayerEquations1D](@ref) is based on the work of +- Y. Xing, X. Zhang (2013) + Positivity-preserving well-balanced discontinuous Galerkin methods for the shallow water equations + on unstructured triangular meshes + [doi: 10.1007/s10915-012-9644-4](https://doi.org/10.1007/s10915-012-9644-4) """ struct PositivityPreservingLimiterShallowWater{N, Variables <: NTuple{N, Any}} variables::Variables @@ -62,7 +76,8 @@ end function limiter_shallow_water!(u, variables::NTuple{N, Any}, mesh, equations::Union{ShallowWaterEquationsWetDry1D, - ShallowWaterEquationsWetDry2D}, + ShallowWaterEquationsWetDry2D, + ShallowWaterMultiLayerEquations1D}, solver, cache) where {N} variable = first(variables) remaining_variables = Base.tail(variables) @@ -77,7 +92,8 @@ end function limiter_shallow_water!(u, variables::Tuple{}, mesh, equations::Union{ShallowWaterEquationsWetDry1D, - ShallowWaterEquationsWetDry2D}, + ShallowWaterEquationsWetDry2D, + ShallowWaterMultiLayerEquations1D}, solver, cache) nothing end diff --git a/src/callbacks_stage/positivity_shallow_water_dg1d.jl b/src/callbacks_stage/positivity_shallow_water_dg1d.jl index 3494662..9f6bbe7 100644 --- a/src/callbacks_stage/positivity_shallow_water_dg1d.jl +++ b/src/callbacks_stage/positivity_shallow_water_dg1d.jl @@ -5,6 +5,8 @@ @muladd begin #! format: noindent +# !!! warning "Experimental code" +# This is an experimental feature and may change in future releases. function limiter_shallow_water!(u, threshold::Real, variable, mesh::Trixi.AbstractMesh{1}, equations::ShallowWaterEquationsWetDry1D, @@ -84,4 +86,89 @@ function limiter_shallow_water!(u, threshold::Real, variable, return nothing end + +# !!! warning "Experimental code" +# This is an experimental feature and may change in future releases. + +# Note that for the `ShallowWaterMultiLayerEquations1D` only the waterheight `h` is limited in +# each layer. Furthermore, a velocity desingularization is applied after the limiting to avoid +# numerical problems near dry states. +function limiter_shallow_water!(u, threshold::Real, variable, + mesh::Trixi.AbstractMesh{1}, + equations::ShallowWaterMultiLayerEquations1D, + dg::DGSEM, cache) + @unpack weights = dg.basis + + Trixi.@threaded for element in eachelement(dg, cache) + # Limit layerwise + for m in eachlayer(equations) + # determine minimum value + value_min = typemax(eltype(u)) + for i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, element) + value_min = min(value_min, variable(u_node, equations)[m]) + end + + # detect if limiting is necessary + value_min < threshold - eps() || continue + + # compute mean value + u_mean = zero(get_node_vars(u, equations, dg, 1, element)) + for i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, element) + u_mean += u_node * weights[i] + end + + # note that the reference element is [-1,1]^ndims(dg), thus the weights sum to 2 + u_mean = u_mean / 2^ndims(mesh) + + # We compute the value directly with the mean values. + # The waterheight `h` is limited independently in each layer. + value_mean = variable(u_mean, equations)[m] + theta = (value_mean - threshold) / (value_mean - value_min) + + for i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, element) + h_node = waterheight(u_node, equations)[m] + h_mean = waterheight(u_mean, equations)[m] + + u[m, i, element] = theta * h_node + (1 - theta) * h_mean + end + end + end + + # "Safety" application of the wet/dry thresholds over all the DG nodes + # on the current `element` after the limiting above in order to avoid dry nodes. + # If the value_mean < threshold before applying limiter, there + # could still be dry nodes afterwards due to logic of the limiting. + # Additionally, a velocity desingularization is applied to avoid numerical + # problems near dry states. + Trixi.@threaded for element in eachelement(dg, cache) + for i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, element) + + h = MVector(waterheight(u_node, equations)) + hv = MVector(momentum(u_node, equations)) + b = u_node[end] + + # Velocity desingularization + hv = h .* (2.0 * h .* hv) ./ + (h .^ 2 .+ max.(h .^ 2, equations.threshold_desingularization)) + + for i in eachlayer(equations) + # Ensure positivity and zero velocity at dry states + if h[i] <= threshold + h[i] = threshold + hv[i] = zero(eltype(u)) + end + end + + u_node = SVector(h..., hv..., b) + + set_node_vars!(u, u_node, equations, dg, i, element) + end + end + + return nothing +end end # @muladd diff --git a/src/callbacks_stage/positivity_shallow_water_dg2d.jl b/src/callbacks_stage/positivity_shallow_water_dg2d.jl index 812f453..69c47e1 100644 --- a/src/callbacks_stage/positivity_shallow_water_dg2d.jl +++ b/src/callbacks_stage/positivity_shallow_water_dg2d.jl @@ -5,6 +5,8 @@ @muladd begin #! format: noindent +# !!! warning "Experimental code" +# This is an experimental feature and may change in future releases. function limiter_shallow_water!(u, threshold::Real, variable, mesh::Trixi.AbstractMesh{2}, equations::ShallowWaterEquationsWetDry2D, dg::DGSEM, diff --git a/src/equations/shallow_water_multilayer_1d.jl b/src/equations/shallow_water_multilayer_1d.jl index 0dec78a..07acbbc 100644 --- a/src/equations/shallow_water_multilayer_1d.jl +++ b/src/equations/shallow_water_multilayer_1d.jl @@ -30,6 +30,12 @@ nonconservative term, which has some benefits for the design of well-balanced sc The additional quantity ``H_0`` is also available to store a reference value for the total water height that is useful to set initial conditions or test the "lake-at-rest" well-balancedness. +Also, there are two thresholds which prevent numerical problems as well as instabilities. Both of +them do not have to be passed, as default values are defined within the struct. The limiters are +used in [`PositivityPreservingLimiterShallowWater`](@ref) on the water height. `threshold_limiter` +acts as a (small) shift on the initial condition and cutoff before the next time step, whereas +`threshold_desingularization` is used in the velocity desingularization. + The bottom topography function ``b(x)`` is set inside the initial condition routine for a particular problem setup. @@ -55,10 +61,14 @@ struct ShallowWaterMultiLayerEquations1D{NVARS, NLAYERS, RealT <: Real} <: AbstractShallowWaterMultiLayerEquations{1, NVARS, NLAYERS} gravity::RealT # gravitational constant H0::RealT # constant "lake-at-rest" total water height + threshold_limiter::RealT # threshold for the positivity-limiter + threshold_desingularization::RealT # threshold for velocity desingularization rhos::SVector{NLAYERS, RealT} # Vector of layer densities function ShallowWaterMultiLayerEquations1D{NVARS, NLAYERS, RealT}(gravity::RealT, H0::RealT, + threshold_limiter::RealT, + threshold_desingularization::RealT, rhos::SVector{NLAYERS, RealT}) where { NVARS, @@ -71,7 +81,7 @@ struct ShallowWaterMultiLayerEquations1D{NVARS, NLAYERS, RealT <: Real} <: throw(ArgumentError("densities must be in increasing order (rhos[1] < rhos[2] < ... < rhos[NLAYERS])")) min(rhos...) > 0 || throw(ArgumentError("densities must be positive")) - new(gravity, H0, rhos) + new(gravity, H0, threshold_limiter, threshold_desingularization, rhos) end end @@ -80,19 +90,32 @@ end # The reference total water height H0 defaults to 0.0 but is used for the "lake-at-rest" # well-balancedness test cases. function ShallowWaterMultiLayerEquations1D(; gravity_constant, - H0 = zero(gravity_constant), rhos) + H0 = zero(gravity_constant), + threshold_limiter = nothing, + threshold_desingularization = nothing, rhos) # Promote all variables to a common type _rhos = promote(rhos...) RealT = promote_type(eltype(_rhos), eltype(gravity_constant), eltype(H0)) __rhos = SVector(map(RealT, _rhos)) + # Set default values for thresholds + if threshold_limiter === nothing + threshold_limiter = 5 * eps(RealT) + end + if threshold_desingularization === nothing + threshold_desingularization = 1e-10 + end + # Extract number of layers and variables NLAYERS = length(rhos) NVARS = 2 * NLAYERS + 1 return ShallowWaterMultiLayerEquations1D{NVARS, NLAYERS, RealT}(gravity_constant, - H0, __rhos) + H0, + threshold_limiter, + threshold_desingularization, + __rhos) end @inline function Base.real(::ShallowWaterMultiLayerEquations1D{NVARS, NLAYERS, RealT}) where { @@ -354,6 +377,79 @@ In the two-layer setting this combination is equivalent to the fluxes in: return SVector(f) end +""" + hydrostatic_reconstruction_ersing_etal(u_ll, u_rr, equations::ShallowWaterMultiLayerEquations1D) + +!!! warning "Experimental code" + This is an experimental feature and may change in future releases. + +A particular type of hydrostatic reconstruction of the water height and bottom topography to +guarantee well-balancedness in the presence of wet/dry transitions and entropy stability for the +[`ShallowWaterMultiLayerEquations1D`](@ref). +The reconstructed solution states `u_ll_star` and `u_rr_star` variables are used to evaluate the +surface numerical flux at the interface. The key idea is a piecewise linear reconstruction of the +bottom topography and water height interfaces using subcells, where the bottom topography is allowed +to be discontinuous. +Use in combination with the generic numerical flux routine [`Trixi.FluxHydrostaticReconstruction`](@extref). +""" +@inline function hydrostatic_reconstruction_ersing_etal(u_ll, u_rr, + equations::ShallowWaterMultiLayerEquations1D) + # Unpack waterheight and bottom topographies + h_ll = waterheight(u_ll, equations) + h_rr = waterheight(u_rr, equations) + b_ll = u_ll[end] + b_rr = u_rr[end] + + # Get the velocities on either side + v_ll = MVector(velocity(u_ll, equations)) + v_rr = MVector(velocity(u_rr, equations)) + + threshold = equations.threshold_limiter + + # Ensure zero velocity at dry states + for i in eachlayer(equations) + if h_ll[i] <= threshold + v_ll[i] = zero(eltype(u_ll)) + end + if h_rr[i] <= threshold + v_rr[i] = zero(eltype(u_rr)) + end + end + + # Calculate interface heights + H_ll = waterheight(cons2prim(u_ll, equations), equations) + H_rr = waterheight(cons2prim(u_rr, equations), equations) + + # Discontinuous reconstruction of the bottom topography + b_ll_star = min(H_ll[1], max(b_rr, b_ll)) + b_rr_star = min(H_rr[1], max(b_rr, b_ll)) + + # Calculate reconstructed interface heights + H_ll_star = max.(H_ll, b_ll_star) + H_rr_star = max.(H_rr, b_rr_star) + + # Initialize reconstructed waterheights + h_ll_star = zero(MVector{nlayers(equations), real(equations)}) + h_rr_star = zero(MVector{nlayers(equations), real(equations)}) + + # Reconstruct the waterheights + for i in eachlayer(equations) + if i == nlayers(equations) # The lowest layer is measured from the bottom topography + h_ll_star[i] = max(H_ll_star[i] - b_ll_star, threshold) + h_rr_star[i] = max(H_rr_star[i] - b_rr_star, threshold) + else + h_ll_star[i] = max(H_ll_star[i] - H_ll_star[i + 1], threshold) + h_rr_star[i] = max(H_rr_star[i] - H_rr_star[i + 1], threshold) + end + end + + # Reconstructed states + u_ll_star = SVector(h_ll_star..., (h_ll_star .* v_ll)..., b_ll_star) + u_rr_star = SVector(h_rr_star..., (h_rr_star .* v_rr)..., b_rr_star) + + return SVector(u_ll_star, u_rr_star) +end + # Specialized `DissipationLocalLaxFriedrichs` to avoid spurious dissipation in the bottom # topography @inline function (dissipation::DissipationLocalLaxFriedrichs)(u_ll, u_rr, @@ -397,7 +493,7 @@ end c_ll = sqrt(equations.gravity * sum(h_ll)) c_rr = sqrt(equations.gravity * sum(h_rr)) - return (max(abs(v_m_ll) + c_ll, abs(v_m_rr) + c_rr)) + return (max(abs(v_m_ll), abs(v_m_rr)) + max(c_ll, c_rr)) end # Convert conservative variables to primitive @@ -422,18 +518,20 @@ end # Convert primitive to conservative variables @inline function Trixi.prim2cons(prim, equations::ShallowWaterMultiLayerEquations1D) - H = prim[1:nlayers(equations)] - v = prim[(nlayers(equations) + 1):(2 * nlayers(equations))] + # To extract the interface height and velocity we reuse the waterheight and momentum functions + # from the conservative variables. + H = waterheight(prim, equations) # For primitive variables this extracts the interface height + v = momentum(prim, equations) # For primitive variables this extracts the velocity b = prim[end] # Calculate waterheight h = MVector{nlayers(equations), real(equations)}(undef) for i in eachlayer(equations) if i < nlayers(equations) - setindex!(h, H[i] - H[i + 1], i) + h[i] = H[i] - H[i + 1] else # The lowest layer is measured from the bottom topography - setindex!(h, H[i] - b, i) + h[i] = H[i] - b end end @@ -535,14 +633,28 @@ end return energy_total(u, equations) - energy_kinetic(u, equations) end -# Calculate the error for the "lake-at-rest" test case where H = ∑h+b should -# be a constant value over time +@inline function Trixi.waterheight_pressure(u, + equations::ShallowWaterMultiLayerEquations1D) + h = waterheight(u, equations) + + return 0.5 * equations.gravity * sum(h.^3) +end + +# Calculate the error for the "lake-at-rest" test case where H = ∑h + b should +# be a constant value over time. +# Note, assumes there is a single reference water height `H0` with which to compare. @inline function Trixi.lake_at_rest_error(u, equations::ShallowWaterMultiLayerEquations1D) h = waterheight(u, equations) b = u[end] - return abs(equations.H0 - (sum(h) + b)) + # For well-balancedness testing with possible wet/dry regions the reference + # water height `H0` accounts for the possibility that the bottom topography + # can emerge out of the water as well as for the threshold offset to avoid + # division by a "hard" zero water heights as well. + H0_wet_dry = max(equations.H0, b + equations.threshold_limiter) + + return abs(H0_wet_dry - (sum(h) + b)) end # Helper function to set the layer values in the flux computation diff --git a/src/solvers/indicators.jl b/src/solvers/indicators.jl index 0ac2a5a..a650546 100644 --- a/src/solvers/indicators.jl +++ b/src/solvers/indicators.jl @@ -40,7 +40,8 @@ end # of the shallow water equations # It modifies the shock-capturing indicator to use full FV method in dry elements # or partially dry elements containing a wet/dry transition. -function IndicatorHennemannGassnerShallowWater(equations::Trixi.AbstractShallowWaterEquations, +function IndicatorHennemannGassnerShallowWater(equations::Union{Trixi.AbstractShallowWaterEquations, + AbstractShallowWaterMultiLayerEquations}, basis; alpha_max = 0.5, alpha_min = 0.001, diff --git a/src/solvers/indicators_1d.jl b/src/solvers/indicators_1d.jl index 8a9f47d..be9c5cd 100644 --- a/src/solvers/indicators_1d.jl +++ b/src/solvers/indicators_1d.jl @@ -5,13 +5,14 @@ @muladd begin #! format: noindent -# Modified indicator for ShallowWaterEquationsWetDry1D to apply full FV method on elements -# containing some "dry" LGL nodes. That is, if an element is partially "wet" then it becomes a -# full FV element. +# Modified indicator for ShallowWaterEquationsWetDry1D and ShallowWaterMultiLayerEquations1D to +# apply full FV method on elements containing some "dry" LGL nodes. That is, if an element is +# partially "wet" then it becomes a full FV element. function (indicator_hg::IndicatorHennemannGassnerShallowWater)(u::AbstractArray{<:Any, 3}, mesh, - equations::ShallowWaterEquationsWetDry1D, + equations::Union{ShallowWaterEquationsWetDry1D, + ShallowWaterMultiLayerEquations1D}, dg::DGSEM, cache; kwargs...) @unpack alpha_max, alpha_min, alpha_smooth, variable = indicator_hg @@ -28,20 +29,22 @@ function (indicator_hg::IndicatorHennemannGassnerShallowWater)(u::AbstractArray{ threshold = 0.5 * 10^(-1.8 * (nnodes(dg))^0.25) parameter_s = log((1 - 0.0001) / 0.0001) - # If the water height `h` at one LGL node is lower than `threshold_partially_wet` - # the indicator sets the element-wise blending factor alpha[element] = 1 - # via the local variable `indicator_wet`. In turn, this ensures that a pure - # FV method is used in partially wet elements and guarantees the well-balanced property. - # - # Hard-coded cut-off value of `threshold_partially_wet = 1e-4` was determined through many numerical experiments. - # Overall idea is to increase robustness when computing the velocity on (nearly) dry elements which - # could be "dangerous" due to division of conservative variables, e.g., v = hv / h. - # Here, the impact of the threshold on the number of elements being updated with FV is not that - # significant. However, its impact on the robustness is very significant. - # The value can be seen as a trade-off between accuracy and stability. - # Well-balancedness of the scheme on partially wet elements with hydrostatic reconstruction - # can only be proven for the FV method (see Chen and Noelle). - # Therefore we set alpha to one regardless of its given maximum value. + #= + If the water height `h` at one LGL node is lower than `threshold_partially_wet` + the indicator sets the element-wise blending factor alpha[element] = 1 + via the local variable `indicator_wet`. In turn, this ensures that a pure + FV method is used in partially wet elements and guarantees the well-balanced property. + + Hard-coded cut-off value of `threshold_partially_wet = 1e-4` was determined through many numerical experiments. + Overall idea is to increase robustness when computing the velocity on (nearly) dry elements which + could be "dangerous" due to division of conservative variables, e.g., v = hv / h. + Here, the impact of the threshold on the number of elements being updated with FV is not that + significant. However, its impact on the robustness is very significant. + The value can be seen as a trade-off between accuracy and stability. + Well-balancedness of the scheme on partially wet elements with hydrostatic reconstruction + can only be proven for the FV method (see Chen and Noelle). + Therefore we set alpha to one regardless of its given maximum value. + =# threshold_partially_wet = 1e-4 Trixi.@threaded for element in eachelement(dg, cache) @@ -54,9 +57,10 @@ function (indicator_hg::IndicatorHennemannGassnerShallowWater)(u::AbstractArray{ # Calculate indicator variables at Gauss-Lobatto nodes for i in eachnode(dg) u_local = get_node_vars(u, equations, dg, i, element) - h, _, _ = u_local + h = waterheight(u_local, equations) - if h <= threshold_partially_wet + # Set indicator to FV if water height is below the threshold + if minimum(h) <= threshold_partially_wet indicator_wet = 0 end From 07a51bad4438a988d66424b6a8778c4921c9e47d Mon Sep 17 00:00:00 2001 From: patrickersing Date: Tue, 9 Apr 2024 11:25:05 +0200 Subject: [PATCH 02/17] add new elixirs to test wet/dry functionality --- .../elixir_shallowwater_multilayer_beach.jl | 125 +++++++++++ ..._shallowwater_multilayer_parabolic_bowl.jl | 122 +++++++++++ ...wwater_multilayer_well_balanced_wet_dry.jl | 194 ++++++++++++++++++ 3 files changed, 441 insertions(+) create mode 100644 examples/tree_1d_dgsem/elixir_shallowwater_multilayer_beach.jl create mode 100644 examples/tree_1d_dgsem/elixir_shallowwater_multilayer_parabolic_bowl.jl create mode 100644 examples/tree_1d_dgsem/elixir_shallowwater_multilayer_well_balanced_wet_dry.jl diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_multilayer_beach.jl b/examples/tree_1d_dgsem/elixir_shallowwater_multilayer_beach.jl new file mode 100644 index 0000000..91767fa --- /dev/null +++ b/examples/tree_1d_dgsem/elixir_shallowwater_multilayer_beach.jl @@ -0,0 +1,125 @@ + +using OrdinaryDiffEq +using Trixi +using TrixiShallowWater + +############################################################################### +# Semidiscretization of the shallow water equations + +equations = ShallowWaterMultiLayerEquations1D(gravity_constant = 9.812, rhos = 1.0) + +""" + initial_condition_beach(x, t, equations:: ShallowWaterMultiLayerEquationsWetDry1D) + +Initial condition to simulate a wave running towards a beach and crashing. Difficult test +including both wetting and drying in the domain using slip wall boundary conditions. +The bottom topography is altered to be differentiable on the domain [0,8] and +differs from the reference below. + +The water height and speed functions used here, are adapted from the initial condition +found in section 5.2 of the paper: + - Andreas Bollermann, Sebastian Noelle, Maria Lukáčová-Medvidová (2011) + Finite volume evolution Galerkin methods for the shallow water equations with dry beds\n + [DOI: 10.4208/cicp.220210.020710a](https://dx.doi.org/10.4208/cicp.220210.020710a) +""" +function initial_condition_beach(x, t, equations::ShallowWaterMultiLayerEquations1D) + D = 1 + delta = 0.02 + gamma = sqrt((3 * delta) / (4 * D)) + x_a = sqrt((4 * D) / (3 * delta)) * acosh(sqrt(20)) + + f = D + 40 * delta * sech(gamma * (8 * x[1] - x_a))^2 + + # steep curved beach + b = 0.01 + 99 / 409600 * 4^x[1] + + if x[1] >= 6 + H = b + v = 0.0 + else + H = f + v = sqrt(equations.gravity / D) * H + end + + # It is mandatory to shift the water level at dry areas to make sure the water height h + # stays positive. The system would not be stable for h set to a hard 0 due to division by h in + # the computation of velocity, e.g., (h v) / h. Therefore, a small dry state threshold + # with a default value of 5*eps() ≈ 1e-15 in double precision, is set in the constructor above + # for the ShallowWaterMultiLayerEquations1D and added to the initial condition if h = 0. + # This default value can be changed within the constructor call depending on the simulation setup. + H = max(H, b + equations.threshold_limiter) + return prim2cons(SVector(H, v, b), equations) +end + +initial_condition = initial_condition_beach +boundary_condition = boundary_condition_slip_wall + +############################################################################### +# Get the DG approximation space + +volume_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal) +surface_flux = (FluxHydrostaticReconstruction(FluxPlusDissipation(flux_ersing_etal, + DissipationLocalLaxFriedrichs()), + hydrostatic_reconstruction_ersing_etal), + FluxHydrostaticReconstruction(flux_nonconservative_ersing_etal, + hydrostatic_reconstruction_ersing_etal)) + +basis = LobattoLegendreBasis(3) + +indicator_sc = IndicatorHennemannGassnerShallowWater(equations, basis, + alpha_max = 0.5, + alpha_min = 0.001, + alpha_smooth = true, + variable = waterheight_pressure) +volume_integral = VolumeIntegralShockCapturingHG(indicator_sc; + volume_flux_dg = volume_flux, + volume_flux_fv = surface_flux) + +solver = DGSEM(basis, surface_flux, volume_integral) + +############################################################################### +# Create the TreeMesh for the domain [0, 8] + +coordinates_min = 0.0 +coordinates_max = 8.0 + +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 7, + n_cells_max = 10_000, + periodicity = false) + +# create the semi discretization object +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_condition) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 10.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 1000 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + save_analysis = false, + extra_analysis_integrals = (energy_kinetic, + energy_internal)) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(dt = 0.5, + save_initial_solution = true, + save_final_solution = true) + +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution) + +stage_limiter! = PositivityPreservingLimiterShallowWater(variables = (Trixi.waterheight,)) + +############################################################################### +# run the simulation + +sol = solve(ode, SSPRK43(stage_limiter!); + ode_default_options()..., callback = callbacks); + +summary_callback() # print the timer summary diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_multilayer_parabolic_bowl.jl b/examples/tree_1d_dgsem/elixir_shallowwater_multilayer_parabolic_bowl.jl new file mode 100644 index 0000000..8257775 --- /dev/null +++ b/examples/tree_1d_dgsem/elixir_shallowwater_multilayer_parabolic_bowl.jl @@ -0,0 +1,122 @@ + +using OrdinaryDiffEq +using Trixi +using TrixiShallowWater + +############################################################################### +# Semidiscretization of the shallow water equations + +equations = ShallowWaterMultiLayerEquations1D(gravity_constant = 9.81, rhos = 1.0) + +""" + initial_condition_parabolic_bowl(x, t, equations:: ShallowWaterMultiLayerEquations1D) + +Well-known initial condition to test the [`hydrostatic_reconstruction_ersing_etal`](@ref) and its +wet-dry mechanics. This test has analytical solutions. The initial condition is defined by the +analytical solution at time t=0. The bottom topography defines a bowl and the water level is given +by an oscillating lake. + +The original test and its analytical solution in two dimensions were first presented in +- William C. Thacker (1981) + Some exact solutions to the nonlinear shallow-water wave equations + [DOI: 10.1017/S0022112081001882](https://doi.org/10.1017/S0022112081001882). + +The particular setup below is taken from Section 6.2 of +- Niklas Wintermeyer, Andrew R. Winters, Gregor J. Gassner and Timothy Warburton (2018) + An entropy stable discontinuous Galerkin method for the shallow water equations on + curvilinear meshes with wet/dry fronts accelerated by GPUs + [DOI: 10.1016/j.jcp.2018.08.038](https://doi.org/10.1016/j.jcp.2018.08.038). +""" +function initial_condition_parabolic_bowl(x, t, equations::ShallowWaterMultiLayerEquations1D) + a = 1 + h_0 = 0.1 + sigma = 0.5 + ω = sqrt(2 * equations.gravity * h_0) / a + + v = -sigma * ω * sin(ω * t) + + b = h_0 * x[1]^2 / a^2 + + H = sigma * h_0 / a^2 * (2 * x[1] * cos(ω * t) - sigma) + h_0 + + #= + It is mandatory to shift the water level at dry areas to make sure the water height h + stays positive. The system would not be stable for h set to a hard 0 due to division by h in + the computation of velocity, e.g., (h v) / h. Therefore, a small dry state threshold + with a default value of 5*eps() ≈ 1e-15 in double precision, is set in the constructor above + for the ShallowWaterMultiLayerEquations1D and added to the initial condition if h = 0. + This default value can be changed within the constructor call depending on the simulation setup. + =# + H = max(H, b + equations.threshold_limiter) + return prim2cons(SVector(H, v, b), equations) +end + +initial_condition = initial_condition_parabolic_bowl + +############################################################################### +# Get the DG approximation space + +volume_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal) +surface_flux = (FluxHydrostaticReconstruction(FluxPlusDissipation(flux_ersing_etal, + DissipationLocalLaxFriedrichs()), + hydrostatic_reconstruction_ersing_etal), + FluxHydrostaticReconstruction(flux_nonconservative_ersing_etal, + hydrostatic_reconstruction_ersing_etal)) + +basis = LobattoLegendreBasis(5) + +indicator_sc = IndicatorHennemannGassnerShallowWater(equations, basis, + alpha_max = 0.5, + alpha_min = 0.001, + alpha_smooth = true, + variable = waterheight_pressure) +volume_integral = VolumeIntegralShockCapturingHG(indicator_sc; + volume_flux_dg = volume_flux, + volume_flux_fv = surface_flux) + +solver = DGSEM(basis, surface_flux, volume_integral) + +############################################################################### +# Create the TreeMesh for the domain [-2, 2] + +coordinates_min = -2.0 +coordinates_max = 2.0 + +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 6, + n_cells_max = 10_000) + +# create the semi discretization object +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 10.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 1000 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + save_analysis = false, + extra_analysis_integrals = (energy_kinetic, + energy_internal)) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 1000, + save_initial_solution = true, + save_final_solution = true) + +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution) + +stage_limiter! = PositivityPreservingLimiterShallowWater(variables = (Trixi.waterheight,)) + +############################################################################### +# run the simulation + +sol = solve(ode, SSPRK43(stage_limiter!); + ode_default_options()..., callback = callbacks); + +summary_callback() # print the timer summary diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_multilayer_well_balanced_wet_dry.jl b/examples/tree_1d_dgsem/elixir_shallowwater_multilayer_well_balanced_wet_dry.jl new file mode 100644 index 0000000..c92c50e --- /dev/null +++ b/examples/tree_1d_dgsem/elixir_shallowwater_multilayer_well_balanced_wet_dry.jl @@ -0,0 +1,194 @@ + +using OrdinaryDiffEq +using Trixi +using TrixiShallowWater + +############################################################################### +# Semidiscretization of the two-layer shallow water equations to test well-balancedness + +equations = ShallowWaterMultiLayerEquations1D(gravity_constant = 1.0, H0 = 2.0, rhos = (1.0, 3.0)) + +""" + initial_condition_twolayer_well_balanced_wet_dry(perturbation, equations:: ShallowWaterMultiLayerEquations1D) + initial_condition_twolayer_well_balanced_wet_dry(x, t, equations:: ShallowWaterMultiLayerEquations1D) + +Initial condition with a complex (discontinuous) bottom topography to test well-balancedness for a +two-layer shallow water system with dry states if `perturbation` is set to `false`. +Additionally, it is possible to set a perturbation in the lower layer to test perturbations from the +lake-at-rest condition by setting the `perturbation` variable to `true`. + +The initial condition is taken from the paper: + - S. Martínez-Aranda, A. Ramos-Pérez, P. García-Navarro (2020) + A 1D shallow-flow model for two-layer flows based on FORCE scheme with wet–dry treatment\n + [DOI:10.2166/hydro.2020.002](https://doi.org/10.2166/hydro.2020.002) +""" +function initial_condition_twolayer_well_balanced_wet_dry(perturbation::Bool, equations::ShallowWaterMultiLayerEquations1D) + return function initial_condition_twolayer_well_balanced_wet_dry(x, t, equations) + # Set interface height for the upper layer + H_upper = 2.0 + + # Set interface height for the lower layer + H_lower = 1.0 + if perturbation == true + if x[1] <= 5.0 + nothing + elseif x[1] <= 15.0 + H_lower = 0.75 + elseif x[1] <= 35.0 + H_lower = 1.5 - (x[1] - 15.0) / 20.0 * 0.5 + elseif x[1] <= 40.0 + H_lower = 1.0 + (x[1] - 35.0) / 5.0 * 0.5 + elseif x[1] <= 55.0 + H_lower = 1.0 + elseif x[1] <= 70.0 + H_lower = 1.6 - (x[1] - 55.0) / 15.0 * 0.4 + elseif x[1] <= 75.0 + nothing + elseif x[1] <= 85.0 + H_lower = 1.0 - (x[1] - 75.0) / 10.0 * 0.3 + elseif x[1] <= 95.0 + H_lower = 1.1 + else + H_lower = 1.5 + end + end + + # Set the bottom topography + if x[1] <= 5.0 + b = 2.5 + elseif x[1] <= 20.0 + b = 0.5 + elseif x[1] <= 40.0 + b = 0.1 + (x[1] - 20.0) / 20.0 * 1.4 + elseif x[1] <= 60.0 + b = 1.0 - (x[1] - 40.0) / 20.0 * 1.3 + elseif x[1] <= 70.0 + b = 0.7 + elseif x[1] <= 75.0 + b = 2.2 + elseif x[1] <= 95.0 + b = 0.4 + else + b = 1.5 + end + + # Set zero initial velocity + v1_upper = 0.0 + v1_lower = 0.0 + + #= + It is mandatory to shift the water level at dry areas to make sure the water height h + stays positive. The system would not be stable for h set to a hard 0 due to division by h in + the computation of velocity, e.g., (h v) / h. Therefore, a small dry state threshold + with a default value of 5*eps() ≈ 1e-15 in double precision, is set in the constructor above + for the ShallowWaterMultiLayerEquations1D and added to the initial condition if h = 0. + This default value can be changed within the constructor call depending on the simulation setup. + =# + H_lower = max(H_lower, b + equations.threshold_limiter) + H_upper = max(H_upper, H_lower + equations.threshold_limiter) + + return prim2cons(SVector(H_upper, H_lower, v1_upper, v1_lower, b), equations) + end +end + + +initial_condition = initial_condition_twolayer_well_balanced_wet_dry(true, equations) + +############################################################################### +# Get the DG approximation space + +volume_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal) +surface_flux = (FluxHydrostaticReconstruction(FluxPlusDissipation(flux_ersing_etal, + DissipationLocalLaxFriedrichs()), + hydrostatic_reconstruction_ersing_etal), + FluxHydrostaticReconstruction(flux_nonconservative_ersing_etal, + hydrostatic_reconstruction_ersing_etal)) +basis = LobattoLegendreBasis(3) + +indicator_sc = IndicatorHennemannGassnerShallowWater(equations, basis, + alpha_max = 0.5, + alpha_min = 0.001, + alpha_smooth = true, + variable = waterheight_pressure) +volume_integral = VolumeIntegralShockCapturingHG(indicator_sc; + volume_flux_dg = volume_flux, + volume_flux_fv = surface_flux) +solver = DGSEM(basis, surface_flux, volume_integral) +############################################################################### +# Get the TreeMesh and setup a periodic mesh + +# The domain of interest [0.0, 100.0] is extended towards -28.0 to match the positions of the +# disctontinuities with TreeMesh. +coordinates_min = -28.0 +coordinates_max = 100.0 +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 7, + n_cells_max = 10_000, + periodicity = true) + +# create the semi discretization object +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 10.0) +ode = semidiscretize(semi, tspan) + +################################################################################# +# 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 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])) + elseif i == nnodes(semi.solver) + x_node = SVector(prevfloat(x_node[1])) + end + u_node = initial_condition(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_callback = AnalysisCallback(semi, interval = analysis_interval, + save_analysis = false, + extra_analysis_integrals = (lake_at_rest_error,)) + +stepsize_callback = StepsizeCallback(cfl = 0.5) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 1000, + save_initial_solution = true, + save_final_solution = true) + +stage_limiter! = PositivityPreservingLimiterShallowWater(variables = (Trixi.waterheight,)) + +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, + stepsize_callback) + +############################################################################### +# run the simulation + +# use a Runge-Kutta method with automatic (error based) time step size control +sol = solve(ode, SSPRK43(stage_limiter!); + ode_default_options()..., callback = callbacks, adaptive = false, dt = 1.0); +summary_callback() # print the timer summary From 1965d9b7b5aaba855402f8ac666ead47ff4ed9e2 Mon Sep 17 00:00:00 2001 From: patrickersing Date: Tue, 9 Apr 2024 11:42:47 +0200 Subject: [PATCH 03/17] add tests --- ...wwater_multilayer_well_balanced_wet_dry.jl | 4 +- test/test_tree_1d.jl | 142 ++++++++++++++++++ test/test_upstream.jl | 1 - 3 files changed, 145 insertions(+), 2 deletions(-) diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_multilayer_well_balanced_wet_dry.jl b/examples/tree_1d_dgsem/elixir_shallowwater_multilayer_well_balanced_wet_dry.jl index c92c50e..1e271e0 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_multilayer_well_balanced_wet_dry.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_multilayer_well_balanced_wet_dry.jl @@ -91,8 +91,10 @@ function initial_condition_twolayer_well_balanced_wet_dry(perturbation::Bool, eq end end +# Default value for the perturbation +perturbation = false -initial_condition = initial_condition_twolayer_well_balanced_wet_dry(true, equations) +initial_condition = initial_condition_twolayer_well_balanced_wet_dry(perturbation, equations) ############################################################################### # Get the DG approximation space diff --git a/test/test_tree_1d.jl b/test/test_tree_1d.jl index 14103c7..6e6e982 100644 --- a/test/test_tree_1d.jl +++ b/test/test_tree_1d.jl @@ -542,6 +542,43 @@ end # 2LSWE end end + @trixi_testset "elixir_shallowwater_multilayer_convergence.jl with hydrostatic_reconstruction_ersing_etal" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_multilayer_convergence.jl"), + l2=[ + 0.0019005747730996517, + 0.0006416926161516576, + 0.0013237483460263738, + 0.0002514843620148158, + 0.00019551931612091513, + 0.00029396289378303574, + 0.0001391263772870812, + ], + linf=[ + 0.0059479250195759725, + 0.0028032845871422407, + 0.0039446424223523735, + 0.0009866378484075788, + 0.000930107485495435, + 0.0012826105967129742, + 0.00021874455861947695, + ], + surface_flux=(FluxHydrostaticReconstruction(FluxPlusDissipation(flux_ersing_etal, + DissipationLocalLaxFriedrichs()), + hydrostatic_reconstruction_ersing_etal), + FluxHydrostaticReconstruction(flux_nonconservative_ersing_etal, + hydrostatic_reconstruction_ersing_etal)), + tspan=(0.0, 0.25)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end + @trixi_testset "elixir_shallowwater_multilayer_well_balanced.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_multilayer_well_balanced.jl"), @@ -574,6 +611,63 @@ end # 2LSWE end end + @trixi_testset "elixir_shallowwater_multilayer_well_balanced_wet_dry.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_multilayer_well_balanced_wet_dry.jl"), + l2=[ + 0.025515633605684252, + 0.020321228314358498, + 2.078418885421747e-16, + 1.5034457079135132e-17, + 0.04746160729122716, + ], + linf=[ + 0.9999999999999989, + 1.0, + 1.8096275952411747e-15, + 4.245490280714591e-16, + 2.0, + ], + tspan=(0.0, 0.25)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end + + @trixi_testset "elixir_shallowwater_multilayer_well_balanced_wet_dry.jl with perturbation" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_multilayer_well_balanced_wet_dry.jl"), + l2=[ + 0.031114402256402142, + 0.02924759823866292, + 0.0021533928763101534, + 0.01412031518684538, + 0.04746160729122716, + ], + linf=[ + 1.2499999999999991, + 0.998672198687742, + 0.04297148454930247, + 0.24675612377725484, + 2.0, + ], + tspan=(0.0, 0.25), + perturbation=true) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end + @trixi_testset "elixir_shallowwater_multilayer_dam_break_ec.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_multilayer_dam_break_ec.jl"), @@ -637,6 +731,54 @@ end # 2LSWE @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 end end + + @trixi_testset "elixir_shallowwater_multilayer_beach.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_multilayer_beach.jl"), + l2=[ + 0.5837610281036327, + 3.7410178422784273, + 6.289710784508112e-8, + ], + linf=[ + 0.9917260088977286, + 7.0962810662023035, + 4.452604027704865e-7, + ], + tspan=(0.0, 0.25)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end + + @trixi_testset "elixir_shallowwater_multilayer_parabolic_bowl.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_multilayer_parabolic_bowl.jl"), + l2=[ + 0.002063065256405533, + 0.00048255459387538827, + 4.0617503413126664e-17, + ], + linf=[ + 0.0036250458403655466, + 0.0008665322585535845, + 1.6653345369377348e-16, + ], + tspan=(0.0, 0.25)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end end # MLSWE end # TreeMesh1D diff --git a/test/test_upstream.jl b/test/test_upstream.jl index 34b354b..f6a0ff6 100644 --- a/test/test_upstream.jl +++ b/test/test_upstream.jl @@ -42,7 +42,6 @@ isdir(outdir) && rm(outdir, recursive = true) @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 end end - # TODO: add upstream tests in 2D and positivity-preserving tests # Shallow water wet/dry 2D # TreeMesh2D From 5d6dafa06c8bd534e2cf372a8a70ffc72c869917 Mon Sep 17 00:00:00 2001 From: patrickersing Date: Tue, 9 Apr 2024 11:43:45 +0200 Subject: [PATCH 04/17] apply formatter --- .../elixir_shallowwater_multilayer_parabolic_bowl.jl | 3 ++- ...ixir_shallowwater_multilayer_well_balanced_wet_dry.jl | 9 ++++++--- src/equations/shallow_water_multilayer_1d.jl | 2 +- src/solvers/indicators_1d.jl | 2 +- 4 files changed, 10 insertions(+), 6 deletions(-) diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_multilayer_parabolic_bowl.jl b/examples/tree_1d_dgsem/elixir_shallowwater_multilayer_parabolic_bowl.jl index 8257775..ccc2dab 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_multilayer_parabolic_bowl.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_multilayer_parabolic_bowl.jl @@ -27,7 +27,8 @@ The particular setup below is taken from Section 6.2 of curvilinear meshes with wet/dry fronts accelerated by GPUs [DOI: 10.1016/j.jcp.2018.08.038](https://doi.org/10.1016/j.jcp.2018.08.038). """ -function initial_condition_parabolic_bowl(x, t, equations::ShallowWaterMultiLayerEquations1D) +function initial_condition_parabolic_bowl(x, t, + equations::ShallowWaterMultiLayerEquations1D) a = 1 h_0 = 0.1 sigma = 0.5 diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_multilayer_well_balanced_wet_dry.jl b/examples/tree_1d_dgsem/elixir_shallowwater_multilayer_well_balanced_wet_dry.jl index 1e271e0..1af8f29 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_multilayer_well_balanced_wet_dry.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_multilayer_well_balanced_wet_dry.jl @@ -6,7 +6,8 @@ using TrixiShallowWater ############################################################################### # Semidiscretization of the two-layer shallow water equations to test well-balancedness -equations = ShallowWaterMultiLayerEquations1D(gravity_constant = 1.0, H0 = 2.0, rhos = (1.0, 3.0)) +equations = ShallowWaterMultiLayerEquations1D(gravity_constant = 1.0, H0 = 2.0, + rhos = (1.0, 3.0)) """ initial_condition_twolayer_well_balanced_wet_dry(perturbation, equations:: ShallowWaterMultiLayerEquations1D) @@ -22,7 +23,8 @@ The initial condition is taken from the paper: A 1D shallow-flow model for two-layer flows based on FORCE scheme with wet–dry treatment\n [DOI:10.2166/hydro.2020.002](https://doi.org/10.2166/hydro.2020.002) """ -function initial_condition_twolayer_well_balanced_wet_dry(perturbation::Bool, equations::ShallowWaterMultiLayerEquations1D) +function initial_condition_twolayer_well_balanced_wet_dry(perturbation::Bool, + equations::ShallowWaterMultiLayerEquations1D) return function initial_condition_twolayer_well_balanced_wet_dry(x, t, equations) # Set interface height for the upper layer H_upper = 2.0 @@ -94,7 +96,8 @@ end # Default value for the perturbation perturbation = false -initial_condition = initial_condition_twolayer_well_balanced_wet_dry(perturbation, equations) +initial_condition = initial_condition_twolayer_well_balanced_wet_dry(perturbation, + equations) ############################################################################### # Get the DG approximation space diff --git a/src/equations/shallow_water_multilayer_1d.jl b/src/equations/shallow_water_multilayer_1d.jl index 07acbbc..99fd4b4 100644 --- a/src/equations/shallow_water_multilayer_1d.jl +++ b/src/equations/shallow_water_multilayer_1d.jl @@ -637,7 +637,7 @@ end equations::ShallowWaterMultiLayerEquations1D) h = waterheight(u, equations) - return 0.5 * equations.gravity * sum(h.^3) + return 0.5 * equations.gravity * sum(h .^ 3) end # Calculate the error for the "lake-at-rest" test case where H = ∑h + b should diff --git a/src/solvers/indicators_1d.jl b/src/solvers/indicators_1d.jl index be9c5cd..0bf149e 100644 --- a/src/solvers/indicators_1d.jl +++ b/src/solvers/indicators_1d.jl @@ -34,7 +34,7 @@ function (indicator_hg::IndicatorHennemannGassnerShallowWater)(u::AbstractArray{ the indicator sets the element-wise blending factor alpha[element] = 1 via the local variable `indicator_wet`. In turn, this ensures that a pure FV method is used in partially wet elements and guarantees the well-balanced property. - + Hard-coded cut-off value of `threshold_partially_wet = 1e-4` was determined through many numerical experiments. Overall idea is to increase robustness when computing the velocity on (nearly) dry elements which could be "dangerous" due to division of conservative variables, e.g., v = hv / h. From 12a794bd1346a55018fddbd381c580294ce71a9c Mon Sep 17 00:00:00 2001 From: patrickersing Date: Tue, 9 Apr 2024 13:16:18 +0200 Subject: [PATCH 05/17] update test values for new max_abs_speeds_naive --- test/test_tree_1d.jl | 48 ++++++++++++++++++++++---------------------- 1 file changed, 24 insertions(+), 24 deletions(-) diff --git a/test/test_tree_1d.jl b/test/test_tree_1d.jl index 6e6e982..81a4e30 100644 --- a/test/test_tree_1d.jl +++ b/test/test_tree_1d.jl @@ -512,21 +512,21 @@ end # 2LSWE @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_multilayer_convergence.jl"), l2=[ - 0.0019005751369889307, - 0.0006416934635148739, - 0.001323748067661235, - 0.0002514848406032196, - 0.00019551914601560465, - 0.00029396199135642574, + 0.0019005749358461095, + 0.0006416934559485677, + 0.0013237478952603332, + 0.00025148428672466024, + 0.00019551927320258115, + 0.0002939621441059124, 0.00013912637728693124, ], linf=[ - 0.005947927726240643, - 0.0028032845734302647, - 0.003944642659854947, - 0.0009866355929912807, - 0.0009301080749861135, - 0.0012826136652659414, + 0.0059479248159197695, + 0.0028032873366074518, + 0.003944640973376878, + 0.0009866380498015914, + 0.0009301097731234664, + 0.0012826130439229644, 0.0002187445586192549, ], surface_flux=(flux_lax_friedrichs, @@ -704,21 +704,21 @@ end # 2LSWE @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_multilayer_dam_break_es.jl"), l2=[ - 0.03715839782687889, - 0.03944992348741886, - 0.047363301719025856, - 0.28686684422935993, - 0.28693864344370795, - 0.21860509590943838, + 0.037157941241555185, + 0.03944955271506525, + 0.0473661084974011, + 0.28686995762998163, + 0.2869427613009073, + 0.21860712413448047, 0.013638618139749504, ], linf=[ - 0.18088706524197873, - 0.22218392524551045, - 0.4463223829540357, - 1.1462905426260082, - 1.14002097899466, - 0.9996443057617685, + 0.18060907942891546, + 0.2221666684315402, + 0.44639752069288585, + 1.146783259827299, + 1.140481663033156, + 0.9997998219394365, 0.4999999999999993, ], tspan=(0.0, 0.25)) From bf748825419ade2a3180d5302280e53a8ecc225a Mon Sep 17 00:00:00 2001 From: patrickersing Date: Tue, 9 Apr 2024 14:02:09 +0200 Subject: [PATCH 06/17] additional wet/dry test with multiple layers --- ...hallowwater_multilayer_dam_break_es_dry.jl | 124 ++++++++++++++++++ test/test_tree_1d.jl | 40 ++++++ 2 files changed, 164 insertions(+) create mode 100644 examples/tree_1d_dgsem/elixir_shallowwater_multilayer_dam_break_es_dry.jl diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_multilayer_dam_break_es_dry.jl b/examples/tree_1d_dgsem/elixir_shallowwater_multilayer_dam_break_es_dry.jl new file mode 100644 index 0000000..e474494 --- /dev/null +++ b/examples/tree_1d_dgsem/elixir_shallowwater_multilayer_dam_break_es_dry.jl @@ -0,0 +1,124 @@ + +using OrdinaryDiffEq +using Trixi +using TrixiShallowWater + +############################################################################### +# Semidiscretization of the multilayer shallow water equations for a dam break +# test over a dry domain with a discontinuous bottom topography function + +equations = ShallowWaterMultiLayerEquations1D(gravity_constant = 9.81, + rhos = (0.8, 0.85, 0.9, 0.95, 1.0)) + +# Initial condition of a dam break with a discontinuous water heights and bottom topography. +# To test the wet/dry functionality this test case considers a dam break over a dry domain with +# multiple layers. +# Works as intended for TreeMesh1D with `initial_refinement_level=5`. If the mesh +# refinement level is changed the initial condition below may need changed as well to +# ensure that the discontinuities lie on an element interface. +function initial_condition_dam_break(x, t, equations::ShallowWaterMultiLayerEquations1D) + # Set the discontinuity + if x[1] <= 10.0 + H = [3.5, 3.0, 2.5, 2.0, 1.5] + b = 0.0 + else + # Right side of the domain is dry + b = 0.5 + H = [b, b, b, b, b] + end + + # It is mandatory to shift the water level at dry areas to make sure the water height h + # stays positive. The system would not be stable for h set to a hard 0 due to division by h in + # the computation of velocity, e.g., (h v) / h. Therefore, a small dry state threshold + # with a default value of 5*eps() ≈ 1e-15 in double precision, is set in the constructor above + # for the ShallowWaterMultiLayerEquations1D and added to the initial condition if h = 0. + # This default value can be changed within the constructor call depending on the simulation setup. + for i in reverse(eachlayer(equations)) + if i == nlayers(equations) + H[i] = max(H[i], b + equations.threshold_limiter) + else + H[i] = max(H[i], H[i + 1] + equations.threshold_limiter) + end + end + + # Initialize zero velocity + v = zero(H) + + return prim2cons(SVector(H..., v..., b), equations) +end + +initial_condition = initial_condition_dam_break + +############################################################################### +# Get the DG approximation space + +volume_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal) +surface_flux = (FluxHydrostaticReconstruction(FluxPlusDissipation(flux_ersing_etal, + DissipationLocalLaxFriedrichs()), + hydrostatic_reconstruction_ersing_etal), + FluxHydrostaticReconstruction(flux_nonconservative_ersing_etal, + hydrostatic_reconstruction_ersing_etal)) + +basis = LobattoLegendreBasis(3) + +indicator_sc = IndicatorHennemannGassnerShallowWater(equations, basis, + alpha_max = 0.5, + alpha_min = 0.001, + alpha_smooth = true, + variable = waterheight_pressure) +volume_integral = VolumeIntegralShockCapturingHG(indicator_sc; + volume_flux_dg = volume_flux, + volume_flux_fv = surface_flux) + +solver = DGSEM(basis, surface_flux, volume_integral) + +############################################################################### +# Get the TreeMesh and setup a non-periodic mesh + +coordinates_min = 0.0 +coordinates_max = 20.0 +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 5, + n_cells_max = 10000, + periodicity = false) + +boundary_condition = boundary_condition_slip_wall + +# create the semidiscretization object +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_condition) + +############################################################################### +# ODE solvers + +tspan = (0.0, 1.0) +ode = semidiscretize(semi, tspan) + +############################################################################### +# Callbacks + +summary_callback = SummaryCallback() + +analysis_interval = 500 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + save_analysis = false) + +stepsize_callback = StepsizeCallback(cfl = 0.2) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 500, + save_initial_solution = true, + save_final_solution = true) + +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, + stepsize_callback, save_solution) + +stage_limiter! = PositivityPreservingLimiterShallowWater(variables = (Trixi.waterheight,)) + +############################################################################### +# run the simulation + +sol = solve(ode, SSPRK43(stage_limiter!), dt = 1.0, + save_everystep = false, callback = callbacks, adaptive = false); +summary_callback() # print the timer summary diff --git a/test/test_tree_1d.jl b/test/test_tree_1d.jl index 81a4e30..e0b81c6 100644 --- a/test/test_tree_1d.jl +++ b/test/test_tree_1d.jl @@ -732,6 +732,46 @@ end # 2LSWE end end + @trixi_testset "elixir_shallowwater_multilayer_dam_break_es_dry.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_multilayer_dam_break_es_dry.jl"), + l2=[ + 0.0539713484191238, + 0.05378401084239003, + 0.05343571027541134, + 0.05294334126567629, + 0.14409266517044308, + 0.2258690603183475, + 0.22395024283599815, + 0.22041584295217526, + 0.21549254848916716, + 0.6274449804372626, + 0.013638618139749523, + ], + linf=[ + 0.27418120077237185, + 0.2738885780879326, + 0.2733530772398547, + 0.2726148157344642, + 0.8243308057901404, + 0.6976533131655047, + 0.6927619903064818, + 0.6836813560099443, + 0.6708713913679563, + 3.0233275007119254, + 0.5, + ], + tspan=(0.0, 0.25)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end + @trixi_testset "elixir_shallowwater_multilayer_beach.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_multilayer_beach.jl"), From 063e0a9b29bb9e2eeba7f1779d191597f3088500 Mon Sep 17 00:00:00 2001 From: patrickersing Date: Tue, 9 Apr 2024 19:43:53 +0200 Subject: [PATCH 07/17] fix typo --- src/callbacks_stage/positivity_shallow_water.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/callbacks_stage/positivity_shallow_water.jl b/src/callbacks_stage/positivity_shallow_water.jl index 32ee722..7adff78 100644 --- a/src/callbacks_stage/positivity_shallow_water.jl +++ b/src/callbacks_stage/positivity_shallow_water.jl @@ -16,7 +16,7 @@ It is applied to all scalar `variables` in their given order using the defined `threshold_limiter` from the equations struct (e.g. in [`ShallowWaterEquationsWetDry1D`](@ref)) to determine the minimal acceptable values. The order of the `variables` is important and might have a strong influence -on the robustness. The limiter is availble for the [`ShallowWaterEquationsWetDry1D`](@ref), +on the robustness. The limiter is available for the [`ShallowWaterEquationsWetDry1D`](@ref), [`ShallowWaterEquationsWetDry2D`](@ref), and [`ShallowWaterMultiLayerEquations1D`](@ref). As opposed to the standard version of the [`Trixi.PositivityPreservingLimiterZhangShu`](@extref), From e64a07fd09258ed4690fedffbe6f5a52d341b3d0 Mon Sep 17 00:00:00 2001 From: patrickersing Date: Wed, 10 Apr 2024 08:59:19 +0200 Subject: [PATCH 08/17] merge branch 'main' into es_reconstruction_1d --- .../elixir_shallowwater_source_terms.jl | 2 +- .../elixir_shallowwater_well_balanced.jl | 2 +- .../tree_1d_dgsem/elixir_shallowwater_ec.jl | 2 +- .../elixir_shallowwater_well_balanced.jl | 2 +- ...ixir_shallowwater_well_balanced_wet_dry.jl | 2 +- .../tree_2d_dgsem/elixir_shallowwater_ec.jl | 2 +- .../elixir_shallowwater_well_balanced.jl | 2 +- .../elixir_shallowwater_well_balanced_wall.jl | 2 +- ...ixir_shallowwater_well_balanced_wet_dry.jl | 2 +- .../elixir_shallowwater_ec.jl | 2 +- .../elixir_shallowwater_ec_shockcapturing.jl | 2 +- .../elixir_shallowwater_well_balanced.jl | 2 +- test/test_structured_2d.jl | 12 +-- test/test_tree_1d.jl | 26 +++--- test/test_tree_2d.jl | 24 ++--- test/test_unstructured_2d.jl | 88 +++++++++---------- test/test_upstream.jl | 24 +++++ 17 files changed, 111 insertions(+), 87 deletions(-) diff --git a/examples/structured_2d_dgsem/elixir_shallowwater_source_terms.jl b/examples/structured_2d_dgsem/elixir_shallowwater_source_terms.jl index 20656cf..bf30def 100644 --- a/examples/structured_2d_dgsem/elixir_shallowwater_source_terms.jl +++ b/examples/structured_2d_dgsem/elixir_shallowwater_source_terms.jl @@ -43,7 +43,7 @@ save_solution = SaveSolutionCallback(interval = 100, save_final_solution = true, solution_variables = cons2prim) -stepsize_callback = StepsizeCallback(cfl = 2.0) +stepsize_callback = StepsizeCallback(cfl = 1.0) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, diff --git a/examples/structured_2d_dgsem/elixir_shallowwater_well_balanced.jl b/examples/structured_2d_dgsem/elixir_shallowwater_well_balanced.jl index 11c8d68..23f218b 100644 --- a/examples/structured_2d_dgsem/elixir_shallowwater_well_balanced.jl +++ b/examples/structured_2d_dgsem/elixir_shallowwater_well_balanced.jl @@ -110,7 +110,7 @@ save_solution = SaveSolutionCallback(interval = 1000, save_initial_solution = true, save_final_solution = true) -stepsize_callback = StepsizeCallback(cfl = 3.0) +stepsize_callback = StepsizeCallback(cfl = 1.0) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, stepsize_callback) diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_ec.jl b/examples/tree_1d_dgsem/elixir_shallowwater_ec.jl index 8046736..1281a15 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_ec.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_ec.jl @@ -79,7 +79,7 @@ save_solution = SaveSolutionCallback(interval = 100, save_initial_solution = true, save_final_solution = true) -stepsize_callback = StepsizeCallback(cfl = 3.0) +stepsize_callback = StepsizeCallback(cfl = 1.0) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, stepsize_callback) diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced.jl b/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced.jl index e78a867..ab46bc4 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced.jl @@ -74,7 +74,7 @@ save_solution = SaveSolutionCallback(interval = 1000, save_initial_solution = true, save_final_solution = true) -stepsize_callback = StepsizeCallback(cfl = 3.0) +stepsize_callback = StepsizeCallback(cfl = 1.0) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, stepsize_callback) 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 3f11b47..44b5931 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 @@ -101,7 +101,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.0) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, stepsize_callback) diff --git a/examples/tree_2d_dgsem/elixir_shallowwater_ec.jl b/examples/tree_2d_dgsem/elixir_shallowwater_ec.jl index 56f0e16..a24b7e3 100644 --- a/examples/tree_2d_dgsem/elixir_shallowwater_ec.jl +++ b/examples/tree_2d_dgsem/elixir_shallowwater_ec.jl @@ -106,7 +106,7 @@ save_solution = SaveSolutionCallback(dt = 0.2, save_initial_solution = true, save_final_solution = true) -stepsize_callback = StepsizeCallback(cfl = 3.0) +stepsize_callback = StepsizeCallback(cfl = 1.0) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, stepsize_callback) diff --git a/examples/tree_2d_dgsem/elixir_shallowwater_well_balanced.jl b/examples/tree_2d_dgsem/elixir_shallowwater_well_balanced.jl index 54b0aa2..3eaaf05 100644 --- a/examples/tree_2d_dgsem/elixir_shallowwater_well_balanced.jl +++ b/examples/tree_2d_dgsem/elixir_shallowwater_well_balanced.jl @@ -106,7 +106,7 @@ save_solution = SaveSolutionCallback(interval = 1000, save_initial_solution = true, save_final_solution = true) -stepsize_callback = StepsizeCallback(cfl = 3.0) +stepsize_callback = StepsizeCallback(cfl = 1.0) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, stepsize_callback) diff --git a/examples/tree_2d_dgsem/elixir_shallowwater_well_balanced_wall.jl b/examples/tree_2d_dgsem/elixir_shallowwater_well_balanced_wall.jl index d92008e..1ff4f3a 100644 --- a/examples/tree_2d_dgsem/elixir_shallowwater_well_balanced_wall.jl +++ b/examples/tree_2d_dgsem/elixir_shallowwater_well_balanced_wall.jl @@ -109,7 +109,7 @@ save_solution = SaveSolutionCallback(interval = 1000, save_initial_solution = true, save_final_solution = true) -stepsize_callback = StepsizeCallback(cfl = 3.0) +stepsize_callback = StepsizeCallback(cfl = 1.0) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, stepsize_callback) diff --git a/examples/tree_2d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl b/examples/tree_2d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl index 8135e28..d451f63 100644 --- a/examples/tree_2d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl +++ b/examples/tree_2d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl @@ -134,7 +134,7 @@ save_solution = SaveSolutionCallback(interval = 1000, save_initial_solution = true, save_final_solution = true) -stepsize_callback = StepsizeCallback(cfl = 2.0) +stepsize_callback = StepsizeCallback(cfl = 1.0) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, stepsize_callback) diff --git a/examples/unstructured_2d_dgsem/elixir_shallowwater_ec.jl b/examples/unstructured_2d_dgsem/elixir_shallowwater_ec.jl index 1415335..0eb9d33 100644 --- a/examples/unstructured_2d_dgsem/elixir_shallowwater_ec.jl +++ b/examples/unstructured_2d_dgsem/elixir_shallowwater_ec.jl @@ -110,7 +110,7 @@ save_solution = SaveSolutionCallback(interval = 1000, save_initial_solution = true, save_final_solution = true) -stepsize_callback = StepsizeCallback(cfl = 3.0) +stepsize_callback = StepsizeCallback(cfl = 1.0) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, stepsize_callback) diff --git a/examples/unstructured_2d_dgsem/elixir_shallowwater_ec_shockcapturing.jl b/examples/unstructured_2d_dgsem/elixir_shallowwater_ec_shockcapturing.jl index 2b550c6..78d6cc7 100644 --- a/examples/unstructured_2d_dgsem/elixir_shallowwater_ec_shockcapturing.jl +++ b/examples/unstructured_2d_dgsem/elixir_shallowwater_ec_shockcapturing.jl @@ -115,7 +115,7 @@ analysis_callback = AnalysisCallback(semi, interval = analysis_interval) alive_callback = AliveCallback(analysis_interval = analysis_interval) -stepsize_callback = StepsizeCallback(cfl = 3.0) +stepsize_callback = StepsizeCallback(cfl = 1.0) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, stepsize_callback) diff --git a/examples/unstructured_2d_dgsem/elixir_shallowwater_well_balanced.jl b/examples/unstructured_2d_dgsem/elixir_shallowwater_well_balanced.jl index e0b76c6..18a497c 100644 --- a/examples/unstructured_2d_dgsem/elixir_shallowwater_well_balanced.jl +++ b/examples/unstructured_2d_dgsem/elixir_shallowwater_well_balanced.jl @@ -110,7 +110,7 @@ save_solution = SaveSolutionCallback(interval = 1000, save_initial_solution = true, save_final_solution = true) -stepsize_callback = StepsizeCallback(cfl = 3.0) +stepsize_callback = StepsizeCallback(cfl = 1.0) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, stepsize_callback) diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl index 506d903..fac1fe6 100644 --- a/test/test_structured_2d.jl +++ b/test/test_structured_2d.jl @@ -20,15 +20,15 @@ isdir(outdir) && rm(outdir, recursive = true) @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"), l2=[ - 0.0017285599436729316, - 0.025584610912606776, - 0.028373834961180594, + 0.0017286908591070864, + 0.025585037307655684, + 0.028374244567802766, 6.274146767730866e-5, ], linf=[ - 0.012972309788264802, - 0.108283714215621, - 0.15831585777928936, + 0.012973752001194772, + 0.10829375385832263, + 0.15832858475438094, 0.00018196759554722775, ], tspan=(0.0, 0.05)) diff --git a/test/test_tree_1d.jl b/test/test_tree_1d.jl index e0b81c6..bfe0608 100644 --- a/test/test_tree_1d.jl +++ b/test/test_tree_1d.jl @@ -19,13 +19,13 @@ isdir(outdir) && rm(outdir, recursive = true) @trixi_testset "elixir_shallowwater_ec.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_ec.jl"), l2=[ - 0.244729018751225, - 0.8583565222389505, + 0.24476140682560343, + 0.8587309324660326, 0.07330427577586297, ], linf=[ - 2.1635021283528504, - 3.8717508164234453, + 2.1636963952308372, + 3.8737770522883115, 1.7711213427919539, ], tspan=(0.0, 0.25)) @@ -42,13 +42,13 @@ isdir(outdir) && rm(outdir, recursive = true) @trixi_testset "elixir_shallowwater_ec.jl with initial_condition_weak_blast_wave" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_ec.jl"), l2=[ - 0.39464782107209717, - 2.03880864210846, + 0.39472828074570576, + 2.0390687947320076, 4.1623084150546725e-10, ], linf=[ - 0.778905801278281, - 3.2409883402608273, + 0.7793741954662221, + 3.2411927977882096, 7.419800190922032e-10, ], initial_condition=initial_condition_weak_blast_wave, @@ -134,14 +134,14 @@ isdir(outdir) && rm(outdir, recursive = true) @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_well_balanced_wet_dry.jl"), l2=[ - 0.00965787167169024, - 5.345454081916856e-14, - 0.03857583749209928, + 0.009657871671690306, + 6.23878532543316e-14, + 0.03857583749209947, ], linf=[ 0.4999999999998892, - 2.2447689894899726e-13, - 1.9999999999999714, + 4.4857500383821557e-13, + 1.9999999999999811, ], tspan=(0.0, 0.25), # Soften the tolerance as test results vary between different CPUs diff --git a/test/test_tree_2d.jl b/test/test_tree_2d.jl index 214e84f..eec03b4 100644 --- a/test/test_tree_2d.jl +++ b/test/test_tree_2d.jl @@ -19,15 +19,15 @@ isdir(outdir) && rm(outdir, recursive = true) @trixi_testset "elixir_shallowwater_ec.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_ec.jl"), l2=[ - 0.991181203601035, - 0.734130029040644, - 0.7447696147162621, + 0.9911802019934329, + 0.7340106828033273, + 0.7446338002084801, 0.5875351036989047, ], linf=[ - 2.0117744577945413, - 2.9962317608172127, - 2.6554999727293653, + 2.0120253138457564, + 2.991158989293406, + 2.6557412817714035, 3.0, ], tspan=(0.0, 0.25)) @@ -344,15 +344,15 @@ isdir(outdir) && rm(outdir, recursive = true) @trixi_testset "elixir_shallowwater_wall.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_wall.jl"), l2=[ - 0.13517233723296504, - 0.20010876311162215, - 0.20010876311162223, + 0.1351723240085936, + 0.20010881416550014, + 0.2001088141654999, 2.719538414346464e-7, ], linf=[ - 0.5303607982988336, - 0.5080989745682338, - 0.5080989745682352, + 0.5303608302490757, + 0.5080987791967457, + 0.5080987791967506, 1.1301675764130437e-6, ], tspan=(0.0, 0.25)) diff --git a/test/test_unstructured_2d.jl b/test/test_unstructured_2d.jl index cd8c5e7..5fcfa32 100644 --- a/test/test_unstructured_2d.jl +++ b/test/test_unstructured_2d.jl @@ -19,15 +19,15 @@ isdir(outdir) && rm(outdir, recursive = true) @trixi_testset "elixir_shallowwater_ec.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_ec.jl"), l2=[ - 0.6106939484178353, - 0.48586236867426724, - 0.48234490854514356, - 0.29467422718511727, + 0.6107326269462766, + 0.48666631722018877, + 0.48309775159067053, + 0.29467422718511704, ], linf=[ - 2.775979948281604, - 3.1721242154451548, - 3.5713448319601393, + 2.776782342826098, + 3.2158378644333707, + 3.652920889487258, 2.052861364219655, ], tspan=(0.0, 0.25)) @@ -133,16 +133,16 @@ isdir(outdir) && rm(outdir, recursive = true) @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"), l2=[ - 0.0011197623982310795, - 0.04456344888447023, - 0.014317376629669337, - 5.089218476758975e-6, + 0.001118134082248467, + 0.044560486817464634, + 0.01430926600634214, + 5.089218476759981e-6, ], linf=[ - 0.007835284004819698, - 0.3486891284278597, - 0.11242778979399048, - 2.6407324614119432e-5, + 0.007798727223654822, + 0.34782952734839157, + 0.11161614702628064, + 2.6407324614341476e-5, ], tspan=(0.0, 0.025)) # Ensure that we do not have excessive memory allocations @@ -159,16 +159,16 @@ isdir(outdir) && rm(outdir, recursive = true) @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"), l2=[ - 0.0011197139793938152, - 0.015430259691310781, - 0.017081031802719724, - 5.089218476758271e-6, + 0.001119719074048042, + 0.015430106961377458, + 0.017081080392855975, + 5.0892184767572974e-6, ], linf=[ - 0.014300809338967824, - 0.12783372461225184, - 0.17625472321992852, - 2.6407324614341476e-5, + 0.014300786288935274, + 0.12782194018738569, + 0.1762517408268396, + 2.640732461456352e-5, ], surface_flux=(FluxHydrostaticReconstruction(FluxHLL(min_max_speed_naive), hydrostatic_reconstruction_audusse_etal), @@ -188,15 +188,15 @@ isdir(outdir) && rm(outdir, recursive = true) @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"), l2=[ - 0.0011196687776346434, - 0.044562672453443995, - 0.014306265289763618, + 0.001118046975499805, + 0.04455969246244461, + 0.014298120235633432, 5.089218476759981e-6, ], linf=[ - 0.007825021762002393, - 0.348550815397918, - 0.1115517935018282, + 0.007776521213640031, + 0.34768318303226353, + 0.11075311228066198, 2.6407324614341476e-5, ], surface_flux=(flux_wintermeyer_etal, @@ -218,16 +218,16 @@ isdir(outdir) && rm(outdir, recursive = true) @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"), l2=[ - 0.0011197139793938727, - 0.015430259691311309, - 0.017081031802719554, - 5.089218476759981e-6, + 0.0011197190740480426, + 0.015430106961377272, + 0.017081080392856225, + 5.0892184767572974e-6, ], linf=[ - 0.014300809338967824, - 0.12783372461224918, - 0.17625472321993918, - 2.6407324614341476e-5, + 0.014300786288934386, + 0.1278219401874079, + 0.1762517408268396, + 2.640732461456352e-5, ], surface_flux=(FluxHLL(min_max_speed_naive), flux_nonconservative_fjordholm_etal), @@ -303,15 +303,15 @@ isdir(outdir) && rm(outdir, recursive = true) @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_ec_shockcapturing.jl"), l2=[ - 0.6124656312639043, - 0.504371951785709, - 0.49180896200746366, - 0.29467422718511727, + 0.612551520607341, + 0.5039173660221961, + 0.49136517934903523, + 0.29467422718511704, ], linf=[ - 2.7639232436274392, - 3.3985508653311767, - 3.3330308209196224, + 2.7636771472622197, + 3.236168963021072, + 3.3363936775653826, 2.052861364219655, ], tspan=(0.0, 0.25)) diff --git a/test/test_upstream.jl b/test/test_upstream.jl index f6a0ff6..1a5a713 100644 --- a/test/test_upstream.jl +++ b/test/test_upstream.jl @@ -18,6 +18,30 @@ isdir(outdir) && rm(outdir, recursive = true) # Run tests for TreeMesh @testset "TreeMesh" begin # Shallow water wet/dry 1D + @trixi_testset "elixir_shallowwater_ec.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "tree_1d_dgsem", + "elixir_shallowwater_ec.jl"), + l2=[ + 0.24476140682560343, + 0.8587309324660326, + 0.07330427577586297, + ], + linf=[ + 2.1636963952308372, + 3.8737770522883115, + 1.7711213427919539, + ], + tspan=(0.0, 0.25)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end + @trixi_testset "TreeMesh1D: elixir_shallowwater_well_balanced_nonperiodic.jl with wall boundary" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "tree_1d_dgsem", "elixir_shallowwater_well_balanced_nonperiodic.jl"), From b9302cd28ae88a0a597f1fad0e28ddad100db9d4 Mon Sep 17 00:00:00 2001 From: patrickersing Date: Wed, 10 Apr 2024 15:48:21 +0200 Subject: [PATCH 09/17] add reference for desingularization --- src/callbacks_stage/positivity_shallow_water.jl | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/callbacks_stage/positivity_shallow_water.jl b/src/callbacks_stage/positivity_shallow_water.jl index 7adff78..253e629 100644 --- a/src/callbacks_stage/positivity_shallow_water.jl +++ b/src/callbacks_stage/positivity_shallow_water.jl @@ -35,7 +35,10 @@ there could still be dry nodes afterwards due to the logic of the limiter. For the [`ShallowWaterMultiLayerEquations1D`](@ref) the implementation differs. In this case the positivity limiter is applied layerwise and only the waterheight `h` is limited within each layer. Furthermore, a velocity desingularization is applied after the limiting to avoid numerical problems -near dry states. +near dry states. Details about the desingularization strategy can be found in Section 2.2 of the paper +- A. Kurganov, G. Petrova (2007) + A second-order well-balanced positivity preserving central-upwind scheme for the Saint-Venant system + [doi: 10.4310/CMS.2007.v5.n1.a6](https://dx.doi.org/10.4310/CMS.2007.v5.n1.a6) This fully-discrete positivity-preserving limiter is based on the work of - Zhang, Shu (2011) From 3d6d89a9381cd544f4ee470550c38ee47df42bf1 Mon Sep 17 00:00:00 2001 From: patrickersing Date: Thu, 11 Apr 2024 14:17:50 +0200 Subject: [PATCH 10/17] address changes from code review --- .../elixir_shallowwater_multilayer_beach.jl | 1 + ..._shallowwater_multilayer_parabolic_bowl.jl | 1 + ...wwater_multilayer_well_balanced_wet_dry.jl | 34 ++----------------- .../positivity_shallow_water.jl | 7 ++-- .../positivity_shallow_water_dg1d.jl | 6 ++-- src/equations/shallow_water_multilayer_1d.jl | 28 +++++++-------- 6 files changed, 26 insertions(+), 51 deletions(-) diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_multilayer_beach.jl b/examples/tree_1d_dgsem/elixir_shallowwater_multilayer_beach.jl index 91767fa..ec3234f 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_multilayer_beach.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_multilayer_beach.jl @@ -6,6 +6,7 @@ using TrixiShallowWater ############################################################################### # Semidiscretization of the shallow water equations +# By passing only a single value for rhos, the system recovers the standard shallow water equations equations = ShallowWaterMultiLayerEquations1D(gravity_constant = 9.812, rhos = 1.0) """ diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_multilayer_parabolic_bowl.jl b/examples/tree_1d_dgsem/elixir_shallowwater_multilayer_parabolic_bowl.jl index ccc2dab..d48557e 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_multilayer_parabolic_bowl.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_multilayer_parabolic_bowl.jl @@ -6,6 +6,7 @@ using TrixiShallowWater ############################################################################### # Semidiscretization of the shallow water equations +# By passing only a single value for rhos, the system recovers the standard shallow water equations equations = ShallowWaterMultiLayerEquations1D(gravity_constant = 9.81, rhos = 1.0) """ diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_multilayer_well_balanced_wet_dry.jl b/examples/tree_1d_dgsem/elixir_shallowwater_multilayer_well_balanced_wet_dry.jl index 1af8f29..b82ddba 100644 --- a/examples/tree_1d_dgsem/elixir_shallowwater_multilayer_well_balanced_wet_dry.jl +++ b/examples/tree_1d_dgsem/elixir_shallowwater_multilayer_well_balanced_wet_dry.jl @@ -75,8 +75,8 @@ function initial_condition_twolayer_well_balanced_wet_dry(perturbation::Bool, end # Set zero initial velocity - v1_upper = 0.0 - v1_lower = 0.0 + v1_upper = zero(H_upper) + v1_lower = zero(H_upper) #= It is mandatory to shift the water level at dry areas to make sure the water height h @@ -140,34 +140,6 @@ semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) tspan = (0.0, 10.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 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])) - elseif i == nnodes(semi.solver) - x_node = SVector(prevfloat(x_node[1])) - end - u_node = initial_condition(x_node, first(tspan), equations) - Trixi.set_node_vars!(u, u_node, equations, semi.solver, i, element) - end -end ############################################################################################# summary_callback = SummaryCallback() @@ -193,7 +165,7 @@ callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, sav ############################################################################### # run the simulation -# use a Runge-Kutta method with automatic (error based) time step size control +# use a Runge-Kutta method with CFL-based time step sol = solve(ode, SSPRK43(stage_limiter!); ode_default_options()..., callback = callbacks, adaptive = false, dt = 1.0); summary_callback() # print the timer summary diff --git a/src/callbacks_stage/positivity_shallow_water.jl b/src/callbacks_stage/positivity_shallow_water.jl index 253e629..8ed1630 100644 --- a/src/callbacks_stage/positivity_shallow_water.jl +++ b/src/callbacks_stage/positivity_shallow_water.jl @@ -8,9 +8,6 @@ """ PositivityPreservingLimiterShallowWater(; variables) -!!! warning "Experimental code" - This is an experimental feature and may change in future releases. - The limiter is specifically designed for the shallow water equations. It is applied to all scalar `variables` in their given order using the defined `threshold_limiter` from the equations struct @@ -50,6 +47,10 @@ The specific implementation for the [`ShallowWaterMultiLayerEquations1D](@ref) i Positivity-preserving well-balanced discontinuous Galerkin methods for the shallow water equations on unstructured triangular meshes [doi: 10.1007/s10915-012-9644-4](https://doi.org/10.1007/s10915-012-9644-4) + + +!!! warning "Experimental code" + This is an experimental feature and may change in future releases. """ struct PositivityPreservingLimiterShallowWater{N, Variables <: NTuple{N, Any}} variables::Variables diff --git a/src/callbacks_stage/positivity_shallow_water_dg1d.jl b/src/callbacks_stage/positivity_shallow_water_dg1d.jl index 9f6bbe7..14b77e4 100644 --- a/src/callbacks_stage/positivity_shallow_water_dg1d.jl +++ b/src/callbacks_stage/positivity_shallow_water_dg1d.jl @@ -87,12 +87,12 @@ function limiter_shallow_water!(u, threshold::Real, variable, return nothing end -# !!! warning "Experimental code" -# This is an experimental feature and may change in future releases. - # Note that for the `ShallowWaterMultiLayerEquations1D` only the waterheight `h` is limited in # each layer. Furthermore, a velocity desingularization is applied after the limiting to avoid # numerical problems near dry states. +# +# !!! warning "Experimental code" +# This is an experimental feature and may change in future releases. function limiter_shallow_water!(u, threshold::Real, variable, mesh::Trixi.AbstractMesh{1}, equations::ShallowWaterMultiLayerEquations1D, diff --git a/src/equations/shallow_water_multilayer_1d.jl b/src/equations/shallow_water_multilayer_1d.jl index a7cc908..6c6dabb 100644 --- a/src/equations/shallow_water_multilayer_1d.jl +++ b/src/equations/shallow_water_multilayer_1d.jl @@ -135,13 +135,13 @@ function Trixi.varnames(::typeof(cons2cons), return (waterheight..., momentum..., "b") end -# We use the interface heights, H = ∑h + b as primitive variables for easier visualization and setting initial +# We use the total layer heights, H = ∑h + b as primitive variables for easier visualization and setting initial # conditions function Trixi.varnames(::typeof(cons2prim), equations::ShallowWaterMultiLayerEquations1D) - interface_height = ntuple(n -> "H" * string(n), Val(nlayers(equations))) + total_layer_height = ntuple(n -> "H" * string(n), Val(nlayers(equations))) velocity = ntuple(n -> "v" * string(n) * "_1", Val(nlayers(equations))) - return (interface_height..., velocity..., "b") + return (total_layer_height..., velocity..., "b") end # Set initial conditions at physical location `x` for time `t` @@ -379,9 +379,6 @@ end """ hydrostatic_reconstruction_ersing_etal(u_ll, u_rr, equations::ShallowWaterMultiLayerEquations1D) -!!! warning "Experimental code" - This is an experimental feature and may change in future releases. - A particular type of hydrostatic reconstruction of the water height and bottom topography to guarantee well-balancedness in the presence of wet/dry transitions and entropy stability for the [`ShallowWaterMultiLayerEquations1D`](@ref). @@ -390,6 +387,9 @@ surface numerical flux at the interface. The key idea is a piecewise linear reco bottom topography and water height interfaces using subcells, where the bottom topography is allowed to be discontinuous. Use in combination with the generic numerical flux routine [`Trixi.FluxHydrostaticReconstruction`](@extref). + +!!! warning "Experimental code" + This is an experimental feature and may change in future releases. """ @inline function hydrostatic_reconstruction_ersing_etal(u_ll, u_rr, equations::ShallowWaterMultiLayerEquations1D) @@ -415,7 +415,7 @@ Use in combination with the generic numerical flux routine [`Trixi.FluxHydrostat end end - # Calculate interface heights + # Calculate total layer heights H_ll = waterheight(cons2prim(u_ll, equations), equations) H_rr = waterheight(cons2prim(u_rr, equations), equations) @@ -423,7 +423,7 @@ Use in combination with the generic numerical flux routine [`Trixi.FluxHydrostat b_ll_star = min(H_ll[1], max(b_rr, b_ll)) b_rr_star = min(H_rr[1], max(b_rr, b_ll)) - # Calculate reconstructed interface heights + # Calculate reconstructed total layer heights H_ll_star = max.(H_ll, b_ll_star) H_rr_star = max.(H_rr, b_rr_star) @@ -501,7 +501,7 @@ end h = waterheight(u, equations) b = u[end] - # Initialize interface height + # Initialize total layer height H = MVector{nlayers(equations), real(equations)}(undef) for i in reverse(eachlayer(equations)) if i == nlayers(equations) @@ -517,9 +517,9 @@ end # Convert primitive to conservative variables @inline function Trixi.prim2cons(prim, equations::ShallowWaterMultiLayerEquations1D) - # To extract the interface height and velocity we reuse the waterheight and momentum functions + # To extract the total layer height and velocity we reuse the waterheight and momentum functions # from the conservative variables. - H = waterheight(prim, equations) # For primitive variables this extracts the interface height + H = waterheight(prim, equations) # For primitive variables this extracts the total layer height v = momentum(prim, equations) # For primitive variables this extracts the velocity b = prim[end] @@ -556,9 +556,9 @@ end # Calculate entropy variables in each layer for i in eachlayer(equations) - # Compute w1[i] = ρ[i] * g * (b + ∑h[k] + ∑σ[k] * h[k]), where σ[k] = ρ[k] / ρ[i] denotes - # the density ratio of different layers - w1 = equations.rhos[i] * (g * b - 0.5 * v[i]^2) + # Compute w1[i] = ρ[i] * g * (b + ∑h[k] + ∑σ[k] * h[k]) - 0.5 * ρ[i] * v[i]^2, + # where σ[k] = ρ[k] / ρ[i] denotes the density ratio of different layers + w1 = equations.rhos[i] * (g * b) - 0.5 * equations.rhos[i] * v[i]^2 for j in eachlayer(equations) if j < i w1 += equations.rhos[i] * g * From 2857753d5373484df8c3da743e1a28e5144b09a5 Mon Sep 17 00:00:00 2001 From: patrickersing Date: Mon, 22 Apr 2024 15:23:06 +0200 Subject: [PATCH 11/17] soften tolerances to fix macOS testing --- test/test_structured_2d.jl | 3 ++- test/test_tree_1d.jl | 8 +++++--- 2 files changed, 7 insertions(+), 4 deletions(-) diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl index fac1fe6..eff0178 100644 --- a/test/test_structured_2d.jl +++ b/test/test_structured_2d.jl @@ -83,7 +83,8 @@ isdir(outdir) && rm(outdir, recursive = true) 4.4849667259339357e-14, 1.9999999999999993, ], - tspan=(0.0, 0.25)) + tspan=(0.0, 0.25), + atol=1e-12) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let diff --git a/test/test_tree_1d.jl b/test/test_tree_1d.jl index bfe0608..8323a71 100644 --- a/test/test_tree_1d.jl +++ b/test/test_tree_1d.jl @@ -145,7 +145,7 @@ isdir(outdir) && rm(outdir, recursive = true) ], tspan=(0.0, 0.25), # Soften the tolerance as test results vary between different CPUs - atol=1000 * eps()) + atol=1e-12) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let @@ -785,7 +785,8 @@ end # 2LSWE 7.0962810662023035, 4.452604027704865e-7, ], - tspan=(0.0, 0.25)) + tspan=(0.0, 0.25), + atol=1e-5) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let @@ -809,7 +810,8 @@ end # 2LSWE 0.0008665322585535845, 1.6653345369377348e-16, ], - tspan=(0.0, 0.25)) + tspan=(0.0, 0.25), + atol=1e-9) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let From 59b8ebd1f7d33abb6997f27848cce07cc07a1c86 Mon Sep 17 00:00:00 2001 From: patrickersing Date: Mon, 22 Apr 2024 16:14:37 +0200 Subject: [PATCH 12/17] rework thresholds --- src/equations/equations.jl | 8 +++++++ src/equations/shallow_water_multilayer_1d.jl | 25 +++++++++++++++----- src/equations/shallow_water_wet_dry_1d.jl | 18 +++++++++++--- src/equations/shallow_water_wet_dry_2d.jl | 18 +++++++++++--- src/solvers/indicators_1d.jl | 2 +- src/solvers/indicators_2d.jl | 2 +- 6 files changed, 59 insertions(+), 14 deletions(-) diff --git a/src/equations/equations.jl b/src/equations/equations.jl index a36f993..b92a6fc 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -41,4 +41,12 @@ Retrieve the number of layers from an equation instance of the `AbstractShallowW } NLAYERS end + +# Provide default thresholds dependend on number format (Currently default thresholds are only provided +# for Float64) +default_threshold_partially_wet(::Type{Float64}) = 1e-4 +default_threshold_partially_wet(catchall) = throw(ArgumentError("threshold_partially_wet must be provided for non-Float64 types")) + +default_threshold_desingularization(::Type{Float64}) = 1e-10 +default_threshold_desingularization(catchall) = throw(ArgumentError("threshold_desingularization must be provided for non-Float64 types")) end # @muladd diff --git a/src/equations/shallow_water_multilayer_1d.jl b/src/equations/shallow_water_multilayer_1d.jl index 6c6dabb..5fe9a90 100644 --- a/src/equations/shallow_water_multilayer_1d.jl +++ b/src/equations/shallow_water_multilayer_1d.jl @@ -30,11 +30,15 @@ nonconservative term, which has some benefits for the design of well-balanced sc The additional quantity ``H_0`` is also available to store a reference value for the total water height that is useful to set initial conditions or test the "lake-at-rest" well-balancedness. -Also, there are two thresholds which prevent numerical problems as well as instabilities. Both of -them do not have to be passed, as default values are defined within the struct. The limiters are +Also, there are two thresholds which prevent numerical problems as well as instabilities. The limiters are used in [`PositivityPreservingLimiterShallowWater`](@ref) on the water height. `threshold_limiter` acts as a (small) shift on the initial condition and cutoff before the next time step, whereas -`threshold_desingularization` is used in the velocity desingularization. +`threshold_desingularization` is used in the velocity desingularization. A third +`threshold_partially_wet` is applied on the water height to define "partially wet" elements in +[`IndicatorHennemannGassnerShallowWater`](@ref), that are then calculated with a pure FV method to +ensure well-balancedness. For `Float64` no threshold needs to be passed, as default values are +defined within the struct. For other number formats `threshold_desingularization` and `threshold_partially_wet` +must be provided. The bottom topography function ``b(x)`` is set inside the initial condition routine for a particular problem setup. @@ -62,12 +66,14 @@ struct ShallowWaterMultiLayerEquations1D{NVARS, NLAYERS, RealT <: Real} <: H0::RealT # constant "lake-at-rest" total water height threshold_limiter::RealT # threshold for the positivity-limiter threshold_desingularization::RealT # threshold for velocity desingularization + threshold_partially_wet::RealT # threshold to define partially wet elements rhos::SVector{NLAYERS, RealT} # Vector of layer densities function ShallowWaterMultiLayerEquations1D{NVARS, NLAYERS, RealT}(gravity::RealT, H0::RealT, threshold_limiter::RealT, threshold_desingularization::RealT, + threshold_partially_wet::RealT, rhos::SVector{NLAYERS, RealT}) where { NVARS, @@ -80,7 +86,8 @@ struct ShallowWaterMultiLayerEquations1D{NVARS, NLAYERS, RealT <: Real} <: throw(ArgumentError("densities must be in increasing order (rhos[1] < rhos[2] < ... < rhos[NLAYERS])")) min(rhos...) > 0 || throw(ArgumentError("densities must be positive")) - new(gravity, H0, threshold_limiter, threshold_desingularization, rhos) + new(gravity, H0, threshold_limiter, threshold_desingularization, + threshold_partially_wet, rhos) end end @@ -91,7 +98,9 @@ end function ShallowWaterMultiLayerEquations1D(; gravity_constant, H0 = zero(gravity_constant), threshold_limiter = nothing, - threshold_desingularization = nothing, rhos) + threshold_desingularization = nothing, + threshold_partially_wet = nothing, + rhos) # Promote all variables to a common type _rhos = promote(rhos...) @@ -103,7 +112,10 @@ function ShallowWaterMultiLayerEquations1D(; gravity_constant, threshold_limiter = 5 * eps(RealT) end if threshold_desingularization === nothing - threshold_desingularization = 1e-10 + threshold_desingularization = default_threshold_desingularization(RealT) + end + if threshold_partially_wet === nothing + threshold_partially_wet = default_threshold_partially_wet(RealT) end # Extract number of layers and variables @@ -114,6 +126,7 @@ function ShallowWaterMultiLayerEquations1D(; gravity_constant, H0, threshold_limiter, threshold_desingularization, + threshold_partially_wet, __rhos) end diff --git a/src/equations/shallow_water_wet_dry_1d.jl b/src/equations/shallow_water_wet_dry_1d.jl index c395786..eacd5ac 100644 --- a/src/equations/shallow_water_wet_dry_1d.jl +++ b/src/equations/shallow_water_wet_dry_1d.jl @@ -28,7 +28,11 @@ Also, there are two thresholds which prevent numerical problems as well as insta have to be passed, as default values are defined within the struct. The first one, `threshold_limiter`, is used in [`PositivityPreservingLimiterShallowWater`](@ref) on the water height, as a (small) shift on the initial condition and cutoff before the next time step. The second one, `threshold_wet`, is applied on the water height to -define when the flow is "wet" before calculating the numerical flux. +define when the flow is "wet" before calculating the numerical flux. A third +`threshold_partially_wet` is applied on the water height to define "partially wet" elements in +[`IndicatorHennemannGassnerShallowWater`](@ref), that are then calculated with a pure FV method to +ensure well-balancedness. For `Float64` no threshold needs to be passed, as default values are +defined within the struct. For other number formats `threshold_partially_wet` must be provided. The bottom topography function ``b(x)`` is set inside the initial condition routine for a particular problem setup. To test the conservative form of the SWE one can set the bottom topography @@ -62,6 +66,10 @@ struct ShallowWaterEquationsWetDry1D{RealT <: Real} <: # before calculating the numerical flux. # Default is 5*eps() which in double precision is ≈1e-15. threshold_wet::RealT + # `threshold_partially_wet` used in `IndicatorHennemannGassnerShallowWater` on the water height + # to define "partially wet" elements. Those elements are calculated with a pure FV method to + # ensure well-balancedness. Default in double precision is 1e-4. + threshold_partially_wet::RealT # Standard shallow water equations for dispatch on Trixi.jl functions basic_swe::ShallowWaterEquations1D{RealT} end @@ -73,7 +81,8 @@ end # Strict default values for thresholds that performed well in many numerical experiments function ShallowWaterEquationsWetDry1D(; gravity_constant, H0 = zero(gravity_constant), threshold_limiter = nothing, - threshold_wet = nothing) + threshold_wet = nothing, + threshold_partially_wet = nothing) T = promote_type(typeof(gravity_constant), typeof(H0)) if threshold_limiter === nothing threshold_limiter = 500 * eps(T) @@ -81,13 +90,16 @@ function ShallowWaterEquationsWetDry1D(; gravity_constant, H0 = zero(gravity_con if threshold_wet === nothing threshold_wet = 5 * eps(T) end + if threshold_partially_wet === nothing + threshold_partially_wet = default_threshold_partially_wet(T) + end # Construct the standard SWE for dispatch. Even though the `basic_swe` already store the # gravity constant and the total water height, we store an extra copy in # `ShallowWaterEquationsWetDry1D` for convenience. basic_swe = ShallowWaterEquations1D(gravity_constant = gravity_constant, H0 = H0) ShallowWaterEquationsWetDry1D(gravity_constant, H0, threshold_limiter, - threshold_wet, basic_swe) + threshold_wet, threshold_partially_wet, basic_swe) end Trixi.have_nonconservative_terms(::ShallowWaterEquationsWetDry1D) = True() diff --git a/src/equations/shallow_water_wet_dry_2d.jl b/src/equations/shallow_water_wet_dry_2d.jl index 32c7e81..1db2a5c 100644 --- a/src/equations/shallow_water_wet_dry_2d.jl +++ b/src/equations/shallow_water_wet_dry_2d.jl @@ -31,7 +31,11 @@ Also, there are two thresholds which prevent numerical problems as well as insta have to be passed, as default values are defined within the struct. The first one, `threshold_limiter`, is used in [`PositivityPreservingLimiterShallowWater`](@ref) on the water height, as a (small) shift on the initial condition and cutoff before the next time step. The second one, `threshold_wet`, is applied on the water height to -define when the flow is "wet" before calculating the numerical flux. +define when the flow is "wet" before calculating the numerical flux. A third +`threshold_partially_wet` is applied on the water height to define "partially wet" elements in +[`IndicatorHennemannGassnerShallowWater`](@ref), that are then calculated with a pure FV method to +ensure well-balancedness. For `Float64` no threshold needs to be passed, as default values are +defined within the struct. For other number formats `threshold_partially_wet` must be provided. The bottom topography function ``b(x,y)`` is set inside the initial condition routine for a particular problem setup. To test the conservative form of the SWE one can set the bottom topography @@ -65,6 +69,10 @@ struct ShallowWaterEquationsWetDry2D{RealT <: Real} <: # before calculating the numerical flux. # Default is 5*eps() which in double precision is ≈1e-15. threshold_wet::RealT + # `threshold_partially_wet` used in `IndicatorHennemannGassnerShallowWater` on the water height + # to define "partially wet" elements. Those elements are calculated with a pure FV method to + # ensure well-balancedness. Default in double precision is 1e-4. + threshold_partially_wet::RealT # Standard shallow water equations for dispatch on Trixi.jl functions basic_swe::ShallowWaterEquations2D{RealT} end @@ -76,7 +84,8 @@ end # Strict default values for thresholds that performed well in many numerical experiments function ShallowWaterEquationsWetDry2D(; gravity_constant, H0 = zero(gravity_constant), threshold_limiter = nothing, - threshold_wet = nothing) + threshold_wet = nothing, + threshold_partially_wet = nothing) T = promote_type(typeof(gravity_constant), typeof(H0)) if threshold_limiter === nothing threshold_limiter = 500 * eps(T) @@ -84,13 +93,16 @@ function ShallowWaterEquationsWetDry2D(; gravity_constant, H0 = zero(gravity_con if threshold_wet === nothing threshold_wet = 5 * eps(T) end + if threshold_partially_wet === nothing + threshold_partially_wet = default_threshold_partially_wet(T) + end # Construct the standard SWE for dispatch. Even though the `basic_swe` already store the # gravity constant and the total water height, we store an extra copy in # `ShallowWaterEquationsWetDry2D` for convenience. basic_swe = ShallowWaterEquations2D(gravity_constant = gravity_constant, H0 = H0) ShallowWaterEquationsWetDry2D(gravity_constant, H0, threshold_limiter, - threshold_wet, basic_swe) + threshold_wet, threshold_partially_wet, basic_swe) end Trixi.have_nonconservative_terms(::ShallowWaterEquationsWetDry2D) = True() diff --git a/src/solvers/indicators_1d.jl b/src/solvers/indicators_1d.jl index 0bf149e..5931e0e 100644 --- a/src/solvers/indicators_1d.jl +++ b/src/solvers/indicators_1d.jl @@ -45,7 +45,7 @@ function (indicator_hg::IndicatorHennemannGassnerShallowWater)(u::AbstractArray{ can only be proven for the FV method (see Chen and Noelle). Therefore we set alpha to one regardless of its given maximum value. =# - threshold_partially_wet = 1e-4 + threshold_partially_wet = equations.threshold_partially_wet Trixi.@threaded for element in eachelement(dg, cache) indicator = indicator_threaded[Threads.threadid()] diff --git a/src/solvers/indicators_2d.jl b/src/solvers/indicators_2d.jl index 1e36414..db9295e 100644 --- a/src/solvers/indicators_2d.jl +++ b/src/solvers/indicators_2d.jl @@ -42,7 +42,7 @@ function (indicator_hg::IndicatorHennemannGassnerShallowWater)(u::AbstractArray{ # Well-balancedness of the scheme on partially wet elements with hydrostatic reconstruction # can only be proven for the FV method (see Chen and Noelle). # Therefore we set alpha to be one regardless of its given value from the modal indicator. - threshold_partially_wet = 1e-4 + threshold_partially_wet = equations.threshold_partially_wet Trixi.@threaded for element in eachelement(dg, cache) indicator = indicator_threaded[Threads.threadid()] From 097492d4204e9c91047e176f48ffba3b58616643 Mon Sep 17 00:00:00 2001 From: patrickersing Date: Mon, 22 Apr 2024 17:14:29 +0200 Subject: [PATCH 13/17] fix typo --- src/equations/equations.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/equations/equations.jl b/src/equations/equations.jl index b92a6fc..9c025bb 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -42,7 +42,7 @@ Retrieve the number of layers from an equation instance of the `AbstractShallowW NLAYERS end -# Provide default thresholds dependend on number format (Currently default thresholds are only provided +# Provide default thresholds dependent on number format (Currently default thresholds are only provided # for Float64) default_threshold_partially_wet(::Type{Float64}) = 1e-4 default_threshold_partially_wet(catchall) = throw(ArgumentError("threshold_partially_wet must be provided for non-Float64 types")) From 039e9d0edecc8bc6e88a5f1d4f16539906d09064 Mon Sep 17 00:00:00 2001 From: patrickersing Date: Mon, 22 Apr 2024 17:28:28 +0200 Subject: [PATCH 14/17] add unit test --- src/equations/equations.jl | 1 + test/test_unit.jl | 5 +++++ 2 files changed, 6 insertions(+) diff --git a/src/equations/equations.jl b/src/equations/equations.jl index 9c025bb..09ff74e 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -42,6 +42,7 @@ Retrieve the number of layers from an equation instance of the `AbstractShallowW NLAYERS end +# TODO: Add suitable default thresholds for Float32 # Provide default thresholds dependent on number format (Currently default thresholds are only provided # for Float64) default_threshold_partially_wet(::Type{Float64}) = 1e-4 diff --git a/test/test_unit.jl b/test/test_unit.jl index 0482e16..ae0e69d 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -239,6 +239,11 @@ end rhos = (0.1, 0.2, 0.3, 0.4)) @test_throws ArgumentError initial_condition_convergence_test(0.0, 0.0, equations) end + +@timed_testset "Exception check for default_threshold functions" begin + @test_throws ArgumentError TrixiShallowWater.default_threshold_partially_wet(Int64) + @test_throws ArgumentError TrixiShallowWater.default_threshold_desingularization(Int64) +end end # Unit tests end # module From 5bda15fc4bae171ac74b9c07e3ef4d2a00783fc4 Mon Sep 17 00:00:00 2001 From: patrickersing Date: Tue, 23 Apr 2024 09:31:35 +0200 Subject: [PATCH 15/17] increase maxiters for single coverage test --- test/test_tree_1d.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/test/test_tree_1d.jl b/test/test_tree_1d.jl index 8323a71..d60dc77 100644 --- a/test/test_tree_1d.jl +++ b/test/test_tree_1d.jl @@ -811,7 +811,8 @@ end # 2LSWE 1.6653345369377348e-16, ], tspan=(0.0, 0.25), - atol=1e-9) + atol=1e-9, + coverage_override=(maxiters = 15,)) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let From 2c597909020e1d9009397801a2d82c813cf21c4a Mon Sep 17 00:00:00 2001 From: Patrick Ersing <114223904+patrickersing@users.noreply.github.com> Date: Sat, 27 Apr 2024 09:54:54 +0200 Subject: [PATCH 16/17] adjust comment Co-authored-by: Andrew Winters --- src/equations/shallow_water_multilayer_1d.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/equations/shallow_water_multilayer_1d.jl b/src/equations/shallow_water_multilayer_1d.jl index 5fe9a90..92fa652 100644 --- a/src/equations/shallow_water_multilayer_1d.jl +++ b/src/equations/shallow_water_multilayer_1d.jl @@ -392,10 +392,10 @@ end """ hydrostatic_reconstruction_ersing_etal(u_ll, u_rr, equations::ShallowWaterMultiLayerEquations1D) -A particular type of hydrostatic reconstruction of the water height and bottom topography to +A particular hydrostatic reconstruction of the water height and bottom topography to guarantee well-balancedness in the presence of wet/dry transitions and entropy stability for the [`ShallowWaterMultiLayerEquations1D`](@ref). -The reconstructed solution states `u_ll_star` and `u_rr_star` variables are used to evaluate the +The reconstructed solution states `u_ll_star` and `u_rr_star` are used to evaluate the surface numerical flux at the interface. The key idea is a piecewise linear reconstruction of the bottom topography and water height interfaces using subcells, where the bottom topography is allowed to be discontinuous. From dcb6abf00b0d66b48e362cf7b0729118096d4af7 Mon Sep 17 00:00:00 2001 From: patrickersing Date: Tue, 30 Apr 2024 14:52:48 +0200 Subject: [PATCH 17/17] update convergence test values --- test/test_tree_1d.jl | 52 ++++++++++++++++++++++---------------------- 1 file changed, 26 insertions(+), 26 deletions(-) diff --git a/test/test_tree_1d.jl b/test/test_tree_1d.jl index 9ccee9e..d7cde02 100644 --- a/test/test_tree_1d.jl +++ b/test/test_tree_1d.jl @@ -512,21 +512,21 @@ end # 2LSWE @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_multilayer_convergence.jl"), l2=[ - 0.0007413378395106811, - 0.0005405760373610748, - 8.772497321219398e-5, - 0.0006969915241438401, - 0.0004924779641660143, - 0.00011816455915431224, + 0.0007413423106723253, + 0.000540576469967569, + 8.772494905937258e-5, + 0.0006969877681271723, + 0.0004924792705658947, + 0.00011816401639126746, 0.000139126377287091, ], linf=[ - 0.002657085141947846, - 0.0028102251075707296, - 0.00022855706065982861, - 0.003310637567035979, - 0.002490976099376596, - 0.0005211086991627756, + 0.0026570823889522366, + 0.0028102271323960926, + 0.0002285569562520129, + 0.003310623577850169, + 0.0024909798619617285, + 0.0005211057995907487, 0.00021874455861881081, ], surface_flux=(flux_lax_friedrichs, @@ -546,22 +546,22 @@ end # 2LSWE @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_multilayer_convergence.jl"), l2=[ - 0.0019005747730996517, - 0.0006416926161516576, - 0.0013237483460263738, - 0.0002514843620148158, - 0.00019551931612091513, - 0.00029396289378303574, - 0.0001391263772870812, + 0.0007413415478415575, + 0.0005405759470722189, + 8.772503316220759e-5, + 0.0006969870846204186, + 0.0004924774643042689, + 0.00011816402423681591, + 0.000139126377287091, ], linf=[ - 0.0059479250195759725, - 0.0028032845871422407, - 0.0039446424223523735, - 0.0009866378484075788, - 0.000930107485495435, - 0.0012826105967129742, - 0.00021874455861947695, + 0.0026571096352370205, + 0.0028102212809490434, + 0.00022855702706781056, + 0.0033106637858648646, + 0.0024909780761511735, + 0.0005211084155695156, + 0.00021874455861881081, ], surface_flux=(FluxHydrostaticReconstruction(FluxPlusDissipation(flux_ersing_etal, DissipationLocalLaxFriedrichs()),