Add comparison with experimental data to docs

Project.toml

@@ -5,6 +5,7 @@ version = "0.8.1-DEV"
 
 [deps]
 BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
+DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
 FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
@@ -27,6 +28,7 @@ TrixiBase = "9a0f1c46-06d5-4909-a5a3-ce25d3fa3284"
 
 [compat]
 BandedMatrices = "1.5"
+DelimitedFiles = "1"
 DiffEqBase = "6.160"
 FastBroadcast = "0.3.5"
 ForwardDiff = "0.10.36, 1"

docs/src/dingemans.md

@@ -43,7 +43,7 @@ sgn = SerreGreenNaghdiEquations1D(bathymetry_type = bathymetry_variable,
 
 # Hyperbolic approximation of Serre-Green-Naghdi equations
 hysgn = HyperbolicSerreGreenNaghdiEquations1D(bathymetry_type = bathymetry_mild_slope,
-                                              lambda = 100.0, gravity = 9.81, eta0 = 0.8)
+                                              lambda = 200.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.
@@ -63,15 +63,15 @@ We create a computational domain that is large enough to contain the entire expe
 ```@example dingemans
 coordinates_min = -138.0
 coordinates_max = 46.0
-N = 512
+N = 1024
 mesh = Mesh1D(coordinates_min, coordinates_max, N)
 nothing # hide
 ```
 
-For the spatial discretization, we use fourth-order accurate [summation-by-parts operators](@ref sbp_operators):
+For the spatial discretization, we use sixth-order accurate [summation-by-parts operators](@ref sbp_operators):
 
 ```@example dingemans
-accuracy_order = 4
+accuracy_order = 6
 solver = Solver(mesh, accuracy_order)
 nothing # hide
 ```
@@ -114,7 +114,7 @@ sol_hysgn = solve(ode_hysgn, Tsit5(); options...)
 nothing # hide
 ```
 
-## Visualization and Comparison
+## Visualization of Temporal Evolution
 
 For proper comparison, we need to account for the fact that the BBM-BBM equations use a different reference level (``\eta_0 = 0``) compared to the other equations. We create a custom conversion function which allows us to easily shift the BBM-BBM results:
 
@@ -147,16 +147,16 @@ for time_val in times
 
     for (i, (semi, sol, label, conversion, linestyle)) in enumerate(models)
         plot!(p, semi => sol,
-            step=step_idx,
-            label=label,
-            conversion=conversion,
-            plot_bathymetry=true,
-            legend=false,
-            title="Dingemans at t = $(time_val)",
-            suptitle="",
-            linewidth=[2 1], # 1 for the bathymetry
-            linestyles=[linestyle :solid], # :solid for the bathymetry
-            color=[i :black], # black for the bathymetry
+              step=step_idx,
+              label=label,
+              conversion=conversion,
+              plot_bathymetry=true,
+              legend=false,
+              title="Dingemans at t = $(time_val)",
+              suptitle="",
+              linewidth=[2 1], # 1 for the bathymetry
+              linestyles=[linestyle :solid], # :solid for the bathymetry
+              color=[i :black], # black for the bathymetry
         )
     end
 
@@ -169,26 +169,80 @@ legend_plot = plot(legend=:top, framestyle=:none, legendfontsize=11)
 for (i, (_, _, label, _, linestyle)) in enumerate(models)
     plot!(legend_plot, [], [], label=label, linestyles=linestyle, linewidth=2, color=i)
 end
-plot!(legend_plot, [], [], label="Bathymetry", color=:black,)
+legend_plot_bathymetry = plot(legend_plot, [], [], label="Bathymetry", color=:black,)
 
 xlims_zoom = [(-25, 15), (0, 40), (5, 45), (-100, -60)]
 snapshot_plots_zoom = [plot(snapshot_plots[i], xlims=xlims_zoom[i], ylims=(0.75, 0.85), title="Zoomed in at t = $(times[i])") for i in 1:4]
 
 # Combine all plots
-all_plots = [snapshot_plots..., legend_plot, snapshot_plots_zoom...]
+all_plots = [snapshot_plots..., legend_plot_bathymetry, snapshot_plots_zoom...]
 plot(all_plots...,
-    size=(900, 1100),
-    layout=@layout([a b; c d; e{0.14h}; f g; h i]),
+     size=(900, 1100),
+     layout=@layout([a b; c d; e{0.14h}; f g; h i]),
 )
 
-savefig("dingemans_comparison.png") # hide
+savefig("dingemans_temporal.png") # hide
 nothing # hide
 ```
 
-![Dingemans comparison](dingemans_comparison.png)
+![Dingemans temporal](dingemans_temporal.png)
 
 The results show how different dispersive wave models capture the wave evolution over the trapezoidal bathymetry.
 
+## Comparison with Experimental Data
+
+During the experiment of Dingemans, the wave evolution was recorded at six gauges along the flume. These measurements provide a reference for validating the numerical models.
+We compare the numerical results with the experimental data at the corresponding gauge locations.
+
+```@example dingemans
+t_values, x_values, experimental_data = data_dingemans()
+
+tlims = [(20, 30), (25, 35), (30, 40), (35, 45), (40, 50), (45, 55)]
+
+snapshot_plots_time = []
+for (j, x_val) in enumerate(x_values)
+    p = plot(t_values, experimental_data[:, j],
+             xlims=tlims[j],
+             ylims=(0.765, 0.865),
+             label="experimental data",
+             linestyle=:dash,
+             color=:gray,
+             markershape=:circle,
+             markercolor=:gray,
+             markersize=1
+        )
+    for (i, (semi, sol, label, conversion, linestyle)) in enumerate(models)
+        plot!(p, semi => sol, x_val,
+              conversion=conversion,
+              label=label,
+              linestyle=linestyle,
+              color=i,
+              suptitle="",
+              legend=false,
+              title="Gauge at x = $(x_val)",
+              linewidth=2
+        )
+
+    end
+
+    push!(snapshot_plots_time, p)
+
+end
+
+legend_plot_data = plot(legend_plot, [], [], label="experimental data", linestyle=:dash, color=:gray, markershape=:circle, markercolor=:gray, markersize=1)
+
+all_plots2 = [snapshot_plots_time..., legend_plot_data]
+plot(all_plots2..., layout=@layout([a b; c d; e f; g{0.16h}]),
+     size=(900, 900), suptitle="")
+
+savefig("dingemans_experimental.png") # hide
+nothing # hide
+```
+
+![Dingemans experimental](dingemans_experimental.png)
+
+In this setup, the numerical results from the Svärd-Kalisch equations can capture the wave evolution more accurately compared to the other models.
+
 ## [Plain program](@id overview-plain-program-dingemans)
 
 Here follows a version of the program without any comments.
@@ -210,7 +264,7 @@ sgn = SerreGreenNaghdiEquations1D(bathymetry_type = bathymetry_variable,
 
 # Hyperbolic approximation of Serre-Green-Naghdi equations
 hysgn = HyperbolicSerreGreenNaghdiEquations1D(bathymetry_type = bathymetry_mild_slope,
-                                              lambda = 100.0, gravity = 9.81, eta0 = 0.8)
+                                              lambda = 200.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.
@@ -220,10 +274,10 @@ boundary_conditions = boundary_condition_periodic
 
 coordinates_min = -138.0
 coordinates_max = 46.0
-N = 512
+N = 1024
 mesh = Mesh1D(coordinates_min, coordinates_max, N)
 
-accuracy_order = 4
+accuracy_order = 6
 solver = Solver(mesh, accuracy_order)
 
 tspan = (0.0, 70.0)
@@ -273,16 +327,16 @@ for time_val in times
 
     for (i, (semi, sol, label, conversion, linestyle)) in enumerate(models)
         plot!(p, semi => sol,
-            step=step_idx,
-            label=label,
-            conversion=conversion,
-            plot_bathymetry=true,
-            legend=false,
-            title="Dingemans at t = $(time_val)",
-            suptitle="",
-            linewidth=[2 1], # 1 for the bathymetry
-            linestyles=[linestyle :solid], # :solid for the bathymetry
-            color=[i :black], # black for the bathymetry
+              step=step_idx,
+              label=label,
+              conversion=conversion,
+              plot_bathymetry=true,
+              legend=false,
+              title="Dingemans at t = $(time_val)",
+              suptitle="",
+              linewidth=[2 1], # 1 for the bathymetry
+              linestyles=[linestyle :solid], # :solid for the bathymetry
+              color=[i :black], # black for the bathymetry
         )
     end
 
@@ -295,17 +349,57 @@ legend_plot = plot(legend=:top, framestyle=:none, legendfontsize=11)
 for (i, (_, _, label, _, linestyle)) in enumerate(models)
     plot!(legend_plot, [], [], label=label, linestyles=linestyle, linewidth=2, color=i)
 end
-plot!(legend_plot, [], [], label="Bathymetry", color=:black,)
+legend_plot_bathymetry = plot(legend_plot, [], [], label="Bathymetry", color=:black,)
 
 xlims_zoom = [(-25, 15), (0, 40), (5, 45), (-100, -60)]
 snapshot_plots_zoom = [plot(snapshot_plots[i], xlims=xlims_zoom[i], ylims=(0.75, 0.85), title="Zoomed in at t = $(times[i])") for i in 1:4]
 
 # Combine all plots
-all_plots = [snapshot_plots..., legend_plot, snapshot_plots_zoom...]
+all_plots = [snapshot_plots..., legend_plot_bathymetry, snapshot_plots_zoom...]
 plot(all_plots...,
-    size=(900, 1100),
-    layout=@layout([a b; c d; e{0.14h}; f g; h i]),
+     size=(900, 1100),
+     layout=@layout([a b; c d; e{0.14h}; f g; h i]),
 )
+
+t_values, x_values, experimental_data = data_dingemans()
+
+tlims = [(20, 30), (25, 35), (30, 40), (35, 45), (40, 50), (45, 55)]
+
+snapshot_plots_time = []
+for (j, x_val) in enumerate(x_values)
+    p = plot(t_values, experimental_data[:, j],
+             xlims=tlims[j],
+             ylims=(0.765, 0.865),
+             label="experimental data",
+             linestyle=:dash,
+             color=:gray,
+             markershape=:circle,
+             markercolor=:gray,
+             markersize=1
+        )
+    for (i, (semi, sol, label, conversion, linestyle)) in enumerate(models)
+        plot!(p, semi => sol, x_val,
+              conversion=conversion,
+              label=label,
+              linestyle=linestyle,
+              color=i,
+              suptitle="",
+              legend=false,
+              title="Gauge at x = $(x_val)",
+              linewidth=2
+        )
+
+    end
+
+    push!(snapshot_plots_time, p)
+
+end
+
+legend_plot_data = plot(legend_plot, [], [], label="experimental data", linestyle=:dash, color=:gray, markershape=:circle, markercolor=:gray, markersize=1)
+
+all_plots2 = [snapshot_plots_time..., legend_plot_data]
+plot(all_plots2..., layout=@layout([a b; c d; e f; g{0.16h}]),
+     size=(900, 900), suptitle="")
 ```
 
 ### References

src/DispersiveShallowWater.jl

@@ -16,6 +16,7 @@ See also: [DispersiveShallowWater.jl](https://github.com/NumericalMathematics/Di
 module DispersiveShallowWater
 
 using BandedMatrices: BandedMatrix
+using DelimitedFiles: readdlm
 using DiffEqBase: DiffEqBase, terminate!
 using FastBroadcast: @..
 using ForwardDiff: ForwardDiff
@@ -56,6 +57,7 @@ using TimerOutputs: TimerOutputs, print_timer, reset_timer!
 @reexport using TrixiBase: trixi_include
 using TrixiBase: TrixiBase, @trixi_timeit, timer
 
+include("experimental_data.jl")
 include("boundary_conditions.jl")
 include("mesh.jl")
 include("equations/equations.jl")
@@ -66,7 +68,9 @@ include("callbacks_step/callbacks_step.jl")
 include("visualization.jl")
 include("util.jl")
 
-export examples_dir, get_examples, default_example, convergence_test
+export examples_dir, get_examples, default_example, convergence_test, data_dir
+
+export data_dingemans
 
 export AbstractShallowWaterEquations,
        KdVEquation1D,

src/experimental_data.jl

@@ -0,0 +1,21 @@
+"""
+    data_dingemans()
+
+Load the experimental data for the Dingemans experiment. Returns the time values, x-coordinates of the six wave gauges,
+and experimental data, which is a matrix of wave heights at each gauge location (columns) over time (rows).
+
+- Dingemans (1994)
+  Comparison of computations with Boussinesq-like models and laboratory measurements
+  [URL: https://resolver.tudelft.nl/uuid:c2091d53-f455-48af-a84b-ac86680455e9](https://resolver.tudelft.nl/uuid:c2091d53-f455-48af-a84b-ac86680455e9)
+- Dingemans (1997):
+  Water Wave Propagation Over Uneven Bottoms (In 2 Parts).
+  [DOI: 10.1142/1241](https://doi.org/10.1142/1241)
+"""
+function data_dingemans()
+    path_dingemans = joinpath(data_dir(), "Dingemans.csv")
+    all_data, _ = readdlm(path_dingemans, ','; header = true)
+    t_values = all_data[:, 1]
+    x_values = (3.04, 9.44, 20.04, 26.04, 30.44, 37.04)
+    experimental_data = all_data[:, 2:end]
+    return t_values, x_values, experimental_data
+end

src/util.jl

@@ -47,6 +47,19 @@ function default_example()
     joinpath(examples_dir(), "bbm_bbm_1d", "bbm_bbm_1d_basic.jl")
 end
 
+"""
+    data_dir()
+
+Return the directory where the data files provided with DispersiveShallowWater.jl are located. If DispersiveShallowWater is
+installed as a regular package (with `]add DispersiveShallowWater`), these files are read-only and should *not* be
+modified. To find out which files are available, use, e.g., `readdir`.
+
+```@example
+readdir(data_dir())
+```
+"""
+data_dir() = pkgdir(DispersiveShallowWater, "data")::String
+
 function convergence_test(example::AbstractString, iterations_or_Ns; kwargs...)
     convergence_test(Main, example::AbstractString, iterations_or_Ns; kwargs...)
 end

test/test_unit.jl

@@ -532,7 +532,11 @@ end
 end
 
 @testitem "util" setup=[Setup] begin
+    @test_nowarn examples_dir()
     @test_nowarn get_examples()
+    @test_nowarn data_dir()
+    t, x, experimental_data = data_dingemans()
+    @test size(experimental_data) == (length(t), length(x))
 
     accuracy_orders = [2, 4, 6]
     for accuracy_order in accuracy_orders