Skip to content

Commit

Permalink
Merge pull request #100 from ranocha/hr/trixi
Browse files Browse the repository at this point in the history
Interfacing with Trixi and general cleanup
  • Loading branch information
ranocha committed Jun 1, 2021
2 parents 5b1ef8b + 01a0b6f commit 8b98546
Show file tree
Hide file tree
Showing 38 changed files with 1,122 additions and 778 deletions.
9 changes: 9 additions & 0 deletions NEWS.md
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,15 @@ used in the Julia ecosystem. Notable changes will be documented in this file
for human readability.


## Changes in the v0.5 lifecycle

#### Deprecated

- The (keyword) argument `parallel::Union{Val{:serial}, Val{:threads}}`
is deprecated in favor of `mode` with possible values
`FastMode()` (default), `SafeMode()`, and `ThreadedMode()`


## Breaking changes from v0.4.x to v0.5

- Switch from British English to American English consistently, e.g.,
Expand Down
3 changes: 1 addition & 2 deletions Project.toml
Original file line number Diff line number Diff line change
Expand Up @@ -5,8 +5,8 @@ version = "0.5.3-pre"

[deps]
ArgCheck = "dce04be8-c92d-5529-be00-80e4d2c0e197"
ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def"
FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
Expand All @@ -23,7 +23,6 @@ Unrolled = "9602ed7d-8fef-5bc8-8597-8f21381861e8"
[compat]
ArgCheck = "1, 2.0"
DiffEqBase = "6"
DiffEqCallbacks = "2.6"
FFTW = "1"
LoopVectorization = "0.12.22"
PolynomialBases = "0.4.5"
Expand Down
37 changes: 15 additions & 22 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -49,14 +49,11 @@ julia> using Plots: plot, plot!

julia> D = periodic_derivative_operator(derivative_order=1, accuracy_order=2,
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 20 nodes,
stencils with 1 nodes to the left, 1 nodes to the right, and coefficients from
Fornberg (1998)
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")
Expand All @@ -78,15 +75,12 @@ 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)
SBP 1st derivative operator of order 2 {T=Float64, Parallel=Val{:serial}}
on a grid in [0.0, 1.0] using 21 nodes
and coefficients given in
Mattsson, Nordström (2004)
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
derivatives.
Journal of Computational Physics 199, pp. 503-540.


julia> x = grid(D); u = exp.(x);

julia> plot(x, D * u, label="numerical")
Expand All @@ -112,18 +106,18 @@ supported.

### Periodic domains

- `periodic_derivative_operator(derivative_order, accuracy_order, xmin, xmax, N)`
- `periodic_derivative_operator(; derivative_order, accuracy_order, xmin, xmax, N)`

These are classical central finite difference operators using `N` nodes on the
interval `[xmin, xmax]`.

- `periodic_derivative_operator(Holoborodko2008(), derivative_order, accuracy_order, xmin, xmax, N)`
- `periodic_derivative_operator(Holoborodko2008(); derivative_order, accuracy_order, xmin, xmax, N)`

