Merge pull request #279 from grahamas/local-rebased

Multidimensional handling for GhostDerivativeOperators (part 3)

Author: ChrisRackauckas

src/DiffEqOperators.jl

@@ -7,6 +7,7 @@ import DiffEqBase: AbstractDiffEqLinearOperator, update_coefficients!, isconstan
 using SparseArrays, ForwardDiff, BandedMatrices, NNlib, LazyArrays, BlockBandedMatrices
 using LazyBandedMatrices, ModelingToolkit
 
+abstract type AbstractDiffEqAffineOperator{T} end
 abstract type AbstractDerivativeOperator{T} <: AbstractDiffEqLinearOperator{T} end
 abstract type AbstractDiffEqCompositeOperator{T} <: AbstractDiffEqLinearOperator{T} end
 abstract type AbstractMatrixFreeOperator{T} <: AbstractDiffEqLinearOperator{T} end
@@ -22,15 +23,14 @@ include("utils.jl")
 include("boundary_padded_arrays.jl")
 
 ### Boundary Operators
-include("derivative_operators/BC_operators.jl")
+include("derivative_operators/bc_operators.jl")
 include("derivative_operators/multi_dim_bc_operators.jl")
 
 ### Derivative Operators
 include("derivative_operators/fornberg.jl")
 include("derivative_operators/derivative_operator.jl")
 include("derivative_operators/abstract_operator_functions.jl")
 include("derivative_operators/convolutions.jl")
-include("derivative_operators/concretization.jl")
 include("derivative_operators/ghost_derivative_operator.jl")
 include("derivative_operators/derivative_operator_functions.jl")
 include("derivative_operators/coefficient_functions.jl")
@@ -42,8 +42,11 @@ include("MOL_discretization.jl")
 
 include("docstrings.jl")
 
+### Concretizations
+include("derivative_operators/concretization.jl")
+
 # The (u,p,t) and (du,u,p,t) interface
-for T in [DiffEqScaledOperator, DiffEqOperatorCombination, DiffEqOperatorComposition]
+for T in [DiffEqScaledOperator, DiffEqOperatorCombination, DiffEqOperatorComposition, GhostDerivativeOperator]
   (L::T)(u,p,t) = (update_coefficients!(L,u,p,t); L * u)
   (L::T)(du,u,p,t) = (update_coefficients!(L,u,p,t); mul!(du,L,u))
 end
@@ -53,8 +56,9 @@ export AnalyticalJacVecOperator, JacVecOperator, getops
 export AbstractDerivativeOperator, DerivativeOperator,
        CenteredDifference, UpwindDifference
 export DirichletBC, Dirichlet0BC, NeumannBC, Neumann0BC, RobinBC, GeneralBC, MultiDimBC, PeriodicBC,
-       MultiDimDirectionalBC, ComposedMultiDimBC,
-       compose, decompose, perpsize
+       MultiDimDirectionalBC, ComposedMultiDimBC
+export compose, decompose, perpsize
+
 export GhostDerivativeOperator
 export MOLFiniteDifference
 end # module

src/boundary_padded_arrays.jl

@@ -18,6 +18,8 @@ end
 
 Base.length(Q::BoundaryPaddedVector) = length(Q.u) + 2
 Base.size(Q::BoundaryPaddedVector) = (length(Q.u) + 2,)
+unpadded_size(Q::BoundaryPaddedVector) = size(Q.u)
+unpadded_size(Q::AbstractArray) = size(Q)
 Base.lastindex(Q::BoundaryPaddedVector) = Base.length(Q)
 
 function Base.getindex(Q::BoundaryPaddedVector,i)
@@ -47,6 +49,7 @@ function Base.size(Q::BoundaryPaddedArray)
     S[getaxis(Q)] += 2
     return Tuple(S)
 end
+unpadded_size(Q::BoundaryPaddedArray) = size(Q.u)
 
 """
 A = compose(padded_arrays::BoundaryPaddedArray...)
@@ -95,6 +98,7 @@ ComposedBoundaryPaddedMatrix{T,V,B} = ComposedBoundaryPaddedArray{T,2,1,V,B}
 ComposedBoundaryPadded3Tensor{T,V,B} = ComposedBoundaryPaddedArray{T,3,2,V,B}
 
 Base.size(Q::ComposedBoundaryPaddedArray) = size(Q.u).+2
+unpadded_size(Q::ComposedBoundaryPaddedArray) = size(Q.u)
 
 """
 Ax, Ay,... = decompose(A::ComposedBoundaryPaddedArray)

