Add BBMEquation1D

.gitignore

@@ -9,3 +9,4 @@ out*/
 run
 run/*
 **/*.code-workspace
+.vscode/

examples/bbm_1d/bbm_1d_basic.jl

@@ -2,9 +2,9 @@ using OrdinaryDiffEq
 using DispersiveShallowWater
 
 ###############################################################################
-# Semidiscretization of the BBM equation
+# Semidiscretization of the BBM equation (conserves the quadratic invariant)
 
-equations = BBMEquation1D()
+equations = BBMEquation1D(split_form = true)
 
 initial_condition = initial_condition_convergence_test
 boundary_conditions = boundary_condition_periodic
@@ -32,7 +32,7 @@ analysis_callback = AnalysisCallback(semi; interval = 10,
                                      extra_analysis_errors = (:conservation_error,),
                                      extra_analysis_integrals = (waterheight_total,
                                                                  entropy_modified,
-                                                                 invariant_cubic))
+                                                                 hamiltonian))
 callbacks = CallbackSet(analysis_callback, summary_callback)
 
 saveat = range(tspan..., length = 100)

examples/bbm_1d/bbm_1d_fourier.jl

@@ -34,7 +34,7 @@ analysis_callback = AnalysisCallback(semi; interval = 10,
                                      extra_analysis_errors = (:conservation_error,),
                                      extra_analysis_integrals = (waterheight_total,
                                                                  entropy_modified,
-                                                                 invariant_cubic))
+                                                                 hamiltonian))
 callbacks = CallbackSet(analysis_callback, summary_callback)
 
 saveat = range(tspan..., length = 100)

examples/bbm_1d/bbm_1d_hamiltonian_relaxation.jl

@@ -0,0 +1,43 @@
+using OrdinaryDiffEq
+using DispersiveShallowWater
+
+###############################################################################
+# Semidiscretization of the BBM equation (conserves the quadratic invariant)
+
+equations = BBMEquation1D(split_form = false)
+
+initial_condition = initial_condition_convergence_test
+boundary_conditions = boundary_condition_periodic
+
+# create homogeneous mesh
+coordinates_min = -90.0
+coordinates_max = 90.0
+N = 512
+mesh = Mesh1D(coordinates_min, coordinates_max, N)
+
+# create solver with periodic SBP operators of accuracy order 4
+accuracy_order = 4
+solver = Solver(mesh, accuracy_order)
+
+# 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, 1000.0)
+ode = semidiscretize(semi, tspan)
+summary_callback = SummaryCallback()
+analysis_callback = AnalysisCallback(semi; interval = 10,
+                                     extra_analysis_errors = (:conservation_error,),
+                                     extra_analysis_integrals = (waterheight_total,
+                                                                 entropy_modified,
+                                                                 hamiltonian))
+
+relaxation_callback = RelaxationCallback(invariant = hamiltonian)
+# Always put relaxation_callback before analysis_callback to guarantee conservation of the invariant
+callbacks = CallbackSet(relaxation_callback, analysis_callback, summary_callback)
+
+saveat = range(tspan..., length = 100)
+sol = solve(ode, Tsit5(), abstol = 1e-7, reltol = 1e-7,
+            save_everystep = false, callback = callbacks, saveat = saveat)

examples/bbm_1d/bbm_1d_relaxation.jl

@@ -32,7 +32,7 @@ analysis_callback = AnalysisCallback(semi; interval = 10,
                                      extra_analysis_errors = (:conservation_error,),
                                      extra_analysis_integrals = (waterheight_total,
                                                                  entropy_modified,
-                                                                 invariant_cubic))
+                                                                 hamiltonian))
 
 relaxation_callback = RelaxationCallback(invariant = entropy_modified)
 # Always put relaxation_callback before analysis_callback to guarantee conservation of the invariant

src/DispersiveShallowWater.jl

@@ -72,7 +72,7 @@ export prim2prim, prim2cons, cons2prim, prim2phys,
        bathymetry, still_water_surface,
        energy_total, entropy, lake_at_rest_error,
        energy_total_modified, entropy_modified,
-       invariant_cubic
+       hamiltonian
 
 export Mesh1D, xmin, xmax, nnodes
 

src/callbacks_step/analysis.jl

@@ -339,4 +339,4 @@ pretty_form_utf(::typeof(energy_total)) = "∫e_total"
 pretty_form_utf(::typeof(entropy_modified)) = "∫U_modified"
 pretty_form_utf(::typeof(energy_total_modified)) = "∫e_modified"
 pretty_form_utf(::typeof(lake_at_rest_error)) = "∫|η-η₀|"
-pretty_form_utf(::typeof(invariant_cubic)) = "∫(η + 1)³"
+pretty_form_utf(::typeof(hamiltonian)) = "∫H"

