Merge pull request #113 from JuliaDiffEq/nonuniform

Generalize the derivative operator

Author: ChrisRackauckas

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/abstract_operator_functions.jl

@@ -20,7 +20,7 @@ checkbounds(A::AbstractDerivativeOperator, k::Integer, j::Colon) =
 @inline function getindex(A::AbstractDerivativeOperator, i::Int, j::Int)
     @boundscheck checkbounds(A, i, j)
     bpc = A.boundary_point_count
-    N = A.dimension
+    N = A.len
     bsl = A.boundary_stencil_length
     slen = A.stencil_length
     if bpc > 0 && 1<=i<=bpc
@@ -104,16 +104,16 @@ end
 
 # Base.length(A::AbstractDerivativeOperator) = A.stencil_length
 Base.ndims(A::AbstractDerivativeOperator) = 2
-Base.size(A::AbstractDerivativeOperator) = (A.dimension, A.dimension + 2)
+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
@@ -164,3 +164,16 @@ end
 function *(A::AbstractDerivativeOperator,B::AbstractDerivativeOperator)
     return BandedMatrix(A)*BandedMatrix(B)
 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)

src/derivative_operators/concretization.jl

@@ -1,4 +1,4 @@
-function LinearAlgebra.Array(A::DerivativeOperator{T}, N::Int=A.dimension) where T
+function LinearAlgebra.Array(A::DerivativeOperator{T}, N::Int=A.len) where T
     L = zeros(T, N, N+2)
     bl = A.boundary_point_count
     stl = A.stencil_length
@@ -16,7 +16,7 @@ function LinearAlgebra.Array(A::DerivativeOperator{T}, N::Int=A.dimension) where
     return L / A.dx^A.derivative_order
 end
 
-function SparseArrays.SparseMatrixCSC(A::DerivativeOperator{T}, N::Int=A.dimension) where T
+function SparseArrays.SparseMatrixCSC(A::DerivativeOperator{T}, N::Int=A.len) where T
     L = spzeros(T, N, N+2)
     bl = A.boundary_point_count
     stl = A.stencil_length