These are central finite difference operators using `N` nodes on the
interval `[xmin, xmax]` and the coefficients of
[Pavel Holoborodko](http://www.holoborodko.com/pavel/numerical-methods/numerical-derivative/smooth-low-noise-differentiators/).

- `fourier_derivative_operator(xmin, xmax, N)`
- `fourier_derivative_operator(; xmin, xmax, N)`

Fourier derivative operators are implemented using the fast Fourier transform of
[FFTW.jl](https://github.com/JuliaMath/FFTW.jl).
Expand All @@ -135,14 +129,14 @@ e.g. to solve elliptic problems of the form `u = (D^2 + I) \ f`.

### Finite (nonperiodic) domains

- `derivative_operator(source_of_coefficients, derivative_order, accuracy_order, xmin, xmax, N)`
- `derivative_operator(source_of_coefficients; derivative_order, accuracy_order, xmin, xmax, N)`

Finite difference SBP operators for first and second derivatives can be obtained
by using `MattssonNordström2004()` as `source_of_coefficients`.
Other sources of coefficients are implemented as well. To obtain a full list
of all operators, use `subtypes(SourceOfCoefficients)`.

- `legendre_derivative_operator(xmin, xmax, N)`
- `legendre_derivative_operator(; xmin, xmax, N)`

Use Lobatto Legendre polynomial collocation schemes on `N`, i.e.
polynomials of degree `N-1`, implemented via
Expand Down Expand Up @@ -181,15 +175,14 @@ of the operators. Therefore, some conversion functions are supplied, e.g.
```julia
julia> using SummationByPartsOperators

julia> D = derivative_operator(MattssonNordström2004(), 1, 2, 0., 1., 5)
SBP 1st derivative operator of order 2 {T=Float64, Parallel=Val{:serial}}
on a grid in [0.0, 1.0] using 5 nodes
and coefficients given in
Mattsson, Nordström (2004)
julia> D = derivative_operator(MattssonNordström2004(),
derivative_order=1, accuracy_order=2,
xmin=0.0, xmax=1.0, N=5)
SBP first-derivative operator of order 2 on a grid in [0.0, 1.0] using 5 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.

Journal of Computational Physics 199, pp. 503-540.

julia> Matrix(D)
5×5 Array{Float64,2}:
Expand Down
14 changes: 4 additions & 10 deletions docs/src/index.md
Original file line number Diff line number Diff line change
Expand Up @@ -32,14 +32,11 @@ julia> using Plots: plot, plot!

julia> D = periodic_derivative_operator(derivative_order=1, accuracy_order=2,
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 20 nodes,
stencils with 1 nodes to the left, 1 nodes to the right, and coefficients from
Fornberg (1998)
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")
Expand All @@ -59,15 +56,12 @@ 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)
SBP 1st derivative operator of order 2 {T=Float64, Parallel=Val{:serial}}
on a grid in [0.0, 1.0] using 21 nodes
and coefficients given in
Mattsson, Nordström (2004)
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
derivatives.
Journal of Computational Physics 199, pp. 503-540.


julia> x = grid(D); u = exp.(x);

julia> plot(x, D * u, label="numerical")
Expand Down
56 changes: 17 additions & 39 deletions docs/src/introduction.md
Original file line number Diff line number Diff line change
Expand Up @@ -44,14 +44,11 @@ julia> using SummationByPartsOperators, LinearAlgebra
julia> D = periodic_derivative_operator(derivative_order=1, accuracy_order=2,
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 20 nodes,
stencils with 1 nodes to the left, 1 nodes to the right, and coefficients from
Fornberg (1998)
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> M = mass_matrix(D)
UniformScaling{Float64}
0.1*I
Expand All @@ -61,8 +58,7 @@ julia> M * Matrix(D) + Matrix(D)' * M |> norm
julia> D = fourier_derivative_operator(xmin=0.0, xmax=2.0, N=20)
Periodic 1st derivative Fourier operator {T=Float64}
on a grid in [0.0, 2.0] using 20 nodes and 11 modes.
on a grid in [0.0, 2.0] using 20 nodes and 11 modes
julia> M = mass_matrix(D)
UniformScaling{Float64}
Expand Down Expand Up @@ -111,15 +107,12 @@ 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 1st derivative operator of order 2 {T=Rational{Int64}, Parallel=Val{:serial}}
on a grid in [0//1, 1//1] using 9 nodes
and coefficients given in
Mattsson, Nordström (2004)
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
Expand Down Expand Up @@ -186,15 +179,12 @@ julia> using SummationByPartsOperators, LinearAlgebra
julia> D = derivative_operator(MattssonNordström2004(), derivative_order=2, accuracy_order=2,
xmin=0//1, xmax=1//1, N=9)
SBP 2nd derivative operator of order 2 {T=Rational{Int64}, Parallel=Val{:serial}}
on a grid in [0//1, 1//1] using 9 nodes
and coefficients given in
Mattsson, Nordström (2004)
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.
julia> M = mass_matrix(D)
9×9 Diagonal{Rational{Int64}, Vector{Rational{Int64}}}:
1//16 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
Expand Down Expand Up @@ -262,25 +252,19 @@ julia> using SummationByPartsOperators, LinearAlgebra
julia> Dp = derivative_operator(Mattsson2017(:plus), derivative_order=1, accuracy_order=2,
xmin=0//1, xmax=1//1, N=9)
SBP 1st derivative operator of order 2 {T=Rational{Int64}, Parallel=Val{:serial}}
on a grid in [0//1, 1//1] using 9 nodes
and coefficients given in
Upwind coefficients (plus) of
Mattsson (2017)
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 1st derivative operator of order 2 {T=Rational{Int64}, Parallel=Val{:serial}}
on a grid in [0//1, 1//1] using 9 nodes
and coefficients given in
Upwind coefficients (minus) of
Mattsson (2017)
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}}}:
Expand Down Expand Up @@ -324,11 +308,8 @@ julia> 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 the Lobatto Legendre nodes in [-1.0, 1.0] using 3 nodes
coupled continuously on the mesh
UniformMesh1D{Float64} with 3 cells in (0.0, 1.0)
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}:
Expand Down Expand Up @@ -363,11 +344,8 @@ 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 the Lobatto Legendre nodes in [-1.0, 1.0] using 3 nodes
coupled discontinuously (upwind: Val{:central}()) on the mesh
UniformPeriodicMesh1D{Float64} with 3 cells in (0.0, 1.0)
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)
julia> M = mass_matrix(D);
Expand Down
40 changes: 22 additions & 18 deletions src/SBP_coefficients/Mattsson2012.jl
Original file line number Diff line number Diff line change
Expand Up @@ -3,37 +3,41 @@
Mattsson2012()
Coefficients of the SBP operators given in
Mattsson (2012)
- Mattsson (2012)
Summation by Parts Operators for Finite Difference Approximations of
Second-Derivatives with Variable Coefficients.
Journal of Scientific Computing 51, pp. 650-682.
"""
struct Mattsson2012 <: SourceOfCoefficients end

function Base.show(io::IO, ::Mattsson2012)
print(io,
" Mattsson (2012) \n",
" Summation by Parts Operators for Finite Difference Approximations of\n",
" Second-Derivatives with Variable Coefficients. \n",
" Journal of Scientific Computing 51, pp. 650-682. \n",
"See also (first derivatives) \n",
" Mattsson, Nordström (2004) \n",
" Summation by parts operators for finite difference approximations of second \n",
" derivatives. \n",
" Journal of Computational Physics 199, pp. 503-540. \n")
function Base.show(io::IO, source::Mattsson2012)
if get(io, :compact, false)
summary(io, source)
else
print(io,
"Mattsson (2012) \n",
" Summation by Parts Operators for Finite Difference Approximations of\n",
" Second-Derivatives with Variable Coefficients. \n",
" Journal of Scientific Computing 51, pp. 650-682. \n",
"See also (first derivatives) \n",
" Mattsson, Nordström (2004) \n",
" Summation by parts operators for finite difference approximations of second \n",
" derivatives. \n",
" Journal of Computational Physics 199, pp. 503-540.")
end
end


@inline function first_derivative_coefficients(source::Mattsson2012, order::Int, T=Float64, parallel=Val{:serial}())
first_derivative_coefficients(MattssonNordström2004(), order, T, parallel)
@inline function first_derivative_coefficients(source::Mattsson2012, order::Int, T=Float64, mode=FastMode())
first_derivative_coefficients(MattssonNordström2004(), order, T, mode)
end

@inline function second_derivative_coefficients(source::Mattsson2012, order::Int, T=Float64, parallel=Val{:serial}())
second_derivative_coefficients(MattssonNordström2004(), order, T, parallel)
@inline function second_derivative_coefficients(source::Mattsson2012, order::Int, T=Float64, mode=FastMode())
second_derivative_coefficients(MattssonNordström2004(), order, T, mode)
end


function var_coef_derivative_coefficients(source::Mattsson2012, derivative_order::Int, accuracy_order::Int, grid, parallel=Val{:serial}())
function var_coef_derivative_coefficients(source::Mattsson2012, derivative_order::Int, accuracy_order::Int, grid, mode=FastMode())
@argcheck derivative_order == 2
T = eltype(grid)
if accuracy_order == 2
Expand Down Expand Up @@ -61,7 +65,7 @@ function var_coef_derivative_coefficients(source::Mattsson2012, derivative_order
end

VarCoefDerivativeCoefficients(coefficient_cache, left_weights, right_weights,
parallel, derivative_order, accuracy_order, source)
mode, derivative_order, accuracy_order, source)
end


Expand Down
Loading

0 comments on commit 8b98546

Please sign in to comment.