Add linear dispersion relation

[JoshuaLampert]

This adds a new feature to compute the dispersion relation and wave speed of the different equations. It also adds a new page in the docs explaining the general concept of dispersion, which uses the new LinearDispersionRelation to compare the dispersive behavior of the different equations.
I have three questions:

• Should we provide built-in support for the full dispersion relation? If yes, I'm not sure about the best option. We could add a dummy equation EulerEquation1D as I do in the docs and dispatch on that. This, however, has the disadvantage that it might lead to the confusion that the Euler equations are supported by DispersiveShallowWater.jl.
• I'm not sure about the hyperbolic Serre-Green-Naghdi equations. There is a dispersion relation given in Favrie and Gavrilyuk, A rapid numerical method for solving Serre-Green-Naghdi equations describing long free surface gravity waves in eq. (19), but I don't know to understand this. Since $\tau = 1/h$, I initially thought, we obtain a classical dispersion relation by $\tau_0 = 1/h_0$, but this doesn't work out dimension-wise. That this doesn't work can also be seen by eq. (18), which differs from the one in the other two papers mentioned in the code. There is an additional third power in the denominator. I think this is related to the fact that the equations are written in Lagrangian form and not in the classical variables, but I'm not sure. Do you know how this dispersion relation would needs to be translated to obtain the classical version, @ranocha?
• I wanted to double-check that the linear dispersion relation for the BBM equation written as $\eta_t + \sqrt{gD}\eta_x + \frac{3}{2}\sqrt{\frac{g}{D}}\eta\eta_x - \frac{1}{6}D^2\eta_{xxt} = 0$ gives indeed the same as the one for the BBM-BBM equations, but I get another factor of $5/2$ in front (linearizing with $\eta = h_0 + h'$ and $D = h_0$ gives $h'_t + 5/2 \sqrt{g h_0} h'_x - 1/6 h_0^2h_{xxt} = 0$). I haven't thought about that when implementing the BBM equation. I would expect the same dispersion relation as for the BBM-BBM equations. Do you have an idea how to deal with that, @ranocha?

Closes #131.

[codecov]

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 98.07%. Comparing base (3c3917e) to head (80f4715).
Report is 2 commits behind head on main.

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #168      +/-   ##
==========================================
+ Coverage   98.02%   98.07%   +0.04%     
==========================================
  Files          19       20       +1     
  Lines        1877     1920      +43     
==========================================
+ Hits         1840     1883      +43     
  Misses         37       37              

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[github-actions]

Benchmark Results

	main	80f4715...	main/80f47157b1d0c3...
bbm_1d/bbm_1d_basic.jl	14.8 ± 0.53 μs	14.5 ± 0.54 μs	1.03
bbm_1d/bbm_1d_fourier.jl	0.223 ± 0.31 ms	0.216 ± 0.019 ms	1.03
bbm_bbm_1d/bbm_bbm_1d_basic_reflecting.jl	0.0806 ± 0.00039 ms	0.0808 ± 0.00042 ms	0.997
bbm_bbm_1d/bbm_bbm_1d_dg.jl	0.034 ± 0.0005 ms	0.0348 ± 0.00059 ms	0.977
bbm_bbm_1d/bbm_bbm_1d_relaxation.jl	27.3 ± 0.42 μs	27.6 ± 0.39 μs	0.989
bbm_bbm_1d/bbm_bbm_1d_upwind_relaxation.jl	0.0488 ± 0.00054 ms	0.0486 ± 0.00085 ms	1
hyperbolic_serre_green_naghdi_1d/hyperbolic_serre_green_naghdi_dingemans.jl	4.11 ± 0.016 μs	4.19 ± 0.014 μs	0.98
serre_green_naghdi_1d/serre_green_naghdi_well_balanced.jl	0.197 ± 0.0082 ms	0.197 ± 0.008 ms	1
svaerd_kalisch_1d/svaerd_kalisch_1d_dingemans_relaxation.jl	0.146 ± 0.0038 ms	0.145 ± 0.0046 ms	1.01
time_to_load	1.88 ± 0.0057 s	1.89 ± 0.015 s	0.995

Benchmark Plots

A plot of the benchmark results have been uploaded as an artifact to the workflow run for this PR.

[ranocha]

I'm not sure about the hyperbolic Serre-Green-Naghdi equations. There is a dispersion relation given in Favrie and Gavrilyuk, A rapid numerical method for solving Serre-Green-Naghdi equations describing long free surface gravity waves in eq. (19), but I don't know to understand this. Since τ = 1 / h, I initially thought, we obtain a classical dispersion relation by τ_0 = 1 / h_0, but this doesn't work out dimension-wise. That this doesn't work can also be seen by eq. (18), which differs from the one in the other two papers mentioned in the code. There is an additional third power in the denominator. I think this is related to the fact that the equations are written in Lagrangian form and not in the classical variables, but I'm not sure. Do you know how this dispersion relation would needs to be translated to obtain the classical version?

Good question... Maybe it would be best to compute the relation from scratch?

[ranocha]

I wanted to double-check that the linear dispersion relation for the BBM equation written as [...] I haven't thought about that when implementing the BBM equation. I would expect the same dispersion relation as for the BBM-BBM equations. Do you have an idea how to deal with that?

Why do you expect the dispersion relations to be the same?

[JoshuaLampert]

I wanted to double-check that the linear dispersion relation for the BBM equation written as [...] I haven't thought about that when implementing the BBM equation. I would expect the same dispersion relation as for the BBM-BBM equations. Do you have an idea how to deal with that?

Why do you expect the dispersion relations to be the same?

I think it would be pretty weird if the phase speed for the BBM equation would be always 5/2 as fast as the phase speed of the BBM-BBM equations. This would also mean that for $k \to 0$, $c$ would not converge to $c_0$, but $to $5/2 c_0$, which would also be strange, I think. But I might have an idea, where this discrepancy comes from. If we linearize the BBM equation around $\eta = 0 + \eta'$, i.e., we set $\eta_0 = 0$ (which is also what we are doing for the BBM-BBM equations), we don't get any linear contribution from the $3/2\sqrt{g/D}\eta\eta_x$ term anymore, i.e., we lose the additional 3/2 part, which in the end leads to the same linear dispersion relation. Thus, I think, the BBM equation also only makes sense for $\eta_0$ as for the BBM-BBM equations. To verify this also numerically, I ran two scripts similar to the example in the new section I added in this PR.

First version, where we use the true $\eta_0$ in the equations and I changed the linear dispersion relation to include the additional factor 5/2.

using DispersiveShallowWater
using Plots
using OrdinaryDiffEqTsit5
using Printf

h0 = eta0 = 0.8
g = 9.81
disp_rel = LinearDispersionRelation(h0)
euler = EulerEquations1D(; gravity_constant = g, eta0 = eta0)
equations = BBMEquation1D(; gravity_constant = g, eta0 = eta0, D = eta0)

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)
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, equations, 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 = wave_speed(disp_rel, equations, k)
anim = @animate for step in eachindex(sol.u)
    t = sol.t[step]
    x_t_euler = x + c_euler * t
    x_t = x + c * t
    index = argmin(abs.(DispersiveShallowWater.grid(semi) .- x_t))
    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], [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" "BBM"],
          color = [:blue :green])
