implement Serre-Green-Naghdi equations for flat bathymetry

Project.toml

@@ -37,7 +37,7 @@ SciMLBase = "2.11"
 SimpleUnPack = "1.1"
 SparseArrays = "1"
 StaticArrays = "1"
-SummationByPartsOperators = "0.5.52"
+SummationByPartsOperators = "0.5.63"
 TimerOutputs = "0.5.7"
 TrixiBase = "0.1.3"
 julia = "1.9"

README.md

@@ -10,10 +10,11 @@
 [![DOI](https://zenodo.org/badge/635090135.svg)](https://zenodo.org/doi/10.5281/zenodo.10034636)
 
 **DispersiveShallowWater.jl** is a [Julia](https://julialang.org/) package that implements structure-preserving numerical methods for dispersive shallow water models.
-To date, it provides provably conservative, entropy-conserving and well-balanced numerical schemes for two dispersive shallow water models:
+To date, it provides provably conservative, entropy-conserving and well-balanced numerical schemes for some dispersive shallow water models:
 
 * the [BBM-BBM equations with varying bottom topography](https://iopscience.iop.org/article/10.1088/1361-6544/ac3c29),
-* the [dispersive shallow water model proposed by Magnus Svärd and Henrik Kalisch](https://arxiv.org/abs/2302.09924).
+* the [dispersive shallow water model proposed by Magnus Svärd and Henrik Kalisch](https://arxiv.org/abs/2302.09924),
+* the [Serre-Green-Naghdi equations](https://arxiv.org/abs/2408.02665).
 
 The semidiscretizations are based on summation by parts (SBP) operators, which are implemented in [SummationByPartsOperators.jl](https://github.com/ranocha/SummationByPartsOperators.jl/).
 In order to obtain fully discrete schemes, the time integration methods from [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl) are used to solve the resulting ordinary differential equations.

examples/serre_green_naghdi_1d/serre_green_naghdi_soliton.jl

@@ -0,0 +1,40 @@
+using OrdinaryDiffEq
+using DispersiveShallowWater
+
+###############################################################################
+# Semidiscretization of the Serre-Green-Naghdi equations
+
+equations = SerreGreenNaghdiEquations1D(gravity_constant = 9.81)
+
+# initial_condition_convergence_test needs periodic boundary conditions
+initial_condition = initial_condition_convergence_test
+boundary_conditions = boundary_condition_periodic
+
+# create homogeneous mesh
+coordinates_min = -50.0
+coordinates_max = 50.0
+N = 512
+mesh = Mesh1D(coordinates_min, coordinates_max, N)
+
+# create solver with periodic SBP operators of accuracy order 4
+accuracy_order = 4
+solver = Solver(mesh, accuracy_order)
+
+# 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, (xmax(mesh) - xmin(mesh)) / sqrt(1.2 * equations.gravity)) # one period
+ode = semidiscretize(semi, tspan)
+summary_callback = SummaryCallback()
+analysis_callback = AnalysisCallback(semi; interval = 100,
+                                     extra_analysis_errors = (:conservation_error,),
+                                     extra_analysis_integrals = (waterheight_total,
+                                                                 entropy_modified))
+callbacks = CallbackSet(analysis_callback, summary_callback)
+
+saveat = range(tspan..., length = 100)
+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_fourier.jl

@@ -0,0 +1,41 @@
+using OrdinaryDiffEq
+using DispersiveShallowWater
+using SummationByPartsOperators: fourier_derivative_operator
+
+###############################################################################
+# Semidiscretization of the Serre-Green-Naghdi equations
+
+equations = SerreGreenNaghdiEquations1D(gravity_constant = 9.81)
+
+# initial_condition_convergence_test needs periodic boundary conditions
+initial_condition = initial_condition_convergence_test
+boundary_conditions = boundary_condition_periodic
+
+# create homogeneous mesh
+coordinates_min = -50.0
+coordinates_max = 50.0
+N = 128
+mesh = Mesh1D(coordinates_min, coordinates_max, N)
+
+# create solver with Fourier pseudospectral collocation method
+D1 = fourier_derivative_operator(xmin(mesh), xmax(mesh), 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, (xmax(mesh) - xmin(mesh)) / sqrt(1.2 * equations.gravity)) # one period
+ode = semidiscretize(semi, tspan)
+summary_callback = SummaryCallback()
+analysis_callback = AnalysisCallback(semi; interval = 100,
+                                     extra_analysis_errors = (:conservation_error,),
+                                     extra_analysis_integrals = (waterheight_total,
+                                                                 entropy_modified))
+callbacks = CallbackSet(analysis_callback, summary_callback)
+
+saveat = range(tspan..., length = 100)
+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_relaxation.jl

@@ -0,0 +1,43 @@
+using OrdinaryDiffEq
+using DispersiveShallowWater
+using SummationByPartsOperators: fourier_derivative_operator
+
+###############################################################################
+# Semidiscretization of the Serre-Green-Naghdi equations
+
+equations = SerreGreenNaghdiEquations1D(gravity_constant = 9.81)
+
+# initial_condition_convergence_test needs periodic boundary conditions
+initial_condition = initial_condition_convergence_test
+boundary_conditions = boundary_condition_periodic
+
+# create homogeneous mesh
+coordinates_min = -50.0
+coordinates_max = 50.0
+N = 128
+mesh = Mesh1D(coordinates_min, coordinates_max, N)
+
+# create solver with Fourier pseudospectral collocation method
+D1 = fourier_derivative_operator(xmin(mesh), xmax(mesh), 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, (xmax(mesh) - xmin(mesh)) / sqrt(1.2 * equations.gravity)) # one period
+ode = semidiscretize(semi, tspan)
+summary_callback = SummaryCallback()
+analysis_callback = AnalysisCallback(semi; interval = 100,
+                                     extra_analysis_errors = (:conservation_error,),
+                                     extra_analysis_integrals = (waterheight_total,
+                                                                 entropy_modified))
+relaxation_callback = RelaxationCallback(invariant = entropy_modified)
+# Always put relaxation_callback before analysis_callback to guarantee conservation of the invariant
+callbacks = CallbackSet(relaxation_callback, analysis_callback, summary_callback)
+
+saveat = range(tspan..., length = 100)
+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_upwind.jl

@@ -0,0 +1,48 @@
+using OrdinaryDiffEq
+using DispersiveShallowWater
+using SummationByPartsOperators: upwind_operators, periodic_derivative_operator
+using SparseArrays: sparse
+
+###############################################################################
+# Semidiscretization of the Serre-Green-Naghdi equations
+
+equations = SerreGreenNaghdiEquations1D(gravity_constant = 9.81)
+
+# initial_condition_convergence_test needs periodic boundary conditions
+initial_condition = initial_condition_convergence_test
+boundary_conditions = boundary_condition_periodic
+
+# create homogeneous mesh
+coordinates_min = -50.0
+coordinates_max = 50.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)
+accuracy_order = 4
+D1 = upwind_operators(periodic_derivative_operator;
+                      derivative_order = 1,
+                      accuracy_order = accuracy_order - 1,
+                      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, (xmax(mesh) - xmin(mesh)) / sqrt(1.2 * equations.gravity)) # one period
+ode = semidiscretize(semi, tspan)
+summary_callback = SummaryCallback()
+analysis_callback = AnalysisCallback(semi; interval = 100,
+                                     extra_analysis_errors = (:conservation_error,),
+                                     extra_analysis_integrals = (waterheight_total,
+                                                                 entropy_modified))
+callbacks = CallbackSet(analysis_callback, summary_callback)
+
+saveat = range(tspan..., length = 100)
+sol = solve(ode, Tsit5(), abstol = 1e-7, reltol = 1e-7,
+            save_everystep = false, callback = callbacks, saveat = saveat)

src/DispersiveShallowWater.jl

@@ -17,14 +17,17 @@ import SciMLBase: u_modified!
 
 @reexport using StaticArrays: SVector
 using SimpleUnPack: @unpack
-using SparseArrays: sparse
-using SummationByPartsOperators: AbstractDerivativeOperator,
+using SparseArrays: sparse, issparse
+using SummationByPartsOperators: SummationByPartsOperators,
+                                 AbstractDerivativeOperator,
                                  PeriodicDerivativeOperator, PeriodicUpwindOperators,
                                  UniformPeriodicCoupledOperator,
                                  DerivativeOperator, UpwindOperators,
                                  UniformCoupledOperator,
+                                 FourierDerivativeOperator,
                                  periodic_derivative_operator,
-                                 derivative_order, integrate, mass_matrix
+                                 derivative_order, integrate, mass_matrix,
+                                 scale_by_mass_matrix!
 import SummationByPartsOperators: grid, xmin, xmax
 using TimerOutputs: TimerOutputs, print_timer, reset_timer!
 @reexport using TrixiBase: trixi_include
@@ -41,11 +44,17 @@ include("util.jl")
 
 export examples_dir, get_examples, default_example, convergence_test
 
-export BBMBBMEquations1D, BBMBBMVariableEquations1D, SvärdKalischEquations1D,
-       SvaerdKalischEquations1D
+export AbstractShallowWaterEquations,
+       BBMBBMEquations1D, BBMBBMVariableEquations1D,
+       SvärdKalischEquations1D, SvaerdKalischEquations1D,
+       SerreGreenNaghdiEquations1D
 
-export prim2prim, prim2cons, cons2prim, waterheight_total, waterheight,
-       velocity, momentum, discharge, energy_total, entropy, lake_at_rest_error,
+export prim2prim, prim2cons, cons2prim,
+       waterheight_total, waterheight,
+       velocity, momentum, discharge,
+       gravity_constant,
+       bathymetry, still_water_surface,
+       energy_total, entropy, lake_at_rest_error,
        energy_total_modified, entropy_modified
 
 export Mesh1D, xmin, xmax, nnodes
@@ -56,6 +65,8 @@ export Semidiscretization, semidiscretize, grid
 
 export boundary_condition_periodic, boundary_condition_reflecting
 
+export bathymetry_flat
+
 export initial_condition_convergence_test,
        initial_condition_manufactured, source_terms_manufactured,
        initial_condition_manufactured_reflecting, source_terms_manufactured_reflecting,

src/callbacks_step/analysis.jl

@@ -116,7 +116,7 @@ end
 
 """
     errors(analysis_callback)
-    
+
 Return the computed errors for each timestep as a named tuple.
 The shape of each entry is (nvariables, ntimesteps).
 """
@@ -330,7 +330,10 @@ function analyze(quantity, q, t, semi::Semidiscretization)
     integrate_quantity(q -> quantity(q, semi.equations), q, semi)
 end
 
-# modified entropy from Svärd-Kalisch equations need to take the whole vector `u` for every point in space
+# The modified entropy/energy of the Svärd-Kalisch and
+# Serre-Green-Naghdi equations
+# takes the whole `q` for every point in space since it requires
+# the derivative of the velocity `v_x`.
 function analyze(quantity::Union{typeof(energy_total_modified), typeof(entropy_modified)},
                  q, t, semi::Semidiscretization)
     integrate_quantity(quantity, q, semi)

src/equations/bbm_bbm_1d.jl

@@ -1,5 +1,5 @@
 @doc raw"""
-    BBMBBMEquations1D(gravity, D, eta0 = 0.0)
+    BBMBBMEquations1D(; gravity_constant, D = 1.0, eta0 = 0.0)
 
 BBM-BBM (Benjamin–Bona–Mahony) system in one spatial dimension. The equations are given by
 ```math
@@ -292,11 +292,4 @@ end
     return waterheight_total(q, equations) - bathymetry(q, equations)
 end
 
-@inline function energy_total(q, equations::BBMBBMEquations1D)
-    eta, v = q
-    D = still_waterdepth(q, equations)
-    e = 0.5 * (equations.gravity * eta^2 + (D + eta - equations.eta0) * v^2)
-    return e
-end
-
 @inline entropy(q, equations::BBMBBMEquations1D) = energy_total(q, equations)

src/equations/bbm_bbm_variable_bathymetry_1d.jl

@@ -1,5 +1,5 @@
 @doc raw"""
-    BBMBBMVariableEquations1D(gravity, eta0 = 1.0)
+    BBMBBMVariableEquations1D(; gravity_constant, eta0 = 0.0)
 
 BBM-BBM (Benjamin–Bona–Mahony) system in one spatial dimension with spatially varying bathymetry. The equations are given by
 ```math
@@ -378,12 +378,6 @@ end
     return waterheight_total(q, equations) - bathymetry(q, equations)
 end
 
-@inline function energy_total(q, equations::BBMBBMVariableEquations1D)
-    eta, v, D = q
-    e = 0.5 * (equations.gravity * eta^2 + (D + eta - equations.eta0) * v^2)
-    return e
-end
-
 @inline entropy(q, equations::BBMBBMVariableEquations1D) = energy_total(q, equations)
 
 # Calculate the error for the "lake-at-rest" test case where eta should