Dingemans

Author: ranocha

.gitignore

@@ -8,3 +8,4 @@ out*/
 /docs/build/
 run
 run/*
+**/*.code-workspace

examples/serre_green_naghdi_1d/serre_green_naghdi_dingemans.jl

@@ -0,0 +1,49 @@
+using OrdinaryDiffEq
+using DispersiveShallowWater
+using SummationByPartsOperators: upwind_operators, periodic_derivative_operator
+
+###############################################################################
+# Semidiscretization of the Svärd-Kalisch equations
+
+bathymetry = bathymetry_variable # or bathymetry_mild_slope
+equations = SerreGreenNaghdiEquations1D(bathymetry;
+                                        gravity_constant = 9.81)
+
+initial_condition = initial_condition_dingemans
+boundary_conditions = boundary_condition_periodic
+
+# create homogeneous mesh
+coordinates_min = -140.0
+coordinates_max = 100.0
+N = 512
+mesh = Mesh1D(coordinates_min, coordinates_max, N)
+
+# create solver with periodic upwind SBP operators of accuracy order 4
+# (note that only central-type with even order are used)
+D1 = upwind_operators(periodic_derivative_operator;
+                      derivative_order = 1, accuracy_order = 3,
+                      xmin = xmin(mesh), xmax = xmax(mesh),
+                      N = nnodes(mesh))
+solver = Solver(D1, nothing)
+
+# semidiscretization holds all the necessary data structures for the spatial discretization
+semi = Semidiscretization(mesh, equations, initial_condition, solver,
+                          boundary_conditions = boundary_conditions)
+
+###############################################################################
+# Create `ODEProblem` and run the simulation
+tspan = (0.0, 70.0)
+ode = semidiscretize(semi, tspan)
+summary_callback = SummaryCallback()
+analysis_callback = AnalysisCallback(semi; interval = 10,
+                                     extra_analysis_errors = (:conservation_error,),
+                                     extra_analysis_integrals = (waterheight_total,
+                                                                 momentum,
+                                                                 entropy,
+                                                                 entropy_modified))
+
+callbacks = CallbackSet(analysis_callback, summary_callback)
+
+saveat = range(tspan..., length = 500)
+sol = solve(ode, Tsit5(), abstol = 1e-7, reltol = 1e-7,
+            save_everystep = false, callback = callbacks, saveat = saveat)

src/equations/serre_green_naghdi_1d.jl

@@ -77,11 +77,19 @@ end
 
 function SerreGreenNaghdiEquations1D(bathymetry = bathymetry_flat;
                                      gravity_constant, eta0 = 0.0)
-    eta0 == 0.0 ||
-        @warn "The still-water surface needs to be 0 for the Serre-Green-Naghdi equations"
     SerreGreenNaghdiEquations1D(bathymetry, gravity_constant, eta0)
 end
 
+function get_name(equations::SerreGreenNaghdiEquations1D)
+    name = equations |> typeof |> nameof |> string
+    bathymetry = equations.bathymetry
+    if bathymetry isa BathymetryFlat
+        return name
+    else # variable bathymetry
+        return name * "-" * string(nameof(typeof(bathymetry)))
+    end
+end
+
 varnames(::typeof(prim2prim), ::SerreGreenNaghdiEquations1D) = ("η", "v", "D")
 varnames(::typeof(prim2cons), ::SerreGreenNaghdiEquations1D) = ("h", "hv", "b")
 
