Merge pull request #1 from JuliaDiffEq/master

rebase

Author: xtalax

src/DiffEqOperators.jl

@@ -25,8 +25,6 @@ include("derivative_operators/BC_operators.jl")
 
 ### Derivative Operators
 include("derivative_operators/fornberg.jl")
-include("derivative_operators/upwind_operator.jl")
-include("derivative_operators/derivative_irreg_operator.jl")
 include("derivative_operators/derivative_operator.jl")
 include("derivative_operators/abstract_operator_functions.jl")
 include("derivative_operators/convolutions.jl")
@@ -44,6 +42,7 @@ end
 
 export MatrixFreeOperator
 export DiffEqScalar, DiffEqArrayOperator, DiffEqIdentity, JacVecOperator, getops
-export AbstractDerivativeOperator, DerivativeOperator, UpwindOperator, FiniteDifference
-export RobinBC
+export AbstractDerivativeOperator, DerivativeOperator,
+       CenteredDifference, UpwindDifference
+export RobinBC, GeneralBC
 end # module

src/derivative_operators/BC_operators.jl

@@ -1,4 +1,4 @@
-abstract type AbstractBC{T} end
+abstract type AbstractBC{T} <: AbstractDiffEqLinearOperator{T} end
 
 """
 Robin, General, and in general Neumann and Dirichlet BCs are all affine opeartors, meaning that they take the form Qx = Qax + Qb.
@@ -10,7 +10,7 @@ struct PeriodicBC{T} <: AbstractBC{T}
 end
 
 """