@@ -34,12 +34,11 @@ function SparseArrays.SparseMatrixCSC(A::DerivativeOperator{T}, N::Int=A.dimensi
     return L / A.dx^A.derivative_order
 end
 
-function SparseArrays.sparse(A::AbstractDerivativeOperator{T}, N::Int=A.dimension) where T
+function SparseArrays.sparse(A::AbstractDerivativeOperator{T}, N::Int=A.len) where T
     SparseMatrixCSC(A,N)
 end
 
-function BandedMatrices.BandedMatrix(A::DerivativeOperator{T}, N::Int=A.dimension) where T
-    N = A.dimension
+function BandedMatrices.BandedMatrix(A::DerivativeOperator{T}, N::Int=A.len) where T
     bl = A.boundary_point_count
     stl = A.stencil_length
     bstl = A.boundary_stencil_length

src/derivative_operators/convolutions.jl

@@ -1,5 +1,5 @@
 # mul! done by convolutions
-function LinearAlgebra.mul!(x_temp::AbstractVector{T}, A::Union{DerivativeOperator{T},UpwindOperator{T}}, x::AbstractVector{T}) where T<:Real
+function LinearAlgebra.mul!(x_temp::AbstractVector{T}, A::DerivativeOperator, x::AbstractVector{T}) where T<:Real
     convolve_BC_left!(x_temp, x, A)
     convolve_interior!(x_temp, x, A)
     convolve_BC_right!(x_temp, x, A)
@@ -9,102 +9,101 @@ end
 ################################################
 
 # Against a standard vector, assume already padded and just apply the stencil
-function convolve_interior!(x_temp::AbstractVector{T}, x::AbstractVector{T}, A::DerivativeOperator{T,S}) where {T<:Real,S<:SVector}
+function convolve_interior!(x_temp::AbstractVector{T}, x::AbstractVector{T}, A::DerivativeOperator) where {T<:Real}
     @assert length(x_temp)+2 == length(x)
-    coeffs = A.stencil_coefs
+    stencil = A.stencil_coefs
+    coeff   = A.coefficients
     mid = div(A.stencil_length,2)
     for i in (1+A.boundary_point_count) : (length(x_temp)-A.boundary_point_count)
         xtempi = zero(T)
+        cur_stencil = eltype(stencil) <: AbstractVector ? stencil[i] : stencil
+        cur_coeff   = eltype(coeff)   <: AbstractVector ? coeff[i] : true
+        cur_stencil = cur_coeff < 0 && A.use_winding ? reverse(cur_stencil) : cur_stencil
         for idx in 1:A.stencil_length
-            xtempi += coeffs[idx] * x[i - mid + idx]
+            xtempi += cur_coeff * cur_stencil[idx] * x[i - mid + idx]
         end
         x_temp[i] = xtempi
     end
 end
 
-function convolve_BC_left!(x_temp::AbstractVector{T}, x::AbstractVector{T}, A::DerivativeOperator{T,S}) where {T<:Real,S<:SVector}
-    coeffs = A.low_boundary_coefs
+function convolve_BC_left!(x_temp::AbstractVector{T}, x::AbstractVector{T}, A::DerivativeOperator) where {T<:Real}
+    stencil = A.low_boundary_coefs
+    coeff   = A.coefficients
     for i in 1 : A.boundary_point_count
-        xtempi = coeffs[i][1]*x[1]
+        xtempi = stencil[i][1]*x[1]
+        cur_stencil = stencil[i]
+        cur_coeff   = eltype(coeff)   <: AbstractVector ? coeff[i] : true
+        cur_stencil = cur_coeff < 0 && A.use_winding ? reverse(cur_stencil) : cur_stencil
         for idx in 2:A.boundary_stencil_length
-            xtempi += coeffs[i][idx] * x[idx]
+            xtempi += cur_coeff * cur_stencil[idx] * x[idx]
         end
         x_temp[i] = xtempi
     end
 end
 
-function convolve_BC_right!(x_temp::AbstractVector{T}, x::AbstractVector{T}, A::DerivativeOperator{T,S}) where {T<:Real,S<:SVector}
-    coeffs = A.high_boundary_coefs
+function convolve_BC_right!(x_temp::AbstractVector{T}, x::AbstractVector{T}, A::DerivativeOperator) where {T<:Real}
+    stencil = A.high_boundary_coefs
+    coeff   = A.coefficients
     for i in 1 : A.boundary_point_count
-        xtempi = coeffs[i][end]*x[end]
+        xtempi = stencil[i][end]*x[end]
+        cur_stencil = stencil[i]
+        cur_coeff   = eltype(coeff)   <: AbstractVector ? coeff[i] : true
+        cur_stencil = cur_coeff < 0 && A.use_winding ? reverse(cur_stencil) : cur_stencil
         for idx in (A.boundary_stencil_length-1):-1:1
-            xtempi += coeffs[i][end-idx] * x[end-idx]
+            xtempi += cur_coeff * cur_stencil[end-idx] * x[end-idx]
         end
         x_temp[end-A.boundary_point_count+i] = xtempi
     end
 end
 
 ###########################################
-# NOT NEEDED ANYMORE
 
 # Against A BC-padded vector, specialize the computation to explicitly use the left, right, and middle parts
-function convolve_interior!(x_temp::AbstractVector{T}, _x::BoundaryPaddedVector, A::DerivativeOperator{T,S}) where {T<:Real,S<:SVector}
-    coeffs = A.stencil_coefs
+function convolve_interior!(x_temp::AbstractVector{T}, _x::BoundaryPaddedVector, A::DerivativeOperator) where {T<:Real}
+    stencil = A.stencil_coefs
+    coeff   = A.coefficients
     x = _x.u
     mid = div(A.stencil_length,2) + 1
     # Just do the middle parts
     for i in (1+A.boundary_length) : (length(x_temp)-A.boundary_length)
         xtempi = zero(T)
+        cur_stencil = eltype(stencil) <: AbstractVector ? stencil[i-A.boundary_length] : stencil
+        cur_coeff   = eltype(coeff)   <: AbstractVector ? coeff[i-A.boundary_length] : true
+        cur_stencil = cur_coeff < 0 && A.use_winding ? reverse(cur_stencil) : cur_stencil
         @inbounds for idx in 1:A.stencil_length
-            xtempi += coeffs[idx] * x[i - (mid-idx) + 1]
+            xtempi += cur_coeff * cur_stencil[idx] * x[i - (mid-idx) + 1]
         end
         x_temp[i] = xtempi
     end
 end
 
-function convolve_BC_left!(x_temp::AbstractVector{T}, _x::BoundaryPaddedVector, A::DerivativeOperator{T,S}) where {T<:Real,S<:SVector}
-    coeffs = A.low_boundary_coefs
+function convolve_BC_left!(x_temp::AbstractVector{T}, _x::BoundaryPaddedVector, A::DerivativeOperator) where {T<:Real}
+    stencil = A.low_boundary_coefs
+    coeff   = A.coefficients
     for i in 1 : A.boundary_length
-        xtempi = coeffs[i][1]*x.l
+        xtempi = stencil[i][1]*x.l
+        cur_stencil = stencil[i]
+        cur_coeff   = eltype(coeff)   <: AbstractVector ? coeff[i-A.boundary_length] : true
+        cur_stencil = cur_coeff < 0 && A.use_winding ? reverse(cur_stencil) : cur_stencil
         @inbounds for idx in 2:A.stencil_length
-            xtempi += coeffs[i][idx] * x.u[idx-1]
+            xtempi += cur_coeff * cur_stencil[idx] * x.u[idx-1]
         end
         x_temp[i] = xtempi
     end
 end
 
-function convolve_BC_right!(x_temp::AbstractVector{T}, _x::BoundaryPaddedVector, A::DerivativeOperator{T,S}) where {T<:Real,S<:SVector}
-    coeffs = A.low_boundary_coefs
+function convolve_BC_right!(x_temp::AbstractVector{T}, _x::BoundaryPaddedVector, A::DerivativeOperator) where {T<:Real}
+    stencil = A.low_boundary_coefs
+    coeff   = A.coefficients
     bc_start = length(x.u) - A.stencil_length
     for i in 1 : A.boundary_length
-        xtempi = coeffs[i][end]*x.r
+        xtempi = stencil[i][end]*x.r
+        cur_stencil = stencil[i]
+        cur_coeff   = eltype(coeff)   <: AbstractVector ? coeff[i-A.boundary_length] : true
+        cur_stencil = cur_coeff < 0 && A.use_winding ? reverse(cur_stencil) : cur_stencil
         @inbounds for idx in A.stencil_length:-1:2
-            xtempi += coeffs[i][end-idx] * x.u[end-idx+1]
+            xtempi += cur_coeff * cur_stencil[end-idx] * x.u[end-idx+1]
         end
         x_temp[i] = xtempi
     end
 end
-
-##########################################
-
-#= INTERIOR CONVOLUTION =#
-function convolve_interior!(x_temp::AbstractVector{T}, x::AbstractVector{T}, A::UpwindOperator{T,S,LBC,RBC}) where {T<:Real,S<:SVector,LBC,RBC}
-    N = length(x)
-    stencil_length = length(A.up_stencil_coefs)
-    stencil_rem = 1 - stencil_length%2
-    Threads.@threads for i in A.boundary_point_count[1]+1 : N-A.boundary_point_count[2]
-        xtempi = zero(T)
-        if A.directions[][i] == false
-            @inbounds for j in 1 : length(A.up_stencil_coefs)
-                xtempi += A.up_stencil_coefs[j] * x[i+j-1-stencil_rem]
-            end
-        else
-            @inbounds for j in -length(A.down_stencil_coefs)+1 : 0
-                xtempi += A.down_stencil_coefs[j+stencil_length] * x[i+j+stencil_rem]
-                # println("i = $i, j = $j, s_idx = $(j+stencil_length), x_idx = $(i+j+stencil_rem), $(A.down_stencil_coefs[j+stencil_length]) * $(x[i+j+stencil_rem]), xtempi = $xtempi")
-            end
-        end
-
-        x_temp[i] = xtempi
-    end
-end