From ca1e56d961902d8545be51a52303ac70d211f532 Mon Sep 17 00:00:00 2001
From: Hendrik Ranocha <ranocha@users.noreply.github.com>
Date: Fri, 8 Apr 2022 13:07:49 +0200
Subject: [PATCH] update docs (#137)

* show versions in benchmarks

* fix typo

* remove doctest filter

* remove doctest filter

* more repl blocks

* mention periodic upwind operators in the introduction

* bump version
---
 Project.toml             |   2 +-
 docs/src/benchmarks.md   |  46 +++++++-
 docs/src/introduction.md | 246 ++++++++++++---------------------------
 3 files changed, 119 insertions(+), 175 deletions(-)

diff --git a/Project.toml b/Project.toml
index b2e8c85c..7fbe67d9 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "SummationByPartsOperators"
 uuid = "9f78cca6-572e-554e-b819-917d2f1cf240"
 author = ["Hendrik Ranocha"]
-version = "0.5.15"
+version = "0.5.16"
 
 [deps]
 ArgCheck = "dce04be8-c92d-5529-be00-80e4d2c0e197"
diff --git a/docs/src/benchmarks.md b/docs/src/benchmarks.md
index 82180e8a..72eadee6 100644
--- a/docs/src/benchmarks.md
+++ b/docs/src/benchmarks.md
@@ -57,6 +57,17 @@ DiffEqOperators.jl.
 doit(D_DEO, "D_DEO:", du, u)
 ```
 
+These results were obtained using the following versions.
+```@example first-derivative-periodic
+using InteractiveUtils
+versioninfo()
+
+using Pkg
+Pkg.status(["SummationByPartsOperators", "DiffEqOperators"],
+           mode=PKGMODE_MANIFEST)
+nothing # hide
+```
+
 
 #### Bounded domains
 
@@ -106,6 +117,17 @@ SummationByPartsOperators.jl.
 doit(D_banded, "D_banded:", du, u)
 ```
 
+These results were obtained using the following versions.
+```@example first-derivative-bounded
+using InteractiveUtils
+versioninfo()
+
+using Pkg
+Pkg.status(["SummationByPartsOperators", "BandedMatrices"],
+           mode=PKGMODE_MANIFEST)
+nothing # hide
+```
+
 
 ## Dissipation operators
 
@@ -160,6 +182,17 @@ to represent the derivative operator. This is obviously a bad idea but 🤷
 doit(Di_full, "Di_full:", du, u)
 ```
 
+These results were obtained using the following versions.
+```@example dissipation
+using InteractiveUtils
+versioninfo()
+
+using Pkg
+Pkg.status(["SummationByPartsOperators", "BandedMatrices"],
+           mode=PKGMODE_MANIFEST)
+nothing # hide
+```
+
 
 ## Structure-of-Arrays (SoA) and Array-of-Structures (AoS)
 
@@ -175,7 +208,7 @@ To demonstrate this, let us set up some benchmark code.
 using BenchmarkTools
 using StaticArrays, StructArrays
 using LinearAlgebra, SparseArrays
-using SummationByPartsOperators, BandedMatrices
+using SummationByPartsOperators
 
 BLAS.set_num_threads(1) # make sure that BLAS is serial to be fair
 
@@ -285,3 +318,14 @@ show(stdout, MIME"text/plain"(), @benchmark mul!($du_soa, $D_sparse, $u_soa))
 println("\nD_full")
 show(stdout, MIME"text/plain"(), @benchmark mul!($du_soa, $D_full, $u_soa))
 ```
+
+These results were obtained using the following versions.
+```@example soa-aos
+using InteractiveUtils
+versioninfo()
+
+using Pkg
+Pkg.status(["SummationByPartsOperators", "StaticArrays", "StructArrays"],
+           mode=PKGMODE_MANIFEST)
+nothing # hide
+```
diff --git a/docs/src/introduction.md b/docs/src/introduction.md
index 6c771226..a24cd003 100644
--- a/docs/src/introduction.md
+++ b/docs/src/introduction.md
@@ -102,45 +102,18 @@ The boundary operators are represented matrix-free via
 [`derivative_left`](@ref) and [`derivative_right`](@ref) for zeroth-order
 derivatives.
 
-```jldoctest; filter = r"((┌.*[\n\r]+)(│ `LoopVectorization.check_args` on your inputs failed; running fallback `@inbounds @fastmath` loop instead\.*[\n\r]+)(│.*[\n\r]+)(└.*[\n\r]+))"
-julia> using SummationByPartsOperators, LinearAlgebra
-
-julia> 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//1
-
-julia> 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//1
-
-julia> M = mass_matrix(D)
-9×9 Diagonal{Rational{Int64}, Vector{Rational{Int64}}}:
- 1//16   ⋅     ⋅     ⋅     ⋅     ⋅     ⋅     ⋅     ⋅
-  ⋅     1//8   ⋅     ⋅     ⋅     ⋅     ⋅     ⋅     ⋅
-  ⋅      ⋅    1//8   ⋅     ⋅     ⋅     ⋅     ⋅     ⋅
-  ⋅      ⋅     ⋅    1//8   ⋅     ⋅     ⋅     ⋅     ⋅
-  ⋅      ⋅     ⋅     ⋅    1//8   ⋅     ⋅     ⋅     ⋅
-  ⋅      ⋅     ⋅     ⋅     ⋅    1//8   ⋅     ⋅     ⋅
-  ⋅      ⋅     ⋅     ⋅     ⋅     ⋅    1//8   ⋅     ⋅
-  ⋅      ⋅     ⋅     ⋅     ⋅     ⋅     ⋅    1//8   ⋅
-  ⋅      ⋅     ⋅     ⋅     ⋅     ⋅     ⋅     ⋅    1//16
-
-julia> M * Matrix(D) + Matrix(D)' * M == tR * tR' - tL * tL'
-true
-
-julia> u = randn(size(grid(D))); derivative_left(D, u, Val(0)) == u[begin]
-true
-
-julia> u = randn(size(grid(D))); derivative_right(D, u, Val(0)) == u[end]
-true
+```@repl
+using SummationByPartsOperators, LinearAlgebra
+
+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)
+
+M * Matrix(D) + Matrix(D)' * M == tR * tR' - tL * tL'
+u = randn(size(grid(D))); derivative_left(D, u, Val(0)) == u[begin]
+u = randn(size(grid(D))); derivative_right(D, u, Val(0)) == u[end]
 ```
 
 Here, we have introduced some additional features. Firstly, exact rational
@@ -174,59 +147,20 @@ linear functionals are available as [`derivative_left`](@ref) and
 `M * D == -A + tR * dR' - tL * dL'`, where `A` is symmetric and positive
 semidefinite.
 
-```jldoctest
-julia> using SummationByPartsOperators, LinearAlgebra
+```@repl
+using SummationByPartsOperators, LinearAlgebra
 
-julia> 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.
+D = derivative_operator(MattssonNordström2004(), derivative_order=2, accuracy_order=2,
+                        xmin=0//1, xmax=1//1, N=9)
 
-julia> M = mass_matrix(D)
-9×9 Diagonal{Rational{Int64}, Vector{Rational{Int64}}}:
- 1//16   ⋅     ⋅     ⋅     ⋅     ⋅     ⋅     ⋅     ⋅
-  ⋅     1//8   ⋅     ⋅     ⋅     ⋅     ⋅     ⋅     ⋅
-  ⋅      ⋅    1//8   ⋅     ⋅     ⋅     ⋅     ⋅     ⋅
-  ⋅      ⋅     ⋅    1//8   ⋅     ⋅     ⋅     ⋅     ⋅
-  ⋅      ⋅     ⋅     ⋅    1//8   ⋅     ⋅     ⋅     ⋅
-  ⋅      ⋅     ⋅     ⋅     ⋅    1//8   ⋅     ⋅     ⋅
-  ⋅      ⋅     ⋅     ⋅     ⋅     ⋅    1//8   ⋅     ⋅
-  ⋅      ⋅     ⋅     ⋅     ⋅     ⋅     ⋅    1//8   ⋅
-  ⋅      ⋅     ⋅     ⋅     ⋅     ⋅     ⋅     ⋅    1//16
-
-julia> tL = derivative_left(D, Val(0)); tL'
-1×9 adjoint(::Vector{Bool}) with eltype Bool:
- 1  0  0  0  0  0  0  0  0
-
-julia> tR = derivative_right(D, Val(0)); tR'
-1×9 adjoint(::Vector{Bool}) with eltype Bool:
- 0  0  0  0  0  0  0  0  1
-
-julia> 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//1
-
-julia> 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//1
-
-julia> 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//1
-
-julia> isposdef(A)
-true
+M = mass_matrix(D)
+tL = derivative_left(D, Val(0)); tL'
+tR = derivative_right(D, Val(0)); tR'
+dL = derivative_left(D, Val(1)); dL'
+dR = derivative_right(D, Val(1)); dR'
+
+A = -M * Matrix(D) + tR * dR' - tL * dL'
+isposdef(A)
 ```
 Usually, there is no need to form `dL, dR` explicitly. Instead, you can use the
 matrix-free variants [`derivative_left`](@ref) and [`derivative_right`](@ref).
