From 62610bb6a858eedfaf8fe61f074aa757441a08d1 Mon Sep 17 00:00:00 2001 From: Patrick Ersing <114223904+patrickersing@users.noreply.github.com> Date: Fri, 23 Feb 2024 05:30:08 +0100 Subject: [PATCH] Separation of TrixiShallowWater.jl (#1809) * remove wet_dry functionality for SWE-1D * remove wet_dry functionality for SWE-2D * remove twolayer equations * remove limiters from swe_quasi_1d * remove export of min_max_speed_chen_noelle * remove elixirs * update news.md * add unit tests to increase coverage * Change news.md according to code review Co-authored-by: Hendrik Ranocha --------- Co-authored-by: Andrew Winters Co-authored-by: Hendrik Ranocha --- NEWS.md | 13 + .../elixir_shallowwater_conical_island.jl | 114 --- .../elixir_shallowwater_parabolic_bowl.jl | 119 --- ...ixir_shallowwater_well_balanced_wet_dry.jl | 207 ----- .../elixir_shallowwater_beach.jl | 123 --- .../elixir_shallowwater_parabolic_bowl.jl | 119 --- ...lixir_shallowwater_twolayer_convergence.jl | 60 -- .../elixir_shallowwater_twolayer_dam_break.jl | 94 -- ...xir_shallowwater_twolayer_well_balanced.jl | 86 -- ...ixir_shallowwater_well_balanced_wet_dry.jl | 172 ---- .../elixir_shallowwater_conical_island.jl | 117 --- .../elixir_shallowwater_parabolic_bowl.jl | 121 --- ...lixir_shallowwater_twolayer_convergence.jl | 60 -- ...xir_shallowwater_twolayer_well_balanced.jl | 81 -- ...ixir_shallowwater_well_balanced_wet_dry.jl | 206 ----- ...ixir_shallowwater_three_mound_dam_break.jl | 133 --- ...lixir_shallowwater_twolayer_convergence.jl | 63 -- .../elixir_shallowwater_twolayer_dam_break.jl | 147 ---- ...xir_shallowwater_twolayer_well_balanced.jl | 81 -- src/Trixi.jl | 12 +- src/callbacks_stage/callbacks_stage.jl | 2 - .../positivity_shallow_water.jl | 89 -- .../positivity_shallow_water_dg1d.jl | 89 -- .../positivity_shallow_water_dg2d.jl | 90 -- src/equations/equations.jl | 2 - src/equations/numerical_fluxes.jl | 23 - src/equations/shallow_water_1d.jl | 187 +--- src/equations/shallow_water_2d.jl | 266 +----- src/equations/shallow_water_quasi_1d.jl | 37 +- src/equations/shallow_water_two_layer_1d.jl | 511 ----------- src/equations/shallow_water_two_layer_2d.jl | 805 ------------------ src/solvers/dgsem_tree/indicators.jl | 76 -- src/solvers/dgsem_tree/indicators_1d.jl | 109 --- src/solvers/dgsem_tree/indicators_2d.jl | 110 --- test/test_structured_2d.jl | 78 -- test/test_tree_1d.jl | 2 - test/test_tree_1d_shallowwater.jl | 75 -- test/test_tree_1d_shallowwater_twolayer.jl | 74 -- test/test_tree_2d_part3.jl | 3 - test/test_tree_2d_shallowwater.jl | 79 -- test/test_tree_2d_shallowwater_twolayer.jl | 88 -- test/test_unit.jl | 34 +- test/test_unstructured_2d.jl | 101 --- 43 files changed, 56 insertions(+), 5002 deletions(-) delete mode 100644 examples/structured_2d_dgsem/elixir_shallowwater_conical_island.jl delete mode 100644 examples/structured_2d_dgsem/elixir_shallowwater_parabolic_bowl.jl delete mode 100644 examples/structured_2d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl delete mode 100644 examples/tree_1d_dgsem/elixir_shallowwater_beach.jl delete mode 100644 examples/tree_1d_dgsem/elixir_shallowwater_parabolic_bowl.jl delete mode 100644 examples/tree_1d_dgsem/elixir_shallowwater_twolayer_convergence.jl delete mode 100644 examples/tree_1d_dgsem/elixir_shallowwater_twolayer_dam_break.jl delete mode 100644 examples/tree_1d_dgsem/elixir_shallowwater_twolayer_well_balanced.jl delete mode 100644 examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl delete mode 100644 examples/tree_2d_dgsem/elixir_shallowwater_conical_island.jl delete mode 100644 examples/tree_2d_dgsem/elixir_shallowwater_parabolic_bowl.jl delete mode 100644 examples/tree_2d_dgsem/elixir_shallowwater_twolayer_convergence.jl delete mode 100644 examples/tree_2d_dgsem/elixir_shallowwater_twolayer_well_balanced.jl delete mode 100644 examples/tree_2d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl delete mode 100644 examples/unstructured_2d_dgsem/elixir_shallowwater_three_mound_dam_break.jl delete mode 100644 examples/unstructured_2d_dgsem/elixir_shallowwater_twolayer_convergence.jl delete mode 100644 examples/unstructured_2d_dgsem/elixir_shallowwater_twolayer_dam_break.jl delete mode 100644 examples/unstructured_2d_dgsem/elixir_shallowwater_twolayer_well_balanced.jl delete mode 100644 src/callbacks_stage/positivity_shallow_water.jl delete mode 100644 src/callbacks_stage/positivity_shallow_water_dg1d.jl delete mode 100644 src/callbacks_stage/positivity_shallow_water_dg2d.jl delete mode 100644 src/equations/shallow_water_two_layer_1d.jl delete mode 100644 src/equations/shallow_water_two_layer_2d.jl delete mode 100644 test/test_tree_1d_shallowwater_twolayer.jl delete mode 100644 test/test_tree_2d_shallowwater_twolayer.jl diff --git a/NEWS.md b/NEWS.md index feccd1f9852..ecc91581e9a 100644 --- a/NEWS.md +++ b/NEWS.md @@ -4,6 +4,19 @@ Trixi.jl follows the interpretation of [semantic versioning (semver)](https://ju used in the Julia ecosystem. Notable changes will be documented in this file for human readability. +## Changes when updating to v0.7 from v0.6.x + +#### Added + +#### Changed + +#### Deprecated + +#### Removed +- Some specialized shallow water specific features are no longer available directly in + Trixi.jl, but are moved to a dedicated repository: [TrixiShallowWater.jl](https://github.com/trixi-framework/TrixiShallowWater.jl). This includes all features related to wetting and drying, as well as the `ShallowWaterTwoLayerEquations1D` and `ShallowWaterTwoLayerEquations2D`. + However, the basic shallow water equations are still part of Trixi.jl. We'll also be updating the TrixiShallowWater.jl documentation with instructions on how to use these relocated features in the future. + ## Changes in the v0.6 lifecycle #### Added diff --git a/examples/structured_2d_dgsem/elixir_shallowwater_conical_island.jl b/examples/structured_2d_dgsem/elixir_shallowwater_conical_island.jl deleted file mode 100644 index e65ed19221e..00000000000 --- a/examples/structured_2d_dgsem/elixir_shallowwater_conical_island.jl +++ /dev/null @@ -1,114 +0,0 @@ - -using OrdinaryDiffEq -using Trixi - -############################################################################### -# Semidiscretization of the shallow water equations -# -# TODO: TrixiShallowWater: wet/dry example elixir - -equations = ShallowWaterEquations2D(gravity_constant = 9.81, H0 = 1.4) - -""" - initial_condition_conical_island(x, t, equations::ShallowWaterEquations2D) - -Initial condition for the [`ShallowWaterEquations2D`](@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::ShallowWaterEquations2D) - # 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 ShallowWaterEquations 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 deleted file mode 100644 index bc198f18835..00000000000 --- a/examples/structured_2d_dgsem/elixir_shallowwater_parabolic_bowl.jl +++ /dev/null @@ -1,119 +0,0 @@ - -using OrdinaryDiffEq -using Trixi - -############################################################################### -# Semidiscretization of the shallow water equations -# -# TODO: TrixiShallowWater: wet/dry example elixir - -equations = ShallowWaterEquations2D(gravity_constant = 9.81) - -""" - initial_condition_parabolic_bowl(x, t, equations:: ShallowWaterEquations2D) - -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::ShallowWaterEquations2D) - 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 ShallowWaterEquations 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_well_balanced_wet_dry.jl b/examples/structured_2d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl deleted file mode 100644 index 8e492b1ba05..00000000000 --- a/examples/structured_2d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl +++ /dev/null @@ -1,207 +0,0 @@ - -using OrdinaryDiffEq -using Trixi -using Printf: @printf, @sprintf - -############################################################################### -# Semidiscretization of the shallow water equations -# -# TODO: TrixiShallowWater: wet/dry example elixir - -equations = ShallowWaterEquations2D(gravity_constant = 9.812) - -""" - initial_condition_well_balanced_chen_noelle(x, t, equations:: ShallowWaterEquations2D) - -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::ShallowWaterEquations2D) - 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 ShallowWaterEquations 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::ShallowWaterEquations2D) - 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 deleted file mode 100644 index 378079ca334..00000000000 --- a/examples/tree_1d_dgsem/elixir_shallowwater_beach.jl +++ /dev/null @@ -1,123 +0,0 @@ - -using OrdinaryDiffEq -using Trixi - -############################################################################### -# Semidiscretization of the shallow water equations -# -# TODO: TrixiShallowWater: wet/dry example elixir - -equations = ShallowWaterEquations1D(gravity_constant = 9.812) - -""" - initial_condition_beach(x, t, equations:: ShallowWaterEquations1D) -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::ShallowWaterEquations1D) - 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 ShallowWaterEquations 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 deleted file mode 100644 index a586562af7e..00000000000 --- a/examples/tree_1d_dgsem/elixir_shallowwater_parabolic_bowl.jl +++ /dev/null @@ -1,119 +0,0 @@ - -using OrdinaryDiffEq -using Trixi - -############################################################################### -# Semidiscretization of the shallow water equations -# -# TODO: TrixiShallowWater: wet/dry example elixir - -equations = ShallowWaterEquations1D(gravity_constant = 9.81) - -""" - initial_condition_parabolic_bowl(x, t, equations:: ShallowWaterEquations1D) - -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::ShallowWaterEquations1D) - 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 ShallowWaterEquations 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_twolayer_convergence.jl b/examples/tree_1d_dgsem/elixir_shallowwater_twolayer_convergence.jl deleted file mode 100644 index e6a01849852..00000000000 --- a/examples/tree_1d_dgsem/elixir_shallowwater_twolayer_convergence.jl +++ /dev/null @@ -1,60 +0,0 @@ - -using OrdinaryDiffEq -using Trixi - -############################################################################### -# Semidiscretization of the two-layer shallow water equations - -equations = ShallowWaterTwoLayerEquations1D(gravity_constant = 10.0, rho_upper = 0.9, - rho_lower = 1.0) - -initial_condition = initial_condition_convergence_test - -############################################################################### -# Get the DG approximation space - -volume_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal) -solver = DGSEM(polydeg = 3, - surface_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal), - volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) - -############################################################################### -# Get the TreeMesh and setup a periodic mesh - -coordinates_min = 0.0 -coordinates_max = 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 = 500, - 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, - save_everystep = false, callback = callbacks); -summary_callback() # print the timer summary diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_twolayer_dam_break.jl b/examples/tree_1d_dgsem/elixir_shallowwater_twolayer_dam_break.jl deleted file mode 100644 index 03b93754d0f..00000000000 --- a/examples/tree_1d_dgsem/elixir_shallowwater_twolayer_dam_break.jl +++ /dev/null @@ -1,94 +0,0 @@ - -using OrdinaryDiffEq -using Trixi - -############################################################################### -# Semidiscretization of the two-layer shallow water equations for a dam break -# test with a discontinuous bottom topography function to test entropy conservation - -equations = ShallowWaterTwoLayerEquations1D(gravity_constant = 9.81, H0 = 2.0, - rho_upper = 0.9, rho_lower = 1.0) - -# Initial condition of a dam break with a discontinuous water heights and bottom topography. -# Works as intended for TreeMesh1D with `initial_refinement_level=5`. 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_dam_break(x, t, equations::ShallowWaterTwoLayerEquations1D) - v1_upper = 0.0 - v1_lower = 0.0 - - # Set the discontinuity - if x[1] <= 10.0 - H_lower = 2.0 - H_upper = 4.0 - b = 0.0 - else - H_lower = 1.5 - H_upper = 3.0 - b = 0.5 - end - - return prim2cons(SVector(H_upper, v1_upper, H_lower, v1_lower, b), equations) -end - -initial_condition = initial_condition_dam_break - -############################################################################### -# Get the DG approximation space - -volume_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal) -solver = DGSEM(polydeg = 3, - surface_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal), - volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) - -############################################################################### -# Get the TreeMesh and setup a non-periodic mesh - -coordinates_min = 0.0 -coordinates_max = 20.0 -mesh = TreeMesh(coordinates_min, coordinates_max, - initial_refinement_level = 5, - n_cells_max = 10000, - periodicity = false) - -boundary_condition = boundary_condition_slip_wall - -# create the semidiscretization object -semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, - boundary_conditions = boundary_condition) - -############################################################################### -# ODE solvers - -tspan = (0.0, 0.4) -ode = semidiscretize(semi, tspan) - -############################################################################### -# Callbacks - -summary_callback = SummaryCallback() - -analysis_interval = 500 -analysis_callback = AnalysisCallback(semi, interval = analysis_interval, - save_analysis = false, - extra_analysis_integrals = (energy_total, - energy_kinetic, - energy_internal)) - -stepsize_callback = StepsizeCallback(cfl = 1.0) - -alive_callback = AliveCallback(analysis_interval = analysis_interval) - -save_solution = SaveSolutionCallback(interval = 500, - 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, - save_everystep = false, callback = callbacks); -summary_callback() # print the timer summary diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_twolayer_well_balanced.jl b/examples/tree_1d_dgsem/elixir_shallowwater_twolayer_well_balanced.jl deleted file mode 100644 index 098e3aaf601..00000000000 --- a/examples/tree_1d_dgsem/elixir_shallowwater_twolayer_well_balanced.jl +++ /dev/null @@ -1,86 +0,0 @@ - -using OrdinaryDiffEq -using Trixi - -############################################################################### -# Semidiscretization of the two-layer shallow water equations to test well-balancedness - -equations = ShallowWaterTwoLayerEquations1D(gravity_constant = 1.0, H0 = 0.6, - rho_upper = 0.9, rho_lower = 1.0) - -""" - initial_condition_fjordholm_well_balanced(x, t, equations::ShallowWaterTwoLayerEquations1D) - -Initial condition to test well balanced with a bottom topography from Fjordholm -- Ulrik Skre Fjordholm (2012) - Energy conservative and stable schemes for the two-layer shallow water equations. - [DOI: 10.1142/9789814417099_0039](https://doi.org/10.1142/9789814417099_0039) -""" -function initial_condition_fjordholm_well_balanced(x, t, - equations::ShallowWaterTwoLayerEquations1D) - inicenter = 0.5 - x_norm = x[1] - inicenter - r = abs(x_norm) - - H_lower = 0.5 - H_upper = 0.6 - v1_upper = 0.0 - v1_lower = 0.0 - b = r <= 0.1 ? 0.2 * (cos(10 * pi * (x[1] - 0.5)) + 1) : 0.0 - return prim2cons(SVector(H_upper, v1_upper, H_lower, v1_lower, b), equations) -end - -initial_condition = initial_condition_fjordholm_well_balanced - -############################################################################### -# Get the DG approximation space - -volume_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal) -solver = DGSEM(polydeg = 3, - surface_flux = (flux_es_ersing_etal, flux_nonconservative_ersing_etal), - volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) - -############################################################################### -# Get the TreeMesh and setup a periodic mesh - -coordinates_min = 0.0 -coordinates_max = 1.0 -mesh = TreeMesh(coordinates_min, coordinates_max, - initial_refinement_level = 4, - n_cells_max = 10_000, - periodicity = true) - -# 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 = (lake_at_rest_error,)) - -stepsize_callback = StepsizeCallback(cfl = 1.0) - -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, - 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_1d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl b/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl deleted file mode 100644 index 26a8960ab46..00000000000 --- a/examples/tree_1d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl +++ /dev/null @@ -1,172 +0,0 @@ - -using OrdinaryDiffEq -using Trixi -using Printf: @printf, @sprintf - -############################################################################### -# Semidiscretization of the shallow water equations -# -# TODO: TrixiShallowWater: wet/dry example elixir - -equations = ShallowWaterEquations1D(gravity_constant = 9.812) - -""" - initial_condition_complex_bottom_well_balanced(x, t, equations:: ShallowWaterEquations1D) - -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::ShallowWaterEquations1D) - 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 ShallowWaterEquations 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::ShallowWaterEquations1D) - 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 deleted file mode 100644 index 349b3741869..00000000000 --- a/examples/tree_2d_dgsem/elixir_shallowwater_conical_island.jl +++ /dev/null @@ -1,117 +0,0 @@ - -using OrdinaryDiffEq -using Trixi - -############################################################################### -# semidiscretization of the shallow water equations -# -# TODO: TrixiShallowWater: wet/dry example elixir - -equations = ShallowWaterEquations2D(gravity_constant = 9.81, H0 = 1.4) - -""" - initial_condition_conical_island(x, t, equations::ShallowWaterEquations2D) - -Initial condition for the [`ShallowWaterEquations2D`](@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::ShallowWaterEquations2D) - # 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 ShallowWaterEquations 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_parabolic_bowl.jl b/examples/tree_2d_dgsem/elixir_shallowwater_parabolic_bowl.jl deleted file mode 100644 index 2008019cc31..00000000000 --- a/examples/tree_2d_dgsem/elixir_shallowwater_parabolic_bowl.jl +++ /dev/null @@ -1,121 +0,0 @@ - -using OrdinaryDiffEq -using Trixi - -############################################################################### -# Semidiscretization of the shallow water equations -# -# TODO: TrixiShallowWater: wet/dry example elixir - -equations = ShallowWaterEquations2D(gravity_constant = 9.81) - -""" - initial_condition_parabolic_bowl(x, t, equations:: ShallowWaterEquations2D) - -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::ShallowWaterEquations2D) - 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 ShallowWaterEquations 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) - -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_twolayer_convergence.jl b/examples/tree_2d_dgsem/elixir_shallowwater_twolayer_convergence.jl deleted file mode 100644 index 790916e4467..00000000000 --- a/examples/tree_2d_dgsem/elixir_shallowwater_twolayer_convergence.jl +++ /dev/null @@ -1,60 +0,0 @@ - -using OrdinaryDiffEq -using Trixi - -############################################################################### -# Semidiscretization of the two-layer shallow water equations - -equations = ShallowWaterTwoLayerEquations2D(gravity_constant = 10.0, rho_upper = 0.9, - rho_lower = 1.0) - -initial_condition = initial_condition_convergence_test - -############################################################################### -# Get the DG approximation space - -volume_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal) -solver = DGSEM(polydeg = 3, - surface_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_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 = 20_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 = 500, - 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, - save_everystep = false, callback = callbacks); -summary_callback() # print the timer summary diff --git a/examples/tree_2d_dgsem/elixir_shallowwater_twolayer_well_balanced.jl b/examples/tree_2d_dgsem/elixir_shallowwater_twolayer_well_balanced.jl deleted file mode 100644 index 264c26390fe..00000000000 --- a/examples/tree_2d_dgsem/elixir_shallowwater_twolayer_well_balanced.jl +++ /dev/null @@ -1,81 +0,0 @@ - -using OrdinaryDiffEq -using Trixi - -############################################################################### -# Semidiscretization of the two-layer shallow water equations with a bottom topography function -# to test well-balancedness - -equations = ShallowWaterTwoLayerEquations2D(gravity_constant = 9.81, H0 = 0.6, - rho_upper = 0.9, rho_lower = 1.0) - -# An initial condition with constant total water height, zero velocities and a bottom topography to -# test well-balancedness -function initial_condition_well_balanced(x, t, equations::ShallowWaterTwoLayerEquations2D) - H_lower = 0.5 - H_upper = 0.6 - v1_upper = 0.0 - v2_upper = 0.0 - v1_lower = 0.0 - v2_lower = 0.0 - b = (((x[1] - 0.5)^2 + (x[2] - 0.5)^2) < 0.04 ? - 0.2 * (cos(4 * pi * sqrt((x[1] - 0.5)^2 + (x[2] + - -0.5)^2)) + 1) : 0.0) - - return prim2cons(SVector(H_upper, v1_upper, v2_upper, H_lower, v1_lower, v2_lower, b), - equations) -end - -initial_condition = initial_condition_well_balanced - -############################################################################### -# Get the DG approximation space - -volume_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal) -surface_flux = (flux_es_ersing_etal, flux_nonconservative_ersing_etal) -solver = DGSEM(polydeg = 3, surface_flux = surface_flux, - volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) - -############################################################################### -# Get the TreeMesh and setup a periodic mesh - -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, - periodicity = true) - -# Create the semi discretization object -semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) - -############################################################################### -# ODE solver - -tspan = (0.0, 10.0) -ode = semidiscretize(semi, tspan) - -summary_callback = SummaryCallback() - -analysis_interval = 1000 -analysis_callback = AnalysisCallback(semi, interval = analysis_interval, - extra_analysis_integrals = (lake_at_rest_error,)) - -stepsize_callback = StepsizeCallback(cfl = 1.0) - -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, - 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 deleted file mode 100644 index 034411c2b54..00000000000 --- a/examples/tree_2d_dgsem/elixir_shallowwater_well_balanced_wet_dry.jl +++ /dev/null @@ -1,206 +0,0 @@ - -using OrdinaryDiffEq -using Trixi -using Printf: @printf, @sprintf - -############################################################################### -# Semidiscretization of the shallow water equations -# -# TODO: TrixiShallowWater: wet/dry example elixir - -equations = ShallowWaterEquations2D(gravity_constant = 9.812) - -""" - initial_condition_well_balanced_chen_noelle(x, t, equations:: ShallowWaterEquations2D) - -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::ShallowWaterEquations2D) - 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 ShallowWaterEquations 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::ShallowWaterEquations2D) - 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_three_mound_dam_break.jl b/examples/unstructured_2d_dgsem/elixir_shallowwater_three_mound_dam_break.jl deleted file mode 100644 index df321aad267..00000000000 --- a/examples/unstructured_2d_dgsem/elixir_shallowwater_three_mound_dam_break.jl +++ /dev/null @@ -1,133 +0,0 @@ - -using OrdinaryDiffEq -using Trixi - -############################################################################### -# semidiscretization of the shallow water equations -# -# TODO: TrixiShallowWater: wet/dry example elixir - -equations = ShallowWaterEquations2D(gravity_constant = 9.81, H0 = 1.875, - threshold_limiter = 1e-12, threshold_wet = 1e-14) - -""" - initial_condition_three_mounds(x, t, equations::ShallowWaterEquations2D) - -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::ShallowWaterEquations2D) - - # 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 ShallowWaterEquations 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::ShallowWaterEquations2D) - # 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) -mesh_file = Trixi.download("https://gist.githubusercontent.com/svengoldberg/c3c87fecb3fc6e46be7f0d1c7cb35f83/raw/e817ecd9e6c4686581d63c46128f9b6468d396d3/mesh_three_mound.mesh", - joinpath(@__DIR__, "mesh_three_mound.mesh")) - -mesh = UnstructuredMesh2D(mesh_file) - -# Create the semi discretization object -semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver; - boundary_conditions = boundary_conditions) - -############################################################################### -# ODE solver - -tspan = (0.0, 20.0) -ode = semidiscretize(semi, tspan) - -############################################################################### -# Callbacks - -summary_callback = SummaryCallback() - -analysis_interval = 1000 -analysis_callback = AnalysisCallback(semi, interval = analysis_interval) - -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_twolayer_convergence.jl b/examples/unstructured_2d_dgsem/elixir_shallowwater_twolayer_convergence.jl deleted file mode 100644 index fcc08b6f991..00000000000 --- a/examples/unstructured_2d_dgsem/elixir_shallowwater_twolayer_convergence.jl +++ /dev/null @@ -1,63 +0,0 @@ - -using OrdinaryDiffEq -using Trixi - -############################################################################### -# Semidiscretization of the two-layer shallow water equations with a periodic -# bottom topography function (set in the initial conditions) - -equations = ShallowWaterTwoLayerEquations2D(gravity_constant = 10.0, rho_upper = 0.9, - rho_lower = 1.0) - -initial_condition = initial_condition_convergence_test - -############################################################################### -# Get the DG approximation space - -volume_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal) -surface_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_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) -mesh_file = Trixi.download("https://gist.githubusercontent.com/andrewwinters5000/8f8cd23df27fcd494553f2a89f3c1ba4/raw/85e3c8d976bbe57ca3d559d653087b0889535295/mesh_alfven_wave_with_twist_and_flip.mesh", - joinpath(@__DIR__, "mesh_alfven_wave_with_twist_and_flip.mesh")) - -mesh = UnstructuredMesh2D(mesh_file, periodicity = true) - -# Create the semidiscretization 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 = 500, - 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_twolayer_dam_break.jl b/examples/unstructured_2d_dgsem/elixir_shallowwater_twolayer_dam_break.jl deleted file mode 100644 index 821f31c52ac..00000000000 --- a/examples/unstructured_2d_dgsem/elixir_shallowwater_twolayer_dam_break.jl +++ /dev/null @@ -1,147 +0,0 @@ - -using OrdinaryDiffEq -using Trixi - -############################################################################### -# Semidiscretization of the two-layer shallow water equations for a dam break test with a -# discontinuous bottom topography function to test energy conservation - -equations = ShallowWaterTwoLayerEquations2D(gravity_constant = 1.0, rho_upper = 0.9, - rho_lower = 1.0) - -# This test case uses a special work around to setup a truly discontinuous bottom topography -# function and initial condition for this academic testcase of entropy conservation. First, a -# dummy initial_condition_dam_break is introduced to create the semidiscretization. Then the initial -# condition is reset with the true discontinuous values from initial_condition_discontinuous_dam_break. - -function initial_condition_dam_break(x, t, equations::ShallowWaterTwoLayerEquations2D) - if x[1] < sqrt(2) / 2 - H_upper = 1.0 - H_lower = 0.6 - b = 0.1 - else - H_upper = 0.9 - H_lower = 0.5 - b = 0.0 - end - - v1_upper = 0.0 - v2_upper = 0.0 - v1_lower = 0.0 - v2_lower = 0.0 - return prim2cons(SVector(H_upper, v1_upper, v2_upper, H_lower, v1_lower, v2_lower, b), - equations) -end - -initial_condition = initial_condition_dam_break - -boundary_condition_constant = BoundaryConditionDirichlet(initial_condition_dam_break) - -############################################################################### -# Get the DG approximation space - -volume_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal) -surface_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal) -solver = DGSEM(polydeg = 6, surface_flux = surface_flux, - volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) - -############################################################################### -# Get the unstructured quad mesh from a file (downloads the file if not available locally) -mesh_file = Trixi.download("https://gist.githubusercontent.com/andrewwinters5000/8f8cd23df27fcd494553f2a89f3c1ba4/raw/85e3c8d976bbe57ca3d559d653087b0889535295/mesh_alfven_wave_with_twist_and_flip.mesh", - joinpath(@__DIR__, "mesh_alfven_wave_with_twist_and_flip.mesh")) - -mesh = UnstructuredMesh2D(mesh_file, periodicity = false) - -# Boundary conditions -boundary_condition = Dict(:Top => boundary_condition_slip_wall, - :Left => boundary_condition_slip_wall, - :Right => boundary_condition_slip_wall, - :Bottom => boundary_condition_slip_wall) - -# Create the semi discretization object -semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, - solver, boundary_conditions = boundary_condition) - -############################################################################### -# ODE solver - -tspan = (0.0, 0.5) -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 -# test case and initial condition. -# In contrast to the usual signature of initial conditions, this one get passed the -# `element_id` explicitly. In particular, this initial conditions works as intended -# only for the specific mesh loaded above! - -function initial_condition_discontinuous_dam_break(x, t, element_id, - equations::ShallowWaterTwoLayerEquations2D) - # Constant values - v1_upper = 0.0 - v2_upper = 0.0 - v1_lower = 0.0 - v2_lower = 0.0 - - # Left side of discontinuity - IDs = [1, 2, 5, 6, 9, 10, 13, 14] - if element_id in IDs - H_upper = 1.0 - H_lower = 0.6 - b = 0.0 - # Right side of discontinuity - else - H_upper = 0.9 - H_lower = 0.5 - b = 0.1 - end - - return prim2cons(SVector(H_upper, v1_upper, v2_upper, H_lower, v1_lower, v2_lower, b), - equations) -end - -# point to the data we want to augment -u = Trixi.wrap_array(ode.u0, semi) -# reset the initial condition -for element in eachelement(semi.solver, semi.cache) - for j in eachnode(semi.solver), i in eachnode(semi.solver) - x_node = Trixi.get_node_coords(semi.cache.elements.node_coordinates, equations, - semi.solver, i, j, element) - u_node = initial_condition_discontinuous_dam_break(x_node, first(tspan), element, - equations) - Trixi.set_node_vars!(u, u_node, equations, semi.solver, i, j, element) - end -end - -############################################################################### -# Callbacks - -summary_callback = SummaryCallback() - -analysis_interval = 500 -analysis_callback = AnalysisCallback(semi, interval = analysis_interval, - save_analysis = false, - extra_analysis_integrals = (energy_total, - energy_kinetic, - energy_internal)) - -alive_callback = AliveCallback(analysis_interval = analysis_interval) - -save_solution = SaveSolutionCallback(interval = 500, - 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_twolayer_well_balanced.jl b/examples/unstructured_2d_dgsem/elixir_shallowwater_twolayer_well_balanced.jl deleted file mode 100644 index ca1f54595bb..00000000000 --- a/examples/unstructured_2d_dgsem/elixir_shallowwater_twolayer_well_balanced.jl +++ /dev/null @@ -1,81 +0,0 @@ - -using OrdinaryDiffEq -using Trixi - -############################################################################### -# Semidiscretization of the two-layer shallow water equations with a discontinuous bottom -# topography to test well-balancedness - -equations = ShallowWaterTwoLayerEquations2D(gravity_constant = 1.0, H0 = 0.6, - rho_upper = 0.9, rho_lower = 1.0) - -# An initial condition with constant total water height, zero velocities and a bottom topography to -# test well-balancedness -function initial_condition_well_balanced(x, t, equations::ShallowWaterTwoLayerEquations2D) - H_lower = 0.5 - H_upper = 0.6 - v1_upper = 0.0 - v2_upper = 0.0 - v1_lower = 0.0 - v2_lower = 0.0 - - # Bottom Topography - b = (((x[1] - 0.5)^2 + (x[2] - 0.5)^2) < 0.04 ? - 0.2 * (cos(4 * pi * sqrt((x[1] - 0.5)^2 + (x[2] + - -0.5)^2)) + 1) : 0.0) - return prim2cons(SVector(H_upper, v1_upper, v2_upper, H_lower, v1_lower, v2_lower, b), - equations) -end - -initial_condition = initial_condition_well_balanced - -############################################################################### -# Get the DG approximation space - -volume_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal) -surface_flux = (flux_es_ersing_etal, flux_nonconservative_ersing_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) -mesh_file = Trixi.download("https://gist.githubusercontent.com/andrewwinters5000/8f8cd23df27fcd494553f2a89f3c1ba4/raw/85e3c8d976bbe57ca3d559d653087b0889535295/mesh_alfven_wave_with_twist_and_flip.mesh", - joinpath(@__DIR__, "mesh_alfven_wave_with_twist_and_flip.mesh")) - -mesh = UnstructuredMesh2D(mesh_file, periodicity = true) - -# Create the semi discretization object -semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) - -############################################################################### -# ODE solver - -tspan = (0.0, 10.0) -ode = semidiscretize(semi, tspan) - -summary_callback = SummaryCallback() - -analysis_interval = 1000 -analysis_callback = AnalysisCallback(semi, interval = analysis_interval, - extra_analysis_integrals = (lake_at_rest_error,)) - -stepsize_callback = StepsizeCallback(cfl = 1.0) - -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, - 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/Trixi.jl b/src/Trixi.jl index b7f7767a9d8..5f8cd9cae8e 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -160,7 +160,6 @@ export AcousticPerturbationEquations2D, InviscidBurgersEquation1D, LatticeBoltzmannEquations2D, LatticeBoltzmannEquations3D, ShallowWaterEquations1D, ShallowWaterEquations2D, - ShallowWaterTwoLayerEquations1D, ShallowWaterTwoLayerEquations2D, ShallowWaterEquationsQuasi1D, LinearizedEulerEquations2D, PolytropicEulerEquations2D, @@ -179,16 +178,12 @@ export flux, flux_central, flux_lax_friedrichs, flux_hll, flux_hllc, flux_hlle, flux_kennedy_gruber, flux_shima_etal, flux_ec, flux_fjordholm_etal, flux_nonconservative_fjordholm_etal, flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal, - flux_es_ersing_etal, flux_nonconservative_ersing_etal, + flux_nonconservative_ersing_etal, flux_chan_etal, flux_nonconservative_chan_etal, flux_winters_etal, hydrostatic_reconstruction_audusse_etal, flux_nonconservative_audusse_etal, -# TODO: TrixiShallowWater: move anything with "chen_noelle" to new file - hydrostatic_reconstruction_chen_noelle, flux_nonconservative_chen_noelle, - flux_hll_chen_noelle, FluxPlusDissipation, DissipationGlobalLaxFriedrichs, DissipationLocalLaxFriedrichs, FluxLaxFriedrichs, max_abs_speed_naive, FluxHLL, min_max_speed_naive, min_max_speed_davis, min_max_speed_einfeldt, - min_max_speed_chen_noelle, FluxLMARS, FluxRotated, flux_shima_etal_turbo, flux_ranocha_turbo, @@ -239,8 +234,6 @@ export DG, VolumeIntegralFluxDifferencing, VolumeIntegralPureLGLFiniteVolume, VolumeIntegralShockCapturingHG, IndicatorHennemannGassner, -# TODO: TrixiShallowWater: move new indicator - IndicatorHennemannGassnerShallowWater, VolumeIntegralUpwind, SurfaceIntegralWeakForm, SurfaceIntegralStrongForm, SurfaceIntegralUpwind, @@ -276,8 +269,7 @@ export load_mesh, load_time, load_timestep, load_timestep!, load_dt, export ControllerThreeLevel, ControllerThreeLevelCombined, IndicatorLöhner, IndicatorLoehner, IndicatorMax -# TODO: TrixiShallowWater: move new limiter -export PositivityPreservingLimiterZhangShu, PositivityPreservingLimiterShallowWater +export PositivityPreservingLimiterZhangShu export trixi_include, examples_dir, get_examples, default_example, default_example_unstructured, ode_default_options diff --git a/src/callbacks_stage/callbacks_stage.jl b/src/callbacks_stage/callbacks_stage.jl index 70d60de7914..d5abc1d227d 100644 --- a/src/callbacks_stage/callbacks_stage.jl +++ b/src/callbacks_stage/callbacks_stage.jl @@ -8,6 +8,4 @@ include("positivity_zhang_shu.jl") include("subcell_limiter_idp_correction.jl") include("subcell_bounds_check.jl") -# TODO: TrixiShallowWater: move specific limiter file -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 deleted file mode 100644 index 36276026fe9..00000000000 --- a/src/callbacks_stage/positivity_shallow_water.jl +++ /dev/null @@ -1,89 +0,0 @@ -# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). -# Since these FMAs can increase the performance of many numerical algorithms, -# we need to opt-in explicitly. -# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. -@muladd begin -#! format: noindent - -# TODO: TrixiShallowWater: generic wet/dry limiter - -""" - 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 [`ShallowWaterEquations1D`](@ref) struct -or the [`ShallowWaterEquations2D`](@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 -`ShallowWaterEquations` 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::AbstractSemidiscretization, - t) - u = wrap_array(u_ode, semi) - @trixi_timeit timer() "positivity-preserving limiter" limiter_shallow_water!(u, - limiter!.variables, - 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{ShallowWaterEquations1D, - ShallowWaterEquations2D}, - 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{ShallowWaterEquations1D, - ShallowWaterEquations2D}, - 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 deleted file mode 100644 index 13c6866e895..00000000000 --- a/src/callbacks_stage/positivity_shallow_water_dg1d.jl +++ /dev/null @@ -1,89 +0,0 @@ -# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). -# Since these FMAs can increase the performance of many numerical algorithms, -# we need to opt-in explicitly. -# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. -@muladd begin -#! format: noindent - -# TODO: TrixiShallowWater: 1D wet/dry limiter should move - -function limiter_shallow_water!(u, threshold::Real, variable, - mesh::AbstractMesh{1}, - equations::ShallowWaterEquations1D, - dg::DGSEM, cache) - @unpack weights = dg.basis - - @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 - @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 deleted file mode 100644 index da3a25fdcf4..00000000000 --- a/src/callbacks_stage/positivity_shallow_water_dg2d.jl +++ /dev/null @@ -1,90 +0,0 @@ -# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). -# Since these FMAs can increase the performance of many numerical algorithms, -# we need to opt-in explicitly. -# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. -@muladd begin -#! format: noindent - -# TODO: TrixiShallowWater: 2D wet/dry limiter should move - -function limiter_shallow_water!(u, threshold::Real, variable, - mesh::AbstractMesh{2}, - equations::ShallowWaterEquations2D, dg::DGSEM, cache) - @unpack weights = dg.basis - - @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 - @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 65875a2a7e5..8f476cf6f16 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -405,8 +405,6 @@ abstract type AbstractShallowWaterEquations{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end include("shallow_water_1d.jl") include("shallow_water_2d.jl") -include("shallow_water_two_layer_1d.jl") -include("shallow_water_two_layer_2d.jl") include("shallow_water_quasi_1d.jl") # CompressibleEulerEquations diff --git a/src/equations/numerical_fluxes.jl b/src/equations/numerical_fluxes.jl index 44d523b6e89..87fcb412244 100644 --- a/src/equations/numerical_fluxes.jl +++ b/src/equations/numerical_fluxes.jl @@ -326,29 +326,6 @@ This is a [`FluxHLL`](@ref)-type two-wave solver with special estimates of the w """ const flux_hlle = FluxHLL(min_max_speed_einfeldt) -# TODO: TrixiShallowWater: move the chen_noelle flux structure to the new package - -# 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_1d.jl` or `shallow_water_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) - """ flux_shima_etal_turbo(u_ll, u_rr, orientation_or_normal_direction, equations) diff --git a/src/equations/shallow_water_1d.jl b/src/equations/shallow_water_1d.jl index 25ce0fa79fe..e348ef946b7 100644 --- a/src/equations/shallow_water_1d.jl +++ b/src/equations/shallow_water_1d.jl @@ -6,7 +6,7 @@ #! format: noindent @doc raw""" - ShallowWaterEquations1D(; gravity, H0 = 0, threshold_limiter = nothing threshold_wet = nothing) + ShallowWaterEquations1D(; gravity, H0 = 0) Shallow water equations (SWE) in one space dimension. The equations are given by ```math @@ -24,12 +24,6 @@ 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)`` 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. @@ -51,35 +45,16 @@ References for the SWE are many but a good introduction is available in Chapter [DOI: 10.1017/CBO9780511791253](https://doi.org/10.1017/CBO9780511791253) """ struct ShallowWaterEquations1D{RealT <: Real} <: AbstractShallowWaterEquations{1, 3} - # TODO: TrixiShallowWater: where should the `threshold_limiter` and `threshold_wet` live? - # how to "properly" export these constants across the two packages? 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 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 ShallowWaterEquations1D(; 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 - ShallowWaterEquations1D(gravity_constant, H0, threshold_limiter, threshold_wet) +function ShallowWaterEquations1D(; gravity_constant, H0 = zero(gravity_constant)) + ShallowWaterEquations1D(gravity_constant, H0) end have_nonconservative_terms(::ShallowWaterEquations1D) = True() @@ -332,54 +307,6 @@ Further details on the hydrostatic reconstruction and its motivation can be foun z) end -# TODO: TrixiShallowWater: move wet/dry specific routine -""" - flux_nonconservative_chen_noelle(u_ll, u_rr, - orientation::Integer, - equations::ShallowWaterEquations1D) - -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. - -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::ShallowWaterEquations1D) - - # 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 - 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, equations::ShallowWaterEquations1D) @@ -521,67 +448,6 @@ Further details on this hydrostatic reconstruction and its motivation can be fou return u_ll_star, u_rr_star end -# TODO: TrixiShallowWater: move wet/dry specific routine -""" - hydrostatic_reconstruction_chen_noelle(u_ll, u_rr, orientation::Integer, - equations::ShallowWaterEquations1D) - -A particular type of hydrostatic reconstruction of the water height to guarantee well-balancedness -for a general bottom topography of the [`ShallowWaterEquations1D`](@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::ShallowWaterEquations1D) - # 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 `ShallowWaterEquations1D` 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 - # Calculate maximum wave speed for local Lax-Friedrichs-type dissipation as the # maximum velocity magnitude plus the maximum speed of sound @inline function max_abs_speed_naive(u_ll, u_rr, orientation::Integer, @@ -646,39 +512,6 @@ end return λ_min, λ_max end -# TODO: TrixiShallowWater: move wet/dry specific routine -""" - 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. - -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::ShallowWaterEquations1D) - # 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) - - 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))) - - return λ_min, λ_max -end - # More refined estimates for minimum and maximum wave speeds for HLL-type fluxes @inline function min_max_speed_davis(u_ll, u_rr, orientation::Integer, equations::ShallowWaterEquations1D) @@ -841,20 +674,10 @@ end 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. -# -# TODO: TrixiShallowWater: where should `threshold_limiter` live? May need -# to modify or have different versions of the `lake_at_rest_error` function +# be a constant value over time. @inline function lake_at_rest_error(u, equations::ShallowWaterEquations1D) 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)) + return abs(equations.H0 - (h + b)) end end # @muladd diff --git a/src/equations/shallow_water_2d.jl b/src/equations/shallow_water_2d.jl index 6728d7d5553..74a299a51e6 100644 --- a/src/equations/shallow_water_2d.jl +++ b/src/equations/shallow_water_2d.jl @@ -6,7 +6,7 @@ #! format: noindent @doc raw""" - ShallowWaterEquations2D(; gravity, H0 = 0, threshold_limiter = nothing, threshold_wet = nothing) + ShallowWaterEquations2D(; gravity, H0 = 0) Shallow water equations (SWE) in two space dimensions. The equations are given by ```math @@ -27,12 +27,6 @@ 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. @@ -54,18 +48,8 @@ References for the SWE are many but a good introduction is available in Chapter [DOI: 10.1017/CBO9780511791253](https://doi.org/10.1017/CBO9780511791253) """ struct ShallowWaterEquations2D{RealT <: Real} <: AbstractShallowWaterEquations{2, 4} - # TODO: TrixiShallowWater: where should the `threshold_limiter` and `threshold_wet` live? - # how to "properly" export these constants across the two packages? 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 end # Allow for flexibility to set the gravitational constant within an elixir depending on the @@ -73,16 +57,8 @@ end # 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 ShallowWaterEquations2D(; 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 - ShallowWaterEquations2D(gravity_constant, H0, threshold_limiter, threshold_wet) +function ShallowWaterEquations2D(; gravity_constant, H0 = zero(gravity_constant)) + ShallowWaterEquations2D(gravity_constant, H0) end have_nonconservative_terms(::ShallowWaterEquations2D) = True() @@ -460,69 +436,6 @@ Further details for the hydrostatic reconstruction and its motivation can be fou return u_ll_star, u_rr_star end -# TODO: TrixiShallowWater: move wet/dry specific routine -""" - hydrostatic_reconstruction_chen_noelle(u_ll, u_rr, orientation::Integer, - equations::ShallowWaterEquations2D) - -A particular type of hydrostatic reconstruction of the water height to guarantee well-balancedness -for a general bottom topography of the [`ShallowWaterEquations2D`](@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::ShallowWaterEquations2D) - # 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 `ShallowWaterEquations2D` 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::ShallowWaterEquations2D) @@ -608,104 +521,6 @@ end return SVector(f1, f2, f3, f4) end -# TODO: TrixiShallowWater: move wet/dry specific routine -""" - flux_nonconservative_chen_noelle(u_ll, u_rr, - orientation::Integer, - equations::ShallowWaterEquations2D) - flux_nonconservative_chen_noelle(u_ll, u_rr, - normal_direction_ll ::AbstractVector, - normal_direction_average ::AbstractVector, - equations::ShallowWaterEquations2D) - -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::ShallowWaterEquations2D) - # 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::ShallowWaterEquations2D) - # 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::ShallowWaterEquations2D) @@ -1020,67 +835,6 @@ end return λ_min, λ_max end -# TODO: TrixiShallowWater: move wet/dry specific routine -""" - 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::ShallowWaterEquations2D) - h_ll = waterheight(u_ll, equations) - v1_ll, v2_ll = velocity(u_ll, equations) - h_rr = 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::ShallowWaterEquations2D) - h_ll = waterheight(u_ll, equations) - v1_ll, v2_ll = velocity(u_ll, equations) - h_rr = 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 - # More refined estimates for minimum and maximum wave speeds for HLL-type fluxes @inline function min_max_speed_davis(u_ll, u_rr, orientation::Integer, equations::ShallowWaterEquations2D) @@ -1327,20 +1081,10 @@ end 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. -# -# TODO: TrixiShallowWater: where should `threshold_limiter` live? May need -# to modify or have different versions of the `lake_at_rest_error` function +# be a constant value over time. @inline function lake_at_rest_error(u, equations::ShallowWaterEquations2D) 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)) + return abs(equations.H0 - (h + b)) end end # @muladd diff --git a/src/equations/shallow_water_quasi_1d.jl b/src/equations/shallow_water_quasi_1d.jl index d52fbab841d..51c360104a7 100644 --- a/src/equations/shallow_water_quasi_1d.jl +++ b/src/equations/shallow_water_quasi_1d.jl @@ -22,12 +22,6 @@ The gravitational constant is denoted by `g`, the (possibly) variable bottom top 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)`` and channel width ``a(x)`` are 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 and ``a`` to one. @@ -47,14 +41,6 @@ struct ShallowWaterEquationsQuasi1D{RealT <: Real} <: AbstractShallowWaterEquations{1, 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 end # Allow for flexibility to set the gravitational constant within an elixir depending on the @@ -62,17 +48,8 @@ end # 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 ShallowWaterEquationsQuasi1D(; 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 - ShallowWaterEquationsQuasi1D(gravity_constant, H0, threshold_limiter, threshold_wet) +function ShallowWaterEquationsQuasi1D(; gravity_constant, H0 = zero(gravity_constant)) + ShallowWaterEquationsQuasi1D(gravity_constant, H0) end have_nonconservative_terms(::ShallowWaterEquationsQuasi1D) = True() @@ -338,18 +315,10 @@ end # be a constant value over time. Note, assumes there is a single reference # water height `H0` with which to compare. # -# TODO: TrixiShallowWater: where should `threshold_limiter` live? May need -# to modify or have different versions of the `lake_at_rest_error` function @inline function lake_at_rest_error(u, equations::ShallowWaterEquationsQuasi1D) _, _, b, _ = u h = waterheight(u, equations) - # 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)) + return abs(equations.H0 - (h + b)) end end # @muladd diff --git a/src/equations/shallow_water_two_layer_1d.jl b/src/equations/shallow_water_two_layer_1d.jl deleted file mode 100644 index 42ff393593e..00000000000 --- a/src/equations/shallow_water_two_layer_1d.jl +++ /dev/null @@ -1,511 +0,0 @@ -# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). -# Since these FMAs can increase the performance of many numerical algorithms, -# we need to opt-in explicitly. -# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. -@muladd begin -#! format: noindent - -# TODO: TrixiShallowWater: 1D two layer equations should move to new package - -@doc raw""" - ShallowWaterTwoLayerEquations1D(gravity, H0, rho_upper, rho_lower) - -Two-Layer Shallow Water equations (2LSWE) in one space dimension. The equations are given by -```math -\begin{alignat*}{4} -&\frac{\partial}{\partial t}h_{upper} -&&+ \frac{\partial}{\partial x}\left(h_{upper} v_{1,upper}\right) -&&= 0 \\ -&\frac{\partial}{\partial t}\left(h_{upper}v_{1,upper}\right) -&&+ \frac{\partial}{\partial x}\left(h_{upper}v_{1,upper}^2 + \dfrac{gh_{upper}^2}{2}\right) -&&= -gh_{upper}\frac{\partial}{\partial x}\left(b+h_{lower}\right)\\ -&\frac{\partial}{\partial t}h_{lower} -&&+ \frac{\partial}{\partial x}\left(h_{lower}v_{1,lower}\right) -&&= 0 \\ -&\frac{\partial}{\partial t}\left(h_{lower}v_{1,lower}\right) -&&+ \frac{\partial}{\partial x}\left(h_{lower}v_{1,lower}^2 + \dfrac{gh_{lower}^2}{2}\right) -&&= -gh_{lower}\frac{\partial}{\partial x}\left(b+\dfrac{\rho_{upper}}{\rho_{lower}}h_{upper}\right). -\end{alignat*} -``` -The unknown quantities of the 2LSWE are the water heights of the {lower} layer ``h_{lower}`` and the -{upper} layer ``h_{upper}`` with respective velocities ``v_{1,upper}`` and ``v_{1,lower}``. The gravitational constant is -denoted by `g`, the layer densitites by ``\rho_{upper}``and ``\rho_{lower}`` and the (possibly) variable -bottom topography function ``b(x)``. The conservative variable water height ``h_{lower}`` is measured -from the bottom topography ``b`` and ``h_{upper}`` relative to ``h_{lower}``, therefore one also defines the -total water heights as ``H_{upper} = h_{upper} + h_{upper} + b`` and ``H_{lower} = h_{lower} + b``. - -The densities must be chosen such that ``\rho_{upper} < \rho_{lower}``, to make sure that the heavier fluid -``\rho_{lower}`` is in the bottom layer and the lighter fluid ``\rho_{upper}`` in the {upper} layer. - -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. - -The bottom topography function ``b(x)`` is set inside the initial condition routine -for a particular problem setup. - -In addition to the unknowns, Trixi 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's visualization tools will visualize the bottom topography by default. - -A good introduction for the 2LSWE is available in Chapter 12 of the book: -- Benoit Cushman-Roisin (2011)\ - Introduction to geophyiscal fluid dynamics: physical and numerical aspects\ - \ - ISBN: 978-0-12-088759-0 -""" -struct ShallowWaterTwoLayerEquations1D{RealT <: Real} <: - AbstractShallowWaterEquations{1, 5} - gravity::RealT # gravitational constant - H0::RealT # constant "lake-at-rest" total water height - rho_upper::RealT # lower layer density - rho_lower::RealT # upper layer density - r::RealT # ratio of rho_upper / rho_lower -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. Densities must be specified such that rho_upper <= rho_lower. -function ShallowWaterTwoLayerEquations1D(; gravity_constant, - H0 = zero(gravity_constant), rho_upper, - rho_lower) - # Assign density ratio if rho_upper <= rho_lower - if rho_upper > rho_lower - error("Invalid input: Densities must be chosen such that rho_upper <= rho_lower") - else - r = rho_upper / rho_lower - end - ShallowWaterTwoLayerEquations1D(gravity_constant, H0, rho_upper, rho_lower, r) -end - -have_nonconservative_terms(::ShallowWaterTwoLayerEquations1D) = True() -function varnames(::typeof(cons2cons), ::ShallowWaterTwoLayerEquations1D) - ("h_upper", "h_v_upper", - "h_lower", "h_v_lower", "b") -end -# Note, we use the total water height, H_lower = h_upper + h_lower + b, and first layer total height -# H_upper = h_upper + b as the first primitive variable for easier visualization and setting initial -# conditions -function varnames(::typeof(cons2prim), ::ShallowWaterTwoLayerEquations1D) - ("H_upper", "v_upper", - "H_lower", "v_lower", "b") -end - -# Set initial conditions at physical location `x` for time `t` -""" - initial_condition_convergence_test(x, t, equations::ShallowWaterTwoLayerEquations1D) - -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 initial_condition_convergence_test(x, t, - equations::ShallowWaterTwoLayerEquations1D) - # some constants are chosen such that the function is periodic on the domain [0,sqrt(2)] - ω = 2.0 * pi * sqrt(2.0) - - H_lower = 2.0 + 0.1 * sin(ω * x[1] + t) - H_upper = 4.0 + 0.1 * cos(ω * x[1] + t) - v_lower = 1.0 - v_upper = 0.9 - b = 1.0 + 0.1 * cos(2.0 * ω * x[1]) - - return prim2cons(SVector(H_upper, v_upper, H_lower, v_lower, b), equations) -end - -""" - source_terms_convergence_test(u, x, t, equations::ShallowWaterTwoLayerEquations1D) - -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). -""" -@inline function source_terms_convergence_test(u, x, t, - equations::ShallowWaterTwoLayerEquations1D) - # Same settings as in `initial_condition_convergence_test`. Some derivative simplify because - # this manufactured solution velocity is taken to be constant - ω = 2 * pi * sqrt(2.0) - - du1 = (-0.1 * cos(t + ω * x[1]) - 0.1 * sin(t + ω * x[1]) - - 0.09 * ω * cos(t + ω * x[1]) + - -0.09 * ω * sin(t + ω * x[1])) - du2 = (5.0 * (-0.1 * ω * cos(t + ω * x[1]) - 0.1 * ω * sin(t + ω * x[1])) * - (4.0 + 0.2 * cos(t + ω * x[1]) + - -0.2 * sin(t + ω * x[1])) + - 0.1 * ω * (20.0 + cos(t + ω * x[1]) - sin(t + ω * x[1])) * - cos(t + - ω * x[1]) - 0.09 * cos(t + ω * x[1]) - 0.09 * sin(t + ω * x[1]) - - 0.081 * ω * cos(t + ω * x[1]) + - -0.081 * ω * sin(t + ω * x[1])) - du3 = 0.1 * cos(t + ω * x[1]) + 0.1 * ω * cos(t + ω * x[1]) + - 0.2 * ω * sin(2.0 * ω * x[1]) - du4 = ((10.0 + sin(t + ω * x[1]) - cos(2ω * x[1])) * - (-0.09 * ω * cos(t + ω * x[1]) - 0.09 * ω * sin(t + - ω * x[1]) - - 0.2 * ω * sin(2 * ω * x[1])) + 0.1 * cos(t + ω * x[1]) + - 0.1 * ω * cos(t + ω * x[1]) + - 5.0 * (0.1 * ω * cos(t + ω * x[1]) + 0.2 * ω * sin(2.0 * ω * x[1])) * - (2.0 + 0.2 * sin(t + ω * x[1]) + - -0.2 * cos(2.0 * ω * x[1])) + 0.2 * ω * sin(2.0 * ω * x[1])) - - return SVector(du1, du2, du3, du4, zero(eltype(u))) -end - -""" - boundary_condition_slip_wall(u_inner, orientation_or_normal, x, t, surface_flux_function, - equations::ShallowWaterTwoLayerEquations1D) - -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 boundary_condition_slip_wall(u_inner, orientation_or_normal, direction, - x, t, surface_flux_function, - equations::ShallowWaterTwoLayerEquations1D) - # create the "external" boundary solution state - u_boundary = SVector(u_inner[1], - -u_inner[2], - u_inner[3], - -u_inner[4], - u_inner[5]) - - # calculate the boundary flux - if iseven(direction) # u_inner is "left" of boundary, u_boundary is "right" of boundary - f = surface_flux_function(u_inner, u_boundary, orientation_or_normal, equations) - else # u_boundary is "left" of boundary, u_inner is "right" of boundary - f = surface_flux_function(u_boundary, u_inner, orientation_or_normal, equations) - end - return f -end - -# Calculate 1D flux for a single point -# Note, the bottom topography has no flux -@inline function flux(u, orientation::Integer, - equations::ShallowWaterTwoLayerEquations1D) - h_upper, h_v_upper, h_lower, h_v_lower, _ = u - - # Calculate velocities - v_upper, v_lower = velocity(u, equations) - # Calculate pressure - p_upper = 0.5 * equations.gravity * h_upper^2 - p_lower = 0.5 * equations.gravity * h_lower^2 - - f1 = h_v_upper - f2 = h_v_upper * v_upper + p_upper - f3 = h_v_lower - f4 = h_v_lower * v_lower + p_lower - - return SVector(f1, f2, f3, f4, zero(eltype(u))) -end - -""" - flux_nonconservative_ersing_etal(u_ll, u_rr, orientation::Integer, - equations::ShallowWaterTwoLayerEquations1D) - -!!! 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 [`ShallowWaterTwoLayerEquations1D`](@ref) and an -additional term that couples the momentum of both layers. - -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 flux_nonconservative_ersing_etal(u_ll, u_rr, - orientation::Integer, - equations::ShallowWaterTwoLayerEquations1D) - # Pull the necessary left and right state information - h_upper_ll, h_lower_ll = waterheight(u_ll, equations) - h_upper_rr, h_lower_rr = waterheight(u_rr, equations) - b_rr = u_rr[5] - b_ll = u_ll[5] - - # Calculate jumps - h_upper_jump = (h_upper_rr - h_upper_ll) - h_lower_jump = (h_lower_rr - h_lower_ll) - b_jump = (b_rr - b_ll) - - z = zero(eltype(u_ll)) - - # Bottom gradient nonconservative term: (0, g*h_upper*(b+h_lower)_x, - # 0, g*h_lower*(b+r*h_upper)_x, 0) - f = SVector(z, - equations.gravity * h_upper_ll * (b_jump + h_lower_jump), - z, - equations.gravity * h_lower_ll * (b_jump + equations.r * h_upper_jump), - z) - return f -end - -""" - flux_wintermeyer_etal(u_ll, u_rr, orientation, - equations::ShallowWaterTwoLayerEquations1D) - -Total energy conservative (mathematical entropy for two-layer shallow water equations) split form. -When the bottom topography is nonzero this scheme will be well-balanced when used with the -nonconservative [`flux_nonconservative_ersing_etal`](@ref). To obtain the flux for the -two-layer shallow water equations the flux that is described in the paper for the normal shallow -water equations is used within each layer. - -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 flux_wintermeyer_etal(u_ll, u_rr, - orientation::Integer, - equations::ShallowWaterTwoLayerEquations1D) - # Unpack left and right state - h_upper_ll, h_v_upper_ll, h_lower_ll, h_v_lower_ll, _ = u_ll - h_upper_rr, h_v_upper_rr, h_lower_rr, h_v_lower_rr, _ = u_rr - - # Get the velocities on either side - v_upper_ll, v_lower_ll = velocity(u_ll, equations) - v_upper_rr, v_lower_rr = velocity(u_rr, equations) - - # Average each factor of products in flux - v_upper_avg = 0.5 * (v_upper_ll + v_upper_rr) - v_lower_avg = 0.5 * (v_lower_ll + v_lower_rr) - p_upper_avg = 0.5 * equations.gravity * h_upper_ll * h_upper_rr - p_lower_avg = 0.5 * equations.gravity * h_lower_ll * h_lower_rr - - # Calculate fluxes - f1 = 0.5 * (h_v_upper_ll + h_v_upper_rr) - f2 = f1 * v_upper_avg + p_upper_avg - f3 = 0.5 * (h_v_lower_ll + h_v_lower_rr) - f4 = f3 * v_lower_avg + p_lower_avg - - return SVector(f1, f2, f3, f4, zero(eltype(u_ll))) -end - -""" - flux_es_ersing_etal(u_ll, u_rr, orientation_or_normal_direction, - equations::ShallowWaterTwoLayerEquations1D) -Entropy stable surface flux for the two-layer shallow water equations. Uses the entropy conservative -[`flux_wintermeyer_etal`](@ref) and adds a Lax-Friedrichs type dissipation dependent on the jump of -entropy variables. - -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 flux_es_ersing_etal(u_ll, u_rr, - orientation::Integer, - equations::ShallowWaterTwoLayerEquations1D) - # Compute entropy conservative flux but without the bottom topography - f_ec = flux_wintermeyer_etal(u_ll, u_rr, - orientation, - equations) - - # Get maximum signal velocity - λ = max_abs_speed_naive(u_ll, u_rr, orientation, equations) - # Get entropy variables but without the bottom topography - q_rr = cons2entropy(u_rr, equations) - q_ll = cons2entropy(u_ll, equations) - - # Average values from left and right - u_avg = (u_ll + u_rr) / 2 - - # Introduce variables for better readability - rho_upper = equations.rho_upper - rho_lower = equations.rho_lower - g = equations.gravity - drho = rho_upper - rho_lower - - # Compute entropy Jacobian coefficients - h11 = -rho_lower / (g * rho_upper * drho) - h12 = -rho_lower * u_avg[2] / (g * rho_upper * u_avg[1] * drho) - h13 = 1.0 / (g * drho) - h14 = u_avg[4] / (g * u_avg[3] * drho) - h21 = -rho_lower * u_avg[2] / (g * rho_upper * u_avg[1] * drho) - h22 = ((g * rho_upper * u_avg[1]^3 - g * rho_lower * u_avg[1]^3 + - -rho_lower * u_avg[2]^2) / (g * rho_upper * u_avg[1]^2 * drho)) - h23 = u_avg[2] / (g * u_avg[1] * drho) - h24 = u_avg[2] * u_avg[4] / (g * u_avg[1] * u_avg[3] * drho) - h31 = 1.0 / (g * drho) - h32 = u_avg[2] / (g * u_avg[1] * drho) - h33 = -1.0 / (g * drho) - h34 = -u_avg[4] / (g * u_avg[3] * drho) - h41 = u_avg[4] / (g * u_avg[3] * drho) - h42 = u_avg[2] * u_avg[4] / (g * u_avg[1] * u_avg[3] * drho) - h43 = -u_avg[4] / (g * u_avg[3] * drho) - h44 = ((g * rho_upper * u_avg[3]^3 - g * rho_lower * u_avg[3]^3 + - -rho_lower * u_avg[4]^2) / (g * rho_lower * u_avg[3]^2 * drho)) - - # Entropy Jacobian matrix - H = @SMatrix [[h11;; h12;; h13;; h14;; 0]; - [h21;; h22;; h23;; h24;; 0]; - [h31;; h32;; h33;; h34;; 0]; - [h41;; h42;; h43;; h44;; 0]; - [0;; 0;; 0;; 0;; 0]] - - # Add dissipation to entropy conservative flux to obtain entropy stable flux - f_es = f_ec - 0.5 * λ * H * (q_rr - q_ll) - - return SVector(f_es[1], f_es[2], f_es[3], f_es[4], zero(eltype(u_ll))) -end - -# Calculate approximation for maximum wave speed for local Lax-Friedrichs-type dissipation as the -# maximum velocity magnitude plus the maximum speed of sound. This function uses approximate -# eigenvalues using the speed of the barotropic mode as there is no simple way to calculate them -# analytically. -# -# A good overview of the derivation is given in: -# - Jonas Nycander, Andrew McC. Hogg, Leela M. Frankcombe (2008) -# Open boundary conditions for nonlinear channel Flows -# [DOI: 10.1016/j.ocemod.2008.06.003](https://doi.org/10.1016/j.ocemod.2008.06.003) -@inline function max_abs_speed_naive(u_ll, u_rr, - orientation::Integer, - equations::ShallowWaterTwoLayerEquations1D) - # Unpack left and right state - h_upper_ll, h_v_upper_ll, h_lower_ll, h_v_lower_ll, _ = u_ll - h_upper_rr, h_v_upper_rr, h_lower_rr, h_v_lower_rr, _ = u_rr - - # Get the averaged velocity - v_m_ll = (h_v_upper_ll + h_v_lower_ll) / (h_upper_ll + h_lower_ll) - v_m_rr = (h_v_upper_rr + h_v_lower_rr) / (h_upper_rr + h_lower_rr) - - # Calculate the wave celerity on the left and right - h_upper_ll, h_lower_ll = waterheight(u_ll, equations) - h_upper_rr, h_lower_rr = waterheight(u_rr, equations) - c_ll = sqrt(equations.gravity * (h_upper_ll + h_lower_ll)) - c_rr = sqrt(equations.gravity * (h_upper_rr + h_lower_rr)) - - return (max(abs(v_m_ll) + c_ll, abs(v_m_rr) + c_rr)) -end - -# Specialized `DissipationLocalLaxFriedrichs` to avoid spurious dissipation in the bottom -# topography -@inline function (dissipation::DissipationLocalLaxFriedrichs)(u_ll, u_rr, - orientation_or_normal_direction, - equations::ShallowWaterTwoLayerEquations1D) - λ = dissipation.max_abs_speed(u_ll, u_rr, orientation_or_normal_direction, - equations) - diss = -0.5 * λ * (u_rr - u_ll) - return SVector(diss[1], diss[2], diss[3], diss[4], zero(eltype(u_ll))) -end - -# Absolute speed of the barotropic mode -@inline function max_abs_speeds(u, equations::ShallowWaterTwoLayerEquations1D) - h_upper, h_v_upper, h_lower, h_v_lower, _ = u - - # Calculate averaged velocity of both layers - v_m = (h_v_upper + h_v_lower) / (h_upper + h_lower) - c = sqrt(equations.gravity * (h_upper + h_lower)) - - return (abs(v_m) + c) -end - -# Helper function to extract the velocity vector from the conservative variables -@inline function velocity(u, equations::ShallowWaterTwoLayerEquations1D) - h_upper, h_v_upper, h_lower, h_v_lower, _ = u - - v_upper = h_v_upper / h_upper - v_lower = h_v_lower / h_lower - return SVector(v_upper, v_lower) -end - -# Convert conservative variables to primitive -@inline function cons2prim(u, equations::ShallowWaterTwoLayerEquations1D) - h_upper, _, h_lower, _, b = u - - H_lower = h_lower + b - H_upper = h_lower + h_upper + b - v_upper, v_lower = velocity(u, equations) - return SVector(H_upper, v_upper, H_lower, v_lower, b) -end - -# Convert conservative variables to entropy variables -# Note, only the first four are the entropy variables, the fifth entry still just carries the -# bottom topography values for convenience -@inline function cons2entropy(u, equations::ShallowWaterTwoLayerEquations1D) - h_upper, _, h_lower, _, b = u - v_upper, v_lower = velocity(u, equations) - - w1 = (equations.rho_upper * - (equations.gravity * (h_upper + h_lower + b) - 0.5 * v_upper^2)) - w2 = equations.rho_upper * v_upper - w3 = (equations.rho_lower * - (equations.gravity * (equations.r * h_upper + h_lower + b) - 0.5 * v_lower^2)) - w4 = equations.rho_lower * v_lower - return SVector(w1, w2, w3, w4, b) -end - -# Convert primitive to conservative variables -@inline function prim2cons(prim, equations::ShallowWaterTwoLayerEquations1D) - H_upper, v_upper, H_lower, v_lower, b = prim - - h_lower = H_lower - b - h_upper = H_upper - h_lower - b - h_v_upper = h_upper * v_upper - h_v_lower = h_lower * v_lower - return SVector(h_upper, h_v_upper, h_lower, h_v_lower, b) -end - -@inline function waterheight(u, equations::ShallowWaterTwoLayerEquations1D) - return SVector(u[1], u[3]) -end - -# Entropy function for the shallow water equations is the total energy -@inline function entropy(cons, equations::ShallowWaterTwoLayerEquations1D) - energy_total(cons, equations) -end - -# Calculate total energy for a conservative state `cons` -@inline function energy_total(cons, equations::ShallowWaterTwoLayerEquations1D) - h_upper, h_v_upper, h_lower, h_v_lower, b = cons - # Set new variables for better readability - g = equations.gravity - rho_upper = equations.rho_upper - rho_lower = equations.rho_lower - - e = (0.5 * rho_upper * (h_v_upper^2 / h_upper + g * h_upper^2) + - 0.5 * rho_lower * (h_v_lower^2 / h_lower + g * h_lower^2) + - g * rho_lower * h_lower * b + g * rho_upper * h_upper * (h_lower + b)) - return e -end - -# Calculate kinetic energy for a conservative state `cons` -@inline function energy_kinetic(u, equations::ShallowWaterTwoLayerEquations1D) - h_upper, h_v_upper, h_lower, h_v_lower, _ = u - return (0.5 * equations.rho_upper * h_v_upper^2 / h_upper + - 0.5 * equations.rho_lower * h_v_lower^2 / h_lower) -end - -# Calculate potential energy for a conservative state `cons` -@inline function energy_internal(cons, equations::ShallowWaterTwoLayerEquations1D) - return energy_total(cons, equations) - energy_kinetic(cons, equations) -end - -# Calculate the error for the "lake-at-rest" test case where H = h_upper+h_lower+b should -# be a constant value over time -@inline function lake_at_rest_error(u, equations::ShallowWaterTwoLayerEquations1D) - h_upper, _, h_lower, _, b = u - return abs(equations.H0 - (h_upper + h_lower + b)) -end -end # @muladd diff --git a/src/equations/shallow_water_two_layer_2d.jl b/src/equations/shallow_water_two_layer_2d.jl deleted file mode 100644 index a31d881f2ef..00000000000 --- a/src/equations/shallow_water_two_layer_2d.jl +++ /dev/null @@ -1,805 +0,0 @@ -# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). -# Since these FMAs can increase the performance of many numerical algorithms, -# we need to opt-in explicitly. -# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. -@muladd begin -#! format: noindent - -# TODO: TrixiShallowWater: 2D two layer equations should move to new package - -@doc raw""" - ShallowWaterTwoLayerEquations2D(gravity, H0, rho_upper, rho_lower) - -Two-Layer Shallow water equations (2LSWE) in two space dimension. The equations are given by -```math -\begin{alignat*}{8} -&\frac{\partial}{\partial t}h_{upper} -&&+ \frac{\partial}{\partial x}\left(h_{upper} v_{1,upper}\right) -&&+ \frac{\partial}{\partial y}\left(h_{upper} v_{2,upper}\right) \quad -&&= \quad 0 \\ -&\frac{\partial}{\partial t}\left(h_{upper} v_{1,upper}\right) -&&+ \frac{\partial}{\partial x}\left(h_{upper} v_{1,upper}^2 + \frac{gh_{upper}^2}{2}\right) -&&+ \frac{\partial}{\partial y}\left(h_{upper} v_{1,upper} v_{2,upper}\right) \quad -&&= -gh_{upper}\frac{\partial}{\partial x}\left(b+h_{lower}\right) \\ -&\frac{\partial}{\partial t}\left(h_{upper} v_{2,upper}\right) -&&+ \frac{\partial}{\partial x}\left(h_{upper} v_{1,upper} v_{2,upper}\right) -&&+ \frac{\partial}{\partial y}\left(h_{upper} v_{2,upper}^2 + \frac{gh_{upper}^2}{2}\right) -&&= -gh_{upper}\frac{\partial}{\partial y}\left(b+h_{lower}\right)\\ -&\frac{\partial}{\partial t}h_{lower} -&&+ \frac{\partial}{\partial x}\left(h_{lower} v_{1,lower}\right) -&&+ \frac{\partial}{\partial y}\left(h_{lower} v_{2,lower}\right) -&&= \quad 0 \\ -&\frac{\partial}{\partial t}\left(h_{lower} v_{1,lower}\right) -&&+ \frac{\partial}{\partial x}\left(h_{lower} v_{1,lower}^2 + \frac{gh_{lower}^2}{2}\right) -&&+ \frac{\partial}{\partial y}\left(h_{lower} v_{1,lower} v_{2,lower}\right) -&&= -gh_{lower}\frac{\partial}{\partial x}\left(b+\frac{\rho_{upper}}{\rho_{lower}} h_{upper}\right)\\ -&\frac{\partial}{\partial t}\left(h_{lower} v_{2,lower}\right) -&&+ \frac{\partial}{\partial x}\left(h_{lower} v_{1,lower} v_{2,lower}\right) -&&+ \frac{\partial}{\partial y}\left(h_{lower} v_{2,lower}^2 + \frac{gh_{lower}^2}{2}\right) -&&= -gh_{lower}\frac{\partial}{\partial y}\left(b+\frac{\rho_{upper}}{\rho_{lower}} h_{upper}\right) -\end{alignat*} -``` -The unknown quantities of the 2LSWE are the water heights of the lower layer ``h_{lower}`` and the -upper -layer ``h_{upper}`` and the respective velocities in x-direction ``v_{1,lower}`` and ``v_{1,upper}`` and in y-direction -``v_{2,lower}`` and ``v_{2,upper}``. The gravitational constant is denoted by `g`, the layer densitites by -``\rho_{upper}``and ``\rho_{lower}`` and the (possibly) variable bottom topography function by ``b(x)``. -Conservative variable water height ``h_{lower}`` is measured from the bottom topography ``b`` and ``h_{upper}`` -relative to ``h_{lower}``, therefore one also defines the total water heights as ``H_{lower} = h_{lower} + b`` and -``H_{upper} = h_{upper} + h_{lower} + b``. - -The densities must be chosen such that ``\rho_{upper} < \rho_{lower}``, to make sure that the heavier fluid -``\rho_{lower}`` is in the bottom layer and the lighter fluid ``\rho_{upper}`` in the upper layer. - -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. - -The bottom topography function ``b(x)`` is set inside the initial condition routine -for a particular problem setup. - -In addition to the unknowns, Trixi 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's visualization tools will visualize the bottom topography by default. - -A good introduction for the 2LSWE is available in Chapter 12 of the book: - - Benoit Cushman-Roisin (2011)\ - Introduction to geophyiscal fluid dynamics: physical and numerical aspects\ - \ - ISBN: 978-0-12-088759-0 -""" -struct ShallowWaterTwoLayerEquations2D{RealT <: Real} <: - AbstractShallowWaterEquations{2, 7} - gravity::RealT # gravitational constant - H0::RealT # constant "lake-at-rest" total water height - rho_upper::RealT # lower layer density - rho_lower::RealT # upper layer density - r::RealT # ratio of rho_upper / rho_lower -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. Densities must be specified such that rho_upper < rho_lower. -function ShallowWaterTwoLayerEquations2D(; gravity_constant, - H0 = zero(gravity_constant), rho_upper, - rho_lower) - # Assign density ratio if rho_upper <= rho_lower - if rho_upper > rho_lower - error("Invalid input: Densities must be chosen such that rho_upper <= rho_lower") - else - r = rho_upper / rho_lower - end - ShallowWaterTwoLayerEquations2D(gravity_constant, H0, rho_upper, rho_lower, r) -end - -have_nonconservative_terms(::ShallowWaterTwoLayerEquations2D) = True() -function varnames(::typeof(cons2cons), ::ShallowWaterTwoLayerEquations2D) - ("h_upper", "h_v1_upper", "h_v2_upper", "h_lower", "h_v1_lower", "h_v2_lower", "b") -end -# Note, we use the total water height, H_upper = h_upper + h_lower + b, and first layer total height -# H_lower = h_lower + b as the first primitive variable for easier visualization and setting initial -# conditions -function varnames(::typeof(cons2prim), ::ShallowWaterTwoLayerEquations2D) - ("H_upper", "v1_upper", "v2_upper", "H_lower", "v1_lower", "v2_lower", "b") -end - -# Set initial conditions at physical location `x` for time `t` -""" - initial_condition_convergence_test(x, t, equations::ShallowWaterTwoLayerEquations2D) - -A smooth initial condition used for convergence tests in combination with -[`source_terms_convergence_test`](@ref). Constants must be set to ``rho_{upper} = 0.9``, -``rho_{lower} = 1.0``, ``g = 10.0``. -""" -function initial_condition_convergence_test(x, t, - equations::ShallowWaterTwoLayerEquations2D) - # some constants are chosen such that the function is periodic on the domain [0,sqrt(2)]^2] - ω = 2.0 * pi * sqrt(2.0) - - H_lower = 2.0 + 0.1 * sin(ω * x[1] + t) * cos(ω * x[2] + t) - H_upper = 4.0 + 0.1 * cos(ω * x[1] + t) * sin(ω * x[2] + t) - v1_lower = 1.0 - v1_upper = 0.9 - v2_lower = 0.9 - v2_upper = 1.0 - b = 1.0 + 0.1 * cos(0.5 * ω * x[1]) * sin(0.5 * ω * x[2]) - - return prim2cons(SVector(H_upper, v1_upper, v2_upper, H_lower, v1_lower, v2_lower, - b), equations) -end - -""" - source_terms_convergence_test(u, x, t, equations::ShallowWaterTwoLayerEquations2D) - -Source terms used for convergence tests in combination with -[`initial_condition_convergence_test`](@ref). -""" -@inline function source_terms_convergence_test(u, x, t, - equations::ShallowWaterTwoLayerEquations2D) - # Same settings as in `initial_condition_convergence_test`. - # some constants are chosen such that the function is periodic on the domain [0,sqrt(2)]^2] - ω = 2.0 * pi * sqrt(2.0) - - # Source terms obtained with SymPy - du1 = 0.01 * ω * cos(t + ω * x[1]) * cos(t + ω * x[2]) + - 0.01 * ω * sin(t + ω * x[1]) * sin(t + ω * x[2]) - du2 = (5.0 * - (-0.1 * ω * cos(t + ω * x[1]) * cos(t + ω * x[2]) - - 0.1 * ω * sin(t + ω * x[1]) * sin(t + - ω * x[2])) * - (4.0 + 0.2cos(t + ω * x[1]) * sin(t + ω * x[2]) - - 0.2 * sin(t + ω * x[1]) * cos(t + - ω * x[2])) + - 0.009 * ω * cos(t + ω * x[1]) * cos(t + ω * x[2]) + - 0.009 * ω * sin(t + ω * x[1]) * sin(t + - ω * x[2]) + - 0.1 * ω * - (20.0 + cos(t + ω * x[1]) * sin(t + ω * x[2]) - - sin(t + ω * x[1]) * cos(t + - ω * x[2])) * cos(t + ω * x[1]) * cos(t + ω * x[2])) - du3 = (5.0 * - (0.1 * ω * cos(t + ω * x[1]) * cos(t + ω * x[2]) + - 0.1 * ω * sin(t + ω * x[1]) * sin(t + - ω * x[2])) * - (4.0 + 0.2 * cos(t + ω * x[1]) * sin(t + ω * x[2]) - - 0.2 * sin(t + ω * x[1]) * cos(t + - ω * x[2])) + - 0.01ω * cos(t + ω * x[1]) * cos(t + ω * x[2]) + - 0.01 * ω * sin(t + ω * x[1]) * sin(t + ω * x[2]) + - -0.1 * ω * - (20.0 + cos(t + ω * x[1]) * sin(t + ω * x[2]) - - sin(t + ω * x[1]) * cos(t + ω * x[2])) * sin(t + - ω * x[1]) * sin(t + ω * x[2])) - du4 = (0.1 * cos(t + ω * x[1]) * cos(t + ω * x[2]) + - 0.1 * ω * cos(t + ω * x[1]) * cos(t + ω * x[2]) + - 0.05 * ω * sin(0.5 * ω * x[1]) * sin(0.5 * ω * x[2]) - - 0.1 * sin(t + ω * x[1]) * sin(t + ω * x[2]) + - -0.045 * ω * cos(0.5 * ω * x[1]) * cos(0.5 * ω * x[2]) - - 0.09 * ω * sin(t + ω * x[1]) * sin(t + ω * x[2])) - du5 = ((10.0 + sin(t + ω * x[1]) * cos(t + ω * x[2]) - - cos(0.5 * ω * x[1]) * sin(0.5 * ω * x[2])) * (-0.09 * ω * cos(t + - ω * x[1]) * cos(t + ω * x[2]) - - 0.09 * ω * sin(t + ω * x[1]) * sin(t + ω * x[2]) + - -0.05 * ω * sin(0.5 * ω * x[1]) * sin(0.5 * ω * x[2])) + - 5.0 * - (0.1 * ω * cos(t + ω * x[1]) * cos(t + ω * x[2]) + - 0.05 * ω * sin(0.5 * ω * x[1]) * sin(0.5 * ω * x[2])) * - (2.0 + 0.2 * sin(t + ω * x[1]) * cos(t + ω * x[2]) + - -0.2 * cos(0.5 * ω * x[1]) * sin(0.5 * ω * x[2])) + - 0.1 * cos(t + ω * x[1]) * cos(t + ω * x[2]) + - 0.1 * ω * cos(t + - ω * x[1]) * cos(t + ω * x[2]) + - 0.05 * ω * sin(0.5 * ω * x[1]) * sin(0.5 * ω * x[2]) - - 0.1 * sin(t + - ω * x[1]) * sin(t + ω * x[2]) - - 0.045 * ω * cos(0.5 * ω * x[1]) * cos(0.5 * ω * x[2]) - - 0.09 * ω * sin(t + - ω * x[1]) * sin(t + ω * x[2])) - du6 = ((10.0 + sin(t + ω * x[1]) * cos(t + ω * x[2]) + - -cos(0.5 * ω * x[1]) * sin(0.5 * ω * x[2])) * - (0.05 * ω * cos(0.5 * ω * x[1]) * cos(0.5 * ω * x[2]) + - 0.09 * ω * cos(t + ω * x[1]) * cos(t + ω * x[2]) + - 0.09 * ω * sin(t + ω * x[1]) * sin(t + ω * x[2])) + - 5.0 * - (-0.05 * ω * cos(0.5 * ω * x[1]) * cos(0.5 * ω * x[2]) - - 0.1 * ω * sin(t + ω * x[1]) * sin(t + - ω * x[2])) * - (2.0 + 0.2 * sin(t + ω * x[1]) * cos(t + ω * x[2]) + - -0.2 * cos(0.5 * ω * x[1]) * sin(0.5 * ω * x[2])) + - 0.09cos(t + ω * x[1]) * cos(t + ω * x[2]) + - 0.09 * ω * cos(t + ω * x[1]) * cos(t + ω * x[2]) + - 0.045 * ω * sin(0.5 * ω * x[1]) * sin(0.5 * ω * x[2]) + - -0.09 * sin(t + ω * x[1]) * sin(t + ω * x[2]) - - 0.0405 * ω * cos(0.5 * ω * x[1]) * cos(0.5 * ω * x[2]) + - -0.081 * ω * sin(t + ω * x[1]) * sin(t + ω * x[2])) - - return SVector(du1, du2, du3, du4, du5, du6, zero(eltype(u))) -end - -""" - boundary_condition_slip_wall(u_inner, normal_direction, x, t, surface_flux_function, - equations::ShallowWaterTwoLayerEquations2D) - -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 boundary_condition_slip_wall(u_inner, normal_direction::AbstractVector, - x, t, surface_flux_function, - equations::ShallowWaterTwoLayerEquations2D) - # normalize the outward pointing direction - normal = normal_direction / norm(normal_direction) - - # compute the normal velocity - v_normal_upper = normal[1] * u_inner[2] + normal[2] * u_inner[3] - v_normal_lower = normal[1] * u_inner[5] + normal[2] * u_inner[6] - - # create the "external" boundary solution state - u_boundary = SVector(u_inner[1], - u_inner[2] - 2.0 * v_normal_upper * normal[1], - u_inner[3] - 2.0 * v_normal_upper * normal[2], - u_inner[4], - u_inner[5] - 2.0 * v_normal_lower * normal[1], - u_inner[6] - 2.0 * v_normal_lower * normal[2], - u_inner[7]) - - # calculate the boundary flux - flux = surface_flux_function(u_inner, u_boundary, normal_direction, equations) - return flux -end - -# Calculate 1D flux for a single point -# Note, the bottom topography has no flux -@inline function flux(u, orientation::Integer, - equations::ShallowWaterTwoLayerEquations2D) - h_upper, h_v1_upper, h_v2_upper, h_lower, h_v1_lower, h_v2_lower, _ = u - - # Calculate velocities - v1_upper, v2_upper, v1_lower, v2_lower = velocity(u, equations) - - # Calculate pressure - p_upper = 0.5 * equations.gravity * h_upper^2 - p_lower = 0.5 * equations.gravity * h_lower^2 - - # Calculate fluxes depending on orientation - if orientation == 1 - f1 = h_v1_upper - f2 = h_v1_upper * v1_upper + p_upper - f3 = h_v1_upper * v2_upper - f4 = h_v1_lower - f5 = h_v1_lower * v1_lower + p_lower - f6 = h_v1_lower * v2_lower - else - f1 = h_v2_upper - f2 = h_v2_upper * v1_upper - f3 = h_v2_upper * v2_upper + p_upper - f4 = h_v2_lower - f5 = h_v2_lower * v1_lower - f6 = h_v2_lower * v2_lower + p_lower - end - return SVector(f1, f2, f3, f4, f5, f6, zero(eltype(u))) -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 flux(u, normal_direction::AbstractVector, - equations::ShallowWaterTwoLayerEquations2D) - h_upper, h_lower = waterheight(u, equations) - v1_upper, v2_upper, v1_lower, v2_lower = velocity(u, equations) - - v_normal_upper = v1_upper * normal_direction[1] + v2_upper * normal_direction[2] - v_normal_lower = v1_lower * normal_direction[1] + v2_lower * normal_direction[2] - h_v_upper_normal = h_upper * v_normal_upper - h_v_lower_normal = h_lower * v_normal_lower - - p_upper = 0.5 * equations.gravity * h_upper^2 - p_lower = 0.5 * equations.gravity * h_lower^2 - - f1 = h_v_upper_normal - f2 = h_v_upper_normal * v1_upper + p_upper * normal_direction[1] - f3 = h_v_upper_normal * v2_upper + p_upper * normal_direction[2] - f4 = h_v_lower_normal - f5 = h_v_lower_normal * v1_lower + p_lower * normal_direction[1] - f6 = h_v_lower_normal * v2_lower + p_lower * normal_direction[2] - - return SVector(f1, f2, f3, f4, f5, f6, zero(eltype(u))) -end - -""" - flux_nonconservative_ersing_etal(u_ll, u_rr, orientation::Integer, - equations::ShallowWaterTwoLayerEquations2D) - flux_nonconservative_ersing_etal(u_ll, u_rr, - normal_direction_ll::AbstractVector, - normal_direction_average::AbstractVector, - equations::ShallowWaterTwoLayerEquations2D) - -!!! 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 [`ShallowWaterTwoLayerEquations2D`](@ref) and an -additional term that couples the momentum of both layers. - -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 flux_nonconservative_ersing_etal(u_ll, u_rr, - orientation::Integer, - equations::ShallowWaterTwoLayerEquations2D) - # Pull the necessary left and right state information - h_upper_ll, h_lower_ll = waterheight(u_ll, equations) - h_upper_rr, h_lower_rr = waterheight(u_rr, equations) - b_rr = u_rr[7] - b_ll = u_ll[7] - - # Calculate jumps - h_upper_jump = (h_upper_rr - h_upper_ll) - h_lower_jump = (h_lower_rr - h_lower_ll) - b_jump = (b_rr - b_ll) - - z = zero(eltype(u_ll)) - - # Bottom gradient nonconservative term: (0, g*h_upper*(b + h_lower)_x, g*h_upper*(b + h_lower)_y , - # 0, g*h_lower*(b + r*h_upper)_x, - # g*h_lower*(b + r*h_upper)_y, 0) - if orientation == 1 - f = SVector(z, - equations.gravity * h_upper_ll * (b_jump + h_lower_jump), - z, z, - equations.gravity * h_lower_ll * - (b_jump + equations.r * h_upper_jump), - z, z) - else # orientation == 2 - f = SVector(z, z, - equations.gravity * h_upper_ll * (b_jump + h_lower_jump), - z, z, - equations.gravity * h_lower_ll * - (b_jump + equations.r * h_upper_jump), - z) - end - - return f -end - -@inline function flux_nonconservative_ersing_etal(u_ll, u_rr, - normal_direction_ll::AbstractVector, - normal_direction_average::AbstractVector, - equations::ShallowWaterTwoLayerEquations2D) - # Pull the necessary left and right state information - h_upper_ll, h_lower_ll = waterheight(u_ll, equations) - h_upper_rr, h_lower_rr = waterheight(u_rr, equations) - b_rr = u_rr[7] - b_ll = u_ll[7] - - # Calculate jumps - h_upper_jump = (h_upper_rr - h_upper_ll) - h_lower_jump = (h_lower_rr - h_lower_ll) - b_jump = (b_rr - b_ll) - - # Note this routine only uses the `normal_direction_average` and the average of the - # bottom topography to get a quadratic split form DG gradient on curved elements - return SVector(zero(eltype(u_ll)), - normal_direction_average[1] * equations.gravity * h_upper_ll * - (b_jump + h_lower_jump), - normal_direction_average[2] * equations.gravity * h_upper_ll * - (b_jump + h_lower_jump), - zero(eltype(u_ll)), - normal_direction_average[1] * equations.gravity * h_lower_ll * - (b_jump + equations.r * h_upper_jump), - normal_direction_average[2] * equations.gravity * h_lower_ll * - (b_jump + equations.r * h_upper_jump), - zero(eltype(u_ll))) -end - -""" - flux_wintermeyer_etal(u_ll, u_rr, orientation, - equations::ShallowWaterTwoLayerEquations2D) - flux_wintermeyer_etal(u_ll, u_rr, - normal_direction::AbstractVector, - equations::ShallowWaterTwoLayerEquations2D) - -Total energy conservative (mathematical entropy for two-layer shallow water equations) split form. -When the bottom topography is nonzero this scheme will be well-balanced when used with the -nonconservative [`flux_nonconservative_ersing_etal`](@ref). To obtain the flux for the -two-layer shallow water equations the flux that is described in the paper for the normal shallow -water equations is used within each layer. - -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 flux_wintermeyer_etal(u_ll, u_rr, - orientation::Integer, - equations::ShallowWaterTwoLayerEquations2D) - # Unpack left and right state - h_upper_ll, h_v1_upper_ll, h_v2_upper_ll, h_lower_ll, h_v1_lower_ll, h_v2_lower_ll, _ = u_ll - h_upper_rr, h_v1_upper_rr, h_v2_upper_rr, h_lower_rr, h_v1_lower_rr, h_v2_lower_rr, _ = u_rr - - # Get the velocities on either side - v1_upper_ll, v2_upper_ll, v1_lower_ll, v2_lower_ll = velocity(u_ll, equations) - v1_upper_rr, v2_upper_rr, v1_lower_rr, v2_lower_rr = velocity(u_rr, equations) - - # Average each factor of products in flux - v1_upper_avg = 0.5 * (v1_upper_ll + v1_upper_rr) - v1_lower_avg = 0.5 * (v1_lower_ll + v1_lower_rr) - v2_upper_avg = 0.5 * (v2_upper_ll + v2_upper_rr) - v2_lower_avg = 0.5 * (v2_lower_ll + v2_lower_rr) - p_upper_avg = 0.5 * equations.gravity * h_upper_ll * h_upper_rr - p_lower_avg = 0.5 * equations.gravity * h_lower_ll * h_lower_rr - - # Calculate fluxes depending on orientation - if orientation == 1 - f1 = 0.5 * (h_v1_upper_ll + h_v1_upper_rr) - f2 = f1 * v1_upper_avg + p_upper_avg - f3 = f1 * v2_upper_avg - f4 = 0.5 * (h_v1_lower_ll + h_v1_lower_rr) - f5 = f4 * v1_lower_avg + p_lower_avg - f6 = f4 * v2_lower_avg - else - f1 = 0.5 * (h_v2_upper_ll + h_v2_upper_rr) - f2 = f1 * v1_upper_avg - f3 = f1 * v2_upper_avg + p_upper_avg - f4 = 0.5 * (h_v2_lower_ll + h_v2_lower_rr) - f5 = f4 * v1_lower_avg - f6 = f4 * v2_lower_avg + p_lower_avg - end - - return SVector(f1, f2, f3, f4, f5, f6, zero(eltype(u_ll))) -end - -@inline function flux_wintermeyer_etal(u_ll, u_rr, - normal_direction::AbstractVector, - equations::ShallowWaterTwoLayerEquations2D) - # Unpack left and right state - h_upper_ll, h_v1_upper_ll, h_v2_upper_ll, h_lower_ll, h_v1_lower_ll, h_v2_lower_ll, _ = u_ll - h_upper_rr, h_v1_upper_rr, h_v2_upper_rr, h_lower_rr, h_v1_lower_rr, h_v2_lower_rr, _ = u_rr - - # Get the velocities on either side - v1_upper_ll, v2_upper_ll, v1_lower_ll, v2_lower_ll = velocity(u_ll, equations) - v1_upper_rr, v2_upper_rr, v1_lower_rr, v2_lower_rr = velocity(u_rr, equations) - - # Average each factor of products in flux - v1_upper_avg = 0.5 * (v1_upper_ll + v1_upper_rr) - v1_lower_avg = 0.5 * (v1_lower_ll + v1_lower_rr) - v2_upper_avg = 0.5 * (v2_upper_ll + v2_upper_rr) - v2_lower_avg = 0.5 * (v2_lower_ll + v2_lower_rr) - p_upper_avg = 0.5 * equations.gravity * h_upper_ll * h_upper_rr - p_lower_avg = 0.5 * equations.gravity * h_lower_ll * h_lower_rr - h_v1_upper_avg = 0.5 * (h_v1_upper_ll + h_v1_upper_rr) - h_v2_upper_avg = 0.5 * (h_v2_upper_ll + h_v2_upper_rr) - h_v1_lower_avg = 0.5 * (h_v1_lower_ll + h_v1_lower_rr) - h_v2_lower_avg = 0.5 * (h_v2_lower_ll + h_v2_lower_rr) - - # Calculate fluxes depending on normal_direction - f1 = h_v1_upper_avg * normal_direction[1] + h_v2_upper_avg * normal_direction[2] - f2 = f1 * v1_upper_avg + p_upper_avg * normal_direction[1] - f3 = f1 * v2_upper_avg + p_upper_avg * normal_direction[2] - f4 = h_v1_lower_avg * normal_direction[1] + h_v2_lower_avg * normal_direction[2] - f5 = f4 * v1_lower_avg + p_lower_avg * normal_direction[1] - f6 = f4 * v2_lower_avg + p_lower_avg * normal_direction[2] - - return SVector(f1, f2, f3, f4, f5, f6, zero(eltype(u_ll))) -end - -""" - flux_es_ersing_etal(u_ll, u_rr, orientation_or_normal_direction, - equations::ShallowWaterTwoLayerEquations2D) - -Entropy stable surface flux for the two-layer shallow water equations. Uses the entropy conservative -[`flux_wintermeyer_etal`](@ref) and adds a Lax-Friedrichs type dissipation dependent on the jump of -entropy variables. - -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 flux_es_ersing_etal(u_ll, u_rr, - orientation_or_normal_direction, - equations::ShallowWaterTwoLayerEquations2D) - # Compute entropy conservative flux but without the bottom topography - f_ec = flux_wintermeyer_etal(u_ll, u_rr, - orientation_or_normal_direction, - equations) - - # Get maximum signal velocity - λ = max_abs_speed_naive(u_ll, u_rr, orientation_or_normal_direction, equations) - - # Get entropy variables but without the bottom topography - q_rr = cons2entropy(u_rr, equations) - q_ll = cons2entropy(u_ll, equations) - - # Average values from left and right - u_avg = (u_ll + u_rr) / 2 - - # Introduce variables for better readability - rho_upper = equations.rho_upper - rho_lower = equations.rho_lower - g = equations.gravity - drho = rho_upper - rho_lower - - # Compute entropy Jacobian coefficients - h11 = -rho_lower / (g * rho_upper * drho) - h12 = -rho_lower * u_avg[2] / (g * rho_upper * u_avg[1] * drho) - h13 = -rho_lower * u_avg[3] / (g * rho_upper * u_avg[1] * drho) - h14 = 1.0 / (g * drho) - h15 = u_avg[5] / (g * u_avg[4] * drho) - h16 = u_avg[6] / (g * u_avg[4] * drho) - h21 = -rho_lower * u_avg[2] / (g * rho_upper * u_avg[1] * drho) - h22 = ((g * rho_upper * u_avg[1]^3 - g * rho_lower * u_avg[1]^3 + - -rho_lower * u_avg[2]^2) / (g * rho_upper * u_avg[1]^2 * drho)) - h23 = -rho_lower * u_avg[2] * u_avg[3] / (g * rho_upper * u_avg[1]^2 * drho) - h24 = u_avg[2] / (g * u_avg[1] * drho) - h25 = u_avg[2] * u_avg[5] / (g * u_avg[1] * u_avg[4] * drho) - h26 = u_avg[2] * u_avg[6] / (g * u_avg[1] * u_avg[4] * drho) - h31 = -rho_lower * u_avg[3] / (g * rho_upper * u_avg[1] * drho) - h32 = -rho_lower * u_avg[2] * u_avg[3] / (g * rho_upper * u_avg[1]^2 * drho) - h33 = ((g * rho_upper * u_avg[1]^3 - g * rho_lower * u_avg[1]^3 + - -rho_lower * u_avg[3]^2) / (g * rho_upper * u_avg[1]^2 * drho)) - h34 = u_avg[3] / (g * u_avg[1] * drho) - h35 = u_avg[3] * u_avg[5] / (g * u_avg[1] * u_avg[4] * drho) - h36 = u_avg[3] * u_avg[6] / (g * u_avg[1] * u_avg[4] * drho) - h41 = 1.0 / (g * drho) - h42 = u_avg[2] / (g * u_avg[1] * drho) - h43 = u_avg[3] / (g * u_avg[1] * drho) - h44 = -1.0 / (g * drho) - h45 = -u_avg[5] / (g * u_avg[4] * drho) - h46 = -u_avg[6] / (g * u_avg[4] * drho) - h51 = u_avg[5] / (g * u_avg[4] * drho) - h52 = u_avg[2] * u_avg[5] / (g * u_avg[1] * u_avg[4] * drho) - h53 = u_avg[3] * u_avg[5] / (g * u_avg[1] * u_avg[4] * drho) - h54 = -u_avg[5] / (g * u_avg[4] * drho) - h55 = ((g * rho_upper * u_avg[4]^3 - g * rho_lower * u_avg[4]^3 + - -rho_lower * u_avg[5]^2) / (g * rho_lower * u_avg[4]^2 * drho)) - h56 = -u_avg[5] * u_avg[6] / (g * u_avg[4]^2 * drho) - h61 = u_avg[6] / (g * u_avg[4] * drho) - h62 = u_avg[2] * u_avg[6] / (g * u_avg[1] * u_avg[4] * drho) - h63 = u_avg[3] * u_avg[6] / (g * u_avg[1] * u_avg[4] * drho) - h64 = -u_avg[6] / (g * u_avg[4] * drho) - h65 = -u_avg[5] * u_avg[6] / (g * u_avg[4]^2 * drho) - h66 = ((g * rho_upper * u_avg[4]^3 - g * rho_lower * u_avg[4]^3 + - -rho_lower * u_avg[6]^2) / (g * rho_lower * u_avg[4]^2 * drho)) - - # Entropy Jacobian matrix - H = @SMatrix [[h11;; h12;; h13;; h14;; h15;; h16;; 0]; - [h21;; h22;; h23;; h24;; h25;; h26;; 0]; - [h31;; h32;; h33;; h34;; h35;; h36;; 0]; - [h41;; h42;; h43;; h44;; h45;; h46;; 0]; - [h51;; h52;; h53;; h54;; h55;; h56;; 0]; - [h61;; h62;; h63;; h64;; h65;; h66;; 0]; - [0;; 0;; 0;; 0;; 0;; 0;; 0]] - - # Add dissipation to entropy conservative flux to obtain entropy stable flux - f_es = f_ec - 0.5 * λ * H * (q_rr - q_ll) - - return SVector(f_es[1], f_es[2], f_es[3], f_es[4], f_es[5], f_es[6], - zero(eltype(u_ll))) -end - -# Calculate approximation for maximum wave speed for local Lax-Friedrichs-type dissipation as the -# maximum velocity magnitude plus the maximum speed of sound. This function uses approximate -# eigenvalues using the speed of the barotropic mode as there is no simple way to calculate them -# analytically. -# -# A good overview of the derivation is given in: -# - Jonas Nycander, Andrew McC. Hogg, Leela M. Frankcombe (2008) -# Open boundary conditions for nonlinear channel Flows -# [DOI: 10.1016/j.ocemod.2008.06.003](https://doi.org/10.1016/j.ocemod.2008.06.003) -@inline function max_abs_speed_naive(u_ll, u_rr, - orientation::Integer, - equations::ShallowWaterTwoLayerEquations2D) - # Unpack left and right state - h_upper_ll, h_v1_upper_ll, h_v2_upper_ll, h_lower_ll, h_v1_lower_ll, h_v2_lower_ll, _ = u_ll - h_upper_rr, h_v1_upper_rr, h_v2_upper_rr, h_lower_rr, h_v1_lower_rr, h_v2_lower_rr, _ = u_rr - - # Calculate averaged velocity of both layers - if orientation == 1 - v_m_ll = (h_v1_upper_ll + h_v1_lower_ll) / (h_upper_ll + h_lower_ll) - v_m_rr = (h_v1_upper_rr + h_v1_lower_rr) / (h_upper_rr + h_lower_rr) - else - v_m_ll = (h_v2_upper_ll + h_v2_lower_ll) / (h_upper_ll + h_lower_ll) - v_m_rr = (h_v2_upper_rr + h_v2_lower_rr) / (h_upper_rr + h_lower_rr) - end - - # Calculate the wave celerity on the left and right - h_upper_ll, h_lower_ll = waterheight(u_ll, equations) - h_upper_rr, h_lower_rr = waterheight(u_rr, equations) - - c_ll = sqrt(equations.gravity * (h_upper_ll + h_lower_ll)) - c_rr = sqrt(equations.gravity * (h_upper_rr + h_lower_rr)) - - return (max(abs(v_m_ll), abs(v_m_rr)) + max(c_ll, c_rr)) -end - -@inline function max_abs_speed_naive(u_ll, u_rr, - normal_direction::AbstractVector, - equations::ShallowWaterTwoLayerEquations2D) - # Unpack left and right state - h_upper_ll, _, _, h_lower_ll, _, _, _ = u_ll - h_upper_rr, _, _, h_lower_rr, _, _, _ = u_rr - - # Extract and compute the velocities in the normal direction - v1_upper_ll, v2_upper_ll, v1_lower_ll, v2_lower_ll = velocity(u_ll, equations) - v1_upper_rr, v2_upper_rr, v1_lower_rr, v2_lower_rr = velocity(u_rr, equations) - - v_upper_dot_n_ll = v1_upper_ll * normal_direction[1] + - v2_upper_ll * normal_direction[2] - v_upper_dot_n_rr = v1_upper_rr * normal_direction[1] + - v2_upper_rr * normal_direction[2] - v_lower_dot_n_ll = v1_lower_ll * normal_direction[1] + - v2_lower_ll * normal_direction[2] - v_lower_dot_n_rr = v1_lower_rr * normal_direction[1] + - v2_lower_rr * normal_direction[2] - - # Calculate averaged velocity of both layers - v_m_ll = (v_upper_dot_n_ll * h_upper_ll + v_lower_dot_n_ll * h_lower_ll) / - (h_upper_ll + h_lower_ll) - v_m_rr = (v_upper_dot_n_rr * h_upper_rr + v_lower_dot_n_rr * h_lower_rr) / - (h_upper_rr + h_lower_rr) - - # Compute the wave celerity on the left and right - h_upper_ll, h_lower_ll = waterheight(u_ll, equations) - h_upper_rr, h_lower_rr = waterheight(u_rr, equations) - - c_ll = sqrt(equations.gravity * (h_upper_ll + h_lower_ll)) - c_rr = sqrt(equations.gravity * (h_upper_rr + h_lower_rr)) - - # The normal velocities are already scaled by the norm - return max(abs(v_m_ll), abs(v_m_rr)) + max(c_ll, c_rr) * norm(normal_direction) -end - -# Specialized `DissipationLocalLaxFriedrichs` to avoid spurious dissipation in the bottom topography -@inline function (dissipation::DissipationLocalLaxFriedrichs)(u_ll, u_rr, - orientation_or_normal_direction, - equations::ShallowWaterTwoLayerEquations2D) - λ = dissipation.max_abs_speed(u_ll, u_rr, orientation_or_normal_direction, - equations) - diss = -0.5 * λ * (u_rr - u_ll) - return SVector(diss[1], diss[2], diss[3], diss[4], diss[5], diss[6], - zero(eltype(u_ll))) -end - -# Absolute speed of the barotropic mode -@inline function max_abs_speeds(u, equations::ShallowWaterTwoLayerEquations2D) - h_upper, h_v1_upper, h_v2_upper, h_lower, h_v1_lower, h_v2_lower, _ = u - - # Calculate averaged velocity of both layers - v1_m = (h_v1_upper + h_v1_lower) / (h_upper + h_lower) - v2_m = (h_v2_upper + h_v2_lower) / (h_upper + h_lower) - - h_upper, h_lower = waterheight(u, equations) - v1_upper, v2_upper, v1_lower, v2_lower = velocity(u, equations) - - c = sqrt(equations.gravity * (h_upper + h_lower)) - return (max(abs(v1_m) + c, abs(v1_upper), abs(v1_lower)), - max(abs(v2_m) + c, abs(v2_upper), abs(v2_lower))) -end - -# Helper function to extract the velocity vector from the conservative variables -@inline function velocity(u, equations::ShallowWaterTwoLayerEquations2D) - h_upper, h_v1_upper, h_v2_upper, h_lower, h_v1_lower, h_v2_lower, _ = u - - v1_upper = h_v1_upper / h_upper - v2_upper = h_v2_upper / h_upper - v1_lower = h_v1_lower / h_lower - v2_lower = h_v2_lower / h_lower - - return SVector(v1_upper, v2_upper, v1_lower, v2_lower) -end - -# Convert conservative variables to primitive -@inline function cons2prim(u, equations::ShallowWaterTwoLayerEquations2D) - h_upper, _, _, h_lower, _, _, b = u - - H_lower = h_lower + b - H_upper = h_lower + h_upper + b - v1_upper, v2_upper, v1_lower, v2_lower = velocity(u, equations) - - return SVector(H_upper, v1_upper, v2_upper, H_lower, v1_lower, v2_lower, b) -end - -# Convert conservative variables to entropy variables -# Note, only the first four are the entropy variables, the fifth entry still just carries the bottom -# topography values for convenience. -# In contrast to general usage the entropy variables are denoted with q instead of w, because w is -# already used for velocity in y-Direction -@inline function cons2entropy(u, equations::ShallowWaterTwoLayerEquations2D) - h_upper, _, _, h_lower, _, _, b = u - # Assign new variables for better readability - rho_upper = equations.rho_upper - rho_lower = equations.rho_lower - v1_upper, v2_upper, v1_lower, v2_lower = velocity(u, equations) - - w1 = (rho_upper * (equations.gravity * (h_upper + h_lower + b) + - -0.5 * (v1_upper^2 + v2_upper^2))) - w2 = rho_upper * v1_upper - w3 = rho_upper * v2_upper - w4 = (rho_lower * (equations.gravity * (equations.r * h_upper + h_lower + b) + - -0.5 * (v1_lower^2 + v2_lower^2))) - w5 = rho_lower * v1_lower - w6 = rho_lower * v2_lower - return SVector(w1, w2, w3, w4, w5, w6, b) -end - -# Convert primitive to conservative variables -@inline function prim2cons(prim, equations::ShallowWaterTwoLayerEquations2D) - H_upper, v1_upper, v2_upper, H_lower, v1_lower, v2_lower, b = prim - - h_lower = H_lower - b - h_upper = H_upper - h_lower - b - h_v1_upper = h_upper * v1_upper - h_v2_upper = h_upper * v2_upper - h_v1_lower = h_lower * v1_lower - h_v2_lower = h_lower * v2_lower - return SVector(h_upper, h_v1_upper, h_v2_upper, h_lower, h_v1_lower, h_v2_lower, b) -end - -@inline function waterheight(u, equations::ShallowWaterTwoLayerEquations2D) - return SVector(u[1], u[4]) -end - -# Entropy function for the shallow water equations is the total energy -@inline function entropy(cons, equations::ShallowWaterTwoLayerEquations2D) - energy_total(cons, equations) -end - -# Calculate total energy for a conservative state `cons` -@inline function energy_total(cons, equations::ShallowWaterTwoLayerEquations2D) - h_upper, h_v1_upper, h_v2_upper, h_lower, h_v2_lower, h_v2_lower, b = cons - g = equations.gravity - rho_upper = equations.rho_upper - rho_lower = equations.rho_lower - - e = (0.5 * rho_upper * - (h_v1_upper^2 / h_upper + h_v2_upper^2 / h_upper + g * h_upper^2) + - 0.5 * rho_lower * - (h_v2_lower^2 / h_lower + h_v2_lower^2 / h_lower + g * h_lower^2) + - g * rho_lower * h_lower * b + g * rho_upper * h_upper * (h_lower + b)) - return e -end - -# Calculate kinetic energy for a conservative state `cons` -@inline function energy_kinetic(u, equations::ShallowWaterTwoLayerEquations2D) - h_upper, h_v1_upper, h_v2_upper, h_lower, h_v2_lower, h_v2_lower, _ = u - - return (0.5 * equations.rho_upper * h_v1_upper^2 / h_upper + - 0.5 * equations.rho_upper * h_v2_upper^2 / h_upper + - 0.5 * equations.rho_lower * h_v2_lower^2 / h_lower + - 0.5 * equations.rho_lower * h_v2_lower^2 / h_lower) -end - -# Calculate potential energy for a conservative state `cons` -@inline function energy_internal(cons, equations::ShallowWaterTwoLayerEquations2D) - return energy_total(cons, equations) - energy_kinetic(cons, equations) -end - -# Calculate the error for the "lake-at-rest" test case where H = h_upper+h_lower+b should -# be a constant value over time -@inline function lake_at_rest_error(u, equations::ShallowWaterTwoLayerEquations2D) - h_upper, _, _, h_lower, _, _, b = u - return abs(equations.H0 - (h_upper + h_lower + b)) -end -end # @muladd diff --git a/src/solvers/dgsem_tree/indicators.jl b/src/solvers/dgsem_tree/indicators.jl index bb9109f2762..9f25a6d2dbb 100644 --- a/src/solvers/dgsem_tree/indicators.jl +++ b/src/solvers/dgsem_tree/indicators.jl @@ -101,82 +101,6 @@ function Base.show(io::IO, ::MIME"text/plain", indicator::IndicatorHennemannGass summary_box(io, "IndicatorHennemannGassner", setup) end -# TODO: TrixiShallowWater: move the new indicator and all associated routines to the new package -""" - 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} <: - 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 cells -function IndicatorHennemannGassnerShallowWater(equations::AbstractShallowWaterEquations, - basis; - alpha_max = 0.5, - alpha_min = 0.001, - alpha_smooth = true, - variable) - alpha_max, alpha_min = promote(alpha_max, alpha_min) - cache = 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"), - ] - summary_box(io, "IndicatorHennemannGassnerShallowWater", setup) - end -end - function (indicator_hg::IndicatorHennemannGassner)(u, mesh, equations, dg::DGSEM, cache; kwargs...) @unpack alpha_smooth = indicator_hg diff --git a/src/solvers/dgsem_tree/indicators_1d.jl b/src/solvers/dgsem_tree/indicators_1d.jl index dff87bfe06c..4796ddcc602 100644 --- a/src/solvers/dgsem_tree/indicators_1d.jl +++ b/src/solvers/dgsem_tree/indicators_1d.jl @@ -24,115 +24,6 @@ function create_cache(typ::Type{IndicatorHennemannGassner}, mesh, create_cache(typ, equations, dg.basis) end -# Modified indicator for ShallowWaterEquations1D to apply full FV method on cells -# containing some "dry" LGL nodes. That is, if an element is partially "wet" then it becomes a -# full FV element. -# -# TODO: TrixiShallowWater: move new indicator type -function (indicator_hg::IndicatorHennemannGassnerShallowWater)(u::AbstractArray{<:Any, - 3}, - mesh, - equations::ShallowWaterEquations1D, - 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 cells 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 cells 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 cells 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 cells 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 - - @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 - 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 - apply_smoothing!(mesh, alpha, alpha_tmp, dg, cache) - end - - return alpha -end - # Use this function barrier and unpack inside to avoid passing closures to Polyester.jl # with @batch (@threaded). # Otherwise, @threaded does not work here with Julia ARM on macOS. diff --git a/src/solvers/dgsem_tree/indicators_2d.jl b/src/solvers/dgsem_tree/indicators_2d.jl index fa8ed481eb9..665d2254e5d 100644 --- a/src/solvers/dgsem_tree/indicators_2d.jl +++ b/src/solvers/dgsem_tree/indicators_2d.jl @@ -28,116 +28,6 @@ function create_cache(typ::Type{IndicatorHennemannGassner}, mesh, create_cache(typ, equations, dg.basis) end -# Modified indicator for ShallowWaterEquations2D to apply full FV method on cells -# containing some "dry" LGL nodes. That is, if an element is partially "wet" then it becomes a -# full FV element. -# -# TODO: TrixiShallowWater: move new indicator type -function (indicator_hg::IndicatorHennemannGassnerShallowWater)(u::AbstractArray{<:Any, - 4}, - mesh, - equations::ShallowWaterEquations2D, - 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 cells 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 cells 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 cells 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 cells 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 - - @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 - 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 - apply_smoothing!(mesh, alpha, alpha_tmp, dg, cache) - end - - return alpha -end - # Use this function barrier and unpack inside to avoid passing closures to Polyester.jl # with @batch (@threaded). # Otherwise, @threaded does not work here with Julia ARM on macOS. diff --git a/test/test_structured_2d.jl b/test/test_structured_2d.jl index 522510a42e3..f5fb169033a 100644 --- a/test/test_structured_2d.jl +++ b/test/test_structured_2d.jl @@ -1,7 +1,5 @@ module TestExamplesStructuredMesh2D -# TODO: TrixiShallowWater: move any wet/dry tests to new package - using Test using Trixi @@ -907,82 +905,6 @@ end 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 - @trixi_testset "elixir_mhd_ec_shockcapturing.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_ec_shockcapturing.jl"), l2=[0.0364192725149364, 0.0426667193422069, 0.04261673001449095, diff --git a/test/test_tree_1d.jl b/test/test_tree_1d.jl index 8b470278ffd..4a25a51a45e 100644 --- a/test/test_tree_1d.jl +++ b/test/test_tree_1d.jl @@ -42,8 +42,6 @@ isdir(outdir) && rm(outdir, recursive = true) # Shallow water include("test_tree_1d_shallowwater.jl") - # Two-layer Shallow Water - include("test_tree_1d_shallowwater_twolayer.jl") # FDSBP methods on the TreeMesh include("test_tree_1d_fdsbp.jl") diff --git a/test/test_tree_1d_shallowwater.jl b/test/test_tree_1d_shallowwater.jl index 2269e858928..f9be63b87fd 100644 --- a/test/test_tree_1d_shallowwater.jl +++ b/test/test_tree_1d_shallowwater.jl @@ -1,7 +1,5 @@ module TestExamples1DShallowWater -# TODO: TrixiShallowWater: move any wet/dry tests to new package - using Test using Trixi @@ -119,32 +117,6 @@ end 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=[ @@ -339,53 +311,6 @@ end 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 - @trixi_testset "elixir_shallow_water_quasi_1d_source_terms.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallow_water_quasi_1d_source_terms.jl"), diff --git a/test/test_tree_1d_shallowwater_twolayer.jl b/test/test_tree_1d_shallowwater_twolayer.jl deleted file mode 100644 index 180fb3ec3b3..00000000000 --- a/test/test_tree_1d_shallowwater_twolayer.jl +++ /dev/null @@ -1,74 +0,0 @@ -module TestExamples1DShallowWaterTwoLayer - -# TODO: TrixiShallowWater: move two layer tests to new package - -using Test -using Trixi - -include("test_trixi.jl") - -EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_1d_dgsem") - -@testset "Shallow Water Two layer" begin - @trixi_testset "elixir_shallowwater_twolayer_convergence.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, - "elixir_shallowwater_twolayer_convergence.jl"), - l2=[0.005012009872109003, 0.002091035326731071, - 0.005049271397924551, - 0.0024633066562966574, 0.0004744186597732739], - linf=[0.0213772149343594, 0.005385752427290447, - 0.02175023787351349, - 0.008212004668840978, 0.0008992474511784199], - 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=[8.949288784402005e-16, 4.0636427176237915e-17, - 0.001002881985401548, - 2.133351105037203e-16, 0.0010028819854016578], - linf=[2.6229018956769323e-15, 1.878451903240623e-16, - 0.005119880996670156, - 8.003199803957679e-16, 0.005119880996670666], - 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.1000774903431289, 0.5670692949571057, 0.08764242501014498, - 0.45412307886094555, 0.013638618139749523], - linf=[0.586718937495144, 2.1215606128311584, 0.5185911311186155, - 1.820382495072612, 0.5], - 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 - -end # module diff --git a/test/test_tree_2d_part3.jl b/test/test_tree_2d_part3.jl index ce9b3bc04f8..0eff564132c 100644 --- a/test/test_tree_2d_part3.jl +++ b/test/test_tree_2d_part3.jl @@ -26,9 +26,6 @@ isdir(outdir) && rm(outdir, recursive = true) # Shallow water include("test_tree_2d_shallowwater.jl") - # Two-Layer Shallow Water - include("test_tree_2d_shallowwater_twolayer.jl") - # FDSBP methods on the TreeMesh include("test_tree_2d_fdsbp.jl") end diff --git a/test/test_tree_2d_shallowwater.jl b/test/test_tree_2d_shallowwater.jl index 1f3dfbf5267..93a8cb63667 100644 --- a/test/test_tree_2d_shallowwater.jl +++ b/test/test_tree_2d_shallowwater.jl @@ -1,7 +1,5 @@ module TestExamples2DShallowWater -# TODO: TrixiShallowWater: move any wet/dry tests to new package - using Test using Trixi @@ -145,32 +143,6 @@ end 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=[ @@ -277,57 +249,6 @@ end 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=[ diff --git a/test/test_tree_2d_shallowwater_twolayer.jl b/test/test_tree_2d_shallowwater_twolayer.jl deleted file mode 100644 index 802bf4e021c..00000000000 --- a/test/test_tree_2d_shallowwater_twolayer.jl +++ /dev/null @@ -1,88 +0,0 @@ -module TestExamples2DShallowWaterTwoLayer - -# TODO: TrixiShallowWater: move two layer tests to new package - -using Test -using Trixi - -include("test_trixi.jl") - -EXAMPLES_DIR = joinpath(examples_dir(), "tree_2d_dgsem") - -@testset "Two-Layer Shallow Water" begin - @trixi_testset "elixir_shallowwater_twolayer_convergence.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, - "elixir_shallowwater_twolayer_convergence.jl"), - l2=[0.0004016779699408397, 0.005466339651545468, - 0.006148841330156112, - 0.0002882339012602492, 0.0030120142442780313, - 0.002680752838455618, - 8.873630921431545e-6], - linf=[0.002788654460984752, 0.01484602033450666, - 0.017572229756493973, - 0.0016010835493927011, 0.009369847995372549, - 0.008407961775489636, - 3.361991620143279e-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=[3.2935164267930016e-16, 4.6800825611195103e-17, - 4.843057532147818e-17, - 0.0030769233188015013, 1.4809161150389857e-16, - 1.509071695038043e-16, - 0.0030769233188014935], - linf=[2.248201624865942e-15, 2.346382070278936e-16, - 2.208565017494899e-16, - 0.026474051138910493, 9.237568031609006e-16, - 7.520758026187046e-16, - 0.026474051138910267], - 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 with flux_lax_friedrichs.jl" begin - @test_trixi_include(joinpath(EXAMPLES_DIR, - "elixir_shallowwater_twolayer_well_balanced.jl"), - l2=[2.0525741072929735e-16, 6.000589392730905e-17, - 6.102759428478984e-17, - 0.0030769233188014905, 1.8421386173122792e-16, - 1.8473184927121752e-16, - 0.0030769233188014935], - linf=[7.355227538141662e-16, 2.960836949170518e-16, - 4.2726562436938764e-16, - 0.02647405113891016, 1.038795478061861e-15, - 1.0401789378532516e-15, - 0.026474051138910267], - 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 - -end # module diff --git a/test/test_unit.jl b/test/test_unit.jl index c1379587cc8..1907a281718 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -420,11 +420,6 @@ end (1.0, 1.0), 1.0) @test_nowarn show(stdout, limiter_idp) - # TODO: TrixiShallowWater: move unit test - indicator_hg_swe = IndicatorHennemannGassnerShallowWater(1.0, 0.0, true, "variable", - "cache") - @test_nowarn show(stdout, indicator_hg_swe) - indicator_loehner = IndicatorLöhner(1.0, "variable", (; cache = nothing)) @test_nowarn show(stdout, indicator_loehner) @@ -543,7 +538,7 @@ end end @timed_testset "Shallow water conversion between conservative/entropy variables" begin - H, v1, v2, b = 3.5, 0.25, 0.1, 0.4 + H, v1, v2, b, a = 3.5, 0.25, 0.1, 0.4, 0.3 let equations = ShallowWaterEquations1D(gravity_constant = 9.8) cons_vars = prim2cons(SVector(H, v1, b), equations) @@ -572,6 +567,14 @@ end entropy_vars = cons2entropy(cons_vars, equations) @test cons_vars ≈ entropy2cons(entropy_vars, equations) end + + let equations = ShallowWaterEquationsQuasi1D(gravity_constant = 9.8) + cons_vars = prim2cons(SVector(H, v1, b, a), equations) + entropy_vars = cons2entropy(cons_vars, equations) + + total_energy = energy_total(cons_vars, equations) + @test entropy(cons_vars, equations) ≈ a * total_energy + end end @timed_testset "boundary_condition_do_nothing" begin @@ -697,6 +700,14 @@ end u = SVector(1, 0.5, 0.0) @test flux_hll(u, u, 1, equations) ≈ flux(u, 1, equations) + u_ll = SVector(0.1, 1.0, 0.0) + u_rr = SVector(0.1, 1.0, 0.0) + @test flux_hll(u_ll, u_rr, 1, equations) ≈ flux(u_ll, 1, equations) + + u_ll = SVector(0.1, -1.0, 0.0) + u_rr = SVector(0.1, -1.0, 0.0) + @test flux_hll(u_ll, u_rr, 1, equations) ≈ flux(u_rr, 1, equations) + equations = ShallowWaterEquations2D(gravity_constant = 9.81) normal_directions = [SVector(1.0, 0.0), SVector(0.0, 1.0), @@ -707,6 +718,17 @@ end @test flux_hll(u, u, normal_direction, equations) ≈ flux(u, normal_direction, equations) end + + normal_direction = SVector(1.0, 0.0, 0.0) + u_ll = SVector(0.1, 1.0, 1.0, 0.0) + u_rr = SVector(0.1, 1.0, 1.0, 0.0) + @test flux_hll(u_ll, u_rr, normal_direction, equations) ≈ + flux(u_ll, normal_direction, equations) + + u_ll = SVector(0.1, -1.0, -1.0, 0.0) + u_rr = SVector(0.1, -1.0, -1.0, 0.0) + @test flux_hll(u_ll, u_rr, normal_direction, equations) ≈ + flux(u_rr, normal_direction, equations) end @timed_testset "Consistency check for HLL flux (naive): MHD" begin diff --git a/test/test_unstructured_2d.jl b/test/test_unstructured_2d.jl index 83b8318c926..04eb9f679aa 100644 --- a/test/test_unstructured_2d.jl +++ b/test/test_unstructured_2d.jl @@ -1,7 +1,5 @@ module TestExamplesUnstructuredMesh2D -# TODO: TrixiShallowWater: move any wet/dry and two layer tests - using Test using Trixi @@ -566,105 +564,6 @@ end 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 - -@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 - # TODO: FD; for now put the unstructured tests for the 2D FDSBP here. @trixi_testset "FDSBP (central): elixir_advection_basic.jl" begin @test_trixi_include(joinpath(pkgdir(Trixi, "examples", "unstructured_2d_fdsbp"),