src/equations/bbm_1d.jl

@@ -1,6 +1,6 @@
 # TODO: Use this nondimensionalized version or use dimensional version including gravity and depth?
 @doc raw"""
-    BBMEquation1D(; bathymetry_type = bathymetry_flat, eta0 = 0.0)
+    BBMEquation1D(; bathymetry_type = bathymetry_flat, eta0 = 0.0, split_form = true)
 
 BBM (Benjamin–Bona–Mahony) equation in one spatial dimension. The equations are given by
 ```math
@@ -15,12 +15,17 @@
 Currently, the equations only support a flat bathymetry ``b = 0``, see [`bathymetry_flat`](@ref).
 
 The BBM equation is first described in Benjamin, Bona, and Mahony (1972).
-The semidiscretization implemented here is developed in Ranocha, Mitsotakis, and Ketcheson (2020).
-It conserves
+The semidiscretization implemented here is developed in Ranocha, Mitsotakis, and Ketcheson (2020)
+for `split_form = true` and in Linders, Ranocha, and Birken (2023) for `split_form = false`.
+If `split_form` is `true` a split form in the semidiscretization is used, which conserves
 - the total water mass (integral of ``h``) as a linear invariant
 - a quadratic invariant (integral of ``\eta^2 + \eta_x^2``), which is called here [`energy_total_modified`](@ref)
   (and [`entropy_modified`](@ref)) because it contains derivatives of the solution
 
+for periodic boundary conditions. If `split_form` is `false` the semidiscretization conserves
+- the total water mass (integral of ``h``) as a linear invariant
+- the Hamiltonian (integral of ``1/6 \eta^3 + 1/2 \eta^2``) (see [`hamiltonian`](@ref))
+
 for periodic boundary conditions.
 
 - Thomas B. Benjamin, Jerry L. Bona and John J. Mahony (1972)
@@ -29,23 +34,27 @@
 - Hendrik Ranocha, Dimitrios Mitsotakis and David I. Ketcheson (2020)
   A Broad Class of Conservative Numerical Methods for Dispersive Wave Equations
   [DOI: 10.4208/cicp.OA-2020-0119](https://doi.org/10.4208/cicp.OA-2020-0119)
+- Viktor Linders, Hendrik Ranocha and Philipp Birken (2023)
+  Resolving entropy growth from iterative methods
+  [DOI: 10.1007/s10543-023-00992-w](https://doi.org/10.1007/s10543-023-00992-w)
 """
 struct BBMEquation1D{Bathymetry <: AbstractBathymetry, RealT <: Real} <:
        AbstractBBMEquation{1, 1}
     bathymetry_type::Bathymetry # type of bathymetry
     eta0::RealT # constant still-water surface
+    split_form::Bool # whether to use a split-form or not
 end
 
-function BBMEquation1D(; bathymetry_type = bathymetry_flat, eta0 = 0.0)
-    BBMEquation1D(bathymetry_type, eta0)
+function BBMEquation1D(; bathymetry_type = bathymetry_flat, eta0 = 0.0, split_form = true)
+    BBMEquation1D(bathymetry_type, eta0, split_form)
 end
 
 """
     initial_condition_convergence_test(x, t, equations::BBMEquation1D, mesh)
 
 A travelling-wave solution used for convergence tests in a periodic domain.
 
-See section 4.1.3 in (there is an error in paper: it should be sech^2 instead of cosh):
+See section 4.1.3 in (there is an error in paper: it should be `sech^2` instead of `cosh`):
 - Hendrik Ranocha, Dimitrios Mitsotakis and David I. Ketcheson (2020)
   A Broad Class of Conservative Numerical Methods for Dispersive Wave Equations
   [DOI: 10.4208/cicp.OA-2020-0119](https://doi.org/10.4208/cicp.OA-2020-0119)
@@ -126,18 +135,24 @@
            solver.D1 isa FourierDerivativeOperator
             mul!(eta2_x, solver.D1, eta2)
             mul!(eta_x, solver.D1, eta)
-            @.. etaeta_x = eta * eta_x
-            # split form
-            @.. deta = -(1 / 3 * (eta2_x + etaeta_x) + eta_x)
+            if equations.split_form
+                @.. etaeta_x = eta * eta_x
+                @.. deta = -(1 / 3 * (eta2_x + etaeta_x) + eta_x)
+            else
+                @.. deta = -(0.5 * eta2_x + eta_x)
+            end
         elseif solver.D1 isa PeriodicUpwindOperators
             mul!(eta2_x, solver.D1.central, eta2)
             mul!(eta_x_upwind, solver.D1.minus, eta)
             mul!(eta_x, solver.D1.central, eta)
-            @.. etaeta_x = eta * eta_x
-            # split form
-            @.. deta = -(1 / 3 * (eta2_x + etaeta_x) + eta_x_upwind)
+            if equations.split_form
+                @.. etaeta_x = eta * eta_x
+                @.. deta = -(1 / 3 * (eta2_x + etaeta_x) + eta_x_upwind)
+            else
+                @.. deta = -(0.5 * eta2_x + eta_x_upwind)
+            end
         else
             @error "unknown type of first-derivative operator: $(typeof(solver.D1))"
         end
     end
 
@@ -179,21 +194,27 @@
     end
 
     @.. e = 0.5 * eta * (eta - eta_xx)
-    return nothing
+    return e
 end
 
 """
