Merge pull request #149 from xtalax/SplitBC-dev

#148 Changes implemented

Author: ChrisRackauckas

src/DiffEqOperators.jl

@@ -14,6 +14,9 @@ abstract type AbstractMatrixFreeOperator{T} <: AbstractDiffEqLinearOperator{T} e
 include("matrixfree_operators.jl")
 include("jacvec_operators.jl")
 
+### Utilities
+include("utils.jl") 
+
 ### Boundary Padded Arrays
 include("boundary_padded_arrays.jl")
 

src/boundary_padded_arrays.jl

@@ -1,6 +1,7 @@
 # Boundary Padded Arrays
 abstract type AbstractBoundaryPaddedArray{T, N} <: AbstractArray{T, N} end
-
+abstract type AbstractDirectionalBoundaryPaddedArray{T, N, D} <: AbstractBoundaryPaddedArray{T, N} end
+abstract type AbstractComposedBoundaryPaddedArray{T, N} <: AbstractBoundaryPaddedArray{T,N} end
 """
 A vector type that extends a vector u with ghost points at either end
 """
@@ -9,6 +10,7 @@ struct BoundaryPaddedVector{T,T2 <: AbstractVector{T}} <: AbstractBoundaryPadded
     r::T
     u::T2
 end
+
 Base.length(Q::BoundaryPaddedVector) = length(Q.u) + 2
 Base.size(Q::BoundaryPaddedVector) = (length(Q.u) + 2,)
 Base.lastindex(Q::BoundaryPaddedVector) = Base.length(Q)
@@ -27,7 +29,7 @@ end
 Higher dimensional generalization of BoundaryPaddedVector, pads an array of dimension N along the dimension D with 2 Arrays of dimension N-1, stored in lower and upper
 
 """
-struct BoundaryPaddedArray{T, D, N, M, V<:AbstractArray{T, N}, B<: AbstractArray{T, M}} <: AbstractBoundaryPaddedArray{T,N}
+struct BoundaryPaddedArray{T, D, N, M, V<:AbstractArray{T, N}, B<: AbstractArray{T, M}} <: AbstractDirectionalBoundaryPaddedArray{T,N, D}
     lower::B #an array of dimension M = N-1, used to extend the lower index boundary
     upper::B #Ditto for the upper index boundary
     u::V
@@ -59,9 +61,9 @@ function compose(padded_arrays::BoundaryPaddedArray...)
     Ds = getaxis.(padded_arrays)
     (length(padded_arrays) == N) || throw(ArgumentError("The padded_arrays must cover every dimension - make sure that the number of padded_arrays is equal to ndims(u)."))
     for D in Ds
-        length(setdiff(Ds, D)) == (N-1) || throw(ArgumentError("There are multiple Arrays that extend along dimension $D - make sure every dimension has a unique extension"))
+        length(setdiff(Ds, D)) < N || throw(ArgumentError("There are multiple Arrays that extend along dimension $D - make sure every dimension has a unique extension"))
     end
-    reduce((|), fill(padded_arrays[1].u, (length(padded_arrays),)) .== getfield.(padded_arrays, :u)) || throw(ArgumentError("The padded_arrays do not all extend the same u!"))
+    any(fill(padded_arrays[1].u, (length(padded_arrays),)) .== getfield.(padded_arrays, :u)) || throw(ArgumentError("The padded_arrays do not all extend the same u!"))
     padded_arrays = padded_arrays[sortperm([Ds...])]
     lower = [padded_array.lower for padded_array in padded_arrays]
     upper = [padded_array.upper for padded_array in padded_arrays]
@@ -71,7 +73,7 @@ end
 
 # Composed BoundaryPaddedArray
 
-struct ComposedBoundaryPaddedArray{T, N, M, V<:AbstractArray{T, N}, B<: AbstractArray{T, M}} <: AbstractBoundaryPaddedArray{T, N}
+struct ComposedBoundaryPaddedArray{T, N, M, V<:AbstractArray{T, N}, B<: AbstractArray{T, M}} <: AbstractComposedBoundaryPaddedArray{T, N}
     lower::Vector{B}
     upper::Vector{B}
     u::V
@@ -96,8 +98,7 @@ Ax, Ay,... = decompose(A::ComposedBoundaryPaddedArray)
 
 Decomposes a ComposedBoundaryPaddedArray in to components that extend along each dimension individually
 """
