Merge pull request #90 from ranocha/hr/DiffEqOperators

clean-up and more benchmarks

docs/src/benchmarks.md

@@ -6,78 +6,105 @@ on virtual machines in the cloud to generate the documentation automatically.
 
 ## First-derivative operators
 
+Periodic domains:
 ```@example
 using BenchmarkTools
-using SummationByPartsOperators
-using BandedMatrices, LinearAlgebra, Random, SparseArrays
+using LinearAlgebra, SparseArrays
+using SummationByPartsOperators, DiffEqOperators
 
 BLAS.set_num_threads(1) # make sure that BLAS is serial to be fair
 
 T = Float64
 xmin, xmax = T(0), T(1)
-N = 10^3
-der_order = 1 # first-derivative operators
-acc_order = 6 # the (interior) order of accuracy is six
-source = MattssonSvärdShoeybi2008()
 
-D_periodic_serial  = periodic_derivative_operator(der_order, acc_order, xmin, xmax, N+1, Val{:serial}())
-D_nonperiodic_serial  = derivative_operator(source, der_order, acc_order, xmin, xmax, N, Val{:serial}())
-D_nonperiodic_sparse  = sparse(D_nonperiodic_serial)
-D_nonperiodic_banded  = BandedMatrix(D_nonperiodic_serial)
+D_SBP = periodic_derivative_operator(derivative_order=1, accuracy_order=2,
+                                     xmin=xmin, xmax=xmax, N=101)
+x = grid(D_SBP)
+D_DEO = CenteredDifference(derivative_order(D_SBP), accuracy_order(D_SBP),
+                           step(x), length(x)) * PeriodicBC(eltype(D_SBP))
 
-Random.seed!(12345)
-u = randn(T, N)
-dest = similar(u)
+D_sparse = sparse(D_SBP)
 
-function doit(D, text, dest, u)
+u = randn(eltype(D_SBP), length(x)); du = similar(u);
+@show D_SBP * u ≈ D_DEO * u ≈ D_sparse * u
+
+function doit(D, text, du, u)
+  println(text)
+  sleep(0.1)
+  show(stdout, MIME"text/plain"(), @benchmark mul!($du, $D, $u))
+  println()
+end
+
+doit(D_SBP, "D_SBP:", du, u)
+doit(D_DEO, "D_DEO:", du, u)
+doit(D_sparse, "D_sparse:", du, u)
+```
+
+
+General domains:
+```@example
+using BenchmarkTools
+using LinearAlgebra, SparseArrays
+using SummationByPartsOperators, BandedMatrices
+
+BLAS.set_num_threads(1) # make sure that BLAS is serial to be fair
+
+T = Float64
+xmin, xmax = T(0), T(1)
+
+D_SBP = derivative_operator(MattssonNordström2004(), derivative_order=1,
+                            accuracy_order=6, xmin=xmin, xmax=xmax, N=10^3)
+D_sparse = sparse(D_SBP)
+D_banded = BandedMatrix(D_SBP)
+
+u = randn(eltype(D_SBP), size(D_SBP, 1)); du = similar(u);
+@show D_SBP * u ≈ D_sparse * u ≈ D_banded * u
+
+function doit(D, text, du, u)
   println(text)
   sleep(0.1)
-  show(stdout, MIME"text/plain"(), @benchmark mul!($dest, $D, $u))
+  show(stdout, MIME"text/plain"(), @benchmark mul!($du, $D, $u))
   println()
 end
 
-doit(D_periodic_serial, "D_periodic_serial:", dest, u)
-doit(D_nonperiodic_serial, "D_nonperiodic_serial:", dest, u)
-doit(D_nonperiodic_sparse, "D_nonperiodic_sparse:", dest, u)
-doit(D_nonperiodic_banded, "D_nonperiodic_banded:", dest, u)
+doit(D_SBP, "D_SBP:", du, u)
+doit(D_sparse, "D_sparse:", du, u)
+doit(D_banded, "D_banded:", du, u)
 ```
 
 
 ## Dissipation operators
 
 ```@example
 using BenchmarkTools
-using SummationByPartsOperators
-using BandedMatrices, LinearAlgebra, Random, SparseArrays
+using LinearAlgebra, SparseArrays
+using SummationByPartsOperators, BandedMatrices
 
 BLAS.set_num_threads(1) # make sure that BLAS is serial to be fair
 
 T = Float64
 xmin, xmax = T(0), T(1)
-N = 10^3
-acc_order = 8
-source_D = MattssonSvärdShoeybi2008()
-source_Di = MattssonSvärdNordström2004()
-
-D_serial  = derivative_operator(source_D, 1, acc_order, xmin, xmax, N, Val{:serial}())
 
-Di_serial  = dissipation_operator(source_Di, D_serial)
-Di_sparse  = sparse(Di_serial)
-Di_full    = Matrix(Di_serial)
+D_SBP = derivative_operator(MattssonNordström2004(), derivative_order=1,
+                            accuracy_order=6, xmin=xmin, xmax=xmax, N=10^3)
+Di_SBP  = dissipation_operator(MattssonSvärdNordström2004(), D_SBP)
+Di_sparse = sparse(Di_SBP)
+Di_banded = BandedMatrix(Di_SBP)
+Di_full   = Matrix(Di_SBP)
 
-Random.seed!(12345)
-u = randn(T, N)
-dest = similar(u)
+u = randn(eltype(D_SBP), size(D_SBP, 1)); du = similar(u);
+@show Di_SBP * u ≈ Di_sparse * u ≈ Di_banded * u ≈ Di_full * u
 
-function doit(D, text, dest, u)
+function doit(D, text, du, u)
   println(text)
   sleep(0.1)
-  show(stdout, MIME"text/plain"(), @benchmark mul!($dest, $D, $u))
+  show(stdout, MIME"text/plain"(), @benchmark mul!($du, $D, $u))
   println()
 end
 
-doit(D_serial, "D_serial:", dest, u)
-doit(Di_serial, "Di_serial:", dest, u)
-doit(Di_sparse, "Di_sparse:", dest, u)
-doit(Di_full, "Di_full:", dest, u)
+doit(D_SBP, "D_SBP:", du, u)
+doit(Di_SBP, "Di_SBP:", du, u)
+doit(Di_sparse, "Di_sparse:", du, u)
+doit(Di_banded, "Di_banded:", du, u)
+doit(Di_full, "Di_full:", du, u)
 ```

