Add reflective boundary conditions for BBMBBMEquations1D

Project.toml

@@ -4,6 +4,7 @@ authors = ["Joshua Lampert <[email protected]>"]
 version = "0.3.3-pre"
 
 [deps]
+BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
 Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -21,6 +22,7 @@ TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
 TrixiBase = "9a0f1c46-06d5-4909-a5a3-ce25d3fa3284"
 
 [compat]
+BandedMatrices = "0.17, 1"
 DiffEqBase = "6.128"
 Interpolations = "0.14.2, 0.15"
 LinearAlgebra = "1"

examples/bbm_bbm_1d/bbm_bbm_1d_basic_reflecting.jl

@@ -0,0 +1,48 @@
+using OrdinaryDiffEq
+using DispersiveShallowWater
+using SummationByPartsOperators: MattssonNordström2004, derivative_operator
+
+###############################################################################
+# Semidiscretization of the BBM-BBM equations
+
+equations = BBMBBMEquations1D(gravity_constant = 9.81, D = 1.0)
+
+initial_condition = initial_condition_manufactured_reflecting
+source_terms = source_terms_manufactured_reflecting
+boundary_conditions = boundary_condition_reflecting
+
+# create homogeneous mesh
+coordinates_min = 0.0
+coordinates_max = 1.0
+N = 512
+mesh = Mesh1D(coordinates_min, coordinates_max, N)
+
+# create solver with SBP operators of accuracy order 4
+accuracy_order = 4
+D1 = derivative_operator(MattssonNordström2004(),
+                         derivative_order = 1, accuracy_order = accuracy_order,
+                         xmin = mesh.xmin, xmax = mesh.xmax, N = mesh.N)
+D2 = derivative_operator(MattssonNordström2004(),
+                         derivative_order = 2, accuracy_order = accuracy_order,
+                         xmin = mesh.xmin, xmax = mesh.xmax, N = mesh.N)
+solver = Solver(D1, D2)
+
+# semidiscretization holds all the necessary data structures for the spatial discretization
+semi = Semidiscretization(mesh, equations, initial_condition, solver,
+                          boundary_conditions = boundary_conditions,
+                          source_terms = source_terms)
+
+###############################################################################
+# Create `ODEProblem` and run the simulation
+tspan = (0.0, 1.0)
+ode = semidiscretize(semi, tspan)
+summary_callback = SummaryCallback()
+analysis_callback = AnalysisCallback(semi; interval = 10,
+                                     extra_analysis_errors = (:conservation_error,),
+                                     extra_analysis_integrals = (waterheight_total,
+                                                                 velocity, entropy))
+callbacks = CallbackSet(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)

src/DispersiveShallowWater.jl

@@ -1,8 +1,9 @@
 module DispersiveShallowWater
 
+using BandedMatrices: BandedMatrix
 using DiffEqBase: DiffEqBase, SciMLBase, terminate!
 using Interpolations: Interpolations, linear_interpolation
-using LinearAlgebra: mul!, I, Diagonal
+using LinearAlgebra: mul!, I, Diagonal, diag
 using PolynomialBases: PolynomialBases
 using Printf: @printf, @sprintf
 using RecipesBase: RecipesBase, @recipe, @series
@@ -18,8 +19,10 @@ using SparseArrays: sparse
 using SummationByPartsOperators: AbstractDerivativeOperator,
                                  PeriodicDerivativeOperator, PeriodicUpwindOperators,
                                  UniformPeriodicCoupledOperator,
+                                 DerivativeOperator, UpwindOperators,
+                                 UniformCoupledOperator,
                                  periodic_derivative_operator,
-                                 derivative_order, integrate
+                                 derivative_order, integrate, mass_matrix
 import SummationByPartsOperators: grid, xmin, xmax
 using TimerOutputs: TimerOutputs, TimerOutput, @timeit, print_timer, reset_timer!
 @reexport using TrixiBase: TrixiBase, trixi_include
@@ -48,10 +51,11 @@ export Solver
 
 export Semidiscretization, semidiscretize, grid
 
