Changes to Dingemans initital condition

docs/src/dingemans.md

@@ -10,7 +10,7 @@ The bathymetry profile consists of:
 
 - A flat section from the wave maker to x = 11.01
 - A linearly increasing slope from x = 11.01 to x = 23.04 (maximum height of 0.6)
-- A flat plateau from x = 23.04 to x = 27.04  
+- A flat plateau from x = 23.04 to x = 27.04
 - A linearly decreasing slope from x = 27.04 to x = 33.07
 - A flat section beyond x = 33.07
 
@@ -37,13 +37,13 @@ bbmbbm = BBMBBMEquations1D(bathymetry_type = bathymetry_variable,
 sk = SvaerdKalischEquations1D(gravity = 9.81, eta0 = 0.8, alpha = 0.0,
                               beta = 0.27946992481203003, gamma = 0.0521077694235589)
 
-# Serre-Green-Naghdi equations with variable bathymetry  
+# Serre-Green-Naghdi equations with variable bathymetry
 sgn = SerreGreenNaghdiEquations1D(bathymetry_type = bathymetry_variable,
-                                  gravity = 9.81)
+                                  eta0 = 0.8, gravity = 9.81)
 
 # Hyperbolic approximation of Serre-Green-Naghdi equations
 hysgn = HyperbolicSerreGreenNaghdiEquations1D(bathymetry_type = bathymetry_mild_slope,
-                                              lambda = 100.0, gravity = 9.81)
+                                              eta0 = 0.8, lambda = 100.0, gravity = 9.81)
                                               # 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.
@@ -135,7 +135,7 @@ y_limits = (-0.03, 0.87)
 models = [
     (semi_bbmbbm, sol_bbmbbm, "BBM-BBM", shifted_waterheight, :solid),
     # Svärd-Kalisch has a phase shift which need to be adjusted for in the initial condition
-    # (semi_sk, sol_sk, "Svärd-Kalisch", waterheight_total, :dashdotdot), 
+    # (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)
 ]
@@ -164,7 +164,7 @@ for time_val in times
     push!(snapshot_plots, p)
 end
 
-# Create legend plot 
+# Create legend plot
 legend_plot = plot(legend=:top, framestyle=:none, legendfontsize=11)
 
 for (i, (_, _, label, _, linestyle)) in enumerate(models)
@@ -206,13 +206,13 @@ bbmbbm = BBMBBMEquations1D(bathymetry_type = bathymetry_variable,
 sk = SvaerdKalischEquations1D(gravity = 9.81, eta0 = 0.8, alpha = 0.0,
                               beta = 0.27946992481203003, gamma = 0.0521077694235589)
 
-# Serre-Green-Naghdi equations with variable bathymetry  
+# Serre-Green-Naghdi equations with variable bathymetry
 sgn = SerreGreenNaghdiEquations1D(bathymetry_type = bathymetry_variable,
-                                  gravity = 9.81)
+                                  eta0 = 0.8, gravity = 9.81)
 
 # Hyperbolic approximation of Serre-Green-Naghdi equations
 hysgn = HyperbolicSerreGreenNaghdiEquations1D(bathymetry_type = bathymetry_mild_slope,
-                                              lambda = 100.0, gravity = 9.81)
+                                              eta0 = 0.8, lambda = 100.0, gravity = 9.81)
                                               # 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.
@@ -263,7 +263,7 @@ y_limits = (-0.03, 0.87)
 models = [
     (semi_bbmbbm, sol_bbmbbm, "BBM-BBM", shifted_waterheight, :solid),
     # Svärd-Kalisch has a phase shift which need to be adjusted for in the initial condition
-    # (semi_sk, sol_sk, "Svärd-Kalisch", waterheight_total, :dashdotdot), 
+    # (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)
 ]
@@ -292,7 +292,7 @@ for time_val in times
     push!(snapshot_plots, p)
 end
 
-# Create legend plot 
+# Create legend plot
 legend_plot = plot(legend=:top, framestyle=:none, legendfontsize=11)
 
 for (i, (_, _, label, _, linestyle)) in enumerate(models)

docs/src/dispersion.md

@@ -88,9 +88,9 @@ function initial_condition_traveling_wave(x, t, equations, mesh)
     omega = frequency(k)
     h0 = reference_height()
     A = 0.02
-    h = A * cos(k * x - omega * t)
-    v = sqrt(equations.gravity / k * tanh(k * h0)) * h / h0
-    eta = h + equations.eta0
+    eta_prime = A * cos(k * x - omega * t)
+    v = sqrt(equations.gravity / k * tanh(k * h0)) * eta_prime / h0
+    eta = eta_prime + equations.eta0
     D = h0
     return SVector(eta, v, D)
 end

examples/hyperbolic_serre_green_naghdi_1d/hyperbolic_serre_green_naghdi_dingemans.jl

@@ -5,8 +5,7 @@ using DispersiveShallowWater
 # Semidiscretization of the hyperbolic Serre-Green-Naghdi equations
 
 equations = HyperbolicSerreGreenNaghdiEquations1D(bathymetry_type = bathymetry_mild_slope,
-                                                  lambda = 500.0,
-                                                  gravity = 9.81)
+                                                  eta0 = 0.8, lambda = 500.0, gravity = 9.81)
 
 initial_condition = initial_condition_dingemans
 boundary_conditions = boundary_condition_periodic

examples/serre_green_naghdi_1d/serre_green_naghdi_dingemans.jl

@@ -7,7 +7,7 @@ using SummationByPartsOperators: upwind_operators, periodic_derivative_operator
 
 # or bathymetry_mild_slope instead of bathymetry_variable
 equations = SerreGreenNaghdiEquations1D(bathymetry_type = bathymetry_variable,
-                                        gravity = 9.81)
+                                        eta0 = 0.8, gravity = 9.81)
 
 initial_condition = initial_condition_dingemans
 boundary_conditions = boundary_condition_periodic

