add support for 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 = -35.0
-coordinates_max = 35.0
+coordinates_min = 0.0
+coordinates_max = 1.0
 N = 512
 mesh = Mesh1D(coordinates_min, coordinates_max, N)
 
@@ -22,11 +22,12 @@ 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)
+                          boundary_conditions = boundary_conditions,
+                          source_terms = source_terms)
 
 ###############################################################################
 # Create `ODEProblem` and run the simulation
-tspan = (0.0, 30.0)
+tspan = (0.0, 1.0)
 ode = semidiscretize(semi, tspan)
 analysis_callback = AnalysisCallback(semi; interval = 10,
                                      extra_analysis_errors = (:conservation_error,),

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,8 @@ export Semidiscretization, semidiscretize, grid
 
 export boundary_condition_periodic
 
-export initial_condition_convergence_test, initial_condition_dingemans
+export initial_condition_convergence_test, source_terms_convergence_test,
+       initial_condition_dingemans
 
 export AnalysisCallback, RelaxationCallback
 export tstops, errors, integrals

src/equations/bbm_bbm_1d.jl

@@ -77,7 +77,7 @@ end
 #   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)
+              ::BoundaryConditionPeriodic, source_terms, solver, cache)
     @unpack invImD2_D, tmp1 = cache
 
     q = wrap_array(u_ode, mesh, equations, solver)
