Basic interface
Here, we discuss the basic interface of SummationByPartsOperators.jl. We assume you are already familiar with the concept of SBP operators in general and the introduction describing how to construct specific operators.
Applying SBP operators
All SBP operators implement the general interface of matrix vector multiplication in Julia. The most simple version is to just use *, e.g.,
julia> using SummationByPartsOperatorsjulia> D = derivative_operator(MattssonNordström2004(), derivative_order = 1, accuracy_order = 2, xmin = 0.0, xmax = 1.0, N = 9)SBP first-derivative operator of order 2 on a grid in [0.0, 1.0] 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> x = grid(D)0.0:0.125:1.0julia> u = @. sin(pi * x)9-element Vector{Float64}: 0.0 0.3826834323650898 0.7071067811865475 0.9238795325112867 1.0 0.9238795325112867 0.7071067811865476 0.3826834323650899 1.2246467991473532e-16julia> D * u9-element Vector{Float64}: 3.0614674589207183 2.82842712474619 2.1647844005847876 1.1715728752538102 0.0 -1.1715728752538097 -2.1647844005847876 -2.82842712474619 -3.0614674589207183julia> @allocated D * u160
As you can see above, calling D * u allocates a new vector for the result. If you want to apply an SBP operator multiple times and need good performance, you should consider using pre-allocating the output and using in-place update instead. This strategy is also described in the performance tips in the Julia manual. Julia provides the function mul! for this purpose.
julia> using LinearAlgebra, InteractiveUtilsjulia> @doc mul!mul!(Y, A, B) -> Y Calculates the matrix-matrix or matrix-vector product AB and stores the result in Y, overwriting the existing value of Y. Note that Y must not be aliased with either A or B. Examples ≡≡≡≡≡≡≡≡≡≡ julia> A=[1.0 2.0; 3.0 4.0]; B=[1.0 1.0; 1.0 1.0]; Y = similar(B); mul!(Y, A, B); julia> Y 2×2 Matrix{Float64}: 3.0 3.0 7.0 7.0 Implementation ≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡ For custom matrix and vector types, it is recommended to implement 5-argument mul! rather than implementing 3-argument mul! directly if possible. mul!(C, A, B, α, β) -> C Combined inplace matrix-matrix or matrix-vector multiply-add A B α + C β. The result is stored in C by overwriting it. Note that C must not be aliased with either A or B. │ Julia 1.3 │ │ Five-argument mul! requires at least Julia 1.3. Examples ≡≡≡≡≡≡≡≡≡≡ julia> A=[1.0 2.0; 3.0 4.0]; B=[1.0 1.0; 1.0 1.0]; C=[1.0 2.0; 3.0 4.0]; julia> mul!(C, A, B, 100.0, 10.0) === C true julia> C 2×2 Matrix{Float64}: 310.0 320.0 730.0 740.0 mul!(du, D::DerivativeOperator, u, α=true, β=false) Efficient in-place version of du = α * D * u + β * du. Note that du must not be aliased with u. Do mul! with the off-diagonal elements of a matrix.
To improve the performance, you can pre-allocate an output vector and call the non-allocating function mul!.
julia> using SummationByPartsOperatorsjulia> D = derivative_operator(MattssonNordström2004(), derivative_order = 1, accuracy_order = 2, xmin = 0.0, xmax = 1.0, N = 9)SBP first-derivative operator of order 2 on a grid in [0.0, 1.0] 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> x = grid(D)0.0:0.125:1.0julia> u = @. sin(pi * x)9-element Vector{Float64}: 0.0 0.3826834323650898 0.7071067811865475 0.9238795325112867 1.0 0.9238795325112867 0.7071067811865476 0.3826834323650899 1.2246467991473532e-16julia> du = similar(u); mul!(du, D, u)julia> du ≈ D * utruejulia> @allocated mul!(du, D, u)0
All operators provided by SummationByPartsOperators.jl implement this 3-argument version of mul!. Most operators also implement the 5-argument version of mul! that can be used to scale the output and add it to some multiple of the result vector.
julia> using SummationByPartsOperatorsjulia> D = derivative_operator(MattssonNordström2004(), derivative_order = 1, accuracy_order = 2, xmin = 0.0, xmax = 1.0, N = 9)SBP first-derivative operator of order 2 on a grid in [0.0, 1.0] 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> x = grid(D); u = @. sin(pi * x); du = similar(u); mul!(du, D, u);julia> mul!(du, D, u, 2) # equivalent to du .= 2 * D * ujulia> du ≈ 2 * D * utruejulia> @allocated mul!(du, D, u, 2)0julia> du_background = rand(length(du)); du .= du_background9-element Vector{Float64}: 0.7517279027973112 0.6015659445719002 0.6744713659019441 0.44092591837629325 0.4060932979564582 0.2514381702985309 0.13676972151180622 0.022916175595679267 0.3085399004917837julia> mul!(du, D, u, 2, 3) # equivalent to du .= 2 * D * u + 3 * dujulia> du ≈ 2 * D * u + 3 * du_backgroundtruejulia> @allocated mul!(du, D, u, 2, 3)0
Integration and the mass/norm matrix
SBP operators come with a mass matrix yielding a quadrature rule. In SummationByPartsOperators.jl, all operators typically have diagonal mass/norm matrices. You can access them via mass_matrix, e.g.,
julia> using SummationByPartsOperatorsjulia> D = derivative_operator(MattssonNordström2004(), derivative_order = 1, accuracy_order = 2, xmin = 0.0, xmax = 1.0, N = 9)SBP first-derivative operator of order 2 on a grid in [0.0, 1.0] 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> mass_matrix(D)9×9 LinearAlgebra.Diagonal{Float64, Vector{Float64}}: 0.0625 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.125 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.125 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.125 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.125 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.125 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.125 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.125 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.0625julia> D = periodic_derivative_operator(derivative_order = 1, accuracy_order = 2, xmin = 0.0, xmax = 1.0, N = 8)Periodic first-derivative operator of order 2 on a grid in [0.0, 1.0] using 8 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> mass_matrix(D)LinearAlgebra.UniformScaling{Float64} 0.125*I
If you want to use the quadrature associated with a mass matrix, you do not need to form it explicitly. Instead, it is recommended to use the function integrate, e.g.,
julia> using SummationByPartsOperatorsjulia> D = derivative_operator(MattssonNordström2004(), derivative_order = 1, accuracy_order = 2, xmin = 0.0, xmax = 1.0, N = 9)SBP first-derivative operator of order 2 on a grid in [0.0, 1.0] 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{Float64, Vector{Float64}}: 0.0625 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.125 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.125 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.125 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.125 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.125 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.125 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.125 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.0625julia> x = grid(D)0.0:0.125:1.0julia> u = x.^29-element Vector{Float64}: 0.0 0.015625 0.0625 0.140625 0.25 0.390625 0.5625 0.765625 1.0julia> integrate(u, D)0.3359375julia> integrate(u, D) ≈ sum(M * u)truejulia> integrate(u, D) ≈ integrate(x -> x^2, x, D)true
For example, you can proceed as follows to compute the error of the SBP operator when computing a derivative as follows.
julia> using SummationByPartsOperatorsjulia> D = derivative_operator(MattssonNordström2004(), derivative_order = 1, accuracy_order = 2, xmin = 0.0, xmax = 1.0, N = 9)SBP first-derivative operator of order 2 on a grid in [0.0, 1.0] 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{Float64, Vector{Float64}}: 0.0625 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.125 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.125 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.125 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.125 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.125 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.125 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.125 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.0625julia> x = grid(D)0.0:0.125:1.0julia> difference = D * x.^3 - 3 * x.^29-element Vector{Float64}: 0.015625 0.015625 0.015625 0.015625 0.015625 0.015625 0.015625 0.015625 -0.359375julia> error_l2 = sqrt(integrate(abs2, difference, D))0.09110862335695782
Multi-dimensional cases or multiple variables
If you want to work with multiple space dimensions, you can still use the 1D operators provided by SummationByPartsOperators.jl if you apply them in a tensor product fashion along each space dimension.
julia> using SummationByPartsOperatorsjulia> D = derivative_operator(MattssonNordström2004(), derivative_order = 1, accuracy_order = 4, xmin = 0.0, xmax = 1.0, N = 9)SBP first-derivative operator of order 4 on a grid in [0.0, 1.0] 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> x = y = grid(D)0.0:0.125:1.0julia> u = x .* y'.^2 # u(x, y) = x y^29×9 Matrix{Float64}: 0.0 0.0 0.0 0.0 … 0.0 0.0 0.0 0.0 0.00195312 0.0078125 0.0175781 0.0703125 0.0957031 0.125 0.0 0.00390625 0.015625 0.0351562 0.140625 0.191406 0.25 0.0 0.00585938 0.0234375 0.0527344 0.210938 0.287109 0.375 0.0 0.0078125 0.03125 0.0703125 0.28125 0.382812 0.5 0.0 0.00976562 0.0390625 0.0878906 … 0.351562 0.478516 0.625 0.0 0.0117188 0.046875 0.105469 0.421875 0.574219 0.75 0.0 0.0136719 0.0546875 0.123047 0.492188 0.669922 0.875 0.0 0.015625 0.0625 0.140625 0.5625 0.765625 1.0julia> let du_dx = zero(u) for j in axes(u, 2) mul!(view(du_dx, :, j), D, view(u, :, j)) end # The derivative of x*y^2 with respect to x is just y^2. # Thus, the result is constant in each column and varies # in the rows. du_dx end9×9 Matrix{Float64}: 0.0 0.015625 0.0625 0.140625 0.25 0.390625 0.5625 0.765625 1.0 0.0 0.015625 0.0625 0.140625 0.25 0.390625 0.5625 0.765625 1.0 0.0 0.015625 0.0625 0.140625 0.25 0.390625 0.5625 0.765625 1.0 0.0 0.015625 0.0625 0.140625 0.25 0.390625 0.5625 0.765625 1.0 0.0 0.015625 0.0625 0.140625 0.25 0.390625 0.5625 0.765625 1.0 0.0 0.015625 0.0625 0.140625 0.25 0.390625 0.5625 0.765625 1.0 0.0 0.015625 0.0625 0.140625 0.25 0.390625 0.5625 0.765625 1.0 0.0 0.015625 0.0625 0.140625 0.25 0.390625 0.5625 0.765625 1.0 0.0 0.015625 0.0625 0.140625 0.25 0.390625 0.5625 0.765625 1.0julia> let du_dy = zero(u) for i in axes(u, 1) mul!(view(du_dy, i, :), D, view(u, i, :)) end # The derivative of x*y^2 with respect to y is 2*x*y. du_dy end9×9 Matrix{Float64}: 0.0 0.0 0.0 0.0 … 0.0 0.0 0.0 0.0 2.1684e-19 0.03125 0.0625 0.09375 0.15625 0.1875 0.21875 0.25 4.33681e-19 0.0625 0.125 0.1875 0.3125 0.375 0.4375 0.5 4.11997e-18 0.09375 0.1875 0.28125 0.46875 0.5625 0.65625 0.75 8.67362e-19 0.125 0.25 0.375 0.625 0.75 0.875 1.0 -1.6263e-17 0.15625 0.3125 0.46875 … 0.78125 0.9375 1.09375 1.25 8.23994e-18 0.1875 0.375 0.5625 0.9375 1.125 1.3125 1.5 1.88651e-17 0.21875 0.4375 0.65625 1.09375 1.3125 1.53125 1.75 1.73472e-18 0.25 0.5 0.75 1.25 1.5 1.75 2.0julia> 2 .* x .* y'9×9 Matrix{Float64}: 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.03125 0.0625 0.09375 0.125 0.15625 0.1875 0.21875 0.25 0.0 0.0625 0.125 0.1875 0.25 0.3125 0.375 0.4375 0.5 0.0 0.09375 0.1875 0.28125 0.375 0.46875 0.5625 0.65625 0.75 0.0 0.125 0.25 0.375 0.5 0.625 0.75 0.875 1.0 0.0 0.15625 0.3125 0.46875 0.625 0.78125 0.9375 1.09375 1.25 0.0 0.1875 0.375 0.5625 0.75 0.9375 1.125 1.3125 1.5 0.0 0.21875 0.4375 0.65625 0.875 1.09375 1.3125 1.53125 1.75 0.0 0.25 0.5 0.75 1.0 1.25 1.5 1.75 2.0
Here, we have used views to interpret parts of the memory of the multi-dimensional arrays as one-diemnsional vectors that can be used together with the operators of SummationByPartsOperators.jl. You can use the same trick if you collect values of multiple variables in a multi-dimensional array.