src/composite_operators.jl

@@ -17,6 +17,7 @@ struct DiffEqScaledOperator{T,F,OpType<:AbstractDiffEqLinearOperator{T}} <: Abst
   op::OpType
 end
 *(α::DiffEqScalar{T,F}, L::AbstractDiffEqLinearOperator{T}) where {T,F} = DiffEqScaledOperator(α, L)
+*(α::Number, L::AbstractDiffEqLinearOperator{T}) where T = DiffEqScaledOperator(DiffEqScalar(convert(T,α)), L)
 -(L::AbstractDiffEqLinearOperator{T}) where {T} = DiffEqScalar(-one(T)) * L
 getops(L::DiffEqScaledOperator) = (L.coeff, L.op)
 Matrix(L::DiffEqScaledOperator) = L.coeff * Matrix(L.op)
@@ -27,15 +28,19 @@ opnorm(L::DiffEqScaledOperator, p::Real=2) = abs(L.coeff) * opnorm(L.op, p)
 getindex(L::DiffEqScaledOperator, i::Int) = L.coeff * L.op[i]
 getindex(L::DiffEqScaledOperator, I::Vararg{Int, N}) where {N} =
   L.coeff * L.op[I...]
-*(L::DiffEqScaledOperator, x::AbstractVecOrMat) = L.coeff * (L.op * x)
-*(x::AbstractVecOrMat, L::DiffEqScaledOperator) = (L.op * x) * L.coeff
-/(L::DiffEqScaledOperator, x::AbstractVecOrMat) = L.coeff * (L.op / x)
-/(x::AbstractVecOrMat, L::DiffEqScaledOperator) = 1/L.coeff * (x / L.op)
-\(L::DiffEqScaledOperator, x::AbstractVecOrMat) = 1/L.coeff * (L.op \ x)
-\(x::AbstractVecOrMat, L::DiffEqScaledOperator) = L.coeff * (x \ L)
-mul!(Y::AbstractVecOrMat, L::DiffEqScaledOperator, B::AbstractVecOrMat) =
-  lmul!(L.coeff, mul!(Y, L.op, B))
-ldiv!(Y::AbstractVecOrMat, L::DiffEqScaledOperator, B::AbstractVecOrMat) =
+*(L::DiffEqScaledOperator, x::AbstractArray) = L.coeff * (L.op * x)
+*(x::AbstractArray, L::DiffEqScaledOperator) = (L.op * x) * L.coeff
+/(L::DiffEqScaledOperator, x::AbstractArray) = L.coeff * (L.op / x)
+/(x::AbstractArray, L::DiffEqScaledOperator) = 1/L.coeff * (x / L.op)
+\(L::DiffEqScaledOperator, x::AbstractArray) = 1/L.coeff * (L.op \ x)
+\(x::AbstractArray, L::DiffEqScaledOperator) = L.coeff * (x \ L)
+for N in (2,3)
+  @eval begin
+    mul!(Y::AbstractArray{T,$N}, L::DiffEqScaledOperator{T}, B::AbstractArray{T,$N}) where {T} =
+        lmul!(Y, L.coeff, mul!(Y, L.op, B))
+  end
+end
+ldiv!(Y::AbstractArray, L::DiffEqScaledOperator, B::AbstractArray) =
   lmul!(1/L.coeff, ldiv!(Y, L.op, B))
 factorize(L::DiffEqScaledOperator) = L.coeff * factorize(L.op)
 for fact in (:lu, :lu!, :qr, :qr!, :cholesky, :cholesky!, :ldlt, :ldlt!,
@@ -46,18 +51,19 @@ end
 
 # Linear Combination
 struct DiffEqOperatorCombination{T,O<:Tuple{Vararg{AbstractDiffEqLinearOperator{T}}},
-  C<:AbstractVector{T}} <: AbstractDiffEqCompositeOperator{T}
-  ops::O
-  cache::C
-  function DiffEqOperatorCombination(ops; cache=nothing)
-    T = eltype(ops[1])
-    if cache == nothing
-      cache = Vector{T}(undef, size(ops[1], 1))
-      fill!(cache,0)
+    C<:AbstractVector{T}} <: AbstractDiffEqCompositeOperator{T}
+    ops::O
+    cache::C
+    function DiffEqOperatorCombination(ops; cache=nothing)
+        T = eltype(ops[1])
+        for i in 2:length(ops)
+            @assert size(ops[i]) == size(ops[1]) "Operators must be of the same size to be combined! Mismatch between $(ops[i]) and $(ops[i-1]), which are operators $i and $(i-1) respectively"
+        end
+        if cache == nothing
+            cache = zeros(T, size(ops[1], 1))
+        end
+        new{T,typeof(ops),typeof(cache)}(ops, cache)
     end
