From d2ecf6040445f59b59b80f1985e0d11a9e1e6698 Mon Sep 17 00:00:00 2001 From: Hendrik Ranocha Date: Wed, 6 Mar 2024 17:12:19 +0100 Subject: [PATCH] WIP: new tutorial on basic interfaces (#254) * WIP * WIP * also describe 5-argument mul! * integration * finish first draft * bump version --- Project.toml | 2 +- docs/make.jl | 1 + docs/src/index.md | 19 +-- docs/src/introduction.md | 7 +- docs/src/tutorials/basic_interface.md | 163 ++++++++++++++++++++++++++ 5 files changed, 180 insertions(+), 12 deletions(-) create mode 100644 docs/src/tutorials/basic_interface.md diff --git a/Project.toml b/Project.toml index eee59594..53e75885 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.53" +version = "0.5.54" [deps] ArgCheck = "dce04be8-c92d-5529-be00-80e4d2c0e197" diff --git a/docs/make.jl b/docs/make.jl index 92f12bb7..ea2b4d13 100644 --- a/docs/make.jl +++ b/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", diff --git a/docs/src/index.md b/docs/src/index.md index e0543afb..b08b52c8 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -46,8 +46,10 @@ 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. @@ -55,9 +57,9 @@ stencils with 1 nodes to the left, 1 nodes to the right, and coefficients of For 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. diff --git a/docs/src/introduction.md b/docs/src/introduction.md index 173fe76c..8f828664 100644 --- a/docs/src/introduction.md +++ b/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) diff --git a/docs/src/tutorials/basic_interface.md b/docs/src/tutorials/basic_interface.md new file mode 100644 index 00000000..a848ac7f --- /dev/null +++ b/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)) +```