@@ -94,6 +94,7 @@ function rhs!(du_ode, u_ode, t, mesh, equations::BBMBBMEquations1D, initial_cond
 
     @. tmp1 = -(equations.gravity * eta + 0.5 * v^2)
     mul!(dv, invImD2_D, tmp1)
+    calc_sources!(dq, q, t, source_terms, equations, solver)
 
     return nothing
 end

src/equations/bbm_bbm_variable_bathymetry_1d.jl

@@ -33,31 +33,44 @@ varnames(::typeof(prim2cons), ::BBMBBMVariableEquations1D) = ("h", "hv", "b")
 """
     initial_condition_convergence_test(x, t, equations::BBMBBMVariableEquations1D, mesh)
 
-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)
+A smooth manufactured solution in combination with `source_terms_convergence_test`.
 """
-function initial_condition_convergence_test(x,
-                                            t,
+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
+    eta = exp(t) * cospi(2 * (x - 2 * t))
+    v = exp(t / 2) * sinpi(2 * (x - t / 2))
+    D = cospi(2 * x)
     return SVector(eta, v, D)
 end
 
+"""
+    source_terms_convergence_test(q, x, t, equations::BBMBBMVariableEquations1D, mesh)
+
+A smooth manufactured solution in combination with `initial_condition_convergence_test`.
+"""
+function source_terms_convergence_test(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 / 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 -
+          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
+    return SVector(dq1, dq2, 0.0)
+end
+
 """
     initial_condition_dingemans(x, t, equations::BBMBBMVariableEquations1D, mesh)
 
@@ -100,7 +113,7 @@ end
 
 function create_cache(mesh,
                       equations::BBMBBMVariableEquations1D,
-                      solver::Solver,
+                      solver,
                       initial_condition,
                       RealT,
                       uEltype)
@@ -133,8 +146,8 @@ 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)
+              initial_condition, ::BoundaryConditionPeriodic, source_terms,
+              solver, cache)
     @unpack invImDKD_D, invImD2K_D, tmp1 = cache
 
     q = wrap_array(u_ode, mesh, equations, solver)
@@ -153,6 +166,7 @@ function rhs!(du_ode, u_ode, t, mesh, equations::BBMBBMVariableEquations1D,
 
     @. tmp1 = -(equations.gravity * eta + 0.5 * v^2)
     mul!(dv, invImD2K_D, tmp1)
+    calc_sources!(dq, q, t, source_terms, equations, solver)
 
     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

src/semidiscretization.jl

@@ -6,7 +6,7 @@ A `struct` containing everything needed to describe a spatial semidiscretization
 of an equation.
 """
 struct Semidiscretization{Mesh, Equations, InitialCondition, BoundaryConditions,
-                          Solver, Cache}
+                          SourceTerms, Solver, Cache}
     mesh::Mesh
     equations::Equations
 
@@ -15,30 +15,34 @@ struct Semidiscretization{Mesh, Equations, InitialCondition, BoundaryConditions,
     initial_condition::InitialCondition
 
     boundary_conditions::BoundaryConditions
+    source_terms::SourceTerms
     solver::Solver
     cache::Cache
 
     function Semidiscretization{Mesh, Equations, InitialCondition, BoundaryConditions,
-                                Solver,
+                                SourceTerms, Solver,
                                 Cache}(mesh::Mesh, equations::Equations,
                                        initial_condition::InitialCondition,
                                        boundary_conditions::BoundaryConditions,
+                                       source_terms::SourceTerms,
                                        solver::Solver,
                                        cache::Cache) where {Mesh, Equations,
                                                             InitialCondition,
-                                                            BoundaryConditions, Solver,
-                                                            Cache}
+                                                            BoundaryConditions, SourceTerms,
+                                                            Solver, Cache}
         @assert ndims(mesh) == ndims(equations)
         @assert xmin(mesh) == xmin(solver.D1)
         @assert xmax(mesh) == xmax(solver.D1)
         @assert nnodes(mesh) == length(grid(solver))
 
-        new(mesh, equations, initial_condition, boundary_conditions, solver, cache)
+        new(mesh, equations, initial_condition, boundary_conditions, source_terms, solver,
+            cache)
     end
 end
 
 """
     Semidiscretization(mesh, equations, initial_condition, solver;
+                       source_terms=nothing,
                        boundary_conditions=boundary_condition_periodic,
                        RealT=real(solver),
                        uEltype=RealT,
@@ -47,6 +51,7 @@ end
 Construct a semidiscretization of a PDE.
 """
 function Semidiscretization(mesh, equations, initial_condition, solver;
+                            source_terms = nothing,
                             boundary_conditions = boundary_condition_periodic,
                             # `RealT` is used as real type for node locations etc.
                             # while `uEltype` is used as element type of solutions etc.
@@ -57,12 +62,10 @@ function Semidiscretization(mesh, equations, initial_condition, solver;
              initial_cache...)
 
     Semidiscretization{typeof(mesh), typeof(equations), typeof(initial_condition),
-                       typeof(boundary_conditions), typeof(solver), typeof(cache)}(mesh,
-                                                                                   equations,
-                                                                                   initial_condition,
-                                                                                   boundary_conditions,
-                                                                                   solver,
-                                                                                   cache)
+                       typeof(boundary_conditions), typeof(source_terms),
+                       typeof(solver), typeof(cache)}(mesh, equations, initial_condition,
+                                                      boundary_conditions, source_terms,
+                                                      solver, cache)
 end
 
 function Base.show(io::IO, semi::Semidiscretization)
@@ -73,6 +76,7 @@ function Base.show(io::IO, semi::Semidiscretization)
     print(io, ", ", semi.equations)
     print(io, ", ", semi.initial_condition)
     print(io, ", ", semi.boundary_conditions)
+    print(io, ", ", semi.source_terms)
     print(io, ", ", semi.solver)
     print(io, ", cache(")
     for (idx, key) in enumerate(keys(semi.cache))
@@ -93,7 +97,8 @@ function Base.show(io::IO, ::MIME"text/plain", semi::Semidiscretization)
         println(io, "    mesh: ", semi.mesh)
         println(io, "    equations: ", get_name(semi.equations))
         println(io, "    initial condition: ", semi.initial_condition)
-        print(io, "    boundary condition: ", semi.boundary_conditions)
+        println(io, "    boundary condition: ", semi.boundary_conditions)
+        print(io, "    source terms: ", semi.source_terms)
     end
 end
 
@@ -196,10 +201,11 @@ function wrap_array(u_ode, semi::Semidiscretization)
 end
 
 function rhs!(du_ode, u_ode, semi::Semidiscretization, t)
-    @unpack mesh, equations, initial_condition, boundary_conditions, solver, cache = semi
+    @unpack mesh, equations, initial_condition, boundary_conditions, solver, source_terms, cache = semi
 
-    rhs!(du_ode, u_ode, t, mesh, equations, initial_condition, boundary_conditions, solver,
-         cache)
+    rhs!(du_ode, u_ode, t, mesh, equations, initial_condition, boundary_conditions,
+         source_terms,
+         solver, cache)
 
     return nothing
 end

src/solver.jl

@@ -114,3 +114,20 @@ function calc_error_norms(u_ode, t, initial_condition, mesh::Mesh1D, equations,
     end
     return l2_error, linf_error
 end
+
+function calc_sources!(dq, q, t, source_terms::Nothing,
+                       equations::AbstractEquations{1}, solver::Solver)
+    return nothing
+end
+
+function calc_sources!(dq, q, t, source_terms,
+                       equations::AbstractEquations{1}, solver::Solver)
+    x = grid(solver)
+    for i in eachnode(solver)
+        local_source = source_terms(view(q, :, i), x[i], t, equations)
+        for v in eachvariable(equations)
+            dq[v, i] += local_source[v]
+        end
+    end
+    return nothing
+end