Merge pull request #91 from ranocha/hr/breaking

breaking changes v0.4.x → v0.5

NEWS.md

@@ -0,0 +1,20 @@
+# Changelog
+
+SummationByPartsOperators.jl follows the interpretation of
+[semantic versioning (semver)](https://julialang.github.io/Pkg.jl/dev/compatibility/#Version-specifier-format-1)
+used in the Julia ecosystem. Notable changes will be documented in this file
+for human readability.
+
+
+## Breaking changes from v0.4.x to v0.5
+
+- Switch from British English to American English consistently, e.g.,
+  `semidiscretise` → `semidiscretize`
+- `add_transpose_derivative_left!` and `add_transpose_derivative_right!`
+  were replaced by the more general functions
+  `mul_transpose_derivative_left!` and `mul_transpose_derivative_right!`,
+  which use the same interface as `mul!`
+- The number of nodes passed to `periodic_central_derivative_operator`, and
+  `periodic_derivative_operator` changed from the number of visualization nodes
+  to the number of compute nodes (= number of visualization nodes minus one),
+  in accordance with `fourier_derivative_operator`

README.md

@@ -32,6 +32,12 @@ is a registered Julia package. Thus, you can install it from the Julia REPL via
 julia> using Pkg; Pkg.add("SummationByPartsOperators")
 ```
 
+If you want to update SummationByPartsOperators.jl, you can use
+```julia
+julia> using Pkg; Pkg.update()
+```
+A brief list of notable changes is available in [`NEWS.md`](NEWS.md).
+
 
 ## Basic examples
 
@@ -42,9 +48,9 @@ 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=21)
+                                        xmin=0.0, xmax=2.0, N=20)
 Periodic 1st derivative operator of order 2 {T=Float64, Parallel=Val{:serial}}
-on a grid in [0.0, 2.0] using 21 nodes,
+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 from
   Fornberg (1998)
   Calculation of Weights in Finite Difference Formulas.

docs/src/benchmarks.md

@@ -18,7 +18,7 @@ T = Float64
 xmin, xmax = T(0), T(1)
 
 D_SBP = periodic_derivative_operator(derivative_order=1, accuracy_order=2,
-                                     xmin=xmin, xmax=xmax, N=101)
+                                     xmin=xmin, xmax=xmax, N=100)
 x = grid(D_SBP)
 D_DEO = CenteredDifference(derivative_order(D_SBP), accuracy_order(D_SBP),
                            step(x), length(x)) * PeriodicBC(eltype(D_SBP))
@@ -41,7 +41,7 @@ doit(D_sparse, "D_sparse:", du, u)
 ```
 
 
-General domains:
+Bounded domains:
 ```@example
 using BenchmarkTools
 using LinearAlgebra, SparseArrays

docs/src/index.md

@@ -31,9 +31,9 @@ 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=21)
+                                        xmin=0.0, xmax=2.0, N=20)
 Periodic 1st derivative operator of order 2 {T=Float64, Parallel=Val{:serial}}
-on a grid in [0.0, 2.0] using 21 nodes,
+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 from
   Fornberg (1998)
   Calculation of Weights in Finite Difference Formulas.

docs/src/introduction.md

@@ -43,9 +43,9 @@ which can be constructed via [`fourier_derivative_operator`](@ref).
 julia> using SummationByPartsOperators, LinearAlgebra
 
 julia> D = periodic_derivative_operator(derivative_order=1, accuracy_order=2,
-                                        xmin=0.0, xmax=2.0, N=21)
+                                        xmin=0.0, xmax=2.0, N=20)
 Periodic 1st derivative operator of order 2 {T=Float64, Parallel=Val{:serial}}
-on a grid in [0.0, 2.0] using 21 nodes,
+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 from
   Fornberg (1998)
   Calculation of Weights in Finite Difference Formulas.
@@ -242,9 +242,9 @@ 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).
 Some procedures imposing boundary conditions weakly require adding the transposed
 boundary derivatives to a grid function, which can be achieved by
-[`add_transpose_derivative_left!`](@ref) and [`add_transpose_derivative_right!`](@ref).
+[`mul_transpose_derivative_left!`](@ref) and [`mul_transpose_derivative_right!`](@ref).
 You can find applications of these operators in the source code of
-[`WaveEquationNonperiodicSemidiscretisation`](@ref).
+[`WaveEquationNonperiodicSemidiscretization`](@ref).
 
 A special case of second-derivative SBP operators are polynomial derivative operators
 on Lobatto-Legendre nodes, implemented in [`legendre_second_derivative_operator`](@ref).
@@ -414,15 +414,15 @@ equations (PDEs). These are shipped with this package and you are encouraged to
 look at their source code to learn more about it.
 
 - Linear scalar advection with variable coefficient:
-  [`VariableLinearAdvectionNonperiodicSemidiscretisation`](@ref)
+  [`VariableLinearAdvectionNonperiodicSemidiscretization`](@ref)
 - Burgers' equation (inviscid):
-  [`BurgersPeriodicSemidiscretisation`](@ref), [`BurgersNonperiodicSemidiscretisation`](@ref)
+  [`BurgersPeriodicSemidiscretization`](@ref), [`BurgersNonperiodicSemidiscretization`](@ref)
 - Scalar conservation law with cubic flux:
-  [`CubicPeriodicSemidiscretisation`](@ref), [`CubicNonperiodicSemidiscretisation`](@ref)
+  [`CubicPeriodicSemidiscretization`](@ref), [`CubicNonperiodicSemidiscretization`](@ref)
 - A scalar conservation law with quartic, non-convex flux:
-  [`QuarticNonconvexPeriodicSemidiscretisation`](@ref)
+  [`QuarticNonconvexPeriodicSemidiscretization`](@ref)
 - The second-order wave equation:
-  [`WaveEquationNonperiodicSemidiscretisation`](@ref)
+  [`WaveEquationNonperiodicSemidiscretization`](@ref)
 
 Some additional examples are included as [Jupyter](https://jupyter.org) notebooks
 in the directory [`notebooks`](https://github.com/ranocha/SummationByPartsOperators.jl/tree/main/notebooks).

docs/src/tutorials/advection_diffusion.md

@@ -48,7 +48,7 @@ Next, we choose first- and second-derivative SBP operators `D1, D2`, evaluate
 the initial data on the grid, and set up the semidiscretization as an ODE problem.
 
 ```@example advection_diffusion
-N = 101 # number of grid points
+N = 100 # number of grid points
 D1 = periodic_derivative_operator(derivative_order=1, accuracy_order=4,
                                   xmin=xmin, xmax=xmax, N=N)
 D2 = periodic_derivative_operator(derivative_order=2, accuracy_order=4,

docs/src/tutorials/variable_linear_advection.md

@@ -18,7 +18,7 @@ is only allowed for inflow boundaries, e.g. ``a(x_{min}) > 0`` at ``x = x_{min}`
 
 [SummationByPartsOperators.jl](https://github.com/ranocha/SummationByPartsOperators.jl)
 includes a pre-built semidiscretization of this equation:
-[`VariableLinearAdvectionNonperiodicSemidiscretisation`](@ref).
+[`VariableLinearAdvectionNonperiodicSemidiscretization`](@ref).
 Have a look at the source code if you want to dig deeper. Below is an example
 demonstrating how to use this semidiscretization.
 
@@ -46,17 +46,17 @@ split_form = Val(false)
 D = derivative_operator(MattssonSvärdShoeybi2008(), 1, interior_order, xmin, xmax, N)
 # whether or not artificial dissipation should be applied: nothing, dissipation_operator(D)
 Di = nothing
-semidisc = VariableLinearAdvectionNonperiodicSemidiscretisation(D, Di, afunc, split_form, left_bc, right_bc)
-ode = semidiscretize(u0func, semidisc, tspan)
+semi = VariableLinearAdvectionNonperiodicSemidiscretization(D, Di, afunc, split_form, left_bc, right_bc)
+ode = semidiscretize(u0func, semi, tspan)
 
 # solve ODE
 sol = solve(ode, SSPRK104(), dt=D.Δx, adaptive=false,
             save_everystep=false)
 
 # visualise the result
 plot(xguide=L"x", yguide=L"u")
-plot!(evaluate_coefficients(sol[1], semidisc), label=L"u_0")
-plot!(evaluate_coefficients(sol[end], semidisc), label=L"u_\mathrm{numerical}")
+plot!(evaluate_coefficients(sol[1], semi), label=L"u_0")
+plot!(evaluate_coefficients(sol[end], semi), label=L"u_\mathrm{numerical}")
 savefig("example_linear_advection.png");
 ```
 

docs/src/tutorials/wave_equation.md

@@ -13,11 +13,11 @@ Consider the linear wave equation
 
 [SummationByPartsOperators.jl](https://github.com/ranocha/SummationByPartsOperators.jl)
 includes a pre-built semidiscretization of this equation:
-[`WaveEquationNonperiodicSemidiscretisation`](@ref).
+[`WaveEquationNonperiodicSemidiscretization`](@ref).
 Have a look at the source code if you want to dig deeper.
 In particular, you can find applications of
 [`derivative_left`](@ref), [`derivative_right`](@ref)
-[`add_transpose_derivative_left!`](@ref), and [`add_transpose_derivative_right!`](@ref).
+[`mul_transpose_derivative_left!`](@ref), and [`mul_transpose_derivative_right!`](@ref).
 Below is an example demonstrating how to use this semidiscretization.
 
 
@@ -38,7 +38,7 @@ right_bc = Val(:HomogeneousDirichlet)
 # setup spatial semidiscretization
 D2 = derivative_operator(MattssonSvärdShoeybi2008(), derivative_order=2,
                          accuracy_order=4, xmin=xmin, xmax=xmax, N=101)
-semi = WaveEquationNonperiodicSemidiscretisation(D2, left_bc, right_bc)
+semi = WaveEquationNonperiodicSemidiscretization(D2, left_bc, right_bc)
 ode = semidiscretize(v0_func, u0_func, semi, tspan)
 
 # solve second-order ODE using a Runge-Kutta-Nyström method
@@ -71,7 +71,7 @@ function create_gif(left_bc::Val{LEFT_BC}, right_bc::Val{RIGHT_BC}) where {LEFT_
 
     D2 = derivative_operator(MattssonSvärdShoeybi2008(), derivative_order=2,
                             accuracy_order=4, xmin=xmin, xmax=xmax, N=101)
-    semi = WaveEquationNonperiodicSemidiscretisation(D2, left_bc, right_bc)
+    semi = WaveEquationNonperiodicSemidiscretization(D2, left_bc, right_bc)
     ode = semidiscretize(v0_func, u0_func, semi, tspan)
 
     sol = solve(ode, DPRKN6(), saveat=range(first(tspan), stop=last(tspan), length=200))

notebooks/Advection_equation.ipynb

@@ -43,20 +43,20 @@
     "# whether a split form should be applied or not\n",
     "split_form = Val(false)\n",
     "\n",
-    "# setup spatial semidiscretisation\n",
+    "# setup spatial semidiscretization\n",
     "D = derivative_operator(MattssonSvärdShoeybi2008(), 1, interior_order, xmin, xmax, N)\n",
     "# whether or not artificial dissipation should be applied: nothing, dissipation_operator(D)\n",
     "Di = nothing\n",
-    "semidisc = VariableLinearAdvectionNonperiodicSemidiscretisation(D, Di, afunc, split_form, left_bc, right_bc)\n",
-    "ode = semidiscretize(u0func, semidisc, tspan)\n",
+    "semi = VariableLinearAdvectionNonperiodicSemidiscretization(D, Di, afunc, split_form, left_bc, right_bc)\n",
+    "ode = semidiscretize(u0func, semi, tspan)\n",
     "\n",
     "# solve ode\n",
     "sol = solve(ode, SSPRK104(), dt=D.Δx, adaptive=false, \n",
     "            saveat=range(first(tspan), stop=last(tspan), length=200))\n",
     "\n",
     "# visualise the result\n",
     "plot(xguide=L\"x\", yguide=L\"u\")\n",
-    "plot!(evaluate_coefficients(sol[end], semidisc), label=\"\")"
+    "plot!(evaluate_coefficients(sol[end], semi), label=\"\")"
    ]
   },
   {
@@ -68,7 +68,7 @@
     "# make a movie\n",
     "anim = Animation()\n",
     "idx = 1\n",
-    "x, u = evaluate_coefficients(sol[idx], semidisc)\n",
+    "x, u = evaluate_coefficients(sol[idx], semi)\n",
     "\n",
     "fig = plot(x, u, xguide=L\"x\", yguide=L\"u\", xlim=extrema(x), ylim=(-1.05, 1.05),\n",
     "           #size=(1024,768), dpi=250,\n",

notebooks/Benchmarks.ipynb

@@ -23,8 +23,8 @@
     "acc_order = 6\n",
     "source = MattssonSvärdShoeybi2008()\n",
     "\n",
-    "D_periodic_serial  = periodic_derivative_operator(der_order, acc_order, xmin, xmax, N+1, Val{:serial}())\n",
-    "D_periodic_threads = periodic_derivative_operator(der_order, acc_order, xmin, xmax, N+1, Val{:threads}())\n",
+    "D_periodic_serial  = periodic_derivative_operator(der_order, acc_order, xmin, xmax, N, Val{:serial}())\n",
+    "D_periodic_threads = periodic_derivative_operator(der_order, acc_order, xmin, xmax, N, Val{:threads}())\n",
     "D_nonperiodic_serial  = derivative_operator(source, der_order, acc_order, xmin, xmax, N, Val{:serial}())\n",
     "D_nonperiodic_sparse  = sparse(D_nonperiodic_serial)\n",
     "D_nonperiodic_banded  = BandedMatrix(D_nonperiodic_serial)\n",

notebooks/Wave_equation.ipynb

@@ -41,19 +41,19 @@
     "interior_order = 4\n",
     "N = 101\n",
     "\n",
-    "# setup spatial semidiscretisation\n",
+    "# setup spatial semidiscretization\n",
     "D2 = derivative_operator(MattssonSvärdShoeybi2008(), 2, interior_order, xmin, xmax, N)\n",
-    "semidisc = WaveEquationNonperiodicSemidiscretisation(D2, left_bc, right_bc)\n",
-    "ode = semidiscretize(du0func, u0func, semidisc, tspan)\n",
+    "semi = WaveEquationNonperiodicSemidiscretization(D2, left_bc, right_bc)\n",
+    "ode = semidiscretize(du0func, u0func, semi, tspan)\n",
     "\n",
     "# solve ode\n",
     "sol = solve(ode, DPRKN6(), dt=0.25*D2.Δx, adaptive=false, \n",
     "            saveat=range(first(tspan), stop=last(tspan), length=200))\n",
     "\n",
     "# visualise the result\n",
     "plot(xguide=L\"x\")\n",
-    "plot!(evaluate_coefficients(sol[end].x[2], semidisc), label=L\"u\")\n",
-    "plot!(evaluate_coefficients(sol[end].x[1], semidisc), label=L\"\\partial_t u\")"
+    "plot!(evaluate_coefficients(sol[end].x[2], semi), label=L\"u\")\n",
+    "plot!(evaluate_coefficients(sol[end].x[1], semi), label=L\"\\partial_t u\")"
    ]
   },
   {
@@ -65,7 +65,7 @@
     "# make a movie\n",
     "anim = Animation()\n",
     "idx = 1\n",
-    "x, u = evaluate_coefficients(sol[idx].x[2], semidisc)\n",
+    "x, u = evaluate_coefficients(sol[idx].x[2], semi)\n",
     "\n",
     "fig = plot(x, u, xguide=L\"x\", yguide=L\"u\", xlim=extrema(x), ylim=(-1., 1.),\n",
     "           #size=(1024,768), dpi=250,\n",
@@ -117,13 +117,13 @@
     "Nx = 151\n",
     "Ny = 151\n",
     "\n",
-    "# setup spatial semidiscretisation\n",
+    "# setup spatial semidiscretization\n",
     "D2x = derivative_operator(MattssonSvärdShoeybi2008(), 2, interior_order, xmin, xmax, Nx)\n",
-    "semidiscx = WaveEquationNonperiodicSemidiscretisation(D2x, xmin_bc, xmax_bc)\n",
+    "semidiscx = WaveEquationNonperiodicSemidiscretization(D2x, xmin_bc, xmax_bc)\n",
     "opx = LinearMap{Float64}((ddu,u)->semidiscx(ddu,zero(u),u,nothing,0.), Nx, ismutating=true) |> sparse\n",
     "\n",
     "D2y = derivative_operator(MattssonSvärdShoeybi2008(), 2, interior_order, ymin, ymax, Ny)\n",
-    "semidiscy = WaveEquationNonperiodicSemidiscretisation(D2y, ymin_bc, ymax_bc)\n",
+    "semidiscy = WaveEquationNonperiodicSemidiscretization(D2y, ymin_bc, ymax_bc)\n",
     "opy = LinearMap{Float64}((ddu,u)->semidiscy(ddu,zero(u),u,nothing,0.), Ny, ismutating=true) |> sparse\n",
     "\n",
     "x = SummationByPartsOperators.grid(D2x)\n",

src/SBP_operators.jl

@@ -355,40 +355,52 @@ end
 
 
 """
-    add_transpose_derivative_left!(u, D::DerivativeOperator, der_order::Val{N}, α)
+    mul!(du, D::DerivativeOperator, u, α=true, β=false)
 
-Add `α` times the transposed `N`-th derivative functional to the grid function `u`
+Efficient in-place version of `du = α * D * u + β * du`. Note that `du` must not
+be aliased with `u`.
+"""
+function mul! end
+
+
+"""
+    mul_transpose_derivative_left!(u, D::DerivativeOperator, der_order::Val{N}, α=true, β=false)
+
+Set the grid function `u` to `α` times the transposed `N`-th derivative functional
+applied to `u` plus `β` times `u` in the domain of the `N`-th derivative functional
 at the left boundary of the grid.
+Thus, the coefficients `α, β` have the same meaning as in [`mul!`](@ref).
 """
-@inline function add_transpose_derivative_left!(u::AbstractVector, D::DerivativeOperator, der_order::Val{N}, α) where {N}
+@inline function mul_transpose_derivative_left!(u::AbstractVector, D::DerivativeOperator, der_order::Val{N}, α=true, β=false) where {N}
     factor = α / D.Δx^N
     coef = D.coefficients.left_boundary_derivatives[N].coef
     for i in eachindex(coef)
-        u[i] += factor * coef[i]
+        u[i] = factor * coef[i] + β * u[i]
     end
 end
 
-@inline function add_transpose_derivative_left!(u::AbstractVector, D::DerivativeOperator, der_order::Val{0}, α)
-    factor = α
-    u[begin] += factor * u[begin]
+@inline function mul_transpose_derivative_left!(u::AbstractVector, D::DerivativeOperator, der_order::Val{0}, α=true, β=false)
+    u[begin] = α * u[begin] + β * u[begin]
     return nothing
 end
 
 """
-    add_transpose_derivative_right!(u, D::DerivativeOperator, der_order::Val{N}, α)
+    mul_transpose_derivative_right!(u, D::DerivativeOperator, der_order::Val{N}, α=true, β=false)
 
-Add `α` times the transposed `N`-th derivative functional to the grid function `u`
+Set the grid function `u` to `α` times the transposed `N`-th derivative functional
+applied to `u` plus `β` times `u` in the domain of the `N`-th derivative functional
 at the right boundary of the grid.
+Thus, the coefficients `α, β` have the same meaning as in [`mul!`](@ref).
 """
-@inline function add_transpose_derivative_right!(u::AbstractVector, D::DerivativeOperator, der_order::Val{N}, α) where {N}
+@inline function mul_transpose_derivative_right!(u::AbstractVector, D::DerivativeOperator, der_order::Val{N}, α=true, β=false) where {N}
     factor = α / D.Δx^N
     coef = D.coefficients.right_boundary_derivatives[N].coef
     for i in eachindex(coef)
         u[end-i+1] += factor * coef[i]
     end
 end
 
-@inline function add_transpose_derivative_right!(u::AbstractVector, D::DerivativeOperator, der_order::Val{0}, α)
+@inline function mul_transpose_derivative_right!(u::AbstractVector, D::DerivativeOperator, der_order::Val{0}, α=true, β=false)
     factor = α
     u[end] += factor * u[end]
     return nothing
@@ -480,7 +492,7 @@ end
 Create a [`DerivativeOperator`](@ref) approximating the `derivative_order`-th derivative
 on a grid between `xmin` and `xmax` with `N` grid points up to order of accuracy
 `accuracy_order`. with coefficients given by `source_of_coefficients`.
-The evaluation of the derivative can be parallised using threads by chosing
+The evaluation of the derivative can be parallized using threads by chosing
 `parallel=Val{:threads}())`.
 """
 function derivative_operator(source_of_coefficients, derivative_order, accuracy_order, xmin, xmax, N, parallel=Val{:serial}())