src/equations/bbm_bbm_1d.jl

@@ -223,6 +223,8 @@ function source_terms_manufactured_reflecting(q, x, t,
     return SVector(s1, s2, zero(s1))
 end
 
+dingemans_calibration(equations::BBMBBMEquations1D) = 2.7
+
 """
     initial_condition_dingemans(x, t, equations::BBMBBMEquations1D, mesh)
 
@@ -249,12 +251,13 @@ function initial_condition_dingemans(x, t, equations::BBMBBMEquations1D, mesh)
     A = 0.02
     # omega = 2*pi/(2.02*sqrt(2))
     k = 0.8406220896381442 # precomputed result of find_zero(k -> omega^2 - g * k * tanh(k * h0), 1.0) using Roots.jl
-    if x < -30.5 * pi / k || x > -8.5 * pi / k
-        h = 0.0
+    offset = dingemans_calibration(equations)
+    if x - offset < -34.5 * pi / k || x - offset > -4.5 * pi / k
+        eta_prime = 0.0
     else
-        h = A * cos(k * x)
+        eta_prime = A * cos(k * (x - offset))
     end
-    v = sqrt(g / k * tanh(k * h0)) * h / h0
+    v = sqrt(g / k * tanh(k * h0)) * eta_prime / 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
@@ -264,9 +267,9 @@ function initial_condition_dingemans(x, t, equations::BBMBBMEquations1D, mesh)
     else
         b = 0.0
     end
-    # Here, we compute eta - h0!! To obtain the original eta, h0 = 0.8 needs to be added again!
+    # Here, we compute eta - eta0!! To obtain the original eta, eta0 = 0.8 needs to be added again!
     # This is because the BBM-BBM equations are only implemented for eta0 = 0
-    eta = h
+    eta = eta_prime
     D = h0 - b
     return SVector(eta, v, D)
 end

src/equations/equations.jl

@@ -614,6 +614,10 @@ function solve_system_matrix!(dv, system_matrix, ::Union{BBMEquation1D, BBMBBMEq
     ldiv!(system_matrix, dv)
 end
 
+# To match the experimental data from Dingemans, the initial condition needs to be shifted
+# to account for the phase shift to the experimental data.
+dingemans_calibration(equations) = 0
+
 """
     initial_condition_dingemans(x, t, equations::AbstractShallowWaterEquations, mesh)
 
@@ -635,12 +639,13 @@ function initial_condition_dingemans(x, t, equations::AbstractShallowWaterEquati
     A = 0.02
     # omega = 2*pi/(2.02*sqrt(2))
     k = 0.8406220896381442 # precomputed result of find_zero(k -> omega^2 - g * k * tanh(k * h0), 1.0) using Roots.jl
-    if x < -30.5 * pi / k || x > -8.5 * pi / k
-        h = 0.0
+    offset = dingemans_calibration(equations)
+    if x - offset < -34.5 * pi / k || x - offset > -4.5 * pi / k
+        eta_prime = 0.0
     else
-        h = A * cos(k * x)
+        eta_prime = A * cos(k * (x - offset))
     end
-    v = sqrt(g / k * tanh(k * h0)) * h / h0
+    v = sqrt(g / k * tanh(k * h0)) * eta_prime / 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
@@ -651,7 +656,7 @@ function initial_condition_dingemans(x, t, equations::AbstractShallowWaterEquati
         b = 0.0
     end
     eta0 = equations.eta0
-    eta = h + h0
+    eta = eta_prime + eta0
     D = eta0 - b
     return SVector(eta, v, D)
 end

src/equations/hyperbolic_serre_green_naghdi_1d.jl

@@ -321,6 +321,8 @@ function source_terms_manufactured_reflecting(q, x, t,
     return SVector(s1, s2, s3, s4, s5)
 end
 
+dingemans_calibration(equations::HyperbolicSerreGreenNaghdiEquations1D) = 2.4
+
 function create_cache(mesh, equations::HyperbolicSerreGreenNaghdiEquations1D,
                       solver, initial_condition,
                       boundary_conditions::Union{BoundaryConditionPeriodic,

src/equations/serre_green_naghdi_1d.jl

@@ -414,7 +414,7 @@ function source_terms_manufactured_reflecting(q, x, t,
     cos3_pix = cospi(x)^3
     cos4_pix = cospi(x)^4
 
-    # Compute s1 
+    # Compute s1
     s1 = 2(2 + cospix + x) + (-1 - 2(2 + cospix + x) * t) * sinpix * t -
          2(1 - pi * sinpix) * sinpix * t2 * x +
          cospix * (-1 - 2(2 + cospix + x) * t) * pi * t * x
@@ -496,6 +496,8 @@ function source_terms_manufactured_reflecting(q, x, t,
     return SVector(s1, s2, zero(s1))
 end
 
+dingemans_calibration(equations::SerreGreenNaghdiEquations1D) = 2.4
+
 # flat bathymetry with periodic or reflecting boundary conditions
 function create_cache(mesh,
                       equations::SerreGreenNaghdiEquations1D{BathymetryFlat},

src/equations/svaerd_kalisch_1d.jl

@@ -197,6 +197,8 @@ function source_terms_manufactured_reflecting(q, x, t, equations::SvaerdKalischE
     return SVector(s1, s2, zero(s1))
 end
 
+dingemans_calibration(equations::SvaerdKalischEquations1D) = 2.2
+
 # For periodic boundary conditions
 function assemble_system_matrix!(cache, h,
                                  ::SvaerdKalischEquations1D,