Split Q Multidimensional BCs [Frozen, pending review]

[xtalax]

Similar to the other PR, but now the MultiDimBC operators only act along a specified axis, stored in the type information. Likewise Q*u now returns a BoundaryPaddedArray that extends only along the axis of Q. There is also support for composed Q, that extends on all dimensions at once, and a ComposedBoundaryPaddedArray.
There is now an additional constructor which makes it easy to extend an atomic BC to act on a whole dimension, as well as a constructor which will generate all operators needed to extend an atomic BC to the whole domain.

Atomic Boundary Conditions

An atomic boundary condition can be thought of as acting on a single vector, applying boundary conditions at each end.

The currently implemented AtomicBCs are as follows:

• RobinBC: Implements a robin boundary condition, to any approximation order for the derivative.
  Usage: RobinBC([αl, βl, γl], [αr, βr, γr], dx, approximation_order::Int) corresponding to a BC of the form αl*u[0] + βl*u'[0] = γl, αr*u[end+1] + βr*u'[end+1] = γr. For a uniform grid, dx should be a value, for nonuniform it should be the vector of step sizes for x.

• DirichletBC and Dirichlet0BC: Implemented as specialised constructors to RobinBC.
  Usage: DirichletBC(αl, αr) corresponding to a BC of the form u[0] = αl, u[end+1] = αr. Dirichlet0BC(T::Type) gives a bc where αl and αr are eual to 0. For a uniform grid, dx should be a value, for nonuniform it should be the vector of step sizes for x.

• NeumannBC and Neumann0BC: Implemented as specialised constructors to RobinBC.
  Usage: NeumannBC([αl, αr], dx, approximation_order) corresponding to a BC of the form u'[0] = αl, u'[end+1] = αr. Neumann0BC(dx, approximation_order) gives a bc where αl and αr are eual to 0. For a uniform grid, dx should be a value, for nonuniform it should be the vector of step sizes for x.

• GeneralBC: A generalization of the RobinBC, including derivative terms up to any order.
  Usage: GeneralBC(αl::Vector{T}, αr::Vector{T}, dx, approximation_order) Represents a condition of the form αl[1] + αl[2]u[0] + αl[3]u'[0] + αl[4]u''[0]+... = 0, for the left boundary, and αr[1] + αr[2]u[end] + αr[3]u'[end] + αr[4]u''[end]+... = 0 for the right boundary. αl and αr can be of any length greater than 1, and don't have to be the same length in general. For a uniform grid, dx should be a value, for nonuniform it should be the vector of step sizes for x.

• PeriodicBC: Wraps the domain around such that u[0] = u[end] and u[end+1] = u[1].
  Usage: PeriodicBC{T}()

Boundary Conditions on Multi-Dimensional Domains

A MultiDimBC applies a boundary condition along the entire boundary orthoganal to a dimension of some discretized domain.
As an example: If you have an array a with size(a) = (n, m, o), then a MultiDimBC that acts on the first dimension of a would create an extension of size(m,o) at both the lower and upper index end of the 1st dimension of a.
There are some constructors that make it easy to apply the same AtomicBC along the entire boundary of a domain:

dx = 0.01
n = 100
m = 120
o = 78

x = range(0.0; step=dx, length=n)
y = range(0.0; step=dx, length=m)
z = range(0.0; step=dx, length=o)

a = rand(n,m, o)

q3 = Dirichlet0BC(Float64) #Create an atomic bc
#Create a BC operator that applies the atomic BC to the whole boundary orthoganal to dim 3
Qz = MultiDimBC(q3, size(a), 3) 

#Extend A along dim 3, creating a BoundaryPaddedArray extended along this dimension.
az = Qz*a

You may want to apply different boundary conditions on different parts of a boundary. That can be achieved by building an Array{AtomicBC{T}, N-1} (where a isa Array{T, N}) of Atomic BCs with size = size(a)[setdiff(1:N, dim)], and passing that to MultiDimBC like so:

#Create atomic BCs
q1 = RobinBC([1.0, 2.0, 3.0], [0.0, -1.0, 2.0], (0.1, 0.1), 4.0)
q2 = PeriodicBC{Float64}()

#Fill an Array so that each BC is applied on half of the boundary of dim 1
BCx = vcat(fill(q1, (div(m,2), o)), fill(q2, (div(m,2), o)))  

#Create a BC that is dependent on x and z with a comprehension
BCy = [NeumannBC([X*Z, X+Z], dx, 3) for X in x, Z in z] 

Qx = MultiDimBC(BCx, 1) #create an operator that extends on dim 1
Qy = MultiDimBC(BCy, 2) #create an operator that extends on dim 2

ax = Qx*a #BoundaryPaddedArray extended on dim 1 
ay = Qy*a #BoundaryPaddedArray extended on dim 2

There is also an easy constructor that extends the same atomic BC to cover all boundaries at once (a must be defined on a uniform grid if you want to use a Neumann, Neumann0, Robin or General boundary condition, these have special constructors for the non uniform grid case):

