Add support for source terms

examples/bbm_bbm_1d/bbm_bbm_1d_manufactured.jl

@@ -0,0 +1,40 @@
+using OrdinaryDiffEq
+using DispersiveShallowWater
+
+###############################################################################
+# Semidiscretization of the BBM-BBM equations
+
+equations = BBMBBMEquations1D(gravity_constant = 9.81, D = 4.0)
+
+initial_condition = initial_condition_manufactured
+source_terms = source_terms_manufactured
+boundary_conditions = boundary_condition_periodic
+
+# create homogeneous mesh
+coordinates_min = 0.0
+coordinates_max = 1.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,
+                          source_terms = source_terms)
+
+###############################################################################
+# Create `ODEProblem` and run the simulation
+tspan = (0.0, 1.0)
+ode = semidiscretize(semi, tspan)
+analysis_callback = AnalysisCallback(semi; interval = 10,
+                                     extra_analysis_errors = (:conservation_error,),
+                                     extra_analysis_integrals = (waterheight_total,
+                                                                 velocity, entropy))
+callbacks = CallbackSet(analysis_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_bbm_variable_bathymetry_1d/bbm_bbm_variable_bathymetry_1d_manufactured.jl

@@ -0,0 +1,40 @@
+using OrdinaryDiffEq
+using DispersiveShallowWater
+
+###############################################################################
+# Semidiscretization of the BBM-BBM equations
+
+equations = BBMBBMVariableEquations1D(gravity_constant = 9.81)
+
+initial_condition = initial_condition_manufactured
+source_terms = source_terms_manufactured
+boundary_conditions = boundary_condition_periodic
+
+# create homogeneous mesh
+coordinates_min = 0.0
+coordinates_max = 1.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,
+                          source_terms = source_terms)
+
+###############################################################################
+# Create `ODEProblem` and run the simulation
+tspan = (0.0, 1.0)
+ode = semidiscretize(semi, tspan)
+analysis_callback = AnalysisCallback(semi; interval = 10,
+                                     extra_analysis_errors = (:conservation_error,),
+                                     extra_analysis_integrals = (waterheight_total,
+                                                                 velocity, entropy))
+callbacks = CallbackSet(analysis_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

@@ -23,8 +23,8 @@ import SummationByPartsOperators: grid, xmin, xmax
 
 include("boundary_conditions.jl")
 include("mesh.jl")
-include("solver.jl")
 include("equations/equations.jl")
+include("solver.jl")
 include("semidiscretization.jl")
 include("callbacks_step/callbacks_step.jl")
 include("visualization.jl")
@@ -47,7 +47,9 @@ export Semidiscretization, semidiscretize, grid
 
 export boundary_condition_periodic
 
-export initial_condition_convergence_test, initial_condition_dingemans
+export initial_condition_convergence_test,
+       initial_condition_manufactured, source_terms_manufactured,
+       initial_condition_dingemans
 
 export AnalysisCallback, RelaxationCallback
 export tstops, errors, integrals

src/equations/bbm_bbm_1d.jl

@@ -53,6 +53,41 @@ function initial_condition_convergence_test(x, t, equations::BBMBBMEquations1D,
     return SVector(eta, v)
 end
 
+"""
+    initial_condition_manufactured(x, t, equations::BBMBBMEquations1D, mesh)
+
+A smooth manufactured solution in combination with `source_terms_manufactured`.
+"""
+function initial_condition_manufactured(x, t,
+                                        equations::BBMBBMEquations1D,
+                                        mesh)
+    eta = exp(t) * cospi(2 * (x - 2 * t))
+    v = exp(t / 2) * sinpi(2 * (x - t / 2))
+    return SVector(eta, v)
+end
+
+"""
+    source_terms_manufactured(q, x, t, equations::BBMBBMEquations1D, mesh)
+
+A smooth manufactured solution in combination with `initial_condition_manufactured`.
+"""
+function source_terms_manufactured(q, x, t, equations::BBMBBMEquations1D)
+    g = equations.gravity
+    D = equations.D
+    a3 = cospi(t - 2 * x)
+    a4 = sinpi(t - 2 * x)
+    a5 = sinpi(2 * t - 4 * x)
+    a6 = sinpi(4 * t - 2 * x)
+    a7 = cospi(4 * t - 2 * x)
+    dq1 = -2 * pi^2 * D^2 * (4 * pi * a6 - a7) * exp(t) / 3 + 2 * pi * D * exp(t / 2) * a3 -
+          2 * pi * exp(3 * t / 2) * a4 * a6 + 2 * pi * exp(3 * t / 2) * a3 * a7 -
+          4 * pi * exp(t) * a6 + exp(t) * a7
+    dq2 = -pi^2 * D^2 * (a4 + 2 * pi * a3) * exp(t / 2) / 3 + 2 * pi * g * exp(t) * a6 -
+          exp(t / 2) * a4 / 2 - pi * exp(t / 2) * a3 - pi * exp(t) * a5
+
+    return SVector(dq1, dq2)
+end
+#
 function create_cache(mesh,
                       equations::BBMBBMEquations1D,
                       solver,
@@ -61,24 +96,23 @@ function create_cache(mesh,
                       uEltype)
     if solver.D1 isa PeriodicDerivativeOperator ||
        solver.D1 isa UniformPeriodicCoupledOperator
-        invImD2_D = (I - 1 / 6 * equations.D^2 * sparse(solver.D2)) \ Matrix(solver.D1)
+        invImD2 = inv(I - 1 / 6 * equations.D^2 * Matrix(solver.D2))
     elseif solver.D1 isa PeriodicUpwindOperators
-        invImD2_D = (I - 1 / 6 * equations.D^2 * sparse(solver.D2)) \
-                    Matrix(solver.D1.central)
+        invImD2 = inv(I - 1 / 6 * equations.D^2 * Matrix(solver.D2))
     else
         @error "unknown type of first-derivative operator"
     end
-    tmp1 = Array{RealT}(undef, nnodes(mesh))
-    return (invImD2_D = invImD2_D, tmp1 = tmp1)
+    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
 # - 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,
-              ::BoundaryConditionPeriodic, solver, cache)
-    @unpack invImD2_D, tmp1 = cache
+              ::BoundaryConditionPeriodic, source_terms, solver, cache)
+    @unpack invImD2 = cache
 
     q = wrap_array(u_ode, mesh, equations, solver)
     dq = wrap_array(du_ode, mesh, equations, solver)
@@ -89,11 +123,21 @@ function rhs!(du_ode, u_ode, t, mesh, equations::BBMBBMEquations1D, initial_cond
     dv = view(dq, 2, :)
 
     # energy and mass conservative semidiscretization
-    @. tmp1 = -(equations.D * v + eta * v)
-    mul!(deta, invImD2_D, tmp1)
+    if solver.D1 isa PeriodicDerivativeOperator ||
+       solver.D1 isa UniformPeriodicCoupledOperator
+        deta[:] = -solver.D1 * (equations.D * v + eta .* v)
+        dv[:] = -solver.D1 * (equations.gravity * eta + 0.5 * v .^ 2)
+    elseif solver.D1 isa PeriodicUpwindOperators
+        deta[:] = -solver.D1.central * (equations.D * v + eta .* v)
+        dv[:] = -solver.D1.central * (equations.gravity * eta + 0.5 * v .^ 2)
+    else
+        @error "unknown type of first-derivative operator"
+    end
+
+    calc_sources!(dq, q, t, source_terms, equations, solver)
 
-    @. tmp1 = -(equations.gravity * eta + 0.5 * v^2)
-    mul!(dv, invImD2_D, tmp1)
+    deta[:] = invImD2 * deta
+    dv[:] = invImD2 * dv
 
     return nothing
 end

src/equations/bbm_bbm_variable_bathymetry_1d.jl

@@ -41,8 +41,7 @@ For details see Example 5 in Section 3 from (here adapted for dimensional equati
   Exact Traveling-Wave Solutions to Bidirectional Wave Equations
   [DOI: 10.1023/A:1026667903256](https://doi.org/10.1023/A:1026667903256)
 """
-function initial_condition_convergence_test(x,
-                                            t,
+function initial_condition_convergence_test(x, t,
                                             equations::BBMBBMVariableEquations1D,
                                             mesh)
     g = equations.gravity
@@ -58,6 +57,49 @@ function initial_condition_convergence_test(x,
     return SVector(eta, v, D)
 end
 
+"""
+    initial_condition_manufactured(x, t, equations::BBMBBMVariableEquations1D, mesh)
+
+A smooth manufactured solution in combination with `source_terms_manufactured`.
+"""
+function initial_condition_manufactured(x, t,
+                                        equations::BBMBBMVariableEquations1D,
+                                        mesh)
+    eta = exp(t) * cospi(2 * (x - 2 * t))
+    v = exp(t / 2) * sinpi(2 * (x - t / 2))
+    D = 5.0 + 2.0 * cospi(2 * x)
+    return SVector(eta, v, D)
+end
+
+"""
+    source_terms_manufactured(q, x, t, equations::BBMBBMVariableEquations1D, mesh)
+
+A smooth manufactured solution in combination with `initial_condition_manufactured`.
+"""
+function source_terms_manufactured(q, x, t, equations::BBMBBMVariableEquations1D)
+    g = equations.gravity
+    a1 = cospi(2 * x)
+    a2 = sinpi(2 * x)
+    a3 = cospi(t - 2 * x)
+    a4 = sinpi(t - 2 * x)
+    a5 = sinpi(2 * t - 4 * x)
+    a6 = sinpi(4 * t - 2 * x)
+    a7 = cospi(4 * t - 2 * x)
+    dq1 = -2 * pi^2 * (4 * pi * a6 - a7) * (2 * a1 + 5)^2 * exp(t) / 3 +
+          8 * pi^2 * (a6 + 4 * pi * a7) * (2 * a1 + 5) * exp(t) * a2 / 3 +
+          2 * pi * (2 * a1 + 5) * exp(t / 2) * a3 - 2 * pi * exp(3 * t / 2) * a4 * a6 +
+          2 * pi * exp(3 * t / 2) * a3 * a7 + 4 * pi * exp(t / 2) * a2 * a4 -
+          4 * pi * exp(t) * a6 + exp(t) * a7
+    dq2 = 2 * pi * g * exp(t) * a6 -
+          pi^2 *
+          (8 * (2 * pi * a4 - a3) * (2 * a1 + 5) * a2 +
+           (a4 + 2 * pi * a3) * (2 * a1 + 5)^2 +
+           4 * (a4 + 2 * pi * a3) * (16 * sinpi(x)^4 - 26 * sinpi(x)^2 + 7)) * exp(t / 2) /
+          3 - exp(t / 2) * a4 / 2 - pi * exp(t / 2) * a3 - pi * exp(t) * a5
+
+    return SVector(dq1, dq2, 0.0)
+end
+
 """
     initial_condition_dingemans(x, t, equations::BBMBBMVariableEquations1D, mesh)
 
@@ -100,7 +142,7 @@ end
 
 function create_cache(mesh,
                       equations::BBMBBMVariableEquations1D,
-                      solver::Solver,
+                      solver,
                       initial_condition,
                       RealT,
                       uEltype)
@@ -113,18 +155,16 @@ function create_cache(mesh,
     K = Diagonal(D .^ 2)
     if solver.D1 isa PeriodicDerivativeOperator ||
        solver.D1 isa UniformPeriodicCoupledOperator
-        invImDKD_D = (I - 1 / 6 * sparse(solver.D1) * K * sparse(solver.D1)) \
-                     Matrix(solver.D1)
-        invImD2K_D = (I - 1 / 6 * sparse(solver.D2) * K) \ Matrix(solver.D1)
+        invImDKD = inv(I - 1 / 6 * Matrix(solver.D1) * K * Matrix(solver.D1))
+        invImD2K = inv(I - 1 / 6 * Matrix(solver.D2) * K)
     elseif solver.D1 isa PeriodicUpwindOperators
-        invImDKD_D = (I - 1 / 6 * sparse(solver.D1.minus) * K * sparse(solver.D1.plus)) \
-                     Matrix(solver.D1.minus)
-        invImD2K_D = (I - 1 / 6 * sparse(solver.D2) * K) \ Matrix(solver.D1.plus)
+        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"
     end
-    tmp1 = Array{RealT}(undef, nnodes(mesh))
-    return (invImDKD_D = invImDKD_D, invImD2K_D = invImD2K_D, tmp1 = tmp1)
+    tmp1 = Array{RealT}(undef, nnodes(mesh)) # tmp1 is needed for the `RelaxationCallback`
+    return (invImDKD = invImDKD, invImD2K = invImD2K, tmp1 = tmp1)
 end
 
 # Discretization that conserves the mass (for eta and v) and the energy for periodic boundary conditions, see
@@ -133,9 +173,9 @@ end
 #   [DOI: 10.4208/cicp.OA-2020-0119](https://doi.org/10.4208/cicp.OA-2020-0119)
 # Here, adapted for spatially varying bathymetry.
 function rhs!(du_ode, u_ode, t, mesh, equations::BBMBBMVariableEquations1D,
-              initial_condition,
-              ::BoundaryConditionPeriodic, solver, cache)
-    @unpack invImDKD_D, invImD2K_D, tmp1 = cache
+              initial_condition, ::BoundaryConditionPeriodic, source_terms,
+              solver, cache)
+    @unpack invImDKD, invImD2K = cache
 
     q = wrap_array(u_ode, mesh, equations, solver)
     dq = wrap_array(du_ode, mesh, equations, solver)
@@ -148,11 +188,21 @@ function rhs!(du_ode, u_ode, t, mesh, equations::BBMBBMVariableEquations1D,
     dD = view(dq, 3, :)
     fill!(dD, zero(eltype(dD)))
 
-    @. tmp1 = -(D * v + eta * v)
-    mul!(deta, invImDKD_D, tmp1)
+    if solver.D1 isa PeriodicDerivativeOperator ||
+       solver.D1 isa UniformPeriodicCoupledOperator
+        deta[:] = -solver.D1 * (D .* v + eta .* v)
+        dv[:] = -solver.D1 * (equations.gravity * eta + 0.5 * v .^ 2)
+    elseif solver.D1 isa PeriodicUpwindOperators
+        deta[:] = -solver.D1.minus * (D .* v + eta .* v)
+        dv[:] = -solver.D1.plus * (equations.gravity * eta + 0.5 * v .^ 2)
+    else
+        @error "unknown type of first-derivative operator"
+    end
+
+    calc_sources!(dq, q, t, source_terms, equations, solver)
 
-    @. tmp1 = -(equations.gravity * eta + 0.5 * v^2)
-    mul!(dv, invImD2K_D, tmp1)
+    deta[:] = invImDKD * deta
+    dv[:] = invImD2K * dv
 
     return nothing
 end

src/equations/svaerd_kalisch_1d.jl

@@ -86,7 +86,7 @@ end
 
 function create_cache(mesh,
                       equations::SvaerdKalischEquations1D,
-                      solver::Solver,
+                      solver,
                       initial_condition,
                       RealT,
                       uEltype)
@@ -123,8 +123,8 @@ end
 
 # Discretization that conserves the mass (for eta and v) and is energy-bounded for periodic boundary conditions
 function rhs!(du_ode, u_ode, t, mesh, equations::SvaerdKalischEquations1D,
-              initial_condition,
-              ::BoundaryConditionPeriodic, solver, cache)
+              initial_condition, ::BoundaryConditionPeriodic, source_terms,
+              solver, cache)
     @unpack hmD1betaD1, D1betaD1, d, h, hv, alpha_hat, beta_hat, gamma_hat, tmp1, tmp2, D1_central = cache
     q = wrap_array(u_ode, mesh, equations, solver)
     dq = wrap_array(du_ode, mesh, equations, solver)
@@ -176,6 +176,7 @@ function rhs!(du_ode, u_ode, t, mesh, equations::SvaerdKalischEquations1D,
     #              0.5 * D1_central * (gamma_hat .* (solver.D2 * v)) -
     #              0.5 * solver.D2 * (gamma_hat .* D1v))
     dv[:] = hmD1betaD1 \ tmp2
+    calc_sources!(dq, q, t, source_terms, equations, solver)
 
     return nothing
 end