From 35126fbfcaa237fb9bdf1ae7dbdadc5e6183b7db Mon Sep 17 00:00:00 2001 From: Patrick Ersing <114223904+patrickersing@users.noreply.github.com> Date: Tue, 2 Apr 2024 15:52:06 +0200 Subject: [PATCH] [WIP] Add 1D multilayer shallow water equations (#30) * add first implementation of the ml-swe in 1D * fix typo and formatting * add well-balanced test and lake_at_rest_error * add additional analysis functions and unit tests * add wall_bc, specialized dissipation and tests * switch to three-layer convergence test * update reference values for tests * Adjust comments * apply formatter * apply formatter * add changes from code review * fix comment * add comment for energy_total * add unit test for initial_condition_convergence --- ...xir_shallowwater_multilayer_convergence.jl | 61 ++ ...ir_shallowwater_multilayer_dam_break_ec.jl | 92 +++ ...ir_shallowwater_multilayer_dam_break_es.jl | 94 +++ ...r_shallowwater_multilayer_well_balanced.jl | 86 +++ src/TrixiShallowWater.jl | 12 +- src/equations/equations.jl | 28 + src/equations/shallow_water_multilayer_1d.jl | 554 ++++++++++++++++++ src/equations/shallow_water_two_layer_1d.jl | 4 +- src/equations/shallow_water_two_layer_2d.jl | 4 +- src/equations/shallow_water_wet_dry_1d.jl | 2 +- src/equations/shallow_water_wet_dry_2d.jl | 2 +- test/test_tree_1d.jl | 164 ++++++ test/test_unit.jl | 68 ++- 13 files changed, 1154 insertions(+), 17 deletions(-) create mode 100644 examples/tree_1d_dgsem/elixir_shallowwater_multilayer_convergence.jl create mode 100644 examples/tree_1d_dgsem/elixir_shallowwater_multilayer_dam_break_ec.jl create mode 100644 examples/tree_1d_dgsem/elixir_shallowwater_multilayer_dam_break_es.jl create mode 100644 examples/tree_1d_dgsem/elixir_shallowwater_multilayer_well_balanced.jl create mode 100644 src/equations/shallow_water_multilayer_1d.jl diff --git a/examples/tree_1d_dgsem/elixir_shallowwater_multilayer_convergence.jl b/examples/tree_1d_dgsem/elixir_shallowwater_multilayer_convergence.jl new file mode 100644 index 0000000..8a9384f --- /dev/null +++ b/examples/tree_1d_dgsem/elixir_shallowwater_multilayer_convergence.jl @@ -0,0 +1,61 @@ + +using OrdinaryDiffEq +using Trixi +using TrixiShallowWater + +############################################################################### +# Semidiscretization of the multilayer shallow water equations with three layers + +equations = ShallowWaterMultiLayerEquations1D(gravity_constant = 10.0, + rhos = (0.9, 1.0, 1.1)) + +initial_condition = initial_condition_convergence_test + +############################################################################### +# Get the DG approximation space + +volume_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal) +solver = DGSEM(polydeg = 3, + surface_flux = (flux_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 = sqrt(2.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, + 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_multilayer_dam_break_ec.jl b/examples/tree_1d_dgsem/elixir_shallowwater_multilayer_dam_break_ec.jl new file mode 100644 index 0000000..7757f13 --- /dev/null +++ b/examples/tree_1d_dgsem/elixir_shallowwater_multilayer_dam_break_ec.jl @@ -0,0 +1,92 @@ + +using OrdinaryDiffEq +using Trixi +using TrixiShallowWater + +############################################################################### +# Semidiscretization of the multilayer shallow water equations for a dam break +# test with a discontinuous bottom topography function to test entropy conservation + +equations = ShallowWaterMultiLayerEquations1D(gravity_constant = 9.81, + rhos = (0.85, 0.9, 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::ShallowWaterMultiLayerEquations1D) + v = [0.0, 0.0, 0.0] + + # Set the discontinuity + if x[1] <= 10.0 + H = [6.0, 4.0, 2.0] + b = 0.0 + else + H = [5.5, 3.5, 1.5] + b = 0.5 + end + + return prim2cons(SVector(H..., v..., b), equations) +end + +initial_condition = initial_condition_dam_break + +############################################################################### +# Get the DG approximation space + +volume_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal) +solver = DGSEM(polydeg = 3, + surface_flux = (flux_ersing_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_multilayer_dam_break_es.jl b/examples/tree_1d_dgsem/elixir_shallowwater_multilayer_dam_break_es.jl new file mode 100644 index 0000000..5f7c202 --- /dev/null +++ b/examples/tree_1d_dgsem/elixir_shallowwater_multilayer_dam_break_es.jl @@ -0,0 +1,94 @@ + +using OrdinaryDiffEq +using Trixi +using TrixiShallowWater + +############################################################################### +# Semidiscretization of the multilayer shallow water equations for a dam break +# test with a discontinuous bottom topography function for an entropy stable flux + +equations = ShallowWaterMultiLayerEquations1D(gravity_constant = 9.81, + rhos = (0.85, 0.9, 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::ShallowWaterMultiLayerEquations1D) + v = [0.0, 0.0, 0.0] + + # Set the discontinuity + if x[1] <= 10.0 + H = [6.0, 4.0, 2.0] + b = 0.0 + else + H = [5.5, 3.5, 1.5] + b = 0.5 + end + + return prim2cons(SVector(H..., v..., b), equations) +end + +initial_condition = initial_condition_dam_break + +############################################################################### +# Get the DG approximation space + +volume_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal) +solver = DGSEM(polydeg = 3, + surface_flux = (FluxPlusDissipation(flux_ersing_etal, + DissipationLocalLaxFriedrichs()), + 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_multilayer_well_balanced.jl b/examples/tree_1d_dgsem/elixir_shallowwater_multilayer_well_balanced.jl new file mode 100644 index 0000000..4bc9c63 --- /dev/null +++ b/examples/tree_1d_dgsem/elixir_shallowwater_multilayer_well_balanced.jl @@ -0,0 +1,86 @@ + +using OrdinaryDiffEq +using Trixi +using TrixiShallowWater + +############################################################################### +# Semidiscretization of the multilayer shallow water equations to test well-balancedness + +equations = ShallowWaterMultiLayerEquations1D(gravity_constant = 1.0, H0 = 0.7, + rhos = (0.8, 0.9, 1.0)) + +""" + initial_condition_fjordholm_well_balanced(x, t, equations::ShallowWaterMultiLayerEquations1D) + +Initial condition to test well balanced with a bottom topography adapted 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::ShallowWaterMultiLayerEquations1D) + inicenter = 0.5 + x_norm = x[1] - inicenter + r = abs(x_norm) + + H = [0.7, 0.6, 0.5] + v = [0.0, 0.0, 0.0] + b = r <= 0.1 ? 0.2 * (cos(10 * pi * (x[1] - 0.5)) + 1) : 0.0 + + return prim2cons(SVector(H..., v..., b), equations) +end + +initial_condition = initial_condition_fjordholm_well_balanced + +############################################################################### +# Get the DG approximation space + +volume_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal) +solver = DGSEM(polydeg = 3, + surface_flux = (flux_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/src/TrixiShallowWater.jl b/src/TrixiShallowWater.jl index 3932671..ac74b97 100644 --- a/src/TrixiShallowWater.jl +++ b/src/TrixiShallowWater.jl @@ -8,7 +8,7 @@ using Trixi # Import additional symbols that are not exported by Trixi.jl using Trixi: get_node_vars, set_node_vars!, waterheight using MuladdMacro: @muladd -using StaticArrays: SVector, @SMatrix +using StaticArrays: SVector, @SMatrix, MVector using Static: True, False using LinearAlgebra: norm @@ -19,15 +19,17 @@ include("solvers/indicators.jl") # Export types/functions that define the public API of TrixiShallowWater.jl export ShallowWaterEquationsWetDry1D, ShallowWaterEquationsWetDry2D, - ShallowWaterTwoLayerEquations1D, ShallowWaterTwoLayerEquations2D + ShallowWaterTwoLayerEquations1D, ShallowWaterTwoLayerEquations2D, + ShallowWaterMultiLayerEquations1D export hydrostatic_reconstruction_chen_noelle, flux_nonconservative_chen_noelle, - min_max_speed_chen_noelle, - flux_hll_chen_noelle + min_max_speed_chen_noelle, flux_hll_chen_noelle, + flux_ersing_etal, flux_es_ersing_etal -export flux_es_ersing_etal +export nlayers, eachlayer export PositivityPreservingLimiterShallowWater + export IndicatorHennemannGassnerShallowWater end diff --git a/src/equations/equations.jl b/src/equations/equations.jl index 076a994..a36f993 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -13,4 +13,32 @@ include("shallow_water_wet_dry_2d.jl") include("shallow_water_two_layer_1d.jl") include("shallow_water_two_layer_2d.jl") + +abstract type AbstractShallowWaterMultiLayerEquations{NDIMS, NVARS, NLAYERS} <: + Trixi.AbstractEquations{NDIMS, NVARS} end +include("shallow_water_multilayer_1d.jl") + +""" + eachlayer(equations::AbstractShallowWaterMultiLayerEquations) + +Return an iterator over the indices that specify the location in relevant data structures +for the layers in `AbstractShallowWaterMultiLayerEquations`. +""" +@inline function eachlayer(equations::AbstractShallowWaterMultiLayerEquations) + Base.OneTo(nlayers(equations)) +end + +""" + nlayers(equations::AbstractShallowWaterMultiLayerEquations) + +Retrieve the number of layers from an equation instance of the `AbstractShallowWaterMultiLayerEquations`. +""" +@inline function nlayers(::AbstractShallowWaterMultiLayerEquations{NDIMS, NVARS, + NLAYERS}) where { + NDIMS, + NVARS, + NLAYERS + } + NLAYERS +end end # @muladd diff --git a/src/equations/shallow_water_multilayer_1d.jl b/src/equations/shallow_water_multilayer_1d.jl new file mode 100644 index 0000000..0dec78a --- /dev/null +++ b/src/equations/shallow_water_multilayer_1d.jl @@ -0,0 +1,554 @@ +# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). +# Since these FMAs can increase the performance of many numerical algorithms, +# we need to opt-in explicitly. +# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. +@muladd begin +#! format: noindent + +@doc raw""" + ShallowWaterMultiLayerEquations1D(gravity, H0, rhos) + +Multi-Layer Shallow Water equations (MLSWE) in one space dimension. The equations are given by +```math +\left\{ + \begin{aligned} + &\partial_t h_m + \partial_x h_mv_m = 0,\\ + &\partial h_mv_m + \partial_x h_mv_m^2 = -gh_m\partial_x \bigg(b + \sum\limits_{k\geq j}h_k + \sum\limits_{k\ + ISBN: 978-0-12-088759-0 +""" + +struct ShallowWaterMultiLayerEquations1D{NVARS, NLAYERS, RealT <: Real} <: + AbstractShallowWaterMultiLayerEquations{1, NVARS, NLAYERS} + gravity::RealT # gravitational constant + H0::RealT # constant "lake-at-rest" total water height + rhos::SVector{NLAYERS, RealT} # Vector of layer densities + + function ShallowWaterMultiLayerEquations1D{NVARS, NLAYERS, RealT}(gravity::RealT, + H0::RealT, + rhos::SVector{NLAYERS, + RealT}) where { + NVARS, + NLAYERS, + RealT <: + Real + } + # Ensure that layer densities are all positive and in increasing order + issorted(rhos) || + throw(ArgumentError("densities must be in increasing order (rhos[1] < rhos[2] < ... < rhos[NLAYERS])")) + min(rhos...) > 0 || throw(ArgumentError("densities must be positive")) + + new(gravity, H0, rhos) + end +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 ShallowWaterMultiLayerEquations1D(; gravity_constant, + H0 = zero(gravity_constant), rhos) + + # Promote all variables to a common type + _rhos = promote(rhos...) + RealT = promote_type(eltype(_rhos), eltype(gravity_constant), eltype(H0)) + __rhos = SVector(map(RealT, _rhos)) + + # Extract number of layers and variables + NLAYERS = length(rhos) + NVARS = 2 * NLAYERS + 1 + + return ShallowWaterMultiLayerEquations1D{NVARS, NLAYERS, RealT}(gravity_constant, + H0, __rhos) +end + +@inline function Base.real(::ShallowWaterMultiLayerEquations1D{NVARS, NLAYERS, RealT}) where { + NVARS, + NLAYERS, + RealT <: + Real + } + RealT +end + +Trixi.have_nonconservative_terms(::ShallowWaterMultiLayerEquations1D) = True() +function Trixi.varnames(::typeof(cons2cons), + equations::ShallowWaterMultiLayerEquations1D) + waterheight = ntuple(n -> "h" * string(n), Val(nlayers(equations))) + momentum = ntuple(n -> "h" * string(n) * "_v", Val(nlayers(equations))) + + return (waterheight..., momentum..., "b") +end + +# We use the interface heights, H = ∑h + b as primitive variables for easier visualization and setting initial +# conditions +function Trixi.varnames(::typeof(cons2prim), + equations::ShallowWaterMultiLayerEquations1D) + interface_height = ntuple(n -> "H" * string(n), Val(nlayers(equations))) + velocity = ntuple(n -> "v" * string(n) * "_1", Val(nlayers(equations))) + return (interface_height..., velocity..., "b") +end + +# Set initial conditions at physical location `x` for time `t` +""" + initial_condition_convergence_test(x, t, equations::ShallowWaterMultiLayerEquations1D) + +A smooth initial condition for a three-layer configuration used for convergence tests in combination with +[`source_terms_convergence_test`](@ref) (and +[`BoundaryConditionDirichlet(initial_condition_convergence_test)`](@ref) in non-periodic domains). +""" +function Trixi.initial_condition_convergence_test(x, t, + equations::ShallowWaterMultiLayerEquations1D) + # Check if the equation has been set up with three layers + if nlayers(equations) != 3 + throw(ArgumentError("This initial condition is only valid for three layers")) + end + + # Some constants are chosen such that the function is periodic on the domain [0,sqrt(2)] + ω = 2.0 * pi * sqrt(2.0) + + H1 = 4.0 + 0.1 * cos(ω * x[1] + t) + H2 = 2.0 + 0.1 * sin(ω * x[1] + t) + H3 = 1.5 + 0.1 * cos(ω * x[1] + t) + v1 = 0.9 + v2 = 1.0 + v3 = 1.1 + b = 1.0 + 0.1 * cos(2.0 * ω * x[1]) + + return prim2cons(SVector(H1, H2, H3, v1, v2, v3, b), equations) +end + +""" + source_terms_convergence_test(u, x, t, equations::ShallowWaterMultiLayerEquations1D) + +Source terms used for convergence tests with a three-layer configuration in combination with +[`initial_condition_convergence_test`](@ref) +(and [`BoundaryConditionDirichlet(initial_condition_convergence_test)`](@ref) +in non-periodic domains). +""" +@inline function Trixi.source_terms_convergence_test(u, x, t, + equations::ShallowWaterMultiLayerEquations1D) + # 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 * sin(t + x[1] * ω) - 0.1 * cos(t + x[1] * ω) + + 0.9 * (-0.1 * sin(t + x[1] * ω) * ω - 0.1 * cos(t + x[1] * ω) * ω) + du2 = 0.1 * sin(t + x[1] * ω) + 0.1 * cos(t + x[1] * ω) + + 0.1 * sin(t + x[1] * ω) * ω + + 0.1 * cos(t + x[1] * ω) * ω + du3 = -0.1 * sin(t + x[1] * ω) + + 1.1 * (-0.1 * sin(t + x[1] * ω) * ω + 0.2 * sin(2.0 * x[1] * ω) * ω) + du4 = 0.9 * (-0.1 * sin(t + x[1] * ω) - 0.1 * cos(t + x[1] * ω)) + + 0.81 * (-0.1 * sin(t + x[1] * ω) * ω - 0.1 * cos(t + x[1] * ω) * ω) + + 10.0 * (2.0 - 0.1 * sin(t + x[1] * ω) + 0.1 * cos(t + x[1] * ω)) * + (-0.1 * sin(t + x[1] * ω) * ω - 0.1 * cos(t + x[1] * ω) * ω) + + (2.0 - 0.1 * sin(t + x[1] * ω) + 0.1 * cos(t + x[1] * ω)) * + cos(t + x[1] * ω) * ω + du5 = 0.1 * sin(t + x[1] * ω) + 0.1 * cos(t + x[1] * ω) + + 0.1 * sin(t + x[1] * ω) * ω + + 0.1 * cos(t + x[1] * ω) * ω + + 10.0 * (0.5 + 0.1 * sin(t + x[1] * ω) - 0.1 * cos(t + x[1] * ω)) * + (-0.1 * sin(t + x[1] * ω) * ω + + 0.9 * (-0.1 * sin(t + x[1] * ω) * ω - 0.1 * cos(t + x[1] * ω) * ω)) + + 10.0 * (0.5 + 0.1 * sin(t + x[1] * ω) - 0.1 * cos(t + x[1] * ω)) * + (0.1 * sin(t + x[1] * ω) * ω + 0.1 * cos(t + x[1] * ω) * ω) + du6 = -0.11 * sin(t + x[1] * ω) + + 1.21 * (-0.1 * sin(t + x[1] * ω) * ω + 0.2 * sin(2.0x[1] * ω) * ω) + + 10.0 * (0.5 - 0.1 * cos(2.0 * x[1] * ω) + 0.1 * cos(t + x[1] * ω)) * + (-0.1 * sin(t + x[1] * ω) * ω + 0.2 * sin(2.0 * x[1] * ω) * ω) + + 10.0 * (0.5 - 0.1 * cos(2.0 * x[1] * ω) + 0.1 * cos(t + x[1] * ω)) * + (0.8181818181818181 * + (-0.1 * sin(t + x[1] * ω) * ω - 0.1 * cos(t + x[1] * ω) * ω) + + 0.9090909090909091 * + (0.1 * sin(t + x[1] * ω) * ω + 0.1 * cos(t + x[1] * ω) * ω) - + 0.2 * sin(2.0 * x[1] * ω) * ω) + + return SVector(du1, du2, du3, du4, du5, du6, zero(eltype(u))) +end + +""" + boundary_condition_slip_wall(u_inner, orientation_or_normal, x, t, surface_flux_function, + equations::ShallowWaterMultiLayerEquations1D) + +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::ShallowWaterMultiLayerEquations1D) + # Create the "external" boundary solution state + h = waterheight(u_inner, equations) + hv = momentum(u_inner, equations) + b = u_inner[end] + + u_boundary = SVector(h..., -hv..., b) + + # 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 advective portion of the flux for a single point +@inline function Trixi.flux(u, orientation::Integer, + equations::ShallowWaterMultiLayerEquations1D) + # Extract waterheight and momentum and compute velocities + hv = momentum(u, equations) + v = velocity(u, equations) + + # Initialize flux vector + f = zero(MVector{2 * nlayers(equations) + 1, real(equations)}) + # Calculate fluxes in each layer + for i in eachlayer(equations) + f_h = hv[i] + # The momentum flux simplifies as the pressure is included in the nonconservative term + f_hv = hv[i] * v[i] + + setlayer!(f, f_h, f_hv, i, equations) + end + return SVector(f) +end + +""" + flux_nonconservative_ersing_etal(u_ll, u_rr, orientation::Integer, + equations::ShallowWaterTwoLayerEquations1D) + +Non-symmetric path-conservative two-point flux discretizing the nonconservative (source) term +that contains the gradients of the bottom topography and waterheights from the coupling between layers +and the nonconservative pressure formulation [`ShallowWaterMultiLayerEquations1D`](@ref). + +When the bottom topography is nonzero this scheme will be well-balanced when used with [`flux_ersing_etal`](@ref). + +In the two-layer setting this combination is equivalent to the fluxes in: +- 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::ShallowWaterMultiLayerEquations1D) + # Pull the necessary left and right state information + h_ll = waterheight(u_ll, equations) + h_rr = waterheight(u_rr, equations) + b_rr = u_rr[end] + b_ll = u_ll[end] + + # Compute the jumps + h_jump = h_rr - h_ll + b_jump = b_rr - b_ll + g = equations.gravity + + # Initialize flux vector + f = zero(MVector{2 * nlayers(equations) + 1, real(equations)}) + + # Compute the nonconservative flux in each layer (0, ..., 0, f_hv[1], ..., f_hv[NLAYERS], 0), + # where f_hv[i] = g * h[i] * (b + ∑h[k] + ∑σ[k] * h[k])_x and σ[k] = ρ[k] / ρ[i] denotes the + # density ratio of different layers + for i in eachlayer(equations) + f_hv = g * h_ll[i] * b_jump + for j in eachlayer(equations) + if j < i + f_hv += g * h_ll[i] * + (equations.rhos[j] / equations.rhos[i] * h_jump[j]) + else # (i