-    invariant_cubic(q, equations::BBMEquation1D)
+    hamiltonian!(H, q_global, equations::BBMEquation1D, cache)
 
-Return the cubic invariant of the primitive variables `q` for the
-[`BBMEquation1D`](@ref). The cubic invariant is given by
+Return the Hamiltonian `H` of the primitive variables `q_global` for the
+[`BBMEquation1D`](@ref). The Hamiltonian is given by
 ```math
-(\\eta + 1)^3.
+\\frac{1}{6}\\eta^3 + \\frac{1}{2}\\eta^2.
 ```
+
+`q_global` is a vector of the primitive variables at ALL nodes.
+
+See also [`hamiltonian`](@ref).
 """
-function invariant_cubic(q, equations::BBMEquation1D)
-    eta = waterheight_total(q, equations)
-    return (eta + 1)^3
+function hamiltonian!(H, q_global, equations::BBMEquation1D, cache)
+    eta, = q_global.x
+    (; eta2, tmp1) = cache
+    @.. eta2 = eta^2
+    @.. tmp1 = eta^3
+    @.. H = 1 / 6 * tmp1 + 1 / 2 * eta2
+    return H
 end
-
-varnames(::typeof(invariant_cubic), equations) = ("(η + 1)³",)

src/equations/equations.jl

@@ -356,7 +356,7 @@ varnames(::typeof(energy_total_modified), equations) = ("e_modified",)
 
 Alias for [`energy_total_modified`](@ref).
 """
-@inline function entropy_modified(q_global, equations::AbstractShallowWaterEquations, cache)
+@inline function entropy_modified(q_global, equations::AbstractEquations, cache)
     e = similar(q_global.x[begin])
     return entropy_modified!(e, q_global, equations, cache)
 end
@@ -373,6 +373,28 @@ In-place version of [`entropy_modified`](@ref).
 
 varnames(::typeof(entropy_modified), equations) = ("U_modified",)
 
+# The Hamiltonian takes the whole `q_global` for every point in space
+"""
+    hamiltonian(q_global, equations, cache)
+
+Return the Hamiltonian of the primitive variables `q_global` for the
+`equations`. The Hamiltonian is a conserved quantity and may contain
+derivatives of the solution.
+
+`q_global` is a vector of the primitive variables at ALL nodes.
+`cache` needs to hold the SBP operators used by the `solver` if non-hydrostatic
+terms are present.
+
+Internally, this function allocates a vector for the output and
+calls [`DispersiveShallowWater.hamiltonian!`](@ref).
+"""
+function hamiltonian(q_global, equations::AbstractEquations, cache)
+    H = similar(q_global.x[begin])
+    return hamiltonian!(H, q_global, equations, cache)
+end
+
+varnames(::typeof(hamiltonian), equations) = ("H",)
+
 # Add methods to show some information on systems of equations.
 function Base.show(io::IO, equations::AbstractEquations)
     # Since this is not performance-critical, we can use `@nospecialize` to reduce latency.

src/semidiscretization.jl

@@ -147,12 +147,13 @@ end
 # Obtain the function, which has an additional `!` appended to the name
 inplace_version(f) = getfield(@__MODULE__, Symbol(string(nameof(f)) * "!"))
 
-# The entropy/energy of the Svärd-Kalisch and Serre-Green-Naghdi equations
-# takes the whole `q` for every point in space since it requires
+# The modified entropy/energy and the Hamiltonian
+# take the whole `q` for every point in space since it requires
 # the derivative of the velocity `v_x`.
 function integrate_quantity!(quantity,
                              func::Union{typeof(energy_total_modified),
-                                         typeof(entropy_modified)}, q,
+                                         typeof(entropy_modified),
+                                         typeof(hamiltonian)}, q,
                              semi::Semidiscretization)
     inplace_version(func)(quantity, q, semi.equations, semi.cache)
     integrate(quantity, semi)

test/test_bbm_1d.jl

