From fd8096a5e4858d54c0114f0be9c4da519bb515a3 Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Wed, 19 Jun 2024 10:27:56 +0200 Subject: [PATCH 1/7] Bump julia-actions/cache from 1 to 2 (#47) Bumps [julia-actions/cache](https://github.com/julia-actions/cache) from 1 to 2. - [Release notes](https://github.com/julia-actions/cache/releases) - [Changelog](https://github.com/julia-actions/cache/blob/main/devdocs/making_a_new_release.md) - [Commits](https://github.com/julia-actions/cache/compare/v1...v2) --- updated-dependencies: - dependency-name: julia-actions/cache dependency-type: direct:production update-type: version-update:semver-major ... Signed-off-by: dependabot[bot] Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com> --- .github/workflows/Documenter.yml | 2 +- .github/workflows/ci.yml | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/.github/workflows/Documenter.yml b/.github/workflows/Documenter.yml index 2e865ec..89ecdc1 100644 --- a/.github/workflows/Documenter.yml +++ b/.github/workflows/Documenter.yml @@ -33,7 +33,7 @@ jobs: with: version: '1' show-versioninfo: true - - uses: julia-actions/cache@v1 + - uses: julia-actions/cache@v2 - uses: julia-actions/julia-buildpkg@v1 env: PYTHON: "" diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 87d8d73..45237c8 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -70,7 +70,7 @@ jobs: version: ${{ matrix.version }} arch: ${{ matrix.arch }} - run: julia -e 'using InteractiveUtils; versioninfo(verbose=true)' - - uses: julia-actions/cache@v1 + - uses: julia-actions/cache@v2 - uses: julia-actions/julia-buildpkg@v1 - name: Run tests without coverage uses: julia-actions/julia-runtest@v1 @@ -106,7 +106,7 @@ jobs: version: ${{ matrix.version }} arch: ${{ matrix.arch }} - run: julia -e 'using InteractiveUtils; versioninfo(verbose=true)' - - uses: julia-actions/cache@v1 + - uses: julia-actions/cache@v2 - uses: julia-actions/julia-buildpkg@v1 - name: Run tests with coverage uses: julia-actions/julia-runtest@v1 From 7beb48faf25ce447e063540471100dcb46b1b553 Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Tue, 2 Jul 2024 06:07:26 +0200 Subject: [PATCH 2/7] Bump crate-ci/typos from 1.21.0 to 1.22.9 (#49) Bumps [crate-ci/typos](https://github.com/crate-ci/typos) from 1.21.0 to 1.22.9. - [Release notes](https://github.com/crate-ci/typos/releases) - [Changelog](https://github.com/crate-ci/typos/blob/master/CHANGELOG.md) - [Commits](https://github.com/crate-ci/typos/compare/v1.21.0...v1.22.9) --- updated-dependencies: - dependency-name: crate-ci/typos dependency-type: direct:production update-type: version-update:semver-minor ... Signed-off-by: dependabot[bot] Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com> --- .github/workflows/SpellCheck.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/SpellCheck.yml b/.github/workflows/SpellCheck.yml index 10bcc22..8f6d121 100644 --- a/.github/workflows/SpellCheck.yml +++ b/.github/workflows/SpellCheck.yml @@ -10,4 +10,4 @@ jobs: - name: Checkout Actions Repository uses: actions/checkout@v4 - name: Check spelling - uses: crate-ci/typos@v1.21.0 + uses: crate-ci/typos@v1.22.9 From 10234bba0c512a02f292da882057c1668ab899f8 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Tue, 2 Jul 2024 08:09:45 +0200 Subject: [PATCH 3/7] Allow Trixi.jl v0.8 (#48) * Update Project.toml * also adapt main Project.toml --- Project.toml | 2 +- test/Project.toml | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index 6ea3d75..1871d4c 100644 --- a/Project.toml +++ b/Project.toml @@ -16,6 +16,6 @@ LinearAlgebra = "1" MuladdMacro = "0.2.2" Static = "0.3, 0.4, 0.5, 0.6, 0.7, 0.8" StaticArrays = "1" -Trixi = "0.7" +Trixi = "0.7, 0.8" julia = "1.8" Printf = "1" diff --git a/test/Project.toml b/test/Project.toml index d3420ab..c45b5c0 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -7,7 +7,7 @@ Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" [compat] Test = "1" -Trixi = "0.7" +Trixi = "0.7, 0.8" OrdinaryDiffEq = "6.49.1" Downloads = "1" Printf = "1" From 883fc00dabdd262a5f5c5aa92559c7a9c31fc670 Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Fri, 2 Aug 2024 09:17:39 +0200 Subject: [PATCH 4/7] Bump crate-ci/typos from 1.22.9 to 1.23.6 (#50) Bumps [crate-ci/typos](https://github.com/crate-ci/typos) from 1.22.9 to 1.23.6. - [Release notes](https://github.com/crate-ci/typos/releases) - [Changelog](https://github.com/crate-ci/typos/blob/master/CHANGELOG.md) - [Commits](https://github.com/crate-ci/typos/compare/v1.22.9...v1.23.6) --- updated-dependencies: - dependency-name: crate-ci/typos dependency-type: direct:production update-type: version-update:semver-minor ... Signed-off-by: dependabot[bot] Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com> --- .github/workflows/SpellCheck.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/SpellCheck.yml b/.github/workflows/SpellCheck.yml index 8f6d121..71eea71 100644 --- a/.github/workflows/SpellCheck.yml +++ b/.github/workflows/SpellCheck.yml @@ -10,4 +10,4 @@ jobs: - name: Checkout Actions Repository uses: actions/checkout@v4 - name: Check spelling - uses: crate-ci/typos@v1.22.9 + uses: crate-ci/typos@v1.23.6 From cbc57dc13d719485d5256cc1301d74c37179eb5f Mon Sep 17 00:00:00 2001 From: Patrick Ersing <114223904+patrickersing@users.noreply.github.com> Date: Thu, 15 Aug 2024 08:13:07 +0200 Subject: [PATCH 5/7] Shallow water exner equations 1D (#51) * add SWE-Exner1D * fix docstrings * add zero velocity workaround to roe flux * fix typo * add roots * fix convergence tests * fix convergence tests 2 * fix convergence tests 3 * update test_tree_1d reference values * add tests to increase coverage * Apply suggestions from code review Co-authored-by: Andrew Winters * Modify docstrings, fix function q_s --------- Co-authored-by: Andrew Winters --- Project.toml | 6 +- .../elixir_shallowwater_exner_channel.jl | 116 ++++ ...r_shallowwater_exner_source_terms_grass.jl | 69 ++ ...xir_shallowwater_exner_source_terms_mpm.jl | 70 ++ ...elixir_shallowwater_exner_well_balanced.jl | 104 +++ src/TrixiShallowWater.jl | 4 + src/equations/equations.jl | 2 + src/equations/shallow_water_exner_1d.jl | 633 ++++++++++++++++++ test/Project.toml | 14 +- test/test_tree_1d.jl | 174 +++++ test/test_unit.jl | 68 ++ 11 files changed, 1252 insertions(+), 8 deletions(-) create mode 100644 examples/tree_1d_dgsem/elixir_shallowwater_exner_channel.jl create mode 100644 examples/tree_1d_dgsem/elixir_shallowwater_exner_source_terms_grass.jl create mode 100644 examples/tree_1d_dgsem/elixir_shallowwater_exner_source_terms_mpm.jl create mode 100644 examples/tree_1d_dgsem/elixir_shallowwater_exner_well_balanced.jl create mode 100644 src/equations/shallow_water_exner_1d.jl diff --git a/Project.toml b/Project.toml index 1871d4c..831c4d4 100644 --- a/Project.toml +++ b/Project.toml @@ -6,16 +6,18 @@ version = "0.1.0-pre" [deps] LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221" +Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" +Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665" Static = "aedffcd0-7271-4cad-89d0-dc628f76c6d3" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Trixi = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb" -Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" [compat] LinearAlgebra = "1" MuladdMacro = "0.2.2" +Printf = "1" +Roots = "2.1.6" Static = "0.3, 0.4, 0.5, 0.6, 0.7, 0.8" StaticArrays = "1" Trixi = "0.7, 0.8" julia = "1.8" -Printf = "1" diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_exner_channel.jl b/examples/tree_1d_dgsem/elixir_shallowwater_exner_channel.jl new file mode 100644 index 0000000..a62bc4d --- /dev/null +++ b/examples/tree_1d_dgsem/elixir_shallowwater_exner_channel.jl @@ -0,0 +1,116 @@ + +using OrdinaryDiffEq +using Trixi +using TrixiShallowWater +using Roots + +############################################################################### +# Semidiscretization of the shallow water exner equations for a channel flow problem +# with sediment transport + +equations = ShallowWaterExnerEquations1D(gravity_constant = 9.81, rho_f = 1.0, rho_s = 1.0, + porosity = 0.4, + sediment_model = GrassModel(A_g = 0.01)) + +# Initial condition for for a channel flow problem over a sediment hump. +# An asymptotic solution based on the method of characteristics was derived under a rigid-lid +# approximation in chapter 3.5.1 of the thesis: +# - Justin Hudson (2001) +# "Numerical Techniques for Morphodynamic Modelling" +# [PhD Thesis, University of Reading](https://citeseerx.ist.psu.edu/document?repid=rep1&type=pdf&doi=f78dbae9cfbb12ae823975c6ce9d2585b40417ba) +function initial_condition_channel(x, t, + equations::ShallowWaterExnerEquations1D) + (; sediment_model, porosity_inv) = equations + + H_ref = 10.0 # Reference water level + hv = 10.0 + + # Use the method of characteristics to find the asymptotic solution for the bed height, see + # Eq. 3.16 in the reference. + # First use an iterative method to predict x0 + f(x0) = x[1] - x0 - + sediment_model.A_g * porosity_inv * 3 * hv^3 * t * + (H_ref - sinpi((x0 - 300) / 200)^2)^(-4) + fx = Roots.ZeroProblem(f, 400.0) + x0 = Roots.solve(fx, atol = 0.0, rtol = 0.0) + + # If the result is outside 300 < x < 500 the result is invalid and instead compute x0 from + if x0 > 500 || x0 < 300 + x0 = x[1] - + sediment_model.A_g * porosity_inv * 3 * hv^3 * t * H_ref^(-4) + end + + # Compute the sediment and water height + 300 < x0 < 500 ? h_b = 0.1 + sinpi((x0 - 300) / 200)^2 : h_b = 0.1 + h = H_ref - h_b + + return SVector(h, hv, h_b) +end + +initial_condition = initial_condition_channel + +############################################################################### +# Get the DG approximation space + +volume_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal) +surface_flux = (FluxPlusDissipation(flux_ersing_etal, dissipation_roe), + flux_nonconservative_ersing_etal) + +basis = LobattoLegendreBasis(6) + +indicator_sc = IndicatorHennemannGassner(equations, basis, + alpha_max = 0.5, + alpha_min = 0.001, + alpha_smooth = true, + variable = water_sediment_height) +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 periodic mesh + +coordinates_min = 0.0 +coordinates_max = 1000.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 solver + +tspan = (0.0, 30_000.0) +ode = semidiscretize(semi, tspan) + +############################################################################### +# Callbacks + +summary_callback = SummaryCallback() + +analysis_interval = 10000 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +stepsize_callback = StepsizeCallback(cfl = 0.8) + +save_solution = SaveSolutionCallback(dt = 10_000.0, + save_initial_solution = true, + save_final_solution = true) + +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, + stepsize_callback, save_solution) + +############################################################################### +# 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_exner_source_terms_grass.jl b/examples/tree_1d_dgsem/elixir_shallowwater_exner_source_terms_grass.jl new file mode 100644 index 0000000..4389ebe --- /dev/null +++ b/examples/tree_1d_dgsem/elixir_shallowwater_exner_source_terms_grass.jl @@ -0,0 +1,69 @@ +using OrdinaryDiffEq +using Trixi +using TrixiShallowWater + +############################################################################### +# Semidiscretization of the SWE-Exner equations with source terms for convergence testing + +# Equations with Grass model +equations = ShallowWaterExnerEquations1D(gravity_constant = 10.0, rho_f = 0.5, + rho_s = 1.0, porosity = 0.5, + friction = ManningFriction(n = 0.0), + sediment_model = GrassModel(A_g = 0.01)) + +initial_condition = initial_condition_convergence_test + +############################################################################### +# Get the DG approximation space + +volume_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal) +surface_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal) + +solver = DGSEM(polydeg = 4, + surface_flux = surface_flux, + 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 = 2, + 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 = 200 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +stepsize_callback = StepsizeCallback(cfl = 1.0) + +save_solution = SaveSolutionCallback(interval = 200, + save_initial_solution = true, + save_final_solution = true) + +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, + stepsize_callback, save_solution) + +############################################################################### +# 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_exner_source_terms_mpm.jl b/examples/tree_1d_dgsem/elixir_shallowwater_exner_source_terms_mpm.jl new file mode 100644 index 0000000..9022b7e --- /dev/null +++ b/examples/tree_1d_dgsem/elixir_shallowwater_exner_source_terms_mpm.jl @@ -0,0 +1,70 @@ +using OrdinaryDiffEq +using Trixi +using TrixiShallowWater + +############################################################################### +# Semidiscretization of the SWE-Exner equations with source terms for convergence testing + +# Equations with Meyer-Peter-Mueller model +equations = ShallowWaterExnerEquations1D(gravity_constant = 10.0, rho_f = 0.5, + rho_s = 1.0, porosity = 0.5, + friction = ManningFriction(n = 0.01), + sediment_model = MeyerPeterMueller(theta_c = 0.0, + d_s = 1e-3)) + +initial_condition = initial_condition_convergence_test + +############################################################################### +# Get the DG approximation space + +volume_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal) +surface_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal) + +solver = DGSEM(polydeg = 4, + surface_flux = surface_flux, + 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 = 2, + 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 = 200 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +stepsize_callback = StepsizeCallback(cfl = 1.0) + +save_solution = SaveSolutionCallback(interval = 200, + save_initial_solution = true, + save_final_solution = true) + +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, + stepsize_callback, save_solution) + +############################################################################### +# 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_exner_well_balanced.jl b/examples/tree_1d_dgsem/elixir_shallowwater_exner_well_balanced.jl new file mode 100644 index 0000000..facc9b6 --- /dev/null +++ b/examples/tree_1d_dgsem/elixir_shallowwater_exner_well_balanced.jl @@ -0,0 +1,104 @@ + +using OrdinaryDiffEq +using Trixi +using TrixiShallowWater + +############################################################################### +# Semidiscretization of the SWE-Exner equations with a discontinuous sediment bed +# to test well-balancedness + +# Equations with Meyer-Peter-Mueller sedimentation model +equations = ShallowWaterExnerEquations1D(gravity_constant = 10.0, H0 = 1.0, + rho_f = 0.5, rho_s = 1.0, porosity = 0.5, + friction = ManningFriction(n = 0.01), + sediment_model = MeyerPeterMueller(theta_c = 0.047, + d_s = 1e-3)) + +function initial_condition_steady_state(x, t, + equations::ShallowWaterExnerEquations1D) + hv = 0.0 + + # discontinuous sediment bed + if -2.0 < x[1] < 0.0 + h_s = 0.3 + else + h_s = 0.1 + end + + h = 1.0 - h_s + + return SVector(h, hv, h_s) +end + +initial_condition = initial_condition_steady_state + +boundary_condition = boundary_condition_slip_wall + +############################################################################### +# Get the DG approximation space + +volume_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal) +surface_flux = (FluxPlusDissipation(flux_ersing_etal, DissipationLocalLaxFriedrichs()), + flux_nonconservative_ersing_etal) + +basis = LobattoLegendreBasis(3) + +indicator_sc = IndicatorHennemannGassner(equations, basis, + alpha_max = 0.5, + alpha_min = 0.001, + alpha_smooth = true, + variable = water_sediment_height) +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 + +coordinates_min = -2.0 +coordinates_max = 2.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, + source_terms = source_term_bottom_friction, + boundary_conditions = boundary_condition) + +############################################################################### +# ODE solver + +tspan = (0.0, 200.0) +ode = semidiscretize(semi, tspan) + +############################################################################### +# Callbacks + +summary_callback = SummaryCallback() + +analysis_interval = 10000 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + extra_analysis_integrals = (lake_at_rest_error,)) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(dt = 100, + save_initial_solution = true, + save_final_solution = true) + +stepsize_callback = StepsizeCallback(cfl = 1.0) + +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, + stepsize_callback) + +############################################################################### +# run the simulation +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); + +summary_callback() # print the timer summary diff --git a/src/TrixiShallowWater.jl b/src/TrixiShallowWater.jl index 5522fc5..86cd7e9 100644 --- a/src/TrixiShallowWater.jl +++ b/src/TrixiShallowWater.jl @@ -19,6 +19,7 @@ include("solvers/indicators.jl") # Export types/functions that define the public API of TrixiShallowWater.jl export ShallowWaterEquationsWetDry1D, ShallowWaterEquationsWetDry2D, + ShallowWaterExnerEquations1D, ShallowWaterTwoLayerEquations1D, ShallowWaterTwoLayerEquations2D, ShallowWaterMultiLayerEquations1D, ShallowWaterMultiLayerEquations2D @@ -26,6 +27,9 @@ export hydrostatic_reconstruction_chen_noelle, flux_nonconservative_chen_noelle, min_max_speed_chen_noelle, flux_hll_chen_noelle, flux_ersing_etal, flux_es_ersing_etal, hydrostatic_reconstruction_ersing_etal +export ManningFriction, MeyerPeterMueller, GrassModel, ShieldsStressModel, + dissipation_roe, water_sediment_height, source_term_bottom_friction + export nlayers, eachlayer export PositivityPreservingLimiterShallowWater diff --git a/src/equations/equations.jl b/src/equations/equations.jl index b2e69a6..b215b9a 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -11,6 +11,8 @@ include("shallow_water_wet_dry_1d.jl") include("shallow_water_wet_dry_2d.jl") +include("shallow_water_exner_1d.jl") + include("shallow_water_two_layer_1d.jl") include("shallow_water_two_layer_2d.jl") diff --git a/src/equations/shallow_water_exner_1d.jl b/src/equations/shallow_water_exner_1d.jl new file mode 100644 index 0000000..cd043a6 --- /dev/null +++ b/src/equations/shallow_water_exner_1d.jl @@ -0,0 +1,633 @@ +# 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 + +# Abstract type for the different bottom friction models +abstract type Friction{RealT} end + +""" + ManningFriction(; n) + +Creates a Manning friction model for the bottom friction with Manning coefficient `n`. +The type is used to dispatch on the respective friction law through the +`shear_stress_coefficient` when computing the `shear_stress`. +""" +struct ManningFriction{RealT} <: Friction{RealT} + n::RealT +end + +function ManningFriction(; n) + ManningFriction(n) +end + +# Abstract type for the different models to compute sediment discharge +abstract type SedimentModel{RealT} end + +@doc raw""" + ShieldsStressModel(; m_1, m_2, m_3, k_1, k_2, k_3, theta_c, d_s) + +Create a Shields stress model to compute the sediment discharge `q_s` based on the generalized +formulation from equation (1.2) in the given reference. + +The choice of the real constants `m_1`, `m_2`, `m_3`, `k_1`, `k_2`, and `k_3` creates +different models. For example, setting `m_1=0`, `m_2=1.5`, `m_3=0`, `k_1=8`, `k_2=1`, and `k_3=0` +yields the sedimentation model of Meyer-Peter and Müller as given in [`MeyerPeterMueller`](@ref) below. +The Shields stress represents the ratio of agitating and stabilizing forces in the sediment bed where +`theta_c` is the critical Shields stress for incipient motion and `d_s` is the mean diameter of +the sediment grain size. + +- E.D. Fernández-Nieto, T.M. de Luna, G. Narbona-Reina and J. de Dieu Zabsonré (2017)\ + Formal deduction of the Saint-Venant–Exner model including arbitrarily sloping sediment beds and + associated energy\ + [DOI: 10.1051/m2an/2016018](https://doi.org/10.1051/m2an/2016018) +""" +struct ShieldsStressModel{RealT} <: SedimentModel{RealT} + m_1::RealT + m_2::RealT + m_3::RealT + k_1::RealT + k_2::RealT + k_3::RealT + theta_c::RealT # critical shields stress + d_s::RealT # grain diameter +end + +@doc raw""" + GrassModel(; A_g, m_g=3) + +Creates a Grass model to compute the sediment discharge `q_s` as +```math +q_s = A_g v^{m_g} +``` +with the coefficients `A_g` and `m_g`. The constant `A_g` lies in the interval ``[0,1]`` +and is a dimensional calibration constant that is usually measured experimental. It expresses +the kind of interaction between the fluid and the sediment, the strength of which +increases as `A_g` approaches to 1. The factor `m_g` lies in the interval ``[1, 4]``. +Typically, one considers an odd integer value for `m_g` such that the sediment discharge +`q_s` can be differentiated and the model remains valid for all values of the velocity `v`. + +An overview of different formulations to compute the sediment discharge can be found in: +- M.J. Castro Díaz, E.D. Fernández-Nieto, A.M. Ferreiro (2008)\ + Sediment transport models in Shallow Water equations and numerical approach by high order + finite volume methods\ + [DOI:10.1016/j.compfluid.2007.07.017](https://doi.org/10.1016/j.compfluid.2007.07.017) +""" +struct GrassModel{RealT} <: SedimentModel{RealT} + A_g::RealT + m_g::RealT +end + +function GrassModel(; A_g, m_g = 3.0) + return GrassModel(A_g, m_g) +end + +@doc raw""" + MeyerPeterMueller(; theta_c, d_s) + +Creates a Meyer-Peter-Mueller model to compute the sediment discharge +`q_s` with the critical Shields stress `theta_c` and the grain diameter `d_s`. + +An overview of different formulations to compute the sediment discharge can be found in: +- M.J. Castro Díaz, E.D. Fernández-Nieto, A.M. Ferreiro (2008)\ + Sediment transport models in Shallow Water equations and numerical approach by high order + finite volume methods\ + [DOI:10.1016/j.compfluid.2007.07.017](https://doi.org/10.1016/j.compfluid.2007.07.017) +""" +function MeyerPeterMueller(; theta_c, d_s) + return ShieldsStressModel(0.0, 1.5, 0.0, 8.0, 1.0, 0.0, theta_c, d_s) +end + +@doc raw""" + ShallowWaterExnerEquations1D(;gravity_constant, H0 = 0.0, + friction = ManningFriction(n = 0.0), + sediment_model, + porosity, + rho_f, rho_s) +Formulation of the Shallow water-Exner equations in one space dimension that possesses a mathematical +entropy inequality. The equations are given by +```math +\begin{cases} +\partial_t h + \partial_x hv = 0, \\ +\partial_t hv + \partial_x (hv^2) + gh\partial_x (h + h_b) + g\frac{1}{r}h_s\partial_x (rh + h_b) + \frac{\tau}{\rho_f} = 0,\\ +\partial_t h_b + \partial_x q_s = 0, +\end{cases} +``` +The unknown quantities are the water and sediment height ``h``, ``h_b`` and the velocity ``v``. +The sediment discharge ``q_s(h, hv)`` is determined by the `sediment_model` and is used to determine +the active sediment height ``h_s = q_s / v``. +Furthermore ``\tau`` denotes the shear stress at the water-sediment interface and is determined by +the `friction` model. +The gravitational constant is denoted by ``g``, and ``\rho_f`` and ``\rho_s`` are the fluid and sediment +densities, respectively. The density ratio is given by ``r = \rho_f / \rho_s``, where ``r`` lies between +``0 < r < 1`` as the fluid density ``\rho_f`` should be smaller than the sediment density ``\rho_s``. + +The conservative variable water height ``h`` is measured from the sediment height ``h_b``, therefore +one also defines the total water height as ``H = 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. + +The entropy conservative formulation has been derived in the paper: +- E.D. Fernández-Nieto, T.M. de Luna, G. Narbona-Reina and J. de Dieu Zabsonré (2017)\ + Formal deduction of the Saint-Venant–Exner model including arbitrarily sloping sediment beds and + associated energy\ + [DOI: 10.1051/m2an/2016018](https://doi.org/10.1051/m2an/2016018) +""" +struct ShallowWaterExnerEquations1D{RealT <: Real, + FrictionT <: Friction{RealT}, + SedimentT <: SedimentModel{RealT}} <: + Trixi.AbstractShallowWaterEquations{1, 3} + # This structure ensures that friction and sediment models are concrete types + # to prevent allocations + gravity::RealT # gravitational constant + H0::RealT # constant "lake-at-rest" total water height + friction::FrictionT + sediment_model::SedimentT + porosity_inv::RealT # 1/(1-porosity) + rho_f::RealT # density of fluid + rho_s::RealT # density of sediment + r::RealT # density ratio +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. +function ShallowWaterExnerEquations1D(; gravity_constant, H0 = zero(gravity_constant), + friction = ManningFriction(n = 0.0), + sediment_model, + porosity, rho_f, rho_s) + # Precompute common expressions for the porosity and density ratio + porosity_inv = inv(1 - porosity) + r = rho_f / rho_s + return ShallowWaterExnerEquations1D(gravity_constant, H0, friction, sediment_model, + porosity_inv, rho_f, rho_s, r) +end + +Trixi.have_nonconservative_terms(::ShallowWaterExnerEquations1D) = True() +Trixi.varnames(::typeof(cons2cons), ::ShallowWaterExnerEquations1D) = ("h", "hv", "h_b") +# Note, we use the total water height, H = h + h_b, as the first primitive variable for easier +# visualization and setting initial conditions +Trixi.varnames(::typeof(cons2prim), ::ShallowWaterExnerEquations1D) = ("H", "v", "h_b") + +@doc raw""" + boundary_condition_slip_wall(u_inner, orientation_or_normal, x, t, surface_flux_function, + equations::ShallowWaterExnerEquations1D) + +Create a boundary state by reflecting the normal velocity component and keep +the tangential velocity component unchanged. The boundary water height is taken from +the internal value. + +For details see Section 9.2.5 of the book: +- Eleuterio F. Toro (2001)\ + Shock-Capturing Methods for Free-Surface Shallow Flows\ + 1st edition\ + ISBN 0471987662 +""" +@inline function Trixi.boundary_condition_slip_wall(u_inner, orientation_or_normal, + direction, + x, t, + surface_flux_function, + equations::ShallowWaterExnerEquations1D) + # create the "external" boundary solution state + u_boundary = SVector(u_inner[1], + -u_inner[2], + u_inner[3]) + + # calculate the boundary flux + if iseven(direction) # u_inner is "left" of boundary, u_boundary is "right" of boundary + flux = surface_flux_function(u_inner, u_boundary, orientation_or_normal, + equations) + else # u_boundary is "left" of boundary, u_inner is "right" of boundary + flux = surface_flux_function(u_boundary, u_inner, orientation_or_normal, + equations) + end + + return flux +end + +# Set initial conditions at physical location `x` for time `t` +""" + initial_condition_convergence_test(x, t, equations::ShallowWaterExnerEquations1D) + +A smooth initial condition used for convergence tests in combination with +[`Trixi.source_terms_convergence_test`](@ref). +""" +@inline function Trixi.initial_condition_convergence_test(x, t, + equations::ShallowWaterExnerEquations1D) + ω = sqrt(2) * pi + + h = 2.0 + cos(ω * x[1]) * cos(ω * t) + v = 0.5 + h_b = 2.0 + sin(ω * x[1]) * cos(ω * t) + + return SVector(h, h * v, h_b) +end + +""" + source_terms_convergence_test(u, x, t, equations::ShallowWaterExnerEquations1D{T, S, GrassModel{T}}) where {T, S} + +Source terms used for convergence tests in combination with [`Trixi.initial_condition_convergence_test`](@ref) +when using the the [`GrassModel`](@ref) model. + +To use this source term the equations must be set to: +```julia +equations = ShallowWaterExnerEquations1D(gravity_constant = 10.0, rho_f = 0.5, + rho_s = 1.0, porosity = 0.5, + friction = ManningFriction(n = 0.0), + sediment_model = GrassModel(A_g = 0.01) +``` +""" +@inline function Trixi.source_terms_convergence_test(u, x, t, + equations::ShallowWaterExnerEquations1D{T, + S, + GrassModel{T}}) where { + T, + S + } + ω = sqrt(2.0) * pi + A_g = equations.sediment_model.A_g + + h = -cos(x[1] * ω) * sin(t * ω) * ω - 0.5 * sin(x[1] * ω) * cos(t * ω) * ω + hv = -0.5 * cos(x[1] * ω) * sin(t * ω) * ω - 0.25 * sin(x[1] * ω) * cos(t * ω) * ω + + 10.0 * A_g * + (cos(x[1] * ω) * cos(t * ω) * ω - 0.5 * sin(x[1] * ω) * cos(t * ω) * ω) + + 10.0 * (2.0 + cos(x[1] * ω) * cos(t * ω)) * + (cos(x[1] * ω) * cos(t * ω) * ω - sin(x[1] * ω) * cos(t * ω) * ω) + h_b = -sin(x[1] * ω) * sin(t * ω) * ω + return SVector(h, hv, h_b) +end + +""" + source_terms_convergence_test(u, x, t, equations::ShallowWaterExnerEquations1D{T, S, ShieldsStressModel{T}}) where {T, S} + +Source terms used for convergence tests in combination with [`Trixi.initial_condition_convergence_test`](@ref) +when using the [`MeyerPeterMueller`](@ref) model. + +To use this source term the equations must be set to: +```julia +equations = ShallowWaterExnerEquations1D(gravity_constant = 10.0, rho_f = 0.5, + rho_s = 1.0, porosity = 0.5, + friction = ManningFriction(n = 0.01), + sediment_model = MeyerPeterMueller(theta_c = 0.0, + d_s = 1e-3)) +``` +""" +@inline function Trixi.source_terms_convergence_test(u, x, t, + equations::ShallowWaterExnerEquations1D{T, + S, + ShieldsStressModel{T}}) where { + T, + S + } + ω = sqrt(2.0) * pi + (; gravity, porosity_inv, rho_f, rho_s, r) = equations + + n = equations.friction.n + + # Constant expression from the MPM model + c = sqrt(gravity * (1 / r - 1)) * 8.0 * porosity_inv * + (rho_f / (rho_s - rho_f))^(3 / 2) * n^3 + + h = -cos(x[1] * ω) * sin(t * ω) * ω - 0.5 * sin(x[1] * ω) * cos(t * ω) * ω + + hv = ((5.0 * c * + (cos(x[1] * ω) * cos(t * ω) * ω - 0.5 * sin(x[1] * ω) * cos(t * ω) * ω)) / + ((2.0 + cos(x[1] * ω) * cos(t * ω))^0.5) - + 0.5 * cos(x[1] * ω) * sin(t * ω) * ω - + 0.25 * sin(x[1] * ω) * cos(t * ω) * ω + + 10.0 * (2.0 + cos(x[1] * ω) * cos(t * ω)) * + (cos(x[1] * ω) * cos(t * ω) * ω - sin(x[1] * ω) * cos(t * ω) * ω)) + + h_b = ((0.5 * ((0.125 * c) / (2.0 + cos(x[1] * ω) * cos(t * ω))) * sin(x[1] * ω) * + cos(t * ω) * ω) / ((2.0 + cos(x[1] * ω) * cos(t * ω))^0.5) - + sin(x[1] * ω) * sin(t * ω) * ω) + + return SVector(h, hv, h_b) +end + +""" + source_term_bottom_friction(u, x, t, equations::ShallowWaterExnerEquations1D) + +Source term that accounts for the bottom friction in the [ShallowWaterExnerEquations1D](@ref). +The actual friction law is determined through the friction model in `equations.friction`. +""" +@inline function source_term_bottom_friction(u, x, t, + equations::ShallowWaterExnerEquations1D) + return SVector(zero(eltype(u)), + -shear_stress(u, equations), + zero(eltype(u))) +end + +# Calculate 1D flux for a single point +@inline function Trixi.flux(u, orientation::Integer, + equations::ShallowWaterExnerEquations1D) + _, hv, _ = u + + v = velocity(u, equations) + + f1 = hv + f2 = hv * v + f3 = q_s(u, equations) + + return SVector(f1, f2, f3) +end + +""" + flux_nonconservative_ersing_etal(u_ll, u_rr, orientation, equations::ShallowWaterExnerEquations1D) + +Non-symmetric path-conservative two-point flux discretizing the nonconservative terms of the +[`ShallowWaterExnerEquations1D`](@ref) which consists of the hydrostatic pressure of the fluid +layer and an additional pressure contribution from the sediment layer to obtain an entropy +inequality. + +This non-conservative flux should be used together with [`flux_ersing_etal`](@ref) to create a +scheme that is entropy conservative and well-balanced. +""" +@inline function Trixi.flux_nonconservative_ersing_etal(u_ll, u_rr, + orientation::Integer, + equations::ShallowWaterExnerEquations1D) + # Pull the necessary left and right state information + h_ll, _, h_b_ll = u_ll + h_rr, _, h_b_rr = u_rr + + # Calculate jumps + h_jump = h_rr - h_ll + h_b_jump = h_b_rr - h_b_ll + + # Workaround to avoid division by zero, when computing the effective sediment height + if velocity(u_ll, equations) < eps() + h_s_ll = 0.0 + else + h_s_ll = q_s(u_ll, equations) / velocity(u_ll, equations) + end + + z = zero(eltype(u_ll)) + + f2 = equations.gravity * h_ll * (h_jump + h_b_jump) + # Additional nonconservative term to obtain entropy conservative formulation + f2 += (equations.gravity / equations.r * h_s_ll * + (equations.r * h_jump + h_b_jump)) + + return SVector(z, f2, z) +end + +""" + flux_ersing_etal(u_ll, u_rr, orientation::Integer, + equations::ShallowWaterMultiLayerEquations1D) + +Entropy conservative split form, without the hydrostatic pressure. This flux should be used +together with the nonconservative [`flux_nonconservative_ersing_etal`](@ref) to create a scheme +that is entropy conservative and well-balanced. + +To obtain an entropy stable formulation the `surface_flux` can be set as +`FluxPlusDissipation(flux_ersing_etal, DissipationLocalLaxFriedrichs()), flux_nonconservative_ersing_etal`. +""" +@inline function flux_ersing_etal(u_ll, u_rr, orientation::Integer, + equations::ShallowWaterExnerEquations1D) + # Unpack left and right state + _, h_v_ll, _ = u_ll + _, h_v_rr, _ = u_rr + + # Get the velocities on either side + v_ll = velocity(u_ll, equations) + v_rr = velocity(u_rr, equations) + + # Average each factor of products in flux + v_avg = 0.5 * (v_ll + v_rr) + + # Calculate fluxes depending on orientation + f1 = 0.5 * (h_v_ll + h_v_rr) + f2 = f1 * v_avg + f3 = 0.5 * (q_s(u_ll, equations) + q_s(u_rr, equations)) + + return SVector(f1, f2, f3) +end + +""" + dissipation_roe(u_ll, u_rr, orientation_or_normal_direction, + equations::ShallowWaterExnerEquations1D) +Roe-type dissipation term for the [`ShallowWaterExnerEquations1D`](@ref) with an approximate Roe average +for the sediment discharge `q_s`. +""" +@inline function dissipation_roe(u_ll, u_rr, orientation_or_normal_direction, + equations::ShallowWaterExnerEquations1D) + r = equations.r + g = equations.gravity + z = zero(eltype(u_ll)) + + # Get entropy variables and velocities + v_ll = velocity(u_ll, equations) + v_rr = velocity(u_rr, equations) + + # Compute approximate Roe averages. + # The actual Roe average for the sediment discharge `q_s` would depend on the sediment and + # friction model and can be difficult to compute analytically. + # Therefore we only use an approximation here. + h_avg = 0.5 * (u_ll[1] + u_rr[1]) + v_avg = (sqrt(u_ll[1]) * v_ll + sqrt(u_rr[1]) * v_rr) / + (sqrt(u_ll[1]) + sqrt(u_rr[1])) + # Workaround to avoid division by zero, when computing the effective sediment height + if v_avg < eps() + h_s_avg = 0.0 + else + h_s_avg = 0.5 * (q_s(u_ll, equations) / v_ll + q_s(u_rr, equations) / v_rr) + end + + # Compute the eigenvalues using Cardano's formula + λ1, λ2, λ3 = eigvals_cardano(SVector(h_avg, h_avg * v_avg, h_s_avg), equations) + + # Precompute some common expressions + c1 = g * (h_avg + h_s_avg) + c2 = g * (h_avg + h_s_avg / r) + + # Eigenvector matrix + r31 = ((v_avg - λ1)^2 - c1) / c2 + r32 = ((v_avg - λ2)^2 - c1) / c2 + r33 = ((v_avg - λ3)^2 - c1) / c2 + R = @SMatrix [[1 1 1]; [λ1 λ2 λ3]; [r31 r32 r33]] + + # Inverse eigenvector matrix + d1 = (λ1 - λ2) * (λ1 - λ3) + d2 = (λ2 - λ1) * (λ2 - λ3) + d3 = (λ3 - λ2) * (λ3 - λ1) + R_inv = @SMatrix [(c1 - v_avg^2 + λ2 * λ3)/d1 (2 * v_avg - λ2 - λ3)/d1 c2/d1; + (c1 - v_avg^2 + λ1 * λ3)/d2 (2 * v_avg - λ1 - λ3)/d2 c2/d2; + (c1 - v_avg^2 + λ1 * λ2)/d3 (2 * v_avg - λ2 - λ1)/d3 c2/d3] + + # Eigenvalue absolute value matrix + Λ_abs = @SMatrix [abs(λ1) z z; z abs(λ2) z; z z abs(λ3)] + + # Compute dissipation + diss = SVector(-0.5 * R * Λ_abs * R_inv * (u_rr - u_ll)) + + return SVector(diss[1], diss[2], diss[3]) +end + +@inline function Trixi.max_abs_speed_naive(u_ll, u_rr, orientation::Integer, + equations::ShallowWaterExnerEquations1D) + return max(eigvals_cardano(u_rr, equations)..., eigvals_cardano(u_ll, equations)...) +end + +@inline function Trixi.max_abs_speeds(u, equations::ShallowWaterExnerEquations1D) + return maximum(eigvals_cardano(u, equations)) +end + +#Helper function to extract the velocity vector from the conservative variables +@inline function Trixi.velocity(u, equations::ShallowWaterExnerEquations1D) + h, hv, _ = u + + v = hv / h + + return v +end + +# Compute the sediment discharge for Shields stress models +@inline function q_s(u, + equations::ShallowWaterExnerEquations1D{T, S, + ShieldsStressModel{T}}) where { + T, + S + } + (; gravity, rho_f, rho_s, porosity_inv) = equations + (; m_1, m_2, m_3, k_1, k_2, k_3, theta_c, d_s) = equations.sediment_model + + theta = rho_f * abs(shear_stress(u, equations)) / (gravity * (rho_s - rho_f) * d_s) # Shields stress + + Q = d_s * sqrt(gravity * (rho_s / rho_f - 1.0) * d_s) # Characteristic discharge + + return (porosity_inv * Q * sign(theta) * k_1 * theta^m_1 * + (max(theta - k_2 * theta_c, 0.0))^m_2 * + (max(sqrt(theta) - k_3 * sqrt(theta_c), 0.0))^m_3) +end + +# Compute the sediment discharge for the Grass model +@inline function q_s(u, + equations::ShallowWaterExnerEquations1D{T, S, GrassModel{T}}) where { + T, + S + } + (; porosity_inv, sediment_model) = equations + return porosity_inv * sediment_model.A_g * velocity(u, equations)^sediment_model.m_g +end + +# Shear stress formulation using a coefficient to take into account different friction models +@inline function shear_stress(u, equations::ShallowWaterExnerEquations1D) + v = velocity(u, equations) + return equations.gravity * shear_stress_coefficient(u, equations.friction) * v * + abs(v) +end + +# Model dependent shear stress coefficient +@inline function shear_stress_coefficient(u, friction::ManningFriction) + h, _, _ = u + return friction.n^2 / h^(1 / 3) +end + +# Convert conservative variables to primitive +@inline function Trixi.cons2prim(u, equations::ShallowWaterExnerEquations1D) + h, _, h_b = u + + H = h + h_b + v = velocity(u, equations) + + return SVector(H, v, h_b) +end + +# Convert conservative variables to entropy variables +@inline function Trixi.cons2entropy(u, equations::ShallowWaterExnerEquations1D) + h, _, h_b = u + (; gravity, r) = equations + + v = velocity(u, equations) + + w1 = r * (gravity * (h + h_b) - 0.5 * v^2) + w2 = r * v + w3 = gravity * (r * h + h_b) + + return SVector(w1, w2, w3) +end + +# Convert primitive to conservative variables +@inline function Trixi.prim2cons(prim, equations::ShallowWaterExnerEquations1D) + H, v, h_b = prim + + h = H - h_b + hv = h * v + + return SVector(h, hv, h_b) +end + +@inline function Trixi.waterheight(u, equations::ShallowWaterExnerEquations1D) + return u[1] +end + +# Indicator variable used for the shock capturing in `IndicatorHennemannGassnerShallowWater` +@inline function water_sediment_height(u, equations::ShallowWaterExnerEquations1D) + return equations.gravity * u[1] * u[3] +end + +# Mathematical entropy function +@inline function Trixi.entropy(u, equations::ShallowWaterExnerEquations1D) + h, _, h_b = u + v = velocity(u, equations) + (; gravity, r) = equations + + return 0.5 * r * h * v^2 + 0.5 * gravity * (r * h^2 + h_b^2) + r * gravity * h * h_b +end + +# Calculate the error for the "lake-at-rest" test case where H = h + h_b should +# be a constant value over time. +@inline function Trixi.lake_at_rest_error(u, equations::ShallowWaterExnerEquations1D) + h, _, h_b = u + return abs(equations.H0 - (h + h_b)) +end + +# Trigonometric version of Cardano's method for the roots of a cubic polynomial +# in order to compute the eigenvalues of the [`ShallowWaterExnerEquations1D`[(@ref)]. +# Note, assumes only real roots. +@inline function eigvals_cardano(u, equations::ShallowWaterExnerEquations1D) + h = waterheight(u, equations) + v = velocity(u, equations) + g = equations.gravity + r = equations.r + + # Workaround to avoid division by zero, when computing the effective sediment height + if v > eps() + h_s = q_s(u, equations) / v + # Compute gradients of q_s using automatic differentiation. + # Introduces a closure to make q_s a function of u only. This is necessary since the + # gradient function only accepts functions of one variable. + dq_s_dh, dq_s_dhv, _ = Trixi.ForwardDiff.gradient(u -> q_s(u, equations), u) + else + h_s = 0.0 + dq_s_dh = 0.0 + dq_s_dhv = 0.0 + end + + # Coefficient for the original cubic equation ax^3 + bx^2 + cx + d + a = -1 + b = 2 * v + c = g * (h + 1 / r * h_s) * dq_s_dhv + g * (h + h_s) - v^2 + d = g * (h + 1 / r * h_s) * dq_s_dh + + # Coefficient of the depressed cubic equation t^3 + pt + q = 0 + p = (3 * a * c - b^2) / (3 * a^2) + q = (2 * b^3 - 9 * a * b * c + 27 * a^2 * d) / (27 * a^3) + + # Roots of the original cubic equation + λ1 = -b / (3 * a) + + 2 * sqrt(-p / 3) * cos(1 / 3 * acos(3 * q / (2 * p) * sqrt(-3 / p))) + λ2 = -b / (3 * a) + + 2 * sqrt(-p / 3) * + cos(1 / 3 * acos(3 * q / (2 * p) * sqrt(-3 / p)) - 2 * π * 1 / 3) + λ3 = -b / (3 * a) + + 2 * sqrt(-p / 3) * + cos(1 / 3 * acos(3 * q / (2 * p) * sqrt(-3 / p)) - 2 * π * 2 / 3) + + return SVector(λ1, λ2, λ3) +end +end # @muladd diff --git a/test/Project.toml b/test/Project.toml index c45b5c0..172eaa4 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -1,16 +1,18 @@ [deps] -Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" -Trixi = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb" -OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6" +OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" +Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +Trixi = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb" [compat] -Test = "1" -Trixi = "0.7, 0.8" -OrdinaryDiffEq = "6.49.1" Downloads = "1" +OrdinaryDiffEq = "6.49.1" Printf = "1" +Roots = "2.1.6" +Test = "1" +Trixi = "0.7, 0.8" [preferences.OrdinaryDiffEq] PrecompileAutoSpecialize = false diff --git a/test/test_tree_1d.jl b/test/test_tree_1d.jl index d7cde02..64752ac 100644 --- a/test/test_tree_1d.jl +++ b/test/test_tree_1d.jl @@ -823,6 +823,180 @@ end # 2LSWE end end end # MLSWE + +@testset "Shallow Water - Exner" begin + @trixi_testset "elixir_shallowwater_exner_source_terms_grass.jl with EC fluxes" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_exner_source_terms_grass.jl"), + l2=[ + 0.0004102960441666415, + 0.0024123111823754154, + 2.855259772927741e-5, + ], + linf=[ + 0.0008005791155958342, + 0.0075032017005062235, + 4.7151297207559395e-5, + ]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end + + @trixi_testset "elixir_shallowwater_exner_source_terms_mpm.jl with Roe dissipation" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_exner_source_terms_mpm.jl"), + l2=[ + 8.314340397306541e-5, + 0.0003737050980420925, + 2.81406288308791e-5, + ], + linf=[ + 0.000319905497986106, + 0.001710420951723135, + 4.494237163465975e-5, + ], + surface_flux=(FluxPlusDissipation(flux_central, + dissipation_roe), + flux_nonconservative_ersing_etal)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end + + @trixi_testset "elixir_shallowwater_exner_source_terms_mpm.jl with LLF dissipation" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_exner_source_terms_mpm.jl"), + l2=[ + 8.494087939853228e-5, + 0.00037012479603853885, + 2.8139644904971178e-5, + ], + linf=[ + 0.0003106352175714644, + 0.0016775444674146378, + 4.4937259775723604e-5, + ], + surface_flux=(FluxPlusDissipation(flux_ersing_etal, + DissipationLocalLaxFriedrichs()), + flux_nonconservative_ersing_etal)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end + + @trixi_testset "elixir_shallowwater_exner_well_balanced.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_exner_well_balanced.jl"), + l2=[ + 0.0204167335836456, + 1.5412805715405497e-15, + 0.020416733583645628, + ], + linf=[ + 0.19999999999999907, + 2.957304143715833e-15, + 0.19999999999999998, + ]) + # 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_exner_well_balanced.jl with Roe dissipation" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_exner_well_balanced.jl"), + l2=[ + 0.010910894511791577, + 1.8877123431891935e-15, + 0.010910894511791757, + ], + linf=[ + 0.19999999999999674, + 3.752915460516524e-15, + 0.19999999999999998, + ], + surface_flux=(FluxPlusDissipation(flux_ersing_etal, + dissipation_roe), + flux_nonconservative_ersing_etal)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end + + @trixi_testset "elixir_shallowwater_exner_channel.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_exner_channel.jl"), + l2=[ + 0.061254947613784645, + 0.0042920880585939165, + 0.06170746499938789, + ], + linf=[ + 0.5552555774807875, + 0.009352028888004682, + 0.549962205546136, + ]) + # 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_exner_channel.jl with EC fluxes" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_shallowwater_exner_channel.jl"), + l2=[ + 0.06511905620487073, + 0.021756041726745886, + 0.06488442680626008, + ], + linf=[ + 0.4644690482223428, + 0.058166305725594114, + 0.4604211411585378, + ], + surface_flux=(flux_ersing_etal, + flux_nonconservative_ersing_etal)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end + end +end # SWE-Exner end # TreeMesh1D # Clean up afterwards: delete TrixiShallowWater.jl output directory diff --git a/test/test_unit.jl b/test/test_unit.jl index 64002cf..e42157f 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -134,6 +134,74 @@ end end end +@timed_testset "SWE-Exner conversion between conservative/entropy variables" begin + h, v, h_b = (1.0, 0.3, 0.1) + + let equations = ShallowWaterExnerEquations1D(gravity_constant = 9.81, rho_f = 0.9, + rho_s = 1.0, porosity = 0.4, + sediment_model = GrassModel(A_g = 0.01)) + # Test conversion between primitive and conservative variables + prim_vars = SVector(h, v, h_b) + cons_vars = prim2cons(prim_vars, equations) + @test prim_vars ≈ cons2prim(cons_vars, equations) + + # Test conversion from conservative to entropy variables + entropy_vars = cons2entropy(cons_vars, equations) + @test entropy_vars ≈ + Trixi.ForwardDiff.gradient(u -> entropy(u, equations), cons_vars) + end +end + +@timed_testset "SWE-Exner eigenvalue / eigenvector computation" begin + h, v, h_b = (1.0, 0.3, 0.1) + + let equations = ShallowWaterExnerEquations1D(gravity_constant = 9.81, rho_f = 0.9, + rho_s = 1.0, porosity = 0.4, + sediment_model = GrassModel(A_g = 0.01)) + u = SVector(h, h * v, h_b) + r = equations.r + g = equations.gravity + + # Compute effective sediment height + h_s = TrixiShallowWater.q_s(SVector(h, h * v, 0.0), equations) / v + dq_s_dh, dq_s_dhv, _ = Trixi.ForwardDiff.gradient(u -> TrixiShallowWater.q_s(u, + equations), + u) + + # flux Jacobian + A = [0 1 0; (g * (h + h_s)-v^2) (2*v) (g*(h + 1 / equations.r * h_s)); + dq_s_dh dq_s_dhv 0] + + # Compute the eigenvalues using Cardano's formula + λ1, λ2, λ3 = TrixiShallowWater.eigvals_cardano(SVector(h, h * v, h_b), + equations) + + # Precompute some common expressions + c1 = g * (h + h_s) + c2 = g * (h + h_s / r) + + # Eigenvector matrix + r31 = ((v - λ1)^2 - c1) / c2 + r32 = ((v - λ2)^2 - c1) / c2 + r33 = ((v - λ3)^2 - c1) / c2 + R = [[1 1 1]; [λ1 λ2 λ3]; [r31 r32 r33]] + + # Inverse eigenvector matrix + d1 = (λ1 - λ2) * (λ1 - λ3) + d2 = (λ2 - λ1) * (λ2 - λ3) + d3 = (λ3 - λ2) * (λ3 - λ1) + R_inv = [(c1 - v^2 + λ2 * λ3)/d1 (2 * v - λ2 - λ3)/d1 c2/d1; + (c1 - v^2 + λ1 * λ3)/d2 (2 * v - λ1 - λ3)/d2 c2/d2; + (c1 - v^2 + λ1 * λ2)/d3 (2 * v - λ2 - λ1)/d3 c2/d3] + + # Eigenvalue vale matrix + Λ = [λ1 0 0; 0 λ2 0; 0 0 λ3] + + @test R * R_inv ≈ [1 0 0; 0 1 0; 0 0 1] + @test A ≈ R * Λ * R_inv + end +end + @timed_testset "Consistency check for waterheight_pressure" begin H, v1, v2, b = 3.5, 0.25, 0.1, 0.4 From 38def646d9d2aad94a048c758e0585cfd8aa183b Mon Sep 17 00:00:00 2001 From: Patrick Ersing <114223904+patrickersing@users.noreply.github.com> Date: Fri, 16 Aug 2024 10:22:00 +0200 Subject: [PATCH 6/7] limit trixi compat bound to latest version (#52) --- Project.toml | 2 +- test/Project.toml | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index 831c4d4..1b79a8e 100644 --- a/Project.toml +++ b/Project.toml @@ -19,5 +19,5 @@ Printf = "1" Roots = "2.1.6" Static = "0.3, 0.4, 0.5, 0.6, 0.7, 0.8" StaticArrays = "1" -Trixi = "0.7, 0.8" +Trixi = "0.7 - 0.8.6" julia = "1.8" diff --git a/test/Project.toml b/test/Project.toml index 172eaa4..a3f6ad8 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -12,7 +12,7 @@ OrdinaryDiffEq = "6.49.1" Printf = "1" Roots = "2.1.6" Test = "1" -Trixi = "0.7, 0.8" +Trixi = "0.7 - 0.8.6" [preferences.OrdinaryDiffEq] PrecompileAutoSpecialize = false From 846f4a1d727b66d796f44cc2e405a9506c48b344 Mon Sep 17 00:00:00 2001 From: Patrick Ersing <114223904+patrickersing@users.noreply.github.com> Date: Mon, 19 Aug 2024 16:02:47 +0200 Subject: [PATCH 7/7] update flux_nonconservative_ersing_etal (#53) --- Project.toml | 2 +- .../tree_2d_dgsem/elixir_shallowwater_wall.jl | 4 +- src/TrixiShallowWater.jl | 3 +- src/equations/shallow_water_exner_1d.jl | 6 +- src/equations/shallow_water_multilayer_1d.jl | 6 +- src/equations/shallow_water_multilayer_2d.jl | 14 +-- src/equations/shallow_water_two_layer_1d.jl | 6 +- src/equations/shallow_water_two_layer_2d.jl | 14 +-- src/equations/shallow_water_wet_dry_1d.jl | 51 +++------ src/equations/shallow_water_wet_dry_2d.jl | 100 ++++-------------- test/Project.toml | 2 +- test/test_tree_1d.jl | 32 +++--- test/test_tree_2d.jl | 26 ++--- test/test_unstructured_2d.jl | 26 ++--- 14 files changed, 97 insertions(+), 195 deletions(-) diff --git a/Project.toml b/Project.toml index 1b79a8e..b6adc02 100644 --- a/Project.toml +++ b/Project.toml @@ -19,5 +19,5 @@ Printf = "1" Roots = "2.1.6" Static = "0.3, 0.4, 0.5, 0.6, 0.7, 0.8" StaticArrays = "1" -Trixi = "0.7 - 0.8.6" +Trixi = "0.8.7" julia = "1.8" diff --git a/examples/tree_2d_dgsem/elixir_shallowwater_wall.jl b/examples/tree_2d_dgsem/elixir_shallowwater_wall.jl index 14bb8c9..39d6009 100644 --- a/examples/tree_2d_dgsem/elixir_shallowwater_wall.jl +++ b/examples/tree_2d_dgsem/elixir_shallowwater_wall.jl @@ -30,8 +30,8 @@ boundary_condition = boundary_condition_slip_wall ############################################################################### # Get the DG approximation space -volume_flux = (flux_wintermeyer_etal, flux_nonconservative_ersing_etal) -surface_flux = (flux_lax_friedrichs, flux_nonconservative_ersing_etal) +volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) +surface_flux = (flux_lax_friedrichs, flux_nonconservative_wintermeyer_etal) solver = DGSEM(polydeg = 3, surface_flux = surface_flux, volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) diff --git a/src/TrixiShallowWater.jl b/src/TrixiShallowWater.jl index 86cd7e9..3fbdec9 100644 --- a/src/TrixiShallowWater.jl +++ b/src/TrixiShallowWater.jl @@ -25,7 +25,8 @@ export ShallowWaterEquationsWetDry1D, ShallowWaterEquationsWetDry2D, export hydrostatic_reconstruction_chen_noelle, flux_nonconservative_chen_noelle, min_max_speed_chen_noelle, flux_hll_chen_noelle, - flux_ersing_etal, flux_es_ersing_etal, hydrostatic_reconstruction_ersing_etal + flux_ersing_etal, flux_nonconservative_ersing_etal, + flux_es_ersing_etal, hydrostatic_reconstruction_ersing_etal export ManningFriction, MeyerPeterMueller, GrassModel, ShieldsStressModel, dissipation_roe, water_sediment_height, source_term_bottom_friction diff --git a/src/equations/shallow_water_exner_1d.jl b/src/equations/shallow_water_exner_1d.jl index cd043a6..3aa34c9 100644 --- a/src/equations/shallow_water_exner_1d.jl +++ b/src/equations/shallow_water_exner_1d.jl @@ -347,9 +347,9 @@ inequality. This non-conservative flux should be used together with [`flux_ersing_etal`](@ref) to create a scheme that is entropy conservative and well-balanced. """ -@inline function Trixi.flux_nonconservative_ersing_etal(u_ll, u_rr, - orientation::Integer, - equations::ShallowWaterExnerEquations1D) +@inline function flux_nonconservative_ersing_etal(u_ll, u_rr, + orientation::Integer, + equations::ShallowWaterExnerEquations1D) # Pull the necessary left and right state information h_ll, _, h_b_ll = u_ll h_rr, _, h_b_rr = u_rr diff --git a/src/equations/shallow_water_multilayer_1d.jl b/src/equations/shallow_water_multilayer_1d.jl index 20067e6..d9e18b2 100644 --- a/src/equations/shallow_water_multilayer_1d.jl +++ b/src/equations/shallow_water_multilayer_1d.jl @@ -310,9 +310,9 @@ In the two-layer setting this combination is equivalent to the fluxes in: curvilinear meshes [DOI: 10.1007/s10915-024-02451-2](https://doi.org/10.1007/s10915-024-02451-2) """ -@inline function Trixi.flux_nonconservative_ersing_etal(u_ll, u_rr, - orientation::Integer, - equations::ShallowWaterMultiLayerEquations1D) +@inline function flux_nonconservative_ersing_etal(u_ll, u_rr, + orientation::Integer, + equations::ShallowWaterMultiLayerEquations1D) # Pull the necessary left and right state information h_ll = waterheight(u_ll, equations) h_rr = waterheight(u_rr, equations) diff --git a/src/equations/shallow_water_multilayer_2d.jl b/src/equations/shallow_water_multilayer_2d.jl index b541b0b..555f18b 100644 --- a/src/equations/shallow_water_multilayer_2d.jl +++ b/src/equations/shallow_water_multilayer_2d.jl @@ -467,9 +467,9 @@ In the two-layer setting this combination is equivalent to the fluxes in: curvilinear meshes [DOI: 10.1007/s10915-024-02451-2](https://doi.org/10.1007/s10915-024-02451-2) """ -@inline function Trixi.flux_nonconservative_ersing_etal(u_ll, u_rr, - orientation::Integer, - equations::ShallowWaterMultiLayerEquations2D) +@inline function flux_nonconservative_ersing_etal(u_ll, u_rr, + orientation::Integer, + equations::ShallowWaterMultiLayerEquations2D) # Pull the necessary left and right state information h_ll = waterheight(u_ll, equations) h_rr = waterheight(u_rr, equations) @@ -508,10 +508,10 @@ In the two-layer setting this combination is equivalent to the fluxes in: return SVector(f) end -@inline function Trixi.flux_nonconservative_ersing_etal(u_ll, u_rr, - normal_direction_ll::AbstractVector, - normal_direction_average::AbstractVector, - equations::ShallowWaterMultiLayerEquations2D) +@inline function flux_nonconservative_ersing_etal(u_ll, u_rr, + normal_direction_ll::AbstractVector, + normal_direction_average::AbstractVector, + equations::ShallowWaterMultiLayerEquations2D) # Pull the necessary left and right state information h_ll = waterheight(u_ll, equations) h_rr = waterheight(u_rr, equations) diff --git a/src/equations/shallow_water_two_layer_1d.jl b/src/equations/shallow_water_two_layer_1d.jl index 7e6276c..a3a0a83 100644 --- a/src/equations/shallow_water_two_layer_1d.jl +++ b/src/equations/shallow_water_two_layer_1d.jl @@ -230,9 +230,9 @@ For further details see: curvilinear meshes [DOI: 10.1007/s10915-024-02451-2](https://doi.org/10.1007/s10915-024-02451-2) """ -@inline function Trixi.flux_nonconservative_ersing_etal(u_ll, u_rr, - orientation::Integer, - equations::ShallowWaterTwoLayerEquations1D) +@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) diff --git a/src/equations/shallow_water_two_layer_2d.jl b/src/equations/shallow_water_two_layer_2d.jl index 7b8fd6f..6f4ccfe 100644 --- a/src/equations/shallow_water_two_layer_2d.jl +++ b/src/equations/shallow_water_two_layer_2d.jl @@ -343,9 +343,9 @@ For further details see: curvilinear meshes [DOI: 10.1007/s10915-024-02451-2](https://doi.org/10.1007/s10915-024-02451-2) """ -@inline function Trixi.flux_nonconservative_ersing_etal(u_ll, u_rr, - orientation::Integer, - equations::ShallowWaterTwoLayerEquations2D) +@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) @@ -381,10 +381,10 @@ For further details see: return f end -@inline function Trixi.flux_nonconservative_ersing_etal(u_ll, u_rr, - normal_direction_ll::AbstractVector, - normal_direction_average::AbstractVector, - equations::ShallowWaterTwoLayerEquations2D) +@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) diff --git a/src/equations/shallow_water_wet_dry_1d.jl b/src/equations/shallow_water_wet_dry_1d.jl index 8cca15f..b9ad94a 100644 --- a/src/equations/shallow_water_wet_dry_1d.jl +++ b/src/equations/shallow_water_wet_dry_1d.jl @@ -214,11 +214,18 @@ end Non-symmetric two-point volume flux discretizing the nonconservative (source) term that contains the gradient of the bottom topography [`ShallowWaterEquationsWetDry1D`](@ref). -Further details are available in the paper:#include("numerical_fluxes.jl") +Gives entropy conservation and well-balancedness on both the volume and surface when combined with +[`flux_wintermeyer_etal`](@ref). + +Further details are available in the papers: - 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) +- Patrick Ersing, Andrew R. Winters (2023) + An entropy stable discontinuous Galerkin method for the two-layer shallow water equations on + curvilinear meshes + [DOI: 10.48550/arXiv.2306.12699](https://doi.org/10.48550/arXiv.2306.12699) """ @inline function Trixi.flux_nonconservative_wintermeyer_etal(u_ll, u_rr, orientation::Integer, @@ -234,11 +241,8 @@ end Non-symmetric two-point surface flux discretizing the nonconservative (source) term of that contains the gradient of the bottom topography [`ShallowWaterEquationsWetDry1D`](@ref). -This contains additional terms compared to [`flux_nonconservative_wintermeyer_etal`](@ref) -that account for possible discontinuities in the bottom topography function. -Thus, this flux should be used in general at interfaces. For flux differencing volume terms, -[`flux_nonconservative_wintermeyer_etal`](@ref) is analytically equivalent but slightly -cheaper. +This flux can be used together with [`flux_fjordholm_etal`](@ref) at interfaces to ensure entropy +conservation and well-balancedness. Further details for the original finite volume formulation are available in - Ulrik S. Fjordholm, Siddhartha Mishr and Eitan Tadmor (2011) @@ -320,40 +324,12 @@ Further details on the hydrostatic reconstruction and its motivation can be foun 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), + -equations.gravity * (h_ll_star + h_ll) * (b_ll - b_star), z) end -""" - flux_nonconservative_ersing_etal(u_ll, u_rr, orientation::Integer, - equations::ShallowWaterEquationsWetDry1D) - -Non-symmetric path-conservative two-point volume flux discretizing the nonconservative (source) term -that contains the gradient of the bottom topography [`ShallowWaterEquationsWetDry1D`](@ref). - -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.1007/s10915-024-02451-2](https://doi.org/10.1007/s10915-024-02451-2) -""" -@inline function Trixi.flux_nonconservative_ersing_etal(u_ll, u_rr, - orientation::Integer, - equations::ShallowWaterEquationsWetDry1D) - return Trixi.flux_nonconservative_ersing_etal(u_ll, u_rr, orientation, - equations.basic_swe) -end - """ flux_fjordholm_etal(u_ll, u_rr, orientation, equations::ShallowWaterEquationsWetDry1D) @@ -379,7 +355,8 @@ end Total energy conservative (mathematical entropy for shallow water equations) split form. When the bottom topography is nonzero this scheme will be well-balanced when used as a `volume_flux`. -The `surface_flux` should still use, e.g., [`flux_fjordholm_etal`](@ref). +For the `surface_flux` either [`flux_wintermeyer_etal`](@ref) or [`flux_fjordholm_etal`](@ref) can +be used to ensure well-balancedness and entropy conservation. Further details are available in Theorem 1 of the paper: - Niklas Wintermeyer, Andrew R. Winters, Gregor J. Gassner and David A. Kopriva (2017) diff --git a/src/equations/shallow_water_wet_dry_2d.jl b/src/equations/shallow_water_wet_dry_2d.jl index 363fceb..76fd761 100644 --- a/src/equations/shallow_water_wet_dry_2d.jl +++ b/src/equations/shallow_water_wet_dry_2d.jl @@ -250,16 +250,18 @@ end Non-symmetric two-point volume flux discretizing the nonconservative (source) term that contains the gradient of the bottom topography [`ShallowWaterEquationsWetDry2D`](@ref). -On curvilinear meshes, this nonconservative flux depends on both the -contravariant vector (normal direction) at the current node and the averaged -one. This is different from numerical fluxes used to discretize conservative -terms. +For the `surface_flux` either [`flux_wintermeyer_etal`](@ref) or [`flux_fjordholm_etal`](@ref) can +be used to ensure well-balancedness and entropy conservation. -Further details are available in the paper: +Further details are available in the papers: - 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) +- Patrick Ersing, Andrew R. Winters (2023) + An entropy stable discontinuous Galerkin method for the two-layer shallow water equations on + curvilinear meshes + [DOI: 10.48550/arXiv.2306.12699](https://doi.org/10.48550/arXiv.2306.12699) """ @inline function Trixi.flux_nonconservative_wintermeyer_etal(u_ll, u_rr, orientation::Integer, @@ -289,16 +291,8 @@ end Non-symmetric two-point surface flux discretizing the nonconservative (source) term of that contains the gradient of the bottom topography [`ShallowWaterEquationsWetDry2D`](@ref). -On curvilinear meshes, this nonconservative flux depends on both the -contravariant vector (normal direction) at the current node and the averaged -one. This is different from numerical fluxes used to discretize conservative -terms. - -This contains additional terms compared to [`flux_nonconservative_wintermeyer_etal`](@ref) -that account for possible discontinuities in the bottom topography function. -Thus, this flux should be used in general at interfaces. For flux differencing volume terms, -[`flux_nonconservative_wintermeyer_etal`](@ref) is analytically equivalent but slightly -cheaper. +This flux can be used together with [`flux_fjordholm_etal`](@ref) at interfaces to ensure entropy +conservation and well-balancedness. Further details for the original finite volume formulation are available in - Ulrik S. Fjordholm, Siddhartha Mishr and Eitan Tadmor (2011) @@ -487,18 +481,15 @@ Further details on the hydrostatic reconstruction and its motivation can be foun 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), + -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), + -g * (h_ll_star + h_ll) * (b_ll - b_star), z) end @@ -524,20 +515,10 @@ end # 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) + f2 = -normal_direction_average[1] * equations.gravity * (h_ll_star + h_ll) * + (b_ll - b_star) + f3 = -normal_direction_average[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)) @@ -545,52 +526,6 @@ end return SVector(f1, f2, f3, f4) end -""" - flux_nonconservative_ersing_etal(u_ll, u_rr, orientation::Integer, - equations::ShallowWaterEquationsWetDry2D) - flux_nonconservative_ersing_etal(u_ll, u_rr, - normal_direction_ll::AbstractVector, - normal_direction_average::AbstractVector, - equations::ShallowWaterEquationsWetDry2D) - -!!! warning "Experimental code" - This numerical flux is experimental and may change in any future release. - -Non-symmetric path-conservative two-point volume flux discretizing the nonconservative (source) term -that contains the gradient of the bottom topography [`ShallowWaterEquationsWetDry2D`](@ref). - -On curvilinear meshes, this nonconservative flux depends on both the -contravariant vector (normal direction) at the current node and the averaged -one. This is different from numerical fluxes used to discretize conservative -terms. - -This is a modified version of [`flux_nonconservative_wintermeyer_etal`](@ref) that gives entropy -conservation and well-balancedness in both the volume and surface when combined with -[`flux_wintermeyer_etal`](@ref). - -For further details see: -- Patrick Ersing, Andrew R. Winters (2023) - An entropy stable discontinuous Galerkin method for the two-layer shallow water equations on - curvilinear meshes - [DOI: 10.1007/s10915-024-02451-2](https://doi.org/10.1007/s10915-024-02451-2) -""" -@inline function Trixi.flux_nonconservative_ersing_etal(u_ll, u_rr, - orientation::Integer, - equations::ShallowWaterEquationsWetDry2D) - return Trixi.flux_nonconservative_ersing_etal(u_ll, u_rr, orientation, - equations.basic_swe) -end - -@inline function Trixi.flux_nonconservative_ersing_etal(u_ll, u_rr, - normal_direction_ll::AbstractVector, - normal_direction_average::AbstractVector, - equations::ShallowWaterEquationsWetDry2D) - Trixi.flux_nonconservative_ersing_etal(u_ll, u_rr, - normal_direction_ll, - normal_direction_average, - equations.basic_swe) -end - """ flux_fjordholm_etal(u_ll, u_rr, orientation_or_normal_direction, equations::ShallowWaterEquationsWetDry2D) @@ -622,7 +557,8 @@ end Total energy conservative (mathematical entropy for shallow water equations) split form. When the bottom topography is nonzero this scheme will be well-balanced when used as a `volume_flux`. -The `surface_flux` should still use, e.g., [`flux_fjordholm_etal`](@ref). +For the `surface_flux` either [`flux_wintermeyer_etal`](@ref) or [`flux_fjordholm_etal`](@ref) can +be used to ensure well-balancedness and entropy conservation. Further details are available in Theorem 1 of the paper: - Niklas Wintermeyer, Andrew R. Winters, Gregor J. Gassner and David A. Kopriva (2017) diff --git a/test/Project.toml b/test/Project.toml index a3f6ad8..ff4b5ba 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -12,7 +12,7 @@ OrdinaryDiffEq = "6.49.1" Printf = "1" Roots = "2.1.6" Test = "1" -Trixi = "0.7 - 0.8.6" +Trixi = "0.8.7" [preferences.OrdinaryDiffEq] PrecompileAutoSpecialize = false diff --git a/test/test_tree_1d.jl b/test/test_tree_1d.jl index 64752ac..1d38aac 100644 --- a/test/test_tree_1d.jl +++ b/test/test_tree_1d.jl @@ -106,7 +106,7 @@ isdir(outdir) && rm(outdir, recursive = true) end end - @trixi_testset "elixir_shallowwater_well_balanced.jl with flux_nonconservative_ersing_etal" begin + @trixi_testset "elixir_shallowwater_well_balanced.jl with flux_nonconservative_wintermeyer_etal" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_well_balanced.jl"), l2=[ @@ -116,9 +116,7 @@ isdir(outdir) && rm(outdir, recursive = true) ], linf=[2.0000000000000004, 3.0610625110157164e-14, 2.0], surface_flux=(flux_wintermeyer_etal, - flux_nonconservative_ersing_etal), - volume_flux=(flux_wintermeyer_etal, - flux_nonconservative_ersing_etal), + flux_nonconservative_wintermeyer_etal), tspan=(0.0, 0.25)) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) @@ -206,7 +204,7 @@ isdir(outdir) && rm(outdir, recursive = true) end end - @trixi_testset "elixir_shallowwater_source_terms.jl with flux_nonconservative_ersing_etal" begin + @trixi_testset "elixir_shallowwater_source_terms.jl with flux_nonconservative_wintermeyer_etal" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"), l2=[ @@ -220,9 +218,7 @@ isdir(outdir) && rm(outdir, recursive = true) 9.098379777450205e-5, ], surface_flux=(flux_wintermeyer_etal, - flux_nonconservative_ersing_etal), - volume_flux=(flux_wintermeyer_etal, - flux_nonconservative_ersing_etal), + flux_nonconservative_wintermeyer_etal), tspan=(0.0, 0.025)) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) @@ -238,13 +234,13 @@ isdir(outdir) && rm(outdir, recursive = true) @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms_dirichlet.jl"), l2=[ - 0.0022851099219788917, - 0.01560453773635554, - 4.43649172558535e-5, + 0.0022667320585353927, + 0.01571629729279524, + 4.4364917255842716e-5, ], linf=[ - 0.008934615705174398, - 0.059403169140869405, + 0.008945234652224965, + 0.059403165802872415, 9.098379777405796e-5, ], tspan=(0.0, 0.025)) @@ -262,13 +258,13 @@ isdir(outdir) && rm(outdir, recursive = true) @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms_dirichlet.jl"), l2=[ - 0.0022956052733432287, - 0.015540053559855601, - 4.43649172558535e-5, + 0.0022774071143995952, + 0.01566214422689219, + 4.4364917255842716e-5, ], linf=[ - 0.008460440313118323, - 0.05720939349382359, + 0.008451721489057373, + 0.05720939055279217, 9.098379777405796e-5, ], surface_flux=(FluxHydrostaticReconstruction(FluxHLL(min_max_speed_naive), diff --git a/test/test_tree_2d.jl b/test/test_tree_2d.jl index c52c6ea..1c6404d 100644 --- a/test/test_tree_2d.jl +++ b/test/test_tree_2d.jl @@ -122,7 +122,7 @@ isdir(outdir) && rm(outdir, recursive = true) end end - @trixi_testset "elixir_shallowwater_well_balanced.jl with flux_nonconservative_ersing_etal" begin + @trixi_testset "elixir_shallowwater_well_balanced.jl with flux_nonconservative_wintermeyer_etal" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_well_balanced.jl"), l2=[ @@ -138,9 +138,7 @@ isdir(outdir) && rm(outdir, recursive = true) 2.1130620376156584, ], surface_flux=(flux_wintermeyer_etal, - flux_nonconservative_ersing_etal), - volume_flux=(flux_wintermeyer_etal, - flux_nonconservative_ersing_etal), + flux_nonconservative_wintermeyer_etal), tspan=(0.0, 0.25)) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) @@ -208,15 +206,15 @@ isdir(outdir) && rm(outdir, recursive = true) @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms_dirichlet.jl"), l2=[ - 0.0018746929418489125, - 0.017332321628469628, - 0.01634953679145536, - 6.274146767717023e-5, + 0.0018596727473552813, + 0.017306217777629147, + 0.016367646997420396, + 6.274146767723934e-5, ], linf=[ - 0.016262353691956388, - 0.08726160620859424, - 0.09043621801418844, + 0.016548007102923368, + 0.08726160568822783, + 0.09043621622245013, 0.0001819675955490041, ], tspan=(0.0, 0.025)) @@ -258,7 +256,7 @@ isdir(outdir) && rm(outdir, recursive = true) end end - @trixi_testset "elixir_shallowwater_source_terms.jl with flux_nonconservative_ersing_etal" begin + @trixi_testset "elixir_shallowwater_source_terms.jl with flux_nonconservative_wintermeyer_etal" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"), l2=[ @@ -274,9 +272,7 @@ isdir(outdir) && rm(outdir, recursive = true) 0.0001819675955490041, ], surface_flux=(flux_wintermeyer_etal, - flux_nonconservative_ersing_etal), - volume_flux=(flux_wintermeyer_etal, - flux_nonconservative_ersing_etal), + flux_nonconservative_wintermeyer_etal), tspan=(0.0, 0.25)) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) diff --git a/test/test_unstructured_2d.jl b/test/test_unstructured_2d.jl index 03e36f5..2c14a5d 100644 --- a/test/test_unstructured_2d.jl +++ b/test/test_unstructured_2d.jl @@ -98,7 +98,7 @@ isdir(outdir) && rm(outdir, recursive = true) end end - @trixi_testset "elixir_shallowwater_well_balanced.jl with flux_nonconservative_ersing_etal" begin + @trixi_testset "elixir_shallowwater_well_balanced.jl with flux_nonconservative_wintermeyer_etal" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_well_balanced.jl"), l2=[ @@ -114,9 +114,7 @@ isdir(outdir) && rm(outdir, recursive = true) 1.513851228231574, ], surface_flux=(flux_wintermeyer_etal, - flux_nonconservative_ersing_etal), - volume_flux=(flux_wintermeyer_etal, - flux_nonconservative_ersing_etal), + flux_nonconservative_wintermeyer_etal), tspan=(0.0, 0.25), atol=1e-10) # Ensure that we do not have excessive memory allocations @@ -184,7 +182,7 @@ isdir(outdir) && rm(outdir, recursive = true) end end - @trixi_testset "elixir_shallowwater_source_terms.jl with flux_nonconservative_ersing_etal" begin + @trixi_testset "elixir_shallowwater_source_terms.jl with flux_nonconservative_wintermeyer_etal" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"), l2=[ @@ -200,9 +198,7 @@ isdir(outdir) && rm(outdir, recursive = true) 2.6407324614341476e-5, ], surface_flux=(flux_wintermeyer_etal, - flux_nonconservative_ersing_etal), - volume_flux=(flux_wintermeyer_etal, - flux_nonconservative_ersing_etal), + flux_nonconservative_wintermeyer_etal), tspan=(0.0, 0.025)) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) @@ -245,15 +241,15 @@ isdir(outdir) && rm(outdir, recursive = true) @trixi_testset "elixir_shallowwater_dirichlet.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_dirichlet.jl"), l2=[ - 1.1577518608940115e-5, - 4.867189932537344e-13, - 4.647273240470541e-13, - 1.1577518608933468e-5, + 1.1577518608928104e-5, + 4.761468345949485e-13, + 4.546054642605431e-13, + 1.157751860893347e-5, ], linf=[ - 8.394063878602864e-5, - 1.1469760027632646e-10, - 1.1146619484429974e-10, + 8.394063878847113e-5, + 1.1211499950000422e-10, + 1.0890394254534975e-10, 8.394063879602065e-5, ], tspan=(0.0, 2.0),