@@ -107,6 +115,46 @@ function initial_condition_convergence_test(x, t, equations::SerreGreenNaghdiEqu
     return SVector(h, v, 0)
 end
 
+"""
+    initial_condition_dingemans(x, t, equations::SerreGreenNaghdiEquations1D, mesh)
+
+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. It is assumed that `equations.eta0 = 0.8`.
+
+References:
+- Magnus Svärd, Henrik Kalisch (2023)
+  A novel energy-bounded Boussinesq model and a well-balanced and stable numerical discretization
+  [arXiv: 2302.09924](https://arxiv.org/abs/2302.09924)
+- Maarten W. Dingemans (1994)
+  Comparison of computations with Boussinesq-like models and laboratory measurements
+  [link](https://repository.tudelft.nl/islandora/object/uuid:c2091d53-f455-48af-a84b-ac86680455e9/datastream/OBJ/download)
+"""
+function initial_condition_dingemans(x, t, equations::SerreGreenNaghdiEquations1D, mesh)
+    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 * h0), 1.0) using Roots.jl
+    if x < -34.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 * 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
+        b = 0.6
+    elseif 27.04 <= x && x < 33.07
+        b = 0.6 * (33.07 - x) / (33.07 - 27.04)
+    else
+        b = 0.0
+    end
+    eta = h + h0
+    D = equations.eta0 - b
+    return SVector(eta, v, D)
+end
+
 # flat bathymetry
 function create_cache(mesh,
                       equations::SerreGreenNaghdiEquations1D{BathymetryFlat},
@@ -245,15 +293,15 @@ function create_cache(mesh,
 
         factorization = cholesky(system_matrix)
 
-        cache = (; h_x, v_x, v_x_upwind, h_hpb_x, b, b_x, hv_x, v2_x,
+        cache = (; h, h_x, v_x, v_x_upwind, h_hpb_x, b, b_x, hv_x, v2_x,
                    h2_v_vx_x, h_vx_x, p_h, p_0, p_x, tmp,
                    M_h_p_h_bx2, M_h3_3, M_h2_bx,
                    D, Dmat_minus, factorization)
     else
         if D isa FourierDerivativeOperator
             Dmat = Matrix(D)
 
-            cache = (; h_x, v_x, h_hpb_x, b, b_x, hv_x, v2_x,
+            cache = (; h, h_x, v_x, h_hpb_x, b, b_x, hv_x, v2_x,
                        h2_v_vx_x, h_vx_x, p_h, p_x, tmp,
                        M_h_p_h_bx2, M_h3_3, M_h2_bx,
                        D, Dmat)
@@ -285,7 +333,7 @@ function create_cache(mesh,
 
             factorization = cholesky(system_matrix)
 
-            cache = (; h_x, v_x, h_hpb_x, b, b_x, hv_x, v2_x,
+            cache = (; h, h_x, v_x, h_hpb_x, b, b_x, hv_x, v2_x,
                        h2_v_vx_x, h_vx_x, p_h, p_x, tmp,
                        M_h_p_h_bx2, M_h3_3, M_h2_bx,
                        D, Dmat, factorization)
@@ -531,20 +579,21 @@ function rhs_sgn_central!(dq, q, equations, source_terms, cache,
 
     # `q` and `dq` are `ArrayPartition`s. They collect the individual
     # arrays for the water height `h` and the velocity `v`.
-    h, v = q.x
+    eta, v = q.x
     dh, dv, dD = dq.x
     fill!(dD, zero(eltype(dD)))
 
     @trixi_timeit timer() "hyperbolic terms" begin
         # Compute all derivatives required below
-        (; h_x, v_x, h_hpb_x, b, b_x, hv_x, v2_x,
+        (; h, h_x, v_x, h_hpb_x, b, b_x, hv_x, v2_x,
           h2_v_vx_x, h_vx_x, p_h, p_x, tmp,
           M_h_p_h_bx2, M_h3_3, M_h2_bx) = cache
         if equations.bathymetry isa BathymetryVariable
             (; ψ) = cache
         end
 
         @. b = equations.eta0 - q.x[3]
+        @. h = eta - b
         mul!(b_x, D, b)
 
         mul!(h_x, D, h)
@@ -662,20 +711,21 @@ function rhs_sgn_upwind!(dq, q, equations, source_terms, cache,
 
     # `q` and `dq` are `ArrayPartition`s. They collect the individual
     # arrays for the water height `h` and the velocity `v`.
-    h, v, b = q.x
+    eta, v, b = q.x
     dh, dv, db = dq.x
     fill!(db, zero(eltype(db)))
 
     @trixi_timeit timer() "hyperbolic terms" begin
         # Compute all derivatives required below
-        (; h_x, v_x, v_x_upwind, h_hpb_x, b, b_x, hv_x, v2_x,
+        (; h, h_x, v_x, v_x_upwind, h_hpb_x, b, b_x, hv_x, v2_x,
           h2_v_vx_x, h_vx_x, p_h, p_0, p_x, tmp,
           M_h_p_h_bx2, M_h3_3, M_h2_bx) = cache
         if equations.bathymetry isa BathymetryVariable
             (; ψ) = cache
         end
 
         @. b = equations.eta0 - q.x[3]
+        @. h = eta - b
         mul!(b_x, D, b)
 
         mul!(h_x, D, h)
@@ -856,14 +906,31 @@ function energy_total_modified(q_global,
     N = length(v)
     e = zeros(eltype(q_global), N)
 
-    # 1/2 g h^2 + 1/2 h v^2 + 1/6 h^3 v_x^2
+    # 1/2 g h^2 + 1/2 h v^2 + 1/6 h^3 w^2
+    # and + 1/8 h (v b_x)^2 for full bathymetry without mild-slope approximation
     if D isa PeriodicUpwindOperators
         mul!(v_x, D.minus, v)
     else
         mul!(v_x, D, v)
     end
 
-    @. e = 1 / 2 * g * h^2 + 1 / 2 * h * v^2 + 1 / 6 * h^3 * v_x^2
+    if equations.bathymetry isa BathymetryFlat
+        b_x = cache.tmp
+        fill!(b_x, zero(eltype(b_x)))
+    else
+        (; b, b_x) = cache
+        @. b = equations.eta0 - q_global.x[3]
+        if D isa PeriodicUpwindOperators
+            mul!(b_x, D.central, b)
+        else
+            mul!(b_x, D, b)
+        end
+    end
+
+    @. e = 1 / 2 * g * h^2 + 1 / 2 * h * v^2 + 1 / 6 * h * (-h * v_x + 1.5 * v * b_x)^2
+    if equations.bathymetry isa BathymetryVariable
+        @. e += 1 / 8 * h * (v * b_x)^2
+    end
 
     return e
 end

test/test_serre_green_naghdi_1d.jl

@@ -179,6 +179,35 @@ EXAMPLES_DIR = joinpath(examples_dir(), "serre_green_naghdi_1d")
 
         @test_allocations(semi, sol, allocs=850_000)
     end
+
+    @trixi_testset "serre_green_naghdi_dingemans.jl" begin
+        @test_trixi_include(joinpath(EXAMPLES_DIR,
+                                     "serre_green_naghdi_dingemans.jl"),
+                            tspan=(0.0, 1.0),
+                            l2=[0.2463295492331802, 0.805357513755897, 0.0],
+                            linf=[0.036265772921736605, 0.11763218152350845, 0.0],
+                            cons_error=[5.684341886080802e-14, 3.509473432346808e-5, 0.0],
+                            change_waterheight=5.684341886080802e-14,
+                            change_entropy=1.9489684518703143e-5,
+                            change_entropy_modified=-1.2177679309388623e-8)
+
+        @test_allocations(semi, sol, allocs=750_000)
+    end
+
+    @trixi_testset "serre_green_naghdi_dingemans.jl with bathymetry_mild_slope" begin
+        @test_trixi_include(joinpath(EXAMPLES_DIR,
+                                     "serre_green_naghdi_dingemans.jl"),
+                            tspan=(0.0, 1.0),
+                            equations = SerreGreenNaghdiEquations1D(bathymetry_mild_slope; gravity_constant = 9.81),
+                            l2=[0.2463295492331802, 0.805357513755897, 0.0],
+                            linf=[0.036265772921736605, 0.11763218152350845, 0.0],
+                            cons_error=[5.684341886080802e-14, 3.509473432346624e-5, 0.0],
+                            change_waterheight=5.684341886080802e-14,
+                            change_entropy=1.9489684518703143e-5,
+                            change_entropy_modified=-1.2177679309388623e-8)
+
+        @test_allocations(semi, sol, allocs=750_000)
+    end
 end
 
 end # module