Skip to content

Commit

Permalink
update docs (#137)
Browse files Browse the repository at this point in the history 
* show versions in benchmarks

* fix typo

* remove doctest filter

* remove doctest filter

* more repl blocks

* mention periodic upwind operators in the introduction

* bump version
  • Loading branch information
@ranocha
ranocha committed Apr 8, 2022
1 parent 4e084b8 commit ca1e56d
Show file tree
Hide file tree
Showing 3 changed files with 119 additions and 175 deletions.
2 changes: 1 addition & 1 deletion Project.toml
View file
Original file line number Diff line number Diff line change
@@ -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"
Expand Down
46 changes: 45 additions & 1 deletion docs/src/benchmarks.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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)

Expand All @@ -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
Expand Down Expand Up @@ -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
```
246 changes: 73 additions & 173 deletions docs/src/introduction.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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).
Expand All @@ -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

Expand All @@ -302,55 +226,31 @@ 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.
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)
Expand Down

2 comments on commit ca1e56d

@ranocha
Copy link
Owner Author

@ranocha ranocha commented on ca1e56d Apr 8, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@JuliaRegistrator register

@JuliaRegistrator
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Registration pull request created: JuliaRegistries/General/58146

After the above pull request is merged, it is recommended that a tag is created on this repository for the registered package version.

This will be done automatically if the Julia TagBot GitHub Action is installed, or can be done manually through the github interface, or via:

git tag -a v0.5.16 -m "<description of version>" ca1e56d961902d8545be51a52303ac70d211f532
git push origin v0.5.16

Please sign in to comment.