end
gif(anim, "traveling_waves_bbm_nonshift.gif", fps = 25)

As a result, we obtain

We see that this exactly gives the numerical phase speed, which is (of course) much too fast.

Second version, where we use `eta0 = 0.0` in the equations and translate the initial condition and the final solution accordingly (we need the same trick for the BBM-BBM equations to reproduce the phase speed predicted by the linear dispersion relation). Here, I kept the dispersion relation for the BBM equation as the same as for the BBM-BBM equations, i.e., without the factor $5/2$.

using DispersiveShallowWater
using Plots
using OrdinaryDiffEqTsit5
using Printf

h0 = eta0 = 0.8
g = 9.81
disp_rel = LinearDispersionRelation(h0)
euler = EulerEquations1D(; gravity_constant = g, eta0 = eta0)
equations = BBMEquation1D(; gravity_constant = g, eta0 = 0.0, D = eta0)
# equations = BBMBBMEquations1D(; gravity_constant = g, eta0 = 0.0)

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 = 0.8
    A = 0.02
    h = A * cos(k * x - omega * t)
    v = sqrt(equations.gravity / k * tanh(k * h0)) * h / h0
    eta = h
    D = h0
    return SVector(eta)
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, equations, 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 = wave_speed(disp_rel, equations, k)
shifted_waterheight(q, equations) = waterheight_total(q, equations) + 0.8
DispersiveShallowWater.varnames(shifted_waterheight, equations) = ("η",)
anim = @animate for step in eachindex(sol.u)
    t = sol.t[step]
    x_t_euler = x + c_euler * t
    x_t = x + c * t
    index = argmin(abs.(DispersiveShallowWater.grid(semi) .- x_t))
    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], [eta[index] + 0.8],
             color = :green, label = nothing)
    plot!(semi => sol, plot_initial = true, plot_bathymetry = false,
          conversion = shifted_waterheight, step = step, legend = :topleft, linewidth = 2,
          plot_title = @sprintf("t = %.3f", t), yrange = (0.77, 0.83),
          linestyles = [:solid :dot], labels = ["Euler" "BBM"],
          color = [:blue :green])
