From 3205c8fed7acdb190f76d333fc07723a040cbecc Mon Sep 17 00:00:00 2001 From: patrickersing Date: Thu, 11 Apr 2024 12:02:04 +0200 Subject: [PATCH] add multilayer SWE in 2D --- ...xir_shallowwater_multilayer_convergence.jl | 61 ++ ...lixir_shallowwater_multilayer_dam_break.jl | 92 ++ ...r_shallowwater_multilayer_well_balanced.jl | 79 ++ ...xir_shallowwater_multilayer_convergence.jl | 62 ++ ...lixir_shallowwater_multilayer_dam_break.jl | 97 ++ ...r_shallowwater_multilayer_well_balanced.jl | 77 ++ src/TrixiShallowWater.jl | 2 +- src/equations/equations.jl | 1 + src/equations/shallow_water_multilayer_2d.jl | 829 ++++++++++++++++++ test/test_tree_2d.jl | 248 ++++++ test/test_unit.jl | 34 +- test/test_unstructured_2d.jl | 160 ++++ 12 files changed, 1739 insertions(+), 3 deletions(-) create mode 100644 examples/tree_2d_dgsem/elixir_shallowwater_multilayer_convergence.jl create mode 100644 examples/tree_2d_dgsem/elixir_shallowwater_multilayer_dam_break.jl create mode 100644 examples/tree_2d_dgsem/elixir_shallowwater_multilayer_well_balanced.jl create mode 100644 examples/unstructured_2d_dgsem/elixir_shallowwater_multilayer_convergence.jl create mode 100644 examples/unstructured_2d_dgsem/elixir_shallowwater_multilayer_dam_break.jl create mode 100644 examples/unstructured_2d_dgsem/elixir_shallowwater_multilayer_well_balanced.jl create mode 100644 src/equations/shallow_water_multilayer_2d.jl diff --git a/examples/tree_2d_dgsem/elixir_shallowwater_multilayer_convergence.jl b/examples/tree_2d_dgsem/elixir_shallowwater_multilayer_convergence.jl new file mode 100644 index 0000000..391eef6 --- /dev/null +++ b/examples/tree_2d_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 = ShallowWaterMultiLayerEquations2D(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, 0.0) +coordinates_max = (sqrt(2.0), 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_2d_dgsem/elixir_shallowwater_multilayer_dam_break.jl b/examples/tree_2d_dgsem/elixir_shallowwater_multilayer_dam_break.jl new file mode 100644 index 0000000..2dbb108 --- /dev/null +++ b/examples/tree_2d_dgsem/elixir_shallowwater_multilayer_dam_break.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 = ShallowWaterMultiLayerEquations2D(gravity_constant = 1.0, + rhos = (0.9, 0.95, 1.0)) + +# This academic testcase sets up a discontinuous bottom topography +# function and initial condition to test entropy conservation. + +function initial_condition_dam_break(x, t, equations::ShallowWaterMultiLayerEquations2D) + # Bottom topography + b = 0.3 * exp(-0.5 * ((x[1])^2 + (x[2])^2)) + + if x[1] < 0.0 + H = SVector(1.0, 0.8, 0.6) + else + H = SVector(0.9, 0.7, 0.5) + b += 0.1 + end + + v1 = zero(H) + v2 = zero(H) + return prim2cons(SVector(H..., v1..., v2..., b), + equations) +end + +initial_condition = initial_condition_dam_break + +boundary_condition_constant = BoundaryConditionDirichlet(initial_condition_dam_break) + +############################################################################### +# Get the DG approximation space + +volume_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal) +surface_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal) +solver = DGSEM(polydeg = 3, surface_flux = surface_flux, + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +############################################################################### +# Get the TreeMesh and setup a periodic mesh + +coordinates_min = (-1.0, -1.0) +coordinates_max = (1.0, 1.0) +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 4, + n_cells_max = 10_000, + periodicity = true) + +# Create the semi discretization object +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) +############################################################################### +# ODE solver + +tspan = (0.0, 2.0) +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)) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 500, + save_initial_solution = true, + save_final_solution = true) + +stepsize_callback = StepsizeCallback(cfl = 1.0) + +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, + stepsize_callback) + +############################################################################### +# run the simulation + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); +summary_callback() # print the timer summary diff --git a/examples/tree_2d_dgsem/elixir_shallowwater_multilayer_well_balanced.jl b/examples/tree_2d_dgsem/elixir_shallowwater_multilayer_well_balanced.jl new file mode 100644 index 0000000..1a944e1 --- /dev/null +++ b/examples/tree_2d_dgsem/elixir_shallowwater_multilayer_well_balanced.jl @@ -0,0 +1,79 @@ + +using OrdinaryDiffEq +using Trixi +using TrixiShallowWater + +############################################################################### +# Semidiscretization of the two-layer shallow water equations with a bottom topography function +# to test well-balancedness + +equations = ShallowWaterMultiLayerEquations2D(gravity_constant = 9.81, H0 = 0.6, + rhos = (0.7, 0.8, 0.9, 1.0)) + +# An initial condition with constant total water height, zero velocities and a bottom topography to +# test well-balancedness +function initial_condition_well_balanced(x, t, equations::ShallowWaterMultiLayerEquations2D) + H = SVector(0.6, 0.55, 0.5, 0.45) + v1 = zero(H) + v2 = zero(H) + b = (((x[1] - 0.5)^2 + (x[2] - 0.5)^2) < 0.04 ? + 0.2 * (cos(4 * pi * sqrt((x[1] - 0.5)^2 + (x[2] + + -0.5)^2)) + 1) : 0.0) + + return prim2cons(SVector(H..., v1..., v2..., b), + equations) +end + +initial_condition = initial_condition_well_balanced + +############################################################################### +# 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 = 3, surface_flux = surface_flux, + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +############################################################################### +# Get the TreeMesh and setup a periodic mesh + +coordinates_min = (0.0, 0.0) +coordinates_max = (1.0, 1.0) +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 3, + n_cells_max = 10_000, + periodicity = true) + +# Create the semi discretization object +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +############################################################################### +# ODE solver + +tspan = (0.0, 10.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 1000 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + extra_analysis_integrals = (lake_at_rest_error,)) + +stepsize_callback = StepsizeCallback(cfl = 1.0) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 1000, + save_initial_solution = true, + save_final_solution = true) + +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, + stepsize_callback) + +############################################################################### +# run the simulation + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); +summary_callback() # print the timer summary diff --git a/examples/unstructured_2d_dgsem/elixir_shallowwater_multilayer_convergence.jl b/examples/unstructured_2d_dgsem/elixir_shallowwater_multilayer_convergence.jl new file mode 100644 index 0000000..708e7a2 --- /dev/null +++ b/examples/unstructured_2d_dgsem/elixir_shallowwater_multilayer_convergence.jl @@ -0,0 +1,62 @@ + +using OrdinaryDiffEq +using Trixi +using TrixiShallowWater + +############################################################################### +# Semidiscretization of the two-layer shallow water equations with a periodic +# bottom topography function (set in the initial conditions) + +equations = ShallowWaterMultiLayerEquations2D(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) +surface_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal) +solver = DGSEM(polydeg = 8, surface_flux = surface_flux, + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +############################################################################### +# This setup is for the curved, split form convergence test on a periodic domain + +# Get the unstructured quad mesh from a file (downloads the file if not available locally) +mesh_file = Trixi.download("https://gist.githubusercontent.com/andrewwinters5000/8f8cd23df27fcd494553f2a89f3c1ba4/raw/85e3c8d976bbe57ca3d559d653087b0889535295/mesh_alfven_wave_with_twist_and_flip.mesh", + joinpath(@__DIR__, "mesh_alfven_wave_with_twist_and_flip.mesh")) + +mesh = UnstructuredMesh2D(mesh_file, periodicity = true) + +# Create the semidiscretization object +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + source_terms = source_terms_convergence_test) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 1.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 500 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 500, + save_initial_solution = true, + save_final_solution = true) + +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/unstructured_2d_dgsem/elixir_shallowwater_multilayer_dam_break.jl b/examples/unstructured_2d_dgsem/elixir_shallowwater_multilayer_dam_break.jl new file mode 100644 index 0000000..2545288 --- /dev/null +++ b/examples/unstructured_2d_dgsem/elixir_shallowwater_multilayer_dam_break.jl @@ -0,0 +1,97 @@ + +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 = ShallowWaterMultiLayerEquations2D(gravity_constant = 1.0, + rhos = (0.9, 0.95, 1.0)) + +# This academic testcase sets up a discontinuous bottom topography +# function and initial condition to test entropy conservation. + +function initial_condition_dam_break(x, t, equations::ShallowWaterMultiLayerEquations2D) + # Bottom topography + b = 0.3 * exp(-0.5 * ((x[1] - sqrt(2) / 2)^2 + (x[2] - sqrt(2) / 2)^2)) + + if x[1] < sqrt(2) / 2 + H = SVector(1.0, 0.8, 0.6) + else + H = SVector(0.9, 0.7, 0.5) + b += 0.1 + end + + v1 = zero(H) + v2 = zero(H) + return prim2cons(SVector(H..., v1..., v2..., b), + equations) +end + +initial_condition = initial_condition_dam_break + +boundary_condition_constant = BoundaryConditionDirichlet(initial_condition_dam_break) + +############################################################################### +# Get the DG approximation space + +volume_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal) +surface_flux = (flux_ersing_etal, flux_nonconservative_ersing_etal) +solver = DGSEM(polydeg = 6, surface_flux = surface_flux, + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +############################################################################### +# Get the unstructured quad mesh from a file (downloads the file if not available locally) +mesh_file = Trixi.download("https://gist.githubusercontent.com/andrewwinters5000/8f8cd23df27fcd494553f2a89f3c1ba4/raw/85e3c8d976bbe57ca3d559d653087b0889535295/mesh_alfven_wave_with_twist_and_flip.mesh", + joinpath(@__DIR__, "mesh_alfven_wave_with_twist_and_flip.mesh")) + +mesh = UnstructuredMesh2D(mesh_file, periodicity = false) + +# Boundary conditions +boundary_condition = Dict(:Top => boundary_condition_slip_wall, + :Left => boundary_condition_slip_wall, + :Right => boundary_condition_slip_wall, + :Bottom => boundary_condition_slip_wall) + +# Create the semi discretization object +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, + solver, boundary_conditions = boundary_condition) + +############################################################################### +# ODE solver + +tspan = (0.0, 0.5) +ode = semidiscretize(semi, tspan) + +############################################################################### +# Callbacks + +summary_callback = SummaryCallback() + +analysis_interval = 500 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + save_analysis = false, + extra_analysis_integrals = (energy_total, + energy_kinetic, + energy_internal)) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 500, + save_initial_solution = true, + save_final_solution = true) + +stepsize_callback = StepsizeCallback(cfl = 1.0) + +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, + stepsize_callback) + +############################################################################### +# run the simulation + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); +summary_callback() # print the timer summary diff --git a/examples/unstructured_2d_dgsem/elixir_shallowwater_multilayer_well_balanced.jl b/examples/unstructured_2d_dgsem/elixir_shallowwater_multilayer_well_balanced.jl new file mode 100644 index 0000000..4c4763b --- /dev/null +++ b/examples/unstructured_2d_dgsem/elixir_shallowwater_multilayer_well_balanced.jl @@ -0,0 +1,77 @@ + +using OrdinaryDiffEq +using Trixi +using TrixiShallowWater + +############################################################################### +# Semidiscretization of the two-layer shallow water equations with a bottom topography function +# to test well-balancedness + +equations = ShallowWaterMultiLayerEquations2D(gravity_constant = 9.81, H0 = 0.6, + rhos = (0.7, 0.8, 0.9, 1.0)) + +# An initial condition with constant total water height, zero velocities and a bottom topography to +# test well-balancedness +function initial_condition_well_balanced(x, t, equations::ShallowWaterMultiLayerEquations2D) + H = SVector(0.6, 0.55, 0.5, 0.45) + v1 = zero(H) + v2 = zero(H) + b = (((x[1] - 0.5)^2 + (x[2] - 0.5)^2) < 0.04 ? + 0.2 * (cos(4 * pi * sqrt((x[1] - 0.5)^2 + (x[2] + + -0.5)^2)) + 1) : 0.0) + + return prim2cons(SVector(H..., v1..., v2..., b), + equations) +end + +initial_condition = initial_condition_well_balanced + +############################################################################### +# 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 = 6, surface_flux = surface_flux, + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +############################################################################### +# This setup is for the curved, split form convergence test on a periodic domain + +# Get the unstructured quad mesh from a file (downloads the file if not available locally) +mesh_file = Trixi.download("https://gist.githubusercontent.com/andrewwinters5000/8f8cd23df27fcd494553f2a89f3c1ba4/raw/85e3c8d976bbe57ca3d559d653087b0889535295/mesh_alfven_wave_with_twist_and_flip.mesh", + joinpath(@__DIR__, "mesh_alfven_wave_with_twist_and_flip.mesh")) + +mesh = UnstructuredMesh2D(mesh_file, periodicity = true) + +# Create the semidiscretization object +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) +############################################################################### +# ODE solver + +tspan = (0.0, 10.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 1000 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + extra_analysis_integrals = (lake_at_rest_error,)) + +stepsize_callback = StepsizeCallback(cfl = 1.0) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 1000, + save_initial_solution = true, + save_final_solution = true) + +callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, + stepsize_callback) + +############################################################################### +# run the simulation + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); +summary_callback() # print the timer summary diff --git a/src/TrixiShallowWater.jl b/src/TrixiShallowWater.jl index ac74b97..c3b0147 100644 --- a/src/TrixiShallowWater.jl +++ b/src/TrixiShallowWater.jl @@ -20,7 +20,7 @@ include("solvers/indicators.jl") # Export types/functions that define the public API of TrixiShallowWater.jl export ShallowWaterEquationsWetDry1D, ShallowWaterEquationsWetDry2D, ShallowWaterTwoLayerEquations1D, ShallowWaterTwoLayerEquations2D, - ShallowWaterMultiLayerEquations1D + ShallowWaterMultiLayerEquations1D, ShallowWaterMultiLayerEquations2D export hydrostatic_reconstruction_chen_noelle, flux_nonconservative_chen_noelle, min_max_speed_chen_noelle, flux_hll_chen_noelle, diff --git a/src/equations/equations.jl b/src/equations/equations.jl index a36f993..89713b9 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -17,6 +17,7 @@ include("shallow_water_two_layer_2d.jl") abstract type AbstractShallowWaterMultiLayerEquations{NDIMS, NVARS, NLAYERS} <: Trixi.AbstractEquations{NDIMS, NVARS} end include("shallow_water_multilayer_1d.jl") +include("shallow_water_multilayer_2d.jl") """ eachlayer(equations::AbstractShallowWaterMultiLayerEquations) diff --git a/src/equations/shallow_water_multilayer_2d.jl b/src/equations/shallow_water_multilayer_2d.jl new file mode 100644 index 0000000..73dea6d --- /dev/null +++ b/src/equations/shallow_water_multilayer_2d.jl @@ -0,0 +1,829 @@ +# 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""" + ShallowWaterMultiLayerEquations2D(gravity, H0, rhos) + +Multi-Layer Shallow Water equations (MLSWE) in two space dimension. The equations are given by +```math +\left\{ + \begin{aligned} + &\partial_t h_m + \partial_x h_mv_m = 0,\\ + &\partial h_mv1_m + \partial_x h_mv1_m^2 + \partial_y h_mv1_mv2_m = -gh_m\partial_x \bigg(b + \sum\limits_{k\geq j}h_k + \sum\limits_{k\ + ISBN: 978-0-12-088759-0 +""" + +struct ShallowWaterMultiLayerEquations2D{NVARS, NLAYERS, RealT <: Real} <: + AbstractShallowWaterMultiLayerEquations{2, NVARS, NLAYERS} + gravity::RealT # gravitational constant + H0::RealT # constant "lake-at-rest" total water height + rhos::SVector{NLAYERS, RealT} # Vector of layer densities + + function ShallowWaterMultiLayerEquations2D{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 ShallowWaterMultiLayerEquations2D(; 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 = 3 * NLAYERS + 1 + + return ShallowWaterMultiLayerEquations2D{NVARS, NLAYERS, RealT}(gravity_constant, + H0, __rhos) +end + +@inline function Base.real(::ShallowWaterMultiLayerEquations2D{NVARS, NLAYERS, RealT}) where { + NVARS, + NLAYERS, + RealT <: + Real + } + RealT +end + +Trixi.have_nonconservative_terms(::ShallowWaterMultiLayerEquations2D) = True() +function Trixi.varnames(::typeof(cons2cons), + equations::ShallowWaterMultiLayerEquations2D) + waterheight = ntuple(n -> "h" * string(n), Val(nlayers(equations))) + momentum_1 = ntuple(n -> "h" * string(n) * "_v1", Val(nlayers(equations))) + momentum_2 = ntuple(n -> "h" * string(n) * "_v2", Val(nlayers(equations))) + return (waterheight..., momentum_1..., momentum_2..., "b") +end + +# We use the total layer heights, H = ∑h + b as primitive variables for easier visualization and setting initial +# conditions +function Trixi.varnames(::typeof(cons2prim), + equations::ShallowWaterMultiLayerEquations2D) + total_layer_height = ntuple(n -> "H" * string(n), Val(nlayers(equations))) + velocity_1 = ntuple(n -> "v" * string(n) * "_1", Val(nlayers(equations))) + velocity_2 = ntuple(n -> "v" * string(n) * "_2", Val(nlayers(equations))) + return (total_layer_height..., velocity_1..., velocity_2..., "b") +end + +# TODO: This needs to be redone +# # Set initial conditions at physical location `x` for time `t` +""" + initial_condition_convergence_test(x, t, equations::ShallowWaterMultiLayerEquations2D) + +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::ShallowWaterMultiLayerEquations2D) + # 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) + + H3 = 1.5 + 0.1 * cos(ω * x[1] + t) + 0.1 * cos(ω * x[2] + t) + H2 = 2.0 + 0.1 * sin(ω * x[1] + t) + 0.1 * sin(ω * x[2] + t) + H1 = 4.0 + 0.1 * cos(ω * x[1] + t) + 0.1 * cos(ω * x[2] + t) + H = (H1, H2, H3) + v1 = (0.8, 0.8, 0.8) + v2 = (1.0, 1.0, 1.0) + + b = 1.0 + 0.1 * cos(ω * x[1]) + 0.1 * cos(ω * x[2]) + + return prim2cons(SVector(H..., v1..., v2..., b), equations) +end + +""" + source_terms_convergence_test(u, x, t, equations::ShallowWaterMultiLayerEquations2D) + +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::ShallowWaterMultiLayerEquations2D) + # 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.1cos(t + x[2] * ω) - 0.1sin(t + x[1] * ω) - 0.1sin(t + x[2] * ω) - + 0.1cos(t + x[1] * ω) - 0.1cos(t + x[2] * ω) * ω + + 0.8(-0.1sin(t + x[1] * ω) * ω - 0.1cos(t + x[1] * ω) * ω) - + 0.1sin(t + x[2] * ω) * ω + + du2 = 0.1cos(t + x[2] * ω) + 0.1sin(t + x[1] * ω) + 0.1sin(t + x[2] * ω) + + 0.1cos(t + x[1] * ω) + 0.1cos(t + x[2] * ω) * ω + + 0.8(0.1sin(t + x[1] * ω) * ω + 0.1cos(t + x[1] * ω) * ω) + + 0.1sin(t + x[2] * ω) * ω + + du3 = -0.1sin(t + x[1] * ω) - 0.1sin(t + x[2] * ω) + 0.1sin(x[2] * ω) * ω + + 0.8(0.1sin(x[1] * ω) * ω - 0.1sin(t + x[1] * ω) * ω) - + 0.1sin(t + x[2] * ω) * ω + + du4 = 0.8(-0.1cos(t + x[2] * ω) - 0.1sin(t + x[1] * ω) - 0.1sin(t + x[2] * ω) - + 0.1cos(t + x[1] * ω)) + + 0.8(-0.1cos(t + x[2] * ω) * ω - 0.1sin(t + x[2] * ω) * ω) + + 0.8^2 * (-0.1sin(t + x[1] * ω) * ω - 0.1cos(t + x[1] * ω) * ω) + + 10.0(2.0 + 0.1cos(t + x[2] * ω) - 0.1sin(t + x[1] * ω) - + 0.1sin(t + x[2] * ω) + 0.1cos(t + x[1] * ω)) * + (-0.1sin(t + x[1] * ω) * ω - 0.1cos(t + x[1] * ω) * ω) + + (2.0 + 0.1cos(t + x[2] * ω) - 0.1sin(t + x[1] * ω) - 0.1sin(t + x[2] * ω) + + 0.1cos(t + x[1] * ω)) * cos(t + x[1] * ω) * ω + + du5 = 0.8(0.1cos(t + x[2] * ω) + 0.1sin(t + x[1] * ω) + 0.1sin(t + x[2] * ω) + + 0.1cos(t + x[1] * ω)) + + 0.8(0.1cos(t + x[2] * ω) * ω + 0.1sin(t + x[2] * ω) * ω) + + 0.8^2 * (0.1sin(t + x[1] * ω) * ω + 0.1cos(t + x[1] * ω) * ω) + + 10.0(0.5 - 0.1cos(t + x[2] * ω) + 0.1sin(t + x[1] * ω) + + 0.1sin(t + x[2] * ω) - 0.1cos(t + x[1] * ω)) * + (0.1sin(t + x[1] * ω) * ω + 0.1cos(t + x[1] * ω) * ω) + + 10.0(0.5 - 0.1cos(t + x[2] * ω) + 0.1sin(t + x[1] * ω) + + 0.1sin(t + x[2] * ω) - 0.1cos(t + x[1] * ω)) * + (-0.1sin(t + x[1] * ω) * ω + + 0.9(-0.1sin(t + x[1] * ω) * ω - 0.1cos(t + x[1] * ω) * ω)) + + du6 = 0.8(-0.1sin(t + x[1] * ω) - 0.1sin(t + x[2] * ω)) + + 0.8(0.1sin(x[2] * ω) * ω - 0.1sin(t + x[2] * ω) * ω) + + 0.8^2 * (0.1sin(x[1] * ω) * ω - 0.1sin(t + x[1] * ω) * ω) + + 10.0(0.5 - 0.1cos(x[1] * ω) + 0.1cos(t + x[2] * ω) - 0.1cos(x[2] * ω) + + 0.1cos(t + x[1] * ω)) * + (0.1sin(x[1] * ω) * ω - 0.1sin(t + x[1] * ω) * ω) + + 10.0(0.5 - 0.1cos(x[1] * ω) + 0.1cos(t + x[2] * ω) - 0.1cos(x[2] * ω) + + 0.1cos(t + x[1] * ω)) * (-0.1sin(x[1] * ω) * ω + + 0.9 / 1.1 * (-0.1sin(t + x[1] * ω) * ω - 0.1cos(t + x[1] * ω) * ω) + + 1.0 / 1.1 * (0.1sin(t + x[1] * ω) * ω + 0.1cos(t + x[1] * ω) * ω)) + + du7 = -0.1cos(t + x[2] * ω) - 0.1sin(t + x[1] * ω) - 0.1sin(t + x[2] * ω) - + 0.1cos(t + x[1] * ω) - 0.1cos(t + x[2] * ω) * ω + + 0.8(-0.1sin(t + x[1] * ω) * ω - 0.1cos(t + x[1] * ω) * ω) - + 0.1sin(t + x[2] * ω) * ω + + 10.0(2.0 + 0.1cos(t + x[2] * ω) - 0.1sin(t + x[1] * ω) - + 0.1sin(t + x[2] * ω) + 0.1cos(t + x[1] * ω)) * + (-0.1cos(t + x[2] * ω) * ω - 0.1sin(t + x[2] * ω) * ω) + + (2.0 + 0.1cos(t + x[2] * ω) - 0.1sin(t + x[1] * ω) - 0.1sin(t + x[2] * ω) + + 0.1cos(t + x[1] * ω)) * cos(t + x[2] * ω) * ω + + du8 = 0.1cos(t + x[2] * ω) + 0.1sin(t + x[1] * ω) + 0.1sin(t + x[2] * ω) + + 0.1cos(t + x[1] * ω) + 0.1cos(t + x[2] * ω) * ω + + 0.8(0.1sin(t + x[1] * ω) * ω + 0.1cos(t + x[1] * ω) * ω) + + 0.1sin(t + x[2] * ω) * ω + + 10.0(0.5 - 0.1cos(t + x[2] * ω) + 0.1sin(t + x[1] * ω) + + 0.1sin(t + x[2] * ω) - 0.1cos(t + x[1] * ω)) * + (0.9(-0.1cos(t + x[2] * ω) * ω - 0.1sin(t + x[2] * ω) * ω) - + 0.1sin(t + x[2] * ω) * ω) + + 10.0(0.5 - 0.1cos(t + x[2] * ω) + 0.1sin(t + x[1] * ω) + + 0.1sin(t + x[2] * ω) - 0.1cos(t + x[1] * ω)) * + (0.1cos(t + x[2] * ω) * ω + 0.1sin(t + x[2] * ω) * ω) + + du9 = -0.1sin(t + x[1] * ω) - 0.1sin(t + x[2] * ω) + 0.1sin(x[2] * ω) * ω + + 0.8(0.1sin(x[1] * ω) * ω - 0.1sin(t + x[1] * ω) * ω) - + 0.1sin(t + x[2] * ω) * ω + + 10.0(0.5 - 0.1cos(x[1] * ω) + 0.1cos(t + x[2] * ω) - 0.1cos(x[2] * ω) + + 0.1cos(t + x[1] * ω)) * + (0.9 / 1.1 * (-0.1cos(t + x[2] * ω) * ω - 0.1sin(t + x[2] * ω) * ω) + + 1.0 / 1.1 * (0.1cos(t + x[2] * ω) * ω + 0.1sin(t + x[2] * ω) * ω) - + 0.1sin(x[2] * ω) * ω) + + 10.0(0.5 - 0.1cos(x[1] * ω) + 0.1cos(t + x[2] * ω) - 0.1cos(x[2] * ω) + + 0.1cos(t + x[1] * ω)) * (0.1sin(x[2] * ω) * ω - 0.1sin(t + x[2] * ω) * ω) + + return SVector(du1, du2, du3, du4, du5, du6, du7, du8, du9, zero(eltype(u))) +end + +""" + boundary_condition_slip_wall(u_inner, normal_direction, x, t, surface_flux_function, + equations::ShallowWaterMultiLayerEquations2D) +Create a boundary state by reflecting the normal velocity component and keep +the tangential velocity component unchanged. The boundary water height is taken from +the internal value. +For details see Section 9.2.5 of the book: +- Eleuterio F. Toro (2001) + Shock-Capturing Methods for Free-Surface Shallow Flows + 1st edition + ISBN 0471987662 +""" +@inline function Trixi.boundary_condition_slip_wall(u_inner, + normal_direction::AbstractVector, + x, t, + surface_flux_function, + equations::ShallowWaterMultiLayerEquations2D) + # normalize the outward pointing direction + normal = normal_direction / norm(normal_direction) + + # Extract internal values + h = waterheight(u_inner, equations) + h_v1, h_v2 = momentum(u_inner, equations) + b = u_inner[end] + + # u_boundary = SVector(h..., -hv..., b) + + # compute the normal velocity + u_normal = normal[1] * h_v1 + normal[2] * h_v2 + + # create the "external" boundary solution state + u_boundary = SVector(h..., + (h_v1 - 2.0 * u_normal * normal[1])..., + (h_v2 - 2.0 * u_normal * normal[2])..., + b) + + # calculate the boundary flux + flux = surface_flux_function(u_inner, u_boundary, normal_direction, equations) + + return flux +end + +""" + boundary_condition_slip_wall(u_inner, orientation, direction, x, t, + surface_flux_function, equations::ShallowWaterMultiLayerEquations2D) + +Should be used together with [`TreeMesh`](@ref). +""" +@inline function Trixi.boundary_condition_slip_wall(u_inner, orientation, + direction, x, t, + surface_flux_function, + equations::ShallowWaterMultiLayerEquations2D) + # Extract internal values + h = waterheight(u_inner, equations) + h_v1, h_v2 = momentum(u_inner, equations) + b = u_inner[end] + + ## get the appropriate normal vector from the orientation + if orientation == 1 + u_boundary = SVector(h..., -h_v1..., h_v2..., b) + else # orientation == 2 + u_boundary = SVector(h..., h_v1..., -h_v2..., b) + end + + # Calculate boundary flux + if iseven(direction) # u_inner is "left" of boundary, u_boundary is "right" of boundary + flux = surface_flux_function(u_inner, u_boundary, orientation, equations) + else # u_boundary is "left" of boundary, u_inner is "right" of boundary + flux = surface_flux_function(u_boundary, u_inner, orientation, equations) + end + + return flux +end + +# Calculate 1D advective portion of the flux for a single point +# Note, the bottom topography has no flux +@inline function Trixi.flux(u, orientation::Integer, + equations::ShallowWaterMultiLayerEquations2D) + + # Extract waterheight and momentum and compute velocities + h_v1, h_v2 = momentum(u, equations) + v1, v2 = velocity(u, equations) + + # Initialize flux vector + f = zero(MVector{3 * nlayers(equations) + 1, real(equations)}) + + # Calculate fluxes in each layer + # The momentum flux simplifies as the pressure is included in the nonconservative term + if orientation == 1 + for i in eachlayer(equations) + f_h = h_v1[i] + f_hv1 = f_h * v1[i] + f_hv2 = f_h * v2[i] + setlayer!(f, f_h, f_hv1, f_hv2, i, equations) + end + else # orientation == 2 + for i in eachlayer(equations) + f_h = h_v2[i] + f_hv1 = f_h * v1[i] + f_hv2 = f_h * v2[i] + setlayer!(f, f_h, f_hv1, f_hv2, i, equations) + end + end + + return SVector(f) +end + +# Calculate 1D flux for a single point in the normal direction +# Note, this directional vector is not normalized and the bottom topography has no flux +@inline function Trixi.flux(u, normal_direction::AbstractVector, + equations::ShallowWaterMultiLayerEquations2D) + # Extract waterheight and momentum and compute velocities + h = waterheight(u, equations) + v1, v2 = velocity(u, equations) + + v_normal = v1 * normal_direction[1] + v2 * normal_direction[2] + h_v_normal = h .* v_normal + + # Initialize flux vector + f = zero(MVector{3 * nlayers(equations) + 1, real(equations)}) + + # Calculate fluxes in each layer + # The momentum flux simplifies as the pressure is included in the nonconservative term + for i in eachlayer(equations) + f_h = h_v_normal[i] + f_hv1 = f_h * v1[i] + f_hv2 = f_h * v2[i] + setlayer!(f, f_h, f_hv1, f_hv2, i, equations) + end + + return SVector(f) +end + +""" + flux_nonconservative_ersing_etal(u_ll, u_rr, orientation::Integer, + equations::ShallowWaterMultiLayerEquations2D) + flux_nonconservative_ersing_etal(u_ll, u_rr, + normal_direction_ll::AbstractVector, + normal_direction_average::AbstractVector, + equations::ShallowWaterMultiLayerEquations2D) + +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 [`ShallowWaterMultiLayerEquations2D`](@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::ShallowWaterMultiLayerEquations2D) + # 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{3 * nlayers(equations) + 1, real(equations)}) + + # Compute the nonconservative flux in each layer + # 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