q = Dirichlet0BC(Float64)
Qx, Qy, Qz = MultiDimBC(q, size(a)) #leave off the dim argument to extend all dimensions
ax = Qx*a #BoundaryPaddedArray extended on dim 1 
ay = Qy*a #BoundaryPaddedArray extended on dim 2
az = Qz*a #BoundaryPaddedArray extended on dim 3

On a non uniform grid, for the above mentioned boundary types, you can cover the whole boundary on all dimensions like this:

Qx, Qy, Qz = RobinBC(l, r, (dx::Vector, dy::Vector, dz::Vector), approximation_order, size(a))

There are analogously defined constructors for the other above mentioned BCs, simply replace the usual dx with a tuple containing vectors of the step sizes for each dimension, and append size(a) to the argument signature.

Another constructor for PeriodicBC:

Qx, Qy, Qz = PeriodicBC{Float64}(size(a))

ComposedBCOperators/BoundaryPaddedArrays

It is also possible to create an operator that extends all boundaries at once - a ComposedMultiDimBC. This supports syntax like (Lxx + Lyy + Lzz)Q as opposed to the split representation which can be used like LxxQx + LyyQy +LzzQz.

One goes between these two representations like this:

Q = compose(Qx,Qy,Qz)
Qx, Qy, Qz = decompose(Q)

Also works with BoundaryPaddedArrays, if Ax, Ay, and Az are padded arrays extended in each of the dimensions respectively:

A = compose(Ax,Ay,Az) #This A is a ComposedBoundaryPaddedArray
Ax, Ay, Az = decompose(A)

compose() can accept padded arrays / bc operators in any order, as long as there is exactly 1 array/operator that extends on each dimension. An array/operator must be supplied for each dimension in order for compose to work properly.

The 2 representations are equivalent, but the composed representation is faster for some types of solve [which?].

julia> compose(Qx,Qy,Qz)*a == compose(Qx*a, Qy*a, Qz*a)
true

getindex works as expected with BoundaryPaddedArrays and ComposedBoundaryPaddedArrays, for any number of dimensions. Colon indexing is supported for the BoundaryPaddedArray, but only on the extended dimension OR the other dimensions. It is not supported for the ComposedBoundaryPaddedArray.
It is recomended when writing fast algorithms that operate on (Composed)BoundaryPaddedArrays to index in to the fields of the padded array directly, rather than relying on the getindex function for this type as it is quite slow.

Here is the BoundaryPaddedArray for reference (D is axis of extension):

struct BoundaryPaddedArray{T<:Number, D, N, M, V<:AbstractArray{T, N}, B<: AbstractArray{T, M}} <: AbstractBoundaryPaddedArray{T,D,N}
    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
end

And the ComposedBoundaryPaddedArray:

struct ComposedBoundaryPaddedArray{T<:Number, N, M, V<:AbstractArray{T, N}, B<: AbstractArray{T, M}} <: AbstractBoundaryPaddedArray{T, N}
    #A vector of dimension M=N-1 arrays that extend the lower index ends of each boundary
    lower::Vector{B} 
    #The same case for the upper index ends.
    upper::Vector{B}
    u::V
end

Concretizations are now implemented for the multi dimensional BC operators.
Per the suggestion from @ajozefiak, the dimensions of the affine parts are permuted so that they point along the dimension of extension for the operator.

Array(Q::SingleLayerMultiDimBC) - the split Q type, returns a tuple, the first element of which is an array of the shape of the boundary, filled with the linear operator parts of the respective Atomic BCs, the second element is an array of affine parts.

Array(Q::ComposedMultiDimBC) returns a tuple with length equal to the number of dimensions, each element of which is a tuple as described in the paragraph above, corresponding to the BCs for each dimension.

There should be a way to unpack the array of Q_Ls in to some huge block banded matrix, I'm having trouble working out how though.

Sparse concretization is also implemented with methods of sparse() and SparseMatrixCSC()

BandedMatrix() concretization is now implemented for all BCs.
This is not recommended for PeriodicBC though since its matrix has entries in the upper right and lower left corners - This could be moved to the affine part, but would include the difficulties of BridgeBC (No longer in this PR)

There is a utility function perpsize(a::AbstractArray, dim::Int) which returns the size of a perpendicular to dimension dim

[ChrisRackauckas]

leave this to AbstractVector, as an optimization down the line may be to use StaticVector here.

[ChrisRackauckas]

No reason to allocate the vectors here. Should be a tuple. Let's mark this as a small improvement to make in a follow up PR.

[ChrisRackauckas]

oh I see why there's a compose function

[ChrisRackauckas]

dim may need to be a type parameter in order for it to infer well for dispatch

[xtalax]

Does something like:

foo(x, T) where T = bar{T}(x)

Work in julia? I thought that the parameters had to be inferred from argument types.

[ChrisRackauckas]

`foo{T}(x) where T = ...

[ChrisRackauckas]

There are still things that I think should get changed, but this is very good and I don't want it to go stale, so let's merge it in and open an issue to change the last few details. Very nicely done!