-  The variables in l are [al, bl, cl], and correspond to a BC of the form al*u(0) + bl*u'(0) = cl
+  The variables in l are [αl, βl, γl], and correspond to a BC of the form al*u(0) + bl*u'(0) = cl
 
   Implements a robin boundary condition operator Q that acts on a vector to give an extended vector as a result
   Referring to (https://github.com/JuliaDiffEq/DiffEqOperators.jl/files/3267835/ghost_node.pdf)
@@ -27,17 +27,17 @@ struct RobinBC{T, V<:AbstractVector{T}} <: AffineBC{T,V}
     a_r::V
     b_r::T
     function RobinBC(l::AbstractArray{T}, r::AbstractArray{T}, dx::AbstractArray{T}, order = one(T)) where {T}
-        cl, al, bl = l
-        cr, ar, br = r
+        αl, βl, γl = l
+        αr, βr, γr = r
         dx_l, dx_r = dx
 
         s = calculate_weights(1, one(T), Array(one(T):convert(T,order+1))) #generate derivative coefficients about the boundary of required approximation order
 
-        a_l = -s[2:end]./(1+al*dx_l*s[1]/bl)
-        a_r = s[end:-1:2]./(1-ar*dx_r*s[1]/br) # for other boundary stencil is flippedlr with *opposite sign*
+        a_l = -s[2:end]./(αl*dx_l/βl + s[1])
+        a_r = s[end:-1:2]./(αr*dx_r/βr - s[1]) # for other boundary stencil is flippedlr with *opposite sign*
 
-        b_l = cl/(al+bl*s[1]/dx_l)
-        b_r = cr/(ar-br*s[1]/dx_r)
+        b_l = γl/(αl+βl*s[1]/dx_l)
+        b_r = γr/(αr-βr*s[1]/dx_r)
 
         return new{T, typeof(a_l)}(a_l, b_l, a_r, b_r)
     end
@@ -86,20 +86,19 @@ struct GeneralBC{T, V<:AbstractVector{T}} <:AffineBC{T,V}
 end
 
 #implement Neumann and Dirichlet as special cases of RobinBC
-NeumannBC(  α::AbstractVector{T}, dx::AbstractVector{T}, order=1) where T = RobinBC([zero(T), one(T), α[1]], [zero(T), one(T), α[2]], dx, order)
-DirichletBC(α::AbstractVector{T}, dx::AbstractVector{T}, order=1) where T = RobinBC([one(T), zero(T), α[1]], [one(T), zero(T), α[2]], dx, order)
-
+NeumannBC(α::AbstractVector{T}, dx::AbstractVector{T}, order = 1) where T = RobinBC([zero(T), one(T), α[1]], [zero(T), one(T), α[2]], dx, order)
+DirichletBC(α::AbstractVector{T}, dx::AbstractVector{T}, order = 1) where T = RobinBC([one(T), zero(T), α[1]], [one(T), zero(T), α[2]], dx, order)
 # other acceptable argument signatures
-RobinBC(al::T, bl::T, cl::T, dx_l::T, ar::T, br::T, cr::T, dx_r::T, order = 1) where T = RobinBC([cl,al,bl], [cr, ar, br], [dx_l, dx_r], order)
+RobinBC(al::T, bl::T, cl::T, dx_l::T, ar::T, br::T, cr::T, dx_r::T, order = 1) where T = RobinBC([al,bl, cl], [ar, br, cr], [dx_l, dx_r], order)
 
 # this  is 'boundary padded vector' as opposed to 'boundary padded array' to distinguish it from the n dimensional implementation that will eventually be neeeded
-struct BoundaryPaddedVector{T,T2 <: AbstractVector{T}}
+struct BoundaryPaddedVector{T,T2 <: AbstractVector{T}} <: AbstractVector{T}
     l::T
     r::T
     u::T2
 end
 
-Base.:*(Q::AffineBC, u) = BoundaryPaddedVector(Q.a_l ⋅ u[1:length(Q.a_l)] + Q.b_l, Q.a_r ⋅ u[(end-length(Q.a_r)+1):end] + Q.b_r, u)
+Base.:*(Q::AffineBC, u::AbstractVector) = BoundaryPaddedVector(Q.a_l ⋅ u[1:length(Q.a_l)] + Q.b_l, Q.a_r ⋅ u[(end-length(Q.a_r)+1):end] + Q.b_r, u)
 
 Base.size(Q::AbstractBC) = (Inf, Inf) #Is this nessecary?
 Base.length(Q::BoundaryPaddedVector) = length(Q.u) + 2
@@ -140,4 +139,14 @@ function LinearAlgebra.Array(Q::BoundaryPaddedVector)
     return [Q.l; Q.u; Q.r]
 end
 
-#TODO: Implement Sparse concretization
+function Base.convert(::Type{Array},A::AbstractBC{T}) where T
+    Array(A)
+end
+
+function Base.convert(::Type{SparseMatrixCSC},A::AbstractBC{T}) where T
+    SparseMatrixCSC(A)
+end
+
+function Base.convert(::Type{AbstractMatrix},A::AbstractBC{T}) where T
+    SparseMatrixCSC(A)
+end

src/derivative_operators/abstract_operator_functions.jl

@@ -17,52 +17,56 @@ checkbounds(A::AbstractDerivativeOperator, k::Colon, j::Integer) =
 checkbounds(A::AbstractDerivativeOperator, k::Integer, j::Colon) =
     (0 < k ≤ size(A, 1) || throw(BoundsError(A, (k,size(A,2)))))
 
-# ~~ getindex ~~
 @inline function getindex(A::AbstractDerivativeOperator, i::Int, j::Int)
     @boundscheck checkbounds(A, i, j)
-    mid = div(A.stencil_length, 2) + 1
-    bpc = A.stencil_length - mid
-    l = max(1, low(j, mid, bpc))
-    h = min(A.stencil_length, high(j, mid, bpc, A.stencil_length, A.dimension))
-    slen = h - l + 1
-    if abs(i - j) > div(slen, 2)
-        return 0
+    bpc = A.boundary_point_count
+    N = A.len
+    bsl = A.boundary_stencil_length
+    slen = A.stencil_length
+    if bpc > 0 && 1<=i<=bpc
+        if j > bsl
+            return 0
+        else
+            return A.low_boundary_coefs[i][j]
+        end
+    elseif bpc > 0 && (N-bpc)<i<=N
+        if j < N+2-bsl
+            return 0
+        else
+            return A.high_boundary_coefs[bpc-(N-i)][bsl-(N+2-j)]
+        end
     else
-        return A.stencil_coefs[mid + j - i]
+        if j < i-bpc || j > i+slen-bpc-1
+            return 0
+        else
+            return A.stencil_coefs[j-i + 1 + bpc]
+        end
     end
 end
 
 # scalar - colon - colon
-@inline getindex(A::AbstractDerivativeOperator, kr::Colon, jr::Colon) = convert(Array,A)
+@inline getindex(A::AbstractDerivativeOperator, ::Colon, ::Colon) = Array(A)
 
-@inline function getindex(A::AbstractDerivativeOperator, rc::Colon, j)
-    T = eltype(A.stencil_coefs)
-    v = zeros(T, A.dimension)
-    v[j] = one(T)
-    copyto!(v, A*v)
-    return v
+@inline function getindex(A::AbstractDerivativeOperator, ::Colon, j)
+    return BandedMatrix(A)[:,j]
 end
 
 
 # symmetric right now
-@inline function getindex(A::AbstractDerivativeOperator, i, cc::Colon)
-    T = eltype(A.stencil_coefs)
-    v = zeros(T, A.dimension)
-    v[i] = one(T)
-    copyto!(v, A*v)
-    return v
+@inline function getindex(A::AbstractDerivativeOperator, i, ::Colon)
+    return BandedMatrix(A)[i,:]
 end
 
 
 # UnitRanges
-@inline function getindex(A::AbstractDerivativeOperator, rng::UnitRange{Int}, cc::Colon)
-    m = convert(Array,A)
+@inline function getindex(A::AbstractDerivativeOperator, rng::UnitRange{Int}, ::Colon)
+    m = BandedMatrix(A)
     return m[rng, cc]
 end
 
 
-@inline function getindex(A::AbstractDerivativeOperator, rc::Colon, rng::UnitRange{Int})
-    m = convert(Array,A)
+@inline function getindex(A::AbstractDerivativeOperator, ::Colon, rng::UnitRange{Int})
+    m = BandedMatrix(A)
     return m[rnd, cc]
 end
 
@@ -79,35 +83,7 @@ end
 
 
 @inline function getindex(A::AbstractDerivativeOperator{T}, rng::UnitRange{Int}, cng::UnitRange{Int}) where T
-    N = A.dimension
-    if (rng[end] - rng[1]) > ((cng[end] - cng[1]))
-        mat = zeros(T, (N, length(cng)))
-        v = zeros(T, N)
-        for i = cng
-            v[i] = one(T)
-            #=
-                calculating the effect on a unit vector to get the matrix of transformation
-                to get the vector in the new vector space.
-            =#
-            mul!(view(mat, :, i - cng[1] + 1), A, v)
-            v[i] = zero(T)
-        end
-        return mat[rng, :]
-
-    else
-        mat = zeros(T, (length(rng), N))
-        v = zeros(T, N)
-        for i = rng
-            v[i] = one(T)
-            #=
-                calculating the effect on a unit vector to get the matrix of transformation
-                to get the vector in the new vector space.
-            =#
-            mul!(view(mat, i - rng[1] + 1, :), A, v)
-            v[i] = zero(T)
-        end
-        return mat[:, cng]
-    end
+    return BandedMatrix(A)[rng,cng]
 end
 
 #=
@@ -126,74 +102,95 @@ function LinearAlgebra.mul!(x_temp::AbstractArray{T,2}, A::AbstractDerivativeOpe
     end
 end
 
-
 # Base.length(A::AbstractDerivativeOperator) = A.stencil_length
 Base.ndims(A::AbstractDerivativeOperator) = 2
-Base.size(A::AbstractDerivativeOperator) = (A.dimension, A.dimension)
+Base.size(A::AbstractDerivativeOperator) = (A.len, A.len + 2)
 Base.size(A::AbstractDerivativeOperator,i::Integer) = size(A)[i]
 Base.length(A::AbstractDerivativeOperator) = reduce(*, size(A))
 
-
 #=
     For the evenly spaced grid we have a symmetric matrix
 =#
-Base.transpose(A::Union{DerivativeOperator,UpwindOperator}) = A
-Base.adjoint(A::Union{DerivativeOperator,UpwindOperator}) = A
-LinearAlgebra.issymmetric(::Union{DerivativeOperator,UpwindOperator}) = true
+Base.transpose(A::DerivativeOperator) = A
+Base.adjoint(A::DerivativeOperator) = A
+LinearAlgebra.issymmetric(::DerivativeOperator) = true
 
 #=
     Fallback methods that use the full representation of the operator
 =#
-Base.exp(A::AbstractDerivativeOperator{T}) where T = exp(convert(Array,A))
+Base.exp(A::AbstractDerivativeOperator{T}) where T = exp(convert(A))
 Base.:\(A::AbstractVecOrMat, B::AbstractDerivativeOperator) = A \ convert(Array,B)
-Base.:\(A::AbstractDerivativeOperator, B::AbstractVecOrMat) = convert(Array,A) \ B
+Base.:\(A::AbstractDerivativeOperator, B::AbstractVecOrMat) = Array(A) \ B
 Base.:/(A::AbstractVecOrMat, B::AbstractDerivativeOperator) = A / convert(Array,B)
-Base.:/(A::AbstractDerivativeOperator, B::AbstractVecOrMat) = convert(Array,A) / B
+Base.:/(A::AbstractDerivativeOperator, B::AbstractVecOrMat) = Array(A) / B
 
-########################################################################
+#=
+    The Inf opnorm can be calculated easily using the stencil coeffiicents, while other opnorms
+    default to compute from the full matrix form.
+=#
+function LinearAlgebra.opnorm(A::DerivativeOperator, p::Real=2)
+    if p == Inf
+        sum(abs.(A.stencil_coefs)) / A.dx^A.derivative_order
+    else
+        opnorm(BandedMatrix(A), p)
+    end
+end
 
-# Are these necessary?
+########################################################################
 
 get_type(::AbstractDerivativeOperator{T}) where {T} = T
 
 function *(A::AbstractDerivativeOperator,x::AbstractVector)
-    #=
-        We will output a vector which is a supertype of the types of A and x
-        to ensure numerical stability
-    =#
-    get_type(A) != eltype(x) ? error("DiffEqOperator and array are not of same type!") : nothing
-    y = zeros(promote_type(eltype(A),eltype(x)), length(x))
+    y = zeros(promote_type(eltype(A),eltype(x)), length(x)-2)
     LinearAlgebra.mul!(y, A::AbstractDerivativeOperator, x::AbstractVector)
     return y
 end
 
 
 function *(A::AbstractDerivativeOperator,M::AbstractMatrix)
-    #=
-        We will output a vector which is a supertype of the types of A and x
-        to ensure numerical stability
-    =#
-    get_type(A) != eltype(M) ? error("DiffEqOperator and array are not of same type!") : nothing
-    y = zeros(promote_type(eltype(A),eltype(M)), size(M))
+    y = zeros(promote_type(eltype(A),eltype(M)), size(A,1), size(M,2))
     LinearAlgebra.mul!(y, A::AbstractDerivativeOperator, M::AbstractMatrix)
     return y
 end
 
 
 function *(M::AbstractMatrix,A::AbstractDerivativeOperator)
-    #=
-        We will output a vector which is a supertype of the types of A and x
-        to ensure numerical stability
-    =#
-    get_type(A) != eltype(M) ? error("DiffEqOperator and array are not of same type!") : nothing
-    y = zeros(promote_type(eltype(A),eltype(M)), size(M))
-    LinearAlgebra.mul!(y, A::AbstractDerivativeOperator, M::AbstractMatrix)
+    y = zeros(promote_type(eltype(A),eltype(M)), size(M,1), size(A,2))
+    LinearAlgebra.mul!(y, M, BandedMatrix(A))
     return y
 end
 
 
 function *(A::AbstractDerivativeOperator,B::AbstractDerivativeOperator)
-    # TODO: it will result in an operator which calculates
-    #       the derivative of order A.dorder + B.dorder of
-    #       approximation_order = min(approx_A, approx_B)
+    return BandedMatrix(A)*BandedMatrix(B)
+end
+
+################################################################################
+
+function *(coeff_func::Function, A::DerivativeOperator{T,N,Wind}) where {T,N,Wind}
+    coefficients = A.coefficients === nothing ? Vector{T}(undef,A.len) : A.coefficients
+    DerivativeOperator{T,N,Wind,typeof(A.dx),typeof(A.stencil_coefs),
+                       typeof(A.low_boundary_coefs),typeof(coefficients),
+                       typeof(coeff_func)}(
+        A.derivative_order, A.approximation_order,
+        A.dx, A.len, A.stencil_length,
+        A.stencil_coefs,
+        A.boundary_stencil_length,
+        A.boundary_point_count,
+        A.low_boundary_coefs,
+        A.high_boundary_coefs,coefficients,coeff_func
+        )
 end
+
+################################################################################
+
+function DiffEqBase.update_coefficients!(A::AbstractDerivativeOperator,u,p,t)
+    if A.coeff_func !== nothing
+        A.coeff_func(A.coefficients,u,p,t)
+    end
+end
+
+################################################################################
+
+(L::DerivativeOperator)(u,p,t) = L*u
+(L::DerivativeOperator)(du,u,p,t) = mul!(du,L,u)