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 new file mode 100644 index 0000000..8604f70 --- /dev/null +++ b/examples/p4est_2d_dgsem/elixir_shallowwater_well_balanced.jl @@ -0,0 +1,207 @@ + +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.0) + +function initial_condition_wb_testing(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))) + + return prim2cons(SVector(H, v1, v2, b), equations) +end + +initial_condition = initial_condition_wb_testing + +# # Wall BCs +# boundary_condition = Dict(:OuterCircle => boundary_condition_slip_wall) + +# # Dirichlet BCs +# boundary_condition_constant = BoundaryConditionDirichlet(initial_condition) +# boundary_condition = Dict(:OuterCircle => boundary_condition_constant) + +############################################################################### +# Get the DG approximation space + +# # ersing flux in both +# 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 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) + + +# # Fjordholms flux for testing +# surface_flux = (flux_fjordholm_etal, flux_nonconservative_fjordholm_etal) + +# Audusse HR with Rusanov flux +volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_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), +# flux_nonconservative_audusse_etal) + +# 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 + +# # Get the unstructured quad mesh from a file (downloads the file if not available locally) + +# default_mesh_file = joinpath(@__DIR__, "abaqus_outer_circle.inp") +# isfile(default_mesh_file) || +# download("https://gist.githubusercontent.com/andrewwinters5000/df92dd4986909927e96af23c37f6db5f/raw/8620823342f98c505a36351b210aab7f3b368041/abaqus_outer_circle.inp", +# default_mesh_file) +# mesh_file = default_mesh_file + +# mesh = P4estMesh{2}(mesh_file) + +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 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 + +# 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")) + +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) + +############################################################################### +# # Workaround to set a discontinuous bottom topography for debugging and testing. + +# alternative version of the initial conditinon used to setup a truly discontinuous +# bottom topography function for this academic testcase. +# The errors from the analysis callback are not important but the error for this lake at rest test case +# `∑|H0-(h+b)|` should be around machine roundoff +# In contrast to the usual signature of initial conditions, this one get passed the +# `element_id` explicitly. In particular, this initial conditions works as intended +# only for the specific mesh loaded above! +function initial_condition_discontinuous_well_balancedness(x, t, element_id, equations::ShallowWaterEquationsWetDry2D) + # Set the background values + 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))) + + # Setup a discontinuous bottom topography using the element id number + if element_id > 200 && 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 + + 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 = 100 +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_initial_solution = true, + save_final_solution = true) + +stepsize_callback = StepsizeCallback(cfl = 1.0) + +callbacks = CallbackSet(summary_callback, + analysis_callback, + alive_callback, + save_solution, + stepsize_callback) + +############################################################################### +# run the simulation +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), +# sol = solve(ode, SSPRK43(), + 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_wet_dry.jl b/examples/p4est_2d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl new file mode 100644 index 0000000..b5e733e --- /dev/null +++ b/examples/p4est_2d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl @@ -0,0 +1,227 @@ + +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 = 1.235) + +function initial_condition_wb_testing(x, t, equations::ShallowWaterEquationsWetDry2D) + # Calculate primitive variables + 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_wb_testing + +# # Wall BCs +# boundary_condition = Dict(:OuterCircle => boundary_condition_slip_wall) + +# # Dirichlet BCs +# boundary_condition_constant = BoundaryConditionDirichlet(initial_condition) +# boundary_condition = Dict(:OuterCircle => boundary_condition_constant) + +############################################################################### +# Get the DG approximation space + +# # ersing flux in both +# 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 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) + +# # 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) + +# Audusse HR with HLL flux +# surface_flux = (FluxHydrostaticReconstruction(flux_hll, hydrostatic_reconstruction_audusse_etal), +# flux_nonconservative_audusse_etal) + +# 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 + +# # Get the unstructured quad mesh from a file (downloads the file if not available locally) + +# default_mesh_file = joinpath(@__DIR__, "abaqus_outer_circle.inp") +# isfile(default_mesh_file) || +# download("https://gist.githubusercontent.com/andrewwinters5000/df92dd4986909927e96af23c37f6db5f/raw/8620823342f98c505a36351b210aab7f3b368041/abaqus_outer_circle.inp", +# default_mesh_file) +# mesh_file = default_mesh_file + +# mesh = P4estMesh{2}(mesh_file) + +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 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 + +# 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")) + +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) + +############################################################################### +# # Workaround to set a discontinuous bottom topography for debugging and testing. + +# alternative version of the initial conditinon used to setup a truly discontinuous +# bottom topography function for this academic testcase. +# The errors from the analysis callback are not important but the error for this lake at rest test case +# `∑|H0-(h+b)|` should be around machine roundoff +# In contrast to the usual signature of initial conditions, this one get passed the +# `element_id` explicitly. In particular, this initial conditions works as intended +# only for the specific mesh loaded above! +function initial_condition_discontinuous_well_balancedness(x, t, element_id, equations::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 + + # 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 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 = 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(dt = 1.0, + save_initial_solution = true, + save_final_solution = true) + +stepsize_callback = StepsizeCallback(cfl = 0.48) + +callbacks = CallbackSet(summary_callback, + analysis_callback, + alive_callback, + save_solution, + stepsize_callback) + +stage_limiter! = PositivityPreservingLimiterShallowWater(variables = (Trixi.waterheight,)) + +############################################################################### +# run the simulation +sol = solve(ode, SSPRK43(stage_limiter!), + 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/src/TrixiShallowWater.jl b/src/TrixiShallowWater.jl index 86cd7e9..a2a4721 100644 --- a/src/TrixiShallowWater.jl +++ b/src/TrixiShallowWater.jl @@ -1,8 +1,8 @@ module TrixiShallowWater -# While `using Trixi` makes all exported symbols available, in order to extend a method from the +# While `using Trixi` makes all exported symbols available, in order to extend a method from the # `Trixi.jl` module, symbols need to be explicitly qualified with `Trixi.function_name`. -# For more information, see +# For more information, see # https://github.com/trixi-framework/TrixiShallowWater.jl/pull/10#discussion_r1433720559 using Trixi # Import additional symbols that are not exported by Trixi.jl @@ -16,6 +16,7 @@ include("equations/equations.jl") include("equations/numerical_fluxes.jl") include("callbacks_stage/callbacks_stage.jl") include("solvers/indicators.jl") +include("solvers/scratch_p4est.jl") # Export types/functions that define the public API of TrixiShallowWater.jl export ShallowWaterEquationsWetDry1D, ShallowWaterEquationsWetDry2D, diff --git a/src/equations/shallow_water_wet_dry_2d.jl b/src/equations/shallow_water_wet_dry_2d.jl index 363fceb..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) @@ -31,10 +33,10 @@ Also, there are two thresholds which prevent numerical problems as well as insta have to be passed, as default values are defined within the struct. The first one, `threshold_limiter`, is used in [`PositivityPreservingLimiterShallowWater`](@ref) on the water height, as a (small) shift on the initial condition and cutoff before the next time step. The second one, `threshold_wet`, is applied on the water height to -define when the flow is "wet" before calculating the numerical flux. A third -`threshold_partially_wet` is applied on the water height to define "partially wet" elements in +define when the flow is "wet" before calculating the numerical flux. A third +`threshold_partially_wet` is applied on the water height to define "partially wet" elements in [`IndicatorHennemannGassnerShallowWater`](@ref), that are then calculated with a pure FV method to -ensure well-balancedness. For `Float64` no threshold needs to be passed, as default values are +ensure well-balancedness. For `Float64` no threshold needs to be passed, as default values are defined within the struct. For other number formats `threshold_partially_wet` must be provided. The bottom topography function ``b(x,y)`` is set inside the initial condition routine @@ -69,11 +71,11 @@ struct ShallowWaterEquationsWetDry2D{RealT <: Real} <: # before calculating the numerical flux. # Default is 5*eps() which in double precision is ≈1e-15. threshold_wet::RealT - # `threshold_partially_wet` used in `IndicatorHennemannGassnerShallowWater` on the water height - # to define "partially wet" elements. Those elements are calculated with a pure FV method to - # ensure well-balancedness. Default in double precision is 1e-4. + # `threshold_partially_wet` used in `IndicatorHennemannGassnerShallowWater` on the water height + # to define "partially wet" elements. Those elements are calculated with a pure FV method to + # ensure well-balancedness. Default in double precision is 1e-4. threshold_partially_wet::RealT - # Standard shallow water equations for dispatch on Trixi.jl functions + # Standard shallow water equations for dispatch on Trixi.jl functions basic_swe::ShallowWaterEquations2D{RealT} end @@ -96,8 +98,8 @@ function ShallowWaterEquationsWetDry2D(; gravity_constant, H0 = zero(gravity_con if threshold_partially_wet === nothing threshold_partially_wet = default_threshold_partially_wet(T) end - # Construct the standard SWE for dispatch. Even though the `basic_swe` already store the - # gravity constant and the total water height, we store an extra copy in + # Construct the standard SWE for dispatch. Even though the `basic_swe` already store the + # gravity constant and the total water height, we store an extra copy in # `ShallowWaterEquationsWetDry2D` for convenience. basic_swe = ShallowWaterEquations2D(gravity_constant = gravity_constant, H0 = H0) @@ -114,10 +116,10 @@ Trixi.varnames(::typeof(cons2cons), ::ShallowWaterEquationsWetDry2D) = ("h", "h_ Trixi.varnames(::typeof(cons2prim), ::ShallowWaterEquationsWetDry2D) = ("H", "v1", "v2", "b") -# This equation set extends the basic ShallowWaterEquations2D from Trixi.jl with additional -# functionality for wet/dry transitions. Since many functions correspond to the fully wet case, we +# This equation set extends the basic ShallowWaterEquations2D from Trixi.jl with additional +# functionality for wet/dry transitions. Since many functions correspond to the fully wet case, we # make use of the existing functionality and introduce a number of wrapper functions, that dispatch -# to the ShallowWaterEquations2D. +# to the ShallowWaterEquations2D. # Set initial conditions at physical location `x` for time `t` """ @@ -564,13 +566,13 @@ contravariant vector (normal direction) at the current node and the averaged one. This is different from numerical fluxes used to discretize conservative terms. -This is a modified version of [`flux_nonconservative_wintermeyer_etal`](@ref) that gives entropy -conservation and well-balancedness in both the volume and surface when combined with +This is a modified version of [`flux_nonconservative_wintermeyer_etal`](@ref) that gives entropy +conservation and well-balancedness in both the volume and surface when combined with [`flux_wintermeyer_etal`](@ref). For further details see: - Patrick Ersing, Andrew R. Winters (2023) - An entropy stable discontinuous Galerkin method for the two-layer shallow water equations on + An entropy stable discontinuous Galerkin method for the two-layer shallow water equations on curvilinear meshes [DOI: 10.1007/s10915-024-02451-2](https://doi.org/10.1007/s10915-024-02451-2) """ diff --git a/src/solvers/scratch_p4est.jl b/src/solvers/scratch_p4est.jl new file mode 100644 index 0000000..f6f8edd --- /dev/null +++ b/src/solvers/scratch_p4est.jl @@ -0,0 +1,411 @@ +# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). +# Since these FMAs can increase the performance of many numerical algorithms, +# we need to opt-in explicitly. +# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. +@muladd begin +#! format: noindent + +# TODO: once working the mortar methods could likely be extended to the other +# equations 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 mortar type +# and called from the basic `create_cache` method at the top. +# Extra storage is necessary for the numerical flux plus nonconservative term +# on either side of the large and small element mortars. We also require a work +# array to compute the flux with the project large element solutions on each mortar +# to ensure we project back the flux penalty. +function Trixi.create_cache(mesh::Union{P4estMesh{2}, T8codeMesh{2}}, + equations::ShallowWaterEquationsWetDry2D, + mortar_l2::Trixi.LobattoLegendreMortarL2, uEltype) + # TODO: Taal performance using different types + MA2d = Trixi.MArray{Tuple{nvariables(equations), nnodes(mortar_l2)}, + uEltype, 2, + nvariables(equations) * nnodes(mortar_l2)} + fstar_primary_upper_threaded = MA2d[MA2d(undef) for _ in 1:Threads.nthreads()] + fstar_primary_lower_threaded = MA2d[MA2d(undef) for _ in 1:Threads.nthreads()] + fstar_secondary_upper_threaded = MA2d[MA2d(undef) for _ in 1:Threads.nthreads()] + fstar_secondary_lower_threaded = MA2d[MA2d(undef) for _ in 1:Threads.nthreads()] + f_upper_threaded = MA2d[MA2d(undef) for _ in 1:Threads.nthreads()] + f_lower_threaded = MA2d[MA2d(undef) for _ in 1:Threads.nthreads()] + u_threaded = MA2d[MA2d(undef) for _ in 1:Threads.nthreads()] + f_threaded = MA2d[MA2d(undef) for _ in 1:Threads.nthreads()] + + (; fstar_primary_upper_threaded, fstar_primary_lower_threaded, + fstar_secondary_upper_threaded, fstar_secondary_lower_threaded, + f_upper_threaded, f_lower_threaded, u_threaded, f_threaded) +end + +function Trixi.prolong2mortars!(cache, u, + mesh::Union{P4estMesh{2}, T8codeMesh{2}}, + equations::ShallowWaterEquationsWetDry2D, + mortar_l2::Trixi.LobattoLegendreMortarL2, + surface_integral, dg::DGSEM) + @unpack neighbor_ids, node_indices = cache.mortars + index_range = eachnode(dg) + + 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. + small_indices = node_indices[1, mortar] + + i_small_start, i_small_step = Trixi.index_to_start_step_2d(small_indices[1], + index_range) + j_small_start, j_small_step = Trixi.index_to_start_step_2d(small_indices[2], + index_range) + + for position in 1:2 + i_small = i_small_start + j_small = j_small_start + element = neighbor_ids[position, mortar] + for i in eachnode(dg) + for v in eachvariable(equations) + cache.mortars.u[1, v, position, i, mortar] = u[v, i_small, j_small, + element] + end + i_small += i_small_step + j_small += j_small_step + end + end + + # Buffer to copy solution values of the large element in the correct orientation + # before interpolating + u_buffer = cache.u_threaded[Threads.threadid()] + + # Copy solution of large element face to buffer in the + # correct orientation + large_indices = node_indices[2, mortar] + + i_large_start, i_large_step = Trixi.index_to_start_step_2d(large_indices[1], + index_range) + j_large_start, j_large_step = Trixi.index_to_start_step_2d(large_indices[2], + index_range) + + i_large = i_large_start + 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) + # 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] + 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 can modify the values + # in the first solution variable to no longer be the sigma variable of + # Benov et al. and instead be the conservative water height variable `h`. + # Basically, unpacking the sigma variable to create the projected local water + # height from Eq. 41 in Benov et al. + # TODO: My main hope was that this avoids allocations + # TODO: My other hope is that such a strategy will make this code extensible + # to the multilayer equations as well. + 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) + 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 +end + +# TODO: The procedure below allocates a lot. Need a better way to avoid this +function Trixi.calc_mortar_flux!(surface_flux_values, + mesh::Union{P4estMesh{2}, T8codeMesh{2}}, + nonconservative_terms, + equations::ShallowWaterEquationsWetDry2D, + mortar_l2::Trixi.LobattoLegendreMortarL2, + surface_integral, dg::DG, cache) + @unpack neighbor_ids, node_indices = cache.mortars + @unpack contravariant_vectors = cache.elements + @unpack fstar_primary_upper_threaded, fstar_primary_lower_threaded, fstar_secondary_upper_threaded, fstar_secondary_lower_threaded, f_upper_threaded, f_lower_threaded = cache + surface_flux, nonconservative_flux = surface_integral.surface_flux + index_range = eachnode(dg) + + Trixi.@threaded for mortar in Trixi.eachmortar(dg, cache) + # Choose thread-specific pre-allocated containers + # for the numerical flux on the small elements as well as + # the physical flux evaluated at the projected solution + # from the large element + fstar_primary = (fstar_primary_lower_threaded[Threads.threadid()], + fstar_primary_upper_threaded[Threads.threadid()]) + + fstar_secondary = (fstar_secondary_lower_threaded[Threads.threadid()], + fstar_secondary_upper_threaded[Threads.threadid()]) + + f = (f_lower_threaded[Threads.threadid()], + f_upper_threaded[Threads.threadid()]) + + # Get index information on the small elements + small_indices = node_indices[1, mortar] + small_direction = Trixi.indices2direction(small_indices) + + i_small_start, i_small_step = Trixi.index_to_start_step_2d(small_indices[1], + index_range) + j_small_start, j_small_step = Trixi.index_to_start_step_2d(small_indices[2], + index_range) + + for position in 1:2 + i_small = i_small_start + j_small = j_small_start + element = neighbor_ids[position, mortar] + for node in eachnode(dg) + # Get the normal direction on the small element. + # Note, contravariant vectors at interfaces in negative coordinate direction + # are pointing inwards. This is handled by `get_normal_direction`. + normal_direction = Trixi.get_normal_direction(small_direction, + contravariant_vectors, + i_small, j_small, element) + + # Compute the numerical flux in the normal direction on the small elements + Trixi.calc_mortar_flux!(fstar_primary, fstar_secondary, + mesh, nonconservative_terms, equations, + surface_integral, dg, cache, + mortar, position, normal_direction, + node) + + ##### + # TODO: Possibly move this code block indicated by "#####" into a new function + # 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) + + # 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 + # has in-built mechanisms to avoid division by zero in dry regions whereas `Trixi.flux` + # does not have such mechanisms to desingularize the velocity computation. + flux = surface_flux(u_rr, u_rr, normal_direction, equations) + + # Compute nonconservative flux and add it to the conservative flux. + # The nonconservative flux is scaled by a factor of 0.5 based on + # the interpretation of global SBP operators coupled discontinuously via + # central fluxes/SATs + noncons = nonconservative_flux(u_rr, u_rr, + normal_direction, normal_direction, + equations) + + flux_plus_noncons = flux + 0.5 * noncons + + # Copy to the physical flux buffer + set_node_vars!(f[position], flux_plus_noncons, equations, dg, node) + ##### + + i_small += i_small_step + j_small += j_small_step + end + end + + # Buffer to interpolate flux values of the large element to before + # copying in the correct orientation + u_buffer = cache.u_threaded[Threads.threadid()] + + # In calc_interface_flux!, the interface flux is computed once over each + # interface using the normal from the "primary" element. The result is then + # passed back to the "secondary" element, flipping the sign to account for the + # change in the normal direction. For mortars, this sign flip occurs in + # "mortar_fluxes_to_elements!" instead. + Trixi.mortar_fluxes_to_elements!(surface_flux_values, + mesh, equations, mortar_l2, dg, cache, + mortar, fstar_primary, fstar_secondary, + f, u_buffer) + end + + return nothing +end + +# Inlined version of the mortar flux computation on small elements for equations with conservative and +# nonconservative terms +@inline function Trixi.calc_mortar_flux!(fstar_primary, fstar_secondary, + mesh::Union{P4estMesh{2}, T8codeMesh{2}}, + nonconservative_terms::True, + equations::ShallowWaterEquationsWetDry2D, + surface_integral, dg::DG, cache, + mortar_index, position_index, normal_direction, + node_index) + @unpack u = cache.mortars + surface_flux, nonconservative_flux = surface_integral.surface_flux + + u_ll, u_rr = Trixi.get_surface_node_vars(u, equations, dg, position_index, node_index, + mortar_index) + + # General idea of the mortar flux plus nonconserative term treatment is the following: + # (1) unpack the sigma to create the water height on the mortar from Eq. 41 in Benov et al. + # (2) perform hydrostatic reconstruction on the mortar solution + # (3) use these HR quantities to compute the flux and nonconservative terms + # The first step is actually done in the `prolong2mortars!`routine. + # The other two steps occur within the surface flux and nonconservative computations. + + # Compute conservative flux + flux = surface_flux(u_ll, u_rr, normal_direction, equations) + + # Compute nonconservative flux and add it to the conservative flux. + # The nonconservative flux is scaled by a factor of 0.5 based on + # the interpretation of global SBP operators coupled discontinuously via + # central fluxes/SATs + noncons_primary = nonconservative_flux(u_ll, u_rr, + normal_direction, normal_direction, + equations) + noncons_secondary = nonconservative_flux(u_rr, u_ll, + normal_direction, normal_direction, + equations) + + flux_plus_noncons_primary = flux + 0.5 * noncons_primary + flux_plus_noncons_secondary = flux + 0.5 * noncons_secondary + + # Copy results to the buffers + set_node_vars!(fstar_primary[position_index], flux_plus_noncons_primary, equations, dg, node_index) + set_node_vars!(fstar_secondary[position_index], flux_plus_noncons_secondary, equations, dg, node_index) +end + +@inline function Trixi.mortar_fluxes_to_elements!(surface_flux_values, + mesh::Union{P4estMesh{2}, T8codeMesh{2}}, + equations::ShallowWaterEquationsWetDry2D, + mortar_l2::Trixi.LobattoLegendreMortarL2, + dg::DGSEM, cache, mortar, + fstar_primary, fstar_secondary, + f_large, u_buffer) + @unpack contravariant_vectors = cache.elements + @unpack neighbor_ids, node_indices = cache.mortars + surface_flux, nonconservative_flux = dg.surface_integral.surface_flux + + # Copy surface flux data from small to small + small_indices = node_indices[1, mortar] + small_direction = Trixi.indices2direction(small_indices) + + for position in 1:2 + element = neighbor_ids[position, mortar] + for i in eachnode(dg) + for v in eachvariable(equations) + surface_flux_values[v, i, small_direction, element] = fstar_primary[position][v, + i] + end + end + end + + # Project small numerical fluxes and physical flux penalty computed on the projected + # large element solution back onto large element. + # This is basically Eq. (46) from Benov et al. where the factor of 1/2 is already + # already included in `reverse_upper` and `reverse_lower` operators. + Trixi.multiply_dimensionwise!(u_buffer, + mortar_l2.reverse_upper, (fstar_secondary[2] .- f_large[2]), + mortar_l2.reverse_lower, (fstar_secondary[1] .- f_large[1])) + + # The flux is calculated in the outward direction of the small elements, + # so the sign must be switched to get the flux in outward direction + # of the large element. + # The contravariant vectors of the large element (and therefore the normal + # vectors of the large element as well) are twice as large as the + # contravariant vectors of the small elements. Therefore, the flux needs + # to be scaled by a factor of 2 to obtain the flux of the large element. + u_buffer .*= -2 + + # Copy interpolated flux values from buffer to large element face in the + # correct orientation. + # Note that the index of the small sides will always run forward but + # the index of the large side might need to run backwards for flipped sides. + large_element = neighbor_ids[3, mortar] + large_indices = node_indices[2, mortar] + large_direction = Trixi.indices2direction(large_indices) + + # TODO: We need to store the unprojected solution in the mortars such that we have access to them + # here when we go to compute the flux on the parent elements and remove the physical flux evaluated + # at the upprojected solution state that is present from the volume integral computation + index_range = eachnode(dg) + i_large_start, i_large_step = Trixi.index_to_start_step_2d(large_indices[1], + index_range) + j_large_start, j_large_step = Trixi.index_to_start_step_2d(large_indices[2], + index_range) + + i_large = i_large_start + j_large = j_large_start + # maybe similar(u_buffer) instead? Would avoid adding to the cache + flux_buffer = cache.f_threaded[Threads.threadid()] + for node in eachnode(dg) + # Get the proper normal_direction now that we are back computing on the large element + normal_direction = Trixi.get_normal_direction(large_direction, + contravariant_vectors, + i_large, j_large, large_element) + + # 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 + # has in-built mechanisms to avoid division by zero in dry regions whereas `Trixi.flux` + # does not have such mechanisms to desingularize the velocity computation. + flux = surface_flux(view(cache.mortars.u, 3, :, 1, node, mortar), + view(cache.mortars.u, 3, :, 1, node, mortar), + normal_direction, equations) + + noncons = nonconservative_flux(view(cache.mortars.u, 3, :, 1, node, mortar), + view(cache.mortars.u, 3, :, 1, node, mortar), + normal_direction, normal_direction, + equations) + + flux_plus_noncons = flux + 0.5 * noncons + + set_node_vars!(flux_buffer, flux_plus_noncons, equations, dg, node) + + i_large += i_large_step + j_large += j_large_step + end + + # Compute the surface flux values on the large element. Note, we have projected + # the flux penalty from the small element mortars back onto the large element side. + # However, the large element has a local contribution to the strong form flux penalty + # that must be removed to ensure consistency. This is done via the values precomputed + # and stored above in the `flux_buffer` which is the physical flux and nonconservative + # term evaluated at the unprojected large element solution and `normal_direction`. + if :i_backward in large_indices + for i in eachnode(dg) + for v in eachvariable(equations) + surface_flux_values[v, end + 1 - i, large_direction, large_element] = u_buffer[v, i] + flux_buffer[v, i] + end + end + else + for i in eachnode(dg) + for v in eachvariable(equations) + surface_flux_values[v, i, large_direction, large_element] = u_buffer[v, i] + flux_buffer[v, i] + end + end + end + + return nothing +end +end # @muladd diff --git a/src/solvers/scratch_tree.jl b/src/solvers/scratch_tree.jl new file mode 100644 index 0000000..ecddb59 --- /dev/null +++ b/src/solvers/scratch_tree.jl @@ -0,0 +1,269 @@ +# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). +# Since these FMAs can increase the performance of many numerical algorithms, +# we need to opt-in explicitly. +# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. +@muladd begin +#! format: noindent + +function prolong2mortars!(cache, u, + mesh::TreeMesh{2}, equations, + mortar_l2::LobattoLegendreMortarL2, surface_integral, + dg::DGSEM) + @threaded for mortar in eachmortar(dg, cache) + large_element = cache.mortars.neighbor_ids[3, mortar] + upper_element = cache.mortars.neighbor_ids[2, mortar] + lower_element = cache.mortars.neighbor_ids[1, mortar] + + # Copy solution small to small + if cache.mortars.large_sides[mortar] == 1 # -> small elements on right side + if cache.mortars.orientations[mortar] == 1 + # L2 mortars in x-direction + for l in eachnode(dg) + for v in eachvariable(equations) + cache.mortars.u_upper[2, v, l, mortar] = u[v, 1, l, + upper_element] + cache.mortars.u_lower[2, v, l, mortar] = u[v, 1, l, + lower_element] + end + end + else + # L2 mortars in y-direction + for l in eachnode(dg) + for v in eachvariable(equations) + cache.mortars.u_upper[2, v, l, mortar] = u[v, l, 1, + upper_element] + cache.mortars.u_lower[2, v, l, mortar] = u[v, l, 1, + lower_element] + end + end + end + else # large_sides[mortar] == 2 -> small elements on left side + if cache.mortars.orientations[mortar] == 1 + # L2 mortars in x-direction + for l in eachnode(dg) + for v in eachvariable(equations) + cache.mortars.u_upper[1, v, l, mortar] = u[v, nnodes(dg), l, + upper_element] + cache.mortars.u_lower[1, v, l, mortar] = u[v, nnodes(dg), l, + lower_element] + end + end + else + # L2 mortars in y-direction + for l in eachnode(dg) + for v in eachvariable(equations) + cache.mortars.u_upper[1, v, l, mortar] = u[v, l, nnodes(dg), + upper_element] + cache.mortars.u_lower[1, v, l, mortar] = u[v, l, nnodes(dg), + lower_element] + end + end + end + end + + # Interpolate large element face data to small interface locations + if cache.mortars.large_sides[mortar] == 1 # -> large element on left side + leftright = 1 + if cache.mortars.orientations[mortar] == 1 + # L2 mortars in x-direction + u_large = view(u, :, nnodes(dg), :, large_element) + element_solutions_to_mortars!(cache.mortars, mortar_l2, leftright, + mortar, u_large) + else + # L2 mortars in y-direction + u_large = view(u, :, :, nnodes(dg), large_element) + element_solutions_to_mortars!(cache.mortars, mortar_l2, leftright, + mortar, u_large) + end + else # large_sides[mortar] == 2 -> large element on right side + leftright = 2 + if cache.mortars.orientations[mortar] == 1 + # L2 mortars in x-direction + u_large = view(u, :, 1, :, large_element) + element_solutions_to_mortars!(cache.mortars, mortar_l2, leftright, + mortar, u_large) + else + # L2 mortars in y-direction + u_large = view(u, :, :, 1, large_element) + element_solutions_to_mortars!(cache.mortars, mortar_l2, leftright, + mortar, u_large) + end + end + end + + return nothing +end + +@inline function element_solutions_to_mortars!(mortars, + mortar_l2::LobattoLegendreMortarL2, + leftright, mortar, + u_large::AbstractArray{<:Any, 2}) + multiply_dimensionwise!(view(mortars.u_upper, leftright, :, :, mortar), + mortar_l2.forward_upper, u_large) + multiply_dimensionwise!(view(mortars.u_lower, leftright, :, :, mortar), + mortar_l2.forward_lower, u_large) + return nothing +end + +function calc_mortar_flux!(surface_flux_values, + mesh::TreeMesh{2}, + nonconservative_terms::True, equations, + mortar_l2::LobattoLegendreMortarL2, + surface_integral, dg::DG, cache) + surface_flux, nonconservative_flux = surface_integral.surface_flux + @unpack u_lower, u_upper, orientations, large_sides = cache.mortars + @unpack fstar_upper_threaded, fstar_lower_threaded = cache + + @threaded for mortar in eachmortar(dg, cache) + # Choose thread-specific pre-allocated container + fstar_upper = fstar_upper_threaded[Threads.threadid()] + fstar_lower = fstar_lower_threaded[Threads.threadid()] + + # Calculate fluxes + orientation = orientations[mortar] + calc_fstar!(fstar_upper, equations, surface_flux, dg, u_upper, mortar, + orientation) + calc_fstar!(fstar_lower, equations, surface_flux, dg, u_lower, mortar, + orientation) + + # Add nonconservative fluxes. + # These need to be adapted on the geometry (left/right) since the order of + # the arguments matters, based on the global SBP operator interpretation. + # The same interpretation (global SBP operators coupled discontinuously via + # central fluxes/SATs) explains why we need the factor 0.5. + # Alternatively, you can also follow the argumentation of Bohm et al. 2018 + # ("nonconservative diamond flux") + if large_sides[mortar] == 1 # -> small elements on right side + for i in eachnode(dg) + # Pull the left and right solutions + u_upper_ll, u_upper_rr = get_surface_node_vars(u_upper, equations, dg, + i, mortar) + u_lower_ll, u_lower_rr = get_surface_node_vars(u_lower, equations, dg, + i, mortar) + # Call pointwise nonconservative term + noncons_upper = nonconservative_flux(u_upper_ll, u_upper_rr, + orientation, equations) + noncons_lower = nonconservative_flux(u_lower_ll, u_lower_rr, + orientation, equations) + # Add to primary and secondary temporary storage + multiply_add_to_node_vars!(fstar_upper, 0.5, noncons_upper, equations, + dg, i) + multiply_add_to_node_vars!(fstar_lower, 0.5, noncons_lower, equations, + dg, i) + end + else # large_sides[mortar] == 2 -> small elements on the left + for i in eachnode(dg) + # Pull the left and right solutions + u_upper_ll, u_upper_rr = get_surface_node_vars(u_upper, equations, dg, + i, mortar) + u_lower_ll, u_lower_rr = get_surface_node_vars(u_lower, equations, dg, + i, mortar) + # Call pointwise nonconservative term + noncons_upper = nonconservative_flux(u_upper_rr, u_upper_ll, + orientation, equations) + noncons_lower = nonconservative_flux(u_lower_rr, u_lower_ll, + orientation, equations) + # Add to primary and secondary temporary storage + multiply_add_to_node_vars!(fstar_upper, 0.5, noncons_upper, equations, + dg, i) + multiply_add_to_node_vars!(fstar_lower, 0.5, noncons_lower, equations, + dg, i) + end + end + + mortar_fluxes_to_elements!(surface_flux_values, + mesh, equations, mortar_l2, dg, cache, + mortar, fstar_upper, fstar_lower) + end + + return nothing +end + +@inline function calc_fstar!(destination::AbstractArray{<:Any, 2}, equations, + surface_flux, dg::DGSEM, + u_interfaces, interface, orientation) + for i in eachnode(dg) + # Call pointwise two-point numerical flux function + u_ll, u_rr = get_surface_node_vars(u_interfaces, equations, dg, i, interface) + flux = surface_flux(u_ll, u_rr, orientation, equations) + + # Copy flux to left and right element storage + set_node_vars!(destination, flux, equations, dg, i) + end + + return nothing +end + +@inline function mortar_fluxes_to_elements!(surface_flux_values, + mesh::TreeMesh{2}, equations, + mortar_l2::LobattoLegendreMortarL2, + dg::DGSEM, cache, + mortar, fstar_upper, fstar_lower) + large_element = cache.mortars.neighbor_ids[3, mortar] + upper_element = cache.mortars.neighbor_ids[2, mortar] + lower_element = cache.mortars.neighbor_ids[1, mortar] + + # Copy flux small to small + if cache.mortars.large_sides[mortar] == 1 # -> small elements on right side + if cache.mortars.orientations[mortar] == 1 + # L2 mortars in x-direction + direction = 1 + else + # L2 mortars in y-direction + direction = 3 + end + else # large_sides[mortar] == 2 -> small elements on left side + if cache.mortars.orientations[mortar] == 1 + # L2 mortars in x-direction + direction = 2 + else + # L2 mortars in y-direction + direction = 4 + end + end + surface_flux_values[:, :, direction, upper_element] .= fstar_upper + surface_flux_values[:, :, direction, lower_element] .= fstar_lower + + # Project small fluxes to large element + if cache.mortars.large_sides[mortar] == 1 # -> large element on left side + if cache.mortars.orientations[mortar] == 1 + # L2 mortars in x-direction + direction = 2 + else + # L2 mortars in y-direction + direction = 4 + end + else # large_sides[mortar] == 2 -> large element on right side + if cache.mortars.orientations[mortar] == 1 + # L2 mortars in x-direction + direction = 1 + else + # L2 mortars in y-direction + direction = 3 + end + end + + # TODO: Taal performance + # for v in eachvariable(equations) + # # The code below is semantically equivalent to + # # surface_flux_values[v, :, direction, large_element] .= + # # (mortar_l2.reverse_upper * fstar_upper[v, :] + mortar_l2.reverse_lower * fstar_lower[v, :]) + # # but faster and does not allocate. + # # Note that `true * some_float == some_float` in Julia, i.e. `true` acts as + # # a universal `one`. Hence, the second `mul!` means "add the matrix-vector + # # product to the current value of the destination". + # @views mul!(surface_flux_values[v, :, direction, large_element], + # mortar_l2.reverse_upper, fstar_upper[v, :]) + # @views mul!(surface_flux_values[v, :, direction, large_element], + # mortar_l2.reverse_lower, fstar_lower[v, :], true, true) + # end + # The code above could be replaced by the following code. However, the relative efficiency + # depends on the types of fstar_upper/fstar_lower and dg.l2mortar_reverse_upper. + # Using StaticArrays for both makes the code above faster for common test cases. + multiply_dimensionwise!(view(surface_flux_values, :, :, direction, large_element), + mortar_l2.reverse_upper, fstar_upper, + mortar_l2.reverse_lower, fstar_lower) + + return nothing +end +end # @muladd