end
gif(anim, "traveling_waves_bbm_shift.gif", fps = 25)

Here, the result is:

Again, we match the numerically observed phase speed, but this time the speed makes much more sense as it is only a bit too slow as expected.

To conclude, I would add a similar warning as for the BBM-BBM equations also for the BBM equation. Under this assumption, both equations have the same dispersion relation.
I am always so confused about these different notations and when we need which assumptions. It seems like no-one is making this explicit, but I think it makes sense now and is consistent.

[JoshuaLampert]

I'm not sure about the hyperbolic Serre-Green-Naghdi equations. There is a dispersion relation given in Favrie and Gavrilyuk, A rapid numerical method for solving Serre-Green-Naghdi equations describing long free surface gravity waves in eq. (19), but I don't know to understand this. Since τ = 1 / h, I initially thought, we obtain a classical dispersion relation by τ_0 = 1 / h_0, but this doesn't work out dimension-wise. That this doesn't work can also be seen by eq. (18), which differs from the one in the other two papers mentioned in the code. There is an additional third power in the denominator. I think this is related to the fact that the equations are written in Lagrangian form and not in the classical variables, but I'm not sure. Do you know how this dispersion relation would needs to be translated to obtain the classical version?

Good question... Maybe it would be best to compute the relation from scratch?

I tried to compute the linear dispersion relation for the hyperbolic Serre-Green-Naghdi equations from scratch, but had some issues to eliminate the $\eta$ component (what we call H here). I don't want to spend more time on this now. I would need to think about it again. I would be fine with just not defining a linear dispersion relation for the hyperbolic SGN equations for now and merge this PR soon. We can always add the linear dispersion relation later once we found a way to express it. What do you think, @ranocha?

[ranocha]

I wanted to double-check that the linear dispersion relation for the BBM equation written as [...] I haven't thought about that when implementing the BBM equation. I would expect the same dispersion relation as for the BBM-BBM equations. Do you have an idea how to deal with that?

Why do you expect the dispersion relations to be the same?

I think it would be pretty weird if the phase speed for the BBM equation would be always 5/2 as fast as the phase speed of the BBM-BBM equations. This would also mean that for k → 0 , c would not converge to c 0 , but $to 5 / 2 c 0 , which would also be strange, I think. But I might have an idea, where this discrepancy comes from. If we linearize the BBM equation around η = 0 + η ′ , i.e., we set η 0 = 0 (which is also what we are doing for the BBM-BBM equations), we don't get any linear contribution from the 3 / 2 g / D η η x term anymore, i.e., we lose the additional 3/2 part, which in the end leads to the same linear dispersion relation. Thus, I think, the BBM equation also only makes sense for η 0 as for the BBM-BBM equations. To verify this also numerically, I ran two scripts similar to the example in the new section I added in this PR.

First version, where we use the true
η
0
in the equations and I changed the linear dispersion relation to include the additional factor 5/2.
Second version, where we use eta0 = 0.0 in the equations and translate the initial condition and the final solution accordingly (we need the same trick for the BBM-BBM equations to reproduce the phase speed predicted by the linear dispersion relation). Here, I kept the dispersion relation for the BBM equation as the same as for the BBM-BBM equations, i.e., without the factor
5
/
2
.
To conclude, I would add a similar warning as for the BBM-BBM equations also for the BBM equation. Under this assumption, both equations have the same dispersion relation. I am always so confused about these different notations and when we need which assumptions. It seems like no-one is making this explicit, but I think it makes sense now and is consistent.

Nice 👍 I also ran into the issue of needing to set eta0 = 0 and translating the data accordingly. Can we automate this or generalize the system so that other values are also allowed? (not in this PR of course)