-export boundary_condition_periodic
+export boundary_condition_periodic, boundary_condition_reflecting
 
 export initial_condition_convergence_test,
        initial_condition_manufactured, source_terms_manufactured,
+       initial_condition_manufactured_reflecting, source_terms_manufactured_reflecting,
        initial_condition_dingemans
 
 export AnalysisCallback, RelaxationCallback, SummaryCallback

src/boundary_conditions.jl

@@ -7,3 +7,15 @@ A singleton struct indicating periodic boundary conditions.
 const boundary_condition_periodic = BoundaryConditionPeriodic()
 
 Base.show(io::IO, ::BoundaryConditionPeriodic) = print(io, "boundary_condition_periodic")
+
+struct BoundaryConditionReflecting end
+
+"""
+    boundary_condition_reflecting = DispersiveShallowWater.BoundaryConditionReflecting()
+A singleton struct indicating reflecting boundary conditions.
+"""
+const boundary_condition_reflecting = BoundaryConditionReflecting()
+
+function Base.show(io::IO, ::BoundaryConditionReflecting)
+    print(io, "boundary_condition_reflecting")
+end

src/equations/bbm_bbm_1d.jl

@@ -89,24 +89,77 @@
 
     return SVector(dq1, dq2)
 end
-#
+
+"""
+    initial_condition_manufactured_reflecting(x, t, equations::BBMBBMEquations1D, mesh)
+
+A smooth manufactured solution for reflecting boundary conditions in combination
+with [`source_terms_manufactured_reflecting`](@ref).
+"""
+function initial_condition_manufactured_reflecting(x, t,
+                                                   equations::BBMBBMEquations1D,
+                                                   mesh)
+    eta = exp(2 * t) * cospi(x)
+    v = exp(t) * x * sinpi(x)
+    return SVector(eta, v)
+end
+
+"""
+    source_terms_manufactured_reflecting(q, x, t, equations::BBMBBMEquations1D, mesh)
+
+A smooth manufactured solution for reflecting boundary conditions in combination
+with [`initial_condition_manufactured_reflecting`](@ref).
+"""
+function source_terms_manufactured_reflecting(q, x, t, equations::BBMBBMEquations1D)
+    g = equations.gravity
+    D = equations.D
+    a1 = cospi(2 * x)
+    a2 = sinpi(2 * x)
+    a8 = cospi(x)
+    a9 = sinpi(x)
+    a10 = exp(t)
+    a11 = exp(2 * t)
+    dq1 = (pi^2 * D^2 * a10 * a8 / 3 + pi * D * x * a8 + D * a9 + pi * x * a11 * a1 +
+           a11 * a2 / 2 + 2 * a10 * a8) * a10
+    dq2 = (pi * D^2 * (pi * x * a9 - 2 * a8) / 6 - pi * g * a10 * a9 +
+           pi * x^2 * a10 * a2 / 2 + x * a10 * a9^2 + x * a9) * a10
+
+    return SVector(dq1, dq2)
+end
+
 function create_cache(mesh,
                       equations::BBMBBMEquations1D,
                       solver,
                       initial_condition,
                       RealT,
                       uEltype)
+    tmp1 = Array{RealT}(undef, nnodes(mesh)) # tmp1 is needed for the `RelaxationCallback`
     D = equations.D
     if solver.D1 isa PeriodicDerivativeOperator ||
-       solver.D1 isa UniformPeriodicCoupledOperator
-        invImD2 = inv(I - 1 / 6 * D^2 * Matrix(solver.D2))
-    elseif solver.D1 isa PeriodicUpwindOperators
+       solver.D1 isa UniformPeriodicCoupledOperator ||
+       solver.D1 isa PeriodicUpwindOperators
         invImD2 = inv(I - 1 / 6 * D^2 * Matrix(solver.D2))