@@ -247,53 +181,43 @@ 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.
 
-```jldoctest
-julia> using SummationByPartsOperators, LinearAlgebra
+```@repl
+using SummationByPartsOperators, LinearAlgebra
 
-julia> 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> 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> M = mass_matrix(Dp)
-9×9 Diagonal{Rational{Int64}, Vector{Rational{Int64}}}:
- 1//32   ⋅      ⋅     ⋅     ⋅     ⋅     ⋅     ⋅      ⋅
-  ⋅     5//32   ⋅     ⋅     ⋅     ⋅     ⋅     ⋅      ⋅
-  ⋅      ⋅     1//8   ⋅     ⋅     ⋅     ⋅     ⋅      ⋅
-  ⋅      ⋅      ⋅    1//8   ⋅     ⋅     ⋅     ⋅      ⋅
-  ⋅      ⋅      ⋅     ⋅    1//8   ⋅     ⋅     ⋅      ⋅
-  ⋅      ⋅      ⋅     ⋅     ⋅    1//8   ⋅     ⋅      ⋅
-  ⋅      ⋅      ⋅     ⋅     ⋅     ⋅    1//8   ⋅      ⋅
-  ⋅      ⋅      ⋅     ⋅     ⋅     ⋅     ⋅    5//32   ⋅
-  ⋅      ⋅      ⋅     ⋅     ⋅     ⋅     ⋅     ⋅     1//32
-
-julia> M * Matrix(Dp) + Matrix(Dm)' * M
-9×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//1
-
-julia> minimum(eigvals(-M * (Matrix(Dp) - Matrix(Dm)))) > -100 * eps() # tolerance for zero eigenvalues
-true
+Dp = derivative_operator(Mattsson2017(:plus), derivative_order=1, accuracy_order=2,
+                         xmin=0//1, xmax=1//1, N=9)
+Matrix(Dp)
+Dm = derivative_operator(Mattsson2017(:minus), derivative_order=1, accuracy_order=2,
+                         xmin=0//1, xmax=1//1, N=9)
+Matrix(Dm)
+
+M = mass_matrix(Dp)
+M * Matrix(Dp) + Matrix(Dm)' * M
+minimum(eigvals(-M * (Matrix(Dp) - Matrix(Dm)))) # > 0 up to floating point tolerances
+```
+
+You can also set up fully periodic upwind operators by setting the argument
+`left_offset` of [`periodic_derivative_operator`](@ref) appropriately. For example,
+
+```@repl
+using SummationByPartsOperators, LinearAlgebra
+
+Dp = periodic_derivative_operator(derivative_order=1, accuracy_order=2, left_offset=0,
+                                  xmin=0//1, xmax=1//1, N=8)
+Matrix(Dp)
+Dm = periodic_derivative_operator(derivative_order=1, accuracy_order=2, left_offset=-2,
+                                  xmin=0//1, xmax=1//1, N=8)
+Matrix(Dm)
+
+M = mass_matrix(Dp)
+M * Matrix(Dp) + Matrix(Dm)' * M |> iszero
+minimum(eigvals(-M * (Matrix(Dp) - Matrix(Dm)))) # > 0 up to floating point tolerances
 ```
 
+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.
+
 
 ## Continuous and discontinuous Galerkin methods
 
@@ -302,34 +226,14 @@ If the underlying SBP operators are [`LegendreDerivativeOperator`](@ref)s,
 these are CG spectral element methods (CGSEM). However, a continuous coupling
 of arbitrary SBP operators is supported.
 
-```jldoctest
-julia> using SummationByPartsOperators, LinearAlgebra
+```@repl
+using SummationByPartsOperators, LinearAlgebra
 
-julia> 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.0
-
-julia> mass_matrix(D)
-7×7 Diagonal{Float64, Vector{Float64}}:
- 0.0555556   ⋅         ⋅         ⋅         ⋅         ⋅         ⋅
-  ⋅         0.222222   ⋅         ⋅         ⋅         ⋅         ⋅
-  ⋅          ⋅        0.111111   ⋅         ⋅         ⋅         ⋅
-  ⋅          ⋅         ⋅        0.222222   ⋅         ⋅         ⋅
-  ⋅          ⋅         ⋅         ⋅        0.111111   ⋅         ⋅
-  ⋅          ⋅         ⋅         ⋅         ⋅        0.222222   ⋅
-  ⋅          ⋅         ⋅         ⋅         ⋅         ⋅        0.0555556
+D = couple_continuously(
+        legendre_derivative_operator(xmin=-1.0, xmax=1.0, N=3),
+        UniformMesh1D(xmin=0.0, xmax=1.0, Nx=3))
+Matrix(D)
+mass_matrix(D)
 ```
 
 SBP operators can also be coupled as in discontinuous Galerkin (DG) methods.
@@ -337,20 +241,16 @@ Using a central numerical flux results in central SBP operators; upwind fluxes
 yield upwind SBP operators. If [`LegendreDerivativeOperator`](@ref)s are used,
 the discontinuous coupling yields DG spectral element methods (DGSEM).
 
-```jldoctest
-julia> using SummationByPartsOperators, LinearAlgebra
+```@repl
+using SummationByPartsOperators, LinearAlgebra
 
-julia> 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)
+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))
 
-julia> M = mass_matrix(D);
-
-julia> M * Matrix(D) + Matrix(D)' * M |> iszero
-true
+M = mass_matrix(D);
+M * Matrix(D) + Matrix(D)' * M |> iszero
 ```
 
 Right now, only uniform meshes [`UniformMesh1D`](@ref) and [`UniformPeriodicMesh1D`](@ref)