Fix initial condition from Dingemans experiment for BBM-BBM equations (…

…#91)

* fix initial condition from Dingemans experiment for BBM-BBM equations

use notation related to still water depth consistently

* apply formatter v1.0.45

* remove gifs again (whoops)

* fix docstring

examples/bbm_bbm_variable_bathymetry_1d/bbm_bbm_variable_bathymetry_1d_dingemans.jl

@@ -4,7 +4,7 @@ using DispersiveShallowWater
 ###############################################################################
 # Semidiscretization of the BBM-BBM equations
 
-equations = BBMBBMVariableEquations1D(gravity_constant = 9.81, eta0 = 0.8)
+equations = BBMBBMVariableEquations1D(gravity_constant = 9.81)
 
 initial_condition = initial_condition_dingemans
 boundary_conditions = boundary_condition_periodic

examples/bbm_bbm_variable_bathymetry_1d/bbm_bbm_variable_bathymetry_1d_well_balanced.jl

@@ -6,7 +6,7 @@ using SparseArrays: sparse
 ###############################################################################
 # Semidiscretization of the BBM-BBM equations
 
-equations = BBMBBMVariableEquations1D(gravity_constant = 1.0, eta0 = 2.0)
+equations = BBMBBMVariableEquations1D(gravity_constant = 1.0)
 
 # Setup a truly discontinuous bottom topography function for this academic
 # testcase of well-balancedness. The errors from the analysis callback are
@@ -18,11 +18,11 @@ function initial_condition_discontinuous_well_balancedness(x, t,
     # Set the background values
     eta = equations.eta0
     v = 0.0
-    D = -1.0
+    D = equations.eta0 - 1.0
 
     # Setup a discontinuous bottom topography
     if x >= 0.5 && x <= 0.75
-        D = -1.5 - 0.5 * sinpi(2.0 * x)
+        D = equations.eta0 - 1.5 - 0.5 * sinpi(2.0 * x)
     end
 
     return SVector(eta, v, D)

examples/svaerd_kalisch_1d/svaerd_kalisch_1d_well_balanced.jl

@@ -19,11 +19,11 @@ function initial_condition_discontinuous_well_balancedness(x, t,
     # Set the background values
     eta = equations.eta0
     v = 0.0
-    D = -1.0
+    D = equations.eta0 - 1.0
 
     # Setup a discontinuous bottom topography
     if x >= 0.5 && x <= 0.75
-        D = -1.5 - 0.5 * sinpi(2.0 * x)
+        D = equations.eta0 - 1.5 - 0.5 * sinpi(2.0 * x)
     end
 
     return SVector(eta, v, D)

src/equations/bbm_bbm_1d.jl

@@ -1,5 +1,5 @@
 @doc raw"""
-    BBMBBMEquations1D(gravity, D)
+    BBMBBMEquations1D(gravity, D, eta0 = 0.0)
 
 BBM-BBM (Benjamin–Bona–Mahony) system in one spatial dimension. The equations are given by
 ```math
@@ -9,8 +9,8 @@ BBM-BBM (Benjamin–Bona–Mahony) system in one spatial dimension. The equation
 \end{aligned}
 ```
 The unknown quantities of the BBM-BBM equations are the total water height ``\eta`` and the velocity ``v``.
-The gravitational constant is denoted by `g` and the constant bottom topography (bathymetry) ``b = -D``. The water height above the bathymetry is therefore given by
-``h = \eta + D``.
+The gravitational constant is denoted by `g` and the constant bottom topography (bathymetry) ``b = \eta_0 - D``. The water height above the bathymetry is therefore given by
+``h = \eta - \eta_0 + D``. The BBM-BBM equations are only implemented for ``\eta_0 = 0``.
 
 One reference for the BBM-BBM system can be found in
 - Jerry L. Bona, Min Chen (1998)
@@ -21,10 +21,12 @@ One reference for the BBM-BBM system can be found in
 struct BBMBBMEquations1D{RealT <: Real} <: AbstractBBMBBMEquations{1, 2}
     gravity::RealT # gravitational constant
     D::RealT       # constant bathymetry
+    eta0::RealT    # constant still-water surface
 end
 
-function BBMBBMEquations1D(; gravity_constant, D = 1.0)
-    BBMBBMEquations1D(gravity_constant, D)
+function BBMBBMEquations1D(; gravity_constant, D = 1.0, eta0 = 0.0)
+    eta0 == 0.0 || @warn "The still-water surface needs to be 0 for the BBM-BBM equations"
+    BBMBBMEquations1D(gravity_constant, D, eta0)
 end
 
 varnames(::typeof(prim2prim), ::BBMBBMEquations1D) = ("η", "v")
@@ -94,11 +96,12 @@ function create_cache(mesh,
                       initial_condition,
                       RealT,
                       uEltype)
+    D = equations.D
     if solver.D1 isa PeriodicDerivativeOperator ||
        solver.D1 isa UniformPeriodicCoupledOperator
-        invImD2 = inv(I - 1 / 6 * equations.D^2 * Matrix(solver.D2))
+        invImD2 = inv(I - 1 / 6 * D^2 * Matrix(solver.D2))
     elseif solver.D1 isa PeriodicUpwindOperators
-        invImD2 = inv(I - 1 / 6 * equations.D^2 * Matrix(solver.D2))
+        invImD2 = inv(I - 1 / 6 * D^2 * Matrix(solver.D2))
     else
         @error "unknown type of first-derivative operator"
     end
@@ -122,15 +125,15 @@ function rhs!(du_ode, u_ode, t, mesh, equations::BBMBBMEquations1D, initial_cond
     deta = view(dq, 1, :)
     dv = view(dq, 2, :)
 
+    D = equations.D
     # energy and mass conservative semidiscretization
     if solver.D1 isa PeriodicDerivativeOperator ||
        solver.D1 isa UniformPeriodicCoupledOperator
-        @timeit timer() "deta hyperbolic" deta[:]=-solver.D1 * (equations.D * v + eta .* v)
+        @timeit timer() "deta hyperbolic" deta[:]=-solver.D1 * (D * v + eta .* v)
         @timeit timer() "dv hyperbolic" dv[:]=-solver.D1 *
                                               (equations.gravity * eta + 0.5 * v .^ 2)
     elseif solver.D1 isa PeriodicUpwindOperators
-        @timeit timer() "deta hyperbolic" deta[:]=-solver.D1.central *
-                                                  (equations.D * v + eta .* v)
+        @timeit timer() "deta hyperbolic" deta[:]=-solver.D1.central * (D * v + eta .* v)
         @timeit timer() "dv hyperbolic" dv[:]=-solver.D1.central *
                                               (equations.gravity * eta + 0.5 * v .^ 2)
     else
@@ -148,15 +151,15 @@ end
 @inline function prim2cons(q, equations::BBMBBMEquations1D)
     eta, v = q
 
-    h = eta + equations.D
+    h = eta - equations.eta0 + equations.D
     hv = h * v
     return SVector(h, hv)
 end
 
 @inline function cons2prim(u, equations::BBMBBMEquations1D)
     h, hv = u
 
-    eta = h - equations.D
+    eta = h + equations.eta0 - equations.D
     v = hv / h
     return SVector(eta, v)
 end
@@ -170,7 +173,7 @@ end
 end
 
 @inline function bathymetry(q, equations::BBMBBMEquations1D)
-    return -equations.D
+    return equations.eta0 - equations.D
 end
 
 @inline function waterheight(q, equations::BBMBBMEquations1D)
@@ -179,7 +182,8 @@ end
 
 @inline function energy_total(q, equations::BBMBBMEquations1D)
     eta, v = q
-    e = 0.5 * (equations.gravity * eta^2 + (equations.D + eta) * v^2)
+    D = still_waterdepth(q, equations)
+    e = 0.5 * (equations.gravity * eta^2 + (D + eta - equations.eta0) * v^2)
     return e
 end
 

src/equations/bbm_bbm_variable_bathymetry_1d.jl

@@ -9,8 +9,8 @@ BBM-BBM (Benjamin–Bona–Mahony) system in one spatial dimension with spatiall
 \end{aligned}
 ```
 The unknown quantities of the BBM-BBM equations are the total water height ``\eta`` and the velocity ``v``.
-The gravitational constant is denoted by `g` and the bottom topography (bathymetry) ``b = -D``. The water height above the bathymetry is therefore given by
-``h = \eta + D``.
+The gravitational constant is denoted by `g` and the bottom topography (bathymetry) ``b = \eta_0 - D``. The water height above the bathymetry is therefore given by
+``h = \eta - \eta_0 + D``. The BBM-BBM equations are only implemented for ``\eta_0 = 0``.
 
 One reference for the BBM-BBM system with spatially varying bathymetry can be found in
 - Samer Israwi, Henrik Kalisch, Theodoros Katsaounis, Dimitrios Mitsotakis (2022)
@@ -20,10 +20,11 @@ One reference for the BBM-BBM system with spatially varying bathymetry can be fo
 """
 struct BBMBBMVariableEquations1D{RealT <: Real} <: AbstractBBMBBMEquations{1, 3}
     gravity::RealT # gravitational constant
-    eta0::RealT    # constant "lake-at-rest" total water height
+    eta0::RealT    # constant still-water surface
 end
 
-function BBMBBMVariableEquations1D(; gravity_constant, eta0 = 1.0)
+function BBMBBMVariableEquations1D(; gravity_constant, eta0 = 0.0)
+    eta0 == 0.0 || @warn "The still-water surface needs to be 0 for the BBM-BBM equations"
     BBMBBMVariableEquations1D(gravity_constant, eta0)
 end
 
@@ -107,6 +108,11 @@ The initial condition that uses the dispersion relation of the Euler equations
 to approximate waves generated by a wave maker as it is done by experiments of
 Dingemans. The topography is a trapezoidal.
 
+!!! warning "Translation of water height"
+    The initial condition for the water height is translated to be around 0, which is
+    needed for the simulation because the `BBMBBMVariableEquations1D` are only implemented
+    for ``\eta_0 = 0``.
+
 References:
 - Magnus Svärd, Henrik Kalisch (2023)
   A novel energy-bounded Boussinesq model and a well-balanced and stable numerical discretization
@@ -116,16 +122,16 @@ References:
   [link](https://repository.tudelft.nl/islandora/object/uuid:c2091d53-f455-48af-a84b-ac86680455e9/datastream/OBJ/download)
 """
 function initial_condition_dingemans(x, t, equations::BBMBBMVariableEquations1D, mesh)
-    eta0 = 0.8
+    h0 = 0.8
     A = 0.02
     # omega = 2*pi/(2.02*sqrt(2))
-    k = 0.8406220896381442 # precomputed result of find_zero(k -> omega^2 - equations.gravity * k * tanh(k * eta0), 1.0) using Roots.jl
+    k = 0.8406220896381442 # precomputed result of find_zero(k -> omega^2 - equations.gravity * k * tanh(k * h0), 1.0) using Roots.jl
     if x < -30.5 * pi / k || x > -8.5 * pi / k
         h = 0.0
     else
         h = A * cos(k * x)
     end
-    v = sqrt(equations.gravity / k * tanh(k * eta0)) * h / eta0
+    v = sqrt(equations.gravity / k * tanh(k * h0)) * h / h0
     if 11.01 <= x && x < 23.04
         b = 0.6 * (x - 11.01) / (23.04 - 11.01)
     elseif 23.04 <= x && x < 27.04
@@ -135,8 +141,10 @@ function initial_condition_dingemans(x, t, equations::BBMBBMVariableEquations1D,
     else
         b = 0.0
     end
-    eta = h + eta0
-    D = -b
+    # Here, we compute eta - h0!! To obtain the original eta, h0 = 0.8 needs to be added again!
+    # This is because the BBM-BBM equations are only implemented for eta0 = 0
+    eta = h
+    D = h0 - b
     return SVector(eta, v, D)
 end
 
@@ -150,7 +158,7 @@ function create_cache(mesh,
     D = Array{RealT}(undef, nnodes(mesh))
     x = grid(solver)
     for i in eachnode(solver)
-        D[i] = initial_condition(x[i], 0.0, equations, mesh)[3]
+        D[i] = still_waterdepth(initial_condition(x[i], 0.0, equations, mesh), equations)
     end
     K = Diagonal(D .^ 2)
     if solver.D1 isa PeriodicDerivativeOperator ||
@@ -164,7 +172,7 @@ function create_cache(mesh,
         @error "unknown type of first-derivative operator"
     end
     tmp1 = Array{RealT}(undef, nnodes(mesh)) # tmp1 is needed for the `RelaxationCallback`
-    return (invImDKD = invImDKD, invImD2K = invImD2K, tmp1 = tmp1)
+    return (invImDKD = invImDKD, invImD2K = invImD2K, D = D, tmp1 = tmp1)
 end
 
 # Discretization that conserves the mass (for eta and v) and the energy for periodic boundary conditions, see
@@ -175,14 +183,13 @@ end
 function rhs!(du_ode, u_ode, t, mesh, equations::BBMBBMVariableEquations1D,
               initial_condition, ::BoundaryConditionPeriodic, source_terms,
               solver, cache)
-    @unpack invImDKD, invImD2K = cache
+    @unpack invImDKD, invImD2K, D = cache
 
     q = wrap_array(u_ode, mesh, equations, solver)
     dq = wrap_array(du_ode, mesh, equations, solver)
 
     eta = view(q, 1, :)
     v = view(q, 2, :)
-    D = view(q, 3, :)
     deta = view(dq, 1, :)
     dv = view(dq, 2, :)
     dD = view(dq, 3, :)
@@ -212,9 +219,9 @@ end
 @inline function prim2cons(q, equations::BBMBBMVariableEquations1D)
     eta, v, D = q
 
-    h = eta + D
+    b = bathymetry(q, equations)
+    h = eta - b
     hv = h * v
-    b = -D
     return SVector(h, hv, b)
 end
 
@@ -223,7 +230,7 @@ end
 
     eta = h + b
     v = hv / h
-    D = -b
+    D = equations.eta0 - b
     return SVector(eta, v, D)
 end
 
@@ -236,7 +243,8 @@ end
 end
 
 @inline function bathymetry(q, equations::BBMBBMVariableEquations1D)
-    return -q[3]
+    D = q[3]
+    return equations.eta0 - D
 end
 
 @inline function waterheight(q, equations::BBMBBMVariableEquations1D)
@@ -245,7 +253,7 @@ end
 
 @inline function energy_total(q, equations::BBMBBMVariableEquations1D)
     eta, v, D = q
-    e = 0.5 * (equations.gravity * eta^2 + (D + eta) * v^2)
+    e = 0.5 * (equations.gravity * eta^2 + (D + eta - equations.eta0) * v^2)
     return e
 end
 

src/equations/equations.jl

@@ -133,6 +133,12 @@ See [`momentum`](@ref).
 
 varnames(::typeof(discharge), equations) = ("P",)
 
+@inline function still_waterdepth(q, equations::AbstractEquations)
+    b = bathymetry(q, equations)
+    D = equations.eta0 - b
+    return D
+end
+
 """
     entropy(q, equations)