Merge branch 'main' into patch-1

Author: JoshuaLampert

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 when updating to v0.8 from v0.7.x
+
+#### Changed
+
+- The initial condition `initial_condition_dingemans` uses equation-specific offsets to account for phase shifts
+  compared to the experimental data and a slightly different range for the initial wave train ([#237]).
 
 ## Changes in the v0.7 lifecycle
 

Project.toml

@@ -1,7 +1,7 @@
 name = "DispersiveShallowWater"
 uuid = "4ab875f6-b79a-4e28-9745-4f0293954817"
 authors = ["Joshua Lampert <joshua.lampert@uni-hamburg.de>"]
-version = "0.7.8-DEV"
+version = "0.8.1-DEV"
 
 [deps]
 BandedMatrices = "aae01518-5342-5314-be14-df237901396f"

benchmark/benchmarks.jl

@@ -22,7 +22,7 @@ for elixir in elixirs
     benchname = joinpath(basename(dirname(elixir)), basename(elixir)) * " - rhs!:"
     println("Running $benchname...")
     redirect_stdout(devnull) do
-        trixi_include(@__MODULE__, elixir, tspan = (0.0, 1e-10))
+        trixi_include(@__MODULE__, elixir, tspan = (0.0, 1e-11))
     end
     SUITE[benchname] = @benchmarkable DispersiveShallowWater.rhs!($(similar(sol.u[end])),
                                                                   $(copy(sol.u[end])),

docs/src/basic_example.md

@@ -72,11 +72,11 @@ Lastly, we define the physical domain as the interval from -130 to 20 and we cho
 
 ## Define numerical solver
 
-In the next step, we build a [`Semidiscretization`](@ref) that bundles all ingredients for the spatial discretization of the model. Especially, we need to define a [`Solver`](@ref). 
+In the next step, we build a [`Semidiscretization`](@ref) that bundles all ingredients for the spatial discretization of the model. Especially, we need to define a [`Solver`](@ref).
 
 The simplest way to define a solver when working with [`boundary_condition_periodic`](@ref) is to call the constructor by providing the mesh and a desired order of accuracy.
 
-In the following example, we use an accuracy order of 4. The default constructor simply creates periodic first-, second-, and third-derivative central finite difference summation-by-parts (SBP) operators of the provided order of accuracy. 
+In the following example, we use an accuracy order of 4. The default constructor simply creates periodic first-, second-, and third-derivative central finite difference summation-by-parts (SBP) operators of the provided order of accuracy.
 
 How to use other summation-by-parts operators, is described in the section on [how to customize the solver](@ref customize_solver). Note that for non-periodic boundary conditions, the solver also needs to be created with non-periodic
 operators, see, e.g. [examples/bbm\_bbm\_1d/bbm\_bbm\_1d\_basic\_reflecting.jl](https://github.com/NumericalMathematics/DispersiveShallowWater.jl/blob/main/examples/bbm_bbm_1d/bbm_bbm_1d_basic_reflecting.jl).
@@ -127,11 +127,10 @@ nothing # hide
 
 ![shoaling solution](shoaling_solution.png)
 
-By default, this will plot the bathymetry, but not the initial (analytical) solution. 
+By default, this will plot the bathymetry, but not the initial (analytical) solution.
 
 You can adjust this by passing the boolean values `plot_bathymetry` (if `true`, always plot bathymetry in the first subplot) and `plot_initial`. Note that `plot_initial = true` will evaluate and plot the initial condition function at the same time `t` as the numerical solution being displayed (the final time by default). This means if your initial condition function represents an analytical solution, setting `plot_initial = true` will plot the analytical solution at that specific time for comparison.
 
-
 Plotting an animation over time can, e.g., be done by the following command, which uses `step` to plot the solution at a specific time step. Here `conversion = waterheight_total` makes it so that we only look at the waterheight ``\eta`` and not also the velocity ``v``. More on tutorials for plotting can be found in the chapter [Plotting Simulation Results](@ref plotting).
 
 ```@example overview
@@ -150,7 +149,6 @@ It is also possible to plot the solution variables at a fixed spatial point over
 
 More examples sorted by the simulated equations can be found in the [examples/](https://github.com/NumericalMathematics/DispersiveShallowWater.jl/tree/main/examples) subdirectory.
 
-
 ## [Plain program](@id overview-plain-program)
 
 Here follows a version of the program without any comments.

docs/src/callbacks.md

@@ -1,14 +1,15 @@
 # Callbacks
 
 Callbacks provide additional functionality during simulations, such as monitoring solution properties, analyzing errors, or ensuring conservation of physical quantities.
-[DispersiveShallowWater.jl](https://github.com/NumericalMathematics/DispersiveShallowWater.jl implements three main callback types that can be used individually or in combination to enhance simulation analysis and performance monitoring.
+[DispersiveShallowWater.jl](https://github.com/NumericalMathematics/DispersiveShallowWater.jl) implements three main callback types that can be used individually or in combination to enhance simulation analysis and performance monitoring.
 
 When using multiple callbacks simultaneously, combine them using a `CallbackSet`:
 
 ```julia
 callbacks = CallbackSet(analysis_callback, summary_callback)
 sol = solve(ode, Tsit5(), callback = callbacks)
 ```
+
 More information on the usage of callbacks within the SciML framework can be found [in the documentation](https://docs.sciml.ai/DiffEqDocs/stable/features/callback_functions/).
 
 ## Summary Callback
@@ -28,7 +29,7 @@ sol = solve(ode, Tsit5(), callback = summary_callback)
 
 At the end of the simulation, the callback will display output similar to:
 
-```
+```julia
 ───────────────────────────────────────────────────────────────────────────────────────────
               DispersiveSWE                       Time                    Allocations
                                          ───────────────────────   ────────────────────────
@@ -46,7 +47,6 @@ analyze solution                      3   13.2ms    3.0%  4.39ms    147KiB    0.
 ───────────────────────────────────────────────────────────────────────────────────────────
 ```
 
-
 ## Analysis Callback
 
 The [`AnalysisCallback`](@ref) monitors solution quality and physical properties during the simulation. It computes error norms and tracks conservation of important physical quantities at specified time intervals.
@@ -189,4 +189,4 @@ For additional information on relaxation, how it works, and why and when it is u
 [^RanochaSayyariDalcinParsaniKetcheson2020]:
     Ranocha, Sayyari, Dalcin, Parsani, Ketcheson (2020):
     Relaxation Runge–Kutta Methods: Fully-Discrete Explicit Entropy-Stable Schemes for the Compressible Euler and Navier–Stokes Equations.
-    [DOI: 10.1137/19M1263480](https://doi.org/10.1137/19M1263480)
+    [DOI: 10.1137/19M1263480](https://doi.org/10.1137/19M1263480)

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
 
@@ -33,17 +33,17 @@ Next, we set up the different equation systems we want to compare:
 bbmbbm = BBMBBMEquations1D(bathymetry_type = bathymetry_variable,
                            gravity = 9.81, eta0 = 0.0)
 
-# Svärd-Kalisch equations with specific parameter set
-sk = SvaerdKalischEquations1D(gravity = 9.81, eta0 = 0.8, alpha = 0.0,
-                              beta = 0.27946992481203003, gamma = 0.0521077694235589)
+# Svärd-Kalisch equations with specific parameter set (optimized for small wave numbers)
+sk = SvaerdKalischEquations1D(gravity = 9.81, eta0 = 0.8, alpha = 0.0004040404040404049,
+                              beta = 0.49292929292929294, gamma = 0.15707070707070708)
 
-# Serre-Green-Naghdi equations with variable bathymetry  
+# Serre-Green-Naghdi equations with variable bathymetry
 sgn = SerreGreenNaghdiEquations1D(bathymetry_type = bathymetry_variable,
-                                  gravity = 9.81)
+                                  gravity = 9.81, eta0 = 0.8)
 
 # Hyperbolic approximation of Serre-Green-Naghdi equations
 hysgn = HyperbolicSerreGreenNaghdiEquations1D(bathymetry_type = bathymetry_mild_slope,
-                                              lambda = 100.0, gravity = 9.81)
+                                              lambda = 100.0, gravity = 9.81, eta0 = 0.8)
                                               # 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.
@@ -134,8 +134,7 @@ y_limits = (-0.03, 0.87)
 # Model configurations: (semidiscretization, solution, label, conversion_function, linestyle)
 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 +163,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)
@@ -190,7 +189,6 @@ nothing # hide
 
 The results show how different dispersive wave models capture the wave evolution over the trapezoidal bathymetry.
 
-
 ## [Plain program](@id overview-plain-program-dingemans)
 
 Here follows a version of the program without any comments.
@@ -202,17 +200,17 @@ using DispersiveShallowWater, OrdinaryDiffEqTsit5, Plots
 bbmbbm = BBMBBMEquations1D(bathymetry_type = bathymetry_variable,
                            gravity = 9.81, eta0 = 0.0)
 
-# Svärd-Kalisch equations with specific parameter set
-sk = SvaerdKalischEquations1D(gravity = 9.81, eta0 = 0.8, alpha = 0.0,
-                              beta = 0.27946992481203003, gamma = 0.0521077694235589)
+# Svärd-Kalisch equations with specific parameter set (optimized for small wave numbers)
+sk = SvaerdKalischEquations1D(gravity = 9.81, eta0 = 0.8, alpha = 0.0004040404040404049,
+                              beta = 0.49292929292929294, gamma = 0.15707070707070708)
 
-# Serre-Green-Naghdi equations with variable bathymetry  
+# Serre-Green-Naghdi equations with variable bathymetry
 sgn = SerreGreenNaghdiEquations1D(bathymetry_type = bathymetry_variable,
-                                  gravity = 9.81)
+                                  gravity = 9.81, eta0 = 0.8)
 
 # Hyperbolic approximation of Serre-Green-Naghdi equations
 hysgn = HyperbolicSerreGreenNaghdiEquations1D(bathymetry_type = bathymetry_mild_slope,
-                                              lambda = 100.0, gravity = 9.81)
+                                              lambda = 100.0, gravity = 9.81, eta0 = 0.8)
                                               # 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.
@@ -262,8 +260,7 @@ y_limits = (-0.03, 0.87)
 # Model configurations: (semidiscretization, solution, label, conversion_function, linestyle)
 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 +289,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)
