Merge pull request #159 from shivin9/New_Upwind_Operators

New upwind operators

Author: ChrisRackauckas

src/derivative_operators/abstract_operator_functions.jl

@@ -19,6 +19,11 @@ checkbounds(A::AbstractDerivativeOperator, k::Integer, j::Colon) =
 
 @inline function getindex(A::AbstractDerivativeOperator, i::Int, j::Int)
     @boundscheck checkbounds(A, i, j)
+    # TODO: Implement the lazy version
+    if use_winding(A)
+        return Array(A)[i,j]
+    end
+
     bpc = A.boundary_point_count
     N = A.len
     bsl = A.boundary_stencil_length
@@ -70,7 +75,11 @@ end
     elseif bpc > 0 && (N-bpc)<i<=N
          v[1:bsl]  .= A.high_boundary_coefs[i-(N-1)]
     else
-        v[i-bpc:i-bpc+slen-1] .= A.stencil_coefs
+        if use_winding(A)
+            v[i-bpc+slen:i-bpc+2*slen-1] .= A.stencil_coefs
+        else
+            v[i-bpc:i-bpc+slen-1] .= A.stencil_coefs
+        end
     end
     return v
 end

src/derivative_operators/concretization.jl

@@ -1,21 +1,29 @@
 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
+    stencil_length = A.stencil_length
     bstl = A.boundary_stencil_length
     coeff   = A.coefficients
-    stl_2 = div(stl,2)
-    for i in 1:A.boundary_point_count
+
+    if use_winding(A)
+        stencil_pivot = 1 + A.stencil_length%2
+    else
+        stencil_pivot = div(stencil_length,2)
+    end
+
+    for i in 1:bl
         cur_coeff   = typeof(coeff)   <: AbstractVector ? coeff[i] : coeff isa Number ? coeff : true
         cur_stencil = use_winding(A) && cur_coeff < 0 ? reverse(A.low_boundary_coefs[i]) : A.low_boundary_coefs[i]
         L[i,1:bstl] = cur_coeff * cur_stencil
     end
+
     for i in bl+1:N-bl
         cur_coeff   = typeof(coeff)   <: AbstractVector ? coeff[i] : coeff isa Number ? coeff : true
         stencil     = eltype(A.stencil_coefs) <: AbstractVector ? A.stencil_coefs[i] : A.stencil_coefs
         cur_stencil = use_winding(A) && cur_coeff < 0 ? reverse(stencil) : stencil
-        L[i,i+1-stl_2:i+1+stl_2] = cur_coeff * cur_stencil
+        L[i,i+1-stencil_pivot:i-stencil_pivot+stencil_length] = cur_coeff * cur_stencil
     end
+
     for i in N-bl+1:N
         cur_coeff   = typeof(coeff)   <: AbstractVector ? coeff[i] : coeff isa Number ? coeff : true
         cur_stencil = use_winding(A) && cur_coeff < 0 ? reverse(A.high_boundary_coefs[i-N+bl]) : A.high_boundary_coefs[i-N+bl]
@@ -27,21 +35,30 @@ end
 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
-    stl_2 = div(stl,2)
+    stencil_length = A.stencil_length
+    stencil_pivot = div(stencil_length,2)
     bstl = A.boundary_stencil_length
     coeff   = A.coefficients
+
+    if use_winding(A)
+        stencil_pivot = 1 + A.stencil_length%2
+    else
+        stencil_pivot = div(stencil_length,2)
+    end
+
     for i in 1:A.boundary_point_count
         cur_coeff   = typeof(coeff)   <: AbstractVector ? coeff[i] : coeff isa Number ? coeff : true
         cur_stencil = use_winding(A) && cur_coeff < 0 ? reverse(A.low_boundary_coefs[i]) : A.low_boundary_coefs[i]
         L[i,1:bstl] = cur_coeff * cur_stencil
     end
+
     for i in bl+1:N-bl
         cur_coeff   = typeof(coeff)   <: AbstractVector ? coeff[i] : coeff isa Number ? coeff : true
         stencil     = eltype(A.stencil_coefs) <: AbstractVector ? A.stencil_coefs[i] : A.stencil_coefs
         cur_stencil = use_winding(A) && cur_coeff < 0 ? reverse(stencil) : stencil
-        L[i,i+1-stl_2:i+1+stl_2] = cur_coeff * cur_stencil
+        L[i,i+1-stencil_pivot:i-stencil_pivot+stencil_length] = cur_coeff * cur_stencil
     end
+
     for i in N-bl+1:N
         cur_coeff   = typeof(coeff)   <: AbstractVector ? coeff[i] : coeff isa Number ? coeff : true
         cur_stencil = use_winding(A) && cur_coeff < 0 ? reverse(A.high_boundary_coefs[i-N+bl]) : A.high_boundary_coefs[i-N+bl]
@@ -56,11 +73,16 @@ end
 
 function BandedMatrices.BandedMatrix(A::DerivativeOperator{T}, N::Int=A.len) where T
     bl = A.boundary_point_count
-    stl = A.stencil_length
+    stencil_length = A.stencil_length
     bstl = A.boundary_stencil_length
     coeff   = A.coefficients
-    stl_2 = div(stl,2)
-    L = BandedMatrix{T}(Zeros(N, N+2), (max(stl-3,0,bstl),max(stl-1,0,bstl)))
+    if use_winding(A)
+        stencil_pivot = 1 + A.stencil_length%2
+    else
+        stencil_pivot = div(stencil_length,2)
+    end
+    L = BandedMatrix{T}(Zeros(N, N+2), (max(stencil_length-3,0,bstl),max(stencil_length-1,0,bstl)))
+
     for i in 1:A.boundary_point_count
         cur_coeff   = typeof(coeff)   <: AbstractVector ? coeff[i] : coeff isa Number ? coeff : true
         cur_stencil = use_winding(A) && cur_coeff < 0 ? reverse(A.low_boundary_coefs[i]) : A.low_boundary_coefs[i]
@@ -70,7 +92,7 @@ function BandedMatrices.BandedMatrix(A::DerivativeOperator{T}, N::Int=A.len) whe
         cur_coeff   = typeof(coeff)   <: AbstractVector ? coeff[i] : coeff isa Number ? coeff : true
         stencil     = eltype(A.stencil_coefs) <: AbstractVector ? A.stencil_coefs[i] : A.stencil_coefs
         cur_stencil = use_winding(A) && cur_coeff < 0 ? reverse(stencil) : stencil
-        L[i,i+1-stl_2:i+1+stl_2] = cur_coeff * cur_stencil
+        L[i,i+1-stencil_pivot:i-stencil_pivot+stencil_length] = cur_coeff * cur_stencil
     end
     for i in N-bl+1:N
         cur_coeff   = typeof(coeff)   <: AbstractVector ? coeff[i] : coeff isa Number ? coeff : true
@@ -361,4 +383,4 @@ function BlockBandedMatrices.BandedBlockBandedMatrix(A::DerivativeOperator{T,N},
         B = Kron(BandedMatrix(Eye(n)), BandedMatrix(A))
     end
     return BandedBlockBandedMatrix(B)
-end
+end