Add linear dispersion relation

.gitignore

@@ -6,6 +6,10 @@
 **/Manifest.toml
 out*/
 /docs/build/
+docs/src/changelog.md
+docs/src/changelog_tmp.md
+docs/src/code_of_conduct.md
+docs/src/contributing.md
 run
 run/*
 **/*.code-workspace

docs/make.jl

@@ -81,6 +81,7 @@ makedocs(;
                                   size_threshold_warn = 1000 * 1024),
          pages = ["Home" => "index.md",
              "Overview" => "overview.md",
+             "Dispersion" => "dispersion.md",
              "Development" => "development.md",
              "Reference" => [
                  "TrixiBase" => "ref-trixibase.md",

docs/src/dispersion.md

@@ -0,0 +1,145 @@
+# Dispersive properties of the equations
+
+The equations implemented in DispersiveShallowWater.jl describe the propagation of waves in a shallow water system.
+The equations are dispersive, meaning that the phase speed of the waves depends on their wavelength. This is in
+contrast to non-dispersive waves, where all waves travel at the same speed, such as the solution to the shallow
+water equations. In this case the phase speed is constantly given by ``c_0 = \sqrt{g h_0}``, where ``g`` is the
+acceleration due to gravity and ``h_0`` is a reference water height.
+
+The linear dispersion relation of a wave equation describes the relationship between the angular frequency
+``\omega`` and the wavenumber ``k`` of a wave. The angular frequency is related to the period ``T`` of the wave
+by ``\omega = 2\pi / T`` and the wavenumber is related to the wavelength ``\lambda`` by ``k = 2\pi / \lambda``.
+The simple case of just one wave traveling with speed ``c`` is a harmonic (or plane wave) solution of the form
+``\eta(x, t) = \mathrm{e}^{\mathrm{i}(k x - \omega t)}`` with ``c = \omega / k``. The linear dispersion relation
+of a wave equation is then derived by substituting this solution into a linearized version of the equation and
+solving for ``\omega`` in terms of ``k``.
+
+Since all the equations implemented in DispersiveShallowWater.jl are approximations to the full wave equations
+given by the Euler equations, their dispersion relations do not describe the exact speed of waves, but only
+an approximation to it, where the approximation usually becomes better for longer wavelengths. The dispersion
+relation of the full wave system is given by
+
+```math
+\omega(k) = \pm \sqrt{g k \tanh(h_0 k)}.
+```
+
+In DispersiveShallowWater.jl, we can investigate the dispersion relations of the different equations. Let us
+plot the wave speeds of the different equations normalized by the shallow water wave speed ``c_0`` as a function
+of the wave number. We pick a reference water height of ``h_0 = 0.8`` and gravitational acceleration of ``g = 9.81``.
+
+```@example dispersion
+using DispersiveShallowWater
+using Plots
+
+h0 = eta0 = 0.8
+g = 9.81
+disp_rel = LinearDispersionRelation(h0)
+k = 0.01:0.01:5.0
+
+struct EulerEquations1D <: DispersiveShallowWater.AbstractEquations{1, 0}
+    gravity::Float64
+    eta0::Float64
+end
+function (disp_rel::LinearDispersionRelation)(equations::EulerEquations1D, k)
+    h0 = disp_rel.ref_height
+    g = gravity_constant(equations)
+    return sqrt(g * k * tanh(h0 * k))
+end
+
+euler = EulerEquations1D(g, eta0)
+c_euler = wave_speed.(disp_rel, euler, k; normalize = true)
+plot(k, c_euler, label = "Euler", xlabel = "k", ylabel = "c / c_0", legend = :topright)
+
+bbm = BBMEquation1D(; gravity_constant = g, eta0 = eta0)
+c_bbm = wave_speed.(disp_rel, bbm, k; normalize = true)
+plot!(k, c_bbm, label = "BBM")
+
+sk = SvaerdKalischEquations1D(; gravity_constant = g, eta0 = eta0)
+c_sk = wave_speed.(disp_rel, sk, k; normalize = true)
+plot!(k, c_sk, label = "Svärd-Kalisch")
+
+sgn = SerreGreenNaghdiEquations1D(; gravity_constant = g, eta0 = eta0)
+c_sgn = wave_speed.(disp_rel, sgn, k; normalize = true)
+plot!(k, c_sgn, label = "Serre-Green-Naghdi")
+
+savefig("dispersion_relations.png") # hide
+nothing # hide
+```
+
+![dispersion relations](dispersion_relations.png)
+
+To verify that the wave speed predicted by the dispersion relation is indeed the observed one, let us
+simulate a simple wave traveling in a domain with periodic boundary conditions. We initialize the
+wave as a traveling wave for the Euler equations by using the dispersion relation. As an example we
+solve the problem with the Svärd-Kalisch equations.
+
+```@example dispersion
+using OrdinaryDiffEqTsit5
+using Printf
+
+wave_number() = 3.0
+frequency(k) = disp_rel(euler, k)
+
+function initial_condition_traveling_wave(x, t, equations, mesh)
+    k = wave_number()
+    omega = frequency(k)
+    h0 = equations.eta0
+    A = 0.02
+    h = A * cos(k * x - omega * t)
+    v = sqrt(equations.gravity / k * tanh(k * h0)) * h / h0
+    eta = h + h0
+    D = equations.eta0
+    return SVector(eta, v, D)
+end
+
+k = wave_number()
+coordinates_min = 0
+coordinates_max = 10 * pi / k # five waves (wave length = 2pi/k)
+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, sk, initial_condition_traveling_wave, solver,
+                          boundary_conditions = boundary_condition_periodic)
+
+tspan = (0.0, 3.0)
+ode = semidiscretize(semi, tspan)
+
+saveat = range(tspan..., length = 201)
+sol = solve(ode, Tsit5(), abstol = 1e-7, reltol = 1e-7,
+            save_everystep = false, saveat = saveat, tstops = saveat)
+
+x = 0.2 * (coordinates_max - coordinates_min)
+c_euler = wave_speed(disp_rel, euler, k)
+c_sk = wave_speed(disp_rel, sk, k)
+anim = @animate for step in eachindex(sol.u)
+    t = sol.t[step]
+    x_t_euler = x + c_euler * t
+    x_t_sk = x + c_sk * t
+    index = argmin(abs.(DispersiveShallowWater.grid(semi) .- x_t_sk))
+    scatter([x_t_euler], [initial_condition_traveling_wave(x_t_euler, t, euler, mesh)],
+            color = :blue, label = nothing)
+    eta, _ = sol.u[step].x
+    scatter!([x_t_sk], [eta[index]],
+             color = :green, label = nothing)
+    plot!(semi => sol, plot_initial = true, plot_bathymetry = false,
+          conversion = waterheight_total, step = step, legend = :topleft, linewidth = 2,
+          plot_title = @sprintf("t = %.3f", t), yrange = (0.77, 0.83),
+          linestyles = [:solid :dot], labels = ["Euler" "Svärd-Kälisch"],
+          color = [:blue :green])
+end
+gif(anim, "traveling_waves.gif", fps = 25)
+
+nothing # hide
+```
+
+![traveling waves](traveling_waves.gif)
+
+The dots in the movie show that the wave speed predicted by the dispersion relation is indeed the same
+as the one obtained by numerically solving the equations.
+As expected from the plot above the wave speed of the Svärd-Kalisch equations is a bit faster than the correct
+one, which is due to the approximation made in the equations.

docs/src/ref.md

@@ -16,6 +16,13 @@ Modules = [DispersiveShallowWater]
 Pages = Main.EQUATIONS_FILES
 ```
 
+## Linear dispersion relations
+
+```@autodocs
+Modules = [DispersiveShallowWater]
+Pages = ["dispersion_relation.jl"]
+```
+
 ## Mesh
 
 ```@autodocs

src/DispersiveShallowWater.jl

@@ -52,6 +52,7 @@ using TrixiBase: TrixiBase, @trixi_timeit, timer
 include("boundary_conditions.jl")
 include("mesh.jl")
 include("equations/equations.jl")
+include("dispersion_relation.jl")
 include("solver.jl")
 include("semidiscretization.jl")
 include("callbacks_step/callbacks_step.jl")
@@ -65,6 +66,8 @@ export AbstractShallowWaterEquations,
        SvärdKalischEquations1D, SvaerdKalischEquations1D,
        SerreGreenNaghdiEquations1D, HyperbolicSerreGreenNaghdiEquations1D
 
+export LinearDispersionRelation, wave_speed
+
 export prim2prim, prim2cons, cons2prim, prim2phys,
        waterheight_total, waterheight,
        velocity, momentum, discharge,

src/dispersion_relation.jl

@@ -0,0 +1,99 @@
+@doc raw"""
+    LinearDispersionRelation(ref_height)
+
+A struct representing a linear dispersion relation ``\omega(k)`` of an equation. The reference
+water height `h0` is given by `ref_height`. A dispersion relation can be called as
+`disp_rel(equations, k)` to compute the wave frequency ``\omega(k)`` for a given wavenumber `k`
+and a set of equations.
+
+See also [`wave_speed`](@ref) for computing the wave speed ``c = \omega(k) / k``.
+"""
+struct LinearDispersionRelation{RealT <: Real}
+    ref_height::RealT
+end
+
+function Base.show(io::IO, disp_rel::LinearDispersionRelation)
+    print(io, "LinearDispersionRelation(h0 = ", disp_rel.ref_height, ")")
+end
+
+Base.broadcastable(disp_rel::LinearDispersionRelation) = Ref(disp_rel)
+
+@doc raw"""
+    wave_speed(disp_rel, equations, k; normalize = false)
+
+Compute the wave speed ``c`` for a given wavenumber ``k`` using the linear dispersion relation ``disp_rel``
+of the `equations`.
+The wave speed is given by ``c = \omega(k) / k``. If `normalize` is `true`, the wave speed is normalized
+by the shallow water wave speed ``\sqrt{g h0}``, where `g` is the gravity constant and `h0` is the reference
+water height of the dispersion relation.
+
+See also [`LinearDispersionRelation`](@ref).
+"""
+function wave_speed(disp_rel::LinearDispersionRelation, equations, k;
+                    normalize = false)
+    omega = disp_rel(equations, k)
+    c = omega / k
+    if normalize
+        c /= sqrt(gravity_constant(equations) * disp_rel.ref_height)
+    end
+    return c
+end
+
+# TODO: Factor of 5/2 in front?
+function (disp_rel::LinearDispersionRelation)(equations::BBMEquation1D, k)
+    h0 = disp_rel.ref_height
+    g = gravity_constant(equations)
+    return sqrt(g * h0) * k / (1 + 1 / 6 * (h0 * k)^2)
+end
+
+# See
+# - Joshua Lampert, Hendrik Ranocha (2024)
+#  Structure-Preserving Numerical Methods for Two Nonlinear Systems of Dispersive Wave Equations
+#  [DOI: 10.48550/arXiv.2402.16669](https://doi.org/10.48550/arXiv.2402.16669)
+function (disp_rel::LinearDispersionRelation)(equations::BBMBBMEquations1D, k)
+    h0 = disp_rel.ref_height
+    g = gravity_constant(equations)
+    return sqrt(g * h0) * k / (1 + 1 / 6 * (h0 * k)^2)
+end
+
+# See
+# - 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)
+function (disp_rel::LinearDispersionRelation)(equations::SvärdKalischEquations1D, k)
+    h0 = disp_rel.ref_height
+    g = gravity_constant(equations)
+    c0 = sqrt(g * h0)
+    alpha = equations.alpha * c0 * h0^2
+    beta = equations.beta * h0^3
+    gamma = equations.gamma * c0 * h0^3
+    a = (1 + beta / h0 * k^2)
+    b = (-alpha - beta * alpha / h0 * k^2 - gamma / h0) * k^3
+    c = -g * h0 * k^2 + gamma * alpha / h0 * k^6
+    return (-b + sqrt(b^2 - 4 * a * c)) / (2 * a)
+end
+
+# See, e.g., eq. (45) (α = 1) in
+# - Maria Kazolea, Andrea Filippini, Mario Ricchiuto (2023)
+#   Low dispersion finite volume/element discretization of the enhanced Green-Naghdi equations for
+#   wave propagation, breaking and runup on unstructured meshes
+#   [DOI: 10.1016/j.ocemod.2022.102157](https://doi.org/10.1016/j.ocemod.2022.102157)
+#
+# or eq. (53) (β = 0) in
+# - Didier Clamond Denys Dutykh, Dimitrios Mitsotakis (2017)
+#   Conservative modified Serre–Green–Naghdi equations with improved dispersion characteristics
+#   [DOI: 10.1016/j.cnsns.2016.10.009](https://doi.org/10.1016/j.cnsns.2016.10.009)
+function (disp_rel::LinearDispersionRelation)(equations::SerreGreenNaghdiEquations1D, k)
+    h0 = disp_rel.ref_height
+    g = gravity_constant(equations)
+    return sqrt(g * h0) * k / sqrt(1.0 + (k * h0)^2 / 3)
+end
+
+# TODO: Make this for general lambda (how to understand eq. (19) in Favrie and Gavrilyuk?)
+function (disp_rel::LinearDispersionRelation)(equations::HyperbolicSerreGreenNaghdiEquations1D,
+                                              k)
+    h0 = disp_rel.ref_height
+    g = gravity_constant(equations)
+    lambda = equations.lambda
+    return sqrt(g * h0) * k / sqrt(1.0 + (k * h0)^2 / 3)
+end

src/equations/equations.jl

@@ -20,6 +20,8 @@ for the variables in `equations`. In particular, not the variables themselves ar
 """
 @inline eachvariable(equations::AbstractEquations) = Base.OneTo(nvariables(equations))
 
+Base.broadcastable(equations::AbstractEquations) = (equations,)
+
 @doc raw"""
     AbstractShallowWaterEquations{NDIMS, NVARS}
 

src/visualization.jl

@@ -63,7 +63,7 @@ end
             @series begin
                 subplot --> i
                 linestyle := :solid
-                label := "initial $(names[i])"
+                label --> "initial $(names[i])"
                 grid(semi), data_exact[i, :]
             end
         end

test/test_unit.jl

@@ -403,6 +403,44 @@ end
     @test_nowarn display(summary_callback)
 end
 
+@testitem "LinearDispersionRelation" setup=[Setup] begin
+    disp_rel = LinearDispersionRelation(3.0)
+    @test_nowarn print(disp_rel)
+    @test_nowarn display(disp_rel)
+    g = 9.81
+    k = 2 * pi
+    frequencies = [
+        0.5660455316649682,
+        0.5660455316649682,
+        7.700912310929906,
+        3.1189522995345467,
+        3.1189522995345467
+    ]
+    wave_speeds = [
+        0.09008894437955965,
+        0.09008894437955965,
+        1.2256382606017253,
+        0.4963966757387569,
+        0.4963966757387569
+    ]
+
+    for (i, equations) in enumerate((BBMEquation1D(gravity_constant = g),
+                                     BBMBBMEquations1D(gravity_constant = g),
+                                     SvärdKalischEquations1D(gravity_constant = g),
+                                     SerreGreenNaghdiEquations1D(gravity_constant = g),
+                                     HyperbolicSerreGreenNaghdiEquations1D(gravity_constant = g,
+                                                                           lambda = 1.0)))
+        @test isapprox(disp_rel(equations, k), frequencies[i])
+        @test isapprox(wave_speed(disp_rel, equations, k), wave_speeds[i])
+        # Add test for correct broadcasting
+        @test isapprox(disp_rel.(equations, [k, k]), [frequencies[i], frequencies[i]])
+        @test isapprox(wave_speed.(disp_rel, equations, [k, k]),
+                       [wave_speeds[i], wave_speeds[i]])
+        # For the normalized wave speed we expect c(0) = 1. Use eps() to avoid division by zero in c = omega / k
+        @test isapprox(wave_speed(disp_rel, equations, eps(), normalize = true), 1.0)
+    end
+end
+
 @testitem "util" setup=[Setup] begin
     @test_nowarn get_examples()