From fecddd2f02c7d8a8ebe9f3323fe4975932c6ffe6 Mon Sep 17 00:00:00 2001 From: Andrew Winters Date: Sun, 26 May 2024 00:10:57 +0200 Subject: [PATCH] wet-dry well balancedness works but perturbation test on wet-drt still crashes --- .../elixir_shallowwater_perturbation.jl | 122 ++++++++++ ...lixir_shallowwater_perturbation_wet_dry.jl | 211 ++++++++++++++++++ .../elixir_shallowwater_well_balanced.jl | 6 +- ...ixir_shallowwater_well_balanced_wet_dry.jl | 52 +++-- src/equations/shallow_water_wet_dry_2d.jl | 2 + src/solvers/scratch_p4est.jl | 39 +++- 6 files changed, 393 insertions(+), 39 deletions(-) create mode 100644 examples/p4est_2d_dgsem/elixir_shallowwater_perturbation.jl create mode 100644 examples/p4est_2d_dgsem/elixir_shallowwater_perturbation_wet_dry.jl diff --git a/examples/p4est_2d_dgsem/elixir_shallowwater_perturbation.jl b/examples/p4est_2d_dgsem/elixir_shallowwater_perturbation.jl new file mode 100644 index 0000000..f6401f3 --- /dev/null +++ b/examples/p4est_2d_dgsem/elixir_shallowwater_perturbation.jl @@ -0,0 +1,122 @@ + +using Downloads: download +using OrdinaryDiffEq +using Trixi +using TrixiShallowWater + +############################################################################### +# semidiscretization of the shallow water equations with a continuous +# bottom topography function + +equations = ShallowWaterEquationsWetDry2D(gravity_constant = 9.812, H0 = 2.1) + +function initial_condition_perturbation(x, t, equations::ShallowWaterEquationsWetDry2D) + # Calculate primitive variables + H = equations.H0 + v1 = 0.0 + v2 = 0.0 + + x1, x2 = x + b = (1.75 / 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))) + + # Waterheight perturbation + H = H + 0.5 * exp(-40.0 * ((x[1])^2 + (x[2])^2)) + + 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_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) + +surface_flux = (FluxHydrostaticReconstruction(flux_hll_chen_noelle, hydrostatic_reconstruction_chen_noelle), + flux_nonconservative_chen_noelle) + +# Create the solver +solver = DGSEM(polydeg = 3, surface_flux = surface_flux, + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +############################################################################### +# This setup is for the curved, split form well-balancedness testing + +# 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.2 is even more extreme. +function mapping_twist(xi, eta) + y = eta + 0.175 * cos(1.5 * pi * xi) * cos(0.5 * pi * eta) + x = xi + 0.175 * 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) + +# Create the semi discretization object +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_condition) + +############################################################################### +# ODE solvers, callbacks, etc. + +tspan = (0.0, 0.5) +ode = semidiscretize(semi, tspan) + +############################################################################### + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(dt = 0.04, + save_initial_solution = true, + save_final_solution = true) + +# Define a better variable to use in the AMR indicator +@inline function total_water_height(u, equations::ShallowWaterEquationsWetDry2D) + return u[1] + u[4] +end + +amr_controller = ControllerThreeLevel(semi, IndicatorMax(semi, variable = total_water_height), + base_level = 0, + med_level = 1, med_threshold = 2.01, + max_level = 4, max_threshold = 2.15) + +amr_callback = AMRCallback(semi, amr_controller, + interval = 1, + adapt_initial_condition = false, # to setup discontinuous bottom easier + adapt_initial_condition_only_refine = false) + +stepsize_callback = StepsizeCallback(cfl = 1.0) + +callbacks = CallbackSet(summary_callback, + analysis_callback, + alive_callback, + save_solution, + amr_callback, + stepsize_callback) + +############################################################################### +# run the simulation +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + # adaptive = false, + save_everystep = false, callback = callbacks); +summary_callback() # print the timer summary diff --git a/examples/p4est_2d_dgsem/elixir_shallowwater_perturbation_wet_dry.jl b/examples/p4est_2d_dgsem/elixir_shallowwater_perturbation_wet_dry.jl new file mode 100644 index 0000000..764ffce --- /dev/null +++ b/examples/p4est_2d_dgsem/elixir_shallowwater_perturbation_wet_dry.jl @@ -0,0 +1,211 @@ + +using Downloads: download +using OrdinaryDiffEq +using Trixi +using TrixiShallowWater + +# TODO: This elixir is for debugging purposes of the AMR + well-balanced mortars. +# Currently this crashes very early in the simulation. ONe possible explanation is that +# the projection / interpolation of the solution as it is refined / coarsened may need +# to use a similar strategy as the mortar projection where we ensure that a constant solution +# is projected on the dry elements. Further investigation is needed to see if this is indeed +# the case. + +############################################################################### +# semidiscretization of the shallow water equations with a continuous +# bottom topography function + +equations = ShallowWaterEquationsWetDry2D(gravity_constant = 9.812, H0 = 1.235) + +function initial_condition_perturbation(x, t, equations::ShallowWaterEquationsWetDry2D) + # Calculate primitive variables + H = equations.H0 + v1 = 0.0 + v2 = 0.0 + + x1, x2 = x + b = (1.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))) + + # for testing + # b = zero(eltype(x)) + + H = max(equations.H0, b + equations.threshold_limiter) + + 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_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) + +surface_flux = (FluxHydrostaticReconstruction(flux_hll_chen_noelle, hydrostatic_reconstruction_chen_noelle), + flux_nonconservative_chen_noelle) + +# # Create the solver (for fully wet!) +# solver = DGSEM(polydeg = 3, surface_flux = surface_flux, +# volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +# Create the solver +basis = LobattoLegendreBasis(3) + +indicator_sc = IndicatorHennemannGassnerShallowWater(equations, basis, + alpha_max = 0.5, + alpha_min = 0.001, + alpha_smooth = true, + variable = Trixi.waterheight) # waterheight_pressure) +volume_integral = VolumeIntegralShockCapturingHG(indicator_sc; + volume_flux_dg = volume_flux, + volume_flux_fv = surface_flux) + +solver = DGSEM(basis, surface_flux, volume_integral) + +############################################################################### +# This setup is for the curved, split form well-balancedness testing + +# 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.2 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 2 +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, 0.5) +ode = semidiscretize(semi, tspan) + +function initial_condition_discontinuous_well_balancedness(x, t, element_id, equations::ShallowWaterEquationsWetDry2D) + # Set the background values for velocity + 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 + if element_id > 207 && element_id < 300 # for the forced hanging node grid with with < 3 in function above + 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 discontinous purturbation using the element number + if element_id in [232, 224, 225, 226, 227, 228, 229, 230] + H = H + 0.3 + end + + # Clip the initialization to avoid negative water heights and division by zero + h = max(equations.threshold_limiter, H - b) + + # Return the conservative variables + return SVector(h, h * v1, h * v2, b) +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 = 1000 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval=20, #dt = 0.04, + save_initial_solution = true, + save_final_solution = true) + +# # Define a better variable to use in the AMR indicator +# @inline function total_water_height(u, equations::ShallowWaterEquationsWetDry2D) +# return u[1] + u[4] +# end + +# amr_controller = ControllerThreeLevel(semi, IndicatorMax(semi, variable = total_water_height), +# base_level = 1, +# med_level = 2, med_threshold = 1.225, # 2.01, +# max_level = 4, max_threshold = 1.245) # 2.15) + +# amr_callback = AMRCallback(semi, amr_controller, +# interval = 1, +# adapt_initial_condition = false, +# adapt_initial_condition_only_refine = false) + +stepsize_callback = StepsizeCallback(cfl = 0.2) + +callbacks = CallbackSet(summary_callback, + analysis_callback, + alive_callback, + save_solution, + # amr_callback, + stepsize_callback) + +stage_limiter! = PositivityPreservingLimiterShallowWater(variables = (Trixi.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); +summary_callback() # print the timer summary diff --git a/examples/p4est_2d_dgsem/elixir_shallowwater_well_balanced.jl b/examples/p4est_2d_dgsem/elixir_shallowwater_well_balanced.jl index 9d7bb8c..8604f70 100644 --- a/examples/p4est_2d_dgsem/elixir_shallowwater_well_balanced.jl +++ b/examples/p4est_2d_dgsem/elixir_shallowwater_well_balanced.jl @@ -49,7 +49,7 @@ initial_condition = initial_condition_wb_testing # surface_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) # surface_flux = (FluxHydrostaticReconstruction(flux_hll_chen_noelle, hydrostatic_reconstruction_chen_noelle), -# flux_nonconservative_chen_noelle) +# flux_nonconservative_chen_noelle) # # Fjordholms flux for testing @@ -59,7 +59,7 @@ initial_condition = initial_condition_wb_testing volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) surface_flux = (FluxHydrostaticReconstruction(flux_lax_friedrichs, hydrostatic_reconstruction_audusse_etal), - flux_nonconservative_audusse_etal) + flux_nonconservative_audusse_etal) # Audusse HR with HLL flux # surface_flux = (FluxHydrostaticReconstruction(flux_hll, hydrostatic_reconstruction_audusse_etal), @@ -85,7 +85,7 @@ solver = DGSEM(polydeg = 3, surface_flux = surface_flux, boundary_condition = Dict(:all => boundary_condition_slip_wall) # Affine type mapping to take the [-1,1]^2 domain from the mesh file -# and warp it with a mapping as described in https://arxiv.org/abs/2012.12040 +# and warp it as described in https://arxiv.org/abs/2012.12040 # Warping with the coefficient 0.2 is even more extreme. function mapping_twist(xi, eta) y = eta + 0.15 * cos(1.5 * pi * xi) * cos(0.5 * pi * eta) diff --git a/examples/p4est_2d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl b/examples/p4est_2d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl index 68d2610..b5e733e 100644 --- a/examples/p4est_2d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl +++ b/examples/p4est_2d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl @@ -8,8 +8,7 @@ using TrixiShallowWater # semidiscretization of the shallow water equations with a continuous # bottom topography function -equations = ShallowWaterEquationsWetDry2D(gravity_constant = 9.812, H0 = 1.235, - threshold_partially_wet = 1e-4) +equations = ShallowWaterEquationsWetDry2D(gravity_constant = 9.812, H0 = 1.235) function initial_condition_wb_testing(x, t, equations::ShallowWaterEquationsWetDry2D) # Calculate primitive variables @@ -47,23 +46,22 @@ initial_condition = initial_condition_wb_testing # volume_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal) # surface_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal) -# # Wintermeyer flux in the volume -# volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) +# Wintermeyer flux in the volume +volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) # Wintermeyer flux in the surface for testing -surface_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) - -# surface_flux = (FluxHydrostaticReconstruction(flux_hll_chen_noelle, hydrostatic_reconstruction_chen_noelle), -# flux_nonconservative_chen_noelle) +# surface_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) +surface_flux = (FluxHydrostaticReconstruction(flux_hll_chen_noelle, hydrostatic_reconstruction_chen_noelle), + flux_nonconservative_chen_noelle) # # Fjordholms flux for testing # surface_flux = (flux_fjordholm_etal, flux_nonconservative_fjordholm_etal) # # Audusse HR with Rusanov flux -surface_flux = (FluxHydrostaticReconstruction(flux_lax_friedrichs, - hydrostatic_reconstruction_audusse_etal), - flux_nonconservative_audusse_etal) +# surface_flux = (FluxHydrostaticReconstruction(flux_lax_friedrichs, +# hydrostatic_reconstruction_audusse_etal), +# flux_nonconservative_audusse_etal) # Audusse HR with HLL flux # surface_flux = (FluxHydrostaticReconstruction(flux_hll, hydrostatic_reconstruction_audusse_etal), @@ -76,7 +74,7 @@ indicator_sc = IndicatorHennemannGassnerShallowWater(equations, basis, alpha_max = 0.5, alpha_min = 0.001, alpha_smooth = true, - variable = waterheight_pressure) + variable = Trixi.waterheight) # waterheight_pressure) volume_integral = VolumeIntegralShockCapturingHG(indicator_sc; volume_flux_dg = volume_flux, volume_flux_fv = surface_flux) @@ -99,7 +97,7 @@ solver = DGSEM(basis, surface_flux, volume_integral) boundary_condition = Dict(:all => boundary_condition_slip_wall) # Affine type mapping to take the [-1,1]^2 domain from the mesh file -# and warp it with a mapping as described in https://arxiv.org/abs/2012.12040 +# and warp it as described in https://arxiv.org/abs/2012.12040 # Warping with the coefficient 0.2 is even more extreme. function mapping_twist(xi, eta) y = eta + 0.15 * cos(1.5 * pi * xi) * cos(0.5 * pi * eta) @@ -118,7 +116,7 @@ mesh = P4estMesh{2}(mesh_file, polydeg = 3, # Refine bottom left quadrant of each tree to level 2 function refine_fn(p4est, which_tree, quadrant) quadrant_obj = unsafe_load(quadrant) - if quadrant_obj.x == 0 && quadrant_obj.y == 0 && quadrant_obj.level < 2 + if quadrant_obj.x == 0 && quadrant_obj.y == 0 && quadrant_obj.level < 3 # return true (refine) return Cint(1) else @@ -127,7 +125,7 @@ function refine_fn(p4est, which_tree, quadrant) end end -# Refine recursively until each bottom left quadrant of a tree has level 2 +# 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}, @@ -156,11 +154,12 @@ ode = semidiscretize(semi, tspan) # only for the specific mesh loaded above! function initial_condition_discontinuous_well_balancedness(x, t, element_id, equations::ShallowWaterEquationsWetDry2D) # Set the background values for velocity - v1 = 0.0 - v2 = 0.0 + H = equations.H0 + v1 = zero(H) + v2 = zero(H) x1, x2 = x - b = (1.5 / exp(0.5 * ((x1 - 1.0)^2 + (x2 + 1.0)^2)) + 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)) - @@ -175,9 +174,14 @@ function initial_condition_discontinuous_well_balancedness(x, t, element_id, equ 0.25 / exp(3.5 * ((x1 - 0.4)^2 + (x2 - 0.325)^2))) end - H = max(equations.H0, b + equations.threshold_limiter) + # H = max(equations.H0, b + equations.threshold_limiter) + + # return prim2cons(SVector(H, v1, v2, b), equations) + # Clip the initialization to avoid negative water heights and division by zero + h = max(equations.threshold_limiter, H - b) - return prim2cons(SVector(H, v1, v2, b), equations) + # Return the conservative variables + return SVector(h, h * v1, h * v2, b) end # point to the data we want to augment @@ -194,17 +198,17 @@ end summary_callback = SummaryCallback() -analysis_interval = 100 +analysis_interval = 10000 analysis_callback = AnalysisCallback(semi, interval = analysis_interval, extra_analysis_integrals = (lake_at_rest_error,)) alive_callback = AliveCallback(analysis_interval = analysis_interval) -save_solution = SaveSolutionCallback(interval = 100, # dt = 1.0, +save_solution = SaveSolutionCallback(dt = 1.0, save_initial_solution = true, save_final_solution = true) -stepsize_callback = StepsizeCallback(cfl = 1.0) +stepsize_callback = StepsizeCallback(cfl = 0.48) callbacks = CallbackSet(summary_callback, analysis_callback, @@ -219,5 +223,5 @@ stage_limiter! = PositivityPreservingLimiterShallowWater(variables = (Trixi.wate 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,); + save_everystep = false, callback = callbacks); summary_callback() # print the timer summary diff --git a/src/equations/shallow_water_wet_dry_2d.jl b/src/equations/shallow_water_wet_dry_2d.jl index fa0ef00..3943d67 100644 --- a/src/equations/shallow_water_wet_dry_2d.jl +++ b/src/equations/shallow_water_wet_dry_2d.jl @@ -5,6 +5,8 @@ @muladd begin #! format: noindent +# TODO: update this docstring (and probably others) to include `threshold_partially_wet` too + @doc raw""" ShallowWaterEquationsWetDry2D(; gravity, H0 = 0, threshold_limiter = nothing, threshold_wet = nothing) diff --git a/src/solvers/scratch_p4est.jl b/src/solvers/scratch_p4est.jl index 9cd2e43..f6f8edd 100644 --- a/src/solvers/scratch_p4est.jl +++ b/src/solvers/scratch_p4est.jl @@ -86,22 +86,29 @@ function Trixi.prolong2mortars!(cache, u, j_large = j_large_start element = neighbor_ids[3, mortar] for i in eachnode(dg) + # TODO: Do this is a better way that is agnostic with respect the number of equations + # Compute and save the sigma variable from Benov et al. (essentially H = h+b), + # momenta, and bottom topography into the buffer for projection. This ensures + # that we only project constant solution data in still water regions of the domain. + # OBS! A small shift is likely required to ensure we catch water heights close to the threshold + # TODO: This strategy from Benov et al. (https://doi.org/10.1016/j.jcp.2018.02.008) assumes that + # we know the constant background water height H0 which we pertub around. Fairly restrictive + # in practice but a good place to start with the development. We may nee to consider more sophisticated + # positivity preserving projections of the solution like those found in the ALE-DG community to remove + # this assumption. Then we might be able to directly project the water height `h` instead... + if u[1, i_large, j_large, element] > equations.threshold_limiter + eps() + 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 # from Benov et al. + # u_buffer[1, i] = u[4, i_large, j_large, element] + equations.threshold_limiter # do we need to project this instead? + end + u_buffer[2:4, i] = u[2:4, i_large, j_large, element] + for v in eachvariable(equations) - # Compute and save the sigma variable from Benov et al. (essentially H = h+b), - # momenta, and bottom topography into the buffer for projection. This ensures - # that we only project constant solution data in still water regions of the domain. - if u[1, i_large, j_large, element] > 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 # from Benov et al. - # u_buffer[1, i] = u[4, i_large, j_large, element] + equations.threshold_limiter # do we need to project this instead? - end - u_buffer[2:4, i] = u[2:4, i_large, j_large, element] # TODO: FIX ME! (or find a better way) # OBS! incredibly hacky way to save a copy of the (unprojected) parent solution on the mortar # where the mortar container has been modified to have extra storage space cache.mortars.u[3, v, 1, i, mortar] = u[v, i_large, j_large, element] - # u_buffer[v, i] = u[v, i_large, j_large, element] end i_large += i_large_step j_large += j_large_step @@ -127,6 +134,14 @@ 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) 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 + + # TODO: Can we call the stage limiter here for saftey / debugging? + # adapt this call?! + # limiter_shallow_water!(u, threshold::Real, variable, + # mesh::Trixi.AbstractMesh{2}, + # equations::ShallowWaterEquationsWetDry2D, dg::DGSEM, + # cache) + end return nothing @@ -192,7 +207,7 @@ function Trixi.calc_mortar_flux!(surface_flux_values, # once things are working # The projected solution from the large element is always stored in `u_rr` _, u_rr = Trixi.get_surface_node_vars(cache.mortars.u, equations, dg, - position, node, mortar) + position, node, mortar) # Compute conservative flux. Note, we use the surface flux here evaluated at the same # solution state to recover the physical flux at this point because the surface flux