diff --git a/examples/tree_2d_dgsem/elixir_shallowwater_multilayer_dam_break_dry.jl b/examples/tree_2d_dgsem/elixir_shallowwater_multilayer_dam_break_dry.jl new file mode 100644 index 0000000..5cc5b6f --- /dev/null +++ b/examples/tree_2d_dgsem/elixir_shallowwater_multilayer_dam_break_dry.jl @@ -0,0 +1,163 @@ + +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 = ShallowWaterMultiLayerEquations2D(gravity_constant = 1.0, + rhos = (0.9, 0.95, 1.0)) + +function initial_condition_dam_break(x, t, equations::ShallowWaterMultiLayerEquations2D) + # Bottom topography + b = 1.4 * exp(-10.0 * (x[1]^2 + x[2]^2)) + if x[1] > 0.0 + b += 0.1 + end + + if x[1] < -0.5 + H = [1.0, 0.8, 0.6] + else + H = [b, b, b] + end + + v1 = zero(H) + v2 = zero(H) + + # 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 ShallowWaterMultiLayerEquations2D 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 + + return prim2cons(SVector(H..., v1..., v2..., b), + equations) +end + +initial_condition = initial_condition_dam_break + +boundary_conditions = 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) + +############################################################################### +# Get the TreeMesh and setup a non-periodic mesh with wall boundary conditions + +coordinates_min = (-1.0, -1.0) +coordinates_max = (1.0, 1.0) +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 3, + n_cells_max = 10_000, + periodicity = false) + +# Create the semi discretization object +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_conditions) +############################################################################### +# ODE solver + +tspan = (0.0, 2.0) +ode = semidiscretize(semi, tspan) +############################################################################### +#= +Workaround for TreeMesh2D to set true discontinuities for debugging and testing. +Essentially, this is a slight augmentation of the `compute_coefficients` where the `x` node values +passed here are slightly perturbed in order to set a true discontinuity that avoids the doubled +value of the LGL nodes at a particular element interface. +=# + +# 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), j in eachnode(semi.solver) + x_node = Trixi.get_node_coords(semi.cache.elements.node_coordinates, equations, + semi.solver, i, j, element) + # Changing the node positions passed to the initial condition by the minimum + # amount possible with the current type of floating point numbers allows setting + # discontinuous initial data in a simple way. In particular, a check like `if x < x_jump` + # works if the jump location `x_jump` is at the position of an interface. + if i == 1 && j == 1 # bottom left corner + x_node = SVector(nextfloat(x_node[1]), nextfloat(x_node[2])) + elseif i == 1 && j == nnodes(semi.solver) # top left corner + x_node = SVector(nextfloat(x_node[1]), prevfloat(x_node[2])) + elseif i == nnodes(semi.solver) && j == 1 # bottom right corner + x_node = SVector(prevfloat(x_node[1]), nextfloat(x_node[2])) + elseif i == nnodes(semi.solver) && j == nnodes(semi.solver) # top right corner + x_node = SVector(prevfloat(x_node[1]), prevfloat(x_node[2])) + elseif i == 1 # left boundary + x_node = SVector(nextfloat(x_node[1]), x_node[2]) + elseif j == 1 # bottom boundary + x_node = SVector(x_node[1], nextfloat(x_node[2])) + elseif i == nnodes(semi.solver) # right boundary + x_node = SVector(prevfloat(x_node[1]), x_node[2]) + elseif j == nnodes(semi.solver) # top boundary + x_node = SVector(x_node[1], prevfloat(x_node[2])) + end + + u_node = initial_condition_dam_break(x_node, first(tspan), equations) + Trixi.set_node_vars!(u, u_node, equations, semi.solver, i, j, element) + end +end + +############################################################################### +# Callbacks + +summary_callback = SummaryCallback() + +analysis_interval = 50 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + save_analysis = false, + extra_analysis_integrals = (energy_total, + energy_kinetic, + energy_internal)) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 100, + save_initial_solution = true, + save_final_solution = true) + +stepsize_callback = StepsizeCallback(cfl = 0.3) + +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, + stepsize_callback) + +stage_limiter! = PositivityPreservingLimiterShallowWater(variables = (Trixi.waterheight,)) + +############################################################################### +# run the simulation + +# use a Runge-Kutta method with CFL-based time step +sol = solve(ode, SSPRK43(stage_limiter!); + ode_default_options()..., callback = callbacks, adaptive = false, dt = 1.0); +summary_callback() # print the timer summary diff --git a/examples/unstructured_2d_dgsem/elixir_shallowwater_multilayer_dam_break_dry.jl b/examples/unstructured_2d_dgsem/elixir_shallowwater_multilayer_dam_break_dry.jl new file mode 100644 index 0000000..3241dc1 --- /dev/null +++ b/examples/unstructured_2d_dgsem/elixir_shallowwater_multilayer_dam_break_dry.jl @@ -0,0 +1,182 @@ + +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 = ShallowWaterMultiLayerEquations2D(gravity_constant = 1.0, + rhos = (0.9, 0.95, 1.0)) + +# This test case uses a special work around to setup a truly discontinuous bottom topography +# function and initial condition. First, a dummy initial_condition_dam_break is introduced to create +# the semidiscretization. Then the initial condition is reset with the true discontinuous values +# from initial_condition_discontinuous_dam_break. + +function initial_condition_dam_break(x, t, equations::ShallowWaterMultiLayerEquations2D) + # Bottom topography + b = 1.4 * exp(-10.0 * ((x[1] - sqrt(2) / 2)^2 + (x[2] - sqrt(2) / 2)^2)) + + if x[1] < sqrt(2) / 2 + H = [1.0, 0.8, 0.6] + else + b += 0.1 + H = [b, b, b] + end + + v1 = zero(H) + v2 = zero(H) + + # 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 ShallowWaterMultiLayerEquations2D 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 + + return prim2cons(SVector(H..., v1..., v2..., 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(6) + +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 unstructured quad mesh from a file (downloads the file if not available locally) +mesh_file = Trixi.download("https://gist.githubusercontent.com/andrewwinters5000/8f8cd23df27fcd494553f2a89f3c1ba4/raw/85e3c8d976bbe57ca3d559d653087b0889535295/mesh_alfven_wave_with_twist_and_flip.mesh", + joinpath(@__DIR__, "mesh_alfven_wave_with_twist_and_flip.mesh")) + +mesh = UnstructuredMesh2D(mesh_file, periodicity = false) + +# Boundary conditions +boundary_condition = Dict(:Top => boundary_condition_slip_wall, + :Left => boundary_condition_slip_wall, + :Right => boundary_condition_slip_wall, + :Bottom => boundary_condition_slip_wall) + +# Create the semi discretization object +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, + solver, boundary_conditions = boundary_condition) + +############################################################################### +# ODE solver + +tspan = (0.0, 2.0) +ode = semidiscretize(semi, tspan) + +############################################################################### +# Workaround to set a discontinuous bottom topography and initial condition. + +# alternative version of the initial conditinon used to setup a truly discontinuous +# test case and initial condition. +# In contrast to the usual signature of initial conditions, this one get passed the +# `element_id` explicitly. In particular, this initial conditions works as intended +# only for the specific mesh loaded above! +function initial_condition_discontinuous_dam_break(x, t, element_id, + equations::ShallowWaterMultiLayerEquations2D) + # Bottom topography + b = 1.4 * exp(-10.0 * ((x[1] - sqrt(2) / 2)^2 + (x[2] - sqrt(2) / 2)^2)) + + # Left side of discontinuity + IDs = [1, 2, 5, 6, 9, 10, 13, 14] + if element_id in IDs + H = [1.0, 0.8, 0.6] + else # Right side of discontinuity + b += 0.1 + H = [b, b, b] + end + + v1 = zero(H) + v2 = zero(H) + + # 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 ShallowWaterMultiLayerEquations2D 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 + + return prim2cons(SVector(H..., v1..., v2..., b), + equations) +end + +# point to the data we want to augment +u = Trixi.wrap_array(ode.u0, semi) +# reset the initial condition +for element in eachelement(semi.solver, semi.cache) + for j in eachnode(semi.solver), i in eachnode(semi.solver) + x_node = Trixi.get_node_coords(semi.cache.elements.node_coordinates, equations, + semi.solver, i, j, element) + u_node = initial_condition_discontinuous_dam_break(x_node, first(tspan), element, + equations) + Trixi.set_node_vars!(u, u_node, equations, semi.solver, i, j, element) + end +end + +############################################################################### +# Callbacks + +summary_callback = SummaryCallback() + +analysis_interval = 10 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + save_analysis = false, + extra_analysis_integrals = (energy_total, + energy_kinetic, + energy_internal)) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 100, + save_initial_solution = true, + save_final_solution = true) + +stepsize_callback = StepsizeCallback(cfl = 0.5) + +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, + stepsize_callback) + +stage_limiter! = PositivityPreservingLimiterShallowWater(variables = (Trixi.waterheight,)) + +############################################################################### +# run the simulation + +sol = solve(ode, SSPRK43(stage_limiter!); + ode_default_options()..., callback = callbacks, adaptive = false, dt = 1.0); +summary_callback() # print the timer summary diff --git a/examples/unstructured_2d_dgsem/elixir_shallowwater_multilayer_well_balanced_wet_dry.jl b/examples/unstructured_2d_dgsem/elixir_shallowwater_multilayer_well_balanced_wet_dry.jl new file mode 100644 index 0000000..b65f16a --- /dev/null +++ b/examples/unstructured_2d_dgsem/elixir_shallowwater_multilayer_well_balanced_wet_dry.jl @@ -0,0 +1,179 @@ + +using Downloads: download +using OrdinaryDiffEq +using Trixi +using TrixiShallowWater + +############################################################################### +# Semidiscretization of the multilayer shallow water equations with a discontinuous +# bottom topography function (set in the initial conditions) and dry domains + +equations = ShallowWaterMultiLayerEquations2D(gravity_constant = 9.81, H0 = 1.5, + rhos = (0.9, 0.95, 1.0)) + +# An initial condition with constant total water height and zero velocities to test well-balancedness. +# Note, this routine is used to compute errors in the analysis callback but the initialization is +# overwritten by `initial_condition_discontinuous_well_balancedness` below. +function initial_condition_well_balancedness(x, t, + equations::ShallowWaterMultiLayerEquations2D) + # Set the background values + H = [1.5, 1.0, 0.5] + v1 = zero(H) + v2 = zero(H) + + # bottom topography taken from Pond.control in [HOHQMesh](https://github.com/trixi-framework/HOHQMesh) + x1, x2 = x + b = (1.5 / exp(0.5 * ((x1 - 1.0)^2 + (x2 - 1.0)^2)) + + 0.75 / exp(0.5 * ((x1 + 1.0)^2 + (x2 + 1.0)^2))) + + # 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 ShallowWaterMultiLayerEquations2D 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 + + return prim2cons(SVector(H..., v1..., v2..., b), equations) +end + +initial_condition = initial_condition_well_balancedness + +############################################################################### +# 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(6) + +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) + +############################################################################### +# This setup is for the curved, split form well-balancedness testing + +# Get the unstructured quad mesh from a file (downloads the file if not available locally) +default_mesh_file = joinpath(@__DIR__, "mesh_alfven_wave_with_twist_and_flip.mesh") +isfile(default_mesh_file) || + download("https://gist.githubusercontent.com/andrewwinters5000/8f8cd23df27fcd494553f2a89f3c1ba4/raw/85e3c8d976bbe57ca3d559d653087b0889535295/mesh_alfven_wave_with_twist_and_flip.mesh", + default_mesh_file) +mesh_file = default_mesh_file + +mesh = UnstructuredMesh2D(mesh_file, periodicity = true) + +# Create the semi discretization object +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +############################################################################### +# ODE solver + +tspan = (0.0, 1.0) +ode = semidiscretize(semi, tspan) + +############################################################################### +# Workaround to set a discontinuous bottom topography for debugging and testing. + +# alternative version of the initial conditinon used to setup a truly discontinuous +# bottom topography function for this academic testcase. +# The errors from the analysis callback are not important but the error for this lake at rest test case +# `∑|H0-(h+b)|` should be around machine roundoff +# In contrast to the usual signature of initial conditions, this one get passed the +# `element_id` explicitly. In particular, this initial conditions works as intended +# only for the specific mesh loaded above! +function initial_condition_discontinuous_well_balancedness(x, t, element_id, + equations::ShallowWaterMultiLayerEquations2D) + # Set the background values + H = [1.5, 1.0, 0.5] + v1 = zero(H) + v2 = zero(H) + + b = 0.2 + 0.1 * sin(2.0 * pi * x[1]) + 0.1 * cos(2.0 * pi * x[2]) + + # Setup a discontinuous bottom topography to create wet/dry transitions using the `element_id` number + if element_id == 7 + b += 0.1 + elseif element_id == 6 + b += 0.5 + elseif element_id == 10 + b += 1.0 + elseif element_id == 11 + b += 1.5 + end + + # It is mandatory to shift the water level at dry areas to make sure the water height h + # stays positive. The system would not be stable for h set to a hard 0 due to division by h in + # the computation of velocity, e.g., (h v) / h. Therefore, a small dry state threshold + # 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 + + return prim2cons(SVector(H..., v1..., v2..., b), equations) +end + +# point to the data we want to augment +u = Trixi.wrap_array(ode.u0, semi) +# reset the initial condition +for element in eachelement(semi.solver, semi.cache) + for j in eachnode(semi.solver), i in eachnode(semi.solver) + x_node = Trixi.get_node_coords(semi.cache.elements.node_coordinates, equations, + semi.solver, i, j, element) + u_node = initial_condition_discontinuous_well_balancedness(x_node, first(tspan), + element, equations) + Trixi.set_node_vars!(u, u_node, equations, semi.solver, i, j, element) + end +end + +############################################################################### +# Callbacks + +summary_callback = SummaryCallback() + +analysis_interval = 500 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + extra_analysis_integrals = (lake_at_rest_error,)) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 1000, + save_initial_solution = true, + save_final_solution = true) + +stepsize_callback = StepsizeCallback(cfl = 1.0) + +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, + stepsize_callback) + +stage_limiter! = PositivityPreservingLimiterShallowWater(variables = (Trixi.waterheight,)) + +############################################################################### +# run the simulation + +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 8ed1630..80a5307 100644 --- a/src/callbacks_stage/positivity_shallow_water.jl +++ b/src/callbacks_stage/positivity_shallow_water.jl @@ -81,7 +81,8 @@ function limiter_shallow_water!(u, variables::NTuple{N, Any}, mesh, equations::Union{ShallowWaterEquationsWetDry1D, ShallowWaterEquationsWetDry2D, - ShallowWaterMultiLayerEquations1D}, + ShallowWaterMultiLayerEquations1D, + ShallowWaterMultiLayerEquations2D}, solver, cache) where {N} variable = first(variables) remaining_variables = Base.tail(variables) @@ -97,7 +98,8 @@ function limiter_shallow_water!(u, variables::Tuple{}, mesh, equations::Union{ShallowWaterEquationsWetDry1D, ShallowWaterEquationsWetDry2D, - ShallowWaterMultiLayerEquations1D}, + ShallowWaterMultiLayerEquations1D, + ShallowWaterMultiLayerEquations2D}, solver, cache) nothing end diff --git a/src/callbacks_stage/positivity_shallow_water_dg2d.jl b/src/callbacks_stage/positivity_shallow_water_dg2d.jl index 69c47e1..9acdc4d 100644 --- a/src/callbacks_stage/positivity_shallow_water_dg2d.jl +++ b/src/callbacks_stage/positivity_shallow_water_dg2d.jl @@ -88,4 +88,86 @@ function limiter_shallow_water!(u, threshold::Real, variable, return nothing end + +# !!! warning "Experimental code" +# This is an experimental feature and may change in future releases. +function limiter_shallow_water!(u, threshold::Real, variable, + mesh::Trixi.AbstractMesh{2}, + equations::ShallowWaterMultiLayerEquations2D, 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 j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, j, 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, 1, element)) + for j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, j, element) + u_mean += u_node * weights[i] * weights[j] + 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 j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, j, element) + h_node = waterheight(u_node, equations)[m] + h_mean = waterheight(u_mean, equations)[m] + + u[m, i, j, 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 + Trixi.@threaded for element in eachelement(dg, cache) + for j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, j, element) + + h = MVector(waterheight(u_node, equations)) + h_v1 = MVector(momentum(u_node, equations)[1]) + h_v2 = MVector(momentum(u_node, equations)[2]) + b = u_node[end] + + # Velocity desingularization + h_v1 = h .* (2.0 * h .* h_v1) ./ + (h .^ 2 .+ max.(h .^ 2, equations.threshold_desingularization)) + h_v2 = h .* (2.0 * h .* h_v2) ./ + (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 + h_v1[i] = zero(eltype(u)) + h_v2[i] = zero(eltype(u)) + end + end + + u_node = SVector(h..., h_v1..., h_v2..., b) + + set_node_vars!(u, u_node, equations, dg, i, j, element) + end + end + + return nothing +end end # @muladd diff --git a/src/equations/numerical_fluxes.jl b/src/equations/numerical_fluxes.jl index a523388..1d282fa 100644 --- a/src/equations/numerical_fluxes.jl +++ b/src/equations/numerical_fluxes.jl @@ -28,4 +28,22 @@ For complete details see Section 2.4 of the following reference [DOI: 10.1137/15M1053074](https://doi.org/10.1137/15M1053074) """ const flux_hll_chen_noelle = FluxHLL(min_max_speed_chen_noelle) + +# Additional version of `FluxHydrostaticReconstruction` to add support for nonconservative fluxes on +# unstructured meshes. These can depend on both the contravariant vectors (normal direction) at +# the current node and the averaged ones. +@inline function (numflux::Trixi.FluxHydrostaticReconstruction)(u_ll, u_rr, + normal_direction_ll, + normal_direction_average, + equations::Trixi.AbstractEquations) + @unpack numerical_flux, hydrostatic_reconstruction = numflux + + # Create the reconstructed left/right solution states in conservative form + u_ll_star, u_rr_star = hydrostatic_reconstruction(u_ll, u_rr, equations) + + # Use the reconstructed states to compute the nonconservative surface flux + return numerical_flux(u_ll_star, u_rr_star, normal_direction_ll, + normal_direction_average, + equations) +end end diff --git a/src/equations/shallow_water_multilayer_2d.jl b/src/equations/shallow_water_multilayer_2d.jl index da26309..b541b0b 100644 --- a/src/equations/shallow_water_multilayer_2d.jl +++ b/src/equations/shallow_water_multilayer_2d.jl @@ -31,6 +31,16 @@ nonconservative term, which has some benefits for the design of well-balanced sc The additional quantity ``H_0`` is also available to store a reference value for the total water height that is useful to set initial conditions or test the "lake-at-rest" well-balancedness. +Also, there are two thresholds which prevent numerical problems as well as instabilities. The limiters are +used in [`PositivityPreservingLimiterShallowWater`](@ref) on the water height. `threshold_limiter` +acts as a (small) shift on the initial condition and cutoff before the next time step, whereas +`threshold_desingularization` is used in the velocity desingularization. A third +`threshold_partially_wet` is applied on the water height to define "partially wet" elements in +[`IndicatorHennemannGassnerShallowWater`](@ref), that are then calculated with a pure FV method to +ensure well-balancedness. For `Float64` no threshold needs to be passed, as default values are +defined within the struct. For other number formats `threshold_desingularization` and `threshold_partially_wet` +must be provided. + The bottom topography function ``b(x)`` is set inside the initial condition routine for a particular problem setup. @@ -55,10 +65,16 @@ struct ShallowWaterMultiLayerEquations2D{NVARS, NLAYERS, RealT <: Real} <: AbstractShallowWaterMultiLayerEquations{2, NVARS, NLAYERS} gravity::RealT # gravitational constant H0::RealT # constant "lake-at-rest" total water height + threshold_limiter::RealT # threshold for the positivity-limiter + threshold_desingularization::RealT # threshold for velocity desingularization + threshold_partially_wet::RealT # threshold to define partially wet elements rhos::SVector{NLAYERS, RealT} # Vector of layer densities function ShallowWaterMultiLayerEquations2D{NVARS, NLAYERS, RealT}(gravity::RealT, H0::RealT, + threshold_limiter::RealT, + threshold_desingularization::RealT, + threshold_partially_wet::RealT, rhos::SVector{NLAYERS, RealT}) where { NVARS, @@ -71,7 +87,8 @@ struct ShallowWaterMultiLayerEquations2D{NVARS, NLAYERS, RealT <: Real} <: throw(ArgumentError("densities must be in increasing order (rhos[1] < rhos[2] < ... < rhos[NLAYERS])")) min(rhos...) > 0 || throw(ArgumentError("densities must be positive")) - new(gravity, H0, rhos) + new(gravity, H0, threshold_limiter, threshold_desingularization, + threshold_partially_wet, rhos) end end @@ -80,19 +97,38 @@ end # The reference total water height H0 defaults to 0.0 but is used for the "lake-at-rest" # well-balancedness test cases. function ShallowWaterMultiLayerEquations2D(; gravity_constant, - H0 = zero(gravity_constant), rhos) + H0 = zero(gravity_constant), + threshold_limiter = nothing, + threshold_desingularization = nothing, + threshold_partially_wet = nothing, + rhos) # Promote all variables to a common type _rhos = promote(rhos...) RealT = promote_type(eltype(_rhos), eltype(gravity_constant), eltype(H0)) __rhos = SVector(map(RealT, _rhos)) + # Set default values for thresholds + if threshold_limiter === nothing + threshold_limiter = 5 * eps(RealT) + end + if threshold_desingularization === nothing + threshold_desingularization = default_threshold_desingularization(RealT) + end + if threshold_partially_wet === nothing + threshold_partially_wet = default_threshold_partially_wet(RealT) + end + # Extract number of layers and variables NLAYERS = length(rhos) NVARS = 3 * NLAYERS + 1 return ShallowWaterMultiLayerEquations2D{NVARS, NLAYERS, RealT}(gravity_constant, - H0, __rhos) + H0, + threshold_limiter, + threshold_desingularization, + threshold_partially_wet, + __rhos) end @inline function Base.real(::ShallowWaterMultiLayerEquations2D{NVARS, NLAYERS, RealT}) where { @@ -505,7 +541,7 @@ end f_hv * normal_direction_average[2], i, equations) end - return (SVector(f)) + return SVector(f) end """ @@ -601,6 +637,85 @@ end return SVector(f) end +""" + hydrostatic_reconstruction_ersing_etal(u_ll, u_rr, equations::ShallowWaterMultiLayerEquations2D) + +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 +[`ShallowWaterMultiLayerEquations2D`](@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). + +!!! 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::ShallowWaterMultiLayerEquations2D) + # Unpack waterheight and bottom topographies + h_ll = MVector(waterheight(u_ll, equations)) + h_rr = MVector(waterheight(u_rr, equations)) + b_ll = u_ll[end] + b_rr = u_rr[end] + + # Get the velocities on either side + v1_ll = MVector(velocity(u_ll, equations)[1]) + v2_ll = MVector(velocity(u_ll, equations)[2]) + v1_rr = MVector(velocity(u_rr, equations)[1]) + v2_rr = MVector(velocity(u_rr, equations)[2]) + + threshold = equations.threshold_limiter + + # Ensure zero velocity at dry states + for i in eachlayer(equations) + if h_ll[i] <= threshold + v1_ll[i] = zero(eltype(u_ll)) + v2_ll[i] = zero(eltype(u_ll)) + end + if h_rr[i] <= threshold + v1_rr[i] = zero(eltype(u_rr)) + v2_rr[i] = zero(eltype(u_rr)) + end + end + + # Calculate total layer heights + H_ll = waterheight(cons2prim(u_ll, equations), equations) + H_rr = waterheight(cons2prim(u_rr, equations), equations) + + # Discontinuous reconstruction of the bottom topography + b_ll_star = min(H_ll[1], max(b_rr, b_ll)) + b_rr_star = min(H_rr[1], max(b_rr, b_ll)) + + # Calculate reconstructed total layer heights + H_ll_star = max.(H_ll, b_ll_star) + H_rr_star = max.(H_rr, b_rr_star) + + # Initialize reconstructed waterheights + h_ll_star = zero(MVector{nlayers(equations), real(equations)}) + h_rr_star = zero(MVector{nlayers(equations), real(equations)}) + + # Reconstruct the waterheights + for i in eachlayer(equations) + if i == nlayers(equations) # The lowest layer is measured from the bottom topography + h_ll_star[i] = max(H_ll_star[i] - b_ll_star, threshold) + h_rr_star[i] = max(H_rr_star[i] - b_rr_star, threshold) + else + h_ll_star[i] = max(H_ll_star[i] - H_ll_star[i + 1], threshold) + h_rr_star[i] = max(H_rr_star[i] - H_rr_star[i + 1], threshold) + end + end + + # Reconstructed states + u_ll_star = SVector(h_ll_star..., (h_ll_star .* v1_ll)..., (h_ll_star .* v2_ll)..., + b_ll_star) + u_rr_star = SVector(h_rr_star..., (h_rr_star .* v1_rr)..., (h_rr_star .* v2_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, @@ -834,14 +949,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::ShallowWaterMultiLayerEquations2D) + 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::ShallowWaterMultiLayerEquations2D) 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_2d.jl b/src/solvers/indicators_2d.jl index db9295e..833c30c 100644 --- a/src/solvers/indicators_2d.jl +++ b/src/solvers/indicators_2d.jl @@ -5,13 +5,14 @@ @muladd begin #! format: noindent -# Modified indicator for ShallowWaterEquationsWetDry2D 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 ShallowWaterEquationsWetDry2D and ShallowWaterMultiLayerEquations2D 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, 4}, mesh, - equations::ShallowWaterEquationsWetDry2D, + equations::Union{ShallowWaterEquationsWetDry2D, + ShallowWaterMultiLayerEquations2D}, dg::DGSEM, cache; kwargs...) @unpack alpha_max, alpha_min, alpha_smooth, variable = indicator_hg @@ -55,9 +56,10 @@ function (indicator_hg::IndicatorHennemannGassnerShallowWater)(u::AbstractArray{ # Calculate indicator variables at Gauss-Lobatto nodes for j in eachnode(dg), i in eachnode(dg) u_local = get_node_vars(u, equations, dg, i, j, element) - h, _, _, _ = u_local + h = waterheight(u_local, equations) - if h <= threshold_partially_wet + # Set indicator to FV if water height is below the threshold + if minimum(h) <= threshold_partially_wet indicator_wet = 0 end diff --git a/test/test_tree_2d.jl b/test/test_tree_2d.jl index 18a780a..c52c6ea 100644 --- a/test/test_tree_2d.jl +++ b/test/test_tree_2d.jl @@ -523,6 +523,48 @@ end # 2LSWE end end + @trixi_testset "elixir_shallowwater_multilayer_convergence.jl with FluxHydrostaticReconstruction" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_multilayer_convergence.jl"), + l2=[ + 0.0015784809121655386, + 0.0010725791806396989, + 0.000361175462962745, + 0.001010776670918376, + 0.0009612673392275332, + 0.0002639057696207499, + 0.0011407118095590378, + 0.0011974691247903164, + 0.00031384421384170354, + 0.00019675440964325044, + ], + linf=[ + 0.006611273262456141, + 0.004798759663784402, + 0.0014761631666105335, + 0.004378481797736367, + 0.004059620398099983, + 0.001072537228128334, + 0.004650692998510397, + 0.00501486495741188, + 0.0013175223929257074, + 0.0004374891172380657, + ], + surface_flux=(FluxHydrostaticReconstruction(flux_ersing_etal, + 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"), @@ -690,6 +732,47 @@ end # 2LSWE @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 end end + + @trixi_testset "elixir_shallowwater_multilayer_dam_break_dry.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_multilayer_dam_break_dry.jl"), + l2=[ + 0.030684387782428005, + 0.03094529104416787, + 0.09023105191044616, + 0.021006479937144513, + 0.02094556081555727, + 0.06145762035878689, + 0.0010401307741441598, + 0.0010346793974473935, + 0.002952821531601029, + 0.005470808030402103, + ], + linf=[ + 0.10807087634534615, + 0.11202631201715663, + 0.33977650047095415, + 0.05694587797245012, + 0.056598515593249694, + 0.17171444846932465, + 0.00827507857831736, + 0.008095427727762589, + 0.026626783423061535, + 0.1016120899921184, + ], + tspan=(0.0, 0.25), + # Increase iterations for coverage testing to trigger the + # positivity limiter + coverage_override=(maxiters = 130, tspan = (0.0, 1.5))) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end end # MLSWE end # TreeMesh2D diff --git a/test/test_unstructured_2d.jl b/test/test_unstructured_2d.jl index 8bdd0bd..03e36f5 100644 --- a/test/test_unstructured_2d.jl +++ b/test/test_unstructured_2d.jl @@ -506,6 +506,48 @@ end # 2LSWE end end + @trixi_testset "elixir_shallowwater_multilayer_convergence.jl with FluxHydrostaticReconstruction" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_multilayer_convergence.jl"), + l2=[ + 2.8451070677919015e-5, + 2.5413328287320148e-5, + 1.3776844295056111e-5, + 3.216416098536785e-5, + 1.964525963539178e-5, + 1.2494077535320725e-5, + 3.842168759685933e-5, + 2.3940963660143302e-5, + 1.553940305228399e-5, + 7.244891628173627e-7, + ], + linf=[ + 0.00016865976829461005, + 0.0001694347594369816, + 0.00010161333860081445, + 0.00018234203943068295, + 0.00014403358429010416, + 7.821014428416317e-5, + 0.0002159224956841399, + 0.00016273539534478187, + 9.798918562298198e-5, + 4.040896422807805e-6, + ], + surface_flux=(FluxHydrostaticReconstruction(flux_ersing_etal, + 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"), @@ -596,6 +638,44 @@ 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.258337877456028, + 0.4536334309305916, + 0.2424975342044823, + 1.861822857108904e-15, + 1.0728615421620665e-15, + 4.893078129679801e-16, + 2.951473837608095e-15, + 9.819198357343179e-16, + 3.666055853677385e-16, + 0.9049896638396252, + ], + linf=[ + 0.5000000000000003, + 0.5000000000000064, + 0.485509867641643, + 1.7168360326077987e-14, + 6.4508928741529934e-15, + 3.663826180637297e-15, + 2.370261384584603e-14, + 6.3907606640361836e-15, + 2.0890952713058902e-15, + 1.3458935664973586, + ], + 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_dam_break.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_multilayer_dam_break.jl"), @@ -673,6 +753,44 @@ end # 2LSWE @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 end end + + @trixi_testset "elixir_shallowwater_multilayer_dam_break_dry.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_multilayer_dam_break_dry.jl"), + l2=[ + 0.025056397690317766, + 0.02413438932652429, + 0.05305193031437963, + 0.01570655272565154, + 0.015434075201987807, + 0.03802395905825049, + 0.0035119966268674424, + 0.003140567605820363, + 0.0047903431628184764, + 0.00400367487891041, + ], + linf=[ + 0.13756514254828095, + 0.13141572813228075, + 0.40279814892528143, + 0.05125099836463482, + 0.05082705371717001, + 0.20100787456634925, + 0.021397757474313436, + 0.019301532073029846, + 0.031247259773652267, + 0.10003205938749304, + ], + 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 # UnstructuredMesh2D