-decompose(A::ComposedBoundaryPaddedArray) = Tuple([BoundaryPaddedArray{gettype(A), ndims(A), ndims(A)-1, typeof(lower[1])}(A.lower[i], A.upper[i], A.u) for i in 1:ndims(A)])
-
+decompose(A::ComposedBoundaryPaddedArray) = Tuple([BoundaryPaddedArray{gettype(A), i, ndims(A), ndims(A)-1, typeof(lower[1])}(A.lower[i], A.upper[i], A.u) for i in 1:ndims(A)])
 
 Base.length(Q::AbstractBoundaryPaddedArray) = reduce((*), size(Q))
 Base.firstindex(Q::AbstractBoundaryPaddedArray, d::Int) = 1
@@ -106,9 +107,6 @@ Base.lastindex(Q::AbstractBoundaryPaddedArray, d::Int) = size(Q)[d]
 gettype(Q::AbstractBoundaryPaddedArray{T,N}) where {T,N} = T
 Base.ndims(Q::AbstractBoundaryPaddedArray{T,N}) where {T,N} = N
 
-add_dim(A::AbstractArray, i) = reshape(A, size(A)...,i)
-add_dim(i) = i
-
 function Base.getindex(Q::BoundaryPaddedArray{T,D,N,M,V,B}, _inds::Vararg{Int,N}) where {T,D,N,M,V,B} #supports range and colon indexing!
     inds = [_inds...]
     S = size(Q)