src/SBP_operators.jl

@@ -75,11 +75,7 @@ function Base.show(io::IO, coefficients::DerivativeCoefficients)
 end
 
 
-"""
-    mul!(dest::AbstractVector, coefficients::DerivativeCoefficients, u::AbstractVector, α, β)
-
-Compute `α*D*u + β*dest` and store the result in `dest`.
-"""
+# Compute `α*D*u + β*dest` and store the result in `dest`.
 function mul!(dest::AbstractVector, coefficients::DerivativeCoefficients, u::AbstractVector, α, β)
     @unpack left_boundary, right_boundary, lower_coef, central_coef, upper_coef, parallel = coefficients
 
@@ -93,11 +89,7 @@ function mul!(dest::AbstractVector, coefficients::DerivativeCoefficients, u::Abs
     convolve_interior_coefficients!(dest, lower_coef, central_coef, upper_coef, u, α, β, length(left_boundary), length(right_boundary), parallel)
 end
 
-"""
-    mul!(dest::AbstractVector, coefficients::DerivativeCoefficients, u::AbstractVector, α)
-
-Compute `α*D*u` and store the result in `dest`.
-"""
+# Compute `α*D*u` and store the result in `dest`.
 function mul!(dest::AbstractVector, coefficients::DerivativeCoefficients, u::AbstractVector, α)
     @unpack left_boundary, right_boundary, lower_coef, central_coef, upper_coef, parallel = coefficients
 
@@ -433,11 +425,7 @@ end
 
 
 
-"""
-    mul!(dest::AbstractVector, D::DerivativeOperator, u::AbstractVector, α, β)
-
-Compute `α*D*u + β*dest` and store the result in `dest`.
-"""
+# Compute `α*D*u + β*dest` and store the result in `dest`.
 Base.@propagate_inbounds function mul!(dest::AbstractVector, D::DerivativeOperator, u::AbstractVector, α, β)
     @boundscheck begin
         @argcheck size(D, 2) == length(u) DimensionMismatch
@@ -446,11 +434,7 @@ Base.@propagate_inbounds function mul!(dest::AbstractVector, D::DerivativeOperat
     @inbounds mul!(dest, D.coefficients, u, α*D.factor, β)
 end
 
-"""
-    mul!(dest::AbstractVector, D::DerivativeOperator, u::AbstractVector, α)
-
-Compute `α*D*u` and store the result in `dest`.
-"""
+# Compute `α*D*u` and store the result in `dest`.
 Base.@propagate_inbounds function mul!(dest::AbstractVector, D::DerivativeOperator, u::AbstractVector, α)
     @boundscheck begin
         @argcheck size(D, 2) == length(u) DimensionMismatch

src/dissipation_operators.jl

@@ -52,11 +52,7 @@ function Base.show(io::IO, coefficients::VarCoefDerivativeCoefficients)
 end
 
 
-"""
-    mul!(dest::AbstractVector, coefficients::VarCoefDerivativeCoefficients, u::AbstractVector, b::AbstractVector, α, β)
-
-Compute `α*D*u + β*dest` using the coefficients `b` and store the result in `dest`.
-"""
+# Compute `α*D*u + β*dest` using the coefficients `b` and store the result in `dest`.
 function mul!(dest::AbstractVector, coefficients::VarCoefDerivativeCoefficients,
                 u::AbstractVector, b::AbstractVector, α, β)
     @unpack coefficient_cache, parallel = coefficients
@@ -71,11 +67,7 @@ function mul!(dest::AbstractVector, coefficients::VarCoefDerivativeCoefficients,
     convolve_interior_coefficients!(dest, coefficient_cache, u, b, α, β, parallel)
 end
 
-"""
-    mul!(dest::AbstractVector, coefficients::VarCoefDerivativeCoefficients, u::AbstractVector, b::AbstractVector, α)
-
-Compute `α*D*u` using the coefficients `b` and store the result in `dest`.
-"""
+# Compute `α*D*u` using the coefficients `b` and store the result in `dest`.
 function mul!(dest::AbstractVector, coefficients::VarCoefDerivativeCoefficients,
                 u::AbstractVector, b::AbstractVector, α)
     @unpack coefficient_cache, parallel = coefficients
@@ -155,21 +147,20 @@ end
 
 A dissipation operator on a nonperiodic finite difference grid.
 """
-struct DissipationOperator{T,CoefficientCache,Parallel,SourceOfCoefficients,Grid} <: AbstractVariableCoefficientNonperiodicDerivativeOperator{T}
+struct DissipationOperator{T,Coefficients<:VarCoefDerivativeCoefficients{T},Grid} <: AbstractVariableCoefficientNonperiodicDerivativeOperator{T}
     factor::T
-    coefficients::VarCoefDerivativeCoefficients{T,CoefficientCache,Parallel,SourceOfCoefficients}
+    coefficients::Coefficients
     grid::Grid
     b::Vector{T}
 
     function DissipationOperator(factor::T,
-                                 coefficients::VarCoefDerivativeCoefficients{T,CoefficientCache,Parallel,SourceOfCoefficients},
-                                 grid::Grid, b::Vector{T}) where {T,CoefficientCache,Parallel,SourceOfCoefficients,Grid}
+                                 coefficients::Coefficients,
+                                 grid::Grid, b::Vector{T}) where {T,Coefficients<:VarCoefDerivativeCoefficients{T},Grid}
         @argcheck checkbounds(Bool, grid, coefficients.coefficient_cache) DimensionMismatch
         @argcheck length(grid) == length(b)
         factor < 0 && @warn("Negative dissipation strength shouldn't be used.")
 
-        new{T,CoefficientCache,Parallel,SourceOfCoefficients,Grid}(
-            factor,coefficients, grid, b)
+        new{T,Coefficients,Grid}(factor,coefficients, grid, b)
     end
 end
 
@@ -188,11 +179,8 @@ function Base.show(io::IO, D::DissipationOperator{T}) where {T}
 end
 
 
-"""
-    mul!(dest::AbstractVector, D::DissipationOperator, u::AbstractVector, α, β)
 
-Compute `α*D*u + β*dest` and store the result in `dest`.
-"""
+# Compute `α*D*u + β*dest` and store the result in `dest`.
 Base.@propagate_inbounds function mul!(dest::AbstractVector, D::DissipationOperator,
                                        u::AbstractVector, α, β)
     @boundscheck begin
@@ -202,11 +190,7 @@ Base.@propagate_inbounds function mul!(dest::AbstractVector, D::DissipationOpera
     @inbounds mul!(dest, D.coefficients, u, D.b, D.factor*α, β)
 end
 
-"""
-    mul!(dest::AbstractVector, D::DissipationOperator, u::AbstractVector, α)
-
-Compute `α*D*u` and store the result in `dest`.
-"""
+# Compute `α*D*u` and store the result in `dest`.
 Base.@propagate_inbounds function mul!(dest::AbstractVector, D::DissipationOperator,
                                        u::AbstractVector, α)
     @boundscheck begin

src/periodic_operators.jl

@@ -30,11 +30,7 @@ struct PeriodicDerivativeCoefficients{T,LowerOffset,UpperOffset,Parallel,SourceO
 end
 
 
-"""
-    mul!(dest::AbstractVector, coefficients::PeriodicDerivativeCoefficients, u::AbstractVector, α, β)
-
-Compute `α*D*u + β*dest` and store the result in `dest`.
-"""
+# Compute `α*D*u + β*dest` and store the result in `dest`.
 function mul!(dest::AbstractVector, coefficients::PeriodicDerivativeCoefficients{T,LowerOffset,UpperOffset}, u::AbstractVector, α, β) where {T,LowerOffset,UpperOffset}
     @boundscheck begin
         @argcheck length(dest) == length(u) DimensionMismatch
@@ -46,11 +42,7 @@ function mul!(dest::AbstractVector, coefficients::PeriodicDerivativeCoefficients
     convolve_interior_coefficients!(dest, lower_coef, central_coef, upper_coef, u, α, β, parallel)
 end
 
-"""
-    mul!(dest::AbstractVector, coefficients::PeriodicDerivativeCoefficients, u::AbstractVector, α)
-
-Compute `α*D*u` and store the result in `dest`.
-"""
+# Compute `α*D*u` and store the result in `dest`.
 function mul!(dest::AbstractVector, coefficients::PeriodicDerivativeCoefficients{T,LowerOffset,UpperOffset}, u::AbstractVector, α) where {T,LowerOffset,UpperOffset}
     @boundscheck begin
         @argcheck length(dest) == length(u) DimensionMismatch
@@ -628,11 +620,7 @@ function Base.show(io::IO, D::PeriodicDerivativeOperator{T,LowerOffset,UpperOffs
 end
 
 
-"""
-    mul!(dest::AbstractVector, D::PeriodicDerivativeOperator, u::AbstractVector, α, β)
-
-Compute `α*D*u + β*dest` and store the result in `dest`.
-"""
+# Compute `α*D*u + β*dest` and store the result in `dest`.
 Base.@propagate_inbounds function mul!(dest::AbstractVector, D::PeriodicDerivativeOperator, u::AbstractVector, α, β)
     @boundscheck begin
         @argcheck size(D, 2) == length(u) DimensionMismatch
@@ -641,11 +629,7 @@ Base.@propagate_inbounds function mul!(dest::AbstractVector, D::PeriodicDerivati
     @inbounds mul!(dest, D.coefficients, u, α*D.factor, β)
 end
 
-"""
-    mul!(dest::AbstractVector, D::PeriodicDerivativeOperator, u::AbstractVector, α)
-
-Compute `α*D*u` and store the result in `dest`.
-"""
+# Compute `α*D*u` and store the result in `dest`.
 Base.@propagate_inbounds function mul!(dest::AbstractVector, D::PeriodicDerivativeOperator, u::AbstractVector, α)
     @boundscheck begin
         @argcheck size(D, 2) == length(u) DimensionMismatch
@@ -862,21 +846,13 @@ function Base.show(io::IO, Di::PeriodicDissipationOperator{T}) where {T}
 end
 
 
-"""
-    mul!(dest::AbstractVector, D::PeriodicDissipationOperator, u::AbstractVector, α, β)
-
-Compute `α*D*u + β*dest` and store the result in `dest`.
-"""
+# Compute `α*D*u + β*dest` and store the result in `dest`.
 Base.@propagate_inbounds function mul!(dest::AbstractVector, Di::PeriodicDissipationOperator,
                                        u::AbstractVector, α, β)
     mul!(dest, Di.Di, u, Di.factor*α, β)
 end
 
-"""
-    mul!(dest::AbstractVector, D::PeriodicDissipationOperator, u::AbstractVector, α)
-
-Compute `α*D*u` and store the result in `dest`.
-"""
+# Compute `α*D*u` and store the result in `dest`.
 Base.@propagate_inbounds function mul!(dest::AbstractVector, Di::PeriodicDissipationOperator,
                                        u::AbstractVector, α)
     mul!(dest, Di.Di, u, Di.factor*α)

src/var_coef_operators.jl

@@ -44,11 +44,8 @@ function Base.show(io::IO, D::VarCoefDerivativeOperator{T}) where {T}
 end
 
 
-"""
-    mul!(dest::AbstractVector, D::VarCoefDerivativeOperator, u::AbstractVector, α, β)
 
-Compute `α*D*u + β*dest` and store the result in `dest`.
-"""
+# Compute `α*D*u + β*dest` and store the result in `dest`.
 Base.@propagate_inbounds function mul!(dest::AbstractVector, D::VarCoefDerivativeOperator,
                                         u::AbstractVector, α, β)
     @boundscheck begin
@@ -58,11 +55,7 @@ Base.@propagate_inbounds function mul!(dest::AbstractVector, D::VarCoefDerivativ
     @inbounds mul!(dest, D.coefficients, u, D.b, D.factor*α, β)
 end
 
-"""
-    mul!(dest::AbstractVector, D::VarCoefDerivativeOperator, u::AbstractVector, α)
-
-Compute `α*D*u` and store the result in `dest`.
-"""
+# Compute `α*D*u` and store the result in `dest`.
 Base.@propagate_inbounds function mul!(dest::AbstractVector, D::VarCoefDerivativeOperator,
                                         u::AbstractVector, α)
     @boundscheck begin