fix calculation of source terms

examples/bbm_bbm_variable_bathymetry_1d/bbm_bbm_variable_bathymetry_1d_basic.jl

@@ -6,13 +6,13 @@ using DispersiveShallowWater
 
 equations = BBMBBMVariableEquations1D(gravity_constant = 9.81)
 
+# initial_condition_convergence_test needs periodic boundary conditions
 initial_condition = initial_condition_convergence_test
-source_terms = source_terms_convergence_test
 boundary_conditions = boundary_condition_periodic
 
 # create homogeneous mesh
-coordinates_min = 0.0
-coordinates_max = 1.0
+coordinates_min = -35.0
+coordinates_max = 35.0
 N = 512
 mesh = Mesh1D(coordinates_min, coordinates_max, N)
 
@@ -22,12 +22,11 @@ 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)
+                          boundary_conditions = boundary_conditions)
 
 ###############################################################################
 # Create `ODEProblem` and run the simulation
-tspan = (0.0, 1.0)
+tspan = (0.0, 30.0)
 ode = semidiscretize(semi, tspan)
 analysis_callback = AnalysisCallback(semi; interval = 10,
                                      extra_analysis_errors = (:conservation_error,),

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

@@ -47,7 +47,8 @@ export Semidiscretization, semidiscretize, grid
 
 export boundary_condition_periodic
 
-export initial_condition_convergence_test, source_terms_convergence_test,
+export initial_condition_convergence_test,
+       initial_condition_manufactured, source_terms_manufactured,
        initial_condition_dingemans
 
 export AnalysisCallback, RelaxationCallback

src/equations/bbm_bbm_variable_bathymetry_1d.jl

@@ -33,23 +33,50 @@ varnames(::typeof(prim2cons), ::BBMBBMVariableEquations1D) = ("h", "hv", "b")
 """
     initial_condition_convergence_test(x, t, equations::BBMBBMVariableEquations1D, mesh)
 
-A smooth manufactured solution in combination with `source_terms_convergence_test`.
+A travelling-wave solution used for convergence tests in a periodic domain.
+The bathymetry is constant.
+
+For details see Example 5 in Section 3 from (here adapted for dimensional equations):
+- Min Chen (1997)
+  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,
                                             equations::BBMBBMVariableEquations1D,
                                             mesh)
+    g = equations.gravity
+    D = 2.0 # constant bathymetry in this case
+    c = 5 / 2 * sqrt(D * g)
+    rho = 18 / 5
+    x_t = mod(x - c * t - xmin(mesh), xmax(mesh) - xmin(mesh)) + xmin(mesh)
+
+    theta = 0.5 * sqrt(rho) * x_t / D
+    eta = -D + c^2 * rho^2 / (81 * g) +
+          5 * c^2 * rho^2 / (108 * g) * (2 * sech(theta)^2 - 3 * sech(theta)^4)
+    v = c * (1 - 5 * rho / 18) + 5 * c * rho / 6 * sech(theta)^2
+    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 = cospi(2 * x)
+    D = 5.0 + 2.0 * cospi(2 * x)
     return SVector(eta, v, D)
 end
 
 """
-    source_terms_convergence_test(q, x, t, equations::BBMBBMVariableEquations1D, mesh)
+    source_terms_manufactured(q, x, t, equations::BBMBBMVariableEquations1D, mesh)
 
-A smooth manufactured solution in combination with `initial_condition_convergence_test`.
+A smooth manufactured solution in combination with `initial_condition_manufactured`.
 """
-function source_terms_convergence_test(q, x, t, equations::BBMBBMVariableEquations1D)
+function source_terms_manufactured(q, x, t, equations::BBMBBMVariableEquations1D)
     g = equations.gravity
     a1 = cospi(2 * x)
     a2 = sinpi(2 * x)
@@ -58,16 +85,18 @@ function source_terms_convergence_test(q, x, t, equations::BBMBBMVariableEquatio
     a5 = sinpi(2 * t - 4 * x)
     a6 = sinpi(4 * t - 2 * x)
     a7 = cospi(4 * t - 2 * x)
-    dq1 = -2 / 3 * pi^2 * (4 * pi * a6 - a7) * exp(t) * a1^2 +
-          4 / 3 * pi^2 * (a6 + 4 * pi * a7) * exp(t) * a2 * a1 -
-          2 * pi * exp(3 * t / 2) * a4 * a6 + 2 * pi * exp(3 * t / 2) * a3 * a7 +
-          2 * pi * exp(t / 2) * a2 * a4 + 2 * pi * exp(t / 2) * a1 * a3 -
+    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 / 3 * (2 * pi * a4 - a3) * a2 * a1 -
-           4 / 3 * (a2^2 - a1^2) * (a4 + 2 * pi * a3) + 2 / 3 * (a4 + 2 * pi * a3) * a1^2) *
-          exp(t / 2) / 2 - exp(t / 2) * a4 / 2 - pi * exp(t / 2) * a3 - pi * exp(t) * a5
+          (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
 
@@ -126,18 +155,17 @@ 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)
+    tmp2 = similar(tmp1)
+    return (invImDKD = invImDKD, invImD2K = invImD2K, tmp1 = tmp1, tmp2 = tmp2)
 end
 
 # Discretization that conserves the mass (for eta and v) and the energy for periodic boundary conditions, see
@@ -148,7 +176,7 @@ end
 function rhs!(du_ode, u_ode, t, mesh, equations::BBMBBMVariableEquations1D,
               initial_condition, ::BoundaryConditionPeriodic, source_terms,
               solver, cache)
-    @unpack invImDKD_D, invImD2K_D, tmp1 = cache
+    @unpack invImDKD, invImD2K, tmp1 = cache
 
     q = wrap_array(u_ode, mesh, equations, solver)
     dq = wrap_array(du_ode, mesh, equations, solver)
@@ -161,13 +189,22 @@ 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
 
-    @. tmp1 = -(equations.gravity * eta + 0.5 * v^2)
-    mul!(dv, invImD2K_D, tmp1)
     calc_sources!(dq, q, t, source_terms, equations, solver)
 
+    deta[:] = invImDKD * deta
+    dv[:] = invImD2K * dv
+
     return nothing
 end