diff --git a/.github/workflows/FormatCheck.yml b/.github/workflows/FormatCheck.yml index 1310d5d..8ecb87a 100644 --- a/.github/workflows/FormatCheck.yml +++ b/.github/workflows/FormatCheck.yml @@ -25,8 +25,10 @@ jobs: # This will use the latest version by default but you can set the version like so: # # julia -e 'using Pkg; Pkg.add(PackageSpec(name = "JuliaFormatter", version = "0.13.0"))' + # + # TODO: remove version restriction once formatting has sorted itself out run: | - julia -e 'using Pkg; Pkg.add(PackageSpec(name = "JuliaFormatter"))' + julia -e 'using Pkg; Pkg.add(PackageSpec(name = "JuliaFormatter", version="1.0.45"))' julia -e 'using JuliaFormatter; format(".")' - name: Format check run: | diff --git a/Project.toml b/Project.toml index 1ec3ece..6ea3d75 100644 --- a/Project.toml +++ b/Project.toml @@ -4,14 +4,18 @@ authors = ["Andrew R. Winters ", "Michael Schlottke- version = "0.1.0-pre" [deps] +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221" Static = "aedffcd0-7271-4cad-89d0-dc628f76c6d3" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Trixi = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb" +Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" [compat] +LinearAlgebra = "1" MuladdMacro = "0.2.2" Static = "0.3, 0.4, 0.5, 0.6, 0.7, 0.8" StaticArrays = "1" -Trixi = "0.5.17, 0.6" +Trixi = "0.7" julia = "1.8" +Printf = "1" diff --git a/examples/structured_2d_dgsem/elixir_shallowwater_conical_island.jl b/examples/structured_2d_dgsem/elixir_shallowwater_conical_island.jl new file mode 100644 index 0000000..6696130 --- /dev/null +++ b/examples/structured_2d_dgsem/elixir_shallowwater_conical_island.jl @@ -0,0 +1,113 @@ + +using OrdinaryDiffEq +using Trixi +using TrixiShallowWater + +############################################################################### +# Semidiscretization of the shallow water equations + +equations = ShallowWaterEquationsWetDry2D(gravity_constant = 9.81, H0 = 1.4) + +""" + initial_condition_conical_island(x, t, equations::ShallowWaterEquationsWetDry2D) + +Initial condition for the [`ShallowWaterEquationsWetDry2D`](@ref) to test the [`hydrostatic_reconstruction_chen_noelle`](@ref) +and its handling of discontinuous water heights at the start in combination with wetting and +drying. The bottom topography is given by a conical island in the middle of the domain. Around that +island, there is a cylindrical water column at t=0 and the rest of the domain is dry. This +discontinuous water height is smoothed by a logistic function. This simulation uses periodic +boundary conditions. +""" +function initial_condition_conical_island(x, t, equations::ShallowWaterEquationsWetDry2D) + # Set the background values + + v1 = 0.0 + v2 = 0.0 + + x1, x2 = x + b = max(0.1, 1.0 - 4.0 * sqrt(x1^2 + x2^2)) + + # use a logistic function to transfer water height value smoothly + L = equations.H0 # maximum of function + x0 = 0.3 # center point of function + k = -25.0 # sharpness of transfer + + H = max(b, L / (1.0 + exp(-k * (sqrt(x1^2 + x2^2) - x0)))) + + # 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 v1) / h. Therefore, a small dry state threshold + # with a default value of 500*eps() ≈ 1e-13 in double precision, is set in the constructor above + # for the ShallowWaterEquationsWetDry and added to the initial condition if h = 0. + # This default value can be changed within the constructor call depending on the simulation setup. + H = max(H, b + equations.threshold_limiter) + return prim2cons(SVector(H, v1, v2, b), equations) +end + +initial_condition = initial_condition_conical_island + +############################################################################### +# 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) + +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 StructuredMesh and setup a periodic mesh + +coordinates_min = (-1.0, -1.0) +coordinates_max = (1.0, 1.0) + +cells_per_dimension = (16, 16) + +mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max) + +# Create the semi discretization object +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +############################################################################### +# ODE solver + +tspan = (0.0, 10.0) +ode = semidiscretize(semi, tspan) + +############################################################################### +# Callbacks + +summary_callback = SummaryCallback() + +analysis_interval = 1000 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 100, + save_initial_solution = true, + save_final_solution = true) + +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution) + +############################################################################### +# run the simulation + +stage_limiter! = PositivityPreservingLimiterShallowWater(variables = (Trixi.waterheight,)) + +sol = solve(ode, SSPRK43(stage_limiter!); + ode_default_options()..., callback = callbacks); + +summary_callback() # print the timer summary diff --git a/examples/structured_2d_dgsem/elixir_shallowwater_parabolic_bowl.jl b/examples/structured_2d_dgsem/elixir_shallowwater_parabolic_bowl.jl new file mode 100644 index 0000000..1d90985 --- /dev/null +++ b/examples/structured_2d_dgsem/elixir_shallowwater_parabolic_bowl.jl @@ -0,0 +1,118 @@ + +using OrdinaryDiffEq +using Trixi +using TrixiShallowWater + +############################################################################### +# Semidiscretization of the shallow water equations + +equations = ShallowWaterEquationsWetDry2D(gravity_constant = 9.81) + +""" + initial_condition_parabolic_bowl(x, t, equations:: ShallowWaterEquationsWetDry2D) + +Well-known initial condition to test the [`hydrostatic_reconstruction_chen_noelle`](@ref) and its +wet-dry mechanics. This test has an analytical solution. The initial condition is defined by the +analytical solution at time t=0. The bottom topography defines a bowl and the water level is given +by an oscillating lake. + +The original test and its analytical solution were first presented in +- William C. Thacker (1981) + Some exact solutions to the nonlinear shallow-water wave equations + [DOI: 10.1017/S0022112081001882](https://doi.org/10.1017/S0022112081001882). + +The particular setup below is taken from Section 6.2 of +- Niklas Wintermeyer, Andrew R. Winters, Gregor J. Gassner and Timothy Warburton (2018) + An entropy stable discontinuous Galerkin method for the shallow water equations on + curvilinear meshes with wet/dry fronts accelerated by GPUs + [DOI: 10.1016/j.jcp.2018.08.038](https://doi.org/10.1016/j.jcp.2018.08.038). +""" +function initial_condition_parabolic_bowl(x, t, equations::ShallowWaterEquationsWetDry2D) + a = 1.0 + h_0 = 0.1 + sigma = 0.5 + ω = sqrt(2 * equations.gravity * h_0) / a + + v1 = -sigma * ω * sin(ω * t) + v2 = sigma * ω * cos(ω * t) + + b = h_0 * ((x[1])^2 + (x[2])^2) / a^2 + + H = sigma * h_0 / a^2 * (2 * x[1] * cos(ω * t) + 2 * x[2] * sin(ω * t) - sigma) + h_0 + + # It is mandatory to shift the water level at dry areas to make sure the water height h + # stays positive. The system would not be stable for h set to a hard 0 due to division by h in + # the computation of velocity, e.g., (h v1) / h. Therefore, a small dry state threshold + # with a default value of 500*eps() ≈ 1e-13 in double precision, is set in the constructor above + # for the ShallowWaterEquationsWetDry and added to the initial condition if h = 0. + # This default value can be changed within the constructor call depending on the simulation setup. + H = max(H, b + equations.threshold_limiter) + return prim2cons(SVector(H, v1, v2, b), equations) +end + +initial_condition = initial_condition_parabolic_bowl + +############################################################################### +# 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) + +basis = LobattoLegendreBasis(4) + +indicator_sc = IndicatorHennemannGassnerShallowWater(equations, basis, + alpha_max = 0.6, + 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) + +############################################################################### + +coordinates_min = (-2.0, -2.0) +coordinates_max = (2.0, 2.0) + +cells_per_dimension = (150, 150) + +mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max) + +# create the semi discretization object +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 1.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 1000 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + save_analysis = false, + extra_analysis_integrals = (energy_kinetic, + energy_internal)) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 100, + save_initial_solution = true, + save_final_solution = true) + +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution) + +stage_limiter! = PositivityPreservingLimiterShallowWater(variables = (Trixi.waterheight,)) + +############################################################################### +# run the simulation + +sol = solve(ode, SSPRK43(stage_limiter!); + ode_default_options()..., callback = callbacks); + +summary_callback() # print the timer summary diff --git a/examples/structured_2d_dgsem/elixir_shallowwater_source_terms.jl b/examples/structured_2d_dgsem/elixir_shallowwater_source_terms.jl new file mode 100644 index 0000000..20656cf --- /dev/null +++ b/examples/structured_2d_dgsem/elixir_shallowwater_source_terms.jl @@ -0,0 +1,59 @@ + +using OrdinaryDiffEq +using Trixi +using TrixiShallowWater + +############################################################################### +# semidiscretization of the shallow water equations + +equations = ShallowWaterEquationsWetDry2D(gravity_constant = 9.81) + +initial_condition = initial_condition_convergence_test + +volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) +solver = DGSEM(polydeg = 3, + surface_flux = (flux_lax_friedrichs, flux_nonconservative_fjordholm_etal), + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +coordinates_min = (0.0, 0.0) +coordinates_max = (sqrt(2.0), sqrt(2.0)) + +cells_per_dimension = (8, 8) + +mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + source_terms = source_terms_convergence_test) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 1.0) +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(interval = 100, + save_initial_solution = true, + save_final_solution = true, + solution_variables = cons2prim) + +stepsize_callback = StepsizeCallback(cfl = 2.0) + +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback, + save_solution, + 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 + save_everystep = false, callback = callbacks); +summary_callback() # print the timer summary diff --git a/examples/structured_2d_dgsem/elixir_shallowwater_well_balanced.jl b/examples/structured_2d_dgsem/elixir_shallowwater_well_balanced.jl new file mode 100644 index 0000000..11c8d68 --- /dev/null +++ b/examples/structured_2d_dgsem/elixir_shallowwater_well_balanced.jl @@ -0,0 +1,124 @@ + +using Downloads: download +using OrdinaryDiffEq +using Trixi +using TrixiShallowWater + +############################################################################### +# semidiscretization of the shallow water equations with a discontinuous +# bottom topography function (set in the initial conditions) + +equations = ShallowWaterEquationsWetDry2D(gravity_constant = 9.81, H0 = 3.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::ShallowWaterEquationsWetDry2D) + # Set the background values + H = equations.H0 + v1 = 0.0 + v2 = 0.0 + # bottom topography taken from Pond.control in [HOHQMesh](https://github.com/trixi-framework/HOHQMesh) + x1, x2 = x + b = (1.5 / exp(0.5 * ((x1 - 1.0)^2 + (x2 - 1.0)^2)) + + 0.75 / exp(0.5 * ((x1 + 1.0)^2 + (x2 + 1.0)^2))) + return prim2cons(SVector(H, v1, v2, b), equations) +end + +initial_condition = initial_condition_well_balancedness + +############################################################################### +# Get the DG approximation space + +volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) +surface_flux = (FluxHydrostaticReconstruction(flux_lax_friedrichs, + hydrostatic_reconstruction_audusse_etal), + flux_nonconservative_audusse_etal) +solver = DGSEM(polydeg = 4, surface_flux = surface_flux, + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +############################################################################### +# This setup a structured periodic mesh + +coordinates_min = (-2.0, -2.0) +coordinates_max = (2.0, 2.0) + +cells_per_dimension = (4, 4) + +mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max) + +# Create the semi discretization object +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +############################################################################### +# ODE solver + +tspan = (0.0, 100.0) +ode = semidiscretize(semi, tspan) + +############################################################################### +# Workaround to set a discontinuous bottom topography for debugging and testing. + +# alternative version of the initial conditinon used to setup a truly discontinuous +# bottom topography function for this academic testcase. +# The errors from the analysis callback are not important but the error for this lake at rest test case +# `∑|H0-(h+b)|` should be around machine roundoff +# In contrast to the usual signature of initial conditions, this one get passed the +# `element_id` explicitly. In particular, this initial conditions works as intended +# only for the specific mesh loaded above! +function initial_condition_discontinuous_well_balancedness(x, t, element_id, + equations::ShallowWaterEquationsWetDry2D) + # Set the background values + H = equations.H0 + v1 = 0.0 + v2 = 0.0 + b = 0.0 + + # Setup a discontinuous bottom topography using the element id number + if element_id == 7 + b = 2.0 + 0.5 * sin(2.0 * pi * x[1]) + 0.5 * cos(2.0 * pi * x[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 + +############################################################################### +# Callbacks + +summary_callback = SummaryCallback() + +analysis_interval = 1000 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + extra_analysis_integrals = (lake_at_rest_error,)) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 1000, + save_initial_solution = true, + save_final_solution = true) + +stepsize_callback = StepsizeCallback(cfl = 3.0) + +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, + 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 + save_everystep = false, callback = callbacks); +summary_callback() # print the timer summary diff --git a/examples/structured_2d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl b/examples/structured_2d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl new file mode 100644 index 0000000..7d4d4c1 --- /dev/null +++ b/examples/structured_2d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl @@ -0,0 +1,206 @@ + +using OrdinaryDiffEq +using Trixi +using TrixiShallowWater +using Printf: @printf, @sprintf + +############################################################################### +# Semidiscretization of the shallow water equations + +equations = ShallowWaterEquationsWetDry2D(gravity_constant = 9.812) + +""" + initial_condition_well_balanced_chen_noelle(x, t, equations:: ShallowWaterEquationsWetDry2D) + +Initial condition with a complex (discontinuous) bottom topography to test the well-balanced +property for the [`hydrostatic_reconstruction_chen_noelle`](@ref) including dry areas within the +domain. 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. + +The initial condition is taken from Section 5.2 of the paper: +- Guoxian Chen and Sebastian Noelle (2017) + A new hydrostatic reconstruction scheme based on subcell reconstructions + [DOI:10.1137/15M1053074](https://dx.doi.org/10.1137/15M1053074) +""" +function initial_condition_complex_bottom_well_balanced(x, t, + equations::ShallowWaterEquationsWetDry2D) + v1 = 0 + v2 = 0 + b = sin(4 * pi * x[1]) + 3 + + if x[1] >= 0.5 + b = sin(4 * pi * x[1]) + 1 + end + + H = max(b, 2.5) + + if x[1] >= 0.5 + H = max(b, 1.5) + end + + # It is mandatory to shift the water level at dry areas to make sure the water height h + # stays positive. The system would not be stable for h set to a hard 0 due to division by h in + # the computation of velocity, e.g., (h v1) / h. Therefore, a small dry state threshold + # with a default value of 500*eps() ≈ 1e-13 in double precision, is set in the constructor above + # for the ShallowWaterEquationsWetDry and added to the initial condition if h = 0. + # This default value can be changed within the constructor call depending on the simulation setup. + H = max(H, b + equations.threshold_limiter) + return prim2cons(SVector(H, v1, v2, b), equations) +end + +initial_condition = initial_condition_complex_bottom_well_balanced + +############################################################################### +# 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) + +basis = LobattoLegendreBasis(3) + +indicator_sc = IndicatorHennemannGassnerShallowWater(equations, basis, + alpha_max = 0.5, + alpha_min = 0.001, + alpha_smooth = true, + variable = waterheight_pressure) +volume_integral = VolumeIntegralShockCapturingHG(indicator_sc; + volume_flux_dg = volume_flux, + volume_flux_fv = surface_flux) + +solver = DGSEM(basis, surface_flux, volume_integral) + +############################################################################### +# Create the StructuredMesh for the domain [0, 1]^2 + +coordinates_min = (0.0, 0.0) +coordinates_max = (1.0, 1.0) + +cells_per_dimension = (16, 16) + +mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max) + +# create the semi discretization object +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 10.0) +ode = semidiscretize(semi, tspan) + +############################################################################### +# Workaround to set a discontinuous water and bottom topography for +# debugging and testing. Essentially, this is a slight augmentation of the +# `compute_coefficients` where the `x` node value passed here is slightly +# perturbed to the left / right in order to set a true discontinuity that avoids +# the doubled value of the LGL nodes at a particular element interface. +# +# Note! The errors from the analysis callback are not important but the error +# for this lake at rest test case `∑|H0-(h+b)|` should be near machine roundoff. + +# point to the data we want to augment +u = Trixi.wrap_array(ode.u0, semi) +# reset the initial condition +for element in eachelement(semi.solver, semi.cache) + for 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) + # We know that the discontinuity is a vertical line. Slightly augment the x value by a factor + # of unit roundoff to avoid the repeted value from the LGL nodes at at interface. + if i == 1 + x_node = SVector(nextfloat(x_node[1]), x_node[2]) + elseif i == nnodes(semi.solver) + x_node = SVector(prevfloat(x_node[1]), x_node[2]) + end + u_node = initial_condition_complex_bottom_well_balanced(x_node, first(tspan), + equations) + Trixi.set_node_vars!(u, u_node, equations, semi.solver, i, j, element) + end +end + +############################################################################### +# Callbacks + +summary_callback = SummaryCallback() + +analysis_interval = 1000 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + save_analysis = false) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 1000, + save_initial_solution = true, + save_final_solution = true) + +stepsize_callback = StepsizeCallback(cfl = 1.0) + +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, + stepsize_callback) + +stage_limiter! = PositivityPreservingLimiterShallowWater(variables = (Trixi.waterheight,)) + +############################################################################### +# run the simulation + +sol = solve(ode, SSPRK43(stage_limiter!); dt = 1.0, + ode_default_options()..., callback = callbacks, adaptive = false); + +summary_callback() # print the timer summary + +############################################################################### +# Workaround to compute the well-balancedness error for this particular problem +# that has two reference water heights. One for a lake to the left of the +# discontinuous bottom topography `H0_upper = 2.5` and another for a lake to the +# right of the discontinuous bottom topography `H0_lower = 1.5`. + +# Declare a special version of the function to compute the lake-at-rest error +# OBS! The reference water height values are hardcoded for convenience. +function lake_at_rest_error_two_level(u, x, equations::ShallowWaterEquationsWetDry2D) + h, _, _, b = u + + # For well-balancedness testing with possible wet/dry regions the reference + # water height `H0` accounts for the possibility that the bottom topography + # can emerge out of the water as well as for the threshold offset to avoid + # division by a "hard" zero water heights as well. + if x[1] < 0.5 + H0_wet_dry = max(2.5, b + equations.threshold_limiter) + else + H0_wet_dry = max(1.5, b + equations.threshold_limiter) + end + + return abs(H0_wet_dry - (h + b)) +end + +# point to the data we want to analyze +u = Trixi.wrap_array(sol[end], semi) +# Perform the actual integration of the well-balancedness error over the domain +l1_well_balance_error = Trixi.integrate_via_indices(u, mesh, equations, semi.solver, + semi.cache; + normalize = true) do u, i, j, element, + equations, solver + x_node = Trixi.get_node_coords(semi.cache.elements.node_coordinates, equations, solver, + i, j, element) + # We know that the discontinuity is a vertical line. Slightly augment the x value by a factor + # of unit roundoff to avoid the repeted value from the LGL nodes at at interface. + if i == 1 + x_node = SVector(nextfloat(x_node[1]), x_node[2]) + elseif i == nnodes(semi.solver) + x_node = SVector(prevfloat(x_node[1]), x_node[2]) + end + u_local = Trixi.get_node_vars(u, equations, solver, i, j, element) + return lake_at_rest_error_two_level(u_local, x_node, equations) +end + +# report the well-balancedness lake-at-rest error to the screen +println("─"^100) +println(" Lake-at-rest error for '", Trixi.get_name(equations), "' with ", summary(solver), + " at final time " * @sprintf("%10.8e", tspan[end])) + +@printf(" %-12s:", Trixi.pretty_form_utf(lake_at_rest_error)) +@printf(" % 10.8e", l1_well_balance_error) +println() +println("─"^100) diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_beach.jl b/examples/tree_1d_dgsem/elixir_shallowwater_beach.jl new file mode 100644 index 0000000..4a87144 --- /dev/null +++ b/examples/tree_1d_dgsem/elixir_shallowwater_beach.jl @@ -0,0 +1,123 @@ + +using OrdinaryDiffEq +using Trixi +using TrixiShallowWater + +############################################################################### +# Semidiscretization of the shallow water equations + +equations = ShallowWaterEquationsWetDry1D(gravity_constant = 9.812) + +""" + initial_condition_beach(x, t, equations:: ShallowWaterEquationsWetDry1D) + +Initial condition to simulate a wave running towards a beach and crashing. Difficult test +including both wetting and drying in the domain using slip wall boundary conditions. +The bottom topography is altered to be differentiable on the domain [0,8] and +differs from the reference below. + +The water height and speed functions used here, are adapted from the initial condition +found in section 5.2 of the paper: + - Andreas Bollermann, Sebastian Noelle, Maria Lukáčová-Medvid’ová (2011) + Finite volume evolution Galerkin methods for the shallow water equations with dry beds\n + [DOI: 10.4208/cicp.220210.020710a](https://dx.doi.org/10.4208/cicp.220210.020710a) +""" +function initial_condition_beach(x, t, equations::ShallowWaterEquationsWetDry1D) + D = 1 + delta = 0.02 + gamma = sqrt((3 * delta) / (4 * D)) + x_a = sqrt((4 * D) / (3 * delta)) * acosh(sqrt(20)) + + f = D + 40 * delta * sech(gamma * (8 * x[1] - x_a))^2 + + # steep curved beach + b = 0.01 + 99 / 409600 * 4^x[1] + + if x[1] >= 6 + H = b + v = 0.0 + else + H = f + v = sqrt(equations.gravity / D) * H + end + + # It is mandatory to shift the water level at dry areas to make sure the water height h + # stays positive. The system would not be stable for h set to a hard 0 due to division by h in + # the computation of velocity, e.g., (h v) / h. Therefore, a small dry state threshold + # with a default value of 500*eps() ≈ 1e-13 in double precision, is set in the constructor above + # for the ShallowWaterEquationsWetDry and added to the initial condition if h = 0. + # This default value can be changed within the constructor call depending on the simulation setup. + H = max(H, b + equations.threshold_limiter) + return prim2cons(SVector(H, v, b), equations) +end + +initial_condition = initial_condition_beach +boundary_condition = boundary_condition_slip_wall + +############################################################################### +# Get the DG approximation space + +volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) +surface_flux = (FluxHydrostaticReconstruction(flux_hll_chen_noelle, + hydrostatic_reconstruction_chen_noelle), + flux_nonconservative_chen_noelle) + +basis = LobattoLegendreBasis(3) + +indicator_sc = IndicatorHennemannGassnerShallowWater(equations, basis, + alpha_max = 0.5, + alpha_min = 0.001, + alpha_smooth = true, + variable = waterheight_pressure) +volume_integral = VolumeIntegralShockCapturingHG(indicator_sc; + volume_flux_dg = volume_flux, + volume_flux_fv = surface_flux) + +solver = DGSEM(basis, surface_flux, volume_integral) + +############################################################################### +# Create the TreeMesh for the domain [0, 8] + +coordinates_min = 0.0 +coordinates_max = 8.0 + +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 7, + n_cells_max = 10_000, + periodicity = false) + +# create the semi discretization object +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_condition) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 10.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 1000 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + save_analysis = false, + extra_analysis_integrals = (energy_kinetic, + energy_internal)) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(dt = 0.5, + save_initial_solution = true, + save_final_solution = true) + +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution) + +stage_limiter! = PositivityPreservingLimiterShallowWater(variables = (Trixi.waterheight,)) + +############################################################################### +# run the simulation + +sol = solve(ode, SSPRK43(stage_limiter!); + ode_default_options()..., callback = callbacks); + +summary_callback() # print the timer summary diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_parabolic_bowl.jl b/examples/tree_1d_dgsem/elixir_shallowwater_parabolic_bowl.jl new file mode 100644 index 0000000..5f40bc4 --- /dev/null +++ b/examples/tree_1d_dgsem/elixir_shallowwater_parabolic_bowl.jl @@ -0,0 +1,118 @@ + +using OrdinaryDiffEq +using Trixi +using TrixiShallowWater + +############################################################################### +# Semidiscretization of the shallow water equations + +equations = ShallowWaterEquationsWetDry1D(gravity_constant = 9.81) + +""" + initial_condition_parabolic_bowl(x, t, equations:: ShallowWaterEquationsWetDry1D) + +Well-known initial condition to test the [`hydrostatic_reconstruction_chen_noelle`](@ref) and its +wet-dry mechanics. This test has analytical solutions. The initial condition is defined by the +analytical solution at time t=0. The bottom topography defines a bowl and the water level is given +by an oscillating lake. + +The original test and its analytical solution in two dimensions were first presented in +- William C. Thacker (1981) + Some exact solutions to the nonlinear shallow-water wave equations + [DOI: 10.1017/S0022112081001882](https://doi.org/10.1017/S0022112081001882). + +The particular setup below is taken from Section 6.2 of +- Niklas Wintermeyer, Andrew R. Winters, Gregor J. Gassner and Timothy Warburton (2018) + An entropy stable discontinuous Galerkin method for the shallow water equations on + curvilinear meshes with wet/dry fronts accelerated by GPUs + [DOI: 10.1016/j.jcp.2018.08.038](https://doi.org/10.1016/j.jcp.2018.08.038). +""" +function initial_condition_parabolic_bowl(x, t, equations::ShallowWaterEquationsWetDry1D) + a = 1 + h_0 = 0.1 + sigma = 0.5 + ω = sqrt(2 * equations.gravity * h_0) / a + + v = -sigma * ω * sin(ω * t) + + b = h_0 * x[1]^2 / a^2 + + H = sigma * h_0 / a^2 * (2 * x[1] * cos(ω * t) - sigma) + h_0 + + # It is mandatory to shift the water level at dry areas to make sure the water height h + # stays positive. The system would not be stable for h set to a hard 0 due to division by h in + # the computation of velocity, e.g., (h v) / h. Therefore, a small dry state threshold + # with a default value of 500*eps() ≈ 1e-13 in double precision, is set in the constructor above + # for the ShallowWaterEquationsWetDry and added to the initial condition if h = 0. + # This default value can be changed within the constructor call depending on the simulation setup. + H = max(H, b + equations.threshold_limiter) + return prim2cons(SVector(H, v, b), equations) +end + +initial_condition = initial_condition_parabolic_bowl + +############################################################################### +# Get the DG approximation space + +volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) +surface_flux = (FluxHydrostaticReconstruction(flux_hll_chen_noelle, + hydrostatic_reconstruction_chen_noelle), + flux_nonconservative_chen_noelle) + +basis = LobattoLegendreBasis(5) + +indicator_sc = IndicatorHennemannGassnerShallowWater(equations, basis, + alpha_max = 0.5, + alpha_min = 0.001, + alpha_smooth = true, + variable = waterheight_pressure) +volume_integral = VolumeIntegralShockCapturingHG(indicator_sc; + volume_flux_dg = volume_flux, + volume_flux_fv = surface_flux) + +solver = DGSEM(basis, surface_flux, volume_integral) + +############################################################################### +# Create the TreeMesh for the domain [-2, 2] + +coordinates_min = -2.0 +coordinates_max = 2.0 + +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 6, + n_cells_max = 10_000) + +# create the semi discretization object +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 10.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 1000 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + save_analysis = false, + extra_analysis_integrals = (energy_kinetic, + energy_internal)) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 1000, + save_initial_solution = true, + save_final_solution = true) + +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution) + +stage_limiter! = PositivityPreservingLimiterShallowWater(variables = (Trixi.waterheight,)) + +############################################################################### +# run the simulation + +sol = solve(ode, SSPRK43(stage_limiter!); + ode_default_options()..., callback = callbacks); + +summary_callback() # print the timer summary diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_shock_capturing.jl b/examples/tree_1d_dgsem/elixir_shallowwater_shock_capturing.jl new file mode 100644 index 0000000..bbb65d0 --- /dev/null +++ b/examples/tree_1d_dgsem/elixir_shallowwater_shock_capturing.jl @@ -0,0 +1,117 @@ + +using OrdinaryDiffEq +using Trixi +using TrixiShallowWater + +############################################################################### +# Semidiscretization of the shallow water equations + +equations = ShallowWaterEquationsWetDry1D(gravity_constant = 9.812, H0 = 1.75) + +# Initial condition with a truly discontinuous velocity and bottom topography. +# Works as intended for TreeMesh1D with `initial_refinement_level=3`. If the mesh +# refinement level is changed the initial condition below may need changed as well to +# ensure that the discontinuities lie on an element interface. +function initial_condition_stone_throw_discontinuous_bottom(x, t, + equations::ShallowWaterEquationsWetDry1D) + + # Calculate primitive variables + + # flat lake + H = equations.H0 + + # Discontinuous velocity + v = 0.0 + if x[1] >= -0.75 && x[1] <= 0.0 + v = -1.0 + elseif x[1] >= 0.0 && x[1] <= 0.75 + v = 1.0 + end + + b = (1.5 / exp(0.5 * ((x[1] - 1.0)^2)) + + 0.75 / exp(0.5 * ((x[1] + 1.0)^2))) + + # Force a discontinuous bottom topography + if x[1] >= -1.5 && x[1] <= 0.0 + b = 0.5 + end + + return prim2cons(SVector(H, v, b), equations) +end + +initial_condition = initial_condition_stone_throw_discontinuous_bottom + +boundary_condition = boundary_condition_slip_wall + +############################################################################### +# Get the DG approximation space + +volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) +surface_flux = (FluxHydrostaticReconstruction(flux_lax_friedrichs, + hydrostatic_reconstruction_audusse_etal), + flux_nonconservative_audusse_etal) +basis = LobattoLegendreBasis(4) + +indicator_sc = IndicatorHennemannGassner(equations, basis, + alpha_max = 0.5, + alpha_min = 0.001, + alpha_smooth = true, + variable = waterheight_pressure) +volume_integral = VolumeIntegralShockCapturingHG(indicator_sc; + volume_flux_dg = volume_flux, + volume_flux_fv = surface_flux) + +solver = DGSEM(basis, surface_flux, volume_integral) + +############################################################################### +# Create the TreeMesh for the domain [-3, 3] + +coordinates_min = -3.0 +coordinates_max = 3.0 +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 3, + n_cells_max = 10_000, + periodicity = false) + +# create the semi discretization object +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_condition) + +############################################################################### +# ODE solver + +tspan = (0.0, 3.0) +ode = semidiscretize(semi, tspan) + +############################################################################### +# Callbacks + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + save_analysis = false, + extra_analysis_integrals = (energy_kinetic, + energy_internal, + lake_at_rest_error)) + +# Enable in-situ visualization with a new plot generated every 50 time steps +# and we explicitly pass that the plot data will be one-dimensional +# visualization = VisualizationCallback(interval=50, plot_data_creator=PlotData1D) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 100, + save_initial_solution = true, + save_final_solution = true) + +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution)#, +# visualization) + +############################################################################### +# run the simulation + +# use a Runge-Kutta method with automatic (error based) time step size control +sol = solve(ode, RDPK3SpFSAL49(); abstol = 1.0e-7, reltol = 1.0e-7, + ode_default_options()..., callback = callbacks); +summary_callback() # print the timer summary diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl b/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl new file mode 100644 index 0000000..3f11b47 --- /dev/null +++ b/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl @@ -0,0 +1,171 @@ + +using OrdinaryDiffEq +using Trixi +using Printf: @printf, @sprintf +using TrixiShallowWater + +############################################################################### +# Semidiscretization of the shallow water equations + +equations = ShallowWaterEquationsWetDry1D(gravity_constant = 9.812) + +""" + initial_condition_complex_bottom_well_balanced(x, t, equations:: ShallowWaterEquationsWetDry1D) + +Initial condition with a complex (discontinuous) bottom topography to test the well-balanced +property for the [`hydrostatic_reconstruction_chen_noelle`](@ref) including dry areas within the +domain. 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. + +The initial condition is taken from Section 5.2 of the paper: +- Guoxian Chen and Sebastian Noelle (2017) + A new hydrostatic reconstruction scheme based on subcell reconstructions + [DOI:10.1137/15M1053074](https://dx.doi.org/10.1137/15M1053074) +""" +function initial_condition_complex_bottom_well_balanced(x, t, + equations::ShallowWaterEquationsWetDry1D) + v = 0.0 + b = sin(4 * pi * x[1]) + 3 + + if x[1] >= 0.5 + b = sin(4 * pi * x[1]) + 1 + end + + H = max(b, 2.5) + + if x[1] >= 0.5 + H = max(b, 1.5) + end + + # It is mandatory to shift the water level at dry areas to make sure the water height h + # stays positive. The system would not be stable for h set to a hard 0 due to division by h in + # the computation of velocity, e.g., (h v) / h. Therefore, a small dry state threshold + # with a default value of 500*eps() ≈ 1e-13 in double precision, is set in the constructor above + # for the ShallowWaterEquationsWetDry and added to the initial condition if h = 0. + # This default value can be changed within the constructor call depending on the simulation setup. + H = max(H, b + equations.threshold_limiter) + return prim2cons(SVector(H, v, b), equations) +end + +initial_condition = initial_condition_complex_bottom_well_balanced + +############################################################################### +# 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) + +basis = LobattoLegendreBasis(3) + +indicator_sc = IndicatorHennemannGassnerShallowWater(equations, basis, + alpha_max = 0.5, + alpha_min = 0.001, + alpha_smooth = true, + variable = waterheight_pressure) +volume_integral = VolumeIntegralShockCapturingHG(indicator_sc; + volume_flux_dg = volume_flux, + volume_flux_fv = surface_flux) + +solver = DGSEM(basis, surface_flux, volume_integral) + +############################################################################### +# Create the TreeMesh for the domain [0, 1] + +coordinates_min = 0.0 +coordinates_max = 1.0 + +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 6, + n_cells_max = 10_000) + +# create the semi discretization object +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 25.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 5000 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + save_analysis = false) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 5000, + save_initial_solution = true, + save_final_solution = true) + +stepsize_callback = StepsizeCallback(cfl = 1.5) + +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, + stepsize_callback) + +stage_limiter! = PositivityPreservingLimiterShallowWater(variables = (Trixi.waterheight,)) + +############################################################################### +# run the simulation + +sol = solve(ode, SSPRK43(stage_limiter!); dt = 1.0, + ode_default_options()..., callback = callbacks, adaptive = false); + +summary_callback() # print the timer summary + +############################################################################### +# Workaround to compute the well-balancedness error for this particular problem +# that has two reference water heights. One for a lake to the left of the +# discontinuous bottom topography `H0_upper = 2.5` and another for a lake to the +# right of the discontinuous bottom topography `H0_lower = 1.5`. + +# Declare a special version of the function to compute the lake-at-rest error +# OBS! The reference water height values are hardcoded for convenience. +function lake_at_rest_error_two_level(u, x, equations::ShallowWaterEquationsWetDry1D) + h, _, b = u + + # For well-balancedness testing with possible wet/dry regions the reference + # water height `H0` accounts for the possibility that the bottom topography + # can emerge out of the water as well as for the threshold offset to avoid + # division by a "hard" zero water heights as well. + if x[1] < 0.5 + H0_wet_dry = max(2.5, b + equations.threshold_limiter) + else + H0_wet_dry = max(1.5, b + equations.threshold_limiter) + end + + return abs(H0_wet_dry - (h + b)) +end + +# point to the data we want to analyze +u = Trixi.wrap_array(sol[end], semi) +# Perform the actual integration of the well-balancedness error over the domain +l1_well_balance_error = Trixi.integrate_via_indices(u, mesh, equations, semi.solver, + semi.cache; + normalize = true) do u, i, element, + equations, solver + x_node = Trixi.get_node_coords(semi.cache.elements.node_coordinates, equations, solver, + i, element) + # We know that the discontinuity is a vertical line. Slightly augment the x value by a factor + # of unit roundoff to avoid the repeted value from the LGL nodes at at interface. + if i == 1 + x_node = SVector(nextfloat(x_node[1])) + elseif i == nnodes(semi.solver) + x_node = SVector(prevfloat(x_node[1])) + end + u_local = Trixi.get_node_vars(u, equations, solver, i, element) + return lake_at_rest_error_two_level(u_local, x_node, equations) +end + +# report the well-balancedness lake-at-rest error to the screen +println("─"^100) +println(" Lake-at-rest error for '", Trixi.get_name(equations), "' with ", summary(solver), + " at final time " * @sprintf("%10.8e", tspan[end])) + +@printf(" %-12s:", Trixi.pretty_form_utf(lake_at_rest_error)) +@printf(" % 10.8e", l1_well_balance_error) +println() +println("─"^100) diff --git a/examples/tree_2d_dgsem/elixir_shallowwater_conical_island.jl b/examples/tree_2d_dgsem/elixir_shallowwater_conical_island.jl new file mode 100644 index 0000000..b688e59 --- /dev/null +++ b/examples/tree_2d_dgsem/elixir_shallowwater_conical_island.jl @@ -0,0 +1,116 @@ + +using OrdinaryDiffEq +using Trixi +using TrixiShallowWater + +############################################################################### +# semidiscretization of the shallow water equations + +equations = ShallowWaterEquationsWetDry2D(gravity_constant = 9.81, H0 = 1.4) + +""" + initial_condition_conical_island(x, t, equations::ShallowWaterEquationsWetDry2D) + +Initial condition for the [`ShallowWaterEquationsWetDry2D`](@ref) to test the [`hydrostatic_reconstruction_chen_noelle`](@ref) +and its handling of discontinuous water heights at the start in combination with wetting and +drying. The bottom topography is given by a conical island in the middle of the domain. Around that +island, there is a cylindrical water column at t=0 and the rest of the domain is dry. This +discontinuous water height is smoothed by a logistic function. This simulation uses a Dirichlet +boundary condition with the initial values. Due to the dry cells at the boundary, this has the +effect of an outflow which can be seen in the simulation. +""" +function initial_condition_conical_island(x, t, equations::ShallowWaterEquationsWetDry2D) + # Set the background values + + v1 = 0.0 + v2 = 0.0 + + x1, x2 = x + b = max(0.1, 1.0 - 4.0 * sqrt(x1^2 + x2^2)) + + # use a logistic function to transfer water height value smoothly + L = equations.H0 # maximum of function + x0 = 0.3 # center point of function + k = -25.0 # sharpness of transfer + + H = max(b, L / (1.0 + exp(-k * (sqrt(x1^2 + x2^2) - x0)))) + + # 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 v1) / h. Therefore, a small dry state threshold + # with a default value of 500*eps() ≈ 1e-13 in double precision, is set in the constructor above + # for the ShallowWaterEquationsWetDry and added to the initial condition if h = 0. + # This default value can be changed within the constructor call depending on the simulation setup. + H = max(H, b + equations.threshold_limiter) + return prim2cons(SVector(H, v1, v2, b), equations) +end + +initial_condition = initial_condition_conical_island +boundary_conditions = BoundaryConditionDirichlet(initial_condition) + +############################################################################### +# 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) + +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 TreeMesh and setup a mesh + +coordinates_min = (-1.0, -1.0) +coordinates_max = (1.0, 1.0) +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 4, + n_cells_max = 10_000, + periodicity = false) + +# Create the semi discretization object +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_conditions) + +############################################################################### +# ODE solver + +tspan = (0.0, 10.0) +ode = semidiscretize(semi, tspan) + +############################################################################### +# Callbacks + +summary_callback = SummaryCallback() + +analysis_interval = 1000 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 100, + save_initial_solution = true, + save_final_solution = true) + +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution) + +############################################################################### +# run the simulation + +stage_limiter! = PositivityPreservingLimiterShallowWater(variables = (Trixi.waterheight,)) + +sol = solve(ode, SSPRK43(stage_limiter!); + ode_default_options()..., callback = callbacks); + +summary_callback() # print the timer summary diff --git a/examples/tree_2d_dgsem/elixir_shallowwater_ec.jl b/examples/tree_2d_dgsem/elixir_shallowwater_ec.jl new file mode 100644 index 0000000..56f0e16 --- /dev/null +++ b/examples/tree_2d_dgsem/elixir_shallowwater_ec.jl @@ -0,0 +1,120 @@ + +using OrdinaryDiffEq +using Trixi +using TrixiShallowWater + +############################################################################### +# semidiscretization of the shallow water equations with a discontinuous +# bottom topography function + +equations = ShallowWaterEquationsWetDry2D(gravity_constant = 9.81) + +# Note, this initial condition is used to compute errors in the analysis callback but the initialization is +# overwritten by `initial_condition_ec_discontinuous_bottom` below. +initial_condition = initial_condition_weak_blast_wave + +############################################################################### +# Get the DG approximation space + +volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) +solver = DGSEM(polydeg = 4, + surface_flux = (flux_fjordholm_etal, flux_nonconservative_fjordholm_etal), + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +############################################################################### +# Get the TreeMesh and setup a periodic mesh + +coordinates_min = (-1.0, -1.0) +coordinates_max = (1.0, 1.0) +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 2, + n_cells_max = 10_000) + +# Create the semi discretization object +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +############################################################################### +# ODE solver + +tspan = (0.0, 2.0) +ode = semidiscretize(semi, tspan) + +############################################################################### +# Workaround to set a discontinuous bottom topography and initial condition for debugging and testing. + +# alternative version of the initial conditinon used to setup a truly discontinuous +# bottom topography function and initial condition for this academic testcase of entropy conservation. +# The errors from the analysis callback are not important but `∑∂S/∂U ⋅ Uₜ` 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 TreeMesh2D with initial_refinement_level=2. +function initial_condition_ec_discontinuous_bottom(x, t, element_id, + equations::ShallowWaterEquationsWetDry2D) + # Set up polar coordinates + inicenter = SVector(0.7, 0.7) + x_norm = x[1] - inicenter[1] + y_norm = x[2] - inicenter[2] + r = sqrt(x_norm^2 + y_norm^2) + phi = atan(y_norm, x_norm) + sin_phi, cos_phi = sincos(phi) + + # Set the background values + H = 4.25 + v1 = 0.0 + v2 = 0.0 + b = 0.0 + + # setup the discontinuous water height and velocities + if element_id == 10 + H = 5.0 + v1 = 0.1882 * cos_phi + v2 = 0.1882 * sin_phi + end + + # Setup a discontinuous bottom topography using the element id number + if element_id == 7 + b = 2.0 + 0.5 * sin(2.0 * pi * x[1]) + 0.5 * cos(2.0 * pi * x[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_ec_discontinuous_bottom(x_node, first(tspan), element, + equations) + Trixi.set_node_vars!(u, u_node, equations, semi.solver, i, j, element) + end +end + +############################################################################### +# Callbacks + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(dt = 0.2, + save_initial_solution = true, + save_final_solution = true) + +stepsize_callback = StepsizeCallback(cfl = 3.0) + +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, + 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 + save_everystep = false, callback = callbacks); +summary_callback() # print the timer summary diff --git a/examples/tree_2d_dgsem/elixir_shallowwater_parabolic_bowl.jl b/examples/tree_2d_dgsem/elixir_shallowwater_parabolic_bowl.jl new file mode 100644 index 0000000..1aae708 --- /dev/null +++ b/examples/tree_2d_dgsem/elixir_shallowwater_parabolic_bowl.jl @@ -0,0 +1,121 @@ + +using OrdinaryDiffEq +using Trixi +using TrixiShallowWater + +############################################################################### +# Semidiscretization of the shallow water equations + +equations = ShallowWaterEquationsWetDry2D(gravity_constant = 9.81) + +""" + initial_condition_parabolic_bowl(x, t, equations:: ShallowWaterEquationsWetDry2D) + +Well-known initial condition to test the [`hydrostatic_reconstruction_chen_noelle`](@ref) and its +wet-dry mechanics. This test has an analytical solution. The initial condition is defined by the +analytical solution at time t=0. The bottom topography defines a bowl and the water level is given +by an oscillating lake. + +The original test and its analytical solution were first presented in +- William C. Thacker (1981) + Some exact solutions to the nonlinear shallow-water wave equations + [DOI: 10.1017/S0022112081001882](https://doi.org/10.1017/S0022112081001882). + +The particular setup below is taken from Section 6.2 of +- Niklas Wintermeyer, Andrew R. Winters, Gregor J. Gassner and Timothy Warburton (2018) + An entropy stable discontinuous Galerkin method for the shallow water equations on + curvilinear meshes with wet/dry fronts accelerated by GPUs + [DOI: 10.1016/j.jcp.2018.08.038](https://doi.org/10.1016/j.jcp.2018.08.038). +""" +function initial_condition_parabolic_bowl(x, t, equations::ShallowWaterEquationsWetDry2D) + a = 1.0 + h_0 = 0.1 + sigma = 0.5 + ω = sqrt(2 * equations.gravity * h_0) / a + + v1 = -sigma * ω * sin(ω * t) + v2 = sigma * ω * cos(ω * t) + + b = h_0 * ((x[1])^2 + (x[2])^2) / a^2 + + H = sigma * h_0 / a^2 * (2 * x[1] * cos(ω * t) + 2 * x[2] * sin(ω * t) - sigma) + h_0 + + # It is mandatory to shift the water level at dry areas to make sure the water height h + # stays positive. The system would not be stable for h set to a hard 0 due to division by h in + # the computation of velocity, e.g., (h v1) / h. Therefore, a small dry state threshold + # with a default value of 500*eps() ≈ 1e-13 in double precision, is set in the constructor above + # for the ShallowWaterEquationsWetDry and added to the initial condition if h = 0. + # This default value can be changed within the constructor call depending on the simulation setup. + H = max(H, b + equations.threshold_limiter) + return prim2cons(SVector(H, v1, v2, b), equations) +end + +initial_condition = initial_condition_parabolic_bowl +############################################################################### +# 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) + +#basis = LobattoLegendreBasis(7) +basis = LobattoLegendreBasis(3) + +indicator_sc = IndicatorHennemannGassnerShallowWater(equations, basis, + alpha_max = 0.6, + alpha_min = 0.001, + alpha_smooth = true, + variable = waterheight_pressure) +volume_integral = VolumeIntegralShockCapturingHG(indicator_sc; + volume_flux_dg = volume_flux, + volume_flux_fv = surface_flux) + +solver = DGSEM(basis, surface_flux, volume_integral) + +############################################################################### +# Create the TreeMesh for the domain [-2, 2]^2 + +coordinates_min = (-2.0, -2.0) +coordinates_max = (2.0, 2.0) + +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 5, + n_cells_max = 10_000) + +# create the semi discretization object +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 1.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 1000 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + save_analysis = false, + extra_analysis_integrals = (energy_kinetic, + energy_internal)) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 100, + save_initial_solution = true, + save_final_solution = true) + +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution) + +stage_limiter! = PositivityPreservingLimiterShallowWater(variables = (Trixi.waterheight,)) + +############################################################################### +# run the simulation + +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution) + +sol = solve(ode, SSPRK43(stage_limiter!); + ode_default_options()..., callback = callbacks); + +summary_callback() # print the timer summary diff --git a/examples/tree_2d_dgsem/elixir_shallowwater_source_terms.jl b/examples/tree_2d_dgsem/elixir_shallowwater_source_terms.jl new file mode 100644 index 0000000..a159acb --- /dev/null +++ b/examples/tree_2d_dgsem/elixir_shallowwater_source_terms.jl @@ -0,0 +1,60 @@ + +using OrdinaryDiffEq +using Trixi +using TrixiShallowWater + +############################################################################### +# semidiscretization of the shallow water equations + +equations = ShallowWaterEquationsWetDry2D(gravity_constant = 9.81) + +initial_condition = initial_condition_convergence_test # MMS EOC test + +############################################################################### +# Get the DG approximation space + +volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) +solver = DGSEM(polydeg = 3, + surface_flux = (flux_lax_friedrichs, flux_nonconservative_fjordholm_etal), + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +############################################################################### +# Get the TreeMesh and setup a periodic mesh + +coordinates_min = (0.0, 0.0) +coordinates_max = (sqrt(2.0), sqrt(2.0)) +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 3, + n_cells_max = 10_000, + periodicity = true) + +# create the semi discretization object +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + source_terms = source_terms_convergence_test) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 1.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 500 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 200, + save_initial_solution = true, + save_final_solution = true) + +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution) + +############################################################################### +# run the simulation + +# use a Runge-Kutta method with automatic (error based) time step size control +sol = solve(ode, RDPK3SpFSAL49(); abstol = 1.0e-8, reltol = 1.0e-8, + ode_default_options()..., callback = callbacks); +summary_callback() # print the timer summary diff --git a/examples/tree_2d_dgsem/elixir_shallowwater_source_terms_dirichlet.jl b/examples/tree_2d_dgsem/elixir_shallowwater_source_terms_dirichlet.jl new file mode 100644 index 0000000..013895f --- /dev/null +++ b/examples/tree_2d_dgsem/elixir_shallowwater_source_terms_dirichlet.jl @@ -0,0 +1,63 @@ + +using OrdinaryDiffEq +using Trixi +using TrixiShallowWater + +############################################################################### +# Semidiscretization of the shallow water equations + +equations = ShallowWaterEquationsWetDry2D(gravity_constant = 9.81) + +initial_condition = initial_condition_convergence_test + +boundary_condition = BoundaryConditionDirichlet(initial_condition) + +############################################################################### +# Get the DG approximation space + +volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) +solver = DGSEM(polydeg = 3, + surface_flux = (flux_lax_friedrichs, flux_nonconservative_fjordholm_etal), + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +############################################################################### +# Get the TreeMesh and setup a periodic mesh + +coordinates_min = (0.0, 0.0) +coordinates_max = (sqrt(2.0), sqrt(2.0)) +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 3, + n_cells_max = 10_000, + periodicity = false) + +# create the semi discretization object +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_condition, + source_terms = source_terms_convergence_test) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 1.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 500 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 200, + save_initial_solution = true, + save_final_solution = true) + +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution) + +############################################################################### +# run the simulation + +# use a Runge-Kutta method with automatic (error based) time step size control +sol = solve(ode, RDPK3SpFSAL49(); abstol = 1.0e-8, reltol = 1.0e-8, + ode_default_options()..., callback = callbacks); +summary_callback() # print the timer summary diff --git a/examples/tree_2d_dgsem/elixir_shallowwater_wall.jl b/examples/tree_2d_dgsem/elixir_shallowwater_wall.jl new file mode 100644 index 0000000..14bb8c9 --- /dev/null +++ b/examples/tree_2d_dgsem/elixir_shallowwater_wall.jl @@ -0,0 +1,83 @@ +using OrdinaryDiffEq +using Trixi +using TrixiShallowWater + +############################################################################### +# Semidiscretization of the shallow water equations + +equations = ShallowWaterEquationsWetDry2D(gravity_constant = 9.81, H0 = 3.25) + +# An initial condition with a bottom topography and a perturbation in the waterheight to test +# boundary_condition_slip_wall +function initial_condition_perturbation(x, t, equations::ShallowWaterEquationsWetDry2D) + # Set the background values + H = equations.H0 + v1 = 0.0 + v2 = 0.0 + + # Bottom topography + b = 1.5 * exp(-0.5 * ((x[1])^2 + (x[2])^2)) + # Waterheight perturbation + H = H + 0.5 * exp(-10.0 * ((x[1])^2 + (x[2])^2)) + + return prim2cons(SVector(H, v1, v2, b), equations) +end + +initial_condition = initial_condition_perturbation + +boundary_condition = boundary_condition_slip_wall + +############################################################################### +# Get the DG approximation space + +volume_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal) +surface_flux = (flux_lax_friedrichs, flux_nonconservative_ersing_etal) +solver = DGSEM(polydeg = 3, surface_flux = surface_flux, + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +############################################################################### +# Get the TreeMesh and setup a non-periodic mesh + +coordinates_min = (-1.0, -1.0) +coordinates_max = (1.0, 1.0) +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 4, + n_cells_max = 10_000, + periodicity = false) + +# create the semi discretization object +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_condition) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 0.25) +ode = semidiscretize(semi, tspan) + +############################################################################### +# Callbacks + +summary_callback = SummaryCallback() + +analysis_interval = 1000 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 1000, + save_initial_solution = true, + save_final_solution = true) + +stepsize_callback = StepsizeCallback(cfl = 1.0) + +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, + stepsize_callback) + +############################################################################### +# 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 + save_everystep = false, callback = callbacks); +summary_callback() # print the timer summary diff --git a/examples/tree_2d_dgsem/elixir_shallowwater_well_balanced.jl b/examples/tree_2d_dgsem/elixir_shallowwater_well_balanced.jl new file mode 100644 index 0000000..54b0aa2 --- /dev/null +++ b/examples/tree_2d_dgsem/elixir_shallowwater_well_balanced.jl @@ -0,0 +1,120 @@ + +using OrdinaryDiffEq +using Trixi +using TrixiShallowWater + +############################################################################### +# semidiscretization of the shallow water equations with a discontinuous +# bottom topography function + +equations = ShallowWaterEquationsWetDry2D(gravity_constant = 9.81, H0 = 3.25) + +# 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::ShallowWaterEquationsWetDry2D) + # Set the background values + H = equations.H0 + v1 = 0.0 + v2 = 0.0 + # bottom topography taken from Pond.control in [HOHQMesh](https://github.com/trixi-framework/HOHQMesh) + x1, x2 = x + b = (1.5 / exp(0.5 * ((x1 - 1.0)^2 + (x2 - 1.0)^2)) + + 0.75 / exp(0.5 * ((x1 + 1.0)^2 + (x2 + 1.0)^2))) + return prim2cons(SVector(H, v1, v2, b), equations) +end + +initial_condition = initial_condition_well_balancedness + +############################################################################### +# Get the DG approximation space + +volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) +surface_flux = (flux_fjordholm_etal, flux_nonconservative_fjordholm_etal) +solver = DGSEM(polydeg = 4, surface_flux = surface_flux, + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +############################################################################### +# Get the TreeMesh and setup a periodic mesh + +coordinates_min = (-1.0, -1.0) +coordinates_max = (1.0, 1.0) +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 2, + n_cells_max = 10_000) + +# Create the semi discretization object +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +############################################################################### +# ODE solver + +tspan = (0.0, 100.0) +ode = semidiscretize(semi, tspan) + +############################################################################### +# Workaround to set a discontinuous bottom topography and initial condition for debugging and testing. + +# alternative version of the initial conditinon used to setup a truly discontinuous +# bottom topography function for this academic testcase of well-balancedness. +# 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 TreeMesh2D with initial_refinement_level=2. +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 + b = 0.0 + + # Setup a discontinuous bottom topography using the element id number + if element_id == 7 + b = 2.0 + 0.5 * sin(2.0 * pi * x[1]) + 0.5 * cos(2.0 * pi * x[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 + +############################################################################### +# Callbacks + +summary_callback = SummaryCallback() + +analysis_interval = 1000 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + extra_analysis_integrals = (lake_at_rest_error,)) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 1000, + save_initial_solution = true, + save_final_solution = true) + +stepsize_callback = StepsizeCallback(cfl = 3.0) + +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, + 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 + save_everystep = false, callback = callbacks); +summary_callback() # print the timer summary diff --git a/examples/tree_2d_dgsem/elixir_shallowwater_well_balanced_wall.jl b/examples/tree_2d_dgsem/elixir_shallowwater_well_balanced_wall.jl new file mode 100644 index 0000000..d92008e --- /dev/null +++ b/examples/tree_2d_dgsem/elixir_shallowwater_well_balanced_wall.jl @@ -0,0 +1,123 @@ +using OrdinaryDiffEq +using Trixi +using TrixiShallowWater + +############################################################################### +# semidiscretization of the shallow water equations with a discontinuous +# bottom topography function + +equations = ShallowWaterEquationsWetDry2D(gravity_constant = 9.81, H0 = 3.25) + +# 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::ShallowWaterEquationsWetDry2D) + # Set the background values + H = equations.H0 + v1 = 0.0 + v2 = 0.0 + # bottom topography taken from Pond.control in [HOHQMesh](https://github.com/trixi-framework/HOHQMesh) + x1, x2 = x + b = (1.5 / exp(0.5 * ((x1 - 1.0)^2 + (x2 - 1.0)^2)) + + 0.75 / exp(0.5 * ((x1 + 1.0)^2 + (x2 + 1.0)^2))) + return prim2cons(SVector(H, v1, v2, b), equations) +end + +initial_condition = initial_condition_well_balancedness + +boundary_condition = boundary_condition_slip_wall + +############################################################################### +# Get the DG approximation space + +volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) +surface_flux = (flux_fjordholm_etal, flux_nonconservative_fjordholm_etal) +solver = DGSEM(polydeg = 4, surface_flux = surface_flux, + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +############################################################################### +# Get the TreeMesh and setup a periodic mesh + +coordinates_min = (-1.0, -1.0) +coordinates_max = (1.0, 1.0) +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 2, + n_cells_max = 10_000, + periodicity = false) + +# create the semi discretization object +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_condition) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 100.0) +ode = semidiscretize(semi, tspan) + +############################################################################### +# Workaround to set a discontinuous bottom topography and initial condition for debugging and testing. + +# alternative version of the initial conditinon used to setup a truly discontinuous +# bottom topography function for this academic testcase of well-balancedness. +# 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 TreeMesh2D with initial_refinement_level=2. +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 + b = 0.0 + + # Setup a discontinuous bottom topography using the element id number + if element_id == 7 + b = 2.0 + 0.5 * sin(2.0 * pi * x[1]) + 0.5 * cos(2.0 * pi * x[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 + +############################################################################### +# Callbacks + +summary_callback = SummaryCallback() + +analysis_interval = 1000 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + extra_analysis_integrals = (lake_at_rest_error,)) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 1000, + save_initial_solution = true, + save_final_solution = true) + +stepsize_callback = StepsizeCallback(cfl = 3.0) + +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, + 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 + save_everystep = false, callback = callbacks); +summary_callback() # print the timer summary diff --git a/examples/tree_2d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl b/examples/tree_2d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl new file mode 100644 index 0000000..8135e28 --- /dev/null +++ b/examples/tree_2d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl @@ -0,0 +1,205 @@ + +using OrdinaryDiffEq +using Trixi +using TrixiShallowWater +using Printf: @printf, @sprintf + +############################################################################### +# Semidiscretization of the shallow water equations + +equations = ShallowWaterEquationsWetDry2D(gravity_constant = 9.812) + +""" + initial_condition_well_balanced_chen_noelle(x, t, equations:: ShallowWaterEquationsWetDry2D) + +Initial condition with a complex (discontinuous) bottom topography to test the well-balanced +property for the [`hydrostatic_reconstruction_chen_noelle`](@ref) including dry areas within the +domain. 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. + +The initial condition is taken from Section 5.2 of the paper: +- Guoxian Chen and Sebastian Noelle (2017) + A new hydrostatic reconstruction scheme based on subcell reconstructions + [DOI:10.1137/15M1053074](https://dx.doi.org/10.1137/15M1053074) +""" +function initial_condition_complex_bottom_well_balanced(x, t, + equations::ShallowWaterEquationsWetDry2D) + v1 = 0 + v2 = 0 + b = sin(4 * pi * x[1]) + 3 + + if x[1] >= 0.5 + b = sin(4 * pi * x[1]) + 1 + end + + H = max(b, 2.5) + if x[1] >= 0.5 + H = max(b, 1.5) + end + + # It is mandatory to shift the water level at dry areas to make sure the water height h + # stays positive. The system would not be stable for h set to a hard 0 due to division by h in + # the computation of velocity, e.g., (h v1) / h. Therefore, a small dry state threshold + # with a default value of 500*eps() ≈ 1e-13 in double precision, is set in the constructor above + # for the ShallowWaterEquationsWetDry and added to the initial condition if h = 0. + # This default value can be changed within the constructor call depending on the simulation setup. + H = max(H, b + equations.threshold_limiter) + return prim2cons(SVector(H, v1, v2, b), equations) +end + +initial_condition = initial_condition_complex_bottom_well_balanced + +############################################################################### +# 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) + +basis = LobattoLegendreBasis(3) + +indicator_sc = IndicatorHennemannGassnerShallowWater(equations, basis, + alpha_max = 0.5, + alpha_min = 0.001, + alpha_smooth = true, + variable = waterheight_pressure) +volume_integral = VolumeIntegralShockCapturingHG(indicator_sc; + volume_flux_dg = volume_flux, + volume_flux_fv = surface_flux) + +solver = DGSEM(basis, surface_flux, volume_integral) + +############################################################################### +# Create the TreeMesh for the domain [0, 1]^2 + +coordinates_min = (0.0, 0.0) +coordinates_max = (1.0, 1.0) + +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 3, + n_cells_max = 10_000) + +# create the semi discretization object +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 50.0) +ode = semidiscretize(semi, tspan) + +############################################################################### +# Workaround to set a discontinuous water and bottom topography for +# debugging and testing. Essentially, this is a slight augmentation of the +# `compute_coefficients` where the `x` node value passed here is slightly +# perturbed to the left / right in order to set a true discontinuity that avoids +# the doubled value of the LGL nodes at a particular element interface. +# +# Note! The errors from the analysis callback are not important but the error +# for this lake at rest test case `∑|H0-(h+b)|` should be near machine roundoff. + +# point to the data we want to augment +u = Trixi.wrap_array(ode.u0, semi) +# reset the initial condition +for element in eachelement(semi.solver, semi.cache) + for 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) + # We know that the discontinuity is a vertical line. Slightly augment the x value by a factor + # of unit roundoff to avoid the repeted value from the LGL nodes at at interface. + if i == 1 + x_node = SVector(nextfloat(x_node[1]), x_node[2]) + elseif i == nnodes(semi.solver) + x_node = SVector(prevfloat(x_node[1]), x_node[2]) + end + u_node = initial_condition_complex_bottom_well_balanced(x_node, first(tspan), + equations) + Trixi.set_node_vars!(u, u_node, equations, semi.solver, i, j, element) + end +end + +############################################################################### +# Callbacks + +summary_callback = SummaryCallback() + +analysis_interval = 1000 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + save_analysis = false) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 1000, + save_initial_solution = true, + save_final_solution = true) + +stepsize_callback = StepsizeCallback(cfl = 2.0) + +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, + stepsize_callback) + +stage_limiter! = PositivityPreservingLimiterShallowWater(variables = (Trixi.waterheight,)) + +############################################################################### +# run the simulation + +sol = solve(ode, SSPRK43(stage_limiter!); dt = 1.0, + ode_default_options()..., callback = callbacks, adaptive = false); + +summary_callback() # print the timer summary + +############################################################################### +# Workaround to compute the well-balancedness error for this particular problem +# that has two reference water heights. One for a lake to the left of the +# discontinuous bottom topography `H0_upper = 2.5` and another for a lake to the +# right of the discontinuous bottom topography `H0_lower = 1.5`. + +# Declare a special version of the function to compute the lake-at-rest error +# OBS! The reference water height values are hardcoded for convenience. +function lake_at_rest_error_two_level(u, x, equations::ShallowWaterEquationsWetDry2D) + h, _, _, b = u + + # For well-balancedness testing with possible wet/dry regions the reference + # water height `H0` accounts for the possibility that the bottom topography + # can emerge out of the water as well as for the threshold offset to avoid + # division by a "hard" zero water heights as well. + + if x[1] < 0.5 + H0_wet_dry = max(2.5, b + equations.threshold_limiter) + else + H0_wet_dry = max(1.5, b + equations.threshold_limiter) + end + + return abs(H0_wet_dry - (h + b)) +end + +# point to the data we want to analyze +u = Trixi.wrap_array(sol[end], semi) +# Perform the actual integration of the well-balancedness error over the domain +l1_well_balance_error = Trixi.integrate_via_indices(u, mesh, equations, semi.solver, + semi.cache; + normalize = true) do u, i, j, element, + equations, solver + x_node = Trixi.get_node_coords(semi.cache.elements.node_coordinates, equations, solver, + i, j, element) + # We know that the discontinuity is a vertical line. Slightly augment the x value by a factor + # of unit roundoff to avoid the repeted value from the LGL nodes at at interface. + if i == 1 + x_node = SVector(nextfloat(x_node[1]), x_node[2]) + elseif i == nnodes(semi.solver) + x_node = SVector(prevfloat(x_node[1]), x_node[2]) + end + u_local = Trixi.get_node_vars(u, equations, solver, i, j, element) + return lake_at_rest_error_two_level(u_local, x_node, equations) +end + +# report the well-balancedness lake-at-rest error to the screen +println("─"^100) +println(" Lake-at-rest error for '", Trixi.get_name(equations), "' with ", summary(solver), + " at final time " * @sprintf("%10.8e", tspan[end])) + +@printf(" %-12s:", Trixi.pretty_form_utf(lake_at_rest_error)) +@printf(" % 10.8e", l1_well_balance_error) +println() +println("─"^100) diff --git a/examples/unstructured_2d_dgsem/elixir_shallowwater_dirichlet.jl b/examples/unstructured_2d_dgsem/elixir_shallowwater_dirichlet.jl new file mode 100644 index 0000000..c987d98 --- /dev/null +++ b/examples/unstructured_2d_dgsem/elixir_shallowwater_dirichlet.jl @@ -0,0 +1,81 @@ + +using Downloads: download +using OrdinaryDiffEq +using Trixi +using TrixiShallowWater + +############################################################################### +# semidiscretization of the shallow water equations with a continuous +# bottom topography function (set in the initial conditions) + +equations = ShallowWaterEquationsWetDry2D(gravity_constant = 1.0, H0 = 3.0) + +# An initial condition with constant total water height and zero velocities to test well-balancedness. +function initial_condition_well_balancedness(x, t, equations::ShallowWaterEquationsWetDry2D) + # Set the background values + H = equations.H0 + v1 = 0.0 + v2 = 0.0 + # bottom topography taken from Pond.control in [HOHQMesh](https://github.com/trixi-framework/HOHQMesh) + x1, x2 = x + b = (1.5 / exp(0.5 * ((x1 - 1.0)^2 + (x2 - 1.0)^2)) + + 0.75 / exp(0.5 * ((x1 + 1.0)^2 + (x2 + 1.0)^2))) + return prim2cons(SVector(H, v1, v2, b), equations) +end + +initial_condition = initial_condition_well_balancedness + +boundary_condition_constant = BoundaryConditionDirichlet(initial_condition) +boundary_condition = Dict(:OuterCircle => boundary_condition_constant) + +############################################################################### +# Get the DG approximation space + +volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) +solver = DGSEM(polydeg = 4, surface_flux = (flux_hll, flux_nonconservative_fjordholm_etal), + 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__, "mesh_outer_circle.mesh") +isfile(default_mesh_file) || + download("https://gist.githubusercontent.com/andrewwinters5000/9beddd9cd00e2a0a15865129eeb24928/raw/be71e67fa48bc4e1e97f5f6cd77c3ed34c6ba9be/mesh_outer_circle.mesh", + default_mesh_file) +mesh_file = default_mesh_file + +mesh = UnstructuredMesh2D(mesh_file) + +# Create the semi discretization object +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_condition) + +############################################################################### +# ODE solvers, callbacks, etc. + +tspan = (0.0, 100.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 1000 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + save_analysis = true, + extra_analysis_integrals = (lake_at_rest_error,)) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 1000, + save_initial_solution = true, + save_final_solution = true) + +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution) + +############################################################################### +# run the simulation + +# use a Runge-Kutta method with automatic (error based) time step size control +sol = solve(ode, RDPK3SpFSAL49(); abstol = 1.0e-11, reltol = 1.0e-11, + ode_default_options()..., callback = callbacks); +summary_callback() # print the timer summary diff --git a/examples/unstructured_2d_dgsem/elixir_shallowwater_ec.jl b/examples/unstructured_2d_dgsem/elixir_shallowwater_ec.jl new file mode 100644 index 0000000..1415335 --- /dev/null +++ b/examples/unstructured_2d_dgsem/elixir_shallowwater_ec.jl @@ -0,0 +1,124 @@ + +using Downloads: download +using OrdinaryDiffEq +using Trixi +using TrixiShallowWater + +############################################################################### +# semidiscretization of the shallow water equations with a discontinuous +# bottom topography function + +equations = ShallowWaterEquationsWetDry2D(gravity_constant = 9.81) + +# Note, this initial condition is used to compute errors in the analysis callback but the initialization is +# overwritten by `initial_condition_ec_discontinuous_bottom` below. +initial_condition = initial_condition_weak_blast_wave + +############################################################################### +# Get the DG approximation space + +volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) +solver = DGSEM(polydeg = 6, + surface_flux = (flux_fjordholm_etal, flux_nonconservative_fjordholm_etal), + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +############################################################################### +# This setup is for the curved, split form entropy conservation testing (needs periodic BCs) + +# Get the unstructured quad mesh from a file (downloads the file if not available locally) +default_mesh_file = joinpath(@__DIR__, "mesh_alfven_wave_with_twist_and_flip.mesh") +isfile(default_mesh_file) || + download("https://gist.githubusercontent.com/andrewwinters5000/8f8cd23df27fcd494553f2a89f3c1ba4/raw/85e3c8d976bbe57ca3d559d653087b0889535295/mesh_alfven_wave_with_twist_and_flip.mesh", + default_mesh_file) +mesh_file = default_mesh_file + +mesh = UnstructuredMesh2D(mesh_file, periodicity = true) + +# Create the semi discretization object +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +############################################################################### +# ODE solver + +tspan = (0.0, 2.0) +ode = semidiscretize(semi, tspan) + +############################################################################### +# Workaround to set a discontinuous bottom topography and initial condition for debugging and testing. + +# alternative version of the initial conditinon used to setup a truly discontinuous +# bottom topography function and initial condition for this academic testcase of entropy conservation. +# The errors from the analysis callback are not important but `∑∂S/∂U ⋅ Uₜ` 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_ec_discontinuous_bottom(x, t, element_id, + equations::ShallowWaterEquationsWetDry2D) + # Set up polar coordinates + inicenter = SVector(0.7, 0.7) + x_norm = x[1] - inicenter[1] + y_norm = x[2] - inicenter[2] + r = sqrt(x_norm^2 + y_norm^2) + phi = atan(y_norm, x_norm) + sin_phi, cos_phi = sincos(phi) + + # Set the background values + H = 3.25 + v1 = 0.0 + v2 = 0.0 + b = 0.0 + + # setup the discontinuous water height and velocities + if element_id == 10 + H = 4.0 + v1 = 0.1882 * cos_phi + v2 = 0.1882 * sin_phi + end + + # Setup a discontinuous bottom topography using the element id number + if element_id == 7 + b = 2.0 + 0.5 * sin(2.0 * pi * x[1]) + 0.5 * cos(2.0 * pi * x[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_ec_discontinuous_bottom(x_node, first(tspan), element, + equations) + Trixi.set_node_vars!(u, u_node, equations, semi.solver, i, j, element) + end +end + +############################################################################### +# Callbacks + +summary_callback = SummaryCallback() + +analysis_interval = 1000 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 1000, + save_initial_solution = true, + save_final_solution = true) + +stepsize_callback = StepsizeCallback(cfl = 3.0) + +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, + 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 + save_everystep = false, callback = callbacks); +summary_callback() # print the timer summary diff --git a/examples/unstructured_2d_dgsem/elixir_shallowwater_ec_shockcapturing.jl b/examples/unstructured_2d_dgsem/elixir_shallowwater_ec_shockcapturing.jl new file mode 100644 index 0000000..2b550c6 --- /dev/null +++ b/examples/unstructured_2d_dgsem/elixir_shallowwater_ec_shockcapturing.jl @@ -0,0 +1,129 @@ + +using Downloads: download +using OrdinaryDiffEq +using Trixi +using TrixiShallowWater + +############################################################################### +# semidiscretization of the shallow water equations with a discontinuous +# bottom topography function + +equations = ShallowWaterEquationsWetDry2D(gravity_constant = 9.81) + +# Note, this initial condition is used to compute errors in the analysis callback but the initialization is +# overwritten by `initial_condition_ec_discontinuous_bottom` below. +initial_condition = initial_condition_weak_blast_wave + +############################################################################### +# Get the DG approximation space + +surface_flux = (flux_fjordholm_etal, flux_nonconservative_fjordholm_etal) +volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) +polydeg = 6 +basis = LobattoLegendreBasis(polydeg) +indicator_sc = IndicatorHennemannGassner(equations, basis, + alpha_max = 0.5, + alpha_min = 0.001, + alpha_smooth = true, + variable = waterheight_pressure) +volume_integral = VolumeIntegralShockCapturingHG(indicator_sc; + volume_flux_dg = volume_flux, + volume_flux_fv = surface_flux) +solver = DGSEM(basis, surface_flux, volume_integral) + +############################################################################### +# This setup is for the curved, split form entropy conservation testing (needs periodic BCs) + +# Get the unstructured quad mesh from a file (downloads the file if not available locally) +default_mesh_file = joinpath(@__DIR__, "mesh_alfven_wave_with_twist_and_flip.mesh") +isfile(default_mesh_file) || + download("https://gist.githubusercontent.com/andrewwinters5000/8f8cd23df27fcd494553f2a89f3c1ba4/raw/85e3c8d976bbe57ca3d559d653087b0889535295/mesh_alfven_wave_with_twist_and_flip.mesh", + default_mesh_file) +mesh_file = default_mesh_file + +mesh = UnstructuredMesh2D(mesh_file, periodicity = true) + +# Create the semi discretization object +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +############################################################################### +# ODE solver + +tspan = (0.0, 2.0) +ode = semidiscretize(semi, tspan) + +############################################################################### +# Workaround to set a discontinuous bottom topography and initial condition for debugging and testing. + +# alternative version of the initial conditinon used to setup a truly discontinuous +# bottom topography function and initial condition for this academic testcase of entropy conservation. +# The errors from the analysis callback are not important but `∑∂S/∂U ⋅ Uₜ` 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_ec_discontinuous_bottom(x, t, element_id, + equations::ShallowWaterEquationsWetDry2D) + # Set up polar coordinates + inicenter = SVector(0.7, 0.7) + x_norm = x[1] - inicenter[1] + y_norm = x[2] - inicenter[2] + r = sqrt(x_norm^2 + y_norm^2) + phi = atan(y_norm, x_norm) + sin_phi, cos_phi = sincos(phi) + + # Set the background values + H = 3.25 + v1 = 0.0 + v2 = 0.0 + b = 0.0 + + # setup the discontinuous water height and velocities + if element_id == 10 + H = 4.0 + v1 = 0.1882 * cos_phi + v2 = 0.1882 * sin_phi + end + + # Setup a discontinuous bottom topography using the element id number + if element_id == 7 + b = 2.0 + 0.5 * sin(2.0 * pi * x[1]) + 0.5 * cos(2.0 * pi * x[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_ec_discontinuous_bottom(x_node, first(tspan), element, + equations) + Trixi.set_node_vars!(u, u_node, equations, semi.solver, i, j, element) + end +end + +############################################################################### +# Callbacks + +summary_callback = SummaryCallback() + +analysis_interval = 1000 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +stepsize_callback = StepsizeCallback(cfl = 3.0) + +callbacks = CallbackSet(summary_callback, analysis_callback, alive_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 + save_everystep = false, callback = callbacks); +summary_callback() # print the timer summary diff --git a/examples/unstructured_2d_dgsem/elixir_shallowwater_source_terms.jl b/examples/unstructured_2d_dgsem/elixir_shallowwater_source_terms.jl new file mode 100644 index 0000000..95342b5 --- /dev/null +++ b/examples/unstructured_2d_dgsem/elixir_shallowwater_source_terms.jl @@ -0,0 +1,67 @@ + +using Downloads: download +using OrdinaryDiffEq +using Trixi +using TrixiShallowWater + +############################################################################### +# semidiscretization of the shallow water equations with a periodic +# bottom topography function (set in the initial conditions) + +equations = ShallowWaterEquationsWetDry2D(gravity_constant = 9.81) + +initial_condition = initial_condition_convergence_test + +############################################################################### +# Get the DG approximation space + +volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) +surface_flux = (flux_fjordholm_etal, flux_nonconservative_fjordholm_etal) +solver = DGSEM(polydeg = 6, surface_flux = surface_flux, + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +############################################################################### +# This setup is for the curved, split form convergence test on a periodic domain + +# Get the unstructured quad mesh from a file (downloads the file if not available locally) +default_mesh_file = joinpath(@__DIR__, "mesh_alfven_wave_with_twist_and_flip.mesh") +isfile(default_mesh_file) || + download("https://gist.githubusercontent.com/andrewwinters5000/8f8cd23df27fcd494553f2a89f3c1ba4/raw/85e3c8d976bbe57ca3d559d653087b0889535295/mesh_alfven_wave_with_twist_and_flip.mesh", + default_mesh_file) +mesh_file = default_mesh_file + +mesh = UnstructuredMesh2D(mesh_file, periodicity = true) + +# Create the semi discretization object +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + source_terms = source_terms_convergence_test) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 1.0) +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(interval = 100, + 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), + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); +summary_callback() # print the timer summary diff --git a/examples/unstructured_2d_dgsem/elixir_shallowwater_three_mound_dam_break.jl b/examples/unstructured_2d_dgsem/elixir_shallowwater_three_mound_dam_break.jl new file mode 100644 index 0000000..54e93b5 --- /dev/null +++ b/examples/unstructured_2d_dgsem/elixir_shallowwater_three_mound_dam_break.jl @@ -0,0 +1,139 @@ + +using Downloads: download +using OrdinaryDiffEq +using Trixi +using TrixiShallowWater + +############################################################################### +# semidiscretization of the shallow water equations + +equations = ShallowWaterEquationsWetDry2D(gravity_constant = 9.81, H0 = 1.875, + threshold_limiter = 1e-12, threshold_wet = 1e-14) + +""" + initial_condition_three_mounds(x, t, equations::ShallowWaterEquationsWetDry2D) + +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::ShallowWaterEquationsWetDry2D) + + # Set the background values + v1 = 0.0 + v2 = 0.0 + + 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 = equations.H0 # 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)))) + + # Avoid division by zero by adjusting the initial condition with a small dry state threshold + # that defaults to 500*eps() ≈ 1e-13 in double precision and is set in the constructor above + # for the ShallowWaterEquationsWetDry struct. + H = max(H, b + equations.threshold_limiter) + 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_function, + equations::ShallowWaterEquationsWetDry2D) + # Impulse and bottom from inside, height from external state + u_outer = SVector(equations.threshold_wet, u_inner[2], u_inner[3], u_inner[4]) + + # calculate the boundary flux + flux = surface_flux_function(u_inner, u_outer, normal_direction, equations) + + return 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_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) +surface_flux = (FluxHydrostaticReconstruction(flux_hll_chen_noelle, + hydrostatic_reconstruction_chen_noelle), + flux_nonconservative_chen_noelle) + +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) + +default_meshfile = joinpath(@__DIR__, "mesh_three_mound.mesh") + +isfile(default_meshfile) || + download("https://gist.githubusercontent.com/svengoldberg/c3c87fecb3fc6e46be7f0d1c7cb35f83/raw/e817ecd9e6c4686581d63c46128f9b6468d396d3/mesh_three_mound.mesh", + default_meshfile) + +meshfile = default_meshfile + +mesh = UnstructuredMesh2D(meshfile) + +# 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) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 100, + save_initial_solution = true, + save_final_solution = true) + +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution) + +############################################################################### +# run the simulation + +stage_limiter! = PositivityPreservingLimiterShallowWater(variables = (Trixi.waterheight,)) + +sol = solve(ode, SSPRK43(stage_limiter!); + ode_default_options()..., callback = callbacks); +summary_callback() # print the timer summary diff --git a/examples/unstructured_2d_dgsem/elixir_shallowwater_wall_bc_shockcapturing.jl b/examples/unstructured_2d_dgsem/elixir_shallowwater_wall_bc_shockcapturing.jl new file mode 100644 index 0000000..c7b80be --- /dev/null +++ b/examples/unstructured_2d_dgsem/elixir_shallowwater_wall_bc_shockcapturing.jl @@ -0,0 +1,100 @@ + +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_stone_throw(x, t, equations::ShallowWaterEquationsWetDry2D) + # Set up polar coordinates + inicenter = SVector(0.15, 0.15) + x_norm = x[1] - inicenter[1] + y_norm = x[2] - inicenter[2] + r = sqrt(x_norm^2 + y_norm^2) + + # Calculate primitive variables + H = equations.H0 + v1 = r < 0.6 ? 2.0 : 0.0 + v2 = r < 0.6 ? -2.0 : 0.0 + # bottom topography taken from Pond.control in [HOHQMesh](https://github.com/trixi-framework/HOHQMesh) + x1, x2 = x + b = (1.5 / exp(0.5 * ((x1 - 1.0)^2 + (x2 - 1.0)^2)) + + 0.75 / exp(0.5 * ((x1 + 1.0)^2 + (x2 + 1.0)^2))) + + return prim2cons(SVector(H, v1, v2, b), equations) +end + +initial_condition = initial_condition_stone_throw + +boundary_condition = Dict(:OuterCircle => boundary_condition_slip_wall) + +############################################################################### +# Get the DG approximation space + +surface_flux = (FluxHydrostaticReconstruction(flux_hll, + hydrostatic_reconstruction_audusse_etal), + flux_nonconservative_audusse_etal) +volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) +polydeg = 6 +basis = LobattoLegendreBasis(polydeg) +indicator_sc = IndicatorHennemannGassner(equations, basis, + alpha_max = 0.5, + alpha_min = 0.001, + alpha_smooth = true, + variable = Trixi.waterheight) +volume_integral = VolumeIntegralShockCapturingHG(indicator_sc; + volume_flux_dg = volume_flux, + volume_flux_fv = surface_flux) + +solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux, + volume_integral = volume_integral) + +############################################################################### +# Get the unstructured quad mesh from a file (downloads the file if not available locally) + +default_mesh_file = joinpath(@__DIR__, "mesh_outer_circle.mesh") +isfile(default_mesh_file) || + download("https://gist.githubusercontent.com/andrewwinters5000/9beddd9cd00e2a0a15865129eeb24928/raw/be71e67fa48bc4e1e97f5f6cd77c3ed34c6ba9be/mesh_outer_circle.mesh", + default_mesh_file) +mesh_file = default_mesh_file + +mesh = UnstructuredMesh2D(mesh_file) + +# Create the semi discretization object +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_condition) + +############################################################################### +# ODE solvers, callbacks, etc. + +tspan = (0.0, 3.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + save_analysis = false, + extra_analysis_integrals = (energy_kinetic, + energy_internal)) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 100, + save_initial_solution = true, + save_final_solution = true) + +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution) + +############################################################################### +# run the simulation + +# use a Runge-Kutta method with automatic (error based) time step size control +sol = solve(ode, RDPK3SpFSAL49(); abstol = 1.0e-8, reltol = 1.0e-8, + ode_default_options()..., callback = callbacks); +summary_callback() # print the timer summary diff --git a/examples/unstructured_2d_dgsem/elixir_shallowwater_well_balanced.jl b/examples/unstructured_2d_dgsem/elixir_shallowwater_well_balanced.jl new file mode 100644 index 0000000..e0b76c6 --- /dev/null +++ b/examples/unstructured_2d_dgsem/elixir_shallowwater_well_balanced.jl @@ -0,0 +1,124 @@ + +using Downloads: download +using OrdinaryDiffEq +using Trixi +using TrixiShallowWater + +############################################################################### +# semidiscretization of the shallow water equations with a discontinuous +# bottom topography function (set in the initial conditions) + +equations = ShallowWaterEquationsWetDry2D(gravity_constant = 9.81, H0 = 3.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::ShallowWaterEquationsWetDry2D) + # Set the background values + H = equations.H0 + v1 = 0.0 + v2 = 0.0 + # bottom topography taken from Pond.control in [HOHQMesh](https://github.com/trixi-framework/HOHQMesh) + x1, x2 = x + b = (1.5 / exp(0.5 * ((x1 - 1.0)^2 + (x2 - 1.0)^2)) + + 0.75 / exp(0.5 * ((x1 + 1.0)^2 + (x2 + 1.0)^2))) + return prim2cons(SVector(H, v1, v2, b), equations) +end + +initial_condition = initial_condition_well_balancedness + +############################################################################### +# Get the DG approximation space + +volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) +surface_flux = (flux_fjordholm_etal, flux_nonconservative_fjordholm_etal) +solver = DGSEM(polydeg = 6, 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__, "mesh_alfven_wave_with_twist_and_flip.mesh") +isfile(default_mesh_file) || + download("https://gist.githubusercontent.com/andrewwinters5000/8f8cd23df27fcd494553f2a89f3c1ba4/raw/85e3c8d976bbe57ca3d559d653087b0889535295/mesh_alfven_wave_with_twist_and_flip.mesh", + default_mesh_file) +mesh_file = default_mesh_file + +mesh = UnstructuredMesh2D(mesh_file, periodicity = true) + +# Create the semi discretization object +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +############################################################################### +# ODE solver + +tspan = (0.0, 100.0) +ode = semidiscretize(semi, tspan) + +############################################################################### +# Workaround to set a discontinuous bottom topography for debugging and testing. + +# alternative version of the initial conditinon used to setup a truly discontinuous +# bottom topography function for this academic testcase. +# The errors from the analysis callback are not important but the error for this lake at rest test case +# `∑|H0-(h+b)|` should be around machine roundoff +# In contrast to the usual signature of initial conditions, this one get passed the +# `element_id` explicitly. In particular, this initial conditions works as intended +# only for the specific mesh loaded above! +function initial_condition_discontinuous_well_balancedness(x, t, element_id, + equations::ShallowWaterEquationsWetDry2D) + # Set the background values + H = equations.H0 + v1 = 0.0 + v2 = 0.0 + b = 0.0 + + # Setup a discontinuous bottom topography using the element id number + if element_id == 7 + b = 2.0 + 0.5 * sin(2.0 * pi * x[1]) + 0.5 * cos(2.0 * pi * x[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 + +############################################################################### +# Callbacks + +summary_callback = SummaryCallback() + +analysis_interval = 1000 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + extra_analysis_integrals = (lake_at_rest_error,)) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 1000, + save_initial_solution = true, + save_final_solution = true) + +stepsize_callback = StepsizeCallback(cfl = 3.0) + +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, + 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 + save_everystep = false, callback = callbacks); +summary_callback() # print the timer summary diff --git a/src/TrixiShallowWater.jl b/src/TrixiShallowWater.jl index 7cfe6da..c065f79 100644 --- a/src/TrixiShallowWater.jl +++ b/src/TrixiShallowWater.jl @@ -1,17 +1,30 @@ module TrixiShallowWater -# We decided to import only Trixi.jl and qualify symbols explicitly with e.g. `Trixi.function_name`. + +# 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 # https://github.com/trixi-framework/TrixiShallowWater.jl/pull/10#discussion_r1433720559 -using Trixi: Trixi +using Trixi +# Import additional symbols that are not exported by Trixi.jl +using Trixi: get_node_vars, set_node_vars! using MuladdMacro: @muladd using StaticArrays: SVector using Static: True, False +using LinearAlgebra: norm include("equations/equations.jl") +include("equations/numerical_fluxes.jl") +include("callbacks_stage/callbacks_stage.jl") +include("solvers/indicators.jl") + +# Export types/functions that define the public API of TrixiShallowWater.jl +export ShallowWaterEquationsWetDry1D, ShallowWaterEquationsWetDry2D + +export hydrostatic_reconstruction_chen_noelle, flux_nonconservative_chen_noelle, + min_max_speed_chen_noelle, + flux_hll_chen_noelle -# export types/functions that define the public API of TrixiShallowWater.jl -export ShallowWaterEquationsWetDry1D -# TODO: These function are currently exported by Trixi.jl. Needs to be uncommented when removed from Trixi.jl -#export hydrostatic_reconstruction_chen_noelle, flux_nonconservative_chen_noelle, min_max_speed_chen_noelle +export PositivityPreservingLimiterShallowWater +export IndicatorHennemannGassnerShallowWater end diff --git a/src/callbacks_stage/callbacks_stage.jl b/src/callbacks_stage/callbacks_stage.jl new file mode 100644 index 0000000..ddc596a --- /dev/null +++ b/src/callbacks_stage/callbacks_stage.jl @@ -0,0 +1,9 @@ +# 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 + +include("positivity_shallow_water.jl") +end # @muladd diff --git a/src/callbacks_stage/positivity_shallow_water.jl b/src/callbacks_stage/positivity_shallow_water.jl new file mode 100644 index 0000000..17b8388 --- /dev/null +++ b/src/callbacks_stage/positivity_shallow_water.jl @@ -0,0 +1,87 @@ +# 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 + +""" + PositivityPreservingLimiterShallowWater(; variables) + +The limiter is specifically designed for the shallow water equations. +It is applied to all scalar `variables` in their given order +using the defined `threshold_limiter` from the [`ShallowWaterEquationsWetDry1D`](@ref) struct +or the [`ShallowWaterEquationsWetDry2D`](@ref) struct to determine the minimal acceptable values. +The order of the `variables` is important and might have a strong influence +on the robustness. + +As opposed to the standard version of the [`PositivityPreservingLimiterZhangShu`](@ref), +nodes with a water height below the `threshold_limiter` are treated in a special way. +To avoid numerical problems caused by velocities close to zero, +the velocity is cut off, such that the node can be identified as "dry". The special feature of the +`ShallowWaterEquationsWetDry` used here is that the bottom topography is stored as an additional +quantity in the solution vector `u`. However, the value of the bottom topography +should not be changed. That is why, it is not limited. + +After the limiting process is applied to all degrees of freedom, for safety reasons, +the `threshold_limiter` is applied again on all the DG nodes in order to avoid water height below. +In the case where the cell mean value is below the threshold before applying the limiter, +there could still be dry nodes afterwards due to the logic of the limiter. + +This fully-discrete positivity-preserving limiter is based on the work of +- Zhang, Shu (2011) + Maximum-principle-satisfying and positivity-preserving high-order schemes + for conservation laws: survey and new developments + [doi: 10.1098/rspa.2011.0153](https://doi.org/10.1098/rspa.2011.0153) +""" +struct PositivityPreservingLimiterShallowWater{N, Variables <: NTuple{N, Any}} + variables::Variables +end + +function PositivityPreservingLimiterShallowWater(; variables) + PositivityPreservingLimiterShallowWater(variables) +end + +function (limiter!::PositivityPreservingLimiterShallowWater)(u_ode, integrator, + semi::Trixi.AbstractSemidiscretization, + t) + u = Trixi.wrap_array(u_ode, semi) + Trixi.@trixi_timeit Trixi.timer() "positivity-preserving limiter" limiter_shallow_water!(u, + limiter!.variables, + Trixi.mesh_equations_solver_cache(semi)...) +end + +# Iterate over tuples in a type-stable way using "lispy tuple programming", +# similar to https://stackoverflow.com/a/55849398: +# Iterating over tuples of different functions isn't type-stable in general +# but accessing the first element of a tuple is type-stable. Hence, it's good +# to process one element at a time and replace iteration by recursion here. +# Note that you shouldn't use this with too many elements per tuple since the +# compile times can increase otherwise - but a handful of elements per tuple +# is definitely fine. +function limiter_shallow_water!(u, variables::NTuple{N, Any}, + mesh, + equations::Union{ShallowWaterEquationsWetDry1D, + ShallowWaterEquationsWetDry2D}, + solver, cache) where {N} + variable = first(variables) + remaining_variables = Base.tail(variables) + + limiter_shallow_water!(u, equations.threshold_limiter, variable, mesh, equations, + solver, cache) + limiter_shallow_water!(u, remaining_variables, mesh, equations, solver, cache) + return nothing +end + +# terminate the type-stable iteration over tuples +function limiter_shallow_water!(u, variables::Tuple{}, + mesh, + equations::Union{ShallowWaterEquationsWetDry1D, + ShallowWaterEquationsWetDry2D}, + solver, cache) + nothing +end + +include("positivity_shallow_water_dg1d.jl") +include("positivity_shallow_water_dg2d.jl") +end # @muladd diff --git a/src/callbacks_stage/positivity_shallow_water_dg1d.jl b/src/callbacks_stage/positivity_shallow_water_dg1d.jl new file mode 100644 index 0000000..3494662 --- /dev/null +++ b/src/callbacks_stage/positivity_shallow_water_dg1d.jl @@ -0,0 +1,87 @@ +# 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 limiter_shallow_water!(u, threshold::Real, variable, + mesh::Trixi.AbstractMesh{1}, + equations::ShallowWaterEquationsWetDry1D, + dg::DGSEM, cache) + @unpack weights = dg.basis + + Trixi.@threaded for element in eachelement(dg, cache) + # determine minimum value + value_min = typemax(eltype(u)) + for i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, element) + value_min = min(value_min, variable(u_node, equations)) + end + + # detect if limiting is necessary + value_min < threshold || continue + + # compute mean value + u_mean = zero(get_node_vars(u, equations, dg, 1, element)) + for i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, element) + u_mean += u_node * weights[i] + end + # note that the reference element is [-1,1]^ndims(dg), thus the weights sum to 2 + u_mean = u_mean / 2^ndims(mesh) + + # We compute the value directly with the mean values, as we assume that + # Jensen's inequality holds (e.g. pressure for compressible Euler equations). + value_mean = variable(u_mean, equations) + theta = (value_mean - threshold) / (value_mean - value_min) + for i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, element) + + # Cut off velocity in case that the waterheight is smaller than the threshold + + h_node, h_v_node, b_node = u_node + h_mean, h_v_mean, _ = u_mean # b_mean is not used as b_node must not be overwritten + + # Set them both to zero to apply linear combination correctly + if h_node <= threshold + h_v_node = zero(eltype(u)) + h_v_mean = zero(eltype(u)) + end + + u_node = SVector(h_node, h_v_node, b_node) + u_mean = SVector(h_mean, h_v_mean, b_node) + + # When velocity is cut off, the only averaged value is the waterheight, + # because the velocity is set to zero and this value is passed. + # Otherwise, the velocity is averaged, as well. + # Note that the auxiliary bottom topography variable `b` is never limited. + set_node_vars!(u, theta * u_node + (1 - theta) * u_mean, + equations, dg, i, element) + end + end + + # "Safety" application of the wet/dry thresholds over all the DG nodes + # on the current `element` after the limiting above in order to avoid dry nodes. + # If the value_mean < threshold before applying limiter, there + # could still be dry nodes afterwards due to logic of the limiting + Trixi.@threaded for element in eachelement(dg, cache) + for i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, element) + + h, hv, b = u_node + + if h <= threshold + h = threshold + hv = zero(eltype(u)) + end + + u_node = SVector(h, hv, b) + + set_node_vars!(u, u_node, equations, dg, i, element) + end + end + + return nothing +end +end # @muladd diff --git a/src/callbacks_stage/positivity_shallow_water_dg2d.jl b/src/callbacks_stage/positivity_shallow_water_dg2d.jl new file mode 100644 index 0000000..812f453 --- /dev/null +++ b/src/callbacks_stage/positivity_shallow_water_dg2d.jl @@ -0,0 +1,89 @@ +# 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 limiter_shallow_water!(u, threshold::Real, variable, + mesh::Trixi.AbstractMesh{2}, + equations::ShallowWaterEquationsWetDry2D, dg::DGSEM, + cache) + @unpack weights = dg.basis + + Trixi.@threaded for element in eachelement(dg, cache) + # determine minimum value + value_min = typemax(eltype(u)) + for j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, j, element) + value_min = min(value_min, variable(u_node, equations)) + end + + # detect if limiting is necessary + value_min < threshold || continue + + # compute mean value + u_mean = zero(get_node_vars(u, equations, dg, 1, 1, element)) + for j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, j, element) + u_mean += u_node * weights[i] * weights[j] + end + # note that the reference element is [-1,1]^ndims(dg), thus the weights sum to 2 + u_mean = u_mean / 2^ndims(mesh) + + # We compute the value directly with the mean values, as we assume that + # Jensen's inequality holds (e.g. pressure for compressible Euler equations). + value_mean = variable(u_mean, equations) + theta = (value_mean - threshold) / (value_mean - value_min) + for j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, j, element) + + # Cut off velocity in case that the water height is smaller than the threshold + + h_node, h_v1_node, h_v2_node, b_node = u_node + h_mean, h_v1_mean, h_v2_mean, _ = u_mean # b_mean is not used as it must not be overwritten + + if h_node <= threshold + h_v1_node = zero(eltype(u)) + h_v2_node = zero(eltype(u)) + h_v1_mean = zero(eltype(u)) + h_v2_mean = zero(eltype(u)) + end + + u_node = SVector(h_node, h_v1_node, h_v2_node, b_node) + u_mean = SVector(h_mean, h_v1_mean, h_v2_mean, b_node) + + # When velocities are cut off, the only averaged value is the water height, + # because the velocities are set to zero and this value is passed. + # Otherwise, the velocities are averaged, as well. + # Note that the auxiliary bottom topography variable `b` is never limited. + set_node_vars!(u, theta * u_node + (1 - theta) * u_mean, + equations, dg, i, j, element) + end + end + + # "Safety" application of the wet/dry thresholds over all the DG nodes + # on the current `element` after the limiting above in order to avoid dry nodes. + # If the value_mean < threshold before applying limiter, there + # could still be dry nodes afterwards due to logic of the limiting + Trixi.@threaded for element in eachelement(dg, cache) + for j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, j, element) + + h, h_v1, h_v2, b = u_node + + if h <= threshold + h = threshold + h_v1 = zero(eltype(u)) + h_v2 = zero(eltype(u)) + end + + u_node = SVector(h, h_v1, h_v2, b) + + set_node_vars!(u, u_node, equations, dg, i, j, element) + end + end + + return nothing +end +end # @muladd diff --git a/src/equations/equations.jl b/src/equations/equations.jl index 4d31451..9fb8148 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -9,4 +9,5 @@ # Include files with actual implementations for different systems of equations. include("shallow_water_wet_dry_1d.jl") +include("shallow_water_wet_dry_2d.jl") end # @muladd diff --git a/src/equations/numerical_fluxes.jl b/src/equations/numerical_fluxes.jl new file mode 100644 index 0000000..098bc69 --- /dev/null +++ b/src/equations/numerical_fluxes.jl @@ -0,0 +1,31 @@ +# 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 + +# This file contains general numerical fluxes that are not specific to certain equations + +# An empty version of the `min_max_speed_chen_noelle` function is declared here +# in order to create a dimension agnostic version of `flux_hll_chen_noelle`. +# The full description of this wave speed estimate can be found in the docstrings +# for `min_max_speed_chen_noelle` in `shallow_water_wet_dry_1d.jl` or `shallow_water_wet_dry_2d.jl`. + +function min_max_speed_chen_noelle end + +""" + flux_hll_chen_noelle = FluxHLL(min_max_speed_chen_noelle) + +An instance of [`FluxHLL`](@ref) specific to the shallow water equations that +uses the wave speed estimates from [`min_max_speed_chen_noelle`](@ref). +This HLL flux is guaranteed to have zero numerical mass flux out of a "dry" element, +maintain positivity of the water height, and satisfy an entropy inequality. + +For complete details see Section 2.4 of the following reference +- Guoxian Chen and Sebastian Noelle (2017) + A new hydrostatic reconstruction scheme based on subcell reconstructions + [DOI: 10.1137/15M1053074](https://doi.org/10.1137/15M1053074) +""" +const flux_hll_chen_noelle = FluxHLL(min_max_speed_chen_noelle) +end diff --git a/src/equations/shallow_water_wet_dry_1d.jl b/src/equations/shallow_water_wet_dry_1d.jl index 03cb118..095041d 100644 --- a/src/equations/shallow_water_wet_dry_1d.jl +++ b/src/equations/shallow_water_wet_dry_1d.jl @@ -62,6 +62,8 @@ struct ShallowWaterEquationsWetDry1D{RealT <: Real} <: # before calculating the numerical flux. # Default is 5*eps() which in double precision is ≈1e-15. threshold_wet::RealT + # Standard shallow water equations for dispatch on Trixi.jl functions + basic_swe::ShallowWaterEquations1D{RealT} end # Allow for flexibility to set the gravitational constant within an elixir depending on the @@ -79,22 +81,25 @@ function ShallowWaterEquationsWetDry1D(; gravity_constant, H0 = zero(gravity_con if threshold_wet === nothing threshold_wet = 5 * eps(T) end + # Construct the standard SWE for dispatch. Even though the `basic_swe` already store the + # gravity constant and the total water height, we store an extra copy in + # `ShallowWaterEquationsWetDry1D` for convenience. + basic_swe = ShallowWaterEquations1D(gravity_constant = gravity_constant, H0 = H0) + ShallowWaterEquationsWetDry1D(gravity_constant, H0, threshold_limiter, - threshold_wet) + threshold_wet, basic_swe) end Trixi.have_nonconservative_terms(::ShallowWaterEquationsWetDry1D) = True() -function Trixi.varnames(::typeof(Trixi.cons2cons), ::ShallowWaterEquationsWetDry1D) +function Trixi.varnames(::typeof(cons2cons), ::ShallowWaterEquationsWetDry1D) ("h", "h_v", "b") end # Note, we use the total water height, H = h + b, as the first primitive variable for easier # visualization and setting initial conditions -function Trixi.varnames(::typeof(Trixi.cons2prim), ::ShallowWaterEquationsWetDry1D) +function Trixi.varnames(::typeof(cons2prim), ::ShallowWaterEquationsWetDry1D) ("H", "v", "b") end -# TODO: Remove thresholds from ShallowWaterEquations1D after they have been moved - # This equation set extends the basic ShallowWaterEquations1D 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 @@ -112,7 +117,7 @@ A smooth initial condition used for convergence tests in combination with function Trixi.initial_condition_convergence_test(x, t, equations::ShallowWaterEquationsWetDry1D) return Trixi.initial_condition_convergence_test(x, t, - Trixi.ShallowWaterEquations1D(gravity_constant = equations.gravity)) + equations.basic_swe) end """ @@ -130,7 +135,7 @@ as defined in [`initial_condition_convergence_test`](@ref). @inline function Trixi.source_terms_convergence_test(u, x, t, equations::ShallowWaterEquationsWetDry1D) return Trixi.source_terms_convergence_test(u, x, t, - Trixi.ShallowWaterEquations1D(gravity_constant = equations.gravity)) + equations.basic_swe) end """ @@ -142,7 +147,7 @@ Note for the shallow water equations to the total energy acts as a mathematical function Trixi.initial_condition_weak_blast_wave(x, t, equations::ShallowWaterEquationsWetDry1D) return Trixi.initial_condition_weak_blast_wave(x, t, - Trixi.ShallowWaterEquations1D(gravity_constant = equations.gravity)) + equations.basic_swe) end """ @@ -164,10 +169,22 @@ For details see Section 9.2.5 of the book: x, t, surface_flux_function, equations::ShallowWaterEquationsWetDry1D) - return Trixi.boundary_condition_slip_wall(u_inner, orientation_or_normal, direction, - x, t, - surface_flux_function, - Trixi.ShallowWaterEquations1D(gravity_constant = equations.gravity)) + # This can not be dispatched, when Flux Hydrostactic reconstruction is used + # create the "external" boundary solution state + u_boundary = SVector(u_inner[1], + -u_inner[2], + u_inner[3]) + + # calculate the boundary flux + if iseven(direction) # u_inner is "left" of boundary, u_boundary is "right" of boundary + flux = surface_flux_function(u_inner, u_boundary, orientation_or_normal, + equations) + else # u_boundary is "left" of boundary, u_inner is "right" of boundary + flux = surface_flux_function(u_boundary, u_inner, orientation_or_normal, + equations) + end + + return flux end # Calculate 1D flux for a single point @@ -175,7 +192,7 @@ end @inline function Trixi.flux(u, orientation::Integer, equations::ShallowWaterEquationsWetDry1D) return Trixi.flux(u, orientation, - Trixi.ShallowWaterEquations1D(gravity_constant = equations.gravity)) + equations.basic_swe) end """ @@ -195,7 +212,7 @@ Further details are available in the paper:#include("numerical_fluxes.jl") orientation::Integer, equations::ShallowWaterEquationsWetDry1D) return Trixi.flux_nonconservative_wintermeyer_etal(u_ll, u_rr, orientation, - Trixi.ShallowWaterEquations1D(gravity_constant = equations.gravity)) + equations.basic_swe) end """ @@ -225,7 +242,7 @@ and for curvilinear 2D case in the paper: orientation::Integer, equations::ShallowWaterEquationsWetDry1D) return Trixi.flux_nonconservative_fjordholm_etal(u_ll, u_rr, orientation, - Trixi.ShallowWaterEquations1D(gravity_constant = equations.gravity)) + equations.basic_swe) end """ @@ -251,56 +268,55 @@ Further details on the hydrostatic reconstruction and its motivation can be foun equations::ShallowWaterEquationsWetDry1D) return Trixi.flux_nonconservative_audusse_etal(u_ll, u_rr, orientation, - Trixi.ShallowWaterEquations1D(gravity_constant = equations.gravity)) + equations.basic_swe) end -# TODO: This function is currently exported by Trixi.jl. Needs to be uncommented when removed from Trixi.jl -# """ -# flux_nonconservative_chen_noelle(u_ll, u_rr, -# orientation::Integer, -# equations::ShallowWaterEquationsWetDry1D) +""" + flux_nonconservative_chen_noelle(u_ll, u_rr, + orientation::Integer, + equations::ShallowWaterEquationsWetDry1D) -# Non-symmetric two-point surface flux that discretizes the nonconservative (source) term. -# The discretization uses the `hydrostatic_reconstruction_chen_noelle` on the conservative -# variables. +Non-symmetric two-point surface flux that discretizes the nonconservative (source) term. +The discretization uses the `hydrostatic_reconstruction_chen_noelle` on the conservative +variables. -# Should be used together with [`FluxHydrostaticReconstruction`](@ref) and -# [`hydrostatic_reconstruction_chen_noelle`](@ref) in the surface flux to ensure consistency. +Should be used together with [`FluxHydrostaticReconstruction`](@ref) and +[`hydrostatic_reconstruction_chen_noelle`](@ref) in the surface flux to ensure consistency. -# Further details on the hydrostatic reconstruction and its motivation can be found in -# - Guoxian Chen and Sebastian Noelle (2017) -# A new hydrostatic reconstruction scheme based on subcell reconstructions -# [DOI:10.1137/15M1053074](https://dx.doi.org/10.1137/15M1053074) -# """ -# @inline function flux_nonconservative_chen_noelle(u_ll, u_rr, -# orientation::Integer, -# equations::ShallowWaterEquationsWetDry1D) +Further details on the hydrostatic reconstruction and its motivation can be found in +- Guoxian Chen and Sebastian Noelle (2017) + A new hydrostatic reconstruction scheme based on subcell reconstructions + [DOI:10.1137/15M1053074](https://dx.doi.org/10.1137/15M1053074) +""" +@inline function flux_nonconservative_chen_noelle(u_ll, u_rr, + orientation::Integer, + equations::ShallowWaterEquationsWetDry1D) -# # Pull the water height and bottom topography on the left -# h_ll, _, b_ll = u_ll -# h_rr, _, b_rr = u_rr + # Pull the water height and bottom topography on the left + h_ll, _, b_ll = u_ll + h_rr, _, b_rr = u_rr -# H_ll = h_ll + b_ll -# H_rr = h_rr + b_rr + H_ll = h_ll + b_ll + H_rr = h_rr + b_rr -# b_star = min(max(b_ll, b_rr), min(H_ll, H_rr)) + b_star = min(max(b_ll, b_rr), min(H_ll, H_rr)) -# # Create the hydrostatic reconstruction for the left solution state -# u_ll_star, _ = hydrostatic_reconstruction_chen_noelle(u_ll, u_rr, equations) + # Create the hydrostatic reconstruction for the left solution state + u_ll_star, _ = hydrostatic_reconstruction_chen_noelle(u_ll, u_rr, equations) -# # Copy the reconstructed water height for easier to read code -# h_ll_star = u_ll_star[1] + # Copy the reconstructed water height for easier to read code + h_ll_star = u_ll_star[1] -# z = zero(eltype(u_ll)) -# # Includes two parts: -# # (i) Diagonal (consistent) term from the volume flux that uses `b_ll` to avoid -# # cross-averaging across a discontinuous bottom topography -# # (ii) True surface part that uses `h_ll` and `h_ll_star` to handle discontinuous bathymetry -# return SVector(z, -# equations.gravity * h_ll * b_ll - -# equations.gravity * (h_ll_star + h_ll) * (b_ll - b_star), -# z) -# end + z = zero(eltype(u_ll)) + # Includes two parts: + # (i) Diagonal (consistent) term from the volume flux that uses `b_ll` to avoid + # cross-averaging across a discontinuous bottom topography + # (ii) True surface part that uses `h_ll` and `h_ll_star` to handle discontinuous bathymetry + return SVector(z, + equations.gravity * h_ll * b_ll - + equations.gravity * (h_ll_star + h_ll) * (b_ll - b_star), + z) +end """ flux_nonconservative_ersing_etal(u_ll, u_rr, orientation::Integer, @@ -323,7 +339,7 @@ For further details see: orientation::Integer, equations::ShallowWaterEquationsWetDry1D) return Trixi.flux_nonconservative_ersing_etal(u_ll, u_rr, orientation, - Trixi.ShallowWaterEquations1D(gravity_constant = equations.gravity)) + equations.basic_swe) end """ @@ -342,7 +358,7 @@ Details are available in Eq. (4.1) in the paper: @inline function Trixi.flux_fjordholm_etal(u_ll, u_rr, orientation::Integer, equations::ShallowWaterEquationsWetDry1D) return Trixi.flux_fjordholm_etal(u_ll, u_rr, orientation, - Trixi.ShallowWaterEquations1D(gravity_constant = equations.gravity)) + equations.basic_swe) end """ @@ -362,7 +378,7 @@ Further details are available in Theorem 1 of the paper: @inline function Trixi.flux_wintermeyer_etal(u_ll, u_rr, orientation::Integer, equations::ShallowWaterEquationsWetDry1D) return Trixi.flux_wintermeyer_etal(u_ll, u_rr, orientation, - Trixi.ShallowWaterEquations1D(gravity_constant = equations.gravity)) + equations.basic_swe) end """ @@ -382,134 +398,168 @@ Further details on this hydrostatic reconstruction and its motivation can be fou @inline function Trixi.hydrostatic_reconstruction_audusse_etal(u_ll, u_rr, equations::ShallowWaterEquationsWetDry1D) return Trixi.hydrostatic_reconstruction_audusse_etal(u_ll, u_rr, - Trixi.ShallowWaterEquations1D(gravity_constant = equations.gravity)) -end - -# TODO: This function is currently exported by Trixi.jl. Needs to be uncommented when removed from Trixi.jl -# """ -# hydrostatic_reconstruction_chen_noelle(u_ll, u_rr, orientation::Integer, -# equations::ShallowWaterEquationsWetDry1D) - -# A particular type of hydrostatic reconstruction of the water height to guarantee well-balancedness -# for a general bottom topography of the [`ShallowWaterEquationsWetDry1D`](@ref). The reconstructed solution states -# `u_ll_star` and `u_rr_star` variables are used to evaluate the surface numerical flux at the interface. -# The key idea is a linear reconstruction of the bottom and water height at the interfaces using subcells. -# Use in combination with the generic numerical flux routine [`FluxHydrostaticReconstruction`](@ref). - -# Further details on this hydrostatic reconstruction and its motivation can be found in -# - Guoxian Chen and Sebastian Noelle (2017) -# A new hydrostatic reconstruction scheme based on subcell reconstructions -# [DOI:10.1137/15M1053074](https://dx.doi.org/10.1137/15M1053074) -# """ -# @inline function hydrostatic_reconstruction_chen_noelle(u_ll, u_rr, -# equations::ShallowWaterEquationsWetDry1D) -# # Unpack left and right water heights and bottom topographies -# h_ll, _, b_ll = u_ll -# h_rr, _, b_rr = u_rr - -# # Get the velocities on either side -# v_ll = velocity(u_ll, equations) -# v_rr = velocity(u_rr, equations) - -# H_ll = b_ll + h_ll -# H_rr = b_rr + h_rr - -# b_star = min(max(b_ll, b_rr), min(H_ll, H_rr)) - -# # Compute the reconstructed water heights -# h_ll_star = min(H_ll - b_star, h_ll) -# h_rr_star = min(H_rr - b_star, h_rr) - -# # Set the water height to be at least the value stored in the variable threshold after -# # the hydrostatic reconstruction is applied and before the numerical flux is calculated -# # to avoid numerical problem with arbitrary small values. Interfaces with a water height -# # lower or equal to the threshold can be declared as dry. -# # The default value for `threshold_wet` is ≈ 5*eps(), or 1e-15 in double precision, is set -# # in the `ShallowWaterEquationsWetDry1D` struct. This threshold value can be changed in the constructor -# # call of this equation struct in an elixir. -# threshold = equations.threshold_wet - -# if (h_ll_star <= threshold) -# h_ll_star = threshold -# v_ll = zero(v_ll) -# end - -# if (h_rr_star <= threshold) -# h_rr_star = threshold -# v_rr = zero(v_rr) -# end - -# # Create the conservative variables using the reconstruted water heights -# u_ll_star = SVector(h_ll_star, h_ll_star * v_ll, b_ll) -# u_rr_star = SVector(h_rr_star, h_rr_star * v_rr, b_rr) - -# return u_ll_star, u_rr_star -# end + equations.basic_swe) +end + +""" + hydrostatic_reconstruction_chen_noelle(u_ll, u_rr, orientation::Integer, + equations::ShallowWaterEquationsWetDry1D) + +A particular type of hydrostatic reconstruction of the water height to guarantee well-balancedness +for a general bottom topography of the [`ShallowWaterEquationsWetDry1D`](@ref). The reconstructed solution states +`u_ll_star` and `u_rr_star` variables are used to evaluate the surface numerical flux at the interface. +The key idea is a linear reconstruction of the bottom and water height at the interfaces using subcells. +Use in combination with the generic numerical flux routine [`FluxHydrostaticReconstruction`](@ref). + +Further details on this hydrostatic reconstruction and its motivation can be found in +- Guoxian Chen and Sebastian Noelle (2017) + A new hydrostatic reconstruction scheme based on subcell reconstructions + [DOI:10.1137/15M1053074](https://dx.doi.org/10.1137/15M1053074) +""" +@inline function hydrostatic_reconstruction_chen_noelle(u_ll, u_rr, + equations::ShallowWaterEquationsWetDry1D) + # Unpack left and right water heights and bottom topographies + h_ll, _, b_ll = u_ll + h_rr, _, b_rr = u_rr + + # Get the velocities on either side + v_ll = velocity(u_ll, equations) + v_rr = velocity(u_rr, equations) + + H_ll = b_ll + h_ll + H_rr = b_rr + h_rr + + b_star = min(max(b_ll, b_rr), min(H_ll, H_rr)) + + # Compute the reconstructed water heights + h_ll_star = min(H_ll - b_star, h_ll) + h_rr_star = min(H_rr - b_star, h_rr) + + # Set the water height to be at least the value stored in the variable threshold after + # the hydrostatic reconstruction is applied and before the numerical flux is calculated + # to avoid numerical problem with arbitrary small values. Interfaces with a water height + # lower or equal to the threshold can be declared as dry. + # The default value for `threshold_wet` is ≈ 5*eps(), or 1e-15 in double precision, is set + # in the `ShallowWaterEquationsWetDry1D` struct. This threshold value can be changed in the constructor + # call of this equation struct in an elixir. + threshold = equations.threshold_wet + + if (h_ll_star <= threshold) + h_ll_star = threshold + v_ll = zero(v_ll) + end + + if (h_rr_star <= threshold) + h_rr_star = threshold + v_rr = zero(v_rr) + end + + # Create the conservative variables using the reconstruted water heights + u_ll_star = SVector(h_ll_star, h_ll_star * v_ll, b_ll) + u_rr_star = SVector(h_rr_star, h_rr_star * v_rr, b_rr) + + return u_ll_star, u_rr_star +end # Specialized `DissipationLocalLaxFriedrichs` to avoid spurious dissipation in the bottom topography -@inline function (dissipation::Trixi.DissipationLocalLaxFriedrichs)(u_ll, u_rr, - orientation_or_normal_direction, - equations::ShallowWaterEquationsWetDry1D) - return (dissipation::Trixi.DissipationLocalLaxFriedrichs)(u_ll, u_rr, +@inline function (dissipation::DissipationLocalLaxFriedrichs)(u_ll, u_rr, orientation_or_normal_direction, - Trixi.ShallowWaterEquations1D(gravity_constant = equations.gravity)) + equations::ShallowWaterEquationsWetDry1D) + return (dissipation::DissipationLocalLaxFriedrichs)(u_ll, u_rr, + orientation_or_normal_direction, + equations.basic_swe) end # Specialized `FluxHLL` to avoid spurious dissipation in the bottom topography -@inline function (numflux::Trixi.FluxHLL)(u_ll, u_rr, orientation_or_normal_direction, - equations::ShallowWaterEquationsWetDry1D) - return (numflux::Trixi.FluxHLL)(u_ll, u_rr, orientation_or_normal_direction, - Trixi.ShallowWaterEquations1D(gravity_constant = equations.gravity)) +@inline function (numflux::FluxHLL)(u_ll, u_rr, orientation_or_normal_direction, + equations::ShallowWaterEquationsWetDry1D) + λ_min, λ_max = numflux.min_max_speed(u_ll, u_rr, orientation_or_normal_direction, + equations) + + if λ_min >= 0 && λ_max >= 0 + return flux(u_ll, orientation_or_normal_direction, equations) + elseif λ_max <= 0 && λ_min <= 0 + return flux(u_rr, orientation_or_normal_direction, equations) + else + f_ll = flux(u_ll, orientation_or_normal_direction, equations) + f_rr = flux(u_rr, orientation_or_normal_direction, equations) + inv_λ_max_minus_λ_min = inv(λ_max - λ_min) + factor_ll = λ_max * inv_λ_max_minus_λ_min + factor_rr = λ_min * inv_λ_max_minus_λ_min + factor_diss = λ_min * λ_max * inv_λ_max_minus_λ_min + diss = u_rr - u_ll + return factor_ll * f_ll - factor_rr * f_rr + + factor_diss * SVector(diss[1], diss[2], zero(eltype(u_ll))) + end end -# TODO: This function is currently exported by Trixi.jl. Needs to be uncommented when removed from Trixi.jl -# """ -# min_max_speed_chen_noelle(u_ll, u_rr, orientation::Integer, -# equations::ShallowWaterEquationsWetDry1D) +""" + min_max_speed_chen_noelle(u_ll, u_rr, orientation::Integer, + equations::ShallowWaterEquations1D) -# The approximated speeds for the HLL type numerical flux used by Chen and Noelle for their -# hydrostatic reconstruction. As they state in the paper, these speeds are chosen for the numerical -# flux to ensure positivity and to satisfy an entropy inequality. +The approximated speeds for the HLL type numerical flux used by Chen and Noelle for their +hydrostatic reconstruction. As they state in the paper, these speeds are chosen for the numerical +flux to ensure positivity and to satisfy an entropy inequality. -# Further details on this hydrostatic reconstruction and its motivation can be found in -# - Guoxian Chen and Sebastian Noelle (2017) -# A new hydrostatic reconstruction scheme based on subcell reconstructions -# [DOI:10.1137/15M1053074](https://dx.doi.org/10.1137/15M1053074) -# """ -# @inline function min_max_speed_chen_noelle(u_ll, u_rr, orientation::Integer, -# equations::ShallowWaterEquationsWetDry1D) -# # Get the velocity quantities -# v_ll = velocity(u_ll, equations) -# v_rr = velocity(u_rr, equations) +Further details on this hydrostatic reconstruction and its motivation can be found in +- Guoxian Chen and Sebastian Noelle (2017) + A new hydrostatic reconstruction scheme based on subcell reconstructions + [DOI:10.1137/15M1053074](https://dx.doi.org/10.1137/15M1053074) +""" +@inline function min_max_speed_chen_noelle(u_ll, u_rr, orientation::Integer, + equations::ShallowWaterEquationsWetDry1D) + # Get the velocity quantities + v_ll = velocity(u_ll, equations) + v_rr = velocity(u_rr, equations) -# # Calculate the wave celerity on the left and right -# h_ll = waterheight(u_ll, equations) -# h_rr = waterheight(u_rr, equations) + # Calculate the wave celerity on the left and right + h_ll = Trixi.waterheight(u_ll, equations) + h_rr = Trixi.waterheight(u_rr, equations) -# a_ll = sqrt(equations.gravity * h_ll) -# a_rr = sqrt(equations.gravity * h_rr) + a_ll = sqrt(equations.gravity * h_ll) + a_rr = sqrt(equations.gravity * h_rr) -# λ_min = min(v_ll - a_ll, v_rr - a_rr, zero(eltype(u_ll))) -# λ_max = max(v_ll + a_ll, v_rr + a_rr, zero(eltype(u_ll))) + λ_min = min(v_ll - a_ll, v_rr - a_rr, zero(eltype(u_ll))) + λ_max = max(v_ll + a_ll, v_rr + a_rr, zero(eltype(u_ll))) -# return λ_min, λ_max -# end + return λ_min, λ_max +end @inline function Trixi.max_abs_speeds(u, equations::ShallowWaterEquationsWetDry1D) return Trixi.max_abs_speeds(u, - Trixi.ShallowWaterEquations1D(gravity_constant = equations.gravity)) + equations.basic_swe) end -# # Helper function to extract the velocity vector from the conservative variables -# @inline function Trixi.velocity(u, equations::ShallowWaterEquationsWetDry1D) -# return Trixi.velocity(u, -# Trixi.ShallowWaterEquations1D(gravity_constant = equations.gravity)) -# end +# Calculate estimates for minimum and maximum wave speeds for HLL-type fluxes +@inline function Trixi.min_max_speed_naive(u_ll, u_rr, orientation::Integer, + equations::ShallowWaterEquationsWetDry1D) + return Trixi.min_max_speed_naive(u_ll, u_rr, orientation, + equations.basic_swe) +end + +# More refined estimates for minimum and maximum wave speeds for HLL-type fluxes +@inline function Trixi.min_max_speed_davis(u_ll, u_rr, orientation::Integer, + equations::ShallowWaterEquationsWetDry1D) + return Trixi.min_max_speed_davis(u_ll, u_rr, orientation, + equations.basic_swe) +end + +@inline function Trixi.min_max_speed_einfeldt(u_ll, u_rr, orientation::Integer, + equations::ShallowWaterEquationsWetDry1D) + return Trixi.min_max_speed_einfeldt(u_ll, u_rr, orientation, + equations.basic_swe) +end + +# Helper function to extract the velocity vector from the conservative variables +@inline function Trixi.velocity(u, equations::ShallowWaterEquationsWetDry1D) + return Trixi.velocity(u, + equations.basic_swe) +end # Convert conservative variables to primitive @inline function Trixi.cons2prim(u, equations::ShallowWaterEquationsWetDry1D) return Trixi.cons2prim(u, - Trixi.ShallowWaterEquations1D(gravity_constant = equations.gravity)) + equations.basic_swe) end # Convert conservative variables to entropy @@ -517,35 +567,35 @@ end # just carries the bottom topography values for convenience @inline function Trixi.cons2entropy(u, equations::ShallowWaterEquationsWetDry1D) return Trixi.cons2entropy(u, - Trixi.ShallowWaterEquations1D(gravity_constant = equations.gravity)) + equations.basic_swe) end # Convert entropy variables to conservative @inline function Trixi.entropy2cons(w, equations::ShallowWaterEquationsWetDry1D) return Trixi.entropy2cons(w, - Trixi.ShallowWaterEquations1D(gravity_constant = equations.gravity)) + equations.basic_swe) end # Convert primitive to conservative variables @inline function Trixi.prim2cons(prim, equations::ShallowWaterEquationsWetDry1D) return Trixi.prim2cons(prim, - Trixi.ShallowWaterEquations1D(gravity_constant = equations.gravity)) + equations.basic_swe) end -# @inline function Trixi.waterheight(u, equations::ShallowWaterEquationsWetDry1D) -# return Trixi.waterheight(u, -# Trixi.ShallowWaterEquations1D(gravity_constant = equations.gravity)) -# end +@inline function Trixi.waterheight(u, equations::ShallowWaterEquationsWetDry1D) + return Trixi.waterheight(u, + equations.basic_swe) +end -# @inline function Trixi.pressure(u, equations::ShallowWaterEquationsWetDry1D) -# return Trixi.pressure(u, -# Trixi.ShallowWaterEquations1D(gravity_constant = equations.gravity)) -# end +@inline function Trixi.pressure(u, equations::ShallowWaterEquationsWetDry1D) + return Trixi.pressure(u, + equations.basic_swe) +end -# @inline function Trixi.waterheight_pressure(u, equations::ShallowWaterEquationsWetDry1D) -# return Trixi.waterheight_pressure(u, -# Trixi.ShallowWaterEquations1D(gravity_constant = equations.gravity)) -# end +@inline function Trixi.waterheight_pressure(u, equations::ShallowWaterEquationsWetDry1D) + return Trixi.waterheight_pressure(u, + equations.basic_swe) +end # Entropy function for the shallow water equations is the total energy @inline function Trixi.entropy(cons, equations::ShallowWaterEquationsWetDry1D) @@ -555,13 +605,13 @@ end # Calculate total energy for a conservative state `cons` @inline function Trixi.energy_total(cons, equations::ShallowWaterEquationsWetDry1D) return Trixi.energy_total(cons, - Trixi.ShallowWaterEquations1D(gravity_constant = equations.gravity)) + equations.basic_swe) end # Calculate kinetic energy for a conservative state `cons` @inline function Trixi.energy_kinetic(u, equations::ShallowWaterEquationsWetDry1D) return Trixi.energy_kinetic(u, - Trixi.ShallowWaterEquations1D(gravity_constant = equations.gravity)) + equations.basic_swe) end # Calculate potential energy for a conservative state `cons` diff --git a/src/equations/shallow_water_wet_dry_2d.jl b/src/equations/shallow_water_wet_dry_2d.jl new file mode 100644 index 0000000..a9edf53 --- /dev/null +++ b/src/equations/shallow_water_wet_dry_2d.jl @@ -0,0 +1,846 @@ +# 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 + +@doc raw""" + ShallowWaterEquationsWetDry2D(; gravity, H0 = 0, threshold_limiter = nothing, threshold_wet = nothing) + +Shallow water equations (SWE) in two space dimensions. The equations are given by +```math +\begin{aligned} + \frac{\partial h}{\partial t} + \frac{\partial}{\partial x}(h v_1) + + \frac{\partial}{\partial y}(h v_2) &= 0 \\ + \frac{\partial}{\partial t}(h v_1) + \frac{\partial}{\partial x}\left(h v_1^2 + \frac{g}{2}h^2\right) + + \frac{\partial}{\partial y}(h v_1 v_2) + g h \frac{\partial b}{\partial x} &= 0 \\ + \frac{\partial}{\partial t}(h v_2) + \frac{\partial}{\partial x}(h v_1 v_2) + + \frac{\partial}{\partial y}\left(h v_2^2 + \frac{g}{2}h^2\right) + g h \frac{\partial b}{\partial y} &= 0. +\end{aligned} +``` +The unknown quantities of the SWE are the water height ``h`` and the velocities ``\mathbf{v} = (v_1, v_2)^T``. +The gravitational constant is denoted by `g` and the (possibly) variable bottom topography function ``b(x,y)``. +Conservative variable water height ``h`` is measured from the bottom topography ``b``, therefore one +also defines the total water height as ``H = h + b``. + +The additional quantity ``H_0`` is also available to store a reference value for the total water height that +is useful to set initial conditions or test the "lake-at-rest" well-balancedness. + +Also, there are two thresholds which prevent numerical problems as well as instabilities. Both of them do not +have to be passed, as default values are defined within the struct. The 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. + +The bottom topography function ``b(x,y)`` is set inside the initial condition routine +for a particular problem setup. To test the conservative form of the SWE one can set the bottom topography +variable `b` to zero. + +In addition to the unknowns, TrixiShallowWater.jl currently stores the bottom topography values at the approximation points +despite being fixed in time. This is done for convenience of computing the bottom topography gradients +on the fly during the approximation as well as computing auxiliary quantities like the total water height ``H`` +or the entropy variables. +This affects the implementation and use of these equations in various ways: +* The flux values corresponding to the bottom topography must be zero. +* The bottom topography values must be included when defining initial conditions, boundary conditions or + source terms. +* [`AnalysisCallback`](@ref) analyzes this variable. +* Trixi.jl's visualization tools will visualize the bottom topography by default. + +References for the SWE are many but a good introduction is available in Chapter 13 of the book: +- Randall J. LeVeque (2002) + Finite Volume Methods for Hyperbolic Problems + [DOI: 10.1017/CBO9780511791253](https://doi.org/10.1017/CBO9780511791253) +""" +struct ShallowWaterEquationsWetDry2D{RealT <: Real} <: + Trixi.AbstractShallowWaterEquations{2, 4} + gravity::RealT # gravitational constant + H0::RealT # constant "lake-at-rest" total water height + # `threshold_limiter` used in `PositivityPreservingLimiterShallowWater` on water height, + # as a (small) shift on the initial condition and cutoff before the next time step. + # Default is 500*eps() which in double precision is ≈1e-13. + threshold_limiter::RealT + # `threshold_wet` applied on water height to define when the flow is "wet" + # before calculating the numerical flux. + # Default is 5*eps() which in double precision is ≈1e-15. + threshold_wet::RealT + # Standard shallow water equations for dispatch on Trixi.jl functions + basic_swe::ShallowWaterEquations2D{RealT} +end + +# Allow for flexibility to set the gravitational constant within an elixir depending on the +# application where `gravity_constant=1.0` or `gravity_constant=9.81` are common values. +# The reference total water height H0 defaults to 0.0 but is used for the "lake-at-rest" +# well-balancedness test cases. +# Strict default values for thresholds that performed well in many numerical experiments +function ShallowWaterEquationsWetDry2D(; gravity_constant, H0 = zero(gravity_constant), + threshold_limiter = nothing, + threshold_wet = nothing) + T = promote_type(typeof(gravity_constant), typeof(H0)) + if threshold_limiter === nothing + threshold_limiter = 500 * eps(T) + end + if threshold_wet === nothing + threshold_wet = 5 * eps(T) + end + # Construct the standard SWE for dispatch. Even though the `basic_swe` already store the + # gravity constant and the total water height, we store an extra copy in + # `ShallowWaterEquationsWetDry2D` for convenience. + basic_swe = ShallowWaterEquations2D(gravity_constant = gravity_constant, H0 = H0) + + ShallowWaterEquationsWetDry2D(gravity_constant, H0, threshold_limiter, + threshold_wet, basic_swe) +end + +Trixi.have_nonconservative_terms(::ShallowWaterEquationsWetDry2D) = True() +Trixi.varnames(::typeof(cons2cons), ::ShallowWaterEquationsWetDry2D) = ("h", "h_v1", + "h_v2", "b") +# Note, we use the total water height, H = h + b, as the first primitive variable for easier +# visualization and setting initial conditions +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 +# make use of the existing functionality and introduce a number of wrapper functions, that dispatch +# to the ShallowWaterEquations2D. + +# Set initial conditions at physical location `x` for time `t` +""" + initial_condition_convergence_test(x, t, equations::ShallowWaterEquationsWetDry2D) + +A smooth initial condition used for convergence tests in combination with +[`source_terms_convergence_test`](@ref) +(and [`BoundaryConditionDirichlet(initial_condition_convergence_test)`](@ref) in non-periodic domains). +""" +function Trixi.initial_condition_convergence_test(x, t, + equations::ShallowWaterEquationsWetDry2D) + return Trixi.initial_condition_convergence_test(x, t, + equations.basic_swe) +end + +""" + source_terms_convergence_test(u, x, t, equations::ShallowWaterEquationsWetDry2D) + +Source terms used for convergence tests in combination with +[`initial_condition_convergence_test`](@ref) +(and [`BoundaryConditionDirichlet(initial_condition_convergence_test)`](@ref) in non-periodic domains). + +This manufactured solution source term is specifically designed for the bottom topography function +`b(x,y) = 2 + 0.5 * sin(sqrt(2)*pi*x) + 0.5 * sin(sqrt(2)*pi*y)` +as defined in [`initial_condition_convergence_test`](@ref). +""" +@inline function Trixi.source_terms_convergence_test(u, x, t, + equations::ShallowWaterEquationsWetDry2D) + return Trixi.source_terms_convergence_test(u, x, t, + equations.basic_swe) +end + +""" + initial_condition_weak_blast_wave(x, t, equations::ShallowWaterEquationsWetDry2D) + +A weak blast wave discontinuity useful for testing, e.g., total energy conservation. +Note for the shallow water equations to the total energy acts as a mathematical entropy function. +""" +function Trixi.initial_condition_weak_blast_wave(x, t, + equations::ShallowWaterEquationsWetDry2D) + return Trixi.initial_condition_weak_blast_wave(x, t, + equations.basic_swe) +end + +""" + boundary_condition_slip_wall(u_inner, normal_direction, x, t, surface_flux_function, + equations::ShallowWaterEquationsWetDry2D) +Create a boundary state by reflecting the normal velocity component and keep +the tangential velocity component unchanged. The boundary water height is taken from +the internal value. +For details see Section 9.2.5 of the book: +- Eleuterio F. Toro (2001) + Shock-Capturing Methods for Free-Surface Shallow Flows + 1st edition + ISBN 0471987662 +""" +@inline function Trixi.boundary_condition_slip_wall(u_inner, + normal_direction::AbstractVector, + x, t, + surface_flux_function, + equations::ShallowWaterEquationsWetDry2D) + # normalize the outward pointing direction + normal = normal_direction / norm(normal_direction) + + # compute the normal velocity + u_normal = normal[1] * u_inner[2] + normal[2] * u_inner[3] + + # create the "external" boundary solution state + u_boundary = SVector(u_inner[1], + u_inner[2] - 2.0 * u_normal * normal[1], + u_inner[3] - 2.0 * u_normal * normal[2], + u_inner[4]) + + # calculate the boundary flux + flux = surface_flux_function(u_inner, u_boundary, normal_direction, equations) + + return flux +end + +""" + boundary_condition_slip_wall(u_inner, orientation, direction, x, t, + surface_flux_function, equations::ShallowWaterEquationsWetDry2D) + +Should be used together with [`TreeMesh`](@ref). +""" +@inline function Trixi.boundary_condition_slip_wall(u_inner, orientation, + direction, x, t, + surface_flux_function, + equations::ShallowWaterEquationsWetDry2D) + ## get the appropriate normal vector from the orientation + if orientation == 1 + u_boundary = SVector(u_inner[1], -u_inner[2], u_inner[3], u_inner[4]) + else # orientation == 2 + u_boundary = SVector(u_inner[1], u_inner[2], -u_inner[3], u_inner[4]) + end + + # Calculate boundary flux + if iseven(direction) # u_inner is "left" of boundary, u_boundary is "right" of boundary + flux = surface_flux_function(u_inner, u_boundary, orientation, equations) + else # u_boundary is "left" of boundary, u_inner is "right" of boundary + flux = surface_flux_function(u_boundary, u_inner, orientation, equations) + end + + return flux +end + +# Calculate 1D flux for a single point +# Note, the bottom topography has no flux +@inline function Trixi.flux(u, orientation::Integer, + equations::ShallowWaterEquationsWetDry2D) + Trixi.flux(u, orientation, equations.basic_swe) +end + +# Calculate 1D flux for a single point in the normal direction +# Note, this directional vector is not normalized and the bottom topography has no flux +@inline function Trixi.flux(u, normal_direction::AbstractVector, + equations::ShallowWaterEquationsWetDry2D) + Trixi.flux(u, normal_direction, + equations.basic_swe) +end + +""" + flux_nonconservative_wintermeyer_etal(u_ll, u_rr, orientation::Integer, + equations::ShallowWaterEquationsWetDry2D) + flux_nonconservative_wintermeyer_etal(u_ll, u_rr, + normal_direction_ll ::AbstractVector, + normal_direction_average::AbstractVector, + equations::ShallowWaterEquationsWetDry2D) + +Non-symmetric two-point volume flux discretizing the nonconservative (source) term +that contains the gradient of the bottom topography [`ShallowWaterEquationsWetDry2D`](@ref). + +On curvilinear meshes, this nonconservative flux depends on both the +contravariant vector (normal direction) at the current node and the averaged +one. This is different from numerical fluxes used to discretize conservative +terms. + +Further details are available in the paper: +- Niklas Wintermeyer, Andrew R. Winters, Gregor J. Gassner and David A. Kopriva (2017) + An entropy stable nodal discontinuous Galerkin method for the two dimensional + shallow water equations on unstructured curvilinear meshes with discontinuous bathymetry + [DOI: 10.1016/j.jcp.2017.03.036](https://doi.org/10.1016/j.jcp.2017.03.036) +""" +@inline function Trixi.flux_nonconservative_wintermeyer_etal(u_ll, u_rr, + orientation::Integer, + equations::ShallowWaterEquationsWetDry2D) + Trixi.flux_nonconservative_wintermeyer_etal(u_ll, u_rr, orientation, + equations.basic_swe) +end + +@inline function Trixi.flux_nonconservative_wintermeyer_etal(u_ll, u_rr, + normal_direction_ll::AbstractVector, + normal_direction_average::AbstractVector, + equations::ShallowWaterEquationsWetDry2D) + return Trixi.flux_nonconservative_wintermeyer_etal(u_ll, u_rr, + normal_direction_ll, + normal_direction_average, + equations.basic_swe) +end + +""" + flux_nonconservative_fjordholm_etal(u_ll, u_rr, orientation::Integer, + equations::ShallowWaterEquationsWetDry2D) + flux_nonconservative_fjordholm_etal(u_ll, u_rr, + normal_direction_ll ::AbstractVector, + normal_direction_average::AbstractVector, + equations::ShallowWaterEquationsWetDry2D) + +Non-symmetric two-point surface flux discretizing the nonconservative (source) term of +that contains the gradient of the bottom topography [`ShallowWaterEquationsWetDry2D`](@ref). + +On curvilinear meshes, this nonconservative flux depends on both the +contravariant vector (normal direction) at the current node and the averaged +one. This is different from numerical fluxes used to discretize conservative +terms. + +This contains additional terms compared to [`flux_nonconservative_wintermeyer_etal`](@ref) +that account for possible discontinuities in the bottom topography function. +Thus, this flux should be used in general at interfaces. For flux differencing volume terms, +[`flux_nonconservative_wintermeyer_etal`](@ref) is analytically equivalent but slightly +cheaper. + +Further details for the original finite volume formulation are available in +- Ulrik S. Fjordholm, Siddhartha Mishr and Eitan Tadmor (2011) + Well-balanced and energy stable schemes for the shallow water equations with discontinuous topography + [DOI: 10.1016/j.jcp.2011.03.042](https://doi.org/10.1016/j.jcp.2011.03.042) +and for curvilinear 2D case in the paper: +- Niklas Wintermeyer, Andrew R. Winters, Gregor J. Gassner and David A. Kopriva (2017) + An entropy stable nodal discontinuous Galerkin method for the two dimensional + shallow water equations on unstructured curvilinear meshes with discontinuous bathymetry + [DOI: 10.1016/j.jcp.2017.03.036](https://doi.org/10.1016/j.jcp.2017.03.036) +""" +@inline function Trixi.flux_nonconservative_fjordholm_etal(u_ll, u_rr, + orientation::Integer, + equations::ShallowWaterEquationsWetDry2D) + Trixi.flux_nonconservative_fjordholm_etal(u_ll, u_rr, orientation, + equations.basic_swe) +end + +@inline function Trixi.flux_nonconservative_fjordholm_etal(u_ll, u_rr, + normal_direction_ll::AbstractVector, + normal_direction_average::AbstractVector, + equations::ShallowWaterEquationsWetDry2D) + Trixi.flux_nonconservative_fjordholm_etal(u_ll, u_rr, + normal_direction_ll, + normal_direction_average, + equations.basic_swe) +end + +""" + hydrostatic_reconstruction_audusse_etal(u_ll, u_rr, orientation_or_normal_direction, + equations::ShallowWaterEquationsWetDry2D) + +A particular type of hydrostatic reconstruction on the water height to guarantee well-balancedness +for a general bottom topography [`ShallowWaterEquationsWetDry2D`](@ref). The reconstructed solution states +`u_ll_star` and `u_rr_star` variables are used to evaluate the surface numerical flux at the interface. +Use in combination with the generic numerical flux routine [`FluxHydrostaticReconstruction`](@ref). + +Further details for the hydrostatic reconstruction and its motivation can be found in +- Emmanuel Audusse, François Bouchut, Marie-Odile Bristeau, Rupert Klein, and Benoit Perthame (2004) + A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows + [DOI: 10.1137/S1064827503431090](https://doi.org/10.1137/S1064827503431090) +""" +@inline function Trixi.hydrostatic_reconstruction_audusse_etal(u_ll, u_rr, + equations::ShallowWaterEquationsWetDry2D) + return Trixi.hydrostatic_reconstruction_audusse_etal(u_ll, u_rr, + equations.basic_swe) +end + +""" + hydrostatic_reconstruction_chen_noelle(u_ll, u_rr, orientation::Integer, + equations::ShallowWaterEquationsWetDry2D) + +A particular type of hydrostatic reconstruction of the water height to guarantee well-balancedness +for a general bottom topography of the [`ShallowWaterEquationsWetDry2D`](@ref). The reconstructed solution states +`u_ll_star` and `u_rr_star` variables are then used to evaluate the surface numerical flux at the interface. +The key idea is a linear reconstruction of the bottom and water height at the interfaces using subcells. +Use in combination with the generic numerical flux routine [`FluxHydrostaticReconstruction`](@ref). + +Further details on this hydrostatic reconstruction and its motivation can be found in +- Guoxian Chen and Sebastian Noelle (2017) + A new hydrostatic reconstruction scheme based on subcell reconstructions + [DOI:10.1137/15M1053074](https://dx.doi.org/10.1137/15M1053074) +""" +@inline function hydrostatic_reconstruction_chen_noelle(u_ll, u_rr, + equations::ShallowWaterEquationsWetDry2D) + # Unpack left and right water heights and bottom topographies + h_ll, _, _, b_ll = u_ll + h_rr, _, _, b_rr = u_rr + + # Get the velocities on either side + v1_ll, v2_ll = velocity(u_ll, equations) + v1_rr, v2_rr = velocity(u_rr, equations) + + H_ll = b_ll + h_ll + H_rr = b_rr + h_rr + + b_star = min(max(b_ll, b_rr), min(H_ll, H_rr)) + + # Compute the reconstructed water heights + h_ll_star = min(H_ll - b_star, h_ll) + h_rr_star = min(H_rr - b_star, h_rr) + + # Set the water height to be at least the value stored in the variable threshold after + # the hydrostatic reconstruction is applied and before the numerical flux is calculated + # to avoid numerical problem with arbitrary small values. Interfaces with a water height + # lower or equal to the threshold can be declared as dry. + # The default value for `threshold_wet` is ≈5*eps(), or 1e-15 in double precision, is set + # in the `ShallowWaterEquationsWetDry2D` struct. This threshold value can be changed in the constructor + # call of this equation struct in an elixir. + threshold = equations.threshold_wet + + if (h_ll_star <= threshold) + h_ll_star = threshold + v1_ll = zero(v1_ll) + v2_ll = zero(v2_ll) + end + + if (h_rr_star <= threshold) + h_rr_star = threshold + v1_rr = zero(v1_rr) + v2_rr = zero(v2_rr) + end + + # Create the conservative variables using the reconstruted water heights + u_ll_star = SVector(h_ll_star, h_ll_star * v1_ll, h_ll_star * v2_ll, b_ll) + u_rr_star = SVector(h_rr_star, h_rr_star * v1_rr, h_rr_star * v2_rr, b_rr) + + return u_ll_star, u_rr_star +end + +""" + flux_nonconservative_audusse_etal(u_ll, u_rr, orientation::Integer, + equations::ShallowWaterEquationsWetDry2D) + flux_nonconservative_audusse_etal(u_ll, u_rr, + normal_direction_ll ::AbstractVector, + normal_direction_average::AbstractVector, + equations::ShallowWaterEquationsWetDry2D) + +Non-symmetric two-point surface flux that discretizes the nonconservative (source) term. +The discretization uses the `hydrostatic_reconstruction_audusse_etal` on the conservative +variables. + +This hydrostatic reconstruction ensures that the finite volume numerical fluxes remain +well-balanced for discontinuous bottom topographies [`ShallowWaterEquationsWetDry2D`](@ref). +Should be used together with [`FluxHydrostaticReconstruction`](@ref) and +[`hydrostatic_reconstruction_audusse_etal`](@ref) in the surface flux to ensure consistency. + +Further details for the hydrostatic reconstruction and its motivation can be found in +- Emmanuel Audusse, François Bouchut, Marie-Odile Bristeau, Rupert Klein, and Benoit Perthame (2004) + A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows + [DOI: 10.1137/S1064827503431090](https://doi.org/10.1137/S1064827503431090) +""" +@inline function Trixi.flux_nonconservative_audusse_etal(u_ll, u_rr, + orientation::Integer, + equations::ShallowWaterEquationsWetDry2D) + return Trixi.flux_nonconservative_audusse_etal(u_ll, u_rr, orientation, + equations.basic_swe) +end + +@inline function Trixi.flux_nonconservative_audusse_etal(u_ll, u_rr, + normal_direction_ll::AbstractVector, + normal_direction_average::AbstractVector, + equations::ShallowWaterEquationsWetDry2D) + return Trixi.flux_nonconservative_audusse_etal(u_ll, u_rr, + normal_direction_ll, + normal_direction_average, + equations.basic_swe) +end + +""" + flux_nonconservative_chen_noelle(u_ll, u_rr, + orientation::Integer, + equations::ShallowWaterEquationsWetDry2D) + flux_nonconservative_chen_noelle(u_ll, u_rr, + normal_direction_ll ::AbstractVector, + normal_direction_average ::AbstractVector, + equations::ShallowWaterEquationsWetDry2D) + +Non-symmetric two-point surface flux that discretizes the nonconservative (source) term. +The discretization uses the [`hydrostatic_reconstruction_chen_noelle`](@ref) on the conservative +variables. + +Should be used together with [`FluxHydrostaticReconstruction`](@ref) and +[`hydrostatic_reconstruction_chen_noelle`](@ref) in the surface flux to ensure consistency. + +Further details on the hydrostatic reconstruction and its motivation can be found in +- Guoxian Chen and Sebastian Noelle (2017) + A new hydrostatic reconstruction scheme based on subcell reconstructions + [DOI:10.1137/15M1053074](https://dx.doi.org/10.1137/15M1053074) +""" +@inline function flux_nonconservative_chen_noelle(u_ll, u_rr, orientation::Integer, + equations::ShallowWaterEquationsWetDry2D) + # Pull the water height and bottom topography on the left + h_ll, _, _, b_ll = u_ll + h_rr, _, _, b_rr = u_rr + + H_ll = h_ll + b_ll + H_rr = h_rr + b_rr + + b_star = min(max(b_ll, b_rr), min(H_ll, H_rr)) + + # Create the hydrostatic reconstruction for the left solution state + u_ll_star, _ = hydrostatic_reconstruction_chen_noelle(u_ll, u_rr, equations) + + # Copy the reconstructed water height for easier to read code + h_ll_star = u_ll_star[1] + + z = zero(eltype(u_ll)) + # Includes two parts: + # (i) Diagonal (consistent) term from the volume flux that uses `b_ll` to avoid + # cross-averaging across a discontinuous bottom topography + # (ii) True surface part that uses `h_ll` and `h_ll_star` to handle discontinuous bathymetry + g = equations.gravity + if orientation == 1 + f = SVector(z, + g * h_ll * b_ll - g * (h_ll_star + h_ll) * (b_ll - b_star), + z, z) + else # orientation == 2 + f = SVector(z, z, + g * h_ll * b_ll - g * (h_ll_star + h_ll) * (b_ll - b_star), + z) + end + + return f +end + +@inline function flux_nonconservative_chen_noelle(u_ll, u_rr, + normal_direction_ll::AbstractVector, + normal_direction_average::AbstractVector, + equations::ShallowWaterEquationsWetDry2D) + # Pull the water height and bottom topography on the left + h_ll, _, _, b_ll = u_ll + h_rr, _, _, b_rr = u_rr + + H_ll = h_ll + b_ll + H_rr = h_rr + b_rr + + b_star = min(max(b_ll, b_rr), min(H_ll, H_rr)) + + # Create the hydrostatic reconstruction for the left solution state + u_ll_star, _ = hydrostatic_reconstruction_chen_noelle(u_ll, u_rr, equations) + + # Copy the reconstructed water height for easier to read code + h_ll_star = u_ll_star[1] + + # Comes in two parts: + # (i) Diagonal (consistent) term from the volume flux that uses `normal_direction_average` + # but we use `b_ll` to avoid cross-averaging across a discontinuous bottom topography + + f2 = normal_direction_average[1] * equations.gravity * h_ll * b_ll + f3 = normal_direction_average[2] * equations.gravity * h_ll * b_ll + + # (ii) True surface part that uses `normal_direction_ll`, `h_ll` and `h_ll_star` + # to handle discontinuous bathymetry + + f2 -= normal_direction_ll[1] * equations.gravity * (h_ll_star + h_ll) * + (b_ll - b_star) + f3 -= normal_direction_ll[2] * equations.gravity * (h_ll_star + h_ll) * + (b_ll - b_star) + + # First and last equations do not have a nonconservative flux + f1 = f4 = zero(eltype(u_ll)) + + return SVector(f1, f2, f3, f4) +end + +""" + flux_nonconservative_ersing_etal(u_ll, u_rr, orientation::Integer, + equations::ShallowWaterEquationsWetDry2D) + flux_nonconservative_ersing_etal(u_ll, u_rr, + normal_direction_ll::AbstractVector, + normal_direction_average::AbstractVector, + equations::ShallowWaterEquationsWetDry2D) + +!!! warning "Experimental code" + This numerical flux is experimental and may change in any future release. + +Non-symmetric path-conservative two-point volume flux discretizing the nonconservative (source) term +that contains the gradient of the bottom topography [`ShallowWaterEquationsWetDry2D`](@ref). + +On curvilinear meshes, this nonconservative flux depends on both the +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 +[`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 + curvilinear meshes + [DOI: 10.48550/arXiv.2306.12699](https://doi.org/10.48550/arXiv.2306.12699) +""" +@inline function Trixi.flux_nonconservative_ersing_etal(u_ll, u_rr, + orientation::Integer, + equations::ShallowWaterEquationsWetDry2D) + return Trixi.flux_nonconservative_ersing_etal(u_ll, u_rr, orientation, + equations.basic_swe) +end + +@inline function Trixi.flux_nonconservative_ersing_etal(u_ll, u_rr, + normal_direction_ll::AbstractVector, + normal_direction_average::AbstractVector, + equations::ShallowWaterEquationsWetDry2D) + Trixi.flux_nonconservative_ersing_etal(u_ll, u_rr, + normal_direction_ll, + normal_direction_average, + equations.basic_swe) +end + +""" + flux_fjordholm_etal(u_ll, u_rr, orientation_or_normal_direction, + equations::ShallowWaterEquationsWetDry2D) + +Total energy conservative (mathematical entropy for shallow water equations). When the bottom topography +is nonzero this should only be used as a surface flux otherwise the scheme will not be well-balanced. +For well-balancedness in the volume flux use [`flux_wintermeyer_etal`](@ref). + +Details are available in Eq. (4.1) in the paper: +- Ulrik S. Fjordholm, Siddhartha Mishr and Eitan Tadmor (2011) + Well-balanced and energy stable schemes for the shallow water equations with discontinuous topography + [DOI: 10.1016/j.jcp.2011.03.042](https://doi.org/10.1016/j.jcp.2011.03.042) +""" +@inline function Trixi.flux_fjordholm_etal(u_ll, u_rr, orientation::Integer, + equations::ShallowWaterEquationsWetDry2D) + Trixi.flux_fjordholm_etal(u_ll, u_rr, orientation, + equations.basic_swe) +end + +@inline function Trixi.flux_fjordholm_etal(u_ll, u_rr, normal_direction::AbstractVector, + equations::ShallowWaterEquationsWetDry2D) + return Trixi.flux_fjordholm_etal(u_ll, u_rr, normal_direction, + equations.basic_swe) +end + +""" + flux_wintermeyer_etal(u_ll, u_rr, orientation_or_normal_direction, + equations::ShallowWaterEquationsWetDry2D) + +Total energy conservative (mathematical entropy for shallow water equations) split form. +When the bottom topography is nonzero this scheme will be well-balanced when used as a `volume_flux`. +The `surface_flux` should still use, e.g., [`flux_fjordholm_etal`](@ref). + +Further details are available in Theorem 1 of the paper: +- Niklas Wintermeyer, Andrew R. Winters, Gregor J. Gassner and David A. Kopriva (2017) + An entropy stable nodal discontinuous Galerkin method for the two dimensional + shallow water equations on unstructured curvilinear meshes with discontinuous bathymetry + [DOI: 10.1016/j.jcp.2017.03.036](https://doi.org/10.1016/j.jcp.2017.03.036) +""" +@inline function Trixi.flux_wintermeyer_etal(u_ll, u_rr, orientation::Integer, + equations::ShallowWaterEquationsWetDry2D) + return Trixi.flux_wintermeyer_etal(u_ll, u_rr, orientation, + equations.basic_swe) +end + +@inline function Trixi.flux_wintermeyer_etal(u_ll, u_rr, + normal_direction::AbstractVector, + equations::ShallowWaterEquationsWetDry2D) + return Trixi.flux_wintermeyer_etal(u_ll, u_rr, normal_direction, + equations.basic_swe) +end + +# Specialized `DissipationLocalLaxFriedrichs` to avoid spurious dissipation in the bottom topography +@inline function (dissipation::DissipationLocalLaxFriedrichs)(u_ll, u_rr, + orientation_or_normal_direction, + equations::ShallowWaterEquationsWetDry2D) + return (dissipation::DissipationLocalLaxFriedrichs)(u_ll, u_rr, + orientation_or_normal_direction, + equations.basic_swe) +end + +# Specialized `FluxHLL` to avoid spurious dissipation in the bottom topography +@inline function (numflux::FluxHLL)(u_ll, u_rr, orientation_or_normal_direction, + equations::ShallowWaterEquationsWetDry2D) + λ_min, λ_max = numflux.min_max_speed(u_ll, u_rr, orientation_or_normal_direction, + equations) + + if λ_min >= 0 && λ_max >= 0 + return flux(u_ll, orientation_or_normal_direction, equations) + elseif λ_max <= 0 && λ_min <= 0 + return flux(u_rr, orientation_or_normal_direction, equations) + else + f_ll = flux(u_ll, orientation_or_normal_direction, equations) + f_rr = flux(u_rr, orientation_or_normal_direction, equations) + inv_λ_max_minus_λ_min = inv(λ_max - λ_min) + factor_ll = λ_max * inv_λ_max_minus_λ_min + factor_rr = λ_min * inv_λ_max_minus_λ_min + factor_diss = λ_min * λ_max * inv_λ_max_minus_λ_min + diss = u_rr - u_ll + return factor_ll * f_ll - factor_rr * f_rr + + factor_diss * SVector(diss[1], diss[2], diss[3], zero(eltype(u_ll))) + end +end + +""" + min_max_speed_chen_noelle(u_ll, u_rr, orientation::Integer, + equations::ShallowWaterEquations2D) + min_max_speed_chen_noelle(u_ll, u_rr, normal_direction::AbstractVector, + equations::ShallowWaterEquations2D) + +Special estimate of the minimal and maximal wave speed of the shallow water equations for +the left and right states `u_ll, u_rr`. These approximate speeds are used for the HLL-type +numerical flux [`flux_hll_chen_noelle`](@ref). These wave speed estimates +together with a particular hydrostatic reconstruction technique guarantee +that the numerical flux is positive and satisfies an entropy inequality. + +Further details on this hydrostatic reconstruction and its motivation can be found in +the reference below. The definition of the wave speeds are given in Equation (2.20). +- Guoxian Chen and Sebastian Noelle (2017) + A new hydrostatic reconstruction scheme based on subcell reconstructions + [DOI:10.1137/15M1053074](https://dx.doi.org/10.1137/15M1053074) +""" +@inline function min_max_speed_chen_noelle(u_ll, u_rr, orientation::Integer, + equations::ShallowWaterEquationsWetDry2D) + h_ll = Trixi.waterheight(u_ll, equations) + v1_ll, v2_ll = velocity(u_ll, equations) + h_rr = Trixi.waterheight(u_rr, equations) + v1_rr, v2_rr = velocity(u_rr, equations) + + a_ll = sqrt(equations.gravity * h_ll) + a_rr = sqrt(equations.gravity * h_rr) + + if orientation == 1 # x-direction + λ_min = min(v1_ll - a_ll, v1_rr - a_rr, zero(eltype(u_ll))) + λ_max = max(v1_ll + a_ll, v1_rr + a_rr, zero(eltype(u_ll))) + else # y-direction + λ_min = min(v2_ll - a_ll, v2_rr - a_rr, zero(eltype(u_ll))) + λ_max = max(v2_ll + a_ll, v2_rr + a_rr, zero(eltype(u_ll))) + end + + return λ_min, λ_max +end + +@inline function min_max_speed_chen_noelle(u_ll, u_rr, normal_direction::AbstractVector, + equations::ShallowWaterEquationsWetDry2D) + h_ll = Trixi.waterheight(u_ll, equations) + v1_ll, v2_ll = velocity(u_ll, equations) + h_rr = Trixi.waterheight(u_rr, equations) + v1_rr, v2_rr = velocity(u_rr, equations) + + v_normal_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2] + v_normal_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2] + + norm_ = norm(normal_direction) + + a_ll = sqrt(equations.gravity * h_ll) * norm_ + a_rr = sqrt(equations.gravity * h_rr) * norm_ + + λ_min = min(v_normal_ll - a_ll, v_normal_rr - a_rr, zero(eltype(u_ll))) + λ_max = max(v_normal_ll + a_ll, v_normal_rr + a_rr, zero(eltype(u_ll))) + + return λ_min, λ_max +end + +@inline function Trixi.max_abs_speeds(u, equations::ShallowWaterEquationsWetDry2D) + return Trixi.max_abs_speeds(u, + equations.basic_swe) +end + +# Calculate estimates for minimum and maximum wave speeds for HLL-type fluxes +@inline function Trixi.min_max_speed_naive(u_ll, u_rr, orientation::Integer, + equations::ShallowWaterEquationsWetDry2D) + return Trixi.min_max_speed_naive(u_ll, u_rr, orientation, + equations.basic_swe) +end + +@inline function Trixi.min_max_speed_naive(u_ll, u_rr, normal_direction::AbstractVector, + equations::ShallowWaterEquationsWetDry2D) + return Trixi.min_max_speed_naive(u_ll, u_rr, normal_direction, + equations.basic_swe) +end + +# More refined estimates for minimum and maximum wave speeds for HLL-type fluxes +@inline function Trixi.min_max_speed_davis(u_ll, u_rr, orientation::Integer, + equations::ShallowWaterEquationsWetDry2D) + return Trixi.min_max_speed_davis(u_ll, u_rr, orientation, + equations.basic_swe) +end + +@inline function Trixi.min_max_speed_davis(u_ll, u_rr, normal_direction::AbstractVector, + equations::ShallowWaterEquationsWetDry2D) + return Trixi.min_max_speed_davis(u_ll, u_rr, normal_direction, + equations.basic_swe) +end + +@inline function Trixi.min_max_speed_einfeldt(u_ll, u_rr, orientation::Integer, + equations::ShallowWaterEquationsWetDry2D) + return Trixi.min_max_speed_einfeldt(u_ll, u_rr, orientation, + equations.basic_swe) +end + +@inline function Trixi.min_max_speed_einfeldt(u_ll, u_rr, + normal_direction::AbstractVector, + equations::ShallowWaterEquationsWetDry2D) + return Trixi.min_max_speed_einfeldt(u_ll, u_rr, normal_direction, + equations.basic_swe) +end + +# Helper function to extract the velocity vector from the conservative variables +@inline function Trixi.velocity(u, equations::ShallowWaterEquationsWetDry2D) + return Trixi.velocity(u, equations.basic_swe) +end + +# Convert conservative variables to primitive +@inline function Trixi.cons2prim(u, equations::ShallowWaterEquationsWetDry2D) + return Trixi.cons2prim(u, equations.basic_swe) +end + +# Convert conservative variables to entropy +# Note, only the first three are the entropy variables, the fourth entry still +# just carries the bottom topography values for convenience +@inline function Trixi.cons2entropy(u, equations::ShallowWaterEquationsWetDry2D) + return Trixi.cons2entropy(u, + equations.basic_swe) +end + +# Convert entropy variables to conservative +@inline function Trixi.entropy2cons(w, equations::ShallowWaterEquationsWetDry2D) + return Trixi.entropy2cons(w, + equations.basic_swe) +end + +# Convert primitive to conservative variables +@inline function Trixi.prim2cons(prim, equations::ShallowWaterEquationsWetDry2D) + return Trixi.prim2cons(prim, + equations.basic_swe) +end + +@inline function Trixi.waterheight(u, equations::ShallowWaterEquationsWetDry2D) + return Trixi.waterheight(u, + equations.basic_swe) +end + +@inline function Trixi.pressure(u, equations::ShallowWaterEquationsWetDry2D) + return Trixi.pressure(u, equations.basic_swe) +end + +@inline function Trixi.waterheight_pressure(u, equations::ShallowWaterEquationsWetDry2D) + return Trixi.waterheight(u, equations) * Trixi.pressure(u, equations) +end + +# Entropy function for the shallow water equations is the total energy +@inline function Trixi.entropy(cons, equations::ShallowWaterEquationsWetDry2D) + Trixi.energy_total(cons, equations) +end + +# Calculate total energy for a conservative state `cons` +@inline function Trixi.energy_total(cons, equations::ShallowWaterEquationsWetDry2D) + Trixi.energy_total(cons, equations.basic_swe) +end + +# Calculate kinetic energy for a conservative state `cons` +@inline function Trixi.energy_kinetic(u, equations::ShallowWaterEquationsWetDry2D) + Trixi.energy_kinetic(u, equations.basic_swe) +end + +# Calculate potential energy for a conservative state `cons` +@inline function Trixi.energy_internal(cons, equations::ShallowWaterEquationsWetDry2D) + return energy_total(cons, equations) - energy_kinetic(cons, equations) +end + +# Calculate the error for the "lake-at-rest" test case where H = h+b should +# be a constant value over time. Note, assumes there is a single reference +# water height `H0` with which to compare. +@inline function Trixi.lake_at_rest_error(u, equations::ShallowWaterEquationsWetDry2D) + h, _, _, b = u + + # For well-balancedness testing with possible wet/dry regions the reference + # water height `H0` accounts for the possibility that the bottom topography + # can emerge out of the water as well as for the threshold offset to avoid + # division by a "hard" zero water heights as well. + H0_wet_dry = max(equations.H0, b + equations.threshold_limiter) + + return abs(H0_wet_dry - (h + b)) +end +end # @muladd diff --git a/src/solvers/indicators.jl b/src/solvers/indicators.jl new file mode 100644 index 0000000..0ac2a5a --- /dev/null +++ b/src/solvers/indicators.jl @@ -0,0 +1,86 @@ +# 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 + +""" + IndicatorHennemannGassnerShallowWater(equations::AbstractEquations, basis; + alpha_max=0.5, + alpha_min=0.001, + alpha_smooth=true, + variable) + +Modified version of the [`IndicatorHennemannGassner`](@ref) +indicator used for shock-capturing for shallow water equations. After +the element-wise values for the blending factors are computed an additional check +is made to see if the element is partially wet. In this case, partially wet elements +are set to use the pure finite volume scheme that is guaranteed to be well-balanced +for this wet/dry transition state of the flow regime. + +See also [`VolumeIntegralShockCapturingHG`](@ref). + +## References + +- Hennemann, Gassner (2020) + "A provably entropy stable subcell shock capturing approach for high order split form DG" + [arXiv: 2008.12044](https://arxiv.org/abs/2008.12044) +""" +struct IndicatorHennemannGassnerShallowWater{RealT <: Real, Variable, Cache} <: + Trixi.AbstractIndicator + alpha_max::RealT + alpha_min::RealT + alpha_smooth::Bool + variable::Variable + cache::Cache +end + +# this method is used when the indicator is constructed as for shock-capturing volume integrals +# of the shallow water equations +# It modifies the shock-capturing indicator to use full FV method in dry elements +# or partially dry elements containing a wet/dry transition. +function IndicatorHennemannGassnerShallowWater(equations::Trixi.AbstractShallowWaterEquations, + basis; + alpha_max = 0.5, + alpha_min = 0.001, + alpha_smooth = true, + variable) + alpha_max, alpha_min = promote(alpha_max, alpha_min) + cache = Trixi.create_cache(IndicatorHennemannGassner, equations, basis) + IndicatorHennemannGassnerShallowWater{typeof(alpha_max), typeof(variable), + typeof(cache)}(alpha_max, alpha_min, + alpha_smooth, variable, cache) +end + +function Base.show(io::IO, indicator::IndicatorHennemannGassnerShallowWater) + @nospecialize indicator # reduce precompilation time + + print(io, "IndicatorHennemannGassnerShallowWater(") + print(io, indicator.variable) + print(io, ", alpha_max=", indicator.alpha_max) + print(io, ", alpha_min=", indicator.alpha_min) + print(io, ", alpha_smooth=", indicator.alpha_smooth) + print(io, ")") +end + +function Base.show(io::IO, ::MIME"text/plain", + indicator::IndicatorHennemannGassnerShallowWater) + @nospecialize indicator # reduce precompilation time + + if get(io, :compact, false) + show(io, indicator) + else + setup = [ + "indicator variable" => indicator.variable, + "max. α" => indicator.alpha_max, + "min. α" => indicator.alpha_min, + "smooth α" => (indicator.alpha_smooth ? "yes" : "no"), + ] + Trixi.summary_box(io, "IndicatorHennemannGassnerShallowWater", setup) + end +end + +include("indicators_1d.jl") +include("indicators_2d.jl") +end # @muladd diff --git a/src/solvers/indicators_1d.jl b/src/solvers/indicators_1d.jl new file mode 100644 index 0000000..8a9f47d --- /dev/null +++ b/src/solvers/indicators_1d.jl @@ -0,0 +1,115 @@ +# 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 + +# Modified indicator for ShallowWaterEquationsWetDry1D to apply full FV method on elements +# containing some "dry" LGL nodes. That is, if an element is partially "wet" then it becomes a +# full FV element. +function (indicator_hg::IndicatorHennemannGassnerShallowWater)(u::AbstractArray{<:Any, + 3}, + mesh, + equations::ShallowWaterEquationsWetDry1D, + dg::DGSEM, cache; + kwargs...) + @unpack alpha_max, alpha_min, alpha_smooth, variable = indicator_hg + @unpack alpha, alpha_tmp, indicator_threaded, modal_threaded = indicator_hg.cache + # TODO: Taal refactor, when to `resize!` stuff changed possibly by AMR? + # Shall we implement `resize!(semi::AbstractSemidiscretization, new_size)` + # or just `resize!` whenever we call the relevant methods as we do now? + resize!(alpha, nelements(dg, cache)) + if alpha_smooth + resize!(alpha_tmp, nelements(dg, cache)) + end + + # magic parameters + threshold = 0.5 * 10^(-1.8 * (nnodes(dg))^0.25) + parameter_s = log((1 - 0.0001) / 0.0001) + + # If the water height `h` at one LGL node is lower than `threshold_partially_wet` + # the indicator sets the element-wise blending factor alpha[element] = 1 + # via the local variable `indicator_wet`. In turn, this ensures that a pure + # FV method is used in partially wet elements and guarantees the well-balanced property. + # + # Hard-coded cut-off value of `threshold_partially_wet = 1e-4` was determined through many numerical experiments. + # Overall idea is to increase robustness when computing the velocity on (nearly) dry elements which + # could be "dangerous" due to division of conservative variables, e.g., v = hv / h. + # Here, the impact of the threshold on the number of elements being updated with FV is not that + # significant. However, its impact on the robustness is very significant. + # The value can be seen as a trade-off between accuracy and stability. + # Well-balancedness of the scheme on partially wet elements with hydrostatic reconstruction + # can only be proven for the FV method (see Chen and Noelle). + # Therefore we set alpha to one regardless of its given maximum value. + threshold_partially_wet = 1e-4 + + Trixi.@threaded for element in eachelement(dg, cache) + indicator = indicator_threaded[Threads.threadid()] + modal = modal_threaded[Threads.threadid()] + + # (Re-)set dummy variable for alpha_dry + indicator_wet = 1 + + # Calculate indicator variables at Gauss-Lobatto nodes + for i in eachnode(dg) + u_local = get_node_vars(u, equations, dg, i, element) + h, _, _ = u_local + + if h <= threshold_partially_wet + indicator_wet = 0 + end + + indicator[i] = indicator_hg.variable(u_local, equations) + end + + # Convert to modal representation + Trixi.multiply_scalar_dimensionwise!(modal, + dg.basis.inverse_vandermonde_legendre, + indicator) + + # Calculate total energies for all modes, without highest, without two highest + total_energy = zero(eltype(modal)) + for i in 1:nnodes(dg) + total_energy += modal[i]^2 + end + total_energy_clip1 = zero(eltype(modal)) + for i in 1:(nnodes(dg) - 1) + total_energy_clip1 += modal[i]^2 + end + total_energy_clip2 = zero(eltype(modal)) + for i in 1:(nnodes(dg) - 2) + total_energy_clip2 += modal[i]^2 + end + + # Calculate energy in higher modes + energy = max((total_energy - total_energy_clip1) / total_energy, + (total_energy_clip1 - total_energy_clip2) / total_energy_clip1) + + alpha_element = 1 / (1 + exp(-parameter_s / threshold * (energy - threshold))) + + # Take care of the case close to pure DG + if alpha_element < alpha_min + alpha_element = zero(alpha_element) + end + + # Take care of the case close to pure FV + if alpha_element > 1 - alpha_min + alpha_element = one(alpha_element) + end + + # Clip the maximum amount of FV allowed or set to one depending on indicator_wet + if indicator_wet == 0 + alpha[element] = 1 + else # Element is not defined as dry but wet + alpha[element] = min(alpha_max, alpha_element) + end + end + + if alpha_smooth + Trixi.apply_smoothing!(mesh, alpha, alpha_tmp, dg, cache) + end + + return alpha +end +end # @muladd diff --git a/src/solvers/indicators_2d.jl b/src/solvers/indicators_2d.jl new file mode 100644 index 0000000..1e36414 --- /dev/null +++ b/src/solvers/indicators_2d.jl @@ -0,0 +1,116 @@ +# 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 + +# Modified indicator for ShallowWaterEquationsWetDry2D to apply full FV method on elements +# containing some "dry" LGL nodes. That is, if an element is partially "wet" then it becomes a +# full FV element. +function (indicator_hg::IndicatorHennemannGassnerShallowWater)(u::AbstractArray{<:Any, + 4}, + mesh, + equations::ShallowWaterEquationsWetDry2D, + dg::DGSEM, cache; + kwargs...) + @unpack alpha_max, alpha_min, alpha_smooth, variable = indicator_hg + @unpack alpha, alpha_tmp, indicator_threaded, modal_threaded, modal_tmp1_threaded = indicator_hg.cache + # TODO: Taal refactor, when to `resize!` stuff changed possibly by AMR? + # Shall we implement `resize!(semi::AbstractSemidiscretization, new_size)` + # or just `resize!` whenever we call the relevant methods as we do now? + resize!(alpha, nelements(dg, cache)) + if alpha_smooth + resize!(alpha_tmp, nelements(dg, cache)) + end + + # magic parameters + threshold = 0.5 * 10^(-1.8 * (nnodes(dg))^0.25) + parameter_s = log((1 - 0.0001) / 0.0001) + + # If the water height `h` at one LGL node is lower than `threshold_partially_wet` + # the indicator sets the element-wise blending factor alpha[element] = 1 + # via the local variable `indicator_wet`. In turn, this ensures that a pure + # FV method is used in partially wet elements and guarantees the well-balanced property. + # + # Hard-coded cut-off value of `threshold_partially_wet = 1e-4` was determined through many numerical experiments. + # Overall idea is to increase robustness when computing the velocity on (nearly) dry elements which + # could be "dangerous" due to division of conservative variables, e.g., v1 = hv1 / h. + # Here, the impact of the threshold on the number of elements being updated with FV is not that + # significant. However, its impact on the robustness is very significant. + # The value can be seen as a trade-off between accuracy and stability. + # Well-balancedness of the scheme on partially wet elements with hydrostatic reconstruction + # can only be proven for the FV method (see Chen and Noelle). + # Therefore we set alpha to be one regardless of its given value from the modal indicator. + threshold_partially_wet = 1e-4 + + Trixi.@threaded for element in eachelement(dg, cache) + indicator = indicator_threaded[Threads.threadid()] + modal = modal_threaded[Threads.threadid()] + modal_tmp1 = modal_tmp1_threaded[Threads.threadid()] + + # (Re-)set dummy variable for alpha_dry + indicator_wet = 1 + + # Calculate indicator variables at Gauss-Lobatto nodes + for j in eachnode(dg), i in eachnode(dg) + u_local = get_node_vars(u, equations, dg, i, j, element) + h, _, _, _ = u_local + + if h <= threshold_partially_wet + indicator_wet = 0 + end + + indicator[i, j] = indicator_hg.variable(u_local, equations) + end + + # Convert to modal representation + Trixi.multiply_scalar_dimensionwise!(modal, + dg.basis.inverse_vandermonde_legendre, + indicator, modal_tmp1) + + # Calculate total energies for all modes, without highest, without two highest + total_energy = zero(eltype(modal)) + for j in 1:nnodes(dg), i in 1:nnodes(dg) + total_energy += modal[i, j]^2 + end + total_energy_clip1 = zero(eltype(modal)) + for j in 1:(nnodes(dg) - 1), i in 1:(nnodes(dg) - 1) + total_energy_clip1 += modal[i, j]^2 + end + total_energy_clip2 = zero(eltype(modal)) + for j in 1:(nnodes(dg) - 2), i in 1:(nnodes(dg) - 2) + total_energy_clip2 += modal[i, j]^2 + end + + # Calculate energy in higher modes + energy = max((total_energy - total_energy_clip1) / total_energy, + (total_energy_clip1 - total_energy_clip2) / total_energy_clip1) + + alpha_element = 1 / (1 + exp(-parameter_s / threshold * (energy - threshold))) + + # Take care of the case close to pure DG + if alpha_element < alpha_min + alpha_element = zero(alpha_element) + end + + # Take care of the case close to pure FV + if alpha_element > 1 - alpha_min + alpha_element = one(alpha_element) + end + + # Clip the maximum amount of FV allowed or set to 1 depending on indicator_wet + if indicator_wet == 0 + alpha[element] = 1 + else # Element is not defined as dry but wet + alpha[element] = min(alpha_max, alpha_element) + end + end + + if alpha_smooth + Trixi.apply_smoothing!(mesh, alpha, alpha_tmp, dg, cache) + end + + return alpha +end +end # @muladd diff --git a/test/Project.toml b/test/Project.toml index cfac638..d3420ab 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -2,11 +2,15 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" Trixi = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" +Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6" +Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" [compat] Test = "1" -Trixi = "0.5, 0.6" +Trixi = "0.7" OrdinaryDiffEq = "6.49.1" +Downloads = "1" +Printf = "1" [preferences.OrdinaryDiffEq] PrecompileAutoSpecialize = false diff --git a/test/runtests.jl b/test/runtests.jl index 92a68d6..d501312 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -2,17 +2,16 @@ using Trixi using TrixiShallowWater using Test -# We run tests in parallel with CI jobs setting the `TRIXI_TEST` environment +# We run tests with CI jobs setting the `TRIXI_TEST` environment # variable to determine the subset of tests to execute. -# By default, we just run the threaded tests since they are relatively cheap -# and test a good amount of different functionality. const TRIXI_TEST = get(ENV, "TRIXI_TEST", "all") -const TRIXI_MPI_NPROCS = clamp(Sys.CPU_THREADS, 2, 3) -const TRIXI_NTHREADS = clamp(Sys.CPU_THREADS, 2, 3) @time @testset "TrixiShallowWater.jl tests" begin @time if TRIXI_TEST == "all" - include("test_tree_1d_shallowwater_wet_dry.jl") + include("test_tree_1d.jl") + include("test_tree_2d.jl") + include("test_unstructured_2d.jl") + include("test_structured_2d.jl") include("test_unit.jl") end diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl new file mode 100644 index 0000000..6312471 --- /dev/null +++ b/test/test_structured_2d.jl @@ -0,0 +1,148 @@ +module TestExamplesStructuredMesh2D + +using Test +using Trixi +using TrixiShallowWater + +include("test_trixi.jl") + +EXAMPLES_DIR = pkgdir(TrixiShallowWater, "examples", "structured_2d_dgsem") + +# Start with a clean environment: remove TrixiShallowWater.jl output directory if it exists +outdir = "out" +isdir(outdir) && rm(outdir, recursive = true) + +@testset "StructuredMesh2D" begin +#! format: noindent + +@trixi_testset "elixir_shallowwater_source_terms.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"), + l2=[ + 0.0017285599436729316, + 0.025584610912606776, + 0.028373834961180594, + 6.274146767730866e-5, + ], + linf=[ + 0.012972309788264802, + 0.108283714215621, + 0.15831585777928936, + 0.00018196759554722775, + ], + tspan=(0.0, 0.05)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end + +@trixi_testset "elixir_shallowwater_well_balanced.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_well_balanced.jl"), + l2=[ + 0.7920927046419308, + 9.92129670988898e-15, + 1.0118635033124588e-14, + 0.7920927046419308, + ], + linf=[ + 2.408429868800133, + 5.5835419986809516e-14, + 5.448874313931364e-14, + 2.4084298688001335, + ], + tspan=(0.0, 0.25)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end + +@trixi_testset "elixir_shallowwater_well_balanced_wet_dry.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_well_balanced_wet_dry.jl"), + l2=[ + 0.019731646454942086, + 1.0694532773278277e-14, + 1.1969913383405568e-14, + 0.0771517260037954, + ], + linf=[ + 0.4999999999998892, + 6.067153702623552e-14, + 4.4849667259339357e-14, + 1.9999999999999993, + ], + tspan=(0.0, 0.25)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end + +@trixi_testset "elixir_shallowwater_conical_island.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_conical_island.jl"), + l2=[ + 0.04593154164306353, + 0.1644534881916908, + 0.16445348819169076, + 0.0011537702354532122, + ], + linf=[ + 0.21100717610846442, + 0.9501592344310412, + 0.950159234431041, + 0.021790250683516296, + ], + tspan=(0.0, 0.025)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end + +@trixi_testset "elixir_shallowwater_parabolic_bowl.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_parabolic_bowl.jl"), + l2=[ + 0.00015285369980313484, + 1.9536806395943226e-5, + 9.936906607758672e-5, + 5.0686313334616055e-15, + ], + linf=[ + 0.003316119030459211, + 0.0005075409427972817, + 0.001986721761060583, + 4.701794509287538e-14, + ], + tspan=(0.0, 0.025), cells_per_dimension=(40, 40)) + # 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 + +# Clean up afterwards: delete TrixiShallowWater.jl output directory +@test_nowarn rm(outdir, recursive = true) + +end # module diff --git a/test/test_tree_1d.jl b/test/test_tree_1d.jl index 7d96552..e5a1dee 100644 --- a/test/test_tree_1d.jl +++ b/test/test_tree_1d.jl @@ -1,27 +1,401 @@ module TestExamplesTree1D using Test +using Trixi using TrixiShallowWater include("test_trixi.jl") EXAMPLES_DIR = pkgdir(TrixiShallowWater, "examples", "tree_1d_dgsem") -# Start with a clean environment: remove Trixi.jl output directory if it exists +# Start with a clean environment: remove TrixiShallowWater.jl output directory if it exists outdir = "out" isdir(outdir) && rm(outdir, recursive = true) @testset "TreeMesh1D" begin #! format: noindent -# Run basic tests -@testset "Examples 1D" begin - # Shallow water - include("test_tree_1d_shallowwater_wet_dry.jl") +@trixi_testset "elixir_shallowwater_ec.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_ec.jl"), + l2=[0.244729018751225, 0.8583565222389505, 0.07330427577586297], + linf=[ + 2.1635021283528504, + 3.8717508164234453, + 1.7711213427919539, + ], + tspan=(0.0, 0.25)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end end -# Clean up afterwards: delete Trixi.jl output directory -@test_nowarn rm(outdir, recursive = true) +@trixi_testset "elixir_shallowwater_ec.jl with initial_condition_weak_blast_wave" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_ec.jl"), + l2=[ + 0.39464782107209717, + 2.03880864210846, + 4.1623084150546725e-10, + ], + linf=[ + 0.778905801278281, + 3.2409883402608273, + 7.419800190922032e-10, + ], + initial_condition=initial_condition_weak_blast_wave, + tspan=(0.0, 0.25)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end + +@trixi_testset "elixir_shallowwater_well_balanced.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_well_balanced.jl"), + l2=[ + 0.10416666834254829, + 1.4352935256803184e-14, + 0.10416666834254838, + ], + linf=[1.9999999999999996, 3.248036646353028e-14, 2.0], + tspan=(0.0, 0.25)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end + +@trixi_testset "elixir_shallowwater_well_balanced.jl with FluxHydrostaticReconstruction" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_well_balanced.jl"), + l2=[ + 0.10416666834254835, + 1.1891029971551825e-14, + 0.10416666834254838, + ], + linf=[2.0000000000000018, 2.4019608337954543e-14, 2.0], + surface_flux=(FluxHydrostaticReconstruction(flux_lax_friedrichs, + hydrostatic_reconstruction_audusse_etal), + flux_nonconservative_audusse_etal), + tspan=(0.0, 0.25)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end + +@trixi_testset "elixir_shallowwater_well_balanced.jl with flux_nonconservative_ersing_etal" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_well_balanced.jl"), + l2=[ + 0.10416666834254838, + 1.6657566141935285e-14, + 0.10416666834254838, + ], + linf=[2.0000000000000004, 3.0610625110157164e-14, 2.0], + surface_flux=(flux_wintermeyer_etal, + flux_nonconservative_ersing_etal), + volume_flux=(flux_wintermeyer_etal, + flux_nonconservative_ersing_etal), + tspan=(0.0, 0.25)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end + +@trixi_testset "elixir_shallowwater_well_balanced_wet_dry.jl with FluxHydrostaticReconstruction" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_well_balanced_wet_dry.jl"), + l2=[ + 0.00965787167169024, + 5.345454081916856e-14, + 0.03857583749209928, + ], + linf=[ + 0.4999999999998892, + 2.2447689894899726e-13, + 1.9999999999999714, + ], + tspan=(0.0, 0.25), + # Soften the tolerance as test results vary between different CPUs + atol=1000 * eps()) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end + +@trixi_testset "elixir_shallowwater_source_terms.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"), + l2=[ + 0.0022363707373868713, + 0.01576799981934617, + 4.436491725585346e-5, + ], + linf=[ + 0.00893601803417754, + 0.05939797350246456, + 9.098379777405796e-5, + ], + tspan=(0.0, 0.025)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end + +@trixi_testset "elixir_shallowwater_source_terms.jl with flux_hll" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"), + l2=[ + 0.0022758146627220154, + 0.015864082886204556, + 4.436491725585346e-5, + ], + linf=[ + 0.008457195427364006, + 0.057201667446161064, + 9.098379777405796e-5, + ], + tspan=(0.0, 0.025), + surface_flux=(FluxHLL(min_max_speed_naive), + flux_nonconservative_fjordholm_etal)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end + +@trixi_testset "elixir_shallowwater_source_terms.jl with flux_nonconservative_ersing_etal" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"), + l2=[ + 0.005774284062933275, + 0.017408601639513584, + 4.43649172561843e-5, + ], + linf=[ + 0.01639116193303547, + 0.05102877460799604, + 9.098379777450205e-5, + ], + surface_flux=(flux_wintermeyer_etal, + flux_nonconservative_ersing_etal), + volume_flux=(flux_wintermeyer_etal, + flux_nonconservative_ersing_etal), + tspan=(0.0, 0.025)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end + +@trixi_testset "elixir_shallowwater_source_terms_dirichlet.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_source_terms_dirichlet.jl"), + l2=[ + 0.0022851099219788917, + 0.01560453773635554, + 4.43649172558535e-5, + ], + linf=[ + 0.008934615705174398, + 0.059403169140869405, + 9.098379777405796e-5, + ], + tspan=(0.0, 0.025)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end + +@trixi_testset "elixir_shallowwater_source_terms_dirichlet.jl with FluxHydrostaticReconstruction" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_source_terms_dirichlet.jl"), + l2=[ + 0.0022956052733432287, + 0.015540053559855601, + 4.43649172558535e-5, + ], + linf=[ + 0.008460440313118323, + 0.05720939349382359, + 9.098379777405796e-5, + ], + surface_flux=(FluxHydrostaticReconstruction(FluxHLL(min_max_speed_naive), + hydrostatic_reconstruction_audusse_etal), + flux_nonconservative_audusse_etal), + tspan=(0.0, 0.025)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end + +@trixi_testset "elixir_shallowwater_well_balanced_nonperiodic.jl with Dirichlet boundary" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_well_balanced_nonperiodic.jl"), + l2=[ + 1.725964362045055e-8, + 5.0427180314307505e-16, + 1.7259643530442137e-8, + ], + linf=[ + 3.844551077492042e-8, + 3.469453422316143e-15, + 3.844551077492042e-8, + ], + tspan=(0.0, 0.25), + surface_flux=(FluxHLL(min_max_speed_naive), + flux_nonconservative_fjordholm_etal),) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end + +@trixi_testset "elixir_shallowwater_well_balanced_nonperiodic.jl with wall boundary" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_well_balanced_nonperiodic.jl"), + l2=[ + 1.7259643614361866e-8, + 3.5519018243195145e-16, + 1.7259643530442137e-8, + ], + linf=[ + 3.844551010878661e-8, + 9.846474508971374e-16, + 3.844551077492042e-8, + ], + tspan=(0.0, 0.25), + surface_flux=(FluxHLL(min_max_speed_naive), + flux_nonconservative_fjordholm_etal), + boundary_condition=boundary_condition_slip_wall) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end + +@trixi_testset "elixir_shallowwater_shock_capturing.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_shock_capturing.jl"), + l2=[0.07424140641160326, 0.2148642632748155, 0.0372579849000542], + linf=[ + 1.1209754279344226, + 1.3230788645853582, + 0.8646939843534251, + ], + tspan=(0.0, 0.05)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end + +@trixi_testset "elixir_shallowwater_beach.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_beach.jl"), + l2=[ + 0.17979210479598923, + 1.2377495706611434, + 6.289818963361573e-8, + ], + linf=[ + 0.845938394800688, + 3.3740800777086575, + 4.4541473087633676e-7, + ], + tspan=(0.0, 0.05), + atol=3e-10) # see https://github.com/trixi-framework/Trixi.jl/issues/1617 + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end + +@trixi_testset "elixir_shallowwater_parabolic_bowl.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_parabolic_bowl.jl"), + l2=[ + 8.965981683033589e-5, + 1.8565707397810857e-5, + 4.1043039226164336e-17, + ], + linf=[ + 0.00041080213807871235, + 0.00014823261488938177, + 2.220446049250313e-16, + ], + tspan=(0.0, 0.05)) + # 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 # TreeMesh1D +# Clean up afterwards: delete TrixiShallowWater.jl output directory +@test_nowarn rm(outdir, recursive = true) + end # module diff --git a/test/test_tree_1d_shallowwater_wet_dry.jl b/test/test_tree_1d_shallowwater_wet_dry.jl deleted file mode 100644 index e2cec1a..0000000 --- a/test/test_tree_1d_shallowwater_wet_dry.jl +++ /dev/null @@ -1,392 +0,0 @@ -module TestExamples1DShallowWaterWetDry - -using Test -using Trixi -using TrixiShallowWater - -include("test_trixi.jl") - -EXAMPLES_DIR = pkgdir(TrixiShallowWater, "examples", "tree_1d_dgsem") - -@testset "Shallow Water Wet Dry" begin -#! format: noindent - -@trixi_testset "elixir_shallowwater_ec.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_ec.jl"), - l2=[0.244729018751225, 0.8583565222389505, 0.07330427577586297], - linf=[ - 2.1635021283528504, - 3.8717508164234453, - 1.7711213427919539, - ], - tspan=(0.0, 0.25)) - # Ensure that we do not have excessive memory allocations - # (e.g., from type instabilities) - let - t = sol.t[end] - u_ode = sol.u[end] - du_ode = similar(u_ode) - @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 - end -end - -@trixi_testset "elixir_shallowwater_ec.jl with initial_condition_weak_blast_wave" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_ec.jl"), - l2=[ - 0.39464782107209717, - 2.03880864210846, - 4.1623084150546725e-10, - ], - linf=[ - 0.778905801278281, - 3.2409883402608273, - 7.419800190922032e-10, - ], - initial_condition=initial_condition_weak_blast_wave, - tspan=(0.0, 0.25)) - # Ensure that we do not have excessive memory allocations - # (e.g., from type instabilities) - let - t = sol.t[end] - u_ode = sol.u[end] - du_ode = similar(u_ode) - @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 - end -end - -@trixi_testset "elixir_shallowwater_well_balanced.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_well_balanced.jl"), - l2=[ - 0.10416666834254829, - 1.4352935256803184e-14, - 0.10416666834254838, - ], - linf=[1.9999999999999996, 3.248036646353028e-14, 2.0], - tspan=(0.0, 0.25)) - # Ensure that we do not have excessive memory allocations - # (e.g., from type instabilities) - let - t = sol.t[end] - u_ode = sol.u[end] - du_ode = similar(u_ode) - @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 - end -end - -@trixi_testset "elixir_shallowwater_well_balanced.jl with FluxHydrostaticReconstruction" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_well_balanced.jl"), - l2=[ - 0.10416666834254835, - 1.1891029971551825e-14, - 0.10416666834254838, - ], - linf=[2.0000000000000018, 2.4019608337954543e-14, 2.0], - surface_flux=(FluxHydrostaticReconstruction(flux_lax_friedrichs, - hydrostatic_reconstruction_audusse_etal), - flux_nonconservative_audusse_etal), - tspan=(0.0, 0.25)) - # Ensure that we do not have excessive memory allocations - # (e.g., from type instabilities) - let - t = sol.t[end] - u_ode = sol.u[end] - du_ode = similar(u_ode) - @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 - end -end - -@trixi_testset "elixir_shallowwater_well_balanced.jl with flux_nonconservative_ersing_etal" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_well_balanced.jl"), - l2=[ - 0.10416666834254838, - 1.6657566141935285e-14, - 0.10416666834254838, - ], - linf=[2.0000000000000004, 3.0610625110157164e-14, 2.0], - surface_flux=(flux_wintermeyer_etal, - flux_nonconservative_ersing_etal), - volume_flux=(flux_wintermeyer_etal, - flux_nonconservative_ersing_etal), - tspan=(0.0, 0.25)) - # Ensure that we do not have excessive memory allocations - # (e.g., from type instabilities) - let - t = sol.t[end] - u_ode = sol.u[end] - du_ode = similar(u_ode) - @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 - end -end - -# TODO: Activate test when Wet_Dry functionality is moved -# @trixi_testset "elixir_shallowwater_well_balanced_wet_dry.jl with FluxHydrostaticReconstruction" begin -# @test_trixi_include(joinpath(EXAMPLES_DIR, -# "elixir_shallowwater_well_balanced_wet_dry.jl"), -# l2=[ -# 0.00965787167169024, -# 5.345454081916856e-14, -# 0.03857583749209928, -# ], -# linf=[ -# 0.4999999999998892, -# 2.2447689894899726e-13, -# 1.9999999999999714, -# ], -# tspan=(0.0, 0.25), -# # Soften the tolerance as test results vary between different CPUs -# atol=1000 * eps()) -# # Ensure that we do not have excessive memory allocations -# # (e.g., from type instabilities) -# let -# t = sol.t[end] -# u_ode = sol.u[end] -# du_ode = similar(u_ode) -# @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 -# end -# end - -@trixi_testset "elixir_shallowwater_source_terms.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"), - l2=[ - 0.0022363707373868713, - 0.01576799981934617, - 4.436491725585346e-5, - ], - linf=[ - 0.00893601803417754, - 0.05939797350246456, - 9.098379777405796e-5, - ], - tspan=(0.0, 0.025)) - # Ensure that we do not have excessive memory allocations - # (e.g., from type instabilities) - let - t = sol.t[end] - u_ode = sol.u[end] - du_ode = similar(u_ode) - @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 - end -end - -@trixi_testset "elixir_shallowwater_source_terms.jl with flux_hll" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"), - l2=[ - 0.0022758146627220154, - 0.015864082886204556, - 4.436491725585346e-5, - ], - linf=[ - 0.008457195427364006, - 0.057201667446161064, - 9.098379777405796e-5, - ], - tspan=(0.0, 0.025), - surface_flux=(flux_hll, flux_nonconservative_fjordholm_etal)) - # Ensure that we do not have excessive memory allocations - # (e.g., from type instabilities) - let - t = sol.t[end] - u_ode = sol.u[end] - du_ode = similar(u_ode) - @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 - end -end - -@trixi_testset "elixir_shallowwater_source_terms.jl with flux_nonconservative_ersing_etal" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"), - l2=[ - 0.005774284062933275, - 0.017408601639513584, - 4.43649172561843e-5, - ], - linf=[ - 0.01639116193303547, - 0.05102877460799604, - 9.098379777450205e-5, - ], - surface_flux=(flux_wintermeyer_etal, - flux_nonconservative_ersing_etal), - volume_flux=(flux_wintermeyer_etal, - flux_nonconservative_ersing_etal), - tspan=(0.0, 0.025)) - # Ensure that we do not have excessive memory allocations - # (e.g., from type instabilities) - let - t = sol.t[end] - u_ode = sol.u[end] - du_ode = similar(u_ode) - @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 - end -end - -@trixi_testset "elixir_shallowwater_source_terms_dirichlet.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, - "elixir_shallowwater_source_terms_dirichlet.jl"), - l2=[ - 0.0022851099219788917, - 0.01560453773635554, - 4.43649172558535e-5, - ], - linf=[ - 0.008934615705174398, - 0.059403169140869405, - 9.098379777405796e-5, - ], - tspan=(0.0, 0.025)) - # Ensure that we do not have excessive memory allocations - # (e.g., from type instabilities) - let - t = sol.t[end] - u_ode = sol.u[end] - du_ode = similar(u_ode) - @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 - end -end - -@trixi_testset "elixir_shallowwater_source_terms_dirichlet.jl with FluxHydrostaticReconstruction" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, - "elixir_shallowwater_source_terms_dirichlet.jl"), - l2=[ - 0.0022956052733432287, - 0.015540053559855601, - 4.43649172558535e-5, - ], - linf=[ - 0.008460440313118323, - 0.05720939349382359, - 9.098379777405796e-5, - ], - surface_flux=(FluxHydrostaticReconstruction(flux_hll, - hydrostatic_reconstruction_audusse_etal), - flux_nonconservative_audusse_etal), - tspan=(0.0, 0.025)) - # Ensure that we do not have excessive memory allocations - # (e.g., from type instabilities) - let - t = sol.t[end] - u_ode = sol.u[end] - du_ode = similar(u_ode) - @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 - end -end - -@trixi_testset "elixir_shallowwater_well_balanced_nonperiodic.jl with Dirichlet boundary" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, - "elixir_shallowwater_well_balanced_nonperiodic.jl"), - l2=[ - 1.725964362045055e-8, - 5.0427180314307505e-16, - 1.7259643530442137e-8, - ], - linf=[ - 3.844551077492042e-8, - 3.469453422316143e-15, - 3.844551077492042e-8, - ], - tspan=(0.0, 0.25)) - # Ensure that we do not have excessive memory allocations - # (e.g., from type instabilities) - let - t = sol.t[end] - u_ode = sol.u[end] - du_ode = similar(u_ode) - @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 - end -end - -@trixi_testset "elixir_shallowwater_well_balanced_nonperiodic.jl with wall boundary" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, - "elixir_shallowwater_well_balanced_nonperiodic.jl"), - l2=[ - 1.7259643614361866e-8, - 3.5519018243195145e-16, - 1.7259643530442137e-8, - ], - linf=[ - 3.844551010878661e-8, - 9.846474508971374e-16, - 3.844551077492042e-8, - ], - tspan=(0.0, 0.25), - boundary_condition=boundary_condition_slip_wall) - # 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 - -# TODO: Activate test when Wet_Dry functionality is moved -# @trixi_testset "elixir_shallowwater_shock_capturing.jl" begin -# @test_trixi_include(joinpath(EXAMPLES_DIR, -# "elixir_shallowwater_shock_capturing.jl"), -# l2=[0.07424140641160326, 0.2148642632748155, 0.0372579849000542], -# linf=[ -# 1.1209754279344226, -# 1.3230788645853582, -# 0.8646939843534251, -# ], -# tspan=(0.0, 0.05)) -# # 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 - -# TODO: Activate test when Wet_Dry functionality is moved -# @trixi_testset "elixir_shallowwater_beach.jl" begin -# @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_beach.jl"), -# l2=[ -# 0.17979210479598923, -# 1.2377495706611434, -# 6.289818963361573e-8, -# ], -# linf=[ -# 0.845938394800688, -# 3.3740800777086575, -# 4.4541473087633676e-7, -# ], -# tspan=(0.0, 0.05), -# atol=3e-10) # see https://github.com/trixi-framework/Trixi.jl/issues/1617 -# # 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 - -# TODO: Activate test when Wet_Dry functionality is moved -# @trixi_testset "elixir_shallowwater_parabolic_bowl.jl" begin -# @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_parabolic_bowl.jl"), -# l2=[ -# 8.965981683033589e-5, -# 1.8565707397810857e-5, -# 4.1043039226164336e-17, -# ], -# linf=[ -# 0.00041080213807871235, -# 0.00014823261488938177, -# 2.220446049250313e-16, -# ], -# tspan=(0.0, 0.05)) -# # 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 -end # module diff --git a/test/test_tree_2d.jl b/test/test_tree_2d.jl new file mode 100644 index 0000000..e0c324c --- /dev/null +++ b/test/test_tree_2d.jl @@ -0,0 +1,364 @@ +module TestExamples2DShallowWaterWetDry + +using Test +using Trixi +using TrixiShallowWater + +include("test_trixi.jl") + +EXAMPLES_DIR = pkgdir(TrixiShallowWater, "examples", "tree_2d_dgsem") + +# Start with a clean environment: remove TrixiShallowWater.jl output directory if it exists +outdir = "out" +isdir(outdir) && rm(outdir, recursive = true) + +@testset "TreeMesh2D" begin +#! format: noindent + +@trixi_testset "elixir_shallowwater_ec.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_ec.jl"), + l2=[ + 0.991181203601035, + 0.734130029040644, + 0.7447696147162621, + 0.5875351036989047, + ], + linf=[ + 2.0117744577945413, + 2.9962317608172127, + 2.6554999727293653, + 3.0, + ], + tspan=(0.0, 0.25)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end + +@trixi_testset "elixir_shallowwater_well_balanced.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_well_balanced.jl"), + l2=[ + 0.9130579602987144, + 1.0602847041965408e-14, + 1.082225645390032e-14, + 0.9130579602987147, + ], + linf=[ + 2.113062037615659, + 4.6613606802974e-14, + 5.4225772771633196e-14, + 2.1130620376156584, + ], + tspan=(0.0, 0.25)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end + +@trixi_testset "elixir_shallowwater_well_balanced_wall.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_well_balanced_wall.jl"), + l2=[ + 0.9130579602987144, + 1.0602847041965408e-14, + 1.082225645390032e-14, + 0.9130579602987147, + ], + linf=[ + 2.113062037615659, + 4.6613606802974e-14, + 5.4225772771633196e-14, + 2.1130620376156584, + ], + tspan=(0.0, 0.25)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end + +@trixi_testset "elixir_shallowwater_well_balanced.jl with FluxHydrostaticReconstruction" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_well_balanced.jl"), + l2=[ + 0.9130579602987147, + 9.68729463970494e-15, + 9.694538537436981e-15, + 0.9130579602987147, + ], + linf=[ + 2.1130620376156584, + 2.3875905654916432e-14, + 2.2492839032269154e-14, + 2.1130620376156584, + ], + surface_flux=(FluxHydrostaticReconstruction(flux_lax_friedrichs, + hydrostatic_reconstruction_audusse_etal), + flux_nonconservative_audusse_etal), + tspan=(0.0, 0.25)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end + +@trixi_testset "elixir_shallowwater_well_balanced.jl with flux_nonconservative_ersing_etal" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_well_balanced.jl"), + l2=[ + 0.9130579602987146, + 1.0323158914614244e-14, + 1.0276096319430528e-14, + 0.9130579602987147, + ], + linf=[ + 2.11306203761566, + 4.063916419044386e-14, + 3.694484044448245e-14, + 2.1130620376156584, + ], + surface_flux=(flux_wintermeyer_etal, + flux_nonconservative_ersing_etal), + volume_flux=(flux_wintermeyer_etal, + flux_nonconservative_ersing_etal), + tspan=(0.0, 0.25)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end + +@trixi_testset "elixir_shallowwater_well_balanced_wet_dry.jl with FluxHydrostaticReconstruction" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_well_balanced_wet_dry.jl"), + l2=[ + 0.030186039395610056, + 2.513287752536758e-14, + 1.3631397744897607e-16, + 0.10911781485920438, + ], + linf=[ + 0.49999999999993505, + 5.5278950497971455e-14, + 7.462550826772548e-16, + 2.0, + ], + tspan=(0.0, 0.25)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end + +@trixi_testset "elixir_shallowwater_source_terms.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"), + l2=[ + 0.001868474306068482, + 0.01731687445878443, + 0.017649083171490863, + 6.274146767717023e-5, + ], + linf=[ + 0.016962486402209986, + 0.08768628853889782, + 0.09038488750767648, + 0.0001819675955490041, + ], + tspan=(0.0, 0.025)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end + +@trixi_testset "elixir_shallowwater_source_terms_dirichlet.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_source_terms_dirichlet.jl"), + l2=[ + 0.0018746929418489125, + 0.017332321628469628, + 0.01634953679145536, + 6.274146767717023e-5, + ], + linf=[ + 0.016262353691956388, + 0.08726160620859424, + 0.09043621801418844, + 0.0001819675955490041, + ], + tspan=(0.0, 0.025)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end + +@trixi_testset "elixir_shallowwater_source_terms.jl with flux_hll" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"), + l2=[ + 0.0018957692481057034, + 0.016943229710439864, + 0.01755623297390675, + 6.274146767717414e-5, + ], + linf=[ + 0.015156105797771602, + 0.07964811135780492, + 0.0839787097210376, + 0.0001819675955490041, + ], + tspan=(0.0, 0.025), + surface_flux=(FluxHLL(min_max_speed_naive), + flux_nonconservative_fjordholm_etal)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end + +@trixi_testset "elixir_shallowwater_source_terms.jl with flux_nonconservative_ersing_etal" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"), + l2=[ + 0.002471853426064005, + 0.05619168608950033, + 0.11844727575152562, + 6.274146767730281e-5, + ], + linf=[ + 0.014332922987500218, + 0.2141204806174546, + 0.5392313755637872, + 0.0001819675955490041, + ], + surface_flux=(flux_wintermeyer_etal, + flux_nonconservative_ersing_etal), + volume_flux=(flux_wintermeyer_etal, + flux_nonconservative_ersing_etal), + tspan=(0.0, 0.25)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end + +@trixi_testset "elixir_shallowwater_conical_island.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_conical_island.jl"), + l2=[ + 0.0459315416430658, + 0.1644534881916991, + 0.16445348819169914, + 0.0011537702354532694, + ], + linf=[ + 0.21100717610846464, + 0.9501592344310412, + 0.9501592344310417, + 0.021790250683516282, + ], + tspan=(0.0, 0.025)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end + +@trixi_testset "elixir_shallowwater_parabolic_bowl.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_parabolic_bowl.jl"), + l2=[ + 0.00025345501281482687, + 4.4525120338817177e-5, + 0.00015991819160294247, + 7.750412064917294e-15, + ], + linf=[ + 0.004664246019836723, + 0.0004972780116736669, + 0.0028735707270457628, + 6.866729407306593e-14, + ], + tspan=(0.0, 0.025), + basis=LobattoLegendreBasis(3)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end + +@trixi_testset "elixir_shallowwater_wall.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_wall.jl"), + l2=[ + 0.13517233723296504, + 0.20010876311162215, + 0.20010876311162223, + 2.719538414346464e-7, + ], + linf=[ + 0.5303607982988336, + 0.5080989745682338, + 0.5080989745682352, + 1.1301675764130437e-6, + ], + 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 # TreeMesh2D + +# Clean up afterwards: delete TrixiShallowWater.jl output directory +@test_nowarn rm(outdir, recursive = true) + +end # module diff --git a/test/test_unit.jl b/test/test_unit.jl index 9517443..9d77370 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -6,7 +6,7 @@ using TrixiShallowWater include("test_trixi.jl") -# Start with a clean environment: remove Trixi.jl output directory if it exists +# Start with a clean environment: remove TrixiShallowWater.jl output directory if it exists outdir = "out" isdir(outdir) && rm(outdir, recursive = true) @@ -14,25 +14,138 @@ isdir(outdir) && rm(outdir, recursive = true) @testset "Unit tests" begin #! format: noindent +@timed_testset "Printing indicators/controllers" begin + # OBS! Constructing indicators/controllers using the parameters below doesn't make sense. It's + # just useful to run basic tests of `show` methods. + + indicator_hg_swe = IndicatorHennemannGassnerShallowWater(1.0, 0.0, true, "variable", + "cache") + @test_nowarn show(stdout, indicator_hg_swe) +end + @timed_testset "Shallow water conversion between conservative/entropy variables" begin - H, v, b = 3.5, 0.25, 0.4 + H, v1, v2, b = 3.5, 0.25, 0.1, 0.4 let equations = ShallowWaterEquationsWetDry1D(gravity_constant = 9.8) - cons_vars = Trixi.prim2cons(SVector(H, v, b), equations) - entropy_vars = Trixi.cons2entropy(cons_vars, equations) - @test cons_vars ≈ Trixi.entropy2cons(entropy_vars, equations) + cons_vars = prim2cons(SVector(H, v1, b), equations) + entropy_vars = cons2entropy(cons_vars, equations) + @test cons_vars ≈ entropy2cons(entropy_vars, equations) - total_energy = Trixi.energy_total(cons_vars, equations) - @test total_energy ≈ Trixi.entropy(cons_vars, equations) + total_energy = energy_total(cons_vars, equations) + @test total_energy ≈ entropy(cons_vars, equations) @test total_energy ≈ - Trixi.energy_internal(cons_vars, equations) + + energy_internal(cons_vars, equations) + energy_kinetic(cons_vars, equations) # test tuple args - cons_vars = Trixi.prim2cons((H, v, b), equations) - entropy_vars = Trixi.cons2entropy(cons_vars, equations) - @test cons_vars ≈ Trixi.entropy2cons(entropy_vars, equations) + cons_vars = prim2cons((H, v1, b), equations) + entropy_vars = cons2entropy(cons_vars, equations) + @test cons_vars ≈ entropy2cons(entropy_vars, equations) + end + + let equations = ShallowWaterEquationsWetDry2D(gravity_constant = 9.8) + cons_vars = prim2cons(SVector(H, v1, v2, b), equations) + entropy_vars = cons2entropy(cons_vars, equations) + @test cons_vars ≈ entropy2cons(entropy_vars, equations) + + total_energy = energy_total(cons_vars, equations) + @test total_energy ≈ entropy(cons_vars, equations) + + # test tuple args + cons_vars = prim2cons((H, v1, v2, b), equations) + entropy_vars = cons2entropy(cons_vars, equations) + @test cons_vars ≈ entropy2cons(entropy_vars, equations) end end + +@timed_testset "Consistency check for waterheight_pressure" begin + H, v1, v2, b = 3.5, 0.25, 0.1, 0.4 + + let equations = ShallowWaterEquationsWetDry1D(gravity_constant = 9.8) + cons_vars = prim2cons(SVector(H, v1, b), equations) + @test waterheight_pressure(cons_vars, equations) ≈ + Trixi.waterheight(cons_vars, equations) * pressure(cons_vars, equations) + end + + let equations = ShallowWaterEquationsWetDry2D(gravity_constant = 9.8) + cons_vars = prim2cons(SVector(H, v1, v2, b), equations) + @test waterheight_pressure(cons_vars, equations) ≈ + Trixi.waterheight(cons_vars, equations) * pressure(cons_vars, equations) + end +end + +# Test consistency of the wrapper functions for the wave speed estimates +@timed_testset "Consistency check for wave speed estimates of the SWE" begin + let equations = ShallowWaterEquationsWetDry1D(gravity_constant = 9.8) + u_rr = SVector(1.2, 0.3, 0.7) + u_ll = SVector(0.25, 0.1, 0.4) + orientation = 1 + + @test min_max_speed_naive(u_ll, u_rr, orientation, + equations) == + min_max_speed_naive(u_ll, u_rr, orientation, + equations.basic_swe) + @test min_max_speed_davis(u_ll, u_rr, orientation, + equations) == + min_max_speed_davis(u_ll, u_rr, orientation, + equations.basic_swe) + @test min_max_speed_einfeldt(u_ll, u_rr, orientation, + equations) == + min_max_speed_einfeldt(u_ll, u_rr, orientation, + equations.basic_swe) + @test min_max_speed_naive(u_ll, u_rr, orientation, + equations) == + min_max_speed_naive(u_ll, u_rr, orientation, + equations.basic_swe) + end + + let equations = ShallowWaterEquationsWetDry2D(gravity_constant = 9.8) + u_rr = SVector(1.2, 0.3, 0.2, 0.7) + u_ll = SVector(0.25, 0.1, 0.3, 0.4) + orientations = [1, 2] + normal_directions = [SVector(1.0, 0.0), + SVector(0.0, 1.0), + SVector(0.5, -0.5), + SVector(-1.2, 0.3)] + + for orientation in orientations + @test min_max_speed_naive(u_ll, u_rr, orientation, + equations) == + min_max_speed_naive(u_ll, u_rr, orientation, + equations.basic_swe) + @test min_max_speed_davis(u_ll, u_rr, orientation, + equations) == + min_max_speed_davis(u_ll, u_rr, orientation, + equations.basic_swe) + @test min_max_speed_einfeldt(u_ll, u_rr, orientation, + equations) == + min_max_speed_einfeldt(u_ll, u_rr, orientation, + equations.basic_swe) + @test min_max_speed_naive(u_ll, u_rr, orientation, + equations) == + min_max_speed_naive(u_ll, u_rr, orientation, + equations.basic_swe) + end + + for normal_direction in normal_directions + @test min_max_speed_naive(u_ll, u_rr, normal_direction, + equations) == + min_max_speed_naive(u_ll, u_rr, normal_direction, + equations.basic_swe) + @test min_max_speed_davis(u_ll, u_rr, normal_direction, + equations) == + min_max_speed_davis(u_ll, u_rr, normal_direction, + equations.basic_swe) + @test min_max_speed_einfeldt(u_ll, u_rr, normal_direction, + equations) == + min_max_speed_einfeldt(u_ll, u_rr, normal_direction, + equations.basic_swe) + @test min_max_speed_naive(u_ll, u_rr, normal_direction, + equations) == + min_max_speed_naive(u_ll, u_rr, normal_direction, + equations.basic_swe) + end + end end +end # Unit tests end # module diff --git a/test/test_unstructured_2d.jl b/test/test_unstructured_2d.jl new file mode 100644 index 0000000..c2ca180 --- /dev/null +++ b/test/test_unstructured_2d.jl @@ -0,0 +1,420 @@ +module TestExamplesUnstructuredMesh2D + +using Test +using Trixi +using TrixiShallowWater + +include("test_trixi.jl") + +EXAMPLES_DIR = pkgdir(TrixiShallowWater, "examples", "unstructured_2d_dgsem") + +# Start with a clean environment: remove TrixiShallowWater.jl output directory if it exists +outdir = "out" +isdir(outdir) && rm(outdir, recursive = true) + +@testset "UnstructuredMesh2D" begin +#! format: noindent + +@trixi_testset "elixir_shallowwater_ec.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_ec.jl"), + l2=[ + 0.6106939484178353, + 0.48586236867426724, + 0.48234490854514356, + 0.29467422718511727, + ], + linf=[ + 2.775979948281604, + 3.1721242154451548, + 3.5713448319601393, + 2.052861364219655, + ], + tspan=(0.0, 0.25)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end + +@trixi_testset "elixir_shallowwater_well_balanced.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_well_balanced.jl"), + l2=[ + 1.2164292510839076, + 2.6118925543469468e-12, + 1.1636046671473883e-12, + 1.2164292510839079, + ], + linf=[ + 1.5138512282315846, + 4.998482888288039e-11, + 2.0246214978154587e-11, + 1.513851228231574, + ], + tspan=(0.0, 0.25)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end + +@trixi_testset "elixir_shallowwater_well_balanced.jl with FluxHydrostaticReconstruction" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_well_balanced.jl"), + l2=[ + 1.2164292510839085, + 1.2643106818778908e-12, + 7.46884905098358e-13, + 1.2164292510839079, + ], + linf=[ + 1.513851228231562, + 1.6287765844373185e-11, + 6.8766999132716964e-12, + 1.513851228231574, + ], + surface_flux=(FluxHydrostaticReconstruction(flux_lax_friedrichs, + hydrostatic_reconstruction_audusse_etal), + flux_nonconservative_audusse_etal), + tspan=(0.0, 0.2)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end + +@trixi_testset "elixir_shallowwater_well_balanced.jl with flux_nonconservative_ersing_etal" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_well_balanced.jl"), + l2=[ + 1.2164292510839083, + 2.590643638636187e-12, + 1.0945471514840143e-12, + 1.2164292510839079, + ], + linf=[ + 1.5138512282315792, + 5.0276441977281156e-11, + 1.9816934589292803e-11, + 1.513851228231574, + ], + surface_flux=(flux_wintermeyer_etal, + flux_nonconservative_ersing_etal), + volume_flux=(flux_wintermeyer_etal, + flux_nonconservative_ersing_etal), + tspan=(0.0, 0.25)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end + +@trixi_testset "elixir_shallowwater_source_terms.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"), + l2=[ + 0.0011197623982310795, + 0.04456344888447023, + 0.014317376629669337, + 5.089218476758975e-6, + ], + linf=[ + 0.007835284004819698, + 0.3486891284278597, + 0.11242778979399048, + 2.6407324614119432e-5, + ], + tspan=(0.0, 0.025)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end + +@trixi_testset "elixir_shallowwater_source_terms.jl with FluxHydrostaticReconstruction" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"), + l2=[ + 0.0011197139793938152, + 0.015430259691310781, + 0.017081031802719724, + 5.089218476758271e-6, + ], + linf=[ + 0.014300809338967824, + 0.12783372461225184, + 0.17625472321992852, + 2.6407324614341476e-5, + ], + surface_flux=(FluxHydrostaticReconstruction(FluxHLL(min_max_speed_naive), + hydrostatic_reconstruction_audusse_etal), + flux_nonconservative_audusse_etal), + tspan=(0.0, 0.025)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end + +@trixi_testset "elixir_shallowwater_source_terms.jl with flux_nonconservative_ersing_etal" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"), + l2=[ + 0.0011196687776346434, + 0.044562672453443995, + 0.014306265289763618, + 5.089218476759981e-6, + ], + linf=[ + 0.007825021762002393, + 0.348550815397918, + 0.1115517935018282, + 2.6407324614341476e-5, + ], + surface_flux=(flux_wintermeyer_etal, + flux_nonconservative_ersing_etal), + volume_flux=(flux_wintermeyer_etal, + flux_nonconservative_ersing_etal), + tspan=(0.0, 0.025)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end + +@trixi_testset "elixir_shallowwater_source_terms.jl with flux_hll" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"), + l2=[ + 0.0011197139793938727, + 0.015430259691311309, + 0.017081031802719554, + 5.089218476759981e-6, + ], + linf=[ + 0.014300809338967824, + 0.12783372461224918, + 0.17625472321993918, + 2.6407324614341476e-5, + ], + surface_flux=(FluxHLL(min_max_speed_naive), + flux_nonconservative_fjordholm_etal), + tspan=(0.0, 0.025)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end + +@trixi_testset "elixir_shallowwater_dirichlet.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_dirichlet.jl"), + l2=[ + 1.1577518608940115e-5, + 4.867189932537344e-13, + 4.647273240470541e-13, + 1.1577518608933468e-5, + ], + linf=[ + 8.394063878602864e-5, + 1.1469760027632646e-10, + 1.1146619484429974e-10, + 8.394063879602065e-5, + ], + tspan=(0.0, 2.0), + surface_flux=(FluxHLL(min_max_speed_naive), + flux_nonconservative_fjordholm_etal)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end + +@trixi_testset "elixir_shallowwater_wall_bc_shockcapturing.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_wall_bc_shockcapturing.jl"), + l2=[ + 0.04444388691670699, + 0.1527771788033111, + 0.1593763537203512, + 6.225080476986749e-8, + ], + linf=[ + 0.6526506870169639, + 1.980765893182952, + 2.4807635459119757, + 3.982097158683473e-7, + ], + tspan=(0.0, 0.05), + surface_flux=(FluxHydrostaticReconstruction(FluxHLL(min_max_speed_naive), + hydrostatic_reconstruction_audusse_etal), + flux_nonconservative_audusse_etal)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end + +@trixi_testset "elixir_shallowwater_ec_shockcapturing.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_ec_shockcapturing.jl"), + l2=[ + 0.6124656312639043, + 0.504371951785709, + 0.49180896200746366, + 0.29467422718511727, + ], + linf=[ + 2.7639232436274392, + 3.3985508653311767, + 3.3330308209196224, + 2.052861364219655, + ], + tspan=(0.0, 0.25)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end + +@trixi_testset "elixir_shallowwater_three_mound_dam_break.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_three_mound_dam_break.jl"), + l2=[ + 0.0892957892027502, + 0.30648836484407915, + 2.28712547616214e-15, + 0.0008778654298684622, + ], + linf=[ + 0.850329472915091, + 2.330631694956507, + 5.783660020252348e-14, + 0.04326237921249021, + ], + basis=LobattoLegendreBasis(3), + 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 + +# TODO: activate when twolayer equations are available +# @trixi_testset "elixir_shallowwater_twolayer_convergence.jl" begin +# @test_trixi_include(joinpath(EXAMPLES_DIR, +# "elixir_shallowwater_twolayer_convergence.jl"), +# l2=[0.0007935561625451243, 0.008825315509943844, +# 0.002429969315645897, +# 0.0007580145888686304, 0.004495741879625235, +# 0.0015758146898767814, +# 6.849532064729749e-6], +# linf=[0.0059205195991136605, 0.08072126590166251, +# 0.03463806075399023, +# 0.005884818649227186, 0.042658506561995546, +# 0.014125956138838602, 2.5829318284764646e-5], +# tspan=(0.0, 0.25)) +# # Ensure that we do not have excessive memory allocations +# # (e.g., from type instabilities) +# let +# t = sol.t[end] +# u_ode = sol.u[end] +# du_ode = similar(u_ode) +# @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 +# end +# end + +# @trixi_testset "elixir_shallowwater_twolayer_well_balanced.jl" begin +# @test_trixi_include(joinpath(EXAMPLES_DIR, +# "elixir_shallowwater_twolayer_well_balanced.jl"), +# l2=[4.706532184998499e-16, 1.1215950712872183e-15, +# 6.7822712922421565e-16, +# 0.002192812926266047, 5.506855295923691e-15, +# 3.3105180099689275e-15, +# 0.0021928129262660085], +# linf=[4.468647674116255e-15, 1.3607872120431166e-14, +# 9.557155049520056e-15, +# 0.024280130945632084, 6.68910907640583e-14, +# 4.7000983997100496e-14, +# 0.024280130945632732], +# tspan=(0.0, 0.25)) +# # Ensure that we do not have excessive memory allocations +# # (e.g., from type instabilities) +# let +# t = sol.t[end] +# u_ode = sol.u[end] +# du_ode = similar(u_ode) +# @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 +# end +# end + +# @trixi_testset "elixir_shallowwater_twolayer_dam_break.jl with flux_lax_friedrichs" begin +# @test_trixi_include(joinpath(EXAMPLES_DIR, +# "elixir_shallowwater_twolayer_dam_break.jl"), +# l2=[0.012447632879122346, 0.012361250464676683, +# 0.0009551519536340908, +# 0.09119400061322577, 0.015276216721920347, +# 0.0012126995108983853, 0.09991983966647647], +# linf=[0.044305765721807444, 0.03279620980615845, +# 0.010754320388190101, +# 0.111309922939555, 0.03663360204931427, +# 0.014332822306649284, +# 0.10000000000000003], +# surface_flux=(flux_lax_friedrichs, +# flux_nonconservative_ersing_etal), +# tspan=(0.0, 0.25)) +# # Ensure that we do not have excessive memory allocations +# # (e.g., from type instabilities) +# let +# t = sol.t[end] +# u_ode = sol.u[end] +# du_ode = similar(u_ode) +# @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 +# end +# end +end + +# Clean up afterwards: delete TrixiShallowWater.jl output directory +@test_nowarn rm(outdir, recursive = true) + +end # module diff --git a/test/test_upstream.jl b/test/test_upstream.jl index 3beb128..34b354b 100644 --- a/test/test_upstream.jl +++ b/test/test_upstream.jl @@ -18,7 +18,7 @@ isdir(outdir) && rm(outdir, recursive = true) # Run tests for TreeMesh @testset "TreeMesh" begin # Shallow water wet/dry 1D - @trixi_testset "1D-Test: elixir_shallowwater_well_balanced_nonperiodic.jl with wall boundary" begin + @trixi_testset "TreeMesh1D: elixir_shallowwater_well_balanced_nonperiodic.jl with wall boundary" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "tree_1d_dgsem", "elixir_shallowwater_well_balanced_nonperiodic.jl"), l2=[ @@ -42,6 +42,90 @@ isdir(outdir) && rm(outdir, recursive = true) @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 end end + # TODO: add upstream tests in 2D and positivity-preserving tests + + # Shallow water wet/dry 2D + # TreeMesh2D + @trixi_testset "TreeMesh2D: elixir_shallowwater_conical_island.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "tree_2d_dgsem", + "elixir_shallowwater_conical_island.jl"), + l2=[ + 0.0459315416430658, + 0.1644534881916991, + 0.16445348819169914, + 0.0011537702354532694, + ], + linf=[ + 0.21100717610846464, + 0.9501592344310412, + 0.9501592344310417, + 0.021790250683516282, + ], + tspan=(0.0, 0.025)) + # 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 + # Unstructured2D + @trixi_testset "Unstructured2D: elixir_shallowwater_wall_bc_shockcapturing.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "unstructured_2d_dgsem", + "elixir_shallowwater_wall_bc_shockcapturing.jl"), + l2=[ + 0.04444388691670699, + 0.1527771788033111, + 0.1593763537203512, + 6.225080476986749e-8, + ], + linf=[ + 0.6526506870169639, + 1.980765893182952, + 2.4807635459119757, + 3.982097158683473e-7, + ], + tspan=(0.0, 0.05), + surface_flux=(FluxHydrostaticReconstruction(FluxHLL(min_max_speed_naive), + hydrostatic_reconstruction_audusse_etal), + flux_nonconservative_audusse_etal)) + # 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 + # Structured2D + @trixi_testset "Structured2D: elixir_shallowwater_conical_island.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "structured_2d_dgsem", + "elixir_shallowwater_conical_island.jl"), + l2=[ + 0.04593154164306353, + 0.1644534881916908, + 0.16445348819169076, + 0.0011537702354532122, + ], + linf=[ + 0.21100717610846442, + 0.9501592344310412, + 0.950159234431041, + 0.021790250683516296, + ], + tspan=(0.0, 0.025)) + # 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 # Clean up afterwards: delete output directory