Introduction
Summation-by-parts (SBP) operators are discrete derivative operators designed to enable (semi-) discrete stability proofs mimicking the energy method from the continuous level. To do so, SBP operators mimic integration-by-parts discretely. Here, we will briefly explain the basic concepts. If you want to learn more about this subject, the classical review articles of [SvärdNordström2014] and [FernándezHickenZingg2014] are good starting points. More recent references and applications of SBP operators from many classes implemented in SummationByPartsOperators.jl are given by [RanochaMitsotakisKetcheson2021].
Since SBP operators are designed to mimic integration-by-parts, they need a notion of derivatives and integrals. Here, derivatives are interpreted as linear operators D (derivative matrices) and integrals are interpreted as discrete inner products, represented by the associated mass/norm matrices M. Thus, the discrete derivative of a grid function u is D * u and the discrete inner product of two grid functions u and v is dot(u, M, v), where M = mass_matrix(D). Here, we have already introduced some basic interfaces provided by SummationByPartsOperators.jl:
- Derivative operators act as linear operators implementing
*(andmul!for more efficient in-place updates avoiding allocations). - The mass matrix associated to an SBP derivative operator can be retrieved via
mass_matrix.
Periodic domains
Periodic (central) SBP operators mimic the properties of differential operators on periodic domains. Hence, they are
- skew-symmetric if they approximate odd derivatives
- symmetric and semi-definite if they approximate even derivatives; second-derivative operators are negative semi-definite, fourth-derivative operators are positive semi-definite etc.
Classical central finite difference operators on periodic domains are periodic SBP operators. They can be constructed via periodic_derivative_operator. Similarly, Fourier collocation methods can be interpreted as periodic SBP operators, which can be constructed via fourier_derivative_operator.
julia> using SummationByPartsOperators, LinearAlgebra
julia> D = periodic_derivative_operator(derivative_order=1, accuracy_order=2,
xmin=0.0, xmax=2.0, N=20)
Periodic first-derivative operator of order 2 on a grid in [0.0, 2.0] using 20 nodes,
stencils with 1 nodes to the left, 1 nodes to the right, and coefficients of Fornberg (1998)
Calculation of Weights in Finite Difference Formulas.
SIAM Rev. 40.3, pp. 685-691.
julia> M = mass_matrix(D)
UniformScaling{Float64}
0.1*I
julia> M * Matrix(D) + Matrix(D)' * M |> norm
0.0
julia> D = fourier_derivative_operator(xmin=0.0, xmax=2.0, N=20)
Periodic 1st derivative Fourier operator {T=Float64}
on a grid in [0.0, 2.0] using 20 nodes and 11 modes
julia> M = mass_matrix(D)
UniformScaling{Float64}
0.1*I
julia> norm(M * Matrix(D) + Matrix(D)' * M) < 10 * eps(eltype(D))
trueAs you have seen above, conversion methods to other common types such as Matrix, sparse from the standard library SparseArrays, and BandedMatrix from BandedMatrices.jl are available.
Non-periodic domains
On non-periodic domains, additional boundary terms appear. Thus, the basic symmetry properties of SBP operators are the same as the ones of periodic SBP operators modulo boundary terms. Note that the correct handling of boundary terms is the basic reason of the success of SBP operators. In particular for hyperbolic problems, other boundary treatments that might appear senseful can result in catastrophic failure.
First-derivative operators
First-derivative SBP operators need to mimic
\[ \int_{x_\mathrm{min}}^{x_\mathrm{max}} u(x) \bigl( \partial_x v(x) \bigr) \mathrm{d}x + \int_{x_\mathrm{min}}^{x_\mathrm{max}} \bigl( \partial_x u(x) \bigr) v(x) \mathrm{d}x = u(x_\mathrm{max}) v(x_\mathrm{max}) - u(x_\mathrm{min}) v(x_\mathrm{min}).\]
Thus, a discrete evaluation at the boundary of the domain is necessary. For SBP operators with a grid including the boundary nodes, this can be achieved by simply picking the first/last nodal coefficient of a grid function u. If boundary nodes are not included, some interpolation is necessary in general. Nevertheless, getting a boundary value is a linear functional that is often represented in the literature using (transposed) vectors tL, tR. Then, an SBP operator has to satisfy M * D + D' * M == tR * tR' - tL * tL'. The boundary operators are represented matrix-free via derivative_left and derivative_right for zeroth-order derivatives.
julia> using SummationByPartsOperators, LinearAlgebrajulia> D = derivative_operator(MattssonNordström2004(), derivative_order = 1, accuracy_order = 2, xmin = 0//1, xmax = 1//1, N = 9)SBP first-derivative operator of order 2 on a grid in [0//1, 1//1] using 9 nodes and coefficients of Mattsson, Nordström (2004) Summation by parts operators for finite difference approximations of second derivatives. Journal of Computational Physics 199, pp. 503-540.julia> tL = zeros(eltype(D), size(D, 1)); tL[1] = 1; tL'1×9 adjoint(::Vector{Rational{Int64}}) with eltype Rational{Int64}: 1//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1julia> tR = zeros(eltype(D), size(D, 1)); tR[end] = 1; tR'1×9 adjoint(::Vector{Rational{Int64}}) with eltype Rational{Int64}: 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 1//1julia> M = mass_matrix(D)9×9 LinearAlgebra.Diagonal{Rational{Int64}, Vector{Rational{Int64}}}: 1//16 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 1//8 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 1//8 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 1//8 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 1//8 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 1//8 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 1//8 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 1//8 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 1//16julia> M * Matrix(D) + Matrix(D)' * M == tR * tR' - tL * tL'truejulia> u = randn(size(grid(D))); derivative_left(D, u, Val(0)) == u[begin]truejulia> u = randn(size(grid(D))); derivative_right(D, u, Val(0)) == u[end]true
Here, we have introduced some additional features. Firstly, exact rational coefficients are provided, based on the type of xmin and xmax (if available). Secondly, a source_of_coefficients has to be provided when constructing the SBP operator. You can list them using
using InteractiveUtils, SummationByPartsOperators
subtypes(SourceOfCoefficients)23-element Vector{Any}:
BeljaddLeFlochMishraParés2017
DienerDorbandSchnetterTiglio2007
Fornberg1998
GlaubitzNordströmÖffner2023
Holoborodko2008
LanczosLowNoise
MadayTadmor1989
Mattsson2012
Mattsson2014
Mattsson2017
⋮
MattssonNordström2004
MattssonSvärdNordström2004
MattssonSvärdShoeybi2008
SharanBradyLivescu2022
Tadmor1989
Tadmor1993
TadmorWaagan2012Convergent
TadmorWaagan2012Standard
WilliamsDuru2024Here and in the following, the order of accuracy of (finite difference) SBP operators refers to the local order of accuracy in the interior, cf. accuracy_order.
A special case of first-derivative SBP operators are polynomial derivative operators on Lobatto-Legendre nodes, implemented in legendre_derivative_operator.
Second-derivative operators
To mimic integration-by-parts of second derivatives,
\[ \int_{x_\mathrm{min}}^{x_\mathrm{max}} u(x) \bigl( \partial_x^2 v(x) \bigr) \mathrm{d}x = - \int_{x_\mathrm{min}}^{x_\mathrm{max}} \bigl( \partial_x u(x) \bigr) \bigl( \partial_x v(x) \bigr) \mathrm{d}x + u(x_\mathrm{max}) \bigl( \partial_x v(x_\mathrm{max}) \bigr) - \bigl( \partial_x u(x_\mathrm{min})) v(x_\mathrm{min}),\]
the evaluation of the first derivative at the boundaries is necessary. These linear functionals are available as derivative_left and derivative_right. In the literature, they are often called dL and dR. Then, a second-derivative SBP operator has to be of the form M * D == -A + tR * dR' - tL * dL', where A is symmetric and positive semidefinite.
julia> using SummationByPartsOperators, LinearAlgebrajulia> D = derivative_operator(MattssonNordström2004(), derivative_order=2, accuracy_order=2, xmin=0//1, xmax=1//1, N=9)SBP second-derivative operator of order 2 on a grid in [0//1, 1//1] using 9 nodes and coefficients of Mattsson, Nordström (2004) Summation by parts operators for finite difference approximations of second derivatives. Journal of Computational Physics 199, pp. 503-540.julia> M = mass_matrix(D)9×9 LinearAlgebra.Diagonal{Rational{Int64}, Vector{Rational{Int64}}}: 1//16 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 1//8 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 1//8 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 1//8 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 1//8 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 1//8 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 1//8 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 1//8 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 1//16julia> tL = derivative_left(D, Val(0)); tL'1×9 adjoint(::Vector{Bool}) with eltype Bool: 1 0 0 0 0 0 0 0 0julia> tR = derivative_right(D, Val(0)); tR'1×9 adjoint(::Vector{Bool}) with eltype Bool: 0 0 0 0 0 0 0 0 1julia> dL = derivative_left(D, Val(1)); dL'1×9 adjoint(::Vector{Rational{Int64}}) with eltype Rational{Int64}: -12//1 16//1 -4//1 0//1 0//1 0//1 0//1 0//1 0//1julia> dR = derivative_right(D, Val(1)); dR'1×9 adjoint(::Vector{Rational{Int64}}) with eltype Rational{Int64}: 0//1 0//1 0//1 0//1 0//1 0//1 4//1 -16//1 12//1julia> A = -M * Matrix(D) + tR * dR' - tL * dL'9×9 Matrix{Rational{Int64}}: 8//1 -8//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 -8//1 16//1 -8//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 -8//1 16//1 -8//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 -8//1 16//1 -8//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 -8//1 16//1 -8//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 -8//1 16//1 -8//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 -8//1 16//1 -8//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 -8//1 16//1 -8//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 -8//1 8//1julia> isposdef(A)true
Usually, there is no need to form dL, dR explicitly. Instead, you can use the matrix-free variants derivative_left and derivative_right. Some procedures imposing boundary conditions weakly require adding the transposed boundary derivatives to a grid function, which can be achieved by mul_transpose_derivative_left! and mul_transpose_derivative_right!. You can find applications of these operators in the source code of WaveEquationNonperiodicSemidiscretization.
A special case of second-derivative SBP operators are polynomial derivative operators on Lobatto-Legendre nodes, implemented in legendre_second_derivative_operator.
Upwind operators
Upwind SBP operators were introduced by Mattsson2017. They combine two derivative operators Dp (:plus) and Dm (:minus) such that M * Dp + Dm' * M == tR * tR' - tL * tL' and M * (Dp - Dm) is negative semidefinite.
julia> using SummationByPartsOperators, LinearAlgebrajulia> Dp = derivative_operator(Mattsson2017(:plus), derivative_order=1, accuracy_order=2, xmin=0//1, xmax=1//1, N=9)SBP first-derivative operator of order 2 on a grid in [0//1, 1//1] using 9 nodes and coefficients of Mattsson (2017) Diagonal-norm upwind SBP operators. Journal of Computational Physics 335, pp. 283-310. (upwind coefficients plus)julia> Matrix(Dp)9×9 Matrix{Rational{Int64}}: -24//1 40//1 -16//1 0//1 0//1 0//1 0//1 0//1 0//1 -8//5 -8//1 64//5 -16//5 0//1 0//1 0//1 0//1 0//1 0//1 0//1 -12//1 16//1 -4//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 -12//1 16//1 -4//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 -12//1 16//1 -4//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 -12//1 16//1 -4//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 -12//1 16//1 -4//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 -8//1 8//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 -8//1 8//1julia> Dm = derivative_operator(Mattsson2017(:minus), derivative_order=1, accuracy_order=2, xmin=0//1, xmax=1//1, N=9)SBP first-derivative operator of order 2 on a grid in [0//1, 1//1] using 9 nodes and coefficients of Mattsson (2017) Diagonal-norm upwind SBP operators. Journal of Computational Physics 335, pp. 283-310. (upwind coefficients minus)julia> Matrix(Dm)9×9 Matrix{Rational{Int64}}: -8//1 8//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 -8//1 8//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 4//1 -16//1 12//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 4//1 -16//1 12//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 4//1 -16//1 12//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 4//1 -16//1 12//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 4//1 -16//1 12//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 16//5 -64//5 8//1 8//5 0//1 0//1 0//1 0//1 0//1 0//1 16//1 -40//1 24//1julia> M = mass_matrix(Dp)9×9 LinearAlgebra.Diagonal{Rational{Int64}, Vector{Rational{Int64}}}: 1//32 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 5//32 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 1//8 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 1//8 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 1//8 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 1//8 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 1//8 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 5//32 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 1//32julia> M * Matrix(Dp) + Matrix(Dm)' * M9×9 Matrix{Rational{Int64}}: -1//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 1//1julia> minimum(eigvals(-M * (Matrix(Dp) - Matrix(Dm)))) # > 0 up to floating point tolerances-6.892084304488098e-16
You can also set up fully periodic upwind operators by setting the argument left_offset of periodic_derivative_operator appropriately. For example,
julia> using SummationByPartsOperators, LinearAlgebrajulia> Dp = periodic_derivative_operator(derivative_order=1, accuracy_order=2, left_offset=0, xmin=0//1, xmax=1//1, N=8)Periodic first-derivative operator of order 2 on a grid in [0//1, 1//1] using 8 nodes, stencils with 0 nodes to the left, 2 nodes to the right, and coefficients of Fornberg (1998) Calculation of Weights in Finite Difference Formulas. SIAM Rev. 40.3, pp. 685-691.julia> Matrix(Dp)8×8 Matrix{Rational{Int64}}: -12//1 16//1 -4//1 0//1 0//1 0//1 0//1 0//1 0//1 -12//1 16//1 -4//1 0//1 0//1 0//1 0//1 0//1 0//1 -12//1 16//1 -4//1 0//1 0//1 0//1 0//1 0//1 0//1 -12//1 16//1 -4//1 0//1 0//1 0//1 0//1 0//1 0//1 -12//1 16//1 -4//1 0//1 0//1 0//1 0//1 0//1 0//1 -12//1 16//1 -4//1 -4//1 0//1 0//1 0//1 0//1 0//1 -12//1 16//1 16//1 -4//1 0//1 0//1 0//1 0//1 0//1 -12//1julia> Dm = periodic_derivative_operator(derivative_order=1, accuracy_order=2, left_offset=-2, xmin=0//1, xmax=1//1, N=8)Periodic first-derivative operator of order 2 on a grid in [0//1, 1//1] using 8 nodes, stencils with 2 nodes to the left, 0 nodes to the right, and coefficients of Fornberg (1998) Calculation of Weights in Finite Difference Formulas. SIAM Rev. 40.3, pp. 685-691.julia> Matrix(Dm)8×8 Matrix{Rational{Int64}}: 12//1 0//1 0//1 0//1 0//1 0//1 4//1 -16//1 -16//1 12//1 0//1 0//1 0//1 0//1 0//1 4//1 4//1 -16//1 12//1 0//1 0//1 0//1 0//1 0//1 0//1 4//1 -16//1 12//1 0//1 0//1 0//1 0//1 0//1 0//1 4//1 -16//1 12//1 0//1 0//1 0//1 0//1 0//1 0//1 4//1 -16//1 12//1 0//1 0//1 0//1 0//1 0//1 0//1 4//1 -16//1 12//1 0//1 0//1 0//1 0//1 0//1 0//1 4//1 -16//1 12//1julia> M = mass_matrix(Dp)LinearAlgebra.UniformScaling{Rational{Int64}} 1//8*Ijulia> M * Matrix(Dp) + Matrix(Dm)' * M |> iszerotruejulia> minimum(eigvals(-M * (Matrix(Dp) - Matrix(Dm)))) # > 0 up to floating point tolerances2.1439294572106007e-16
Note that we used N=8 here, i.e., one node less than for the non-periodic example. This is necessary since the additional node at the right boundary is identified with the left boundary node for periodic operators.
To create all upwind operators for a single setup, you can use upwind_operators.
julia> using SummationByPartsOperatorsjulia> D = upwind_operators(Mattsson2017, derivative_order=1, accuracy_order=2, xmin=0, xmax=1//1, N=9)Upwind SBP first-derivative operators of order 2 on a grid in [0//1, 1//1] using 9 nodes and coefficients of Mattsson2017julia> Matrix(D.plus)9×9 Matrix{Rational{Int64}}: -24//1 40//1 -16//1 0//1 0//1 0//1 0//1 0//1 0//1 -8//5 -8//1 64//5 -16//5 0//1 0//1 0//1 0//1 0//1 0//1 0//1 -12//1 16//1 -4//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 -12//1 16//1 -4//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 -12//1 16//1 -4//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 -12//1 16//1 -4//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 -12//1 16//1 -4//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 -8//1 8//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 -8//1 8//1julia> Matrix(D.minus)9×9 Matrix{Rational{Int64}}: -8//1 8//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 -8//1 8//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 4//1 -16//1 12//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 4//1 -16//1 12//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 4//1 -16//1 12//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 4//1 -16//1 12//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 4//1 -16//1 12//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 16//5 -64//5 8//1 8//5 0//1 0//1 0//1 0//1 0//1 0//1 16//1 -40//1 24//1
This also works with periodic upwind operators.
julia> using SummationByPartsOperatorsjulia> D = upwind_operators(periodic_derivative_operator, accuracy_order = 2, xmin = 0, xmax = 1//1, N = 10)Upwind SBP first-derivative operators of order 2 on a grid in [0//1, 9//10] using 10 nodes and coefficients of Fornberg1998julia> Matrix(D.plus)10×10 Matrix{Rational{Int64}}: -15//1 20//1 -5//1 0//1 0//1 … 0//1 0//1 0//1 0//1 0//1 -15//1 20//1 -5//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 -15//1 20//1 -5//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 -15//1 20//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 -15//1 -5//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 … 20//1 -5//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 -15//1 20//1 -5//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 -15//1 20//1 -5//1 -5//1 0//1 0//1 0//1 0//1 0//1 0//1 -15//1 20//1 20//1 -5//1 0//1 0//1 0//1 0//1 0//1 0//1 -15//1julia> Matrix(D.minus)10×10 Matrix{Rational{Int64}}: 15//1 0//1 0//1 0//1 0//1 … 0//1 0//1 5//1 -20//1 -20//1 15//1 0//1 0//1 0//1 0//1 0//1 0//1 5//1 5//1 -20//1 15//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 5//1 -20//1 15//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 5//1 -20//1 15//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 5//1 -20//1 … 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 5//1 15//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 -20//1 15//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 5//1 -20//1 15//1 0//1 0//1 0//1 0//1 0//1 0//1 0//1 5//1 -20//1 15//1
You can also couple upwind operators continuously across elements using couple_continuously to obtain global upwind operators, see below and Theorem 2.4 of [RanochaMitsotakisKetcheson2021].
Similarly, you can couple classical and upwind operators discontinuously across elements using couple_discontinuously to obtain global upwind operators, see below and Theorem 2.2 of [RanochaMitsotakisKetcheson2021].
Continuous Galerkin methods
SBP operators can be coupled to obtain (nodal) continuous Galerkin (CG) methods. If the underlying SBP operators are LegendreDerivativeOperators, these are CG spectral element methods (CGSEM). However, a continuous coupling of arbitrary SBP operators is supported.
julia> using SummationByPartsOperators, LinearAlgebrajulia> D = couple_continuously( legendre_derivative_operator(xmin=-1.0, xmax=1.0, N=3), UniformMesh1D(xmin=0.0, xmax=1.0, Nx=3))First derivative operator {T=Float64} on 3 Lobatto Legendre nodes in [-1.0, 1.0] coupled continuously on UniformMesh1D{Float64} with 3 cells in (0.0, 1.0)julia> Matrix(D)7×7 Matrix{Float64}: -9.0 12.0 -3.0 0.0 0.0 0.0 0.0 -3.0 0.0 3.0 0.0 0.0 0.0 0.0 1.5 -6.0 0.0 6.0 -1.5 0.0 0.0 0.0 0.0 -3.0 0.0 3.0 0.0 0.0 0.0 0.0 1.5 -6.0 0.0 6.0 -1.5 0.0 0.0 0.0 0.0 -3.0 0.0 3.0 0.0 0.0 0.0 0.0 3.0 -12.0 9.0julia> mass_matrix(D)7×7 LinearAlgebra.Diagonal{Float64, Vector{Float64}}: 0.0555556 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.222222 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.111111 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.222222 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.111111 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.222222 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.0555556
Discontinuous Galerkin methods
SBP operators can also be coupled as in discontinuous Galerkin (DG) methods. Using a central numerical flux results in central SBP operators; upwind fluxes yield upwind SBP operators. If LegendreDerivativeOperators are used, the discontinuous coupling yields DG spectral element methods (DGSEM).
julia> using SummationByPartsOperators, LinearAlgebrajulia> D = couple_discontinuously( legendre_derivative_operator(xmin=-1.0, xmax=1.0, N=3), UniformPeriodicMesh1D(xmin=0.0, xmax=1.0, Nx=3), Val(:central))First derivative operator {T=Float64} on 3 Lobatto Legendre nodes in [-1.0, 1.0] coupled discontinuously (upwind: Val{:central}()) on UniformPeriodicMesh1D{Float64} with 3 cells in (0.0, 1.0)julia> M = mass_matrix(D);julia> M * Matrix(D) + Matrix(D)' * M |> iszerotrue
Right now, only uniform meshes UniformMesh1D and UniformPeriodicMesh1D are implemented.
You can also specify a different coupling than Val(:central) to obtain upwind operators.
julia> using SummationByPartsOperators, LinearAlgebrajulia> Dp = couple_discontinuously( legendre_derivative_operator(xmin=-1.0, xmax=1.0, N=3), UniformPeriodicMesh1D(xmin=0.0, xmax=1.0, Nx=3), Val(:plus))First derivative operator {T=Float64} on 3 Lobatto Legendre nodes in [-1.0, 1.0] coupled discontinuously (upwind: Val{:plus}()) on UniformPeriodicMesh1D{Float64} with 3 cells in (0.0, 1.0)julia> Matrix(Dp)9×9 Matrix{Float64}: -9.0 12.0 -3.0 0.0 0.0 0.0 0.0 0.0 0.0 -3.0 0.0 3.0 0.0 0.0 0.0 0.0 0.0 0.0 3.0 -12.0 -9.0 18.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -9.0 12.0 -3.0 0.0 0.0 0.0 0.0 0.0 0.0 -3.0 0.0 3.0 0.0 0.0 0.0 0.0 0.0 0.0 3.0 -12.0 -9.0 18.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -9.0 12.0 -3.0 0.0 0.0 0.0 0.0 0.0 0.0 -3.0 0.0 3.0 18.0 0.0 0.0 0.0 0.0 0.0 3.0 -12.0 -9.0julia> Dm = couple_discontinuously( legendre_derivative_operator(xmin=-1.0, xmax=1.0, N=3), UniformPeriodicMesh1D(xmin=0.0, xmax=1.0, Nx=3), Val(:minus))First derivative operator {T=Float64} on 3 Lobatto Legendre nodes in [-1.0, 1.0] coupled discontinuously (upwind: Val{:minus}()) on UniformPeriodicMesh1D{Float64} with 3 cells in (0.0, 1.0)julia> Matrix(Dm)9×9 Matrix{Float64}: 9.0 12.0 -3.0 0.0 0.0 0.0 0.0 0.0 -18.0 -3.0 0.0 3.0 0.0 0.0 0.0 0.0 0.0 0.0 3.0 -12.0 9.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -18.0 9.0 12.0 -3.0 0.0 0.0 0.0 0.0 0.0 0.0 -3.0 0.0 3.0 0.0 0.0 0.0 0.0 0.0 0.0 3.0 -12.0 9.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -18.0 9.0 12.0 -3.0 0.0 0.0 0.0 0.0 0.0 0.0 -3.0 0.0 3.0 0.0 0.0 0.0 0.0 0.0 0.0 3.0 -12.0 9.0
Basic interfaces and additional features
To actually compute and plot the discrete grid functions, a few additional ingredients are necessary.
- The discrete coefficients of a function on the
gridof an SBP operator can usually be computed asx = grid(D); u = u_function.(x), at least for nodal bases. In general,compute_coefficients(or the in-place versioncompute_coefficients!) can also be used for this task. - To get a grid and discrete values suitable for plotting, you can use
evaluate_coefficients(or the in-place versionevaluate_coefficients!). The plot nodes returned fromevaluate_coefficientscan be different from the nodes of thegridassociated to an SBP operator. - To implement boundary procedures, the weights of the mass matrix at the boundary are often needed. These can be obtained without forming
M = mass_matrix(D)explicitly vialeft_boundary_weightandright_boundary_weight. - Instead of forming a mass matrix explicitly, discrete integrals can be evaluated efficiently using
integrateand vectors can be scaled by the mass matrix or its inverse usingscale_by_mass_matrix!andscale_by_inverse_mass_matrix!, respectively. - Dissipation operators based on the same discrete inner product as SBP derivative operators can be obtained via
dissipation_operator.
Next steps
If you are familiar with SBP operators in general, this introduction might already be enough for you to apply SummationByPartsOperators.jl to your problems. Otherwise, you might want to have a look at the references, the tutorials coming next, or some ready-to-use semidiscretizations of the following partial differential equations (PDEs). These are shipped with this package and you are encouraged to look at their source code to learn more about it.
- Linear scalar advection with variable coefficient:
VariableLinearAdvectionPeriodicSemidiscretization,VariableLinearAdvectionNonperiodicSemidiscretization - Burgers' equation (inviscid):
BurgersPeriodicSemidiscretization,BurgersNonperiodicSemidiscretization - Scalar conservation law with cubic flux:
CubicPeriodicSemidiscretization,CubicNonperiodicSemidiscretization - A scalar conservation law with quartic, non-convex flux:
QuarticNonconvexPeriodicSemidiscretization - The second-order wave equation:
WaveEquationNonperiodicSemidiscretization
Some additional examples are included as Jupyter notebooks in the directory notebooks. Even more examples and research articles making use of SummationByPartsOperators.jl are listed in the section Applications. If you want to know even more, you can have a look at the test.
References
- SvärdNordström2014Svärd, Nordström (2014). Review of summation-by-parts schemes for initial–boundary-value problems. DOI: 10.1016/j.jcp.2014.02.031
- FernándezHickenZingg2014Fernández, Hicken, Zingg (2014). Review of summation-by-parts operators with simultaneous approximation terms for the numerical solution of partial differential equations. DOI: 10.1016/j.compfluid.2014.02.016
- RanochaMitsotakisKetcheson2021Ranocha, Mitsotakis, Ketcheson (2021). A Broad Class of Conservative Numerical Methods for Dispersive Wave Equations. DOI: 10.4208/cicp.OA-2020-0119