-    # TODO: safecheck dimensions
-    new{T,typeof(ops),typeof(cache)}(ops, cache)
-  end
 end
 +(ops::AbstractDiffEqLinearOperator...) = DiffEqOperatorCombination(ops)
 +(L1::DiffEqOperatorCombination, L2::AbstractDiffEqLinearOperator) = DiffEqOperatorCombination((L1.ops..., L2))
@@ -72,10 +78,10 @@ convert(::Type{AbstractMatrix}, L::DiffEqOperatorCombination) =
 size(L::DiffEqOperatorCombination, args...) = size(L.ops[1], args...)
 getindex(L::DiffEqOperatorCombination, i::Int) = sum(op -> op[i], L.ops)
 getindex(L::DiffEqOperatorCombination, I::Vararg{Int, N}) where {N} = sum(op -> op[I...], L.ops)
-*(L::DiffEqOperatorCombination, x::AbstractVecOrMat) = sum(op -> op * x, L.ops)
-*(x::AbstractVecOrMat, L::DiffEqOperatorCombination) = sum(op -> x * op, L.ops)
-/(L::DiffEqOperatorCombination, x::AbstractVecOrMat) = sum(op -> op / x, L.ops)
-\(x::AbstractVecOrMat, L::DiffEqOperatorCombination) = sum(op -> x \ op, L.ops)
+*(L::DiffEqOperatorCombination, x::AbstractArray) = sum(op -> op * x, L.ops)
+*(x::AbstractArray, L::DiffEqOperatorCombination) = sum(op -> x * op, L.ops)
+/(L::DiffEqOperatorCombination, x::AbstractArray) = sum(op -> op / x, L.ops)
+\(x::AbstractArray, L::DiffEqOperatorCombination) = sum(op -> x \ op, L.ops)
 function mul!(y::AbstractVector, L::DiffEqOperatorCombination, b::AbstractVector)
   mul!(y, L.ops[1], b)
   for op in L.ops[2:end]
