From c82880653c7339dd1097951fc75f79fb3e86a2df Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Tue, 13 May 2025 15:00:37 +0200 Subject: [PATCH 01/37] initial commit. Well-balanced working on fully wet mortars --- src/callbacks_step/amr_dg2d.jl | 4 +- src/equations/shallow_water_multilayer_2d.jl | 58 ++++--- src/solvers/dgsem_p4est/containers.jl | 2 +- src/solvers/dgsem_p4est/dg_2d.jl | 165 ++++++++++++++++++- 4 files changed, 194 insertions(+), 35 deletions(-) diff --git a/src/callbacks_step/amr_dg2d.jl b/src/callbacks_step/amr_dg2d.jl index 74938205..c10e0ede 100644 --- a/src/callbacks_step/amr_dg2d.jl +++ b/src/callbacks_step/amr_dg2d.jl @@ -19,7 +19,7 @@ # !!! warning "Experimental code" # This is an experimental feature and may change in future releases. function Trixi.refine!(u_ode::AbstractVector, adaptor, mesh::P4estMesh{2}, - equations::ShallowWaterEquationsWetDry2D, + equations::Union{ShallowWaterEquationsWetDry2D,ShallowWaterMultiLayerEquations2D}, dg::DGSEM, cache, elements_to_refine) # Return early if there is nothing to do if isempty(elements_to_refine) @@ -123,7 +123,7 @@ end # !!! warning "Experimental code" # This is an experimental feature and may change in future releases. function Trixi.coarsen!(u_ode::AbstractVector, adaptor, mesh::P4estMesh{2}, - equations::ShallowWaterEquationsWetDry2D, + equations::Union{ShallowWaterEquationsWetDry2D,ShallowWaterMultiLayerEquations2D}, dg::DGSEM, cache, elements_to_remove) # Return early if there is nothing to do if isempty(elements_to_remove) diff --git a/src/equations/shallow_water_multilayer_2d.jl b/src/equations/shallow_water_multilayer_2d.jl index 2012ea39..1103ab38 100644 --- a/src/equations/shallow_water_multilayer_2d.jl +++ b/src/equations/shallow_water_multilayer_2d.jl @@ -11,7 +11,7 @@ Multi-Layer Shallow Water equations (MLSWE) in two space dimension. The equations are given by ```math \left\{ - \begin{aligned} + \begin{aligned} &\partial_t h_m + \partial_x h_mv_m = 0,\\ &\partial h_mv1_m + \partial_x h_mv1_m^2 + \partial_y h_mv1_mv2_m = -gh_m\partial_x \bigg(b + \sum\limits_{k\geq j}h_k + \sum\limits_{k= + 2 * (equations.threshold_limiter) + u_buffer[1, i] = u[1, i_large, j_large, element] + + u[4, i_large, j_large, element] + else + u_buffer[1, i] = equations.H0 + end + for v in 2:4 + u_buffer[v, i] = u[v, i_large, j_large, element] + end + + i_large += i_large_step + j_large += j_large_step + end + + # Interpolate large element face data from buffer to small face locations + Trixi.multiply_dimensionwise!(view(cache.mortars.u, 2, :, 1, :, mortar), + mortar_l2.forward_lower, + u_buffer) + Trixi.multiply_dimensionwise!(view(cache.mortars.u, 2, :, 2, :, mortar), + mortar_l2.forward_upper, + u_buffer) + + # After the projection of the constant solution we modify the values + # in the first solution variable to no longer be the total water height H = h+b + # and instead recover the conservative water height variable `h`. + # Basically, unpacking the total water height variable to create the projected local water + # height `h` from Eq. 41 in Benov et al. + for i in eachnode(dg) + cache.mortars.u[2, 1, 1, i, mortar] = max(cache.mortars.u[2, 1, 1, i, + mortar] - + cache.mortars.u[2, 4, 1, i, + mortar], + equations.threshold_limiter) + cache.mortars.u[2, 1, 2, i, mortar] = max(cache.mortars.u[2, 1, 2, i, + mortar] - + cache.mortars.u[2, 4, 2, i, + mortar], + equations.threshold_limiter) + + # Safety application of velocity desingularization and water height cutoff on the mortars. + # TODO: Experiment with `threshold_desingularization` and `threshold_limiter` + # + # For details on the motivation of the velocity desingularization see + # - A. Chertock, S. Cui, A. Kurganov, T. Wu (2015) + # Well-balanced positivity preserving central-upwind scheme for + # the shallow water system with friction terms + # [DOI: 10.1002/fld.4023](https://doi.org/10.1002/fld.4023) + + # Mortars with the copied or projected solution + tol = equations.threshold_desingularization + for child in 1:2, side in 1:2 + if cache.mortars.u[side, 1, child, i, mortar] <= + equations.threshold_limiter + cache.mortars.u[side, 1, child, i, mortar] = equations.threshold_limiter + cache.mortars.u[side, 2, child, i, mortar] = zero(eltype(u)) + cache.mortars.u[side, 3, child, i, mortar] = zero(eltype(u)) + + else + h = cache.mortars.u[side, 1, child, i, mortar] + cache.mortars.u[side, 2, child, i, mortar] = h * (2 * h * + cache.mortars.u[side, + 2, + child, + i, + mortar]) / + (h^2 + max(h^2, tol)) + cache.mortars.u[side, 3, child, i, mortar] = h * (2 * h * + cache.mortars.u[side, + 3, + child, + i, + mortar]) / + (h^2 + max(h^2, tol)) + end + end + end + end + + return nothing +end + # !!! warning "Experimental code" # This is an experimental feature and may change in future releases. function Trixi.calc_mortar_flux!(surface_flux_values, - mesh::Union{P4estMesh{2}, T8codeMesh{2}}, + mesh::P4estMesh{2}, nonconservative_terms, equations::ShallowWaterEquationsWetDry2D, mortar_l2::Trixi.LobattoLegendreMortarL2, @@ -291,8 +447,7 @@ end # !!! warning "Experimental code" # This is an experimental feature and may change in future releases. @inline function Trixi.mortar_fluxes_to_elements!(surface_flux_values, - mesh::Union{P4estMesh{2}, - T8codeMesh{2}}, + mesh::P4estMesh{2}, equations::ShallowWaterEquationsWetDry2D, mortar_l2::Trixi.LobattoLegendreMortarL2, dg::DGSEM, cache, mortar, From bb25645e0b26758db30d52ea3f6b6ef4cc0a31ab Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Tue, 13 May 2025 15:02:20 +0200 Subject: [PATCH 02/37] forgot to add new elixir --- ...r_shallowwater_multilayer_well_balanced.jl | 228 ++++++++++++++++++ 1 file changed, 228 insertions(+) create mode 100644 examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_well_balanced.jl diff --git a/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_well_balanced.jl b/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_well_balanced.jl new file mode 100644 index 00000000..0052bee3 --- /dev/null +++ b/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_well_balanced.jl @@ -0,0 +1,228 @@ + +using OrdinaryDiffEqSSPRK, OrdinaryDiffEqLowStorageRK +using Trixi +using TrixiShallowWater + +############################################################################### +# Semidiscretization of the shallow water equations with a discontinuous +# bottom topography function and wet/dry transition elements on a nonconforming mesh + +############################################################################### +# Semidiscretization of the multilayer shallow water equations with one layers + +equations = ShallowWaterMultiLayerEquations2D(gravity = 9.812, + # H0 = 2.0, # for fully wet testing + H0 = 1.235, # for wet/dry transition testing + # To have the same thresholds as the `ShallowWater2D` + threshold_desingularization = 1e-4, + threshold_limiter = 500*eps(), + rhos = (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 = [equations.H0] + v1 = zero(H) + v2 = zero(H) + + x1, x2 = x + b = (1.5 / exp(0.5 * ((x1 - 1)^2 + (x2 + 1)^2)) + + + 0.8 / exp(0.5 * ((x1 + 1)^2 + (x2 - 1)^2)) + - + 0.5 / exp(3.5 * ((x1 - 0.4)^2 + (x2 - 0.325)^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 + +boundary_condition = Dict(:all => 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)) + +# Create the solver +basis = LobattoLegendreBasis(3) + +indicator_sc = IndicatorHennemannGassnerShallowWater(equations, basis, + alpha_max = 0.5, + alpha_min = 0.001, + alpha_smooth = true, + # OBS! One cannot simply use `waterheight` + # here for multilayer equations. + # Throws a conversion error. + # Would need another function that returns sum(h) + # after the `waterheight` call. + 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) + +# Unstructured mesh with 24 cells of the square domain [-1, 1]^n +mesh_file = Trixi.download("https://gist.githubusercontent.com/efaulhaber/63ff2ea224409e55ee8423b3a33e316a/raw/7db58af7446d1479753ae718930741c47a3b79b7/square_unstructured_2.inp", + joinpath(@__DIR__, "square_unstructured_2.inp")) + +# Affine type mapping to take the [-1,1]^2 domain from the mesh file +# and warp it as described in https://arxiv.org/abs/2012.12040 +# Warping with the coefficient 0.15 is even more extreme. +function mapping_twist(xi, eta) + y = eta + 0.15 * cos(1.5 * pi * xi) * cos(0.5 * pi * eta) + x = xi + 0.15 * cos(0.5 * pi * xi) * cos(2.0 * pi * y) + return SVector(x, y) +end + +mesh = P4estMesh{2}(mesh_file, polydeg = 3, + mapping = mapping_twist, + initial_refinement_level = 0) + +# Refine bottom left quadrant of each tree to level 3 +function refine_fn(p4est, which_tree, quadrant) + quadrant_obj = unsafe_load(quadrant) + if quadrant_obj.x == 0 && quadrant_obj.y == 0 && quadrant_obj.level < 3 + # return true (refine) + return Cint(1) + else + # return false (don't refine) + return Cint(0) + end +end + +# Refine recursively until each bottom left quadrant of a tree has level 3 +# The mesh will be rebalanced before the simulation starts +refine_fn_c = @cfunction(refine_fn, Cint, + (Ptr{Trixi.p4est_t}, Ptr{Trixi.p4est_topidx_t}, + Ptr{Trixi.p4est_quadrant_t})) +Trixi.refine_p4est!(mesh.p4est, true, refine_fn_c, C_NULL) + +# Create the semi discretization object +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_condition) + +############################################################################### +# ODE solvers, callbacks, etc. + +tspan = (0.0, 50.0) +ode = semidiscretize(semi, tspan) + +############################################################################### +# Workaround to set a discontinuous bottom topography for debugging and testing. +# 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 = [equations.H0] + v1 = zero(H) + v2 = zero(H) + + x1, x2 = x + b = (1.75 / exp(0.6 * ((x1 - 1.0)^2 + (x2 + 1.0)^2)) + + + 0.8 / exp(0.5 * ((x1 + 1.0)^2 + (x2 - 1.0)^2)) + - + 0.5 / exp(3.5 * ((x1 - 0.4)^2 + (x2 - 0.325)^2))) + + # Setup a discontinuous bottom topography using the element id number + IDs = [collect(114:133); collect(138:141); collect(156:164); collect(208:300)] + if element_id in IDs + b = (0.75 / exp(0.5 * ((x1 - 1.0)^2 + (x2 + 1.0)^2)) + + + 0.4 / exp(0.5 * ((x1 + 1.0)^2 + (x2 - 1.0)^2)) + - + 0.25 / exp(3.5 * ((x1 - 0.4)^2 + (x2 - 0.325)^2))) + 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 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_well_balancedness(x_node, first(tspan), + element, equations) + Trixi.set_node_vars!(u, u_node, equations, semi.solver, i, j, element) + end +end +############################################################################### + +summary_callback = SummaryCallback() + +analysis_interval = 10000 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + extra_analysis_errors = (:conservation_error,), + extra_analysis_integrals = (lake_at_rest_error,)) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(dt = 10.0, + save_initial_solution = true, + save_final_solution = true) + +stepsize_callback = StepsizeCallback(cfl = 0.5) + +callbacks = CallbackSet(summary_callback, + analysis_callback, + alive_callback, + save_solution, + stepsize_callback) + +stage_limiter! = PositivityPreservingLimiterShallowWater(variables = (waterheight,)) + +############################################################################### +# run the simulation +sol = solve(ode, SSPRK43(stage_limiter!), + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + adaptive = false, + save_everystep = false, callback = callbacks, maxiters=1e9); From f6ddaaffad5f1c9b3f366452783d763ca14cf48b Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Tue, 13 May 2025 20:31:13 +0200 Subject: [PATCH 03/37] comment out mortar desingularization to experiment --- ...r_shallowwater_multilayer_well_balanced.jl | 5 +- src/solvers/dgsem_p4est/dg_2d.jl | 82 +++++++++---------- 2 files changed, 45 insertions(+), 42 deletions(-) diff --git a/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_well_balanced.jl b/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_well_balanced.jl index 0052bee3..4db9318e 100644 --- a/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_well_balanced.jl +++ b/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_well_balanced.jl @@ -70,6 +70,9 @@ surface_flux = (FluxHydrostaticReconstruction(FluxPlusDissipation(flux_ersing_et # Create the solver basis = LobattoLegendreBasis(3) +# TODO: Write a generic function so that I can use `waterheight` in the limiter. +# That way this part of the elixir is identical to the `ShallowWater2D` counterpart + indicator_sc = IndicatorHennemannGassnerShallowWater(equations, basis, alpha_max = 0.5, alpha_min = 0.001, @@ -206,7 +209,7 @@ analysis_callback = AnalysisCallback(semi, interval = analysis_interval, alive_callback = AliveCallback(analysis_interval = analysis_interval) -save_solution = SaveSolutionCallback(dt = 10.0, +save_solution = SaveSolutionCallback(dt = 0.1, # 10.0, save_initial_solution = true, save_final_solution = true) diff --git a/src/solvers/dgsem_p4est/dg_2d.jl b/src/solvers/dgsem_p4est/dg_2d.jl index 078b9ef7..d7423ae3 100644 --- a/src/solvers/dgsem_p4est/dg_2d.jl +++ b/src/solvers/dgsem_p4est/dg_2d.jl @@ -290,50 +290,50 @@ function Trixi.prolong2mortars!(cache, u, cache.mortars.u[2, 1, 1, i, mortar] = max(cache.mortars.u[2, 1, 1, i, mortar] - cache.mortars.u[2, 4, 1, i, - mortar], - equations.threshold_limiter) + mortar], 0) + # equations.threshold_limiter) cache.mortars.u[2, 1, 2, i, mortar] = max(cache.mortars.u[2, 1, 2, i, mortar] - cache.mortars.u[2, 4, 2, i, - mortar], - equations.threshold_limiter) - - # Safety application of velocity desingularization and water height cutoff on the mortars. - # TODO: Experiment with `threshold_desingularization` and `threshold_limiter` - # - # For details on the motivation of the velocity desingularization see - # - A. Chertock, S. Cui, A. Kurganov, T. Wu (2015) - # Well-balanced positivity preserving central-upwind scheme for - # the shallow water system with friction terms - # [DOI: 10.1002/fld.4023](https://doi.org/10.1002/fld.4023) - - # Mortars with the copied or projected solution - tol = equations.threshold_desingularization - for child in 1:2, side in 1:2 - if cache.mortars.u[side, 1, child, i, mortar] <= - equations.threshold_limiter - cache.mortars.u[side, 1, child, i, mortar] = equations.threshold_limiter - cache.mortars.u[side, 2, child, i, mortar] = zero(eltype(u)) - cache.mortars.u[side, 3, child, i, mortar] = zero(eltype(u)) - - else - h = cache.mortars.u[side, 1, child, i, mortar] - cache.mortars.u[side, 2, child, i, mortar] = h * (2 * h * - cache.mortars.u[side, - 2, - child, - i, - mortar]) / - (h^2 + max(h^2, tol)) - cache.mortars.u[side, 3, child, i, mortar] = h * (2 * h * - cache.mortars.u[side, - 3, - child, - i, - mortar]) / - (h^2 + max(h^2, tol)) - end - end + mortar], 0) + # equations.threshold_limiter) + + # # Safety application of velocity desingularization and water height cutoff on the mortars. + # # TODO: Experiment with `threshold_desingularization` and `threshold_limiter` + # # + # # For details on the motivation of the velocity desingularization see + # # - A. Chertock, S. Cui, A. Kurganov, T. Wu (2015) + # # Well-balanced positivity preserving central-upwind scheme for + # # the shallow water system with friction terms + # # [DOI: 10.1002/fld.4023](https://doi.org/10.1002/fld.4023) + + # # Mortars with the copied or projected solution + # tol = equations.threshold_desingularization + # for child in 1:2, side in 1:2 + # if cache.mortars.u[side, 1, child, i, mortar] <= + # equations.threshold_limiter + # cache.mortars.u[side, 1, child, i, mortar] = equations.threshold_limiter + # cache.mortars.u[side, 2, child, i, mortar] = zero(eltype(u)) + # cache.mortars.u[side, 3, child, i, mortar] = zero(eltype(u)) + + # else + # h = cache.mortars.u[side, 1, child, i, mortar] + # cache.mortars.u[side, 2, child, i, mortar] = h * (2 * h * + # cache.mortars.u[side, + # 2, + # child, + # i, + # mortar]) / + # (h^2 + max(h^2, tol)) + # cache.mortars.u[side, 3, child, i, mortar] = h * (2 * h * + # cache.mortars.u[side, + # 3, + # child, + # i, + # mortar]) / + # (h^2 + max(h^2, tol)) + # end + # end end end From 40208bb4a580ea4739122a6a9bedd0474db8fb75 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Fri, 16 May 2025 07:02:12 +0200 Subject: [PATCH 04/37] cleanup new prolong2mortars function --- src/solvers/dgsem_p4est/dg_2d.jl | 69 ++++++++------------------------ 1 file changed, 17 insertions(+), 52 deletions(-) diff --git a/src/solvers/dgsem_p4est/dg_2d.jl b/src/solvers/dgsem_p4est/dg_2d.jl index d7423ae3..0da5c082 100644 --- a/src/solvers/dgsem_p4est/dg_2d.jl +++ b/src/solvers/dgsem_p4est/dg_2d.jl @@ -197,7 +197,7 @@ end # This version is for the `ShallowWaterMultiLayerEquations2D` which "peels" the pressure # contribution into the nonconservative term for easier well-balancing, # see Ersing et al. (https://doi.org/10.1016/j.jcp.2025.113802). -# Thus, this can use the original `P4estMortarContainer`, `calc_mortar_flux!`, +# Thus, this uses the original `P4estMortarContainer`, `calc_mortar_flux!` # and `mortar_fluxes_to_elements!` from Trixi.jl. # # !!! warning "Experimental code" @@ -257,10 +257,18 @@ function Trixi.prolong2mortars!(cache, u, # in the ALE-DG community for the Euler equations with gravity to remove this assumption. # That is, we might be able to directly project the water height `h` instead while maintaining # important steady-state solution behavior. - # Note, a small shift is required to ensure we catch water heights close to the threshold - if u[1, i_large, j_large, element] >= - 2 * (equations.threshold_limiter) - u_buffer[1, i] = u[1, i_large, j_large, element] + + # + # Further, worth noting that for fully wet portions of the domain, the logic below + # is unnecessary and the "classic" mortar method is well-balanced and fully conservative + # for the `ShallowWaterMultiLayerEquations2D`. Could be exploited to design a better + # approach to the AMR refinement / coarsening. + # + # Note, some shift is required to ensure we catch water heights close to the threshold + # and maintain well-balancedness. The larger this shift we also maintain the conservation + # errors to be close to double precision roundoff for longer. + # if u[1, i_large, j_large, element] >= 500 * equations.threshold_limiter # ≈ 5.5e-13 by default + if u[1, i_large, j_large, element] >= equations.threshold_desingularization # 1e-10 by default + u_buffer[1, i] = u[1, i_large, j_large, element] + u[4, i_large, j_large, element] else u_buffer[1, i] = equations.H0 @@ -286,54 +294,11 @@ function Trixi.prolong2mortars!(cache, u, # and instead recover the conservative water height variable `h`. # Basically, unpacking the total water height variable to create the projected local water # height `h` from Eq. 41 in Benov et al. + # Note, no desingularization or water height cutoff is needed here as these steps are done + # later in the hydrostatic reconstruction and stage limiter, respectively. for i in eachnode(dg) - cache.mortars.u[2, 1, 1, i, mortar] = max(cache.mortars.u[2, 1, 1, i, - mortar] - - cache.mortars.u[2, 4, 1, i, - mortar], 0) - # equations.threshold_limiter) - cache.mortars.u[2, 1, 2, i, mortar] = max(cache.mortars.u[2, 1, 2, i, - mortar] - - cache.mortars.u[2, 4, 2, i, - mortar], 0) - # equations.threshold_limiter) - - # # Safety application of velocity desingularization and water height cutoff on the mortars. - # # TODO: Experiment with `threshold_desingularization` and `threshold_limiter` - # # - # # For details on the motivation of the velocity desingularization see - # # - A. Chertock, S. Cui, A. Kurganov, T. Wu (2015) - # # Well-balanced positivity preserving central-upwind scheme for - # # the shallow water system with friction terms - # # [DOI: 10.1002/fld.4023](https://doi.org/10.1002/fld.4023) - - # # Mortars with the copied or projected solution - # tol = equations.threshold_desingularization - # for child in 1:2, side in 1:2 - # if cache.mortars.u[side, 1, child, i, mortar] <= - # equations.threshold_limiter - # cache.mortars.u[side, 1, child, i, mortar] = equations.threshold_limiter - # cache.mortars.u[side, 2, child, i, mortar] = zero(eltype(u)) - # cache.mortars.u[side, 3, child, i, mortar] = zero(eltype(u)) - - # else - # h = cache.mortars.u[side, 1, child, i, mortar] - # cache.mortars.u[side, 2, child, i, mortar] = h * (2 * h * - # cache.mortars.u[side, - # 2, - # child, - # i, - # mortar]) / - # (h^2 + max(h^2, tol)) - # cache.mortars.u[side, 3, child, i, mortar] = h * (2 * h * - # cache.mortars.u[side, - # 3, - # child, - # i, - # mortar]) / - # (h^2 + max(h^2, tol)) - # end - # end + cache.mortars.u[2, 1, 1, i, mortar] = cache.mortars.u[2, 1, 1, i, mortar] - cache.mortars.u[2, 4, 1, i, mortar] + cache.mortars.u[2, 1, 2, i, mortar] = cache.mortars.u[2, 1, 2, i, mortar] - cache.mortars.u[2, 4, 2, i, mortar] end end From 3101afdb371c1339be69585880b7f51270defcc5 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Fri, 16 May 2025 07:35:26 +0200 Subject: [PATCH 05/37] cleanup new well-balancedness test elixir --- ...r_shallowwater_multilayer_well_balanced.jl | 33 +++++++------------ 1 file changed, 12 insertions(+), 21 deletions(-) diff --git a/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_well_balanced.jl b/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_well_balanced.jl index 4db9318e..1b4859f4 100644 --- a/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_well_balanced.jl +++ b/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_well_balanced.jl @@ -4,18 +4,11 @@ using Trixi using TrixiShallowWater ############################################################################### -# Semidiscretization of the shallow water equations with a discontinuous -# bottom topography function and wet/dry transition elements on a nonconforming mesh +# Semidiscretization of the multilayer shallow water equations with one layer +# to fallback to the standard shallow water equations. +# Test case for well-balancedness with wet/dry transition elements on a nonconforming mesh. -############################################################################### -# Semidiscretization of the multilayer shallow water equations with one layers - -equations = ShallowWaterMultiLayerEquations2D(gravity = 9.812, - # H0 = 2.0, # for fully wet testing - H0 = 1.235, # for wet/dry transition testing - # To have the same thresholds as the `ShallowWater2D` - threshold_desingularization = 1e-4, - threshold_limiter = 500*eps(), +equations = ShallowWaterMultiLayerEquations2D(gravity = 9.812, H0 = 1.235, rhos = (1.0)) # An initial condition with constant total water height and zero velocities to test well-balancedness. @@ -70,19 +63,17 @@ surface_flux = (FluxHydrostaticReconstruction(FluxPlusDissipation(flux_ersing_et # Create the solver basis = LobattoLegendreBasis(3) -# TODO: Write a generic function so that I can use `waterheight` in the limiter. -# That way this part of the elixir is identical to the `ShallowWater2D` counterpart +# Cannot simply use `waterheight` here for multilayer equations. +# Need a helper function to extract the relevant variable. +@inline function main_waterheight(u, equations) + return waterheight(u, equations)[1] +end indicator_sc = IndicatorHennemannGassnerShallowWater(equations, basis, alpha_max = 0.5, alpha_min = 0.001, alpha_smooth = true, - # OBS! One cannot simply use `waterheight` - # here for multilayer equations. - # Throws a conversion error. - # Would need another function that returns sum(h) - # after the `waterheight` call. - variable = waterheight_pressure) + variable = main_waterheight) volume_integral = VolumeIntegralShockCapturingHG(indicator_sc; volume_flux_dg = volume_flux, volume_flux_fv = surface_flux) @@ -209,11 +200,11 @@ analysis_callback = AnalysisCallback(semi, interval = analysis_interval, alive_callback = AliveCallback(analysis_interval = analysis_interval) -save_solution = SaveSolutionCallback(dt = 0.1, # 10.0, +save_solution = SaveSolutionCallback(dt = 10.0, save_initial_solution = true, save_final_solution = true) -stepsize_callback = StepsizeCallback(cfl = 0.5) +stepsize_callback = StepsizeCallback(cfl = 1.0) callbacks = CallbackSet(summary_callback, analysis_callback, From 9b01c6253fe0cae47c754adc4359e9525c80d9d0 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Fri, 16 May 2025 07:36:41 +0200 Subject: [PATCH 06/37] better elixir file name --- ...hallowwater_multilayer_well_balanced_wet_dry_nonconforming.jl} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename examples/p4est_2d_dgsem/{elixir_shallowwater_multilayer_well_balanced.jl => elixir_shallowwater_multilayer_well_balanced_wet_dry_nonconforming.jl} (100%) diff --git a/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_well_balanced.jl b/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_well_balanced_wet_dry_nonconforming.jl similarity index 100% rename from examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_well_balanced.jl rename to examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_well_balanced_wet_dry_nonconforming.jl From 0b227679efd294f3cb6f2d9e4262cb29fd256e9f Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Fri, 16 May 2025 08:12:11 +0200 Subject: [PATCH 07/37] cleanup and add AMR pertubation test file --- ...ter_multilayer_perturbation_wet_dry_amr.jl | 238 ++++++++++++++++++ ...yer_well_balanced_wet_dry_nonconforming.jl | 6 +- 2 files changed, 241 insertions(+), 3 deletions(-) create mode 100644 examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_perturbation_wet_dry_amr.jl diff --git a/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_perturbation_wet_dry_amr.jl b/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_perturbation_wet_dry_amr.jl new file mode 100644 index 00000000..da6067ee --- /dev/null +++ b/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_perturbation_wet_dry_amr.jl @@ -0,0 +1,238 @@ + +using OrdinaryDiffEqSSPRK, OrdinaryDiffEqLowStorageRK +using Trixi +using TrixiShallowWater + +# Elixir useful for stress testing the AMR + wet/dry + well-balanced mortars. + +############################################################################### +# Semidiscretization of the multilayer shallow water equations with one layer +# to actually be the standard shallow water equations with a discontinuous +# bottom topography function for a perturbed water height +# on a nonconforming mesh with AMR + +equations = ShallowWaterMultiLayerEquations2D(gravity = 9.812, H0 = 1.235, + rhos = (1.0)) + +function initial_condition_perturbation(x, t, equations::ShallowWaterMultiLayerEquations2D) + # Calculate primitive variables + H = [equations.H0] + v1 = zero(H) + v2 = zero(H) + + x1, x2 = x + b = (2.5 / exp(0.5 * ((x1 - 1.0)^2 + (x2 + 1.0)^2)) + + + 0.8 / exp(0.5 * ((x1 + 1.0)^2 + (x2 - 1.0)^2)) + - + 0.5 / exp(3.5 * ((x1 - 0.4)^2 + (x2 - 0.325)^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_perturbation + +# Wall BCs +boundary_condition = Dict(:all => 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)) + +# Create the solver +basis = LobattoLegendreBasis(4) + +# Cannot simply use `waterheight` here for multilayer equations. +# Need a helper function to extract the relevant variable. +@inline function main_waterheight(u, equations) + return waterheight(u, equations)[1] +end + +indicator_sc = IndicatorHennemannGassnerShallowWater(equations, basis, + alpha_max = 0.5, + alpha_min = 0.001, + alpha_smooth = true, + variable = main_waterheight) + +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) + +# Unstructured mesh with 24 cells of the square domain [-1, 1]^2 +mesh_file = Trixi.download("https://gist.githubusercontent.com/efaulhaber/63ff2ea224409e55ee8423b3a33e316a/raw/7db58af7446d1479753ae718930741c47a3b79b7/square_unstructured_2.inp", + joinpath(@__DIR__, "square_unstructured_2.inp")) + +# Affine type mapping to take the [-1,1]^2 domain from the mesh file +# and warp it as described in https://arxiv.org/abs/2012.12040 +# Warping with the coefficient 0.15 is even more extreme. +function mapping_twist(xi, eta) + y = eta + 0.15 * cos(1.5 * pi * xi) * cos(0.5 * pi * eta) + x = xi + 0.15 * cos(0.5 * pi * xi) * cos(2.0 * pi * y) + return SVector(x, y) +end + +mesh = P4estMesh{2}(mesh_file, polydeg = 4, + mapping = mapping_twist, + initial_refinement_level = 0) + +# Refine bottom left quadrant of each tree to level 3 +function refine_fn(p4est, which_tree, quadrant) + quadrant_obj = unsafe_load(quadrant) + if quadrant_obj.x == 0 && quadrant_obj.y == 0 && quadrant_obj.level < 3 + # return true (refine) + return Cint(1) + else + # return false (don't refine) + return Cint(0) + end +end + +# Refine recursively until each bottom left quadrant of a tree has level 3 +# The mesh will be rebalanced before the simulation starts +refine_fn_c = @cfunction(refine_fn, Cint, + (Ptr{Trixi.p4est_t}, Ptr{Trixi.p4est_topidx_t}, + Ptr{Trixi.p4est_quadrant_t})) +Trixi.refine_p4est!(mesh.p4est, true, refine_fn_c, C_NULL) + +# Create the semi discretization object +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_condition) + +############################################################################### +# ODE solvers, callbacks, etc. + +tspan = (0.0, 4.0) +ode = semidiscretize(semi, tspan) + +function initial_condition_discontinuous_perturbation(x, t, element_id, + equations::ShallowWaterMultiLayerEquations2D) + # Set the background values for velocity + H = [equations.H0] + v1 = zero(H) + v2 = zero(H) + + x1, x2 = x + + # For the mesh file version of the testing + b = (1.75 / exp(0.6 * ((x1 - 1.0)^2 + (x2 + 1.0)^2)) + + + 0.8 / exp(0.5 * ((x1 + 1.0)^2 + (x2 - 1.0)^2)) + - + 0.5 / exp(3.5 * ((x1 - 0.4)^2 + (x2 - 0.325)^2))) + + # Setup a discontinuous bottom topography using the element id number + # Note that this requires to have the mesh exactly as created above, + # any additional refinement changes the initial condition, + IDs = [collect(114:133); collect(138:141); collect(156:164); collect(208:300)] + if element_id in IDs + b = (0.75 / exp(0.5 * ((x1 - 1.0)^2 + (x2 + 1.0)^2)) + + + 0.4 / exp(0.5 * ((x1 + 1.0)^2 + (x2 - 1.0)^2)) + - + 0.25 / exp(3.5 * ((x1 - 0.4)^2 + (x2 - 0.325)^2))) + end + + # Put in a discontinuous perturbation using the element number + if element_id in [232, 224, 225, 226, 227, 228, 229, 230] + H = H .+ 1.6 + 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 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_perturbation(x_node, first(tspan), element, + equations) + Trixi.set_node_vars!(u, u_node, equations, semi.solver, i, j, element) + end +end + +############################################################################### + +summary_callback = SummaryCallback() + +analysis_interval = 1000 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + extra_analysis_errors = (:conservation_error,), + extra_analysis_integrals = (lake_at_rest_error,)) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(dt = 0.2, + save_initial_solution = true, + save_final_solution = true) + +amr_controller = ControllerThreeLevel(semi, + IndicatorMax(semi, variable = main_waterheight), + base_level = 0, + med_level = 1, med_threshold = 0.3, + max_level = 4, max_threshold = 1.15) + +amr_callback = AMRCallback(semi, amr_controller, + interval = 5, + adapt_initial_condition = false, + adapt_initial_condition_only_refine = false) + +stepsize_callback = StepsizeCallback(cfl = 0.5) + +callbacks = CallbackSet(summary_callback, + analysis_callback, + alive_callback, + save_solution, + amr_callback, + stepsize_callback) + +stage_limiter! = PositivityPreservingLimiterShallowWater(variables = (waterheight,)) + +############################################################################### +# run the simulation +sol = solve(ode, SSPRK43(stage_limiter!), + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + adaptive = false, + save_everystep = false, callback = callbacks); diff --git a/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_well_balanced_wet_dry_nonconforming.jl b/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_well_balanced_wet_dry_nonconforming.jl index 1b4859f4..455db40d 100644 --- a/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_well_balanced_wet_dry_nonconforming.jl +++ b/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_well_balanced_wet_dry_nonconforming.jl @@ -126,7 +126,7 @@ semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, ############################################################################### # ODE solvers, callbacks, etc. -tspan = (0.0, 50.0) +tspan = (0.0, 100.0) ode = semidiscretize(semi, tspan) ############################################################################### @@ -209,7 +209,7 @@ stepsize_callback = StepsizeCallback(cfl = 1.0) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, - save_solution, + # save_solution, stepsize_callback) stage_limiter! = PositivityPreservingLimiterShallowWater(variables = (waterheight,)) @@ -219,4 +219,4 @@ stage_limiter! = PositivityPreservingLimiterShallowWater(variables = (waterheigh sol = solve(ode, SSPRK43(stage_limiter!), dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback adaptive = false, - save_everystep = false, callback = callbacks, maxiters=1e9); + save_everystep = false, callback = callbacks); From 84bd0d9d7fa1d535bd37e95a4b325ab22968a404 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Fri, 16 May 2025 08:16:45 +0200 Subject: [PATCH 08/37] add Monai flood with AMR elixir --- ...shallowwater_multilayer_monai_flood_amr.jl | 259 ++++++++++++++++++ ...ter_multilayer_perturbation_wet_dry_amr.jl | 2 +- 2 files changed, 260 insertions(+), 1 deletion(-) create mode 100644 examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_monai_flood_amr.jl diff --git a/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_monai_flood_amr.jl b/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_monai_flood_amr.jl new file mode 100644 index 00000000..315db98f --- /dev/null +++ b/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_monai_flood_amr.jl @@ -0,0 +1,259 @@ + +using OrdinaryDiffEqSSPRK, OrdinaryDiffEqLowStorageRK +using Trixi +using TrixiShallowWater +using TrixiBottomTopography + +# See the tutorial +# https://trixi-framework.github.io/TrixiShallowWater.jl/stable/tutorials/elixir_shallowwater_monai_tsunami/ +# for a thorough explanation of this problem setup. +# +# This elixir allows for experimentation with AMR with this complex test case. + +############################################################################### +# Semidiscretization of the multilayer shallow water equations with one layer +# to fallback to be the standard shallow water equations + +equations = ShallowWaterMultiLayerEquations2D(gravity = 9.81, H0 = 0.0, + rhos = (1.0)) + +# Setup the bottom topography data and cubic spline interpolant + +# Create a bicubic B-spline interpolation of the bathymetry data, then create a function +# to evaluate the resulting spline at a given point $(x,y)$. +spline_bathymetry_file = Trixi.download("https://gist.githubusercontent.com/andrewwinters5000/21255c980c4eda5294f91e8dfe6c7e33/raw/1afb73928892774dc3a902e0c46ffd882ef03ee3/monai_bathymetry_data.txt", + joinpath(@__DIR__, "monai_bathymetry_data.txt")); + +# B-spline interpolation of the underlying data +bath_spline_struct = BicubicBSpline(spline_bathymetry_file, end_condition="not-a-knot") +bathymetry(x, y) = spline_interpolation(bath_spline_struct, x, y) + +# Initial condition is a constant background state +function initial_condition_monai_tsunami(x, t, equations::ShallowWaterMultiLayerEquations2D) + # Set the background water height + H = [equations.H0] + + # Initially water is at rest + v1 = zero(H) + v2 = zero(H) + + # Bottom topography values are computed from the bicubic spline created above + x1, x2 = x + b = bathymetry(x1, x2) + + # 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 the conservative variables + return prim2cons(SVector(H..., v1..., v2..., b), equations) +end + +initial_condition = initial_condition_monai_tsunami + +# Tsunami test case uses a specialized wave maker type of boundary condition. +# It is used to model an incident wave that approaches from off-shore +# with a water depth of $h = 13.535\,\text{cm}$. To create the incident wave information +# that is valid over the time interval $t \in [0\,s, 22.5\,s]$. + +# We download the incident wave data that has been preprocessed to be in the format +# required by TrixiBottomTopography. +water_height_data = Trixi.download("https://gist.githubusercontent.com/andrewwinters5000/5b11f5f175bddb326d11d8e28398127e/raw/64980e0e4526e0fcd49589b34ee5458b9a1cebff/monai_wavemaker_bc.txt", + joinpath(@__DIR__, "monai_wavemaker_bc.txt")); + +# Similar to the bathymetry approximation, we construct a cubic B-spline interpolation +# of the data, then create a function to evaluate the resulting spline at a given $t$ value. +h_spline_struct = CubicBSpline(water_height_data; end_condition = "not-a-knot") +H_from_wave_maker(t) = spline_interpolation(h_spline_struct, t) + +# Now we can define the specialized boundary condition for the incident wave maker. +function boundary_condition_wave_maker(u_inner, normal_direction::AbstractVector, x, t, + surface_flux_functions, equations::ShallowWaterMultiLayerEquations2D) + + surface_flux_function, nonconservative_flux_function = surface_flux_functions + + # Compute the water height from the wave maker input file data + # and then clip to avoid negative water heights and division by zero + h_ext = max(equations.threshold_limiter, H_from_wave_maker(t) - u_inner[4]) + + # Compute the incoming velocity as in Eq. (10) of the paper + # - S. Vater, N. Beisiegel, and J. Behrens (2019) + # A limiter-based well-balanced discontinuous Galerkin method for shallow-water flows + # with wetting and drying: Triangular grids + # [DOI: 10.1002/fld.4762](https://doi.org/10.1002/fld.4762) + h0 = 0.13535 # reference incident water height converted to meters + v1_ext = 2 * (sqrt(equations.gravity * h_ext) - sqrt(equations.gravity * h0)) + + # Create the external solution state in the conservative variables + u_outer = SVector([h_ext]..., [h_ext * v1_ext]..., [zero(eltype(x))]..., u_inner[4]) + + # Calculate the boundary flux and nonconservative contributions + flux = surface_flux_function(u_inner, u_outer, normal_direction, equations) + + noncons_flux = nonconservative_flux_function(u_inner, u_outer, normal_direction, + equations) + + return flux, noncons_flux +end + +boundary_condition = Dict(:Bottom => boundary_condition_slip_wall, + :Top => boundary_condition_slip_wall, + :Right => boundary_condition_slip_wall, + :Left => boundary_condition_wave_maker) + +# Manning friction source term +@inline function source_terms_manning_friction(u, x, t, equations::ShallowWaterMultiLayerEquations2D) + # Grab the conservative variables + h, hv_1, hv_2, b = u + + # Desingularization + sh = (2.0 * h[1]) / (h[1]^2 + max(h[1]^2, 1e-8)) + + # friction coefficient + n = 0.001 + + # Compute the common friction term + Sf = -equations.gravity * n^2 * sh^(7.0 / 3.0) * sqrt(hv_1[1]^2 + hv_2[1]^2) + + return SVector(zero(h)..., + Sf * hv_1[1]..., + Sf * hv_2[1]..., + zero(b)) +end + +############################################################################### +# 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) + +# Cannot simply use `waterheight` here for multilayer equations. +# Need a helper function to extract the relevant variable. +@inline function main_waterheight(u, equations) + return waterheight(u, equations)[1] +end + +indicator_sc = IndicatorHennemannGassnerShallowWater(equations, basis, + alpha_max=0.5, + alpha_min=0.001, + alpha_smooth=true, + variable=main_waterheight) +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/d26356bdc8e8d742a3035b3f71c71a68/raw/9d6ceedb844e92313d1dac2318a28c87ffbb9de2/mesh_monai_shore.inp", + joinpath(@__DIR__, "mesh_monai_shore.inp")); + +mesh = P4estMesh{2}(mesh_file) + +# create the semi discretization object +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver; + boundary_conditions=boundary_condition, + source_terms=source_terms_manning_friction) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 22.5) +ode = semidiscretize(semi, tspan) + +############################################################################### +# Callbacks + +summary_callback = SummaryCallback() + +analysis_interval = 1000 +analysis_callback = AnalysisCallback(semi, interval=analysis_interval, + extra_analysis_errors = (:conservation_error,), + save_analysis = false) + +alive_callback = AliveCallback(analysis_interval=analysis_interval) + +save_solution = SaveSolutionCallback(dt = 0.2, + save_initial_solution=true, + save_final_solution=true) + +# TimeSeries Callback is currently unavailable on `P4estMesh` +# time_series = TimeSeriesCallback(semi, [(4.521, 1.196), # G5 +# (4.521, 1.696), # G7 +# (4.521, 2.196)]; # G9 +# interval = 2, +# solution_variables=cons2prim) + +stepsize_callback = StepsizeCallback(cfl = 0.5) + +# Another possible AMR indicator function could be the velocity, such that it only fires +# in regions where the water is moving +@inline function momentum_norm(u, equations::ShallowWaterMultiLayerEquations2D) + _, h_v1, h_v2, _ = u + return sqrt(h_v1[1]^2 + h_v2[1]^2) + # h, h_v1, h_v2, _ = u + # v1 = (2 * h[1] * h_v1[1]) / (h[1]^2 + max(h[1]^2, 1e-8)) + # v2 = (2 * h[1] * h_v2[1]) / (h[1]^2 + max(h[1]^2, 1e-8)) + # if h[1] <= 1e-8 + # v1 = 0.0 + # v2 = 0.0 + # end + # if v1 > 0.525 || v2 > 0.525 + # v1 = 0.0 + # v2 = 0.0 + # end + # return sqrt(v1^2 + v2^2) +end + +amr_controller = ControllerThreeLevel(semi, + IndicatorMax(semi, variable = momentum_norm), + base_level = 0, + med_level = 1, med_threshold = 0.012, + max_level = 3, max_threshold = 0.015) + # velocity is very sensitive to use as an indicator + # med_level = 1, med_threshold = 0.25, + # max_level = 2, max_threshold = 0.475) + +amr_callback = AMRCallback(semi, amr_controller, + interval = 5, + adapt_initial_condition = false, + adapt_initial_condition_only_refine = false) + +callbacks = CallbackSet(summary_callback, + analysis_callback, + alive_callback, + # time_series, + amr_callback, + save_solution, + stepsize_callback) + +############################################################################### +# run the simulation + +stage_limiter! = PositivityPreservingLimiterShallowWater(variables = (waterheight,)) + +sol = solve(ode, SSPRK43(stage_limiter!); + ode_default_options()..., + callback = callbacks, + adaptive = false, + dt = 1.0 # solve needs some value here; overwritten by the `stepsize_callback` + ); diff --git a/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_perturbation_wet_dry_amr.jl b/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_perturbation_wet_dry_amr.jl index da6067ee..60df89e8 100644 --- a/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_perturbation_wet_dry_amr.jl +++ b/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_perturbation_wet_dry_amr.jl @@ -127,7 +127,7 @@ semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, ############################################################################### # ODE solvers, callbacks, etc. -tspan = (0.0, 4.0) +tspan = (0.0, 5.0) ode = semidiscretize(semi, tspan) function initial_condition_discontinuous_perturbation(x, t, element_id, From 3170e460d172ed60841b94d83a5ca7c9ed1d6955 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Fri, 16 May 2025 09:53:05 +0200 Subject: [PATCH 09/37] add tests --- test/test_p4est_2d.jl | 84 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 84 insertions(+) diff --git a/test/test_p4est_2d.jl b/test/test_p4est_2d.jl index fb2e9959..f484ba1e 100644 --- a/test/test_p4est_2d.jl +++ b/test/test_p4est_2d.jl @@ -126,6 +126,90 @@ isdir(outdir) && rm(outdir, recursive = true) end end end # SWE + +@testset "Multilayer Shallow Water" begin + @trixi_testset "elixir_shallowwater_multilayer_well_balanced_wet_dry_nonconforming.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_multilayer_well_balanced_wet_dry_nonconforming.jl"), + l2=[ + 0.17389058166483545, + 1.3442539891317369e-15, + 1.390606878462557e-15, + 0.18415049792609314], + linf=[ + 0.4160052818017864, + 1.0894983713455818e-14, + 1.132829841145215e-14, + 0.41600528180178625], + 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 + + # Note, these values may change as the functionality of well-balanced mortars + # with AMR and wet/dry are further developed according to the issue + # https://github.com/trixi-framework/TrixiShallowWater.jl/issues/77 + @trixi_testset "elixir_shallowwater_multilayer_perturbation_wet_dry_amr.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_multilayer_perturbation_wet_dry_amr.jl"), + l2=[ + 0.3997030663722126, + 0.38644026433480094, + 0.4362923911447016, + 0.4564350570616195], + linf=[ + 1.372326981328246, + 3.3262654799249893, + 3.6541077514853653, + 0.7495177590247986], + tspan=(0.0, 0.025), + coverage_override=(maxiters = 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 + + # Note, complex example that uses TrixiBottomTopography.jl to approximate + # bathymetry and wave maker boundary condition. Tests several components + # of the TrixiShallowWater.jl toolchain. Note, does not run long enough + # for the AMR to fire. + @trixi_testset "elixir_shallowwater_multilayer_perturbation_wet_dry_amr.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_multilayer_perturbation_wet_dry_amr.jl"), + l2=[ + 0.0003068569689525136, + 3.869268640941544e-6, + 6.345700409103784e-11, + 0.0003087516514830762], + linf=[ + 0.004387570753720829, + 3.562731640300165e-5, + 1.190336978580838e-9, + 0.004387569562413193], + 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 # P4estMesh2D # Clean up afterwards: delete TrixiShallowWater.jl output directory From 9a724d5018cbc3ba0b13db2f9b3cab74db5f3e59 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Fri, 16 May 2025 10:29:14 +0200 Subject: [PATCH 10/37] fix incorrect file name in new test --- test/test_p4est_2d.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_p4est_2d.jl b/test/test_p4est_2d.jl index f484ba1e..6272af76 100644 --- a/test/test_p4est_2d.jl +++ b/test/test_p4est_2d.jl @@ -185,9 +185,9 @@ end # SWE # bathymetry and wave maker boundary condition. Tests several components # of the TrixiShallowWater.jl toolchain. Note, does not run long enough # for the AMR to fire. - @trixi_testset "elixir_shallowwater_multilayer_perturbation_wet_dry_amr.jl" begin + @trixi_testset "elixir_shallowwater_multilayer_monai_flood_amr.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, - "elixir_shallowwater_multilayer_perturbation_wet_dry_amr.jl"), + "elixir_shallowwater_multilayer_monai_flood_amr.jl"), l2=[ 0.0003068569689525136, 3.869268640941544e-6, From 74e89a27d180c89e979b20c61775a90fc4cc1c06 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Fri, 16 May 2025 18:41:07 +0200 Subject: [PATCH 11/37] add TrixiBottomTopography to testing Project.toml --- test/Project.toml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/test/Project.toml b/test/Project.toml index 805e9a2d..9e167ec0 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -5,6 +5,7 @@ Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" Trixi = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb" +TrixiBottomTopography = "86af9953-43df-404b-8eaa-d20d82623a82" TrixiTest = "0a316866-cbd0-4425-8bcb-08103b2c1f26" [compat] @@ -14,4 +15,5 @@ Printf = "1" Roots = "2.1.6" Test = "1" Trixi = "0.11.10" +TrixiBottomTopography = "0.1" TrixiTest = "0.1" From 297ad01e63e543ab8e8912b9ec59401d35fee32d Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Fri, 16 May 2025 18:41:15 +0200 Subject: [PATCH 12/37] apply formatter --- src/callbacks_step/amr_dg2d.jl | 6 ++++-- src/solvers/dgsem_p4est/dg_2d.jl | 10 ++++++---- 2 files changed, 10 insertions(+), 6 deletions(-) diff --git a/src/callbacks_step/amr_dg2d.jl b/src/callbacks_step/amr_dg2d.jl index c10e0ede..473a08c1 100644 --- a/src/callbacks_step/amr_dg2d.jl +++ b/src/callbacks_step/amr_dg2d.jl @@ -19,7 +19,8 @@ # !!! warning "Experimental code" # This is an experimental feature and may change in future releases. function Trixi.refine!(u_ode::AbstractVector, adaptor, mesh::P4estMesh{2}, - equations::Union{ShallowWaterEquationsWetDry2D,ShallowWaterMultiLayerEquations2D}, + equations::Union{ShallowWaterEquationsWetDry2D, + ShallowWaterMultiLayerEquations2D}, dg::DGSEM, cache, elements_to_refine) # Return early if there is nothing to do if isempty(elements_to_refine) @@ -123,7 +124,8 @@ end # !!! warning "Experimental code" # This is an experimental feature and may change in future releases. function Trixi.coarsen!(u_ode::AbstractVector, adaptor, mesh::P4estMesh{2}, - equations::Union{ShallowWaterEquationsWetDry2D,ShallowWaterMultiLayerEquations2D}, + equations::Union{ShallowWaterEquationsWetDry2D, + ShallowWaterMultiLayerEquations2D}, dg::DGSEM, cache, elements_to_remove) # Return early if there is nothing to do if isempty(elements_to_remove) diff --git a/src/solvers/dgsem_p4est/dg_2d.jl b/src/solvers/dgsem_p4est/dg_2d.jl index 0da5c082..6d9f56f0 100644 --- a/src/solvers/dgsem_p4est/dg_2d.jl +++ b/src/solvers/dgsem_p4est/dg_2d.jl @@ -268,8 +268,8 @@ function Trixi.prolong2mortars!(cache, u, # errors to be close to double precision roundoff for longer. # if u[1, i_large, j_large, element] >= 500 * equations.threshold_limiter # ≈ 5.5e-13 by default if u[1, i_large, j_large, element] >= equations.threshold_desingularization # 1e-10 by default - u_buffer[1, i] = u[1, i_large, j_large, element] + - u[4, i_large, j_large, element] + u_buffer[1, i] = (u[1, i_large, j_large, element] + + u[4, i_large, j_large, element]) else u_buffer[1, i] = equations.H0 end @@ -297,8 +297,10 @@ function Trixi.prolong2mortars!(cache, u, # Note, no desingularization or water height cutoff is needed here as these steps are done # later in the hydrostatic reconstruction and stage limiter, respectively. for i in eachnode(dg) - cache.mortars.u[2, 1, 1, i, mortar] = cache.mortars.u[2, 1, 1, i, mortar] - cache.mortars.u[2, 4, 1, i, mortar] - cache.mortars.u[2, 1, 2, i, mortar] = cache.mortars.u[2, 1, 2, i, mortar] - cache.mortars.u[2, 4, 2, i, mortar] + cache.mortars.u[2, 1, 1, i, mortar] = (cache.mortars.u[2, 1, 1, i, mortar] - + cache.mortars.u[2, 4, 1, i, mortar]) + cache.mortars.u[2, 1, 2, i, mortar] = (cache.mortars.u[2, 1, 2, i, mortar] - + cache.mortars.u[2, 4, 2, i, mortar]) end end From d410365b6da586e1a396eebcb9402a473fb1f6b1 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Fri, 16 May 2025 23:27:30 +0200 Subject: [PATCH 13/37] add three mound with AMR elixir --- ...shallowwater_multilayer_monai_flood_amr.jl | 518 +++++++++--------- ...wwater_multilayer_three_mound_dam_break.jl | 182 ++++++ ...yer_well_balanced_wet_dry_nonconforming.jl | 3 +- 3 files changed, 443 insertions(+), 260 deletions(-) create mode 100644 examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_three_mound_dam_break.jl diff --git a/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_monai_flood_amr.jl b/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_monai_flood_amr.jl index 315db98f..c8e599f3 100644 --- a/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_monai_flood_amr.jl +++ b/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_monai_flood_amr.jl @@ -1,259 +1,259 @@ - -using OrdinaryDiffEqSSPRK, OrdinaryDiffEqLowStorageRK -using Trixi -using TrixiShallowWater -using TrixiBottomTopography - -# See the tutorial -# https://trixi-framework.github.io/TrixiShallowWater.jl/stable/tutorials/elixir_shallowwater_monai_tsunami/ -# for a thorough explanation of this problem setup. -# -# This elixir allows for experimentation with AMR with this complex test case. - -############################################################################### -# Semidiscretization of the multilayer shallow water equations with one layer -# to fallback to be the standard shallow water equations - -equations = ShallowWaterMultiLayerEquations2D(gravity = 9.81, H0 = 0.0, - rhos = (1.0)) - -# Setup the bottom topography data and cubic spline interpolant - -# Create a bicubic B-spline interpolation of the bathymetry data, then create a function -# to evaluate the resulting spline at a given point $(x,y)$. -spline_bathymetry_file = Trixi.download("https://gist.githubusercontent.com/andrewwinters5000/21255c980c4eda5294f91e8dfe6c7e33/raw/1afb73928892774dc3a902e0c46ffd882ef03ee3/monai_bathymetry_data.txt", - joinpath(@__DIR__, "monai_bathymetry_data.txt")); - -# B-spline interpolation of the underlying data -bath_spline_struct = BicubicBSpline(spline_bathymetry_file, end_condition="not-a-knot") -bathymetry(x, y) = spline_interpolation(bath_spline_struct, x, y) - -# Initial condition is a constant background state -function initial_condition_monai_tsunami(x, t, equations::ShallowWaterMultiLayerEquations2D) - # Set the background water height - H = [equations.H0] - - # Initially water is at rest - v1 = zero(H) - v2 = zero(H) - - # Bottom topography values are computed from the bicubic spline created above - x1, x2 = x - b = bathymetry(x1, x2) - - # 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 the conservative variables - return prim2cons(SVector(H..., v1..., v2..., b), equations) -end - -initial_condition = initial_condition_monai_tsunami - -# Tsunami test case uses a specialized wave maker type of boundary condition. -# It is used to model an incident wave that approaches from off-shore -# with a water depth of $h = 13.535\,\text{cm}$. To create the incident wave information -# that is valid over the time interval $t \in [0\,s, 22.5\,s]$. - -# We download the incident wave data that has been preprocessed to be in the format -# required by TrixiBottomTopography. -water_height_data = Trixi.download("https://gist.githubusercontent.com/andrewwinters5000/5b11f5f175bddb326d11d8e28398127e/raw/64980e0e4526e0fcd49589b34ee5458b9a1cebff/monai_wavemaker_bc.txt", - joinpath(@__DIR__, "monai_wavemaker_bc.txt")); - -# Similar to the bathymetry approximation, we construct a cubic B-spline interpolation -# of the data, then create a function to evaluate the resulting spline at a given $t$ value. -h_spline_struct = CubicBSpline(water_height_data; end_condition = "not-a-knot") -H_from_wave_maker(t) = spline_interpolation(h_spline_struct, t) - -# Now we can define the specialized boundary condition for the incident wave maker. -function boundary_condition_wave_maker(u_inner, normal_direction::AbstractVector, x, t, - surface_flux_functions, equations::ShallowWaterMultiLayerEquations2D) - - surface_flux_function, nonconservative_flux_function = surface_flux_functions - - # Compute the water height from the wave maker input file data - # and then clip to avoid negative water heights and division by zero - h_ext = max(equations.threshold_limiter, H_from_wave_maker(t) - u_inner[4]) - - # Compute the incoming velocity as in Eq. (10) of the paper - # - S. Vater, N. Beisiegel, and J. Behrens (2019) - # A limiter-based well-balanced discontinuous Galerkin method for shallow-water flows - # with wetting and drying: Triangular grids - # [DOI: 10.1002/fld.4762](https://doi.org/10.1002/fld.4762) - h0 = 0.13535 # reference incident water height converted to meters - v1_ext = 2 * (sqrt(equations.gravity * h_ext) - sqrt(equations.gravity * h0)) - - # Create the external solution state in the conservative variables - u_outer = SVector([h_ext]..., [h_ext * v1_ext]..., [zero(eltype(x))]..., u_inner[4]) - - # Calculate the boundary flux and nonconservative contributions - flux = surface_flux_function(u_inner, u_outer, normal_direction, equations) - - noncons_flux = nonconservative_flux_function(u_inner, u_outer, normal_direction, - equations) - - return flux, noncons_flux -end - -boundary_condition = Dict(:Bottom => boundary_condition_slip_wall, - :Top => boundary_condition_slip_wall, - :Right => boundary_condition_slip_wall, - :Left => boundary_condition_wave_maker) - -# Manning friction source term -@inline function source_terms_manning_friction(u, x, t, equations::ShallowWaterMultiLayerEquations2D) - # Grab the conservative variables - h, hv_1, hv_2, b = u - - # Desingularization - sh = (2.0 * h[1]) / (h[1]^2 + max(h[1]^2, 1e-8)) - - # friction coefficient - n = 0.001 - - # Compute the common friction term - Sf = -equations.gravity * n^2 * sh^(7.0 / 3.0) * sqrt(hv_1[1]^2 + hv_2[1]^2) - - return SVector(zero(h)..., - Sf * hv_1[1]..., - Sf * hv_2[1]..., - zero(b)) -end - -############################################################################### -# 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) - -# Cannot simply use `waterheight` here for multilayer equations. -# Need a helper function to extract the relevant variable. -@inline function main_waterheight(u, equations) - return waterheight(u, equations)[1] -end - -indicator_sc = IndicatorHennemannGassnerShallowWater(equations, basis, - alpha_max=0.5, - alpha_min=0.001, - alpha_smooth=true, - variable=main_waterheight) -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/d26356bdc8e8d742a3035b3f71c71a68/raw/9d6ceedb844e92313d1dac2318a28c87ffbb9de2/mesh_monai_shore.inp", - joinpath(@__DIR__, "mesh_monai_shore.inp")); - -mesh = P4estMesh{2}(mesh_file) - -# create the semi discretization object -semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver; - boundary_conditions=boundary_condition, - source_terms=source_terms_manning_friction) - -############################################################################### -# ODE solvers, callbacks etc. - -tspan = (0.0, 22.5) -ode = semidiscretize(semi, tspan) - -############################################################################### -# Callbacks - -summary_callback = SummaryCallback() - -analysis_interval = 1000 -analysis_callback = AnalysisCallback(semi, interval=analysis_interval, - extra_analysis_errors = (:conservation_error,), - save_analysis = false) - -alive_callback = AliveCallback(analysis_interval=analysis_interval) - -save_solution = SaveSolutionCallback(dt = 0.2, - save_initial_solution=true, - save_final_solution=true) - -# TimeSeries Callback is currently unavailable on `P4estMesh` -# time_series = TimeSeriesCallback(semi, [(4.521, 1.196), # G5 -# (4.521, 1.696), # G7 -# (4.521, 2.196)]; # G9 -# interval = 2, -# solution_variables=cons2prim) - -stepsize_callback = StepsizeCallback(cfl = 0.5) - -# Another possible AMR indicator function could be the velocity, such that it only fires -# in regions where the water is moving -@inline function momentum_norm(u, equations::ShallowWaterMultiLayerEquations2D) - _, h_v1, h_v2, _ = u - return sqrt(h_v1[1]^2 + h_v2[1]^2) - # h, h_v1, h_v2, _ = u - # v1 = (2 * h[1] * h_v1[1]) / (h[1]^2 + max(h[1]^2, 1e-8)) - # v2 = (2 * h[1] * h_v2[1]) / (h[1]^2 + max(h[1]^2, 1e-8)) - # if h[1] <= 1e-8 - # v1 = 0.0 - # v2 = 0.0 - # end - # if v1 > 0.525 || v2 > 0.525 - # v1 = 0.0 - # v2 = 0.0 - # end - # return sqrt(v1^2 + v2^2) -end - -amr_controller = ControllerThreeLevel(semi, - IndicatorMax(semi, variable = momentum_norm), - base_level = 0, - med_level = 1, med_threshold = 0.012, - max_level = 3, max_threshold = 0.015) - # velocity is very sensitive to use as an indicator - # med_level = 1, med_threshold = 0.25, - # max_level = 2, max_threshold = 0.475) - -amr_callback = AMRCallback(semi, amr_controller, - interval = 5, - adapt_initial_condition = false, - adapt_initial_condition_only_refine = false) - -callbacks = CallbackSet(summary_callback, - analysis_callback, - alive_callback, - # time_series, - amr_callback, - save_solution, - stepsize_callback) - -############################################################################### -# run the simulation - -stage_limiter! = PositivityPreservingLimiterShallowWater(variables = (waterheight,)) - -sol = solve(ode, SSPRK43(stage_limiter!); - ode_default_options()..., - callback = callbacks, - adaptive = false, - dt = 1.0 # solve needs some value here; overwritten by the `stepsize_callback` - ); + +using OrdinaryDiffEqSSPRK, OrdinaryDiffEqLowStorageRK +using Trixi +using TrixiShallowWater +using TrixiBottomTopography + +# See the tutorial +# https://trixi-framework.github.io/TrixiShallowWater.jl/stable/tutorials/elixir_shallowwater_monai_tsunami/ +# for a thorough explanation of this problem setup. +# +# This elixir allows for experimentation with AMR with this complex test case. + +############################################################################### +# Semidiscretization of the multilayer shallow water equations with one layer +# to fallback to be the standard shallow water equations + +equations = ShallowWaterMultiLayerEquations2D(gravity = 9.81, H0 = 0.0, + rhos = (1.0)) + +# Setup the bottom topography data and cubic spline interpolant + +# Create a bicubic B-spline interpolation of the bathymetry data, then create a function +# to evaluate the resulting spline at a given point $(x,y)$. +spline_bathymetry_file = Trixi.download("https://gist.githubusercontent.com/andrewwinters5000/21255c980c4eda5294f91e8dfe6c7e33/raw/1afb73928892774dc3a902e0c46ffd882ef03ee3/monai_bathymetry_data.txt", + joinpath(@__DIR__, "monai_bathymetry_data.txt")); + +# B-spline interpolation of the underlying data +bath_spline_struct = BicubicBSpline(spline_bathymetry_file, end_condition = "not-a-knot") +bathymetry(x, y) = spline_interpolation(bath_spline_struct, x, y) + +# Initial condition is a constant background state +function initial_condition_monai_tsunami(x, t, equations::ShallowWaterMultiLayerEquations2D) + # Set the background water height + H = [equations.H0] + + # Initially water is at rest + v1 = zero(H) + v2 = zero(H) + + # Bottom topography values are computed from the bicubic spline created above + x1, x2 = x + b = bathymetry(x1, x2) + + # 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 the conservative variables + return prim2cons(SVector(H..., v1..., v2..., b), equations) +end + +initial_condition = initial_condition_monai_tsunami + +# Tsunami test case uses a specialized wave maker type of boundary condition. +# It is used to model an incident wave that approaches from off-shore +# with a water depth of $h = 13.535\,\text{cm}$. To create the incident wave information +# that is valid over the time interval $t \in [0\,s, 22.5\,s]$. + +# We download the incident wave data that has been preprocessed to be in the format +# required by TrixiBottomTopography. +water_height_data = Trixi.download("https://gist.githubusercontent.com/andrewwinters5000/5b11f5f175bddb326d11d8e28398127e/raw/64980e0e4526e0fcd49589b34ee5458b9a1cebff/monai_wavemaker_bc.txt", + joinpath(@__DIR__, "monai_wavemaker_bc.txt")); + +# Similar to the bathymetry approximation, we construct a cubic B-spline interpolation +# of the data, then create a function to evaluate the resulting spline at a given $t$ value. +h_spline_struct = CubicBSpline(water_height_data; end_condition = "not-a-knot") +H_from_wave_maker(t) = spline_interpolation(h_spline_struct, t) + +# Now we can define the specialized boundary condition for the incident wave maker. +function boundary_condition_wave_maker(u_inner, normal_direction::AbstractVector, x, t, + surface_flux_functions, + equations::ShallowWaterMultiLayerEquations2D) + surface_flux_function, nonconservative_flux_function = surface_flux_functions + + # Compute the water height from the wave maker input file data + # and then clip to avoid negative water heights and division by zero + h_ext = max(equations.threshold_limiter, H_from_wave_maker(t) - u_inner[4]) + + # Compute the incoming velocity as in Eq. (10) of the paper + # - S. Vater, N. Beisiegel, and J. Behrens (2019) + # A limiter-based well-balanced discontinuous Galerkin method for shallow-water flows + # with wetting and drying: Triangular grids + # [DOI: 10.1002/fld.4762](https://doi.org/10.1002/fld.4762) + h0 = 0.13535 # reference incident water height converted to meters + v1_ext = 2 * (sqrt(equations.gravity * h_ext) - sqrt(equations.gravity * h0)) + + # Create the external solution state in the conservative variables + u_outer = SVector([h_ext]..., [h_ext * v1_ext]..., [zero(eltype(x))]..., u_inner[4]) + + # Calculate the boundary flux and nonconservative contributions + flux = surface_flux_function(u_inner, u_outer, normal_direction, equations) + + noncons_flux = nonconservative_flux_function(u_inner, u_outer, normal_direction, + equations) + + return flux, noncons_flux +end + +boundary_condition = Dict(:Bottom => boundary_condition_slip_wall, + :Top => boundary_condition_slip_wall, + :Right => boundary_condition_slip_wall, + :Left => boundary_condition_wave_maker) + +# Manning friction source term +@inline function source_terms_manning_friction(u, x, t, + equations::ShallowWaterMultiLayerEquations2D) + # Grab the conservative variables + h, hv_1, hv_2, b = u + + # Desingularization + sh = (2.0 * h[1]) / (h[1]^2 + max(h[1]^2, 1e-8)) + + # friction coefficient + n = 0.001 + + # Compute the common friction term + Sf = -equations.gravity * n^2 * sh^(7.0 / 3.0) * sqrt(hv_1[1]^2 + hv_2[1]^2) + + return SVector(zero(h)..., + Sf * hv_1[1]..., + Sf * hv_2[1]..., + zero(b)) +end + +############################################################################### +# 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) + +# Cannot simply use `waterheight` here for multilayer equations. +# Need a helper function to extract the relevant variable. +@inline function main_waterheight(u, equations) + return waterheight(u, equations)[1] +end + +indicator_sc = IndicatorHennemannGassnerShallowWater(equations, basis, + alpha_max = 0.5, + alpha_min = 0.001, + alpha_smooth = true, + variable = main_waterheight) +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/d26356bdc8e8d742a3035b3f71c71a68/raw/9d6ceedb844e92313d1dac2318a28c87ffbb9de2/mesh_monai_shore.inp", + joinpath(@__DIR__, "mesh_monai_shore.inp")); + +mesh = P4estMesh{2}(mesh_file) + +# create the semi discretization object +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver; + boundary_conditions = boundary_condition, + source_terms = source_terms_manning_friction) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 22.5) +ode = semidiscretize(semi, tspan) + +############################################################################### +# Callbacks + +summary_callback = SummaryCallback() + +analysis_interval = 1000 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + extra_analysis_errors = (:conservation_error,), + save_analysis = false) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(dt = 0.2, + save_initial_solution = true, + save_final_solution = true) + +# TimeSeries Callback is currently unavailable on `P4estMesh` +# time_series = TimeSeriesCallback(semi, [(4.521, 1.196), # G5 +# (4.521, 1.696), # G7 +# (4.521, 2.196)]; # G9 +# interval = 2, +# solution_variables=cons2prim) + +stepsize_callback = StepsizeCallback(cfl = 0.5) + +# Another possible AMR indicator function could be the velocity, such that it only fires +# in regions where the water is moving +@inline function momentum_norm(u, equations::ShallowWaterMultiLayerEquations2D) + _, h_v1, h_v2, _ = u + return sqrt(h_v1[1]^2 + h_v2[1]^2) + # h, h_v1, h_v2, _ = u + # v1 = (2 * h[1] * h_v1[1]) / (h[1]^2 + max(h[1]^2, 1e-8)) + # v2 = (2 * h[1] * h_v2[1]) / (h[1]^2 + max(h[1]^2, 1e-8)) + # if h[1] <= 1e-8 + # v1 = 0.0 + # v2 = 0.0 + # end + # if v1 > 0.525 || v2 > 0.525 + # v1 = 0.0 + # v2 = 0.0 + # end + # return sqrt(v1^2 + v2^2) +end + +amr_controller = ControllerThreeLevel(semi, + IndicatorMax(semi, variable = momentum_norm), + base_level = 0, + med_level = 1, med_threshold = 0.012, + max_level = 3, max_threshold = 0.015) +# velocity is very sensitive to use as an indicator +# med_level = 1, med_threshold = 0.25, +# max_level = 2, max_threshold = 0.475) + +amr_callback = AMRCallback(semi, amr_controller, + interval = 5, + adapt_initial_condition = false, + adapt_initial_condition_only_refine = false) + +callbacks = CallbackSet(summary_callback, + analysis_callback, + alive_callback, + # time_series, + amr_callback, + save_solution, + stepsize_callback) + +############################################################################### +# run the simulation + +stage_limiter! = PositivityPreservingLimiterShallowWater(variables = (waterheight,)) + +sol = solve(ode, SSPRK43(stage_limiter!); + ode_default_options()..., + callback = callbacks, + adaptive = false, + dt = 1.0); diff --git a/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_three_mound_dam_break.jl b/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_three_mound_dam_break.jl new file mode 100644 index 00000000..c2c13999 --- /dev/null +++ b/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_three_mound_dam_break.jl @@ -0,0 +1,182 @@ + +using OrdinaryDiffEqSSPRK, OrdinaryDiffEqLowStorageRK +using Trixi +using TrixiShallowWater + +############################################################################### +# Semidiscretization of the multilayer shallow water equations with one layer +# to fallback to be the standard shallow water equations + +equations = ShallowWaterMultiLayerEquations2D(gravity = 9.81, H0 = 0.0, + threshold_limiter = 1e-12, + rhos = (1.0)) + +""" + initial_condition_three_mounds(x, t, equations::ShallowWaterMultiLayerEquations2D) + +Initial condition simulating a dam break. The bottom topography is given by one large and two smaller +mounds. The mounds are flooded by the water for t > 0. To smooth the discontinuity, a logistic function +is applied. + +The initial conditions is taken from Section 6.3 of the paper: +- Niklas Wintermeyer, Andrew R. Winters, Gregor J. Gassner and Timothy Warburton (2018) + An entropy stable discontinuous Galerkin method for the shallow water equations on + curvilinear meshes with wet/dry fronts accelerated by GPUs\n + [DOI: 10.1016/j.jcp.2018.08.038](https://doi.org/10.1016/j.jcp.2018.08.038) +""" +function initial_condition_three_mounds(x, t, equations::ShallowWaterMultiLayerEquations2D) + + x1, x2 = x + M_1 = 1 - 0.1 * sqrt((x1 - 30.0)^2 + (x2 - 22.5)^2) + M_2 = 1 - 0.1 * sqrt((x1 - 30.0)^2 + (x2 - 7.5)^2) + M_3 = 2.8 - 0.28 * sqrt((x1 - 47.5)^2 + (x2 - 15.0)^2) + + b = max(0.0, M_1, M_2, M_3) + + # use a logistic function to transfer water height value smoothly + L = 1.875 # maximum of function + x0 = 8 # center point of function + k = -75.0 # sharpness of transfer + + H = [max(b, L / (1.0 + exp(-k * (x1 - x0))))] + + # Everything is initially not moving + 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 the conservative variables + return prim2cons(SVector(H..., v1..., v2..., b), equations) +end + +initial_condition = initial_condition_three_mounds + +function boundary_condition_outflow(u_inner, normal_direction::AbstractVector, x, t, + surface_flux_functions, + equations::ShallowWaterMultiLayerEquations2D) + surface_flux_function, nonconservative_flux_function = surface_flux_functions + # Impulse and bottom from inside, height from external state + u_outer = SVector(equations.threshold_limiter, u_inner[2], u_inner[3], u_inner[4]) + + # calculate the boundary flux + flux = surface_flux_function(u_inner, u_outer, normal_direction, equations) + noncons_flux = nonconservative_flux_function(u_inner, u_outer, normal_direction, + equations) + return flux, noncons_flux +end + +boundary_conditions = Dict(:Bottom => boundary_condition_slip_wall, + :Top => boundary_condition_slip_wall, + :Right => boundary_condition_outflow, + :Left => 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(4) + +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/7a34a0116dc84b69449ea25054b7585a/raw/e098000588ac5d8b543fc807b7fcc49df63dc7a0/mesh_three_mound.inp", + joinpath(@__DIR__, "mesh_three_mound.inp")); + +mesh = P4estMesh{2}(mesh_file) + +# Create the semi discretization object +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver; + boundary_conditions = boundary_conditions) + +############################################################################### +# ODE solver + +tspan = (0.0, 20.0) +ode = semidiscretize(semi, tspan) + +############################################################################### +# Callbacks + +summary_callback = SummaryCallback() + +analysis_interval = 1000 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + extra_analysis_errors = (:conservation_error,), + save_analysis = false) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(dt = 0.5, + save_initial_solution = true, + save_final_solution = true) + +stepsize_callback = StepsizeCallback(cfl = 0.5) + +# Cannot simply use `waterheight` here for multilayer equations. +# Need a helper function to extract the relevant variable. +# TODO: This AMR indicator fires too often in the dry regions +@inline function main_waterheight(u, equations) + return waterheight(u, equations)[1] +end + +amr_indicator = IndicatorLöhner(semi, variable = main_waterheight) + +amr_controller = ControllerThreeLevel(semi, amr_indicator, + base_level = 0, + med_level = 1, med_threshold = 0.1, + max_level = 3, max_threshold = 0.25) + +amr_callback = AMRCallback(semi, amr_controller, + interval = 5, + adapt_initial_condition = false, + adapt_initial_condition_only_refine = false) + +callbacks = CallbackSet(summary_callback, + analysis_callback, + alive_callback, + amr_callback, + save_solution, + stepsize_callback) + +############################################################################### +# run the simulation + +stage_limiter! = PositivityPreservingLimiterShallowWater(variables = (waterheight,)) + +sol = solve(ode, SSPRK43(stage_limiter!); + ode_default_options()..., + callback = callbacks, + adaptive = false, + dt = 1.0 # solve needs some value here; overwritten by the `stepsize_callback` + ); diff --git a/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_well_balanced_wet_dry_nonconforming.jl b/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_well_balanced_wet_dry_nonconforming.jl index 455db40d..2bfe02df 100644 --- a/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_well_balanced_wet_dry_nonconforming.jl +++ b/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_well_balanced_wet_dry_nonconforming.jl @@ -14,7 +14,8 @@ equations = ShallowWaterMultiLayerEquations2D(gravity = 9.812, H0 = 1.235, # 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) +function initial_condition_well_balancedness(x, t, + equations::ShallowWaterMultiLayerEquations2D) # Set the background values H = [equations.H0] From dd1783195c39d9af6458216866abc98cec9e71ca Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Fri, 16 May 2025 23:28:04 +0200 Subject: [PATCH 14/37] update file name --- ...> elixir_shallowwater_multilayer_three_mound_dam_break_amr.jl} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename examples/p4est_2d_dgsem/{elixir_shallowwater_multilayer_three_mound_dam_break.jl => elixir_shallowwater_multilayer_three_mound_dam_break_amr.jl} (100%) diff --git a/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_three_mound_dam_break.jl b/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_three_mound_dam_break_amr.jl similarity index 100% rename from examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_three_mound_dam_break.jl rename to examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_three_mound_dam_break_amr.jl From c485695c96d4cc30e92c7860f5f63422705593bc Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Sat, 17 May 2025 02:34:33 +0200 Subject: [PATCH 15/37] formatter --- ...lixir_shallowwater_multilayer_three_mound_dam_break_amr.jl | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_three_mound_dam_break_amr.jl b/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_three_mound_dam_break_amr.jl index c2c13999..639faa15 100644 --- a/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_three_mound_dam_break_amr.jl +++ b/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_three_mound_dam_break_amr.jl @@ -25,7 +25,6 @@ The initial conditions is taken from Section 6.3 of the paper: [DOI: 10.1016/j.jcp.2018.08.038](https://doi.org/10.1016/j.jcp.2018.08.038) """ function initial_condition_three_mounds(x, t, equations::ShallowWaterMultiLayerEquations2D) - x1, x2 = x M_1 = 1 - 0.1 * sqrt((x1 - 30.0)^2 + (x2 - 22.5)^2) M_2 = 1 - 0.1 * sqrt((x1 - 30.0)^2 + (x2 - 7.5)^2) @@ -178,5 +177,4 @@ sol = solve(ode, SSPRK43(stage_limiter!); ode_default_options()..., callback = callbacks, adaptive = false, - dt = 1.0 # solve needs some value here; overwritten by the `stepsize_callback` - ); + dt = 1.0); From 2d0c5ab3b8f2d5b42a2e1455f8b6d8d1a598bbc8 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Sat, 17 May 2025 09:19:08 +0200 Subject: [PATCH 16/37] add three mound elixir to testing --- test/test_p4est_2d.jl | 26 ++++++++++++++++++++++++++ 1 file changed, 26 insertions(+) diff --git a/test/test_p4est_2d.jl b/test/test_p4est_2d.jl index 6272af76..e6c3fc6c 100644 --- a/test/test_p4est_2d.jl +++ b/test/test_p4est_2d.jl @@ -209,6 +209,32 @@ end # SWE @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 end end + + @trixi_testset "elixir_shallowwater_multilayer_three_mound_dam_break_amr.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_multilayer_three_mound_dam_break_amr.jl"), + l2=[ + 0.13676802860182155, + 0.38852715519450065, + 1.8293078549048694e-13, + 0.0019606480801663806], + linf=[ + 1.0877121083544241, + 2.5590836440248874, + 1.7456605305379413e-11, + 0.044079775502972124], + tspan=(0.0, 0.025), + coverage_override=(maxiters = 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 # P4estMesh2D From 2459e4953197f190a52c7d333469e19688818e07 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Sat, 17 May 2025 09:21:39 +0200 Subject: [PATCH 17/37] better comments in new elixirs --- .../elixir_shallowwater_multilayer_monai_flood_amr.jl | 4 ++-- ...lixir_shallowwater_multilayer_three_mound_dam_break_amr.jl | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_monai_flood_amr.jl b/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_monai_flood_amr.jl index c8e599f3..be8e2fc3 100644 --- a/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_monai_flood_amr.jl +++ b/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_monai_flood_amr.jl @@ -255,5 +255,5 @@ stage_limiter! = PositivityPreservingLimiterShallowWater(variables = (waterheigh sol = solve(ode, SSPRK43(stage_limiter!); ode_default_options()..., callback = callbacks, - adaptive = false, - dt = 1.0); + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + adaptive = false); diff --git a/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_three_mound_dam_break_amr.jl b/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_three_mound_dam_break_amr.jl index 639faa15..071e7699 100644 --- a/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_three_mound_dam_break_amr.jl +++ b/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_three_mound_dam_break_amr.jl @@ -176,5 +176,5 @@ stage_limiter! = PositivityPreservingLimiterShallowWater(variables = (waterheigh sol = solve(ode, SSPRK43(stage_limiter!); ode_default_options()..., callback = callbacks, - adaptive = false, - dt = 1.0); + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + adaptive = false); From 53b7dc8c1fbdaa86fa9ee7271006b110f435232f Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Sat, 17 May 2025 20:33:50 +0200 Subject: [PATCH 18/37] adjust amr controller values in three mound test case --- ...xir_shallowwater_multilayer_three_mound_dam_break_amr.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_three_mound_dam_break_amr.jl b/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_three_mound_dam_break_amr.jl index 071e7699..891bfd67 100644 --- a/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_three_mound_dam_break_amr.jl +++ b/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_three_mound_dam_break_amr.jl @@ -140,7 +140,7 @@ save_solution = SaveSolutionCallback(dt = 0.5, save_initial_solution = true, save_final_solution = true) -stepsize_callback = StepsizeCallback(cfl = 0.5) +stepsize_callback = StepsizeCallback(cfl = 0.4) # Cannot simply use `waterheight` here for multilayer equations. # Need a helper function to extract the relevant variable. @@ -153,8 +153,8 @@ amr_indicator = IndicatorLöhner(semi, variable = main_waterheight) amr_controller = ControllerThreeLevel(semi, amr_indicator, base_level = 0, - med_level = 1, med_threshold = 0.1, - max_level = 3, max_threshold = 0.25) + med_level = 1, med_threshold = 0.25, + max_level = 3, max_threshold = 0.5) amr_callback = AMRCallback(semi, amr_controller, interval = 5, From 4657abc6a228ad486634f6d109a65b268c4e86c2 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Sat, 17 May 2025 23:46:09 +0200 Subject: [PATCH 19/37] update three mound test values --- test/test_p4est_2d.jl | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/test/test_p4est_2d.jl b/test/test_p4est_2d.jl index e6c3fc6c..16b9d83c 100644 --- a/test/test_p4est_2d.jl +++ b/test/test_p4est_2d.jl @@ -214,16 +214,16 @@ end # SWE @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_multilayer_three_mound_dam_break_amr.jl"), l2=[ - 0.13676802860182155, - 0.38852715519450065, - 1.8293078549048694e-13, - 0.0019606480801663806], + 0.13142461271166225, + 0.35364621879644254, + 4.778520542281107e-13, + 0.0019606480801664917], linf=[ - 1.0877121083544241, - 2.5590836440248874, - 1.7456605305379413e-11, + 1.1338438044790615, + 2.416846787140322, + 6.256251416535554e-11, 0.044079775502972124], - tspan=(0.0, 0.025), + tspan=(0.0, 0.3), coverage_override=(maxiters = 5,)) # Ensure that we do not have excessive memory allocations From a266921d8cdb4a0692da62d3a9ec76a666936c84 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Mon, 19 May 2025 15:42:35 +0200 Subject: [PATCH 20/37] fix type stability dispatch on spline functions in Monai elixir --- ...shallowwater_multilayer_monai_flood_amr.jl | 23 +++++++++++-------- 1 file changed, 13 insertions(+), 10 deletions(-) diff --git a/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_monai_flood_amr.jl b/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_monai_flood_amr.jl index be8e2fc3..db8894a9 100644 --- a/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_monai_flood_amr.jl +++ b/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_monai_flood_amr.jl @@ -17,16 +17,18 @@ using TrixiBottomTopography equations = ShallowWaterMultiLayerEquations2D(gravity = 9.81, H0 = 0.0, rhos = (1.0)) -# Setup the bottom topography data and cubic spline interpolant +# Get the precision type for type stable dispatch on the helper spline functions +RealT = typeof(equations.gravity) # Create a bicubic B-spline interpolation of the bathymetry data, then create a function # to evaluate the resulting spline at a given point $(x,y)$. spline_bathymetry_file = Trixi.download("https://gist.githubusercontent.com/andrewwinters5000/21255c980c4eda5294f91e8dfe6c7e33/raw/1afb73928892774dc3a902e0c46ffd882ef03ee3/monai_bathymetry_data.txt", joinpath(@__DIR__, "monai_bathymetry_data.txt")); -# B-spline interpolation of the underlying data -bath_spline_struct = BicubicBSpline(spline_bathymetry_file, end_condition = "not-a-knot") -bathymetry(x, y) = spline_interpolation(bath_spline_struct, x, y) +# B-spline interpolation of the underlying data. +# The type of struct is fixed to be `BicubicBSpline` for type stability +const bath_spline_struct = BicubicBSpline(spline_bathymetry_file, end_condition = "not-a-knot") +bathymetry(x::RealT, y::RealT) = spline_interpolation(bath_spline_struct, x, y) # Initial condition is a constant background state function initial_condition_monai_tsunami(x, t, equations::ShallowWaterMultiLayerEquations2D) @@ -73,13 +75,14 @@ water_height_data = Trixi.download("https://gist.githubusercontent.com/andrewwin # Similar to the bathymetry approximation, we construct a cubic B-spline interpolation # of the data, then create a function to evaluate the resulting spline at a given $t$ value. -h_spline_struct = CubicBSpline(water_height_data; end_condition = "not-a-knot") -H_from_wave_maker(t) = spline_interpolation(h_spline_struct, t) +# The type of struct is fixed to be `CubicBSpline` for type stability +const h_spline_struct = CubicBSpline(water_height_data; end_condition = "not-a-knot") +H_from_wave_maker(t::RealT) = spline_interpolation(h_spline_struct, t) # Now we can define the specialized boundary condition for the incident wave maker. -function boundary_condition_wave_maker(u_inner, normal_direction::AbstractVector, x, t, - surface_flux_functions, - equations::ShallowWaterMultiLayerEquations2D) +@inline function boundary_condition_wave_maker(u_inner, normal_direction::AbstractVector, x, t, + surface_flux_functions, + equations::ShallowWaterMultiLayerEquations2D) surface_flux_function, nonconservative_flux_function = surface_flux_functions # Compute the water height from the wave maker input file data @@ -95,7 +98,7 @@ function boundary_condition_wave_maker(u_inner, normal_direction::AbstractVector v1_ext = 2 * (sqrt(equations.gravity * h_ext) - sqrt(equations.gravity * h0)) # Create the external solution state in the conservative variables - u_outer = SVector([h_ext]..., [h_ext * v1_ext]..., [zero(eltype(x))]..., u_inner[4]) + u_outer = SVector(h_ext, h_ext * v1_ext, zero(eltype(x)), u_inner[4]) # Calculate the boundary flux and nonconservative contributions flux = surface_flux_function(u_inner, u_outer, normal_direction, equations) From 65d5f6e3c2532b454151d7a2a5f5fae88c90b510 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Mon, 19 May 2025 15:46:51 +0200 Subject: [PATCH 21/37] adjust the spline constructors in the monai tutorial as well --- .../tutorials/elixir_shallowwater_monai_tsunami.jl | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/docs/src/tutorials/elixir_shallowwater_monai_tsunami.jl b/docs/src/tutorials/elixir_shallowwater_monai_tsunami.jl index ab534e23..0801920a 100644 --- a/docs/src/tutorials/elixir_shallowwater_monai_tsunami.jl +++ b/docs/src/tutorials/elixir_shallowwater_monai_tsunami.jl @@ -156,7 +156,7 @@ generate_mesh(monai); # and specify the gravitational acceleration to `gravity = 9.812` # as well as a background water height `H0 = 0.0`. # In contrast to the [`Trixi.ShallowWaterEquations2D`](@extref Trixi.ShallowWaterEquations2D) type, -# this equation type allows contains additional parameters and methods needed to handle wetting and drying. +# this equation type contains additional parameters and methods needed to handle wetting and drying. equations = ShallowWaterEquationsWetDry2D(gravity = 9.81, H0 = 0.0) # Next, we construct an approximation to the bathymetry with TrixiBottomTopography.jl using @@ -169,8 +169,9 @@ spline_bathymetry_file = Trixi.download("https://gist.githubusercontent.com/andr # Create a bicubic B-spline interpolation of the bathymetry data, then create a function # to evaluate the resulting spline at a given point $(x,y)$. -bath_spline_struct = BicubicBSpline(spline_bathymetry_file, end_condition = "not-a-knot") -bathymetry(x, y) = spline_interpolation(bath_spline_struct, x, y); +# The type of this struct is fixed to be `BicubicBSpline`. +const bath_spline_struct = BicubicBSpline(spline_bathymetry_file, end_condition = "not-a-knot") +bathymetry(x::Float64, y::Float64) = spline_interpolation(bath_spline_struct, x, y); # We then create a function to supply the initial condition for the simulation. @inline function initial_condition_monai_tsunami(x, t, @@ -211,8 +212,9 @@ wavemaker_bc_file = Trixi.download("https://gist.githubusercontent.com/andrewwin # Similar to the bathymetry approximation, we construct a cubic B-spline interpolation # of the data, then create a function to evaluate the resulting spline at a given $t$ value. -h_spline_struct = CubicBSpline(wavemaker_bc_file; end_condition = "not-a-knot") -H_from_wave_maker(t) = spline_interpolation(h_spline_struct, t); +# The type of this struct is fixed to be `CubicBSpline`. +const h_spline_struct = CubicBSpline(wavemaker_bc_file; end_condition = "not-a-knot") +H_from_wave_maker(t::Float64) = spline_interpolation(h_spline_struct, t); # Now we are equipped to define the specialized boundary condition for the incident # wave maker. From 88814af5f1ff3060434425d10e793a4a80aa0109 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Mon, 19 May 2025 15:52:02 +0200 Subject: [PATCH 22/37] adjust comments --- docs/src/tutorials/elixir_shallowwater_monai_tsunami.jl | 4 ++-- .../elixir_shallowwater_multilayer_monai_flood_amr.jl | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/docs/src/tutorials/elixir_shallowwater_monai_tsunami.jl b/docs/src/tutorials/elixir_shallowwater_monai_tsunami.jl index 0801920a..f4087b79 100644 --- a/docs/src/tutorials/elixir_shallowwater_monai_tsunami.jl +++ b/docs/src/tutorials/elixir_shallowwater_monai_tsunami.jl @@ -169,7 +169,7 @@ spline_bathymetry_file = Trixi.download("https://gist.githubusercontent.com/andr # Create a bicubic B-spline interpolation of the bathymetry data, then create a function # to evaluate the resulting spline at a given point $(x,y)$. -# The type of this struct is fixed to be `BicubicBSpline`. +# The type of this struct is fixed as `BicubicBSpline`. const bath_spline_struct = BicubicBSpline(spline_bathymetry_file, end_condition = "not-a-knot") bathymetry(x::Float64, y::Float64) = spline_interpolation(bath_spline_struct, x, y); @@ -212,7 +212,7 @@ wavemaker_bc_file = Trixi.download("https://gist.githubusercontent.com/andrewwin # Similar to the bathymetry approximation, we construct a cubic B-spline interpolation # of the data, then create a function to evaluate the resulting spline at a given $t$ value. -# The type of this struct is fixed to be `CubicBSpline`. +# The type of this struct is fixed as `CubicBSpline`. const h_spline_struct = CubicBSpline(wavemaker_bc_file; end_condition = "not-a-knot") H_from_wave_maker(t::Float64) = spline_interpolation(h_spline_struct, t); diff --git a/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_monai_flood_amr.jl b/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_monai_flood_amr.jl index db8894a9..9569ba76 100644 --- a/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_monai_flood_amr.jl +++ b/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_monai_flood_amr.jl @@ -26,7 +26,7 @@ spline_bathymetry_file = Trixi.download("https://gist.githubusercontent.com/andr joinpath(@__DIR__, "monai_bathymetry_data.txt")); # B-spline interpolation of the underlying data. -# The type of struct is fixed to be `BicubicBSpline` for type stability +# The type of this struct is fixed as `BicubicBSpline`. const bath_spline_struct = BicubicBSpline(spline_bathymetry_file, end_condition = "not-a-knot") bathymetry(x::RealT, y::RealT) = spline_interpolation(bath_spline_struct, x, y) @@ -75,7 +75,7 @@ water_height_data = Trixi.download("https://gist.githubusercontent.com/andrewwin # Similar to the bathymetry approximation, we construct a cubic B-spline interpolation # of the data, then create a function to evaluate the resulting spline at a given $t$ value. -# The type of struct is fixed to be `CubicBSpline` for type stability +# The type of this struct is fixed as `CubicBSpline`. const h_spline_struct = CubicBSpline(water_height_data; end_condition = "not-a-knot") H_from_wave_maker(t::RealT) = spline_interpolation(h_spline_struct, t) From 16f0cfbdca3171fb80b1de4a89fa70d03e151933 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Mon, 19 May 2025 19:44:54 +0200 Subject: [PATCH 23/37] Apply suggestions from code review Co-authored-by: Patrick Ersing <114223904+patrickersing@users.noreply.github.com> --- ...shallowwater_multilayer_monai_flood_amr.jl | 25 ++++++----------- ...ter_multilayer_perturbation_wet_dry_amr.jl | 28 ++++++------------- ...er_multilayer_three_mound_dam_break_amr.jl | 14 +++------- ...yer_well_balanced_wet_dry_nonconforming.jl | 26 +++++------------ src/solvers/dgsem_p4est/dg_2d.jl | 8 ++++-- 5 files changed, 34 insertions(+), 67 deletions(-) diff --git a/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_monai_flood_amr.jl b/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_monai_flood_amr.jl index 9569ba76..5f089a4d 100644 --- a/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_monai_flood_amr.jl +++ b/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_monai_flood_amr.jl @@ -15,7 +15,7 @@ using TrixiBottomTopography # to fallback to be the standard shallow water equations equations = ShallowWaterMultiLayerEquations2D(gravity = 9.81, H0 = 0.0, - rhos = (1.0)) + rhos = 1.0) # Get the precision type for type stable dispatch on the helper spline functions RealT = typeof(equations.gravity) @@ -33,7 +33,7 @@ bathymetry(x::RealT, y::RealT) = spline_interpolation(bath_spline_struct, x, y) # Initial condition is a constant background state function initial_condition_monai_tsunami(x, t, equations::ShallowWaterMultiLayerEquations2D) # Set the background water height - H = [equations.H0] + H = equations.H0 # Initially water is at rest v1 = zero(H) @@ -49,16 +49,10 @@ function initial_condition_monai_tsunami(x, t, equations::ShallowWaterMultiLayer # 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 + H = max(H, b + equations.threshold_limiter) # Return the conservative variables - return prim2cons(SVector(H..., v1..., v2..., b), equations) + return prim2cons(SVector(H, v1, v2, b), equations) end initial_condition = initial_condition_monai_tsunami @@ -118,20 +112,20 @@ boundary_condition = Dict(:Bottom => boundary_condition_slip_wall, @inline function source_terms_manning_friction(u, x, t, equations::ShallowWaterMultiLayerEquations2D) # Grab the conservative variables - h, hv_1, hv_2, b = u + h, hv_1, hv_2, _ = u # Desingularization - sh = (2.0 * h[1]) / (h[1]^2 + max(h[1]^2, 1e-8)) + h = (h^2 + max(h^2, 1e-8)) / (2.0 * h) # friction coefficient n = 0.001 # Compute the common friction term - Sf = -equations.gravity * n^2 * sh^(7.0 / 3.0) * sqrt(hv_1[1]^2 + hv_2[1]^2) + Sf = -equations.gravity * n^2 * h^(-7 / 3) * sqrt(hv_1^2 + hv_2^2) return SVector(zero(h)..., - Sf * hv_1[1]..., - Sf * hv_2[1]..., + Sf * hv_1, + Sf * hv_2, zero(b)) end @@ -245,7 +239,6 @@ amr_callback = AMRCallback(semi, amr_controller, callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, - # time_series, amr_callback, save_solution, stepsize_callback) diff --git a/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_perturbation_wet_dry_amr.jl b/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_perturbation_wet_dry_amr.jl index 60df89e8..9545f474 100644 --- a/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_perturbation_wet_dry_amr.jl +++ b/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_perturbation_wet_dry_amr.jl @@ -12,11 +12,11 @@ using TrixiShallowWater # on a nonconforming mesh with AMR equations = ShallowWaterMultiLayerEquations2D(gravity = 9.812, H0 = 1.235, - rhos = (1.0)) + rhos = 1.0) function initial_condition_perturbation(x, t, equations::ShallowWaterMultiLayerEquations2D) # Calculate primitive variables - H = [equations.H0] + H = equations.H0 v1 = zero(H) v2 = zero(H) @@ -33,15 +33,9 @@ function initial_condition_perturbation(x, t, equations::ShallowWaterMultiLayerE # 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 + H = max(H, b + equations.threshold_limiter) - return prim2cons(SVector(H..., v1..., v2..., b), equations) + return prim2cons(SVector(H, v1, v2, b), equations) end initial_condition = initial_condition_perturbation @@ -133,7 +127,7 @@ ode = semidiscretize(semi, tspan) function initial_condition_discontinuous_perturbation(x, t, element_id, equations::ShallowWaterMultiLayerEquations2D) # Set the background values for velocity - H = [equations.H0] + H = equations.H0 v1 = zero(H) v2 = zero(H) @@ -160,7 +154,7 @@ function initial_condition_discontinuous_perturbation(x, t, element_id, # Put in a discontinuous perturbation using the element number if element_id in [232, 224, 225, 226, 227, 228, 229, 230] - H = H .+ 1.6 + H = H + 1.6 end # It is mandatory to shift the water level at dry areas to make sure the water height h @@ -169,15 +163,9 @@ function initial_condition_discontinuous_perturbation(x, t, element_id, # 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 + H = max(H, b + equations.threshold_limiter) - return prim2cons(SVector(H..., v1..., v2..., b), equations) + return prim2cons(SVector(H, v1, v2, b), equations) end # point to the data we want to augment diff --git a/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_three_mound_dam_break_amr.jl b/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_three_mound_dam_break_amr.jl index 891bfd67..def8891d 100644 --- a/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_three_mound_dam_break_amr.jl +++ b/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_three_mound_dam_break_amr.jl @@ -9,7 +9,7 @@ using TrixiShallowWater equations = ShallowWaterMultiLayerEquations2D(gravity = 9.81, H0 = 0.0, threshold_limiter = 1e-12, - rhos = (1.0)) + rhos = 1.0) """ initial_condition_three_mounds(x, t, equations::ShallowWaterMultiLayerEquations2D) @@ -37,7 +37,7 @@ function initial_condition_three_mounds(x, t, equations::ShallowWaterMultiLayerE x0 = 8 # center point of function k = -75.0 # sharpness of transfer - H = [max(b, L / (1.0 + exp(-k * (x1 - x0))))] + H = max(b, L / (1.0 + exp(-k * (x1 - x0)))) # Everything is initially not moving v1 = zero(H) @@ -49,16 +49,10 @@ function initial_condition_three_mounds(x, t, equations::ShallowWaterMultiLayerE # 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 + H = max(H, b + equations.threshold_limiter) # Return the conservative variables - return prim2cons(SVector(H..., v1..., v2..., b), equations) + return prim2cons(SVector(H, v1, v2, b), equations) end initial_condition = initial_condition_three_mounds diff --git a/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_well_balanced_wet_dry_nonconforming.jl b/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_well_balanced_wet_dry_nonconforming.jl index 2bfe02df..0407a660 100644 --- a/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_well_balanced_wet_dry_nonconforming.jl +++ b/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_well_balanced_wet_dry_nonconforming.jl @@ -9,7 +9,7 @@ using TrixiShallowWater # Test case for well-balancedness with wet/dry transition elements on a nonconforming mesh. equations = ShallowWaterMultiLayerEquations2D(gravity = 9.812, H0 = 1.235, - rhos = (1.0)) + rhos = 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 @@ -18,7 +18,7 @@ function initial_condition_well_balancedness(x, t, equations::ShallowWaterMultiLayerEquations2D) # Set the background values - H = [equations.H0] + H = equations.H0 v1 = zero(H) v2 = zero(H) @@ -35,15 +35,9 @@ function initial_condition_well_balancedness(x, t, # 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 + H = max(H, b + equations.threshold_limiter) - return prim2cons(SVector(H..., v1..., v2..., b), equations) + return prim2cons(SVector(H, v1, v2, b), equations) end initial_condition = initial_condition_well_balancedness @@ -140,7 +134,7 @@ ode = semidiscretize(semi, tspan) function initial_condition_discontinuous_well_balancedness(x, t, element_id, equations::ShallowWaterMultiLayerEquations2D) # Set the background values - H = [equations.H0] + H = equations.H0 v1 = zero(H) v2 = zero(H) @@ -167,15 +161,9 @@ function initial_condition_discontinuous_well_balancedness(x, t, element_id, # 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 + H = max(H, b + equations.threshold_limiter) - return prim2cons(SVector(H..., v1..., v2..., b), equations) + return prim2cons(SVector(H, v1, v2, b), equations) end # point to the data we want to augment diff --git a/src/solvers/dgsem_p4est/dg_2d.jl b/src/solvers/dgsem_p4est/dg_2d.jl index 6d9f56f0..4fb4ea24 100644 --- a/src/solvers/dgsem_p4est/dg_2d.jl +++ b/src/solvers/dgsem_p4est/dg_2d.jl @@ -194,7 +194,7 @@ end # Specialized strategy to project the solution to the mortars and preserve # well-balancedness from the work of Benov et al. (https://doi.org/10.1016/j.jcp.2018.02.008). -# This version is for the `ShallowWaterMultiLayerEquations2D` which "peels" the pressure +# This version is for the `ShallowWaterMultiLayerEquations2D` with a single layer, which "peels" the pressure # contribution into the nonconservative term for easier well-balancing, # see Ersing et al. (https://doi.org/10.1016/j.jcp.2025.113802). # Thus, this uses the original `P4estMortarContainer`, `calc_mortar_flux!` @@ -210,6 +210,11 @@ function Trixi.prolong2mortars!(cache, u, @unpack neighbor_ids, node_indices = cache.mortars index_range = eachnode(dg) + # Raise an error if `prolong2mortars!` is called with multiple layers. + if nlayers(equations) != 1 + error("Non-conforming meshes are only supported for a single layer!") + end + Trixi.@threaded for mortar in Trixi.eachmortar(dg, cache) # Copy solution data from the small elements using "delayed indexing" with # a start value and a step size to get the correct face and orientation. @@ -266,7 +271,6 @@ function Trixi.prolong2mortars!(cache, u, # Note, some shift is required to ensure we catch water heights close to the threshold # and maintain well-balancedness. The larger this shift we also maintain the conservation # errors to be close to double precision roundoff for longer. - # if u[1, i_large, j_large, element] >= 500 * equations.threshold_limiter # ≈ 5.5e-13 by default if u[1, i_large, j_large, element] >= equations.threshold_desingularization # 1e-10 by default u_buffer[1, i] = (u[1, i_large, j_large, element] + u[4, i_large, j_large, element]) From 2485b1a772a3eb11443cd05052b980fa95ed76d3 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Mon, 19 May 2025 19:57:17 +0200 Subject: [PATCH 24/37] adjust Monai elixir and do not allow 0 water heights on the mortars --- ...shallowwater_multilayer_monai_flood_amr.jl | 33 +++---------------- src/solvers/dgsem_p4est/dg_2d.jl | 16 ++++++--- test/test_p4est_2d.jl | 16 ++++----- 3 files changed, 24 insertions(+), 41 deletions(-) diff --git a/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_monai_flood_amr.jl b/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_monai_flood_amr.jl index 5f089a4d..89e46233 100644 --- a/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_monai_flood_amr.jl +++ b/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_monai_flood_amr.jl @@ -115,7 +115,7 @@ boundary_condition = Dict(:Bottom => boundary_condition_slip_wall, h, hv_1, hv_2, _ = u # Desingularization - h = (h^2 + max(h^2, 1e-8)) / (2.0 * h) + h = (h^2 + max(h^2, 1e-8)) / (2 * h) # friction coefficient n = 0.001 @@ -123,10 +123,10 @@ boundary_condition = Dict(:Bottom => boundary_condition_slip_wall, # Compute the common friction term Sf = -equations.gravity * n^2 * h^(-7 / 3) * sqrt(hv_1^2 + hv_2^2) - return SVector(zero(h)..., + return SVector(zero(h), Sf * hv_1, Sf * hv_2, - zero(b)) + zero(h)) end ############################################################################### @@ -194,32 +194,12 @@ save_solution = SaveSolutionCallback(dt = 0.2, save_initial_solution = true, save_final_solution = true) -# TimeSeries Callback is currently unavailable on `P4estMesh` -# time_series = TimeSeriesCallback(semi, [(4.521, 1.196), # G5 -# (4.521, 1.696), # G7 -# (4.521, 2.196)]; # G9 -# interval = 2, -# solution_variables=cons2prim) - stepsize_callback = StepsizeCallback(cfl = 0.5) -# Another possible AMR indicator function could be the velocity, such that it only fires -# in regions where the water is moving +# Another possible AMR indicator function could be the water height with IndicatorLoehner @inline function momentum_norm(u, equations::ShallowWaterMultiLayerEquations2D) _, h_v1, h_v2, _ = u - return sqrt(h_v1[1]^2 + h_v2[1]^2) - # h, h_v1, h_v2, _ = u - # v1 = (2 * h[1] * h_v1[1]) / (h[1]^2 + max(h[1]^2, 1e-8)) - # v2 = (2 * h[1] * h_v2[1]) / (h[1]^2 + max(h[1]^2, 1e-8)) - # if h[1] <= 1e-8 - # v1 = 0.0 - # v2 = 0.0 - # end - # if v1 > 0.525 || v2 > 0.525 - # v1 = 0.0 - # v2 = 0.0 - # end - # return sqrt(v1^2 + v2^2) + return sqrt(h_v1^2 + h_v2^2) end amr_controller = ControllerThreeLevel(semi, @@ -227,9 +207,6 @@ amr_controller = ControllerThreeLevel(semi, base_level = 0, med_level = 1, med_threshold = 0.012, max_level = 3, max_threshold = 0.015) -# velocity is very sensitive to use as an indicator -# med_level = 1, med_threshold = 0.25, -# max_level = 2, max_threshold = 0.475) amr_callback = AMRCallback(semi, amr_controller, interval = 5, diff --git a/src/solvers/dgsem_p4est/dg_2d.jl b/src/solvers/dgsem_p4est/dg_2d.jl index 4fb4ea24..95ce782a 100644 --- a/src/solvers/dgsem_p4est/dg_2d.jl +++ b/src/solvers/dgsem_p4est/dg_2d.jl @@ -214,7 +214,7 @@ function Trixi.prolong2mortars!(cache, u, if nlayers(equations) != 1 error("Non-conforming meshes are only supported for a single layer!") end - + Trixi.@threaded for mortar in Trixi.eachmortar(dg, cache) # Copy solution data from the small elements using "delayed indexing" with # a start value and a step size to get the correct face and orientation. @@ -301,10 +301,16 @@ function Trixi.prolong2mortars!(cache, u, # Note, no desingularization or water height cutoff is needed here as these steps are done # later in the hydrostatic reconstruction and stage limiter, respectively. for i in eachnode(dg) - cache.mortars.u[2, 1, 1, i, mortar] = (cache.mortars.u[2, 1, 1, i, mortar] - - cache.mortars.u[2, 4, 1, i, mortar]) - cache.mortars.u[2, 1, 2, i, mortar] = (cache.mortars.u[2, 1, 2, i, mortar] - - cache.mortars.u[2, 4, 2, i, mortar]) + cache.mortars.u[2, 1, 1, i, mortar] = max(cache.mortars.u[2, 1, 1, i, + mortar] - + cache.mortars.u[2, 4, 1, i, + mortar], + equations.threshold_limiter) + cache.mortars.u[2, 1, 2, i, mortar] = max(cache.mortars.u[2, 1, 2, i, + mortar] - + cache.mortars.u[2, 4, 2, i, + mortar], + equations.threshold_limiter) end end diff --git a/test/test_p4est_2d.jl b/test/test_p4est_2d.jl index 16b9d83c..71723627 100644 --- a/test/test_p4est_2d.jl +++ b/test/test_p4est_2d.jl @@ -189,15 +189,15 @@ end # SWE @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_multilayer_monai_flood_amr.jl"), l2=[ - 0.0003068569689525136, - 3.869268640941544e-6, - 6.345700409103784e-11, - 0.0003087516514830762], + 0.0003068569689525133, + 3.869268640941133e-6, + 6.345700395163902e-11, + 0.00030875165148307556], linf=[ - 0.004387570753720829, - 3.562731640300165e-5, - 1.190336978580838e-9, - 0.004387569562413193], + 0.004387570753720843, + 3.562731640298174e-5, + 1.1903369989050032e-9, + 0.004387569562413221], tspan=(0.0, 0.25)) # Ensure that we do not have excessive memory allocations From fb8cfbc562f6a80c7caddb6f2c46acea9c0f1b70 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Mon, 19 May 2025 20:00:50 +0200 Subject: [PATCH 25/37] keep the save solution callback in new WB elixir --- ...allowwater_multilayer_well_balanced_wet_dry_nonconforming.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_well_balanced_wet_dry_nonconforming.jl b/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_well_balanced_wet_dry_nonconforming.jl index 0407a660..e69dd8fa 100644 --- a/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_well_balanced_wet_dry_nonconforming.jl +++ b/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_well_balanced_wet_dry_nonconforming.jl @@ -198,7 +198,7 @@ stepsize_callback = StepsizeCallback(cfl = 1.0) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, - # save_solution, + save_solution, stepsize_callback) stage_limiter! = PositivityPreservingLimiterShallowWater(variables = (waterheight,)) From c7b6b3879ec6347eab850f1ab653c2a1842d5b70 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Mon, 19 May 2025 20:23:59 +0200 Subject: [PATCH 26/37] fix merge conflicts in Monai elixir --- ...r_shallowwater_multilayer_three_mound_dam_break_amr.jl | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_three_mound_dam_break_amr.jl b/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_three_mound_dam_break_amr.jl index def8891d..d62d5e87 100644 --- a/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_three_mound_dam_break_amr.jl +++ b/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_three_mound_dam_break_amr.jl @@ -130,7 +130,7 @@ analysis_callback = AnalysisCallback(semi, interval = analysis_interval, alive_callback = AliveCallback(analysis_interval = analysis_interval) -save_solution = SaveSolutionCallback(dt = 0.5, +save_solution = SaveSolutionCallback(dt = 0.04, #dt = 0.5, save_initial_solution = true, save_final_solution = true) @@ -147,8 +147,10 @@ amr_indicator = IndicatorLöhner(semi, variable = main_waterheight) amr_controller = ControllerThreeLevel(semi, amr_indicator, base_level = 0, - med_level = 1, med_threshold = 0.25, - max_level = 3, max_threshold = 0.5) + # med_level = 1, med_threshold = 0.25, + # max_level = 3, max_threshold = 0.5) + med_level = 1, med_threshold = 0.1, + max_level = 3, max_threshold = 0.25) amr_callback = AMRCallback(semi, amr_controller, interval = 5, From ecc3f61bdef6386cca2623495b892a76ee3bf30a Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Tue, 20 May 2025 15:59:00 +0200 Subject: [PATCH 27/37] adjust AMR parameters in three mound test case --- .../elixir_shallowwater_monai_tsunami.jl | 2 +- ...r_shallowwater_multilayer_monai_flood_amr.jl | 2 +- ...ater_multilayer_three_mound_dam_break_amr.jl | 7 ++----- test/test_p4est_2d.jl | 17 +++++++++-------- 4 files changed, 13 insertions(+), 15 deletions(-) diff --git a/docs/src/tutorials/elixir_shallowwater_monai_tsunami.jl b/docs/src/tutorials/elixir_shallowwater_monai_tsunami.jl index f4087b79..4063a5bf 100644 --- a/docs/src/tutorials/elixir_shallowwater_monai_tsunami.jl +++ b/docs/src/tutorials/elixir_shallowwater_monai_tsunami.jl @@ -266,7 +266,7 @@ boundary_condition = Dict(:Bottom => boundary_condition_slip_wall, h, hv_1, hv_2, _ = u n = 0.001 # friction coefficient - h = (h^2 + max(h^2, 1e-8)) / (2.0 * h) # desingularization procedure + h = (h^2 + max(h^2, 1e-8)) / (2 * h) # desingularization procedure ## Compute the common friction term Sf = -equations.gravity * n^2 * h^(-7 / 3) * sqrt(hv_1^2 + hv_2^2) diff --git a/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_monai_flood_amr.jl b/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_monai_flood_amr.jl index 89e46233..0280c0ee 100644 --- a/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_monai_flood_amr.jl +++ b/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_monai_flood_amr.jl @@ -196,7 +196,7 @@ save_solution = SaveSolutionCallback(dt = 0.2, stepsize_callback = StepsizeCallback(cfl = 0.5) -# Another possible AMR indicator function could be the water height with IndicatorLoehner +# Another possible AMR indicator function could be the water height with `IndicatorLoehner` @inline function momentum_norm(u, equations::ShallowWaterMultiLayerEquations2D) _, h_v1, h_v2, _ = u return sqrt(h_v1^2 + h_v2^2) diff --git a/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_three_mound_dam_break_amr.jl b/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_three_mound_dam_break_amr.jl index d62d5e87..6bbea7ca 100644 --- a/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_three_mound_dam_break_amr.jl +++ b/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_three_mound_dam_break_amr.jl @@ -130,7 +130,7 @@ analysis_callback = AnalysisCallback(semi, interval = analysis_interval, alive_callback = AliveCallback(analysis_interval = analysis_interval) -save_solution = SaveSolutionCallback(dt = 0.04, #dt = 0.5, +save_solution = SaveSolutionCallback(dt = 0.5, save_initial_solution = true, save_final_solution = true) @@ -138,17 +138,14 @@ stepsize_callback = StepsizeCallback(cfl = 0.4) # Cannot simply use `waterheight` here for multilayer equations. # Need a helper function to extract the relevant variable. -# TODO: This AMR indicator fires too often in the dry regions @inline function main_waterheight(u, equations) return waterheight(u, equations)[1] end -amr_indicator = IndicatorLöhner(semi, variable = main_waterheight) +amr_indicator = IndicatorLoehner(semi, variable = main_waterheight) amr_controller = ControllerThreeLevel(semi, amr_indicator, base_level = 0, - # med_level = 1, med_threshold = 0.25, - # max_level = 3, max_threshold = 0.5) med_level = 1, med_threshold = 0.1, max_level = 3, max_threshold = 0.25) diff --git a/test/test_p4est_2d.jl b/test/test_p4est_2d.jl index 71723627..f33cf631 100644 --- a/test/test_p4est_2d.jl +++ b/test/test_p4est_2d.jl @@ -214,17 +214,18 @@ end # SWE @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_multilayer_three_mound_dam_break_amr.jl"), l2=[ - 0.13142461271166225, - 0.35364621879644254, - 4.778520542281107e-13, - 0.0019606480801664917], + 0.1307018038119991, + 0.35570420514229945, + 9.300123179370889e-14, + 0.0019606480801663984], linf=[ - 1.1338438044790615, - 2.416846787140322, - 6.256251416535554e-11, + 1.1241023021698624, + 2.4980720394823885, + 1.324711724284584e-11, 0.044079775502972124], tspan=(0.0, 0.3), - coverage_override=(maxiters = 5,)) + coverage_override=(maxiters = 5,), + atol=1e-10) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) From f0384fbf581cc954e8124335efc8066763c761a6 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Tue, 20 May 2025 20:16:43 +0200 Subject: [PATCH 28/37] adjust comments in new mortar routine --- src/solvers/dgsem_p4est/dg_2d.jl | 36 ++++++++++++++------------------ 1 file changed, 16 insertions(+), 20 deletions(-) diff --git a/src/solvers/dgsem_p4est/dg_2d.jl b/src/solvers/dgsem_p4est/dg_2d.jl index 95ce782a..f9bed8d2 100644 --- a/src/solvers/dgsem_p4est/dg_2d.jl +++ b/src/solvers/dgsem_p4est/dg_2d.jl @@ -5,11 +5,6 @@ @muladd begin #! format: noindent -# TODO: Once working the mortar methods could likely be extended to the other -# equation types available in the package. Although for the multilayer equations -# care must be taken because the pressure term is separated from the physical flux -# and directly placed in the nonconservative flux - # The methods below are specialized on the `P4estShallowWaterMortarContainer` # mortar type. Extra storage is necessary for the numerical flux plus nonconservative term # on either side of the parent (large) element mortars. It also requires a work @@ -194,8 +189,8 @@ end # Specialized strategy to project the solution to the mortars and preserve # well-balancedness from the work of Benov et al. (https://doi.org/10.1016/j.jcp.2018.02.008). -# This version is for the `ShallowWaterMultiLayerEquations2D` with a single layer, which "peels" the pressure -# contribution into the nonconservative term for easier well-balancing, +# This version is for the `ShallowWaterMultiLayerEquations2D` with a single layer +# that "peels" the pressure contribution into the nonconservative term for easier well-balancing, # see Ersing et al. (https://doi.org/10.1016/j.jcp.2025.113802). # Thus, this uses the original `P4estMortarContainer`, `calc_mortar_flux!` # and `mortar_fluxes_to_elements!` from Trixi.jl. @@ -263,9 +258,9 @@ function Trixi.prolong2mortars!(cache, u, # That is, we might be able to directly project the water height `h` instead while maintaining # important steady-state solution behavior. # - # Further, worth noting that for fully wet portions of the domain, the logic below + # Further, it is worth noting that for fully wet portions of the domain, the logic below # is unnecessary and the "classic" mortar method is well-balanced and fully conservative - # for the `ShallowWaterMultiLayerEquations2D`. Could be exploited to design a better + # for the `ShallowWaterMultiLayerEquations2D`. This could be exploited to design a better # approach to the AMR refinement / coarsening. # # Note, some shift is required to ensure we catch water heights close to the threshold @@ -298,19 +293,20 @@ function Trixi.prolong2mortars!(cache, u, # and instead recover the conservative water height variable `h`. # Basically, unpacking the total water height variable to create the projected local water # height `h` from Eq. 41 in Benov et al. - # Note, no desingularization or water height cutoff is needed here as these steps are done + # Note, no desingularization or water height cutoff is done here as these steps are performed # later in the hydrostatic reconstruction and stage limiter, respectively. + # + # In testing, spurious waves are generated if one instead uses + # `max(H - b, equations.threshold_limiter)` on the mortars and using a velocity + # desingularization on the mortars here causes unpredictable behavior in + # the `LoehnerIndicator` for AMR in the dry regions. This is likely due to the gradient + # estimation in the dry "noisy" regions being sensitive to common indicator variables + # like the local water height. for i in eachnode(dg) - cache.mortars.u[2, 1, 1, i, mortar] = max(cache.mortars.u[2, 1, 1, i, - mortar] - - cache.mortars.u[2, 4, 1, i, - mortar], - equations.threshold_limiter) - cache.mortars.u[2, 1, 2, i, mortar] = max(cache.mortars.u[2, 1, 2, i, - mortar] - - cache.mortars.u[2, 4, 2, i, - mortar], - equations.threshold_limiter) + cache.mortars.u[2, 1, 1, i, mortar] = (cache.mortars.u[2, 1, 1, i, mortar] - + cache.mortars.u[2, 4, 1, i, mortar]) + cache.mortars.u[2, 1, 2, i, mortar] = (cache.mortars.u[2, 1, 2, i, mortar] - + cache.mortars.u[2, 4, 2, i, mortar]) end end From e47cef610e15198ecfa614a8424be73ab18a50dd Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Tue, 20 May 2025 20:18:43 +0200 Subject: [PATCH 29/37] apply formatter --- src/solvers/dgsem_p4est/dg_2d.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/solvers/dgsem_p4est/dg_2d.jl b/src/solvers/dgsem_p4est/dg_2d.jl index f9bed8d2..de9d2d51 100644 --- a/src/solvers/dgsem_p4est/dg_2d.jl +++ b/src/solvers/dgsem_p4est/dg_2d.jl @@ -304,9 +304,9 @@ function Trixi.prolong2mortars!(cache, u, # like the local water height. for i in eachnode(dg) cache.mortars.u[2, 1, 1, i, mortar] = (cache.mortars.u[2, 1, 1, i, mortar] - - cache.mortars.u[2, 4, 1, i, mortar]) + cache.mortars.u[2, 4, 1, i, mortar]) cache.mortars.u[2, 1, 2, i, mortar] = (cache.mortars.u[2, 1, 2, i, mortar] - - cache.mortars.u[2, 4, 2, i, mortar]) + cache.mortars.u[2, 4, 2, i, mortar]) end end From 707b049b25a4c1cf3582457c91a91e0d6d7652bf Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Tue, 20 May 2025 20:28:40 +0200 Subject: [PATCH 30/37] adjust test values --- test/test_p4est_2d.jl | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/test/test_p4est_2d.jl b/test/test_p4est_2d.jl index f33cf631..10ee3103 100644 --- a/test/test_p4est_2d.jl +++ b/test/test_p4est_2d.jl @@ -141,7 +141,8 @@ end # SWE 1.0894983713455818e-14, 1.132829841145215e-14, 0.41600528180178625], - tspan=(0.0, 0.25)) + tspan=(0.0, 0.25), + atol=1e-10) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) let @@ -189,14 +190,14 @@ end # SWE @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_multilayer_monai_flood_amr.jl"), l2=[ - 0.0003068569689525133, - 3.869268640941133e-6, - 6.345700395163902e-11, + 0.00030685696894658095, + 3.869268491914826e-6, + 6.345699712019463e-11, 0.00030875165148307556], linf=[ - 0.004387570753720843, - 3.562731640298174e-5, - 1.1903369989050032e-9, + 0.0043875707537211345, + 3.562731345127931e-5, + 1.1903363999889186e-9, 0.004387569562413221], tspan=(0.0, 0.25)) From 1209e1d4fabbdfd708e0d91e390b693da2e88004 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Tue, 20 May 2025 21:00:08 +0200 Subject: [PATCH 31/37] update header comments in the four new elixirs to summarize the current state of affairs --- .../elixir_shallowwater_multilayer_monai_flood_amr.jl | 7 +++++++ ...shallowwater_multilayer_perturbation_wet_dry_amr.jl | 7 +++++++ ...hallowwater_multilayer_three_mound_dam_break_amr.jl | 9 +++++++++ ...r_multilayer_well_balanced_wet_dry_nonconforming.jl | 10 +++++++++- 4 files changed, 32 insertions(+), 1 deletion(-) diff --git a/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_monai_flood_amr.jl b/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_monai_flood_amr.jl index 0280c0ee..b536c102 100644 --- a/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_monai_flood_amr.jl +++ b/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_monai_flood_amr.jl @@ -9,6 +9,13 @@ using TrixiBottomTopography # for a thorough explanation of this problem setup. # # This elixir allows for experimentation with AMR with this complex test case. +# In the current state, AMR on wet/dry fronts remains sensitive to the limiting of the water height, +# desingularization of the velocities, and the design of the AMR indicator. +# Further, there remains an appreciable loss of conservation, with conservation errors +# of ≈1e-5 in the water height for this test case as the flow evolves. +# More investigation is needed to identify exactly why each aspect of this strange behavior occurs. +# For instance, the conservation loss may be related to the choice of the CFL that guarantees +# a positive water height as well as a positive element-wise mean necessary in the limiting. ############################################################################### # Semidiscretization of the multilayer shallow water equations with one layer diff --git a/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_perturbation_wet_dry_amr.jl b/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_perturbation_wet_dry_amr.jl index 9545f474..1b9893ca 100644 --- a/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_perturbation_wet_dry_amr.jl +++ b/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_perturbation_wet_dry_amr.jl @@ -4,6 +4,13 @@ using Trixi using TrixiShallowWater # Elixir useful for stress testing the AMR + wet/dry + well-balanced mortars. +# In the current state, AMR on wet/dry fronts remains sensitive to the limiting of the water height, +# desingularization of the velocities, and the design of the AMR indicator. +# Further, there remains an appreciable loss of conservation, with conservation errors +# of ≈1e-5 in the water height for this test case as the flow evolves. +# More investigation is needed to identify exactly why each aspect of this strange behavior occurs. +# For instance, the conservation loss may be related to the choice of the CFL that guarantees +# a positive water height as well as a positive element-wise mean necessary in the limiting. ############################################################################### # Semidiscretization of the multilayer shallow water equations with one layer diff --git a/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_three_mound_dam_break_amr.jl b/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_three_mound_dam_break_amr.jl index 6bbea7ca..815ec3b2 100644 --- a/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_three_mound_dam_break_amr.jl +++ b/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_three_mound_dam_break_amr.jl @@ -3,6 +3,15 @@ using OrdinaryDiffEqSSPRK, OrdinaryDiffEqLowStorageRK using Trixi using TrixiShallowWater +# This elixir allows for experimentation with AMR with this complex test case. +# In the current state, AMR on wet/dry fronts remains sensitive to the limiting of the water height, +# desingularization of the velocities, and the design of the AMR indicator. +# Further, there remains an appreciable loss of conservation, with conservation errors +# of ≈1e-5 in the water height for this test case as the flow evolves. +# More investigation is needed to identify exactly why each aspect of this strange behavior occurs. +# For instance, the conservation loss may be related to the choice of the CFL that guarantees +# a positive water height as well as a positive element-wise mean necessary in the limiting. + ############################################################################### # Semidiscretization of the multilayer shallow water equations with one layer # to fallback to be the standard shallow water equations diff --git a/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_well_balanced_wet_dry_nonconforming.jl b/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_well_balanced_wet_dry_nonconforming.jl index e69dd8fa..990b39e9 100644 --- a/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_well_balanced_wet_dry_nonconforming.jl +++ b/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_well_balanced_wet_dry_nonconforming.jl @@ -3,10 +3,18 @@ using OrdinaryDiffEqSSPRK, OrdinaryDiffEqLowStorageRK using Trixi using TrixiShallowWater +# Stress test case for well-balancedness with wet/dry transition elements +# on a nonconforming mesh with a discontinuous bottom topography where +# this elixir allows for experimentation for well-balancedness with this complex test case. + +# In the current state, recovering a well-balancedness error up to double precision machine +# round-off error relies on a special projection procedure at wet/dry fronts. +# However, there are no conservation issues and these conservation errors remain around +# double precision unit round-off for all time. + ############################################################################### # Semidiscretization of the multilayer shallow water equations with one layer # to fallback to the standard shallow water equations. -# Test case for well-balancedness with wet/dry transition elements on a nonconforming mesh. equations = ShallowWaterMultiLayerEquations2D(gravity = 9.812, H0 = 1.235, rhos = 1.0) From 21bd02bb431e434e7dbc494d63b280cf4b8f7c46 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Tue, 20 May 2025 21:05:57 +0200 Subject: [PATCH 32/37] apply formatter to examples --- .../elixir_shallowwater_multilayer_monai_flood_amr.jl | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_monai_flood_amr.jl b/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_monai_flood_amr.jl index b536c102..08c06518 100644 --- a/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_monai_flood_amr.jl +++ b/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_monai_flood_amr.jl @@ -34,7 +34,8 @@ spline_bathymetry_file = Trixi.download("https://gist.githubusercontent.com/andr # B-spline interpolation of the underlying data. # The type of this struct is fixed as `BicubicBSpline`. -const bath_spline_struct = BicubicBSpline(spline_bathymetry_file, end_condition = "not-a-knot") +const bath_spline_struct = BicubicBSpline(spline_bathymetry_file, + end_condition = "not-a-knot") bathymetry(x::RealT, y::RealT) = spline_interpolation(bath_spline_struct, x, y) # Initial condition is a constant background state @@ -81,8 +82,8 @@ const h_spline_struct = CubicBSpline(water_height_data; end_condition = "not-a-k H_from_wave_maker(t::RealT) = spline_interpolation(h_spline_struct, t) # Now we can define the specialized boundary condition for the incident wave maker. -@inline function boundary_condition_wave_maker(u_inner, normal_direction::AbstractVector, x, t, - surface_flux_functions, +@inline function boundary_condition_wave_maker(u_inner, normal_direction::AbstractVector, + x, t, surface_flux_functions, equations::ShallowWaterMultiLayerEquations2D) surface_flux_function, nonconservative_flux_function = surface_flux_functions From 92a196440959f961918ccd59b04b01d5be4b6bb7 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Tue, 20 May 2025 21:12:09 +0200 Subject: [PATCH 33/37] adjust sensitive perturbation test values --- test/test_p4est_2d.jl | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/test/test_p4est_2d.jl b/test/test_p4est_2d.jl index 10ee3103..238468ba 100644 --- a/test/test_p4est_2d.jl +++ b/test/test_p4est_2d.jl @@ -160,16 +160,16 @@ end # SWE @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_multilayer_perturbation_wet_dry_amr.jl"), l2=[ - 0.3997030663722126, - 0.38644026433480094, - 0.4362923911447016, - 0.4564350570616195], + 0.4177955194235321, + 0.2818773677696145, + 0.33178075151092923, + 0.45643505879802193], linf=[ - 1.372326981328246, - 3.3262654799249893, - 3.6541077514853653, + 2.129878820442652, + 3.292431218004944, + 3.3701108685790135, 0.7495177590247986], - tspan=(0.0, 0.025), + tspan=(0.0, 0.01), coverage_override=(maxiters = 5,)) # Ensure that we do not have excessive memory allocations From b54b742542366f01499ae1c3ea823e7fe4b2d3c2 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Tue, 20 May 2025 21:20:21 +0200 Subject: [PATCH 34/37] apply formatter --- .../elixir_shallowwater_multilayer_monai_flood_amr.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_monai_flood_amr.jl b/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_monai_flood_amr.jl index 08c06518..e4131209 100644 --- a/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_monai_flood_amr.jl +++ b/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_monai_flood_amr.jl @@ -78,7 +78,8 @@ water_height_data = Trixi.download("https://gist.githubusercontent.com/andrewwin # Similar to the bathymetry approximation, we construct a cubic B-spline interpolation # of the data, then create a function to evaluate the resulting spline at a given $t$ value. # The type of this struct is fixed as `CubicBSpline`. -const h_spline_struct = CubicBSpline(water_height_data; end_condition = "not-a-knot") +const h_spline_struct = CubicBSpline(water_height_data; + end_condition = "not-a-knot") H_from_wave_maker(t::RealT) = spline_interpolation(h_spline_struct, t) # Now we can define the specialized boundary condition for the incident wave maker. From ce125304ad8dcdf3a019423c801b0172d561a60f Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Tue, 20 May 2025 22:15:42 +0200 Subject: [PATCH 35/37] apply formatter --- docs/src/tutorials/elixir_shallowwater_monai_tsunami.jl | 3 ++- .../elixir_shallowwater_multilayer_monai_flood_amr.jl | 3 +-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/src/tutorials/elixir_shallowwater_monai_tsunami.jl b/docs/src/tutorials/elixir_shallowwater_monai_tsunami.jl index 4063a5bf..a803a680 100644 --- a/docs/src/tutorials/elixir_shallowwater_monai_tsunami.jl +++ b/docs/src/tutorials/elixir_shallowwater_monai_tsunami.jl @@ -170,7 +170,8 @@ spline_bathymetry_file = Trixi.download("https://gist.githubusercontent.com/andr # Create a bicubic B-spline interpolation of the bathymetry data, then create a function # to evaluate the resulting spline at a given point $(x,y)$. # The type of this struct is fixed as `BicubicBSpline`. -const bath_spline_struct = BicubicBSpline(spline_bathymetry_file, end_condition = "not-a-knot") +const bath_spline_struct = BicubicBSpline(spline_bathymetry_file, + end_condition = "not-a-knot") bathymetry(x::Float64, y::Float64) = spline_interpolation(bath_spline_struct, x, y); # We then create a function to supply the initial condition for the simulation. diff --git a/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_monai_flood_amr.jl b/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_monai_flood_amr.jl index e4131209..08c06518 100644 --- a/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_monai_flood_amr.jl +++ b/examples/p4est_2d_dgsem/elixir_shallowwater_multilayer_monai_flood_amr.jl @@ -78,8 +78,7 @@ water_height_data = Trixi.download("https://gist.githubusercontent.com/andrewwin # Similar to the bathymetry approximation, we construct a cubic B-spline interpolation # of the data, then create a function to evaluate the resulting spline at a given $t$ value. # The type of this struct is fixed as `CubicBSpline`. -const h_spline_struct = CubicBSpline(water_height_data; - end_condition = "not-a-knot") +const h_spline_struct = CubicBSpline(water_height_data; end_condition = "not-a-knot") H_from_wave_maker(t::RealT) = spline_interpolation(h_spline_struct, t) # Now we can define the specialized boundary condition for the incident wave maker. From e06d7ff751413dbc793bc787271fc80e139527dc Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Tue, 20 May 2025 22:19:09 +0200 Subject: [PATCH 36/37] slightly soften abs tol for Monai test to 1e-12 --- test/test_p4est_2d.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/test/test_p4est_2d.jl b/test/test_p4est_2d.jl index 238468ba..40fbad01 100644 --- a/test/test_p4est_2d.jl +++ b/test/test_p4est_2d.jl @@ -199,7 +199,8 @@ end # SWE 3.562731345127931e-5, 1.1903363999889186e-9, 0.004387569562413221], - 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) From b94b73d291e66222a4e9c907dd1efd391fc42f1b Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Tue, 20 May 2025 22:35:25 +0200 Subject: [PATCH 37/37] further soften tolenace for Monai --- test/test_p4est_2d.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_p4est_2d.jl b/test/test_p4est_2d.jl index 40fbad01..14a68a3b 100644 --- a/test/test_p4est_2d.jl +++ b/test/test_p4est_2d.jl @@ -200,7 +200,7 @@ end # SWE 1.1903363999889186e-9, 0.004387569562413221], tspan=(0.0, 0.25), - atol=1e-12) + atol=5e-12) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities)