From 935081c434ae0135b8a53d466da3a1f30c2f8e41 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Sat, 17 Aug 2024 20:24:35 +0200 Subject: [PATCH] variable bathymetry Serre-Green-Naghdi equations (#135) * variable bathymetry Serre-Green-Naghdi equations * Dingemans * format * various fixes * Apply suggestions from code review Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> * move function waterheight * bathymetry -> bathymetry_type * add bathymetry_type to elixirs and well-balanced test * move comments * Apply suggestions from code review Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> * format * Apply suggestions from code review Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> * Apply suggestions from code review Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> * Update src/equations/serre_green_naghdi_1d.jl Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> * use D1 for SBP op. * (h + b) -> eta * explain mild slope * -= * format * update Dingemans test values * fix bug * add new test used for benchmarks * fix typo * Apply suggestions from code review Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> * Update examples/serre_green_naghdi_1d/serre_green_naghdi_well_balanced.jl Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> * assemble_system_matrix! * solve_system_matrix! * Update src/equations/serre_green_naghdi_1d.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> * update Dingemans * Update src/equations/serre_green_naghdi_1d.jl Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> * Update src/equations/serre_green_naghdi_1d.jl Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> * fix docstring of energy_total_modified * adapt tolerances for macOS CI * not only integrals --------- Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- .gitignore | 1 + ...bm_variable_bathymetry_1d_well_balanced.jl | 1 - .../serre_green_naghdi_conservation.jl | 69 ++ .../serre_green_naghdi_dingemans.jl | 49 ++ .../serre_green_naghdi_soliton.jl | 4 +- .../serre_green_naghdi_soliton_fourier.jl | 4 +- .../serre_green_naghdi_soliton_relaxation.jl | 4 +- .../serre_green_naghdi_soliton_upwind.jl | 4 +- .../serre_green_naghdi_well_balanced.jl | 66 ++ src/DispersiveShallowWater.jl | 2 +- src/equations/bbm_bbm_1d.jl | 4 - .../bbm_bbm_variable_bathymetry_1d.jl | 11 - src/equations/equations.jl | 53 +- src/equations/serre_green_naghdi_1d.jl | 737 ++++++++++++++---- src/equations/svaerd_kalisch_1d.jl | 11 - test/test_serre_green_naghdi_1d.jl | 207 +++++ 16 files changed, 1041 insertions(+), 186 deletions(-) create mode 100644 examples/serre_green_naghdi_1d/serre_green_naghdi_conservation.jl create mode 100644 examples/serre_green_naghdi_1d/serre_green_naghdi_dingemans.jl create mode 100644 examples/serre_green_naghdi_1d/serre_green_naghdi_well_balanced.jl diff --git a/.gitignore b/.gitignore index c1dd5bb5..ffeb3567 100644 --- a/.gitignore +++ b/.gitignore @@ -8,3 +8,4 @@ out*/ /docs/build/ run run/* +**/*.code-workspace diff --git a/examples/bbm_bbm_variable_bathymetry_1d/bbm_bbm_variable_bathymetry_1d_well_balanced.jl b/examples/bbm_bbm_variable_bathymetry_1d/bbm_bbm_variable_bathymetry_1d_well_balanced.jl index ab640e15..08b8d26c 100644 --- a/examples/bbm_bbm_variable_bathymetry_1d/bbm_bbm_variable_bathymetry_1d_well_balanced.jl +++ b/examples/bbm_bbm_variable_bathymetry_1d/bbm_bbm_variable_bathymetry_1d_well_balanced.jl @@ -57,7 +57,6 @@ analysis_callback = AnalysisCallback(semi; interval = 10, extra_analysis_integrals = (waterheight_total, velocity, entropy, lake_at_rest_error)) -# Always put relaxation_callback before analysis_callback to guarantee conservation of the invariant callbacks = CallbackSet(analysis_callback, summary_callback) dt = 0.5 diff --git a/examples/serre_green_naghdi_1d/serre_green_naghdi_conservation.jl b/examples/serre_green_naghdi_1d/serre_green_naghdi_conservation.jl new file mode 100644 index 00000000..709f3958 --- /dev/null +++ b/examples/serre_green_naghdi_1d/serre_green_naghdi_conservation.jl @@ -0,0 +1,69 @@ +# This elixir contains an artificial setup that can be used to check the +# conservation properties of the equations and numerical methods as well as +# a possible directional bias (if the velocity is set to zero). See +# - Hendrik Ranocha and Mario Ricchiuto (2024) +# Structure-preserving approximations of the Serre-Green-Naghdi +# equations in standard and hyperbolic form +# [arXiv: 2408.02665](https://arxiv.org/abs/2408.02665) + +using OrdinaryDiffEq +using DispersiveShallowWater +using SummationByPartsOperators: upwind_operators, periodic_derivative_operator + +############################################################################### +# Semidiscretization of the Serre-Green-Naghdi equations + +bathymetry_type = bathymetry_variable # or bathymetry_mild_slope +equations = SerreGreenNaghdiEquations1D(bathymetry_type; + gravity_constant = 9.81, + eta0 = 1.0) + +function initial_condition_conservation_test(x, t, + equations::SerreGreenNaghdiEquations1D, + mesh) + eta = 1 + exp(-x^2) + v = 1.0e-2 # set this to zero to test a directional bias + b = 0.25 * cospi(x / 75) + + D = equations.eta0 - b + return SVector(eta, v, D) +end + +# create homogeneous mesh +coordinates_min = -150.0 +coordinates_max = +150.0 +N = 1_000 +mesh = Mesh1D(coordinates_min, coordinates_max, N) + +# create solver with periodic SBP operators of accuracy order 2 +accuracy_order = 2 +solver = Solver(mesh, accuracy_order) + +# semidiscretization holds all the necessary data structures for the spatial discretization +semi = Semidiscretization(mesh, equations, + initial_condition_conservation_test, solver; + boundary_conditions = boundary_condition_periodic) + +############################################################################### +# Create `ODEProblem` and run the simulation +tspan = (0.0, 35.0) +ode = semidiscretize(semi, tspan) + +# The callbacks support an additional `io` argument to write output to a file +# or any other IO stream. The default is stdout. We use this here to enable +# setting it to `devnull` to benchmark the full simulation including the time +# to compute the errors etc. but without the time to write the output to the +# terminal. +io = stdout +summary_callback = SummaryCallback(io) +analysis_callback = AnalysisCallback(semi; interval = 10, io, + extra_analysis_errors = (:conservation_error,), + extra_analysis_integrals = (waterheight_total, + momentum, + entropy, + entropy_modified)) + +callbacks = CallbackSet(analysis_callback, summary_callback) + +sol = solve(ode, Tsit5(); + save_everystep = false, callback = callbacks) diff --git a/examples/serre_green_naghdi_1d/serre_green_naghdi_dingemans.jl b/examples/serre_green_naghdi_1d/serre_green_naghdi_dingemans.jl new file mode 100644 index 00000000..ab4e0692 --- /dev/null +++ b/examples/serre_green_naghdi_1d/serre_green_naghdi_dingemans.jl @@ -0,0 +1,49 @@ +using OrdinaryDiffEq +using DispersiveShallowWater +using SummationByPartsOperators: upwind_operators, periodic_derivative_operator + +############################################################################### +# Semidiscretization of the Serre-Green-Naghdi equations + +bathymetry_type = bathymetry_variable # or bathymetry_mild_slope +equations = SerreGreenNaghdiEquations1D(bathymetry_type; + gravity_constant = 9.81) + +initial_condition = initial_condition_dingemans +boundary_conditions = boundary_condition_periodic + +# create homogeneous mesh +coordinates_min = -138.0 +coordinates_max = 46.0 +N = 512 +mesh = Mesh1D(coordinates_min, coordinates_max, N) + +# create solver with periodic upwind SBP operators of accuracy order 4 +# (note that only central-type with even order are used) +D1 = upwind_operators(periodic_derivative_operator; + derivative_order = 1, accuracy_order = 3, + xmin = xmin(mesh), xmax = xmax(mesh), + N = nnodes(mesh)) +solver = Solver(D1, nothing) + +# semidiscretization holds all the necessary data structures for the spatial discretization +semi = Semidiscretization(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_conditions) + +############################################################################### +# Create `ODEProblem` and run the simulation +tspan = (0.0, 70.0) +ode = semidiscretize(semi, tspan) +summary_callback = SummaryCallback() +analysis_callback = AnalysisCallback(semi; interval = 10, + extra_analysis_errors = (:conservation_error,), + extra_analysis_integrals = (waterheight_total, + momentum, + entropy, + entropy_modified)) + +callbacks = CallbackSet(analysis_callback, summary_callback) + +saveat = range(tspan..., length = 500) +sol = solve(ode, Tsit5(), abstol = 1e-7, reltol = 1e-7, + save_everystep = false, callback = callbacks, saveat = saveat) diff --git a/examples/serre_green_naghdi_1d/serre_green_naghdi_soliton.jl b/examples/serre_green_naghdi_1d/serre_green_naghdi_soliton.jl index 033197d7..bfcaaddf 100644 --- a/examples/serre_green_naghdi_1d/serre_green_naghdi_soliton.jl +++ b/examples/serre_green_naghdi_1d/serre_green_naghdi_soliton.jl @@ -4,7 +4,9 @@ using DispersiveShallowWater ############################################################################### # Semidiscretization of the Serre-Green-Naghdi equations -equations = SerreGreenNaghdiEquations1D(gravity_constant = 9.81) +bathymetry_type = bathymetry_flat +equations = SerreGreenNaghdiEquations1D(bathymetry_type; + gravity_constant = 9.81) # initial_condition_convergence_test needs periodic boundary conditions initial_condition = initial_condition_convergence_test diff --git a/examples/serre_green_naghdi_1d/serre_green_naghdi_soliton_fourier.jl b/examples/serre_green_naghdi_1d/serre_green_naghdi_soliton_fourier.jl index d22196f9..b8f4bf7f 100644 --- a/examples/serre_green_naghdi_1d/serre_green_naghdi_soliton_fourier.jl +++ b/examples/serre_green_naghdi_1d/serre_green_naghdi_soliton_fourier.jl @@ -5,7 +5,9 @@ using SummationByPartsOperators: fourier_derivative_operator ############################################################################### # Semidiscretization of the Serre-Green-Naghdi equations -equations = SerreGreenNaghdiEquations1D(gravity_constant = 9.81) +bathymetry_type = bathymetry_flat +equations = SerreGreenNaghdiEquations1D(bathymetry_type; + gravity_constant = 9.81) # initial_condition_convergence_test needs periodic boundary conditions initial_condition = initial_condition_convergence_test diff --git a/examples/serre_green_naghdi_1d/serre_green_naghdi_soliton_relaxation.jl b/examples/serre_green_naghdi_1d/serre_green_naghdi_soliton_relaxation.jl index e095c9ff..d484a733 100644 --- a/examples/serre_green_naghdi_1d/serre_green_naghdi_soliton_relaxation.jl +++ b/examples/serre_green_naghdi_1d/serre_green_naghdi_soliton_relaxation.jl @@ -5,7 +5,9 @@ using SummationByPartsOperators: fourier_derivative_operator ############################################################################### # Semidiscretization of the Serre-Green-Naghdi equations -equations = SerreGreenNaghdiEquations1D(gravity_constant = 9.81) +bathymetry_type = bathymetry_flat +equations = SerreGreenNaghdiEquations1D(bathymetry_type; + gravity_constant = 9.81) # initial_condition_convergence_test needs periodic boundary conditions initial_condition = initial_condition_convergence_test diff --git a/examples/serre_green_naghdi_1d/serre_green_naghdi_soliton_upwind.jl b/examples/serre_green_naghdi_1d/serre_green_naghdi_soliton_upwind.jl index 8bb14375..cd1c3eb7 100644 --- a/examples/serre_green_naghdi_1d/serre_green_naghdi_soliton_upwind.jl +++ b/examples/serre_green_naghdi_1d/serre_green_naghdi_soliton_upwind.jl @@ -6,7 +6,9 @@ using SparseArrays: sparse ############################################################################### # Semidiscretization of the Serre-Green-Naghdi equations -equations = SerreGreenNaghdiEquations1D(gravity_constant = 9.81) +bathymetry_type = bathymetry_flat +equations = SerreGreenNaghdiEquations1D(bathymetry_type; + gravity_constant = 9.81) # initial_condition_convergence_test needs periodic boundary conditions initial_condition = initial_condition_convergence_test diff --git a/examples/serre_green_naghdi_1d/serre_green_naghdi_well_balanced.jl b/examples/serre_green_naghdi_1d/serre_green_naghdi_well_balanced.jl new file mode 100644 index 00000000..aac19ffb --- /dev/null +++ b/examples/serre_green_naghdi_1d/serre_green_naghdi_well_balanced.jl @@ -0,0 +1,66 @@ +using OrdinaryDiffEq +using DispersiveShallowWater +using SummationByPartsOperators: upwind_operators, periodic_derivative_operator + +############################################################################### +# Semidiscretization of the Serre-Green-Naghdi equations + +bathymetry_type = bathymetry_variable # or bathymetry_mild_slope +equations = SerreGreenNaghdiEquations1D(bathymetry_type; + gravity_constant = 1.0, eta0 = 2.0) + +# Setup a truly discontinuous bottom topography function for this academic +# testcase of well-balancedness. The errors from the analysis callback are +# not important but the error for this lake-at-rest test case +# `∫|η-η₀|` should be around machine roundoff. +function initial_condition_discontinuous_well_balancedness(x, t, + equations::SerreGreenNaghdiEquations1D, + mesh) + # Set the background values + eta = equations.eta0 + v = 0.0 + D = equations.eta0 - 1.0 + + # Setup a discontinuous bottom topography + if x >= 0.5 && x <= 0.75 + D = equations.eta0 - 1.5 - 0.5 * sinpi(2.0 * x) + end + + return SVector(eta, v, D) +end + +initial_condition = initial_condition_discontinuous_well_balancedness +boundary_conditions = boundary_condition_periodic + +# create homogeneous mesh +coordinates_min = -1.0 +coordinates_max = 1.0 +N = 512 +mesh = Mesh1D(coordinates_min, coordinates_max, N) + +# create solver with periodic upwind SBP operators +D1 = upwind_operators(periodic_derivative_operator, + derivative_order = 1, + accuracy_order = 3, # the resulting central operators have order 4 + xmin = xmin(mesh), xmax = xmax(mesh), + N = nnodes(mesh)) +solver = Solver(D1, nothing) + +# semidiscretization holds all the necessary data structures for the spatial discretization +semi = Semidiscretization(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_conditions) + +############################################################################### +# Create `ODEProblem` and run the simulation +tspan = (0.0, 30.0) +ode = semidiscretize(semi, tspan) +summary_callback = SummaryCallback() +analysis_callback = AnalysisCallback(semi; interval = 10, + extra_analysis_errors = (:conservation_error,), + extra_analysis_integrals = (waterheight_total, + momentum, + entropy_modified, + lake_at_rest_error)) +callbacks = CallbackSet(analysis_callback, summary_callback) + +sol = solve(ode, Tsit5(), dt = 0.5, adaptive = false, callback = callbacks) diff --git a/src/DispersiveShallowWater.jl b/src/DispersiveShallowWater.jl index 294bb2e5..f5b139a2 100644 --- a/src/DispersiveShallowWater.jl +++ b/src/DispersiveShallowWater.jl @@ -80,7 +80,7 @@ export Semidiscretization, semidiscretize, grid export boundary_condition_periodic, boundary_condition_reflecting -export bathymetry_flat +export bathymetry_flat, bathymetry_mild_slope, bathymetry_variable export initial_condition_convergence_test, initial_condition_manufactured, source_terms_manufactured, diff --git a/src/equations/bbm_bbm_1d.jl b/src/equations/bbm_bbm_1d.jl index 594607db..18e1ffc2 100644 --- a/src/equations/bbm_bbm_1d.jl +++ b/src/equations/bbm_bbm_1d.jl @@ -288,8 +288,4 @@ end return equations.eta0 - equations.D end -@inline function waterheight(q, equations::BBMBBMEquations1D) - return waterheight_total(q, equations) - bathymetry(q, equations) -end - @inline entropy(q, equations::BBMBBMEquations1D) = energy_total(q, equations) diff --git a/src/equations/bbm_bbm_variable_bathymetry_1d.jl b/src/equations/bbm_bbm_variable_bathymetry_1d.jl index e90ee896..11b9f3d7 100644 --- a/src/equations/bbm_bbm_variable_bathymetry_1d.jl +++ b/src/equations/bbm_bbm_variable_bathymetry_1d.jl @@ -374,15 +374,4 @@ end return equations.eta0 - D end -@inline function waterheight(q, equations::BBMBBMVariableEquations1D) - return waterheight_total(q, equations) - bathymetry(q, equations) -end - @inline entropy(q, equations::BBMBBMVariableEquations1D) = energy_total(q, equations) - -# Calculate the error for the "lake-at-rest" test case where eta should -# be a constant value over time -@inline function lake_at_rest_error(q, equations::BBMBBMVariableEquations1D) - eta, _, _ = q - return abs(equations.eta0 - eta) -end diff --git a/src/equations/equations.jl b/src/equations/equations.jl index 3ceae561..90be0fdc 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -96,8 +96,8 @@ function cons2prim end waterheight_total(q, equations) Return the total waterheight of the primitive variables `q` for a given set of -`equations`, i.e., the [`waterheight`](@ref) plus the -[`bathymetry`](@ref). +`equations`, i.e., the [`waterheight`](@ref) ``h`` plus the +[`bathymetry`](@ref) ``b``. `q` is a vector of the primitive variables at a single node, i.e., a vector of the correct length `nvariables(equations)`. @@ -107,15 +107,19 @@ function waterheight_total end varnames(::typeof(waterheight_total), equations) = ("η",) """ - waterheight(q, equations) + waterheight(q, equations::AbstractShallowWaterEquations) Return the waterheight of the primitive variables `q` for a given set of -`equations`, i.e., the waterheight above the bathymetry. +`equations`, i.e., the waterheight ``h`` above the bathymetry ``b``. `q` is a vector of the primitive variables at a single node, i.e., a vector of the correct length `nvariables(equations)`. + +See also [`waterheight_total`](@ref), [`bathymetry`](@ref). """ -function waterheight end +@inline function waterheight(q, equations::AbstractShallowWaterEquations) + return waterheight_total(q, equations) - bathymetry(q, equations) +end varnames(::typeof(waterheight), equations) = ("h",) @@ -184,6 +188,19 @@ of the correct length `nvariables(equations)`. """ function bathymetry end +""" + lake_at_rest_error(q, equations::AbstractShallowWaterEquations) + +Calculate the error for the "lake-at-rest" test case where the +[`waterheight_total`](@ref) ``\\eta = h + b`` should +be a constant value over time (given by the value ``\\eta_0`` passed to the +`equations` when constructing them). +""" +@inline function lake_at_rest_error(q, equations::AbstractShallowWaterEquations) + eta = waterheight_total(q, equations) + return abs(equations.eta0 - eta) +end + """ entropy(q, equations) @@ -325,14 +342,40 @@ Default analysis integrals used by the [`AnalysisCallback`](@ref). default_analysis_integrals(::AbstractEquations) = Symbol[] abstract type AbstractBathymetry end + struct BathymetryFlat <: AbstractBathymetry end """ bathymetry_flat = DispersiveShallowWater.BathymetryFlat() A singleton struct indicating a flat bathymetry. + +See also [`bathymetry_mild_slope`](@ref) and [`bathymetry_variable`](@ref). """ const bathymetry_flat = BathymetryFlat() +struct BathymetryMildSlope <: AbstractBathymetry end +""" + bathymetry_mild_slope = DispersiveShallowWater.BathymetryMildSlope() + +A singleton struct indicating a variable bathymetry with mild-slope +approximation. Typically, this means that some terms like ``b_x^2`` +are neglected. + +See also [`bathymetry_flat`](@ref) and [`bathymetry_variable`](@ref). +""" +const bathymetry_mild_slope = BathymetryMildSlope() + +struct BathymetryVariable <: AbstractBathymetry end +""" + bathymetry_variable = DispersiveShallowWater.BathymetryVariable() + +A singleton struct indicating a variable bathymetry (without +mild-slope approximation). + +See also [`bathymetry_flat`](@ref) and [`bathymetry_mild_slope`](@ref). +""" +const bathymetry_variable = BathymetryVariable() + # BBM-BBM equations abstract type AbstractBBMBBMEquations{NDIMS, NVARS} <: AbstractShallowWaterEquations{NDIMS, NVARS} end diff --git a/src/equations/serre_green_naghdi_1d.jl b/src/equations/serre_green_naghdi_1d.jl index 0858679f..a3387ab8 100644 --- a/src/equations/serre_green_naghdi_1d.jl +++ b/src/equations/serre_green_naghdi_1d.jl @@ -1,5 +1,6 @@ @doc raw""" - SerreGreenNaghdiEquations1D(; gravity_constant, eta0 = 0.0) + SerreGreenNaghdiEquations1D(bathymetry_type = bathymetry_flat; + gravity_constant, eta0 = 0.0) Serre-Green-Naghdi system in one spatial dimension. The equations for flat bathymetry are given by @@ -18,6 +19,36 @@ is therefore given by ``h = \eta - \eta_0 + D``. The Serre-Green-Naghdi equations are only implemented for ``\eta_0 = 0``. The total water height is therefore given by ``\eta = h + b``. +Three types of variable `bathymetry_type` are supported: +- [`bathymetry_flat`](@ref): flat bathymetry (typically ``b = 0`` everywhere) +- [`bathymetry_mild_slope`](@ref): variable bathymetry with mild-slope approximation +- [`bathymetry_variable`](@ref): general variable bathymetry + +For the mild-slope approximation, the Serre-Green-Naghdi equations are +```math +\begin{aligned} + h_t + (h v)_x &= 0,\\ + h v_t - \frac{1}{3} (h^3 v_{tx})_x + \frac{1}{2} (h^2 b_x u_t)_x - \frac{1}{2} h^2 b_x u_{tx} + \frac{3}{4} h b_x^2 u_t + + \frac{1}{2} g (h^2)_x + g h b_x + \frac{1}{2} h (v^2)_x + + p_x + \frac{3}{2} \frac{p}{h} b_x &= 0,\\ + p &= \frac{1}{3} h^3 v_{x}^2 - \frac{1}{3} h^3 v v_{xx} + + \frac{1}{2} h^2 v (b_x v)_x. +\end{aligned} +``` +For the general case of variable vathymetry without mild-slope +approximation, the Serre-Green-Naghdi equations are +```math +\begin{aligned} + h_t + (h v)_x &= 0,\\ + h v_t - \frac{1}{3} (h^3 v_{tx})_x + \frac{1}{2} (h^2 b_x u_t)_x - \frac{1}{2} h^2 b_x u_{tx} + h b_x^2 u_t + + \frac{1}{2} g (h^2)_x + g h b_x + \frac{1}{2} h (v^2)_x + + p_x + \frac{3}{2} \frac{p}{h} b_x + \psi b_x &= 0,\\ + p &= \frac{1}{3} h^3 v_{x}^2 - \frac{1}{3} h^3 v v_{xx} + + \frac{1}{2} h^2 v (b_x v)_x, + \psi &= \frac{1}{4} h v (b_x v)_x. +\end{aligned} +``` + References for the Serre-Green-Naghdi system can be found in - Serre (1953) Contribution â l'étude des écoulements permanents et variables dans les canaux @@ -28,10 +59,11 @@ References for the Serre-Green-Naghdi system can be found in The semidiscretization implemented here conserves - the total water mass (integral of h) as a linear invariant -- the total momentum (integral of h v) as a nonlinear invariant -- the total energy +- the total momentum (integral of h v) as a nonlinear invariant if the bathymetry is constant +- the total modified energy -for periodic boundary conditions, see +for periodic boundary conditions (see Ranocha and Ricchiuto (2024)). +Additionally, it is well-balanced for the lake-at-rest stationary solution, see - Hendrik Ranocha and Mario Ricchiuto (2024) Structure-preserving approximations of the Serre-Green-Naghdi equations in standard and hyperbolic form @@ -39,15 +71,24 @@ for periodic boundary conditions, see """ struct SerreGreenNaghdiEquations1D{Bathymetry <: AbstractBathymetry, RealT <: Real} <: AbstractSerreGreenNaghdiEquations{1, 3} - bathymetry::Bathymetry # type of bathymetry + bathymetry_type::Bathymetry # type of bathymetry gravity::RealT # gravitational constant eta0::RealT # constant still-water surface end -function SerreGreenNaghdiEquations1D(; gravity_constant, eta0 = 0.0) - eta0 == 0.0 || - @warn "The still-water surface needs to be 0 for the Serre-Green-Naghdi equations" - SerreGreenNaghdiEquations1D(bathymetry_flat, gravity_constant, eta0) +function SerreGreenNaghdiEquations1D(bathymetry_type = bathymetry_flat; + gravity_constant, eta0 = 0.0) + SerreGreenNaghdiEquations1D(bathymetry_type, gravity_constant, eta0) +end + +function get_name(equations::SerreGreenNaghdiEquations1D) + name = equations |> typeof |> nameof |> string + bathymetry_type = equations.bathymetry_type + if bathymetry_type isa BathymetryFlat + return name + else # variable bathymetry_type + return name * "-" * string(nameof(typeof(bathymetry_type))) + end end varnames(::typeof(prim2prim), ::SerreGreenNaghdiEquations1D) = ("η", "v", "D") @@ -75,16 +116,59 @@ function initial_condition_convergence_test(x, t, equations::SerreGreenNaghdiEqu return SVector(h, v, 0) end +""" + initial_condition_dingemans(x, t, equations::SerreGreenNaghdiEquations1D, mesh) + +The initial condition that uses the dispersion relation of the Euler equations +to approximate waves generated by a wave maker as it is done by experiments of +Dingemans for a flow over a trapezoidal bathymetry. +It is assumed that `equations.eta0 = 0.8`. + +References: +- Magnus Svärd, Henrik Kalisch (2023) + A novel energy-bounded Boussinesq model and a well-balanced and stable numerical discretization + [arXiv: 2302.09924](https://arxiv.org/abs/2302.09924) +- Maarten W. Dingemans (1994) + Comparison of computations with Boussinesq-like models and laboratory measurements + [link](https://repository.tudelft.nl/islandora/object/uuid:c2091d53-f455-48af-a84b-ac86680455e9/datastream/OBJ/download) +""" +function initial_condition_dingemans(x, t, equations::SerreGreenNaghdiEquations1D, mesh) + h0 = 0.8 + A = 0.02 + # omega = 2*pi/(2.02*sqrt(2)) + k = 0.8406220896381442 # precomputed result of find_zero(k -> omega^2 - equations.gravity * k * tanh(k * h0), 1.0) using Roots.jl + if x < -34.5 * pi / k || x > -8.5 * pi / k + h = 0.0 + else + h = A * cos(k * x) + end + v = sqrt(equations.gravity / k * tanh(k * h0)) * h / h0 + if 11.01 <= x && x < 23.04 + b = 0.6 * (x - 11.01) / (23.04 - 11.01) + elseif 23.04 <= x && x < 27.04 + b = 0.6 + elseif 27.04 <= x && x < 33.07 + b = 0.6 * (33.07 - x) / (33.07 - 27.04) + else + b = 0.0 + end + eta = h + h0 + D = equations.eta0 - b + return SVector(eta, v, D) +end + +# flat bathymetry function create_cache(mesh, equations::SerreGreenNaghdiEquations1D{BathymetryFlat}, solver, initial_condition, ::BoundaryConditionPeriodic, RealT, uEltype) - D = solver.D1 + D1 = solver.D1 # create temporary storage h = ones(RealT, nnodes(mesh)) + b = zero(h) h_x = zero(h) v_x = zero(h) h2_x = zero(h) @@ -97,62 +181,203 @@ function create_cache(mesh, M_h = zero(h) M_h3_3 = zero(h) - if D isa PeriodicUpwindOperators + if D1 isa PeriodicUpwindOperators v_x_upwind = zero(h) - Dmat_minus = sparse(D.minus) - - # Floating point errors accumulate a bit and the system matrix - # is not necessarily perfectly symmetric but only up to - # round-off errors. We wrap it here to avoid issues with the - # factorization. - @. M_h = h - scale_by_mass_matrix!(M_h, D) - @. M_h3_3 = (1 / 3) * h^3 - scale_by_mass_matrix!(M_h3_3, D) - system_matrix = Symmetric(Diagonal(M_h) - + - Dmat_minus' * Diagonal(M_h3_3) * Dmat_minus) + D1mat_minus = sparse(D1.minus) + + system_matrix = let cache = (; M_h, M_h3_3) + assemble_system_matrix!(cache, h, + D1, D1mat_minus, equations) + end factorization = cholesky(system_matrix) - cache = (; h_x, v_x, v_x_upwind, h2_x, hv_x, v2_x, + cache = (; h, b, h_x, v_x, v_x_upwind, h2_x, hv_x, v2_x, h2_v_vx_x, h_vx_x, p_x, tmp, M_h, M_h3_3, - D, Dmat_minus, factorization) + D1, D1mat_minus, factorization) else - if D isa FourierDerivativeOperator - Dmat = Matrix(D) + if D1 isa FourierDerivativeOperator + D1mat = Matrix(D1) - cache = (; h_x, v_x, h2_x, hv_x, v2_x, + cache = (; h, b, h_x, v_x, h2_x, hv_x, v2_x, h2_v_vx_x, h_vx_x, p_x, tmp, M_h, M_h3_3, - D, Dmat) + D1, D1mat) else - Dmat = sparse(D) - - # Floating point errors accumulate a bit and the system matrix - # is not necessarily perfectly symmetric but only up to - # round-off errors. We wrap it here to avoid issues with the - # factorization. - @. M_h = h - scale_by_mass_matrix!(M_h, D) - @. M_h3_3 = (1 / 3) * h^3 - scale_by_mass_matrix!(M_h3_3, D) - system_matrix = Symmetric(Diagonal(M_h) - + - Dmat' * Diagonal(M_h3_3) * Dmat) + D1mat = sparse(D1) + + system_matrix = let cache = (; M_h, M_h3_3) + assemble_system_matrix!(cache, h, + D1, D1mat, equations) + end factorization = cholesky(system_matrix) - cache = (; h_x, v_x, h2_x, hv_x, v2_x, + cache = (; h, b, h_x, v_x, h2_x, hv_x, v2_x, h2_v_vx_x, h_vx_x, p_x, tmp, M_h, M_h3_3, - D, Dmat, factorization) + D1, D1mat, factorization) end end return cache end +# variable bathymetry +function create_cache(mesh, + equations::SerreGreenNaghdiEquations1D, + solver, + initial_condition, + ::BoundaryConditionPeriodic, + RealT, uEltype) + D1 = solver.D1 + + # create temporary storage + h = ones(RealT, nnodes(mesh)) + h_x = zero(h) + b = zero(h) + b_x = zero(h) + v_x = zero(h) + h_hpb_x = zero(h) + hv_x = zero(h) + v2_x = zero(h) + h2_v_vx_x = zero(h) + h_vx_x = zero(h) + p_h = zero(h) + p_x = zero(h) + tmp = zero(h) + M_h_p_h_bx2 = zero(h) + M_h3_3 = zero(h) + M_h2_bx = zero(h) + + # b_x appears in the system matrix, so it must not be zero to get the + # correct pattern of the factorization. Just setting it to unity does + # also not work since some terms cancel. Thus, b_x needs to be + # "variable enough" to get the correct pattern of non-zero entries in + # the system matrix and the factorization. + # TODO: This is a hack and should be improved. It would be ideal if we + # had access to the initial condition and the exact value of the + # bathymetry here. + let x = grid(D1) + @. b_x = x^3 + end + + if D1 isa PeriodicUpwindOperators + v_x_upwind = zero(h) + p_0 = zero(h) + + D1mat_minus = sparse(D1.minus) + + system_matrix = let cache = (; M_h_p_h_bx2, M_h3_3, M_h2_bx) + assemble_system_matrix!(cache, h, b_x, + D1, D1mat_minus, equations) + end + factorization = cholesky(system_matrix) + + cache = (; h, h_x, v_x, v_x_upwind, h_hpb_x, b, b_x, hv_x, v2_x, + h2_v_vx_x, h_vx_x, p_h, p_0, p_x, tmp, + M_h_p_h_bx2, M_h3_3, M_h2_bx, + D1, D1mat_minus, factorization) + else + if D1 isa FourierDerivativeOperator + D1mat = Matrix(D1) + + cache = (; h, h_x, v_x, h_hpb_x, b, b_x, hv_x, v2_x, + h2_v_vx_x, h_vx_x, p_h, p_x, tmp, + M_h_p_h_bx2, M_h3_3, M_h2_bx, + D1, D1mat) + else + D1mat = sparse(D1) + + system_matrix = let cache = (; M_h_p_h_bx2, M_h3_3, M_h2_bx) + assemble_system_matrix!(cache, h, b_x, + D1, D1mat, equations) + end + factorization = cholesky(system_matrix) + + cache = (; h, h_x, v_x, h_hpb_x, b, b_x, hv_x, v2_x, + h2_v_vx_x, h_vx_x, p_h, p_x, tmp, + M_h_p_h_bx2, M_h3_3, M_h2_bx, + D1, D1mat, factorization) + end + end + + if equations.bathymetry_type isa BathymetryVariable + psi = zero(h) + cache = (; cache..., psi) + end + + return cache +end + +function assemble_system_matrix!(cache, h, D1, D1mat, + equations::SerreGreenNaghdiEquations1D{BathymetryFlat}) + (; M_h, M_h3_3) = cache + + @. M_h = h + scale_by_mass_matrix!(M_h, D1) + @. M_h3_3 = (1 / 3) * h^3 + scale_by_mass_matrix!(M_h3_3, D1) + + # Floating point errors accumulate a bit and the system matrix + # is not necessarily perfectly symmetric but only up to + # round-off errors. We wrap it here to avoid issues with the + # factorization. + return Symmetric(Diagonal(M_h) + D1mat' * Diagonal(M_h3_3) * D1mat) +end + +# variable bathymetry +function assemble_system_matrix!(cache, h, b_x, D1, D1mat, + equations::SerreGreenNaghdiEquations1D) + (; M_h_p_h_bx2, M_h3_3, M_h2_bx) = cache + + if equations.bathymetry_type isa BathymetryMildSlope + factor = 0.75 + elseif equations.bathymetry_type isa BathymetryVariable + factor = 1.0 + end + @. M_h_p_h_bx2 = h + factor * h * b_x^2 + scale_by_mass_matrix!(M_h_p_h_bx2, D1) + inv3 = 1 / 3 + @. M_h3_3 = inv3 * h^3 + scale_by_mass_matrix!(M_h3_3, D1) + @. M_h2_bx = 0.5 * h^2 * b_x + scale_by_mass_matrix!(M_h2_bx, D1) + # Floating point errors accumulate a bit and the system matrix + # is not necessarily perfectly symmetric but only up to + # round-off errors. We wrap it here to avoid issues with the + # factorization. + return Symmetric(Diagonal(M_h_p_h_bx2) + + + D1mat' * (Diagonal(M_h3_3) * D1mat + - + Diagonal(M_h2_bx)) + - + Diagonal(M_h2_bx) * D1mat) +end + +function solve_system_matrix!(dv, system_matrix, rhs, + ::SerreGreenNaghdiEquations1D, + D1, cache) + if issparse(system_matrix) + (; factorization) = cache + cholesky!(factorization, system_matrix; check = false) + if issuccess(factorization) + scale_by_mass_matrix!(rhs, D1) + dv .= factorization \ rhs + else + # The factorization may fail if the time step is too large + # and h becomes negative. + fill!(dv, NaN) + end + else + factorization = cholesky!(system_matrix) + scale_by_mass_matrix!(rhs, D1) + ldiv!(dv, factorization, rhs) + end + return nothing +end + # Discretization that conserves # - the total water mass (integral of h) as a linear invariant # - the total momentum (integral of h v) as a nonlinear invariant @@ -163,58 +388,61 @@ end # equations in standard and hyperbolic form # [arXiv: 2408.02665](https://arxiv.org/abs/2408.02665) # TODO: Implement source terms -# TODO: Implement variable bathymetry function rhs!(dq, q, t, mesh, equations::SerreGreenNaghdiEquations1D, initial_condition, ::BoundaryConditionPeriodic, source_terms::Nothing, solver, cache) - if cache.D isa PeriodicUpwindOperators - rhs_sgn_upwind!(dq, q, equations, source_terms, cache, equations.bathymetry) + if cache.D1 isa PeriodicUpwindOperators + rhs_sgn_upwind!(dq, q, equations, source_terms, cache, equations.bathymetry_type) else - rhs_sgn_central!(dq, q, equations, source_terms, cache, equations.bathymetry) + rhs_sgn_central!(dq, q, equations, source_terms, cache, equations.bathymetry_type) end return nothing end function rhs_sgn_central!(dq, q, equations, source_terms, cache, ::BathymetryFlat) - # Unpack physical parameters and SBP operator `D` as well as the - # SBP operator in sparse matrix form `Dmat` + # Unpack physical parameters and SBP operator `D1` as well as the + # SBP operator in sparse matrix form `D1mat` g = gravity_constant(equations) - (; D, Dmat) = cache + (; D1, D1mat) = cache # `q` and `dq` are `ArrayPartition`s. They collect the individual # arrays for the water height `h` and the velocity `v`. - h, v = q.x - dh, dv = dq.x + eta, v, D = q.x + dh, dv, dD = dq.x # dh = deta since b is constant in time + fill!(dD, zero(eltype(dD))) @trixi_timeit timer() "hyperbolic terms" begin # Compute all derivatives required below - (; h_x, v_x, h2_x, hv_x, v2_x, h2_v_vx_x, + (; h, b, h_x, v_x, h2_x, hv_x, v2_x, h2_v_vx_x, h_vx_x, p_x, tmp, M_h, M_h3_3) = cache - mul!(h_x, D, h) - mul!(v_x, D, v) + @. b = equations.eta0 - D + @. h = eta - b + + mul!(h_x, D1, h) + mul!(v_x, D1, v) @. tmp = h^2 - mul!(h2_x, D, tmp) + mul!(h2_x, D1, tmp) @. tmp = h * v - mul!(hv_x, D, tmp) + mul!(hv_x, D1, tmp) @. tmp = v^2 - mul!(v2_x, D, tmp) + mul!(v2_x, D1, tmp) @. tmp = h^2 * v * v_x - mul!(h2_v_vx_x, D, tmp) + mul!(h2_v_vx_x, D1, tmp) @. tmp = h * v_x - mul!(h_vx_x, D, tmp) + mul!(h_vx_x, D1, tmp) inv6 = 1 / 6 @. tmp = (0.5 * h^2 * (h * v_x + h_x * v) * v_x - inv6 * h * h2_v_vx_x - inv6 * h^2 * v * h_vx_x) - mul!(p_x, D, tmp) + mul!(p_x, D1, tmp) # Plain: h_t + (h v)_x = 0 # @@ -238,87 +466,69 @@ function rhs_sgn_central!(dq, q, equations, source_terms, cache, ::BathymetryFla p_x) end + # The code below is equivalent to + # dv .= (Diagonal(h) - D1mat * Diagonal(1/3 .* h.^3) * D1mat) \ tmp + # but faster since the symbolic factorization is reused. @trixi_timeit timer() "assembling elliptic operator" begin - # The code below is equivalent to - # dv .= (Diagonal(h) - Dmat * Diagonal(1/3 .* h.^3) * Dmat) \ tmp - # but faster since the symbolic factorization is reused. - # Floating point errors accumulate a bit and the system matrix - # is not necessarily perfectly symmetric but only up to round-off - # errors. We wrap it here to avoid issues with the factorization. - @. M_h = h - scale_by_mass_matrix!(M_h, D) - inv3 = 1 / 3 - @. M_h3_3 = inv3 * h^3 - scale_by_mass_matrix!(M_h3_3, D) - system_matrix = Symmetric(Diagonal(M_h) - + - Dmat' * Diagonal(M_h3_3) * Dmat) + system_matrix = assemble_system_matrix!(cache, h, + D1, D1mat, + equations) end @trixi_timeit timer() "solving elliptic system" begin - if issparse(system_matrix) - (; factorization) = cache - cholesky!(factorization, system_matrix; check = false) - if issuccess(factorization) - scale_by_mass_matrix!(tmp, D) - dv .= factorization \ tmp - else - # The factorization may fail if the time step is too large - # and h becomes negative. - fill!(dv, NaN) - end - else - factorization = cholesky!(system_matrix) - scale_by_mass_matrix!(tmp, D) - ldiv!(dv, factorization, tmp) - end + solve_system_matrix!(dv, system_matrix, tmp, + equations, D1, cache) end return nothing end function rhs_sgn_upwind!(dq, q, equations, source_terms, cache, ::BathymetryFlat) - # Unpack physical parameters and SBP operator `D` as well as the - # SBP upwind operator in sparse matrix form `Dmat_minus` + # Unpack physical parameters and SBP operator `D1` as well as the + # SBP upwind operator in sparse matrix form `D1mat_minus` g = gravity_constant(equations) - (; Dmat_minus) = cache - D_upwind = cache.D - D = D_upwind.central + (; D1mat_minus) = cache + D1_upwind = cache.D1 + D1 = D1_upwind.central # `q` and `dq` are `ArrayPartition`s. They collect the individual # arrays for the water height `h` and the velocity `v`. - h, v = q.x - dh, dv = dq.x + eta, v, D = q.x + dh, dv, dD = dq.x # dh = deta since b is constant in time + fill!(dD, zero(eltype(dD))) @trixi_timeit timer() "hyperbolic terms" begin # Compute all derivatives required below - (; h_x, v_x, v_x_upwind, h2_x, hv_x, v2_x, + (; h, b, h_x, v_x, v_x_upwind, h2_x, hv_x, v2_x, h2_v_vx_x, h_vx_x, p_x, tmp, M_h, M_h3_3) = cache - mul!(h_x, D, h) - mul!(v_x, D, v) - mul!(v_x_upwind, D_upwind.minus, v) + @. b = equations.eta0 - D + @. h = eta - b + + mul!(h_x, D1, h) + mul!(v_x, D1, v) + mul!(v_x_upwind, D1_upwind.minus, v) @. tmp = h^2 - mul!(h2_x, D, tmp) + mul!(h2_x, D1, tmp) @. tmp = h * v - mul!(hv_x, D, tmp) + mul!(hv_x, D1, tmp) @. tmp = v^2 - mul!(v2_x, D, tmp) + mul!(v2_x, D1, tmp) @. tmp = h^2 * v * v_x - mul!(h2_v_vx_x, D, tmp) + mul!(h2_v_vx_x, D1, tmp) @. tmp = h * v_x - mul!(h_vx_x, D, tmp) + mul!(h_vx_x, D1, tmp) # p_+ @. tmp = 0.5 * h^2 * (h * v_x + h_x * v) * v_x_upwind - mul!(p_x, D_upwind.plus, tmp) + mul!(p_x, D1_upwind.plus, tmp) # p_0 minv6 = -1 / 6 @. tmp = minv6 * (h * h2_v_vx_x + h^2 * v * h_vx_x) - mul!(p_x, D, tmp, 1.0, 1.0) + mul!(p_x, D1, tmp, 1.0, 1.0) # Plain: h_t + (h v)_x = 0 # @@ -342,36 +552,240 @@ function rhs_sgn_upwind!(dq, q, equations, source_terms, cache, ::BathymetryFlat p_x) end + # The code below is equivalent to + # dv .= (Diagonal(h) - D1mat_plus * Diagonal(1/3 .* h.^3) * D1mat_minus) \ tmp + # but faster since the symbolic factorization is reused. @trixi_timeit timer() "assembling elliptic operator" begin - # The code below is equivalent to - # dv .= (Diagonal(h) - Dmat_plus * Diagonal(1/3 .* h.^3) * Dmat_minus) \ tmp - # but faster since the symbolic factorization is reused. - # Floating point errors accumulate a bit and the system matrix - # is not necessarily perfectly symmetric but only up to round-off errors. - # We wrap it here to avoid issues with the factorization. - @. M_h = h - scale_by_mass_matrix!(M_h, D) - inv3 = 1 / 3 - @. M_h3_3 = inv3 * h^3 - scale_by_mass_matrix!(M_h3_3, D) - system_matrix = Symmetric(Diagonal(M_h) - + - Dmat_minus' * Diagonal(M_h3_3) * Dmat_minus) + system_matrix = assemble_system_matrix!(cache, h, + D1, D1mat_minus, + equations) end @trixi_timeit timer() "solving elliptic system" begin - (; factorization) = cache - cholesky!(factorization, system_matrix; check = false) - if issuccess(factorization) - scale_by_mass_matrix!(tmp, D) - dv .= factorization \ tmp - else - # The factorization may fail if the time step is too large - # and h becomes negative. - fill!(dv, NaN) + solve_system_matrix!(dv, system_matrix, tmp, + equations, D1, cache) + end + + return nothing +end + +function rhs_sgn_central!(dq, q, equations, source_terms, cache, + ::Union{BathymetryMildSlope, BathymetryVariable}) + # Unpack physical parameters and SBP operator `D1` as well as the + # SBP operator in sparse matrix form `D1mat` + g = gravity_constant(equations) + (; D1, D1mat) = cache + + # `q` and `dq` are `ArrayPartition`s. They collect the individual + # arrays for the water height `h` and the velocity `v`. + eta, v, D = q.x + dh, dv, dD = dq.x # dh = deta since b is constant in time + fill!(dD, zero(eltype(dD))) + + @trixi_timeit timer() "hyperbolic terms" begin + # Compute all derivatives required below + (; h, h_x, v_x, h_hpb_x, b, b_x, hv_x, v2_x, + h2_v_vx_x, h_vx_x, p_h, p_x, tmp, + M_h_p_h_bx2, M_h3_3, M_h2_bx) = cache + if equations.bathymetry_type isa BathymetryVariable + (; psi) = cache + end + + @. b = equations.eta0 - D + @. h = eta - b + mul!(b_x, D1, b) + + mul!(h_x, D1, h) + mul!(v_x, D1, v) + @. tmp = h * eta + mul!(h_hpb_x, D1, tmp) + @. tmp = h * v + mul!(hv_x, D1, tmp) + @. tmp = v^2 + mul!(v2_x, D1, tmp) + + @. tmp = h^2 * v * v_x + mul!(h2_v_vx_x, D1, tmp) + @. tmp = h * v_x + mul!(h_vx_x, D1, tmp) + inv6 = 1 / 6 + @. p_h = (0.5 * h * (h * v_x + h_x * v) * v_x + - + inv6 * h2_v_vx_x + - + inv6 * h * v * h_vx_x) + @. tmp = h * b_x * v^2 + mul!(p_x, D1, tmp) + @. p_h += 0.25 * p_x + if equations.bathymetry_type isa BathymetryVariable + @. psi = 0.125 * p_x + end + @. tmp = b_x * v + mul!(p_x, D1, tmp) + @. p_h += 0.25 * h * v * p_x + if equations.bathymetry_type isa BathymetryVariable + @. psi += 0.125 * h * v * p_x + end + @. p_h = p_h - 0.25 * (h_x * v + h * v_x) * b_x * v + if equations.bathymetry_type isa BathymetryVariable + @. psi -= 0.125 * (h_x * v + h * v_x) * b_x * v + end + @. tmp = p_h * h + mul!(p_x, D1, tmp) + + # Plain: h_t + (h v)_x = 0 + # + # Split form for energy conservation: + # h_t + h_x v + h v_x = 0 + @. dh = -(h_x * v + h * v_x) + + # Plain: h v_t + ... = 0 + # + # Split form for energy conservation: + @. tmp = -(g * h_hpb_x - g * eta * h_x + + + 0.5 * h * v2_x + - + 0.5 * v^2 * h_x + + + 0.5 * hv_x * v + - + 0.5 * h * v * v_x + + p_x + + 1.5 * p_h * b_x) + if equations.bathymetry_type isa BathymetryVariable + @. tmp = tmp - psi * b_x end end + # The code below is equivalent to + # dv .= (Diagonal(h .+ factor .* h .* b_x.^2) - D1mat * (Diagonal(1/3 .* h.^3) * D1mat - Diagonal(0.5 .* h.^2 .* b_x) * D1mat) \ tmp + # but faster since the symbolic factorization is reused. + @trixi_timeit timer() "assembling elliptic operator" begin + system_matrix = assemble_system_matrix!(cache, h, b_x, + D1, D1mat, + equations) + end + + @trixi_timeit timer() "solving elliptic system" begin + solve_system_matrix!(dv, system_matrix, tmp, + equations, D1, cache) + end + + return nothing +end + +function rhs_sgn_upwind!(dq, q, equations, source_terms, cache, + ::Union{BathymetryMildSlope, BathymetryVariable}) + # Unpack physical parameters and SBP operator `D1` as well as the + # SBP operator in sparse matrix form `D1mat` + g = gravity_constant(equations) + (; D1mat_minus) = cache + D1_upwind = cache.D1 + D1 = D1_upwind.central + + # `q` and `dq` are `ArrayPartition`s. They collect the individual + # arrays for the water height `h` and the velocity `v`. + eta, v, D = q.x + dh, dv, dD = dq.x # dh = deta since b is constant in time + fill!(dD, zero(eltype(dD))) + + @trixi_timeit timer() "hyperbolic terms" begin + # Compute all derivatives required below + (; h, h_x, v_x, v_x_upwind, h_hpb_x, b, b_x, hv_x, v2_x, + h2_v_vx_x, h_vx_x, p_h, p_0, p_x, tmp, + M_h_p_h_bx2, M_h3_3, M_h2_bx) = cache + if equations.bathymetry_type isa BathymetryVariable + (; psi) = cache + end + + @. b = equations.eta0 - D + @. h = eta - b + mul!(b_x, D1, b) + + mul!(h_x, D1, h) + mul!(v_x, D1, v) + mul!(v_x_upwind, D1_upwind.minus, v) + @. tmp = h * eta + mul!(h_hpb_x, D1, tmp) + @. tmp = h * v + mul!(hv_x, D1, tmp) + @. tmp = v^2 + mul!(v2_x, D1, tmp) + + @. tmp = h^2 * v * v_x + mul!(h2_v_vx_x, D1, tmp) + @. tmp = h * v_x + mul!(h_vx_x, D1, tmp) + # p_0 + minv6 = -1 / 6 + @. p_h = minv6 * (h2_v_vx_x + + + h * v * h_vx_x) + @. tmp = h * b_x * v^2 + mul!(p_x, D1, tmp) + @. p_h += 0.25 * p_x + if equations.bathymetry_type isa BathymetryVariable + @. psi = 0.125 * p_x + end + @. tmp = b_x * v + mul!(p_x, D1, tmp) + @. p_h += 0.25 * h * v * p_x + if equations.bathymetry_type isa BathymetryVariable + @. psi += 0.125 * h * v * p_x + end + @. p_0 = p_h * h + mul!(p_x, D1, p_0) + # p_+ + @. tmp = (0.5 * h * (h * v_x + h_x * v) * v_x_upwind + - + 0.25 * (h_x * v + h * v_x) * b_x * v) + if equations.bathymetry_type isa BathymetryVariable + @. psi -= 0.125 * (h_x * v + h * v_x) * b_x * v + end + @. p_h = p_h + tmp + @. tmp = tmp * h + mul!(p_x, D1_upwind.plus, tmp, 1.0, 1.0) + + # Plain: h_t + (h v)_x = 0 + # + # Split form for energy conservation: + # h_t + h_x v + h v_x = 0 + @. dh = -(h_x * v + h * v_x) + + # Plain: h v_t + ... = 0 + # + # Split form for energy conservation: + @. tmp = -(g * h_hpb_x - g * eta * h_x + + + 0.5 * h * v2_x + - + 0.5 * v^2 * h_x + + + 0.5 * hv_x * v + - + 0.5 * h * v * v_x + + p_x + + 1.5 * p_h * b_x) + if equations.bathymetry_type isa BathymetryVariable + @. tmp = tmp - psi * b_x + end + end + + # The code below is equivalent to + # dv .= (Diagonal(h .+ factor .* h .* b_x.^2) - D1mat * (Diagonal(1/3 .* h.^3) * D1mat - Diagonal(0.5 .* h.^2 .* b_x)) - Diagonal(0.5 .* h.^2 .* b_x) * D1mat) \ tmp + # but faster since the symbolic factorization is reused. + @trixi_timeit timer() "assembling elliptic operator" begin + system_matrix = assemble_system_matrix!(cache, h, b_x, + D1, D1mat_minus, + equations) + end + + @trixi_timeit timer() "solving elliptic system" begin + solve_system_matrix!(dv, system_matrix, tmp, + equations, D1, cache) + end + return nothing end @@ -406,10 +820,6 @@ end return equations.eta0 - D end -@inline function waterheight(q, equations::SerreGreenNaghdiEquations1D) - return waterheight_total(q, equations) - bathymetry(q, equations) -end - # The entropy/energy takes the whole `q` for every point in space """ energy_total_modified(q_global, equations::SerreGreenNaghdiEquations1D, cache) @@ -421,9 +831,17 @@ derivative compared to the usual [`energy_total`](@ref) modeling non-hydrostatic contributions. The [`energy_total_modified`](@ref) is a conserved quantity (for periodic boundary conditions). -For a [`bathymetry_flat`](@ref) the total energy is given by +For a [`bathymetry_flat`](@ref) the total modified energy is given by +```math +\\frac{1}{2} g \\eta^2 + \\frac{1}{2} h v^2 + \\frac{1}{6} h^3 v_x^2. +``` +For a [`bathymetry_mild_slope`](@ref) the total modified energy is given by ```math -\\frac{1}{2} g h^2 + \\frac{1}{2} h v^2 + \\frac{1}{6} h^3 v_x^2. +\\frac{1}{2} g \\eta^2 + \\frac{1}{2} h v^2 + \\frac{1}{6} h (-h v_x + 1.5 v b_x)^2. +``` +For a [`bathymetry_variable`](@ref) the total modified energy has the additional term +```math ++ \\frac{1}{8} h (v b_x)^2. ``` `q_global` is a vector of the primitive variables at ALL nodes. @@ -432,25 +850,46 @@ For a [`bathymetry_flat`](@ref) the total energy is given by function energy_total_modified(q_global, equations::SerreGreenNaghdiEquations1D, cache) - # unpack physical parameters and SBP operator `D` + # unpack physical parameters and SBP operator `D1` g = equations.gravity - (; D, v_x) = cache + (; D1, h, b, v_x) = cache # `q_global` is an `ArrayPartition`. It collects the individual arrays for - # the water height `h` and the velocity `v`. - h, v = q_global.x + # the total water height `eta = h + b` and the velocity `v`. + eta, v = q_global.x + let D = q_global.x[3] + @. b = equations.eta0 - D + @. h = eta - b + end N = length(v) e = zeros(eltype(q_global), N) - # 1/2 g h^2 + 1/2 h v^2 + 1/6 h^3 v_x^2 - if D isa PeriodicUpwindOperators - mul!(v_x, D.minus, v) + # 1/2 g eta^2 + 1/2 h v^2 + 1/6 h^3 w^2 + # and + 1/8 h (v b_x)^2 for full bathymetry without mild-slope approximation + if D1 isa PeriodicUpwindOperators + mul!(v_x, D1.minus, v) else - mul!(v_x, D, v) + mul!(v_x, D1, v) end - @. e = 1 / 2 * g * h^2 + 1 / 2 * h * v^2 + 1 / 6 * h^3 * v_x^2 + if equations.bathymetry_type isa BathymetryFlat + b_x = cache.tmp + fill!(b_x, zero(eltype(b_x))) + else + (; b, b_x) = cache + @. b = equations.eta0 - q_global.x[3] + if D1 isa PeriodicUpwindOperators + mul!(b_x, D1.central, b) + else + mul!(b_x, D1, b) + end + end + + @. e = 1 / 2 * g * eta^2 + 1 / 2 * h * v^2 + 1 / 6 * h * (-h * v_x + 1.5 * v * b_x)^2 + if equations.bathymetry_type isa BathymetryVariable + @. e += 1 / 8 * h * (v * b_x)^2 + end return e end diff --git a/src/equations/svaerd_kalisch_1d.jl b/src/equations/svaerd_kalisch_1d.jl index 0a5949ea..bd8f2793 100644 --- a/src/equations/svaerd_kalisch_1d.jl +++ b/src/equations/svaerd_kalisch_1d.jl @@ -326,10 +326,6 @@ end return equations.eta0 - D end -@inline function waterheight(q, equations::SvaerdKalischEquations1D) - return waterheight_total(q, equations) - bathymetry(q, equations) -end - @inline entropy(u, equations::SvaerdKalischEquations1D) = energy_total(u, equations) # The modified entropy/total energy takes the whole `q` for every point in space @@ -368,10 +364,3 @@ is a conserved quantity of the Svärd-Kalisch equations given by end return e_modified end - -# Calculate the error for the "lake-at-rest" test case where eta should -# be a constant value over time -@inline function lake_at_rest_error(u, equations::SvaerdKalischEquations1D) - eta, _, _ = u - return abs(equations.eta0 - eta) -end diff --git a/test/test_serre_green_naghdi_1d.jl b/test/test_serre_green_naghdi_1d.jl index c5bb1e55..d0679264 100644 --- a/test/test_serre_green_naghdi_1d.jl +++ b/test/test_serre_green_naghdi_1d.jl @@ -21,6 +21,36 @@ EXAMPLES_DIR = joinpath(examples_dir(), "serre_green_naghdi_1d") @test_allocations(semi, sol, 550_000) end + @trixi_testset "serre_green_naghdi_soliton.jl with bathymetry_mild_slope" begin + # same values as serre_green_naghdi_soliton.jl but with more allocations + @test_trixi_include(joinpath(EXAMPLES_DIR, + "serre_green_naghdi_soliton.jl"), + bathymetry_type=bathymetry_mild_slope, + tspan=(0.0, 0.1), + l2=[9.994998669268741e-7, 1.4703973445698635e-6, 0.0], + linf=[6.5496216650196e-7, 1.027617322124641e-6, 0.0], + cons_error=[0.0, 8.174581012099225e-10, 0.0], + change_waterheight=0.0, + change_entropy_modified=-3.1093350116861984e-11) + + @test_allocations(semi, sol, 800_000) + end + + @trixi_testset "serre_green_naghdi_soliton.jl with bathymetry_variable" begin + # same values as serre_green_naghdi_soliton.jl but with more allocations + @test_trixi_include(joinpath(EXAMPLES_DIR, + "serre_green_naghdi_soliton.jl"), + bathymetry_type=bathymetry_variable, + tspan=(0.0, 0.1), + l2=[9.994998669268741e-7, 1.4703973445698635e-6, 0.0], + linf=[6.5496216650196e-7, 1.027617322124641e-6, 0.0], + cons_error=[0.0, 8.174581012099225e-10, 0.0], + change_waterheight=0.0, + change_entropy_modified=-3.1093350116861984e-11) + + @test_allocations(semi, sol, 800_000) + end + @trixi_testset "serre_green_naghdi_soliton_fourier.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "serre_green_naghdi_soliton_fourier.jl"), @@ -34,6 +64,36 @@ EXAMPLES_DIR = joinpath(examples_dir(), "serre_green_naghdi_1d") @test_allocations(semi, sol, allocs=450_000) end + @trixi_testset "serre_green_naghdi_soliton_fourier.jl with bathymetry_mild_slope" begin + # same values as serre_green_naghdi_soliton_fourier.jl but with more allocations + @test_trixi_include(joinpath(EXAMPLES_DIR, + "serre_green_naghdi_soliton_fourier.jl"), + bathymetry_type=bathymetry_mild_slope, + tspan=(0.0, 0.1), + l2=[8.252225014546995e-8, 6.724994492548714e-7, 0.0], + linf=[2.672093302180656e-8, 9.642725156897014e-8, 0.0], + cons_error=[2.842170943040401e-14, 4.627409566637652e-13, 0.0], + change_waterheight=2.842170943040401e-14, + change_entropy_modified=-3.097966327914037e-11) + + @test_allocations(semi, sol, allocs=850_000) + end + + @trixi_testset "serre_green_naghdi_soliton_fourier.jl with bathymetry_variable" begin + # same values as serre_green_naghdi_soliton_fourier.jl but with more allocations + @test_trixi_include(joinpath(EXAMPLES_DIR, + "serre_green_naghdi_soliton_fourier.jl"), + bathymetry_type=bathymetry_variable, + tspan=(0.0, 0.1), + l2=[8.252225014546995e-8, 6.724994492548714e-7, 0.0], + linf=[2.672093302180656e-8, 9.642725156897014e-8, 0.0], + cons_error=[2.842170943040401e-14, 4.627409566637652e-13, 0.0], + change_waterheight=2.842170943040401e-14, + change_entropy_modified=-3.097966327914037e-11) + + @test_allocations(semi, sol, allocs=850_000) + end + @trixi_testset "serre_green_naghdi_soliton_upwind.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "serre_green_naghdi_soliton_upwind.jl"), @@ -47,6 +107,36 @@ EXAMPLES_DIR = joinpath(examples_dir(), "serre_green_naghdi_1d") @test_allocations(semi, sol, allocs=500_000) end + @trixi_testset "serre_green_naghdi_soliton_upwind.jl with bathymetry_mild_slope" begin + # same as serre_green_naghdi_soliton_upwind.jl but with more allocations + @test_trixi_include(joinpath(EXAMPLES_DIR, + "serre_green_naghdi_soliton_upwind.jl"), + tspan=(0.0, 0.1), + bathymetry_type=bathymetry_mild_slope, + l2=[1.4876412924488654e-6, 5.9888097605442856e-6, 0.0], + linf=[1.0863034516361836e-6, 4.105927902009476e-6, 0.0], + cons_error=[4.263256414560601e-14, 4.483030568991353e-8, 0.0], + change_waterheight=4.263256414560601e-14, + change_entropy_modified=-3.1036506698001176e-11) + + @test_allocations(semi, sol, allocs=750_000) + end + + @trixi_testset "serre_green_naghdi_soliton_upwind.jl with bathymetry_variable" begin + # same as serre_green_naghdi_soliton_upwind.jl but with more allocations + @test_trixi_include(joinpath(EXAMPLES_DIR, + "serre_green_naghdi_soliton_upwind.jl"), + tspan=(0.0, 0.1), + bathymetry_type=bathymetry_variable, + l2=[1.4876412924488654e-6, 5.9888097605442856e-6, 0.0], + linf=[1.0863034516361836e-6, 4.105927902009476e-6, 0.0], + cons_error=[4.263256414560601e-14, 4.483030568991353e-8, 0.0], + change_waterheight=4.263256414560601e-14, + change_entropy_modified=-3.1036506698001176e-11) + + @test_allocations(semi, sol, allocs=750_000) + end + @trixi_testset "serre_green_naghdi_soliton_relaxation.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "serre_green_naghdi_soliton_relaxation.jl"), @@ -59,6 +149,123 @@ EXAMPLES_DIR = joinpath(examples_dir(), "serre_green_naghdi_1d") @test_allocations(semi, sol, allocs=450_000) end + + @trixi_testset "serre_green_naghdi_soliton_relaxation.jl with bathymetry_mild_slope" begin + # same as serre_green_naghdi_soliton_relaxation.jl but with more allocations + @test_trixi_include(joinpath(EXAMPLES_DIR, + "serre_green_naghdi_soliton_relaxation.jl"), + tspan=(0.0, 0.1), + bathymetry_type=bathymetry_mild_slope, + l2=[8.252225169608892e-8, 6.724994488577288e-7, 0.0], + linf=[2.6716495016287922e-8, 9.642466235713909e-8, 0.0], + cons_error=[2.842170943040401e-14, 4.649614027130156e-13, 0.0], + change_waterheight=2.842170943040401e-14, + change_entropy_modified=0.0) + + @test_allocations(semi, sol, allocs=850_000) + end + + @trixi_testset "serre_green_naghdi_soliton_relaxation.jl with bathymetry_variable" begin + # same as serre_green_naghdi_soliton_relaxation.jl but with more allocations + @test_trixi_include(joinpath(EXAMPLES_DIR, + "serre_green_naghdi_soliton_relaxation.jl"), + tspan=(0.0, 0.1), + bathymetry_type=bathymetry_variable, + l2=[8.252225169608892e-8, 6.724994488577288e-7, 0.0], + linf=[2.6716495016287922e-8, 9.642466235713909e-8, 0.0], + cons_error=[2.842170943040401e-14, 4.649614027130156e-13, 0.0], + change_waterheight=2.842170943040401e-14, + change_entropy_modified=0.0) + + @test_allocations(semi, sol, allocs=850_000) + end + + @trixi_testset "serre_green_naghdi_well_balanced.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "serre_green_naghdi_well_balanced.jl"), + tspan=(0.0, 2.0), + l2=[0, 0, 0], + linf=[0, 0, 0], + cons_error=[0, 0, 0], + change_waterheight=0.0, + change_momentum=0.0, + change_entropy_modified=0.0, + lake_at_rest=0.0) + + @test_allocations(semi, sol, allocs=750_000) + end + + @trixi_testset "serre_green_naghdi_well_balanced.jl with bathymetry_mild_slope" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "serre_green_naghdi_well_balanced.jl"), + bathymetry_type=bathymetry_mild_slope, + tspan=(0.0, 2.0), + l2=[0, 0, 0], + linf=[0, 0, 0], + cons_error=[0, 0, 0], + change_waterheight=0.0, + change_momentum=0.0, + change_entropy_modified=0.0, + lake_at_rest=0.0) + + @test_allocations(semi, sol, allocs=750_000) + end + + @trixi_testset "serre_green_naghdi_dingemans.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "serre_green_naghdi_dingemans.jl"), + tspan=(0.0, 1.0), + l2=[0.24638385233253557, 0.8055038575889898, 0.0], + linf=[0.03627087716483823, 0.11899056073196576, 0.0], + cons_error=[5.684341886080802e-14, 3.57950443142954e-5, 0.0], + change_waterheight=5.684341886080802e-14, + change_entropy=2.4420619411102962e-5, + change_entropy_modified=-1.0404164640931413e-8) + + @test_allocations(semi, sol, allocs=750_000) + end + + @trixi_testset "serre_green_naghdi_dingemans.jl with bathymetry_mild_slope" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "serre_green_naghdi_dingemans.jl"), + tspan=(0.0, 1.0), + bathymetry_type=bathymetry_mild_slope, + l2=[0.24638385233253557, 0.8055038575889898, 0.0], + linf=[0.03627087716483823, 0.11899056073196576, 0.0], + cons_error=[5.684341886080802e-14, 3.579504431429478e-5, 0.0], + change_waterheight=5.684341886080802e-14, + change_entropy=2.4420619411102962e-5, + change_entropy_modified=-1.0404164640931413e-8) + + @test_allocations(semi, sol, allocs=750_000) + end + + @trixi_testset "serre_green_naghdi_conservation.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "serre_green_naghdi_conservation.jl"), + l2=[1.3655498085989206, 2.3967486930606716, 0.0], + linf=[1.001076318001934, 0.8052527556023067, 0.0], + cons_error=[0.0, 0.0002674927404067162, 0.0], + change_entropy=-0.0584189861183404, + change_entropy_modified=0.059273537492344985, + atol=1e-11, # to make CI pass + atol_ints=4e-9) # to make CI pass + + @test_allocations(semi, sol, allocs=900_000) + end + + @trixi_testset "serre_green_naghdi_conservation.jl with bathymetry_mild_slope" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "serre_green_naghdi_conservation.jl"), + bathymetry_type=bathymetry_mild_slope, + l2=[1.3655493671985637, 2.3967828251339003, 0.0], + linf=[1.001075913983051, 0.8052680970114169, 0.0], + cons_error=[1.1368683772161603e-13, 0.00026407261543415217, 0.0], + change_entropy=-0.058352284294869605, + change_entropy_modified=0.05927339747017868) + + @test_allocations(semi, sol, allocs=900_000) + end end end # module