@@ -92,7 +98,10 @@ struct DiffEqOperatorComposition{T,O<:Tuple{Vararg{AbstractDiffEqLinearOperator{
   caches::C
   function DiffEqOperatorComposition(ops; caches=nothing)
     T = eltype(ops[1])
-    # TODO: safecheck dimensions
+    for i in 2:length(ops)
+      @assert size(ops[i-1], 1) == size(ops[i], 2) "Operations do not have compatable sizes! Mismatch between $(ops[i]) and $(ops[i-1]), which are operators $i and $(i-1) respectively."
+    end
+
     if caches == nothing
       # Construct a list of caches to be used by mul! and ldiv!
       caches = []
@@ -122,12 +131,12 @@ convert(::Type{AbstractMatrix}, L::DiffEqOperatorComposition) =
 size(L::DiffEqOperatorComposition) = (size(L.ops[end], 1), size(L.ops[1], 2))
 size(L::DiffEqOperatorComposition, m::Integer) = size(L)[m]
 opnorm(L::DiffEqOperatorComposition) = prod(opnorm, L.ops)
-*(L::DiffEqOperatorComposition, x::AbstractVecOrMat) = foldl((acc, op) -> op*acc, L.ops; init=x)
-*(x::AbstractVecOrMat, L::DiffEqOperatorComposition) = foldl((acc, op) -> acc*op, reverse(L.ops); init=x)
-/(L::DiffEqOperatorComposition, x::AbstractVecOrMat) = foldl((acc, op) -> op/acc, L.ops; init=x)
-/(x::AbstractVecOrMat, L::DiffEqOperatorComposition) = foldl((acc, op) -> acc/op, L.ops; init=x)
-\(L::DiffEqOperatorComposition, x::AbstractVecOrMat) = foldl((acc, op) -> op\acc, reverse(L.ops); init=x)
-\(x::AbstractVecOrMat, L::DiffEqOperatorComposition) = foldl((acc, op) -> acc\op, reverse(L.ops); init=x)
+*(L::DiffEqOperatorComposition, x::AbstractArray) = foldl((acc, op) -> op*acc, L.ops; init=x)
+*(x::AbstractArray, L::DiffEqOperatorComposition) = foldl((acc, op) -> acc*op, reverse(L.ops); init=x)
+/(L::DiffEqOperatorComposition, x::AbstractArray) = foldl((acc, op) -> op/acc, L.ops; init=x)
+/(x::AbstractArray, L::DiffEqOperatorComposition) = foldl((acc, op) -> acc/op, L.ops; init=x)
+\(L::DiffEqOperatorComposition, x::AbstractArray) = foldl((acc, op) -> op\acc, reverse(L.ops); init=x)
+\(x::AbstractArray, L::DiffEqOperatorComposition) = foldl((acc, op) -> acc\op, reverse(L.ops); init=x)
 function mul!(y::AbstractVector, L::DiffEqOperatorComposition, b::AbstractVector)
   mul!(L.caches[1], L.ops[1], b)
   for i in 2:length(L.ops) - 1

src/derivative_operators/BC_operators.jl src/derivative_operators/bc_operators.jlsrc/derivative_operators/BC_operators.jl renamed to src/derivative_operators/bc_operators.jl

@@ -1,4 +1,4 @@
-abstract type AbstractBC{T} <: AbstractDiffEqLinearOperator{T} end
+abstract type AbstractBC{T} <: AbstractDiffEqAffineOperator{T} end
 
 
 abstract type AtomicBC{T} <: AbstractBC{T} end
@@ -14,9 +14,32 @@ struct NeumannBC{N} end
 struct Neumann0BC{N} end
 struct DirichletBC{N} end
 struct Dirichlet0BC{N} end
+
+"""
+q = PeriodicBC{T}()
+Qx, Qy, ... = PeriodicBC{T}(size(u)) #When all dimensions are to be extended with a periodic boundary condition.
+-------------------------------------------------------------------------------------
+Creates a periodic boundary condition, where the lower index end of some u is extended with the upper index end and vice versa.
+It is not reccomended to concretize this BC type in to a BandedMatrix, since the vast majority of bands will be all 0s. SpatseMatrix concretization is reccomended.
+"""
 struct PeriodicBC{T} <: AtomicBC{T}
     PeriodicBC(T::Type) = new{T}()
 end
+
+"""
+  q = RobinBC(left_coefficients, right_coefficients, dx::T, approximation_order) where T # When this BC extends a dimension with a uniform step size
+  q = RobinBC(left_coefficients, right_coefficients, dx::Vector{T}, approximation_order) where T # When this BC extends a dimension with a non uniform step size. dx should be the vector of step sizes for the whole dimension
+-------------------------------------------------------------------------------------
+  The variables in l are [αl, βl, γl], and correspond to a BC of the form αl*u(0) + βl*u'(0) = γl imposed on the lower index boundary.
+  The variables in r are [αl, βl, γl], and correspond to an analagous boundary on the higher index end.
+  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)
+  Write vector b̄₁ as a vertical concatanation with b0 and the rest of the elements of b̄ ₁, denoted b̄`₁, the same with ū into u0 and ū`. b̄`₁ = b̄`_2 = fill(β/Δx, length(stencil)-1)
+  Pull out the product of u0 and b0 from the dot product. The stencil used to approximate u` is denoted s. b0 = α+(β/Δx)*s[1]
+  Rearrange terms to find a general formula for u0:= -b̄`₁̇⋅ū`/b0 + γ/b0, which is dependent on ū` the robin coefficients and Δx.
+  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, V<:AbstractVector{T}} <: AffineBC{T}
     a_l::V
     b_l::T
@@ -57,8 +80,11 @@ struct RobinBC{T, V<:AbstractVector{T}} <: AffineBC{T}
     end
 end
 
+stencil(q::AffineBC{T}, N::Int) where T = ([transpose(q.a_l) transpose(zeros(T, N-length(q.a_l)))], [transpose(zeros(T, N-length(q.a_r))) transpose(q.a_r)])
+affine(q::AffineBC) = (q.b_l, q.b_r)
 
-
+stencil(q::PeriodicBC{T}, N::Int) where T= ([transpose(zeros(T, N-1)) one(T)], [one(T) transpose(zeros(T, N-1))])
+affine(q::PeriodicBC{T}) where T = (zero(T), zero(T))
 """
 q = GeneralBC(α_leftboundary, α_rightboundary, dx::T, approximation_order)