@@ -139,14 +137,14 @@ function Base.getindex(Q::ComposedBoundaryPaddedArray{T, N, M, V, B} , inds::Var
     for (dim, index) in enumerate(inds)
         if index == 1
             _inds = inds[setdiff(1:N, dim)]
-            if (1 ∈ _inds) | reduce((|), S[setdiff(1:N, dim)] .== _inds)
+            if (1 ∈ _inds) | any(S[setdiff(1:N, dim)] .== _inds)
                 return zero(T)
             else
                 return Q.lower[dim][(_inds.-1)...]
             end
         elseif index == S[dim]
             _inds = inds[setdiff(1:N, dim)]
-            if (1 ∈ _inds) | reduce((|), S[setdiff(1:N, dim)] .== _inds)
+            if (1 ∈ _inds) | any(S[setdiff(1:N, dim)] .== _inds)
                 return zero(T)
             else
                 return Q.upper[dim][(_inds.-1)...]

src/derivative_operators/BC_operators.jl

@@ -37,12 +37,12 @@ end
   The non identity part of Qa is qa:= -b`₁/b0 = -β.*s[2:end]/(α+β*s[1]/Δx). The constant part is Qb = γ/(α+β*s[1]/Δx)
   do the same at the other boundary (amounts to a flip of s[2:end], with the other set of boundary coeffs)
 """
-struct RobinBC{T} <: AffineBC{T}
-    a_l::Vector{T}
+struct RobinBC{T, V<:AbstractVector{T}} <: AffineBC{T}
+    a_l::V
     b_l::T
-    a_r::Vector{T}
+    a_r::V
     b_r::T
-    function RobinBC(l::AbstractVector{T}, r::AbstractVector{T}, dx::T, order = 1) where {T}
+    function RobinBC(l::NTuple{3,T}, r::NTuple{3,T}, dx::T, order = 1) where {T}
         αl, βl, γl = l
         αr, βr, γr = r
 
@@ -54,7 +54,7 @@ struct RobinBC{T} <: AffineBC{T}
         b_l = γl/(αl+βl*s[1]/dx)
         b_r = γr/(αr-βr*s[1]/dx)
 
-        return new{T}(a_l, b_l, a_r, b_r)
+        return new{T, typeof(a_l)}(a_l, b_l, a_r, b_r)
     end
     function RobinBC(l::AbstractVector{T}, r::AbstractVector{T}, dx::AbstractVector{T}, order = 1) where {T}
         αl, βl, γl = l
@@ -73,7 +73,7 @@ struct RobinBC{T} <: AffineBC{T}
         b_l = γl/denom_l
         b_r = γr/denom_r
 
-        return new{T}(a_l, b_l, a_r, b_r)
+        return new{T, typeof(a_l)}(a_l, b_l, a_r, b_r)
     end
 end
 
@@ -91,10 +91,10 @@ This time there are multiple stencils for multiple derivative orders - these can
 All components that multiply u(0) are factored out, turns out to only involve the first colum of S, s̄0. The rest of S is denoted S`. the coeff of u(0) is s̄0⋅ᾱ[3:end] + α[2].
 the remaining components turn out to be ᾱ[3:end]⋅(S`ū`) or equivalantly (transpose(ᾱ[3:end])*S`)⋅ū`. Rearranging, a stencil q_a to be dotted with ū` upon extension can readily be found, along with a constant component q_b
 """
-struct GeneralBC{T} <:AffineBC{T}
-    a_l::Vector{T}
+struct GeneralBC{T, L<:AbstractVector{T}, R<:AbstractVector{T}} <:AffineBC{T}
+    a_l::L
     b_l::T
-    a_r::Vector{T}
+    a_r::R
     b_r::T
     function GeneralBC(αl::AbstractVector{T}, αr::AbstractVector{T}, dx::T, order = 1) where {T}
         nl = length(αl)
@@ -116,11 +116,11 @@ struct GeneralBC{T} <:AffineBC{T}
         denomr = αr[2] .+ αr[3:end] ⋅ s0_r
 
         a_l = -transpose(transpose(αl[3:end]) * Sl) ./denoml
-        a_r = -transpose(transpose(αr[3:end]) * Sr) ./denomr
+        a_r = reverse(-transpose(transpose(αr[3:end]) * Sr) ./denomr)
 
         b_l = -αl[1]/denoml
         b_r = -αr[1]/denomr
-        new{T}(a_l,b_l,reverse!(a_r),b_r)
+        new{T, typeof(a_l), typeof(a_r)}(a_l,b_l,a_r,b_r)
     end
 
     function GeneralBC(αl::AbstractVector{T}, αr::AbstractVector{T}, dx::AbstractVector{T}, order = 1) where {T}
@@ -145,22 +145,22 @@ struct GeneralBC{T} <:AffineBC{T}
         denomr = αr[2] .+ αr[3:end] ⋅ s0_r
 
         a_l = -transpose(transpose(αl[3:end]) * Sl) ./denoml
-        a_r = -transpose(transpose(αr[3:end]) * Sr) ./denomr
+        a_r = reverse(-transpose(transpose(αr[3:end]) * Sr) ./denomr)
 
         b_l = -αl[1]/denoml
         b_r = -αr[1]/denomr
-        new{T}(a_l,b_l,reverse!(a_r),b_r)
+        new{T,  typeof(a_l), typeof(a_r)}(a_l,b_l,a_r,b_r)
     end
 end
 
 
 
 #implement Neumann and Dirichlet as special cases of RobinBC
-NeumannBC(α::AbstractVector{T}, dx::Union{AbstractVector{T}, T}, order = 1) where T = RobinBC([zero(T), one(T), α[1]], [zero(T), one(T), α[2]], dx, order)
-DirichletBC(αl::T, αr::T) where T = RobinBC([one(T), zero(T), αl], [one(T), zero(T), αr], 1.0, 2.0 )
+NeumannBC(α::NTuple{2,T}, dx::Union{AbstractVector{T}, T}, order = 1) where T = RobinBC((zero(T), one(T), α[1]), (zero(T), one(T), α[2]), dx, order)
+DirichletBC(αl::T, αr::T) where T = RobinBC((one(T), zero(T), αl), (one(T), zero(T), αr), 1.0, 2.0 )
 #specialized constructors for Neumann0 and Dirichlet0
 Dirichlet0BC(T::Type) = DirichletBC(zero(T), zero(T))
-Neumann0BC(dx::Union{AbstractVector{T}, T}, order = 1) where T = NeumannBC([zero(T), zero(T)], dx, order)
+Neumann0BC(dx::Union{AbstractVector{T}, T}, order = 1) where T = NeumannBC((zero(T), zero(T)), 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([al,bl, cl], [ar, br, cr], dx_l, order)

src/derivative_operators/multi_dim_bc_operators.jl

@@ -11,6 +11,9 @@ abstract type MultiDimensionalBC{T, N} <: AbstractBC{T} end
     u_temp
 end
 
+"""
+slice_rmul lets you multiply each vector like strip of an array `u` with a linear operator `A`, sliced along dimension `dim`
+"""
 function slice_rmul(A::AbstractDiffEqLinearOperator, u::AbstractArray{T,N}, dim::Int) where {T,N}
     @assert N != 1
     u_temp = similar(u)
@@ -42,9 +45,6 @@ function slice_rmul(A::AbstractArray{B,M}, u::AbstractArray{T,N}, dim::Int) wher
     return (lower, upper)
 end
 
-"""
-slicemul is the only limitation on the BCs here being used up to arbitrary dimension, an N dimensional implementation is needed.
-"""
 
 struct MultiDimDirectionalBC{T<:Number, B<:AtomicBC{T}, D, N, M} <: MultiDimensionalBC{T, N}
     BCs::Array{B,M} #dimension M=N-1 array of BCs to extend dimension D
@@ -94,17 +94,16 @@ MultiDimBC(BC::B, s) where {B<:AtomicBC} = Tuple([MultiDimDirectionalBC{gettype(
 
 PeriodicBC{T}(s) where T = MultiDimBC(PeriodicBC{T}(), s)
 
-NeumannBC(α::AbstractVector{T}, dxyz, order, s) where T = RobinBC([zero(T), one(T), α[1]], [zero(T), one(T), α[2]], dxyz, order, s)
-DirichletBC(αl::T, αr::T, s) where T = RobinBC([one(T), zero(T), αl], [one(T), zero(T), αr], [ones(T, si) for si in s], 2.0, s)
+NeumannBC(α::NTuple{2,T}, dxyz, order, s) where T = RobinBC((zero(T), one(T), α[1]), (zero(T), one(T), α[2]), dxyz, order, s)
+DirichletBC(αl::T, αr::T, s) where T = RobinBC((one(T), zero(T), αl), (one(T), zero(T), αr), [ones(T, si) for si in s], 2.0, s)
 
 Dirichlet0BC(T::Type, s) = DirichletBC(zero(T), zero(T), s)
-Neumann0BC(T::Type, dxyz, order, s) = NeumannBC([zero(T), zero(T)], dxyz, order, s)
+Neumann0BC(T::Type, dxyz, order, s) = NeumannBC((zero(T), zero(T)), dxyz, order, s)
 
-RobinBC(l::AbstractVector{T}, r::AbstractVector{T}, dxyz, order, s) where {T} = Tuple([MultiDimDirectionalBC{T, RobinBC{T}, dim, length(s), length(s)-1}(fill(RobinBC(l, r, dxyz[dim], order), s[setdiff(1:length(s), dim)])) for dim in 1:length(s)])
-GeneralBC(αl::AbstractVector{T}, αr::AbstractVector{T}, dxyz, order, s) where {T} = Tuple([MultiDimDirectionalBC{T, GeneralBC{T}, dim, length(s), length(s)-1}(fill(GeneralBC(αl, αr, dxyz[dim], order), s[setdiff(1:length(s), dim)])) for dim in 1:length(s)])
+RobinBC(l::NTuple{3,T}, r::NTuple{3,T}, dxyz, order, s) where {T} = Tuple([MultiDimDirectionalBC{T, RobinBC{T}, dim, length(s), length(s)-1}(fill(RobinBC(l, r, dxyz[dim], order), perpindex(s,dim))) for dim in 1:length(s)])
+GeneralBC(αl::AbstractVector{T}, αr::AbstractVector{T}, dxyz, order, s) where {T} = Tuple([MultiDimDirectionalBC{T, GeneralBC{T}, dim, length(s), length(s)-1}(fill(GeneralBC(αl, αr, dxyz[dim], order),perpindex(s,dim))) for dim in 1:length(s)])
 
 
-perpsize(A::AbstractArray{T,N}, dim::Integer) where {T,N} = size(A)[setdiff(1:N, dim)] #the size of A perpendicular to dim
 
 """
 Q = compose(BCs...)

src/utils.jl

@@ -0,0 +1,24 @@
+
+"""
+A function that creates a tuple of CartesianIndices of unit length and `N` dimensions, one pointing along each dimension
+"""
+function unit_indices(N::Int) #create unit CartesianIndex for each dimension
+    out = Vector{CartesianIndex{N}}(undef, N)
+    null = zeros(Int64, N)
+    for i in 1:N
+        unit_i = copy(null)
+        unit_i[i] = 1
+        out[i] = CartesianIndex(Tuple(unit_i))
+    end
+    Tuple(out)
+end
+
+add_dims(A::AbstractArray, n::Int) = cat(ndims(a) + n, a)
+
+""
+perpindex(A, dim::Integer) = A[setdiff(1:length(A), dim)]
+
+"""
+the size of A perpendicular to dim
+"""
+perpsize(A::AbstractArray{T,N}, dim::Integer) where {T,N} = size(A)[setdiff(1:N, dim)]

test/MultiDimBC_test.jl

@@ -9,7 +9,7 @@ m = 15
 A = rand(n,m)
 
 #Create atomic BC
-q1 = RobinBC([1.0, 2.0, 3.0], [0.0, -1.0, 2.0], 0.1, 4.0)
+q1 = RobinBC((1.0, 2.0, 3.0), (0.0, -1.0, 2.0), 0.1, 4.0)
 q2 = PeriodicBC{Float64}()
 
 BCx = vcat(fill(q1, div(m,2)), fill(q2, m-div(m,2)))  #The size of BCx has to be all size components *except* for x
@@ -44,7 +44,7 @@ o = 12
 A = rand(n,m, o)
 
 #Create atomic BC
-q1 = RobinBC([1.0, 2.0, 3.0], [0.0, -1.0, 2.0], 0.1, 4.0)
+q1 = RobinBC((1.0, 2.0, 3.0), (0.0, -1.0, 2.0), 0.1, 4.0)
 q2 = PeriodicBC{Float64}()
 
 BCx = vcat(fill(q1, (div(m,2), o)), fill(q2, (m-div(m,2), o)))  #The size of BCx has to be all size components *except* for x

test/bc_coeff_compositions.jl

@@ -42,7 +42,7 @@ end
     cr = rand()
     dx = rand()
 
-    Q = RobinBC([al, bl, cl], [ar, br, cr], dx)
+    Q = RobinBC((al, bl, cl), (ar, br, cr), dx)
     N = 20
     L = CenteredDifference(4,4, 1.0, N)
     L2 = CenteredDifference(2,4, 1.0, N)
@@ -112,7 +112,7 @@ end
     N = length(x)
 
     L = CenteredDifference(2, 2, dx, N)
-    Q = RobinBC([1.0, 0.0, 0.0], [1.0, 0.0, 0.0], dx)
+    Q = RobinBC((1.0, 0.0, 0.0), (1.0, 0.0, 0.0), dx)
     A = L*Q
 
     analytic_L = second_derivative_stencil(N) ./ dx^2
@@ -155,7 +155,7 @@ end
     N = length(x)
 
     L = CenteredDifference(2, 2, dx, N)
-    Q = RobinBC([1.0, 0.0, 4.0], [1.0, 0.0, 4.0], dx)
+    Q = RobinBC((1.0, 0.0, 4.0), (1.0, 0.0, 4.0), dx)
     A = L*Q
 
     analytic_L = second_derivative_stencil(N) ./ dx^2
@@ -203,7 +203,7 @@ end
     u = sin.(x)
 
     L = CenteredDifference(4, 4, dx, N)
-    Q = RobinBC([1.0, 0.0, sin(0.0)], [1.0, 0.0, sin(0.2+dx)], dx)
+    Q = RobinBC((1.0, 0.0, sin(0.0)), (1.0, 0.0, sin(0.2+dx)), dx)
     A = L*Q
 
     analytic_L = fourth_deriv_approx_stencil(N) ./ dx^4

test/robin.jl

@@ -12,7 +12,7 @@ cr = rand(5)
 
 # Construct 5 arbitrary RobinBC operators
 for i in 1:5
-    Q = RobinBC([al[i], bl[i], cl[i]], [ar[i], br[i], cr[i]], dx[i])
+    Q = RobinBC((al[i], bl[i], cl[i]), (ar[i], br[i], cr[i]), dx[i])
 
     Q_L, Q_b = Array(Q,5i)
 
@@ -46,7 +46,7 @@ end
 u0 = -4/10
 uend = 125/12
 u = Vector(1.0:10.0)
-Q = RobinBC([1.0, 6.0, 10.0], [1.0, 6.0, 10.0], 1.0, 3)
+Q = RobinBC((1.0, 6.0, 10.0), (1.0, 6.0, 10.0), 1.0, 3)
 urobinextended = Q*u
 @test urobinextended.l ≈ u0
 @test urobinextended.r ≈ uend