+        return (invImD2 = invImD2, tmp1 = tmp1)
+    elseif solver.D1 isa DerivativeOperator ||
+           solver.D1 isa UniformCoupledOperator ||
+           solver.D1 isa UpwindOperators
+        D1_b = BandedMatrix(solver.D1)
+        M = mass_matrix(solver.D1)
+        Pd = BandedMatrix((-1 => fill(one(eltype(D1_b)), size(D1_b, 1) - 2),),
+                          (size(D1_b, 1), size(D1_b, 1) - 2))
+        D2d = (sparse(solver.D2) * Pd)[2:(end - 1), :]
+        # homogeneous Dirichtlet boundary conditions
+        invImD2d = inv(I - 1 / 6 * D^2 * D2d)
+        m = diag(M)
+        m[1] = 0
+        m[end] = 0
+        PdM = Diagonal(m)
+        # homogeneous Neumann boundary conditions
+        invImD2n = inv(I + 1 / 6 * D^2 * inv(M) * D1_b' * PdM * D1_b)
+        return (invImD2d = invImD2d, invImD2n = invImD2n, tmp1 = tmp1)
     else
-        @error "unknown type of first-derivative operator"
+        @error "unknown type of first-derivative operator: $(typeof(solver.D1))"
     end
-    tmp1 = Array{RealT}(undef, nnodes(mesh)) # tmp1 is needed for the `RelaxationCallback`
-    return (invImD2 = invImD2, tmp1 = tmp1)
 end
 
 # Discretization that conserves the mass (for eta and v) and the energy for periodic boundary conditions, see
@@ -137,7 +190,7 @@
         @timeit timer() "dv hyperbolic" dv[:]=-solver.D1.central *
                                               (equations.gravity * eta + 0.5 * v .^ 2)
     else
-        @error "unknown type of first-derivative operator"
+        @error "unknown type of first-derivative operator: $(typeof(solver.D1))"
     end
 
     @timeit timer() "source terms" calc_sources!(dq, q, t, source_terms, equations, solver)
@@ -148,6 +201,48 @@
     return nothing
 end
 
+# Discretization that conserves the mass (for eta) and the energy for periodic boundary conditions, see
+# - 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)
+function rhs!(du_ode, u_ode, t, mesh, equations::BBMBBMEquations1D, initial_condition,
+              ::BoundaryConditionReflecting, source_terms, solver, cache)
+    @unpack invImD2d, invImD2n = cache
+
+    q = wrap_array(u_ode, mesh, equations, solver)
+    dq = wrap_array(du_ode, mesh, equations, solver)
+
+    eta = view(q, 1, :)
+    v = view(q, 2, :)
+    deta = view(dq, 1, :)
+    dv = view(dq, 2, :)
+
+    D = equations.D
+    # energy and mass conservative semidiscretization
+    if solver.D1 isa DerivativeOperator ||
+       solver.D1 isa UniformCoupledOperator
+        @timeit timer() "deta hyperbolic" deta[:]=-solver.D1 * (D * v + eta .* v)
+        @timeit timer() "dv hyperbolic" dv[:]=-solver.D1 *
+                                              (equations.gravity * eta + 0.5 * v .^ 2)
+    elseif solver.D1 isa UpwindOperators
+        @timeit timer() "deta hyperbolic" deta[:]=-solver.D1.central * (D * v + eta .* v)
+        @timeit timer() "dv hyperbolic" dv[:]=-solver.D1.central *
+                                              (equations.gravity * eta + 0.5 * v .^ 2)
+    else
+        @error "unknown type of first-derivative operator: $(typeof(solver.D1))"
+    end
+
+    @timeit timer() "source terms" calc_sources!(dq, q, t, source_terms, equations, solver)
+
+    @timeit timer() "deta elliptic" deta[:]=invImD2n * deta
+    @timeit timer() "dv elliptic" begin
+        dv[2:(end - 1)] = invImD2d * dv[2:(end - 1)]
+        dv[1] = dv[end] = zero(eltype(dv))
+    end
+
+    return nothing
+end
+
 @inline function prim2cons(q, equations::BBMBBMEquations1D)
     eta, v = q
 

src/equations/bbm_bbm_variable_bathymetry_1d.jl

@@ -169,7 +169,7 @@
         invImDKD = inv(I - 1 / 6 * Matrix(solver.D1.minus) * K * Matrix(solver.D1.plus))
         invImD2K = inv(I - 1 / 6 * Matrix(solver.D2) * K)
     else
-        @error "unknown type of first-derivative operator"
+        @error "unknown type of first-derivative operator: $(typeof(solver.D1))"
     end
     tmp1 = Array{RealT}(undef, nnodes(mesh)) # tmp1 is needed for the `RelaxationCallback`
     return (invImDKD = invImDKD, invImD2K = invImD2K, D = D, tmp1 = tmp1)
@@ -204,7 +204,7 @@
         @timeit timer() "dv hyperbolic" dv[:]=-solver.D1.plus *
                                               (equations.gravity * eta + 0.5 * v .^ 2)
     else
-        @error "unknown type of first-derivative operator"
+        @error "unknown type of first-derivative operator: $(typeof(solver.D1))"
     end
 
     @timeit timer() "source terms" calc_sources!(dq, q, t, source_terms, equations, solver)

src/equations/svaerd_kalisch_1d.jl

@@ -190,7 +190,7 @@
         D1_central = solver.D1.central
         D1betaD1 = sparse(solver.D1.plus) * Diagonal(beta_hat) * sparse(solver.D1.minus)
     else
-        @error "unknown type of first-derivative operator"
+        @error "unknown type of first-derivative operator: $(typeof(solver.D1))"
     end
     return (hmD1betaD1 = hmD1betaD1, D1betaD1 = D1betaD1, D = D, h = h, hv = hv,
             alpha_hat = alpha_hat, beta_hat = beta_hat, gamma_hat = gamma_hat,
@@ -238,7 +238,7 @@
             yD1v = tmp1 .* (solver.D1.plus * v)
             deta[:] = solver.D1.minus * tmp1 - D1_central * hv
         else
-            @error "unknown type of first derivative operator"
+            @error "unknown type of first derivative operator: $(typeof(solver.D1))"
         end
     end
 

src/semidiscretization.jl

@@ -156,7 +156,7 @@ function integrate_quantity!(quantity, func, u_ode, semi::Semidiscretization; wr
     integrate(quantity, semi; wrap = false)
 end
 
-# modified entropy from Svärd-Kalisch equations need to take the whole vector `u` for every point in space
+# modified entropy from Svärd-Kalisch equations need to take the whole vector `q` for every point in space
 function integrate_quantity(func::Union{typeof(energy_total_modified),
                                         typeof(entropy_modified)}, u_ode,
                             semi::Semidiscretization; wrap = true)

test/test_bbm_bbm_1d.jl

@@ -54,6 +54,17 @@ EXAMPLES_DIR = joinpath(examples_dir(), "bbm_bbm_1d")
                             atol=1e-10,
                             atol_ints=1e-10) # in order to make CI pass
     end
+
+    @trixi_testset "bbm_bbm_1d_basic_reflecting" begin
+        @test_trixi_include(joinpath(EXAMPLES_DIR, "bbm_bbm_1d_basic_reflecting.jl"),
+                            tspan=(0.0, 1.0),
+                            l2=[1.44633857650787e-6 2.5981093955688672e-8],
+                            linf=[2.6974310287641856e-6 4.159028832440015e-8],
+                            cons_error=[4.2697991261385974e-11 0.5469460931577577],
+                            change_waterheight=4.2697991261385974e-11,
+                            change_velocity=0.5469460931577577,
+                            change_entropy=130.69415963528576)
+    end
 end
 
 end # module

test/test_unit.jl

@@ -93,6 +93,9 @@ using SparseArrays: sparse, SparseMatrixCSC
         boundary_conditions = boundary_condition_periodic
         @test_nowarn print(boundary_conditions)
         @test_nowarn display(boundary_conditions)
+        boundary_conditions = boundary_condition_reflecting
+        @test_nowarn print(boundary_conditions)
+        @test_nowarn display(boundary_conditions)
     end
 
     @testset "BBMBBMEquations1D" begin