@@ -16,7 +16,7 @@ EXAMPLES_DIR = joinpath(examples_dir(), "bbm_1d")
                             cons_error=[7.105427357601002e-15],
                             change_waterheight=-7.105427357601002e-15,
                             change_entropy_modified=-4.4274395039067826e-7,
-                            change_invariant_cubic=-3.198066679033218e-6)
+                            change_hamiltonian=-5.33011109249415e-7)
 
         @test_allocations(semi, sol, allocs=5_000)
 
@@ -52,6 +52,51 @@ EXAMPLES_DIR = joinpath(examples_dir(), "bbm_1d")
         @test_allocations(semi, sol, allocs=5_000)
     end
 
+    @trixi_testset "bbm_1d_basic with split_form = false" begin
+        @test_trixi_include(joinpath(EXAMPLES_DIR, "bbm_1d_basic.jl"),
+                            tspan=(0.0, 100.0),
+                            split_form=false,
+                            l2=[0.00023679720962102973],
+                            linf=[9.291484780826753e-5],
+                            cons_error=[0.0],
+                            change_waterheight=0.0,
+                            change_entropy_modified=-4.4114923292148944e-7,
+                            change_hamiltonian=-5.313053028643822e-7)
+
+        @test_allocations(semi, sol, allocs=5_000)
+
+        # test upwind operators
+        using SummationByPartsOperators: upwind_operators, periodic_derivative_operator
+        using SparseArrays: sparse
+        using OrdinaryDiffEq: solve
+        D1 = upwind_operators(periodic_derivative_operator; derivative_order = 1,
+                              accuracy_order = accuracy_order, xmin = mesh.xmin,
+                              xmax = mesh.xmax,
+                              N = mesh.N)
+        D2 = sparse(D1.minus) * sparse(D1.plus)
+        solver = Solver(D1, D2)
+        semi = Semidiscretization(mesh, equations, initial_condition, solver,
+                                  boundary_conditions = boundary_conditions)
+        ode = semidiscretize(semi, (0.0, 100.0))
+        sol = solve(ode, Tsit5(), abstol = 1e-7, reltol = 1e-7,
+                    save_everystep = false, callback = callbacks, saveat = saveat)
+        atol = 1e-12
+        rtol = 1e-12
+        errs = errors(analysis_callback)
+        l2 = [0.00016247570402673777]
+        l2_measured = errs.l2_error[:, end]
+        for (l2_expected, l2_actual) in zip(l2, l2_measured)
+            @test isapprox(l2_expected, l2_actual, atol = atol, rtol = rtol)
+        end
+        linf = [6.495267031642049e-5]
+        linf_measured = errs.linf_error[:, end]
+        for (linf_expected, linf_actual) in zip(linf, linf_measured)
+            @test isapprox(linf_expected, linf_actual, atol = atol, rtol = rtol)
+        end
+
+        @test_allocations(semi, sol, allocs=5_000)
+    end
+
     @trixi_testset "bbm_1d_relaxation" begin
         @test_trixi_include(joinpath(EXAMPLES_DIR, "bbm_1d_relaxation.jl"),
                             tspan=(0.0, 100.0),
@@ -60,7 +105,20 @@ EXAMPLES_DIR = joinpath(examples_dir(), "bbm_1d")
                             cons_error=[5.329070518200751e-15],
                             change_waterheight=5.329070518200751e-15,
                             change_entropy_modified=0.0,
-                            change_invariant_cubic=-1.0308752962373546e-8)
+                            change_hamiltonian=-1.7181174261082788e-9)
+
+        @test_allocations(semi, sol, allocs=5_000)
+    end
+
+    @trixi_testset "bbm_1d_hamiltonian_relaxation" begin
+        @test_trixi_include(joinpath(EXAMPLES_DIR, "bbm_1d_hamiltonian_relaxation.jl"),
+                            tspan=(0.0, 100.0),
+                            l2=[0.00023637369561147064],
+                            linf=[9.280778692660752e-5],
+                            cons_error=[1.7763568394002505e-14],
+                            change_waterheight=-1.7763568394002505e-14,
+                            change_entropy_modified=1.6049757078917537e-9,
+                            change_hamiltonian=0.0)
 
         @test_allocations(semi, sol, allocs=5_000)
     end
@@ -73,7 +131,7 @@ EXAMPLES_DIR = joinpath(examples_dir(), "bbm_1d")
                             cons_error=[1.1546319456101628e-14],
                             change_waterheight=-1.1546319456101628e-14,
                             change_entropy_modified=-4.4266038123907947e-7,
-                            change_invariant_cubic=-3.1871549026618595e-6)
+                            change_hamiltonian=-5.311924822226644e-7)
 
         @test_allocations(semi, sol, allocs=5_000)
     end