WIP: new tutorial on basic interfaces (#254)

* WIP

* WIP

* also describe 5-argument mul!

* integration

* finish first draft

* bump version

Project.toml

@@ -1,7 +1,7 @@
 name = "SummationByPartsOperators"
 uuid = "9f78cca6-572e-554e-b819-917d2f1cf240"
 author = ["Hendrik Ranocha"]
-version = "0.5.53"
+version = "0.5.54"
 
 [deps]
 ArgCheck = "dce04be8-c92d-5529-be00-80e4d2c0e197"

docs/make.jl

@@ -53,6 +53,7 @@ makedocs(
     "Home" => "index.md",
     "Introduction" => "introduction.md",
     "Tutorials" => [
+      "tutorials/basic_interface.md",
       "tutorials/constant_linear_advection.md",
       "tutorials/advection_diffusion.md",
       "tutorials/variable_linear_advection.md",

docs/src/index.md

@@ -46,18 +46,20 @@ julia> using SummationByPartsOperators
 
 julia> using Plots: plot, plot!
 
-julia> D = periodic_derivative_operator(derivative_order=1, accuracy_order=2,
-                                        xmin=0.0, xmax=2.0, N=20)
+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> x = grid(D); u = sinpi.(x);
 
-julia> plot(x, D * u, label="numerical")
+julia> plot(x, D * u, label = "numerical")
 
-julia> plot!(x, π .* cospi.(x), label="analytical")
+julia> plot!(x, π .* cospi.(x), label = "analytical")
 ```
 You should see a plot like the following.
 
@@ -70,8 +72,9 @@ julia> using SummationByPartsOperators
 
 julia> using Plots: plot, plot!
 
-julia> D = derivative_operator(MattssonNordström2004(), derivative_order=1, accuracy_order=2,
-                               xmin=0.0, xmax=1.0, N=21)
+julia> D = derivative_operator(MattssonNordström2004(),
+                               derivative_order = 1, accuracy_order = 2,
+                               xmin = 0.0, xmax = 1.0, N = 21)
 SBP first-derivative operator of order 2 on a grid in [0.0, 1.0] using 21 nodes
 and coefficients of Mattsson, Nordström (2004)
   Summation by parts operators for finite difference approximations of second
@@ -80,9 +83,9 @@ and coefficients of Mattsson, Nordström (2004)
 
 julia> x = grid(D); u = exp.(x);
 
-julia> plot(x, D * u, label="numerical")
+julia> plot(x, D * u, label = "numerical")
 
-julia> plot!(x, exp.(x), label="analytical")
+julia> plot!(x, exp.(x), label = "analytical")
 ```
 You should see a plot like the following.
 

docs/src/introduction.md

@@ -1,4 +1,4 @@
-# Introduction
+# [Introduction](@id intro-introduction)
 
 Summation-by-parts (SBP) operators are discrete derivative operators designed to
 enable (semi-) discrete stability proofs mimicking the energy method from the
@@ -105,8 +105,9 @@ derivatives.
 ```@repl
 using SummationByPartsOperators, LinearAlgebra
 
-D = derivative_operator(MattssonNordström2004(), derivative_order=1, accuracy_order=2,
-                        xmin=0//1, xmax=1//1, N=9)
+D = derivative_operator(MattssonNordström2004(),
+                        derivative_order = 1, accuracy_order = 2,
+                        xmin = 0//1, xmax = 1//1, N = 9)
 tL = zeros(eltype(D), size(D, 1)); tL[1] = 1; tL'
 tR = zeros(eltype(D), size(D, 1)); tR[end] = 1; tR'
 M = mass_matrix(D)

docs/src/tutorials/basic_interface.md

@@ -0,0 +1,163 @@
+# Basic interface
+
+Here, we discuss the basic interface of
+[SummationByPartsOperators.jl](https://github.com/ranocha/SummationByPartsOperators.jl).
+We assume you are already familiar with the concept of SBP operators
+in general and the [introduction](@ref intro-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.,
+
+```@repl
+using SummationByPartsOperators
+
+D = derivative_operator(MattssonNordström2004(),
+                        derivative_order = 1, accuracy_order = 2,
+                        xmin = 0.0, xmax = 1.0, N = 9)
+
+x = grid(D)
+
+u = @. sin(pi * x)
+
+D * u
+
+@allocated D * u
+```
+
+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 an in-place update instead.
+Julia provides the function `mul!` for this purpose.
+
+```@repl
+using LinearAlgebra, InteractiveUtils
+
+@doc mul!
+```
+
+To improve the performance, you can pre-allocate an output vector
+and call the non-allocating function `mul!`.
+
+```@repl
+using SummationByPartsOperators
+
+D = derivative_operator(MattssonNordström2004(),
+                        derivative_order = 1, accuracy_order = 2,
+                        xmin = 0.0, xmax = 1.0, N = 9)
+
+x = grid(D)
+
+u = @. sin(pi * x)
+
+du = similar(u); mul!(du, D, u)
+
+du ≈ D * u
+
+@allocated mul!(du, D, u)
+```
+
+All operators provided by
+[SummationByPartsOperators.jl](https://github.com/ranocha/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.
+
+```@repl
+using SummationByPartsOperators
+
+D = derivative_operator(MattssonNordström2004(),
+                        derivative_order = 1, accuracy_order = 2,
+                        xmin = 0.0, xmax = 1.0, N = 9)
+
+x = grid(D); u = @. sin(pi * x); du = similar(u); mul!(du, D, u);
+
+mul!(du, D, u, 2) # equivalent to du .= 2 * D * u
+
+du ≈ 2 * D * u
+
+@allocated mul!(du, D, u, 2)
+
+du_background = rand(length(du)); du .= du_background
+
+mul!(du, D, u, 2, 3) # equivalent to du .= 2 * D * u + 3 * du
+
+du ≈ 2 * D * u + 3 * du_background
+
+@allocated mul!(du, D, u, 2, 3)
+```
+
+
+## Integration and the mass/norm matrix
+
+SBP operators come with a mass matrix yielding a quadrature rule. In
+[SummationByPartsOperators.jl](https://github.com/ranocha/SummationByPartsOperators.jl),
+all operators typically have diagonal mass/norm matrices.
+You can access them via [`mass_matrix`](@ref), e.g.,
+
+```@repl
+using SummationByPartsOperators
+
+D = derivative_operator(MattssonNordström2004(),
+                        derivative_order = 1, accuracy_order = 2,
+                        xmin = 0.0, xmax = 1.0, N = 9)
+
+mass_matrix(D)
+
+D = periodic_derivative_operator(derivative_order = 1,
+                                 accuracy_order = 2,
+                                 xmin = 0.0, xmax = 1.0,
+                                 N = 8)
+
+mass_matrix(D)
+```
+
+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`](@ref), e.g.,
+
+```@repl
+using SummationByPartsOperators
+
+D = derivative_operator(MattssonNordström2004(),
+                        derivative_order = 1, accuracy_order = 2,
+                        xmin = 0.0, xmax = 1.0, N = 9)
+
+M = mass_matrix(D)
+
+x = grid(D)
+
+u = x.^2
+
+integrate(u, D)
+
+integrate(u, D) ≈ sum(M * u)
+
+integrate(u, D) ≈ integrate(x -> x^2, x, D)
+```
+
+For example, you can proceed as follows to compute the error of the
+SBP operator when computing a derivative as follows.
+
+
+
+```@repl
+using SummationByPartsOperators
+
+D = derivative_operator(MattssonNordström2004(),
+                        derivative_order = 1, accuracy_order = 2,
+                        xmin = 0.0, xmax = 1.0, N = 9)
+
+M = mass_matrix(D)
+
+x = grid(D)
+
+difference = D * x.^3 - 3 * x.^2
+
+error_l2 = sqrt(integrate(abs2, difference, D))
+```