[ranocha]

I'm not sure about the hyperbolic Serre-Green-Naghdi equations. There is a dispersion relation given in Favrie and Gavrilyuk, A rapid numerical method for solving Serre-Green-Naghdi equations describing long free surface gravity waves in eq. (19), but I don't know to understand this. Since τ = 1 / h, I initially thought, we obtain a classical dispersion relation by τ_0 = 1 / h_0, but this doesn't work out dimension-wise. That this doesn't work can also be seen by eq. (18), which differs from the one in the other two papers mentioned in the code. There is an additional third power in the denominator. I think this is related to the fact that the equations are written in Lagrangian form and not in the classical variables, but I'm not sure. Do you know how this dispersion relation would needs to be translated to obtain the classical version?

Good question... Maybe it would be best to compute the relation from scratch?

I tried to compute the linear dispersion relation for the hyperbolic Serre-Green-Naghdi equations from scratch, but had some issues to eliminate the η component (what we call H here). I don't want to spend more time on this now. I would need to think about it again. I would be fine with just not defining a linear dispersion relation for the hyperbolic SGN equations for now and merge this PR soon. We can always add the linear dispersion relation later once we found a way to express it. What do you think, @ranocha?

Sure 👍

[JoshuaLampert]

I wanted to double-check that the linear dispersion relation for the BBM equation written as [...] I haven't thought about that when implementing the BBM equation. I would expect the same dispersion relation as for the BBM-BBM equations. Do you have an idea how to deal with that?

Why do you expect the dispersion relations to be the same?

I think it would be pretty weird if the phase speed for the BBM equation would be always 5/2 as fast as the phase speed of the BBM-BBM equations. This would also mean that for k → 0 , c would not converge to c 0 , but $to 5 / 2 c 0 , which would also be strange, I think. But I might have an idea, where this discrepancy comes from. If we linearize the BBM equation around η = 0 + η ′ , i.e., we set η 0 = 0 (which is also what we are doing for the BBM-BBM equations), we don't get any linear contribution from the 3 / 2 g / D η η x term anymore, i.e., we lose the additional 3/2 part, which in the end leads to the same linear dispersion relation. Thus, I think, the BBM equation also only makes sense for η 0 as for the BBM-BBM equations. To verify this also numerically, I ran two scripts similar to the example in the new section I added in this PR.
First version, where we use the true
η
0
in the equations and I changed the linear dispersion relation to include the additional factor 5/2.
Second version, where we use eta0 = 0.0 in the equations and translate the initial condition and the final solution accordingly (we need the same trick for the BBM-BBM equations to reproduce the phase speed predicted by the linear dispersion relation). Here, I kept the dispersion relation for the BBM equation as the same as for the BBM-BBM equations, i.e., without the factor
5
/
2
.
To conclude, I would add a similar warning as for the BBM-BBM equations also for the BBM equation. Under this assumption, both equations have the same dispersion relation. I am always so confused about these different notations and when we need which assumptions. It seems like no-one is making this explicit, but I think it makes sense now and is consistent.

Nice 👍 I also ran into the issue of needing to set eta0 = 0 and translating the data accordingly. Can we automate this or generalize the system so that other values are also allowed? (not in this PR of course)

I'm thinking about this again. What we could do for the linear dispersion relation is to set it for the BBM equation as $k\sqrt{g h_0}(1 + 3/2\eta_0/h_0)/(1 + 1/6(h_0k)^2)$. This way we always match the wave speed for any $\eta_0$ and we obtain the classical one by setting $\eta_0 = 0$. I wonder if in principal we want to allow $\eta_0 \neq 0$ or if we always want to explicitly set it 0 for the BBM and the BBM-BBM equations. I would need to think about it again. IIRC, we used $\eta_0 \neq 0$ for the well-balancedness test. I would need to check whether we really need $\eta_0 \neq 0$ there.

[JoshuaLampert]

I just pushed a commit to generalize the linear dispersions relations for both the BBM and the BBM-BBM equations for general $\eta_0$ and consolidated a more consistent naming of the different variables involved, see 824be7e. This way we can use any $\eta_0$ with any equations and we numerically reproduce the wave speed given by the linear dispersion relation without needing to shift anything. I still think that $\eta_0 = 0$ is the best choice in most applications for the BBM and BBM-BBM equations, but I also think that general $\eta_0$ can be supported.