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 <[email protected]>

* 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 <[email protected]>

* format

* Apply suggestions from code review

Co-authored-by: Joshua Lampert <[email protected]>

* Apply suggestions from code review

Co-authored-by: Joshua Lampert <[email protected]>

* Update src/equations/serre_green_naghdi_1d.jl

Co-authored-by: Joshua Lampert <[email protected]>

* 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 <[email protected]>

* Update examples/serre_green_naghdi_1d/serre_green_naghdi_well_balanced.jl

Co-authored-by: Joshua Lampert <[email protected]>

* 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 <[email protected]>

* Update src/equations/serre_green_naghdi_1d.jl

Co-authored-by: Joshua Lampert <[email protected]>

* fix docstring of energy_total_modified

* adapt tolerances for macOS CI

* not only integrals

---------

Co-authored-by: Joshua Lampert <[email protected]>
Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com>

.gitignore

@@ -8,3 +8,4 @@ out*/
 /docs/build/
 run
 run/*
+**/*.code-workspace

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

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)

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)

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

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

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

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

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)

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,

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)

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

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