Implement hyperbolized Sainte-Marie equations developed by Escalante et al.

NEWS.md

@@ -5,6 +5,12 @@ DispersiveShallowWater.jl follows the interpretation of
 used in the Julia ecosystem. Notable changes will be documented in this file
 for human readability.
 
+## Changes in the v0.10 lifecycle
+
+#### Added
+
+- The `HyperbolicSainteMarieEquations1D` were added for different types of bathymetry ([#288]).
+
 ## Changes when updating to v0.10 from v0.9.x
 
 #### Changed

README.md

@@ -21,7 +21,8 @@ To date, it provides provably conservative, entropy-conserving, and well-balance
 * the [Benjamin-Bona-Mahony (BBM) equation, also known as regularized long-wave equation](https://doi.org/10.4208/cicp.OA-2020-0119),
 * the [BBM-BBM equations with varying bottom topography](https://doi.org/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 [Serre-Green-Naghdi equations in standard and hyperbolic form](https://doi.org/10.1002/num.70016).
+* the [Serre-Green-Naghdi equations in standard and hyperbolic form](https://doi.org/10.1002/num.70016),
+* the [Sainte-Marie equations in hyperbolic form](https://doi.org/10.1016/j.jcp.2019.05.035).
 
 The semidiscretizations are based on summation-by-parts (SBP) operators, which are implemented in [SummationByPartsOperators.jl](https://github.com/ranocha/SummationByPartsOperators.jl/).
 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.

docs/src/dingemans.md

@@ -47,6 +47,10 @@ hysgn = HyperbolicSerreGreenNaghdiEquations1D(bathymetry_type = bathymetry_mild_
                                               # for actual simulations a higher lambda (~500) is recommended
                                               # it is chosen so low to be able to see the difference between it
                                               # and the SGN equation.
+
+# Hyperbolic approximation of the Sainte-Marie equations
+hysm = HyperbolicSainteMarieEquations1D(bathymetry_type = bathymetry_mild_slope,
+                                        gravity = 9.81, eta0 = 0.8, h0 = 0.8)
 nothing # hide
 ```
 
@@ -106,6 +110,8 @@ semi_sgn = Semidiscretization(mesh, sgn, initial_condition, solver_upwind,
                               boundary_conditions = boundary_conditions)
 semi_hysgn = Semidiscretization(mesh, hysgn, initial_condition, solver_central,
                                 boundary_conditions = boundary_conditions)
+semi_hysm = Semidiscretization(mesh, hysm, initial_condition, solver_central,
+                               boundary_conditions = boundary_conditions)
 nothing # hide
 ```
 
@@ -116,12 +122,13 @@ ode_bbmbbm = semidiscretize(semi_bbmbbm, tspan)
 ode_sk = semidiscretize(semi_sk, tspan)
 ode_sgn = semidiscretize(semi_sgn, tspan)
 ode_hysgn = semidiscretize(semi_hysgn, tspan)
-
+ode_hysm = semidiscretize(semi_hysm, tspan)
 options = (; abstol = 1e-7, reltol = 1e-7, save_everystep = false, saveat = saveat)
 sol_bbmbbm = solve(ode_bbmbbm, Tsit5(); options...)
 sol_sk = solve(ode_sk, Tsit5(); options...)
 sol_sgn = solve(ode_sgn, Tsit5(); options...)
 sol_hysgn = solve(ode_hysgn, Tsit5(); options...)
+sol_hysm = solve(ode_hysm, Tsit5(); options...)
 nothing # hide
 ```
 
@@ -147,7 +154,8 @@ models = [
     (semi_bbmbbm, sol_bbmbbm, "BBM-BBM", shifted_waterheight, :solid),
     (semi_sk, sol_sk, "Svärd-Kalisch", waterheight_total, :dashdotdot),
     (semi_sgn, sol_sgn, "Serre-Green-Naghdi", waterheight_total, :dot),
-    (semi_hysgn, sol_hysgn, "Hyperbolic Serre-Green-Naghdi", waterheight_total, :dashdot)
+    (semi_hysgn, sol_hysgn, "Hyperbolic Serre-Green-Naghdi", waterheight_total, :dashdot),
+    (semi_hysm, sol_hysm, "Hyperbolic Sainte-Marie", waterheight_total, :dash),
 ]
 
 # Create snapshot plots for each time
@@ -280,6 +288,9 @@ hysgn = HyperbolicSerreGreenNaghdiEquations1D(bathymetry_type = bathymetry_mild_
                                               # it is chosen so low to be able to see the difference between it
                                               # and the SGN equation.
 
+hysm = HyperbolicSainteMarieEquations1D(bathymetry_type = bathymetry_mild_slope,
+                                        gravity = 9.81, eta0 = 0.8, h0 = 0.8)
+
 initial_condition = initial_condition_dingemans
 boundary_conditions = boundary_condition_periodic
 
@@ -310,17 +321,21 @@ semi_sgn = Semidiscretization(mesh, sgn, initial_condition, solver_upwind,
                               boundary_conditions = boundary_conditions)
 semi_hysgn = Semidiscretization(mesh, hysgn, initial_condition, solver_central,
                                 boundary_conditions = boundary_conditions)
+semi_hysm = Semidiscretization(mesh, hysm, initial_condition, solver_central,
+                               boundary_conditions = boundary_conditions)
 
 ode_bbmbbm = semidiscretize(semi_bbmbbm, tspan)
 ode_sk = semidiscretize(semi_sk, tspan)
 ode_sgn = semidiscretize(semi_sgn, tspan)
 ode_hysgn = semidiscretize(semi_hysgn, tspan)
+ode_hysm = semidiscretize(semi_hysm, tspan)
 
 options = (; abstol = 1e-7, reltol = 1e-7, save_everystep = false, saveat = saveat)
 sol_bbmbbm = solve(ode_bbmbbm, Tsit5(); options...)
 sol_sk = solve(ode_sk, Tsit5(); options...)
 sol_sgn = solve(ode_sgn, Tsit5(); options...)
 sol_hysgn = solve(ode_hysgn, Tsit5(); options...)
+sol_hysm = solve(ode_hysm, Tsit5(); options...)
 
 # BBM-BBM equations need to be translated vertically for comparison
 shifted_waterheight(q, equations) = waterheight_total(q, equations) + 0.8
@@ -335,7 +350,8 @@ models = [
     (semi_bbmbbm, sol_bbmbbm, "BBM-BBM", shifted_waterheight, :solid),
     (semi_sk, sol_sk, "Svärd-Kalisch", waterheight_total, :dashdotdot),
     (semi_sgn, sol_sgn, "Serre-Green-Naghdi", waterheight_total, :dot),
-    (semi_hysgn, sol_hysgn, "Hyperbolic Serre-Green-Naghdi", waterheight_total, :dashdot)
+    (semi_hysgn, sol_hysgn, "Hyperbolic Serre-Green-Naghdi", waterheight_total, :dashdot),
+    (semi_hysm, sol_hysm, "Hyperbolic Sainte-Marie", waterheight_total, :dash),
 ]
 
 # Create snapshot plots for each time

docs/src/dispersion.md

@@ -48,21 +48,29 @@ plot(k, c_euler, label = "Euler", xlabel = "k", ylabel = "c / c_0", legend = :to
 
 kdv = KdVEquation1D(; gravity = g, eta0 = eta0, D = h0)
 c_kdv = wave_speed.(disp_rel, kdv, k; normalize = true)
-plot!(k, c_kdv, label = "KdV")
+plot!(k, c_kdv, label = "KdV", linestyle = :auto)
 
 bbm = BBMEquation1D(; gravity = g, eta0 = eta0, D = h0)
 c_bbm = wave_speed.(disp_rel, bbm, k; normalize = true)
-plot!(k, c_bbm, label = "BBM")
+plot!(k, c_bbm, label = "BBM", linestyle = :auto)
 
 # Optimized set 4 in the preprint
 sk = SvaerdKalischEquations1D(; gravity = g, eta0 = eta0,
                               alpha = 0.0, beta = 0.2308939393939394, gamma = 0.04034343434343434)
 c_sk = wave_speed.(disp_rel, sk, k; normalize = true)
-plot!(k, c_sk, label = "Svärd-Kalisch")
+plot!(k, c_sk, label = "Svärd-Kalisch", linestyle = :auto)
 
 sgn = SerreGreenNaghdiEquations1D(; gravity = g, eta0 = eta0)
 c_sgn = wave_speed.(disp_rel, sgn, k; normalize = true)
-plot!(k, c_sgn, label = "Serre-Green-Naghdi")
+plot!(k, c_sgn, label = "Serre-Green-Naghdi", linestyle = :auto)
+
+sm = SainteMarieEquations1D(; gravity = g, eta0 = eta0)
+c_sm = wave_speed.(disp_rel, sm, k; normalize = true)
+plot!(k, c_sm, label = "Sainte-Marie", linestyle = :auto)
+
+hsm = HyperbolicSainteMarieEquations1D(; gravity = g, eta0 = eta0, h0 = h0, alpha = 3.0)
+c_hsm = wave_speed.(disp_rel, hsm, k; normalize = true)
+plot!(k, c_hsm, label = "Hyperbolic Sainte-Marie (\$\\alpha = 3\$)", linestyle = :auto)
 
 savefig("dispersion_relations.png") # hide
 nothing # hide

docs/src/index.md

@@ -23,7 +23,8 @@ To date, it provides provably conservative, entropy-conserving, and well-balance
 * the [Benjamin-Bona-Mahony (BBM) equation, also known as regularized long-wave equation](https://doi.org/10.4208/cicp.OA-2020-0119),
 * the [BBM-BBM equations with varying bottom topography](https://doi.org/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 [Serre-Green-Naghdi equations in standard and hyperbolic form](https://doi.org/10.1002/num.70016).
+* the [Serre-Green-Naghdi equations in standard and hyperbolic form](https://doi.org/10.1002/num.70016),
+* the [Sainte-Marie equations in hyperbolic form](https://doi.org/10.1016/j.jcp.2019.05.035).
 
 The semidiscretizations are based on summation-by-parts (SBP) operators, which are implemented in [SummationByPartsOperators.jl](https://github.com/ranocha/SummationByPartsOperators.jl/).
 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.

docs/src/overview.md

@@ -1,6 +1,6 @@
 # [Equations](@id equations)
 
-[DispersiveShallowWater.jl](https://github.com/NumericalMathematics/DispersiveShallowWater.jl) provides six different dispersive shallow water equation systems for modeling water waves. Each equation system offers different levels of physical accuracy, computational complexity, and supports various boundary conditions and bathymetry types.
+[DispersiveShallowWater.jl](https://github.com/NumericalMathematics/DispersiveShallowWater.jl) provides seven different dispersive shallow water equation systems for modeling water waves. Each equation system offers different levels of physical accuracy, computational complexity, and supports various boundary conditions and bathymetry types.
 
 ## [Supported Models and Features](@id eq_overview)
 
@@ -14,6 +14,7 @@ The following table provides an overview of all available equation systems and t
 | [`Svärd-Kalisch`](@ref SvaerdKalischEquations1D) | ``(\eta, v, D)`` | ✅ | ✅ᵃ | ❌ | ❌ | ✅ | ✅ | ✅ |
 | [`Serre-Green-Naghdi`](@ref SerreGreenNaghdiEquations1D) | ``(\eta, v, D)`` | ✅ | ✅ | ✅ | ✅ | ✅ | ✅ | ✅ |
 | [`Hyperbolic SGN`](@ref HyperbolicSerreGreenNaghdiEquations1D) |``(\eta, v, D, w, H)`` | ✅ | ✅ | ✅ | ✅ | ❌ | ✅ | ✅ |
+| [`Hyperbolic Sainte-Marie`](@ref HyperbolicSainteMarieEquations1D) |``(\eta, v, D, w, p)`` | ✅ | ❌ | ✅ | ✅ | ❌ | ✅ | ✅ |
 
 *ᵃReflecting boundary conditions for Svärd-Kalisch equations require `alpha = gamma = 0`*
 
@@ -26,10 +27,11 @@ The following table provides an overview of all available equation systems and t
 - ``D``: Still-water depth
 - ``w``: Auxiliary variable in hyperbolic approximation (``\approx -h v_x``)
 - ``H``: Auxiliary variable in hyperbolic approximation (``\approx h``)
+- ``p``: Non-hydrostatic pressure in hyperbolic Sainte-Marie approximation
 
 ## Abstract Shallow Water Equations Interface
 
-Several equation systems in [DispersiveShallowWater.jl](https://github.com/NumericalMathematics/DispersiveShallowWater.jl) (`BBMBBMEquations1D`, `SvaerdKalischEquations1D`, `SerreGreenNaghdiEquations1D`, and `HyperbolicSerreGreenNaghdiEquations1D`) are subtypes of [`AbstractShallowWaterEquations`](@ref). This design reflects that these systems all contain the classical shallow water equations as a subsystem, extended with additional dispersive terms.
+Several equation systems in [DispersiveShallowWater.jl](https://github.com/NumericalMathematics/DispersiveShallowWater.jl) ([`BBMBBMEquations1D`](@ref), [`SvaerdKalischEquations1D`](@ref), [`SerreGreenNaghdiEquations1D`](@ref), [`HyperbolicSerreGreenNaghdiEquations1D`](@ref), and [`HyperbolicSainteMarieEquations1D`](@ref)) are subtypes of [`AbstractShallowWaterEquations`](@ref). This design reflects that these systems all contain the classical shallow water equations as a subsystem, extended with additional dispersive terms.
 
 The common interface provides shared functionality like [`waterheight`](@ref), [`velocity`](@ref), [`energy_total`](@ref), and [`entropy`](@ref). This enables consistent analysis and visualization across different dispersive models while maintaining the underlying shallow water physics.
 
@@ -47,9 +49,9 @@ The common interface provides shared functionality like [`waterheight`](@ref), [
 
 ## Hyperbolic Approximations
 
-Some equations are hyperbolic approximations of other systems (e.g., [`HyperbolicSerreGreenNaghdiEquations1D`](@ref) approximates [`SerreGreenNaghdiEquations1D`](@ref)). These systems support two approaches for initial conditions:
+Some equations are hyperbolic approximations of other systems (e.g., [`HyperbolicSerreGreenNaghdiEquations1D`](@ref) approximates [`SerreGreenNaghdiEquations1D`](@ref), and [`HyperbolicSainteMarieEquations1D`](@ref) approximates [`SainteMarieEquations1D`](@ref)). These systems support two approaches for initial conditions:
 
-1. **Full variables**: Specify all primitive variables including auxiliary variables ``(\eta, v, D, w, H)``
+1. **Full variables**: Specify all primitive variables including auxiliary variables (e.g., ``(\eta, v, D, w, H)`` or ``(\eta, v, D, w, p)``)
 2. **Physical variables**: Specify only the physical variables ``(\eta, v, D)`` from the limit system. The auxiliary variables are automatically initialized using appropriate approximations.
 
 This flexibility allows using the same initial conditions for both the original and hyperbolic systems, simplifying comparisons and testing.
@@ -94,3 +96,13 @@ DispersiveShallowWater.SerreGreenNaghdiEquations1D
 ```@docs
 DispersiveShallowWater.HyperbolicSerreGreenNaghdiEquations1D
 ```
+
+## Sainte-Marie and Hyperbolic Sainte-Marie
+
+```@docs
+DispersiveShallowWater.SainteMarieEquations1D
+```
+
+```@docs
+DispersiveShallowWater.HyperbolicSainteMarieEquations1D
+```

docs/src/ref.md

@@ -19,7 +19,9 @@ Filter = t -> !(t in [DispersiveShallowWater.KdVEquation1D,
                      DispersiveShallowWater.BBMBBMEquations1D,
                      DispersiveShallowWater.SvaerdKalischEquations1D,
                      DispersiveShallowWater.SerreGreenNaghdiEquations1D,
-                     DispersiveShallowWater.HyperbolicSerreGreenNaghdiEquations1D])
+                     DispersiveShallowWater.HyperbolicSerreGreenNaghdiEquations1D,
+                     DispersiveShallowWater.SainteMarieEquations1D,
+                     DispersiveShallowWater.HyperbolicSainteMarieEquations1D])
 ```
 
 ## Linear dispersion relations

examples/hyperbolic_sainte_marie_1d/hyperbolic_sainte_marie_conservation.jl

@@ -0,0 +1,92 @@
+# 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 OrdinaryDiffEqLowStorageRK
+using DispersiveShallowWater
+
+#=
+You can run this example manually and check the conservation properties as follows:
+
+julia> using Revise; using DispersiveShallowWater, OrdinaryDiffEqVerner
+
+julia> trixi_include("examples/hyperbolic_sainte_marie_1d/hyperbolic_sainte_marie_conservation.jl", tol = 1.0e-12, alg = Vern9());
+
+julia> integrals(analysis_callback).waterheight_total |> x -> (x[end] - x[1]) / x[1]
+0.0
+
+julia> integrals(analysis_callback).momentum |> x -> (x[end] - x[1]) / x[1]
+2.9432057444810533e-15 # for flat bathymetry
+
+julia> integrals(analysis_callback).entropy_modified |> x -> (x[end] - x[1]) / x[1]
+-1.064590038909088e-15
+=#
+
+###############################################################################
+# Semidiscretization of the hyperbolic Sainte-Marie equations
+
+equations = HyperbolicSainteMarieEquations1D(bathymetry_type = bathymetry_mild_slope,
+                                             gravity = 9.81,
+                                             h0 = 1.0)
+
+function initial_condition_conservation_test(x, t,
+                                             equations::HyperbolicSainteMarieEquations1D,
+                                             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)
+
+    # We use the feature that we can only return the physical variables
+    # used by the `hyperbolic_approximation_limit`, i.e., the
+    # `SainteMarieEquations1D.`
+    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 = 50, io,
+                                     extra_analysis_errors = (:conservation_error,),
+                                     extra_analysis_integrals = (waterheight_total,
+                                                                 momentum,
+                                                                 entropy,
+                                                                 entropy_modified))
+
+callbacks = CallbackSet(analysis_callback, summary_callback)
+
+# optimized time integration methods like this one are much more efficient
+# for stiff problems (α big) than standard methods like Tsit5()
+alg = RDPK3SpFSAL35()
+tol = 1.0e-4
+sol = solve(ode, alg;
+            reltol = tol, abstol = tol,
+            save_everystep = false, callback = callbacks)