@@ -321,4 +318,4 @@ plot(all_plots...,
 [^Dingemans1997]:
     Dingemans (1997):
     Water Wave Propagation Over Uneven Bottoms (In 2 Parts).
-    [DOI: 10.1142/1241](https://doi.org/10.1142/1241)
+    [DOI: 10.1142/1241](https://doi.org/10.1142/1241)

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

docs/src/overview.md

@@ -22,12 +22,11 @@ The following table provides an overview of all available equation systems and t
 ![water height and bathymetry](bathymetry.png)
 
 - ``\eta``: Total water height
-- ``v``: Velocity in horizontal direction  
+- ``v``: Velocity in horizontal direction
 - ``D``: Still-water depth
 - ``w``: Auxiliary variable in hyperbolic approximation (``\approx -h v_x``)
 - ``H``: Auxiliary variable in hyperbolic approximation (``\approx h``)
 
-
 ## 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.

docs/src/plotting.md

@@ -41,7 +41,7 @@ nothing # hide
 
 # [Plotting Simulation Results](@id plotting)
 
-[DispersiveShallowWater.jl](https://github.com/NumericalMathematics/DispersiveShallowWater.jl) provides flexible plotting capabilities through [Plots.jl](https://github.com/JuliaPlots/Plots.jl) recipes. The plotting system supports various conversion functions, visualization options, and analysis tools. 
+[DispersiveShallowWater.jl](https://github.com/NumericalMathematics/DispersiveShallowWater.jl) provides flexible plotting capabilities through [Plots.jl](https://github.com/JuliaPlots/Plots.jl) recipes. The plotting system supports various conversion functions, visualization options, and analysis tools.
 
 [Makie.jl](https://docs.makie.org/stable/) is not supported yet. [Contributions are welcome](https://github.com/NumericalMathematics/DispersiveShallowWater.jl/issues/220).
 
@@ -59,7 +59,7 @@ using Plots
 
 t = 13.37 # plot solution at (roughly) t = 13.37s
 step_idx = argmin(abs.(saveat .- t)) # get the closest point to 13.37
-p1 = plot(semi => sol, conversion = prim2prim, plot_bathymetry = false, 
+p1 = plot(semi => sol, conversion = prim2prim, plot_bathymetry = false,
           suptitle = "Primitive Variables", step = step_idx)
 p2 = plot(semi => sol, conversion = prim2cons, plot_bathymetry = false,
           suptitle = "Conservative Variables", step = step_idx)

docs/src/solvers.md

@@ -42,7 +42,8 @@ solver = Solver(D1, D2, D3)
 ```
 
 where:
-- `D1` is always required and must be an `AbstractDerivativeOperator` 
+
+- `D1` is always required and must be an `AbstractDerivativeOperator`
 - `D2` and `D3` are optional and can be either `AbstractDerivativeOperator`s, `AbstractMatrix`es, or `nothing`
 
 ## Reflecting Boundary Conditions
@@ -68,6 +69,7 @@ Other possible choices for sources can be found in the [documentation of Summati
 For equations that benefit from upwind discretizations (such as the Serre-Green-Naghdi equations), you can use [upwind SBP Operators](@ref upwind_sbp):
 
 **For periodic boundary conditions:**
+
 ```julia
 using SummationByPartsOperators: upwind_operators, periodic_derivative_operator
 
@@ -81,6 +83,7 @@ solver = Solver(D1)
 ```
 
 **For reflecting boundary conditions:**
+
 ```julia
 using SummationByPartsOperators: Mattsson2017, upwind_operators
 
@@ -197,6 +200,7 @@ semi = Semidiscretization(mesh, equations, initial_condition, solver,
 ```
 
 The solver you choose should be compatible with your boundary conditions:
+
 - Periodic operators (created with `periodic_derivative_operator` or `fourier_derivative_operator`) require `boundary_condition_periodic`
 - Non-periodic operators (created with `MattssonNordström2004`, etc.) are needed for `boundary_condition_reflecting`
 
@@ -205,6 +209,7 @@ The solver you choose should be compatible with your boundary conditions:
 ### Mismatched Operators and Boundary Conditions
 
 **Problem**: Using periodic operators with reflecting boundary conditions or vice versa.
+
 ```julia
 # ❌ This will fail
 D1 = periodic_derivative_operator(1, 4, mesh.xmin, mesh.xmax, mesh.N)
@@ -214,9 +219,10 @@ semi = Semidiscretization(mesh, equations, initial_condition, solver,
 ```
 
 **Solution**: Match operator type to boundary conditions:
+
 ```julia
 # ✅ Correct approach
-D1 = derivative_operator(MattssonNordström2004(), derivative_order = 1, 
+D1 = derivative_operator(MattssonNordström2004(), derivative_order = 1,
                          accuracy_order = 4, xmin = mesh.xmin, xmax = mesh.xmax, N = mesh.N)
 solver = Solver(D1)
 semi = Semidiscretization(mesh, equations, initial_condition, solver,
@@ -226,6 +232,7 @@ semi = Semidiscretization(mesh, equations, initial_condition, solver,
 ### Incorrect Grid Size for DG Methods
 
 **Problem**: Grid size not divisible by polynomial degree + 1 for DG methods.
+
 ```julia
 # ❌ N = 101, p = 3, but 101 is not divisible by (3+1) = 4
 N = 101
@@ -234,6 +241,7 @@ mesh = Mesh1D(coordinates_min, coordinates_max, N)  # Error in DG setup!
 ```
 
 **Solution**: Ensure `N` is divisible by `p + 1`:
+
 ```julia
 # ✅ Adjust N to be divisible by p + 1
 p = 3