[WIP] Correct concretization BCs and ghost derivative operators

[jagot]

Fixes #182.

So far, the different concretizations of boundary conditions have been implemented, ghost derivative operators are up next.

Please check u for PeriodicBCs, I was unsure if they should be padded with zero(T) or one(T); they were of incorrect length N, need to be N+2.

[ChrisRackauckas]

IMO, this should act like sparse and give "the best choice"

[jagot]

Well, none of these convert methods actually work at all, since the matrix constructors require a length argument N to concretize a boundary condition, and there is no possible default length.

[ChrisRackauckas]

Then maybe it's best to just delete them?

[jagot]

Will do.

[ChrisRackauckas]

Suggested change

	Base.convert(::Type{AbstractMatrix},A::AbstractBC) = SparseMatrixCSC(A)
	Base.convert(::Type{AbstractMatrix},A::AbstractBC) = sparse(A)

[xtalax]

@jagot

Just a heads up, The GhostDerivativeOperators are going to change quite a lot once #170 is merged, as they currently assume that the domain is one dimensional - this PR sets them up to work in multiple dimensions.
The concretizations for these are already implemented there - just needing a concretization for the multiple dimensional boundary conditions to banded.

[jagot]

@xtalax Thanks for the heads-up! As you see in the latest commit to this PR, I did not do anything big, just cleaned the concretization code.

However, I did not yet figure out how to shrink the bandwidths properly. You can remove bstl from this line, which will give the output below, but will break some of the tests, so some more thinking on my part is to be done.

julia> L = 10.0
10.0

julia> N = 9
9

julia> δx = L/(N+1)
1.0

julia> Δ = CenteredDifference(2, 2, δx, N)
DerivativeOperator{Float64,1,false,Float64,StaticArrays.SArray{Tuple{3},Float64,1,3},StaticArrays.SArray{Tuple{0},StaticArrays.SArray{Tuple{4},Float64,1,4},1,0},Nothing,Nothing}(2, 2, 1.0, 9, 3, [1.0, -2.0, 1.0], 4, 0, StaticArrays.SArray{Tuple{4},Float64,1,4}[], StaticArrays.SArray{Tuple{4},Float64,1,4}[], nothing, nothing)

julia> Q = Dirichlet0BC(typeof(δx))
RobinBC{Float64,Array{Float64,1}}([-0.0, 0.0], 0.0, [-0.0, 0.0], 0.0)

julia> A = Δ*Q;

julia> first(sparse(A))
9×9 BandedMatrix{Float64,Array{Float64,2},Base.OneTo{Int64}}:
 -2.0   1.0    ⋅     ⋅     ⋅     ⋅     ⋅     ⋅     ⋅
  1.0  -2.0   1.0    ⋅     ⋅     ⋅     ⋅     ⋅     ⋅
   ⋅    1.0  -2.0   1.0    ⋅     ⋅     ⋅     ⋅     ⋅
   ⋅     ⋅    1.0  -2.0   1.0    ⋅     ⋅     ⋅     ⋅
   ⋅     ⋅     ⋅    1.0  -2.0   1.0    ⋅     ⋅     ⋅
   ⋅     ⋅     ⋅     ⋅    1.0  -2.0   1.0    ⋅     ⋅
   ⋅     ⋅     ⋅     ⋅     ⋅    1.0  -2.0   1.0    ⋅
   ⋅     ⋅     ⋅     ⋅     ⋅     ⋅    1.0  -2.0   1.0
   ⋅     ⋅     ⋅     ⋅     ⋅     ⋅     ⋅    1.0  -2.0

However, Neumann BCs (et al.) do not generate optimal bandwidths, since banded–banded multiplication has no way of knowing that most of the sub-/superpdiagonals are actually zero:

julia> Q = Neumann0BC(δx)
RobinBC{Float64,Array{Float64,1}}([1.0], -0.0, [1.0], 0.0)

julia> A = Δ*Q;

julia> first(sparse(A))
9×9 BandedMatrix{Float64,Array{Float64,2},Base.OneTo{Int64}}:
 -1.0   1.0   0.0    ⋅     ⋅     ⋅     ⋅     ⋅     ⋅
  1.0  -2.0   1.0   0.0    ⋅     ⋅     ⋅     ⋅     ⋅
  0.0   1.0  -2.0   1.0   0.0    ⋅     ⋅     ⋅     ⋅
   ⋅    0.0   1.0  -2.0   1.0   0.0    ⋅     ⋅     ⋅
   ⋅     ⋅    0.0   1.0  -2.0   1.0   0.0    ⋅     ⋅
   ⋅     ⋅     ⋅    0.0   1.0  -2.0   1.0   0.0    ⋅
   ⋅     ⋅     ⋅     ⋅    0.0   1.0  -2.0   1.0   0.0
   ⋅     ⋅     ⋅     ⋅     ⋅    0.0   1.0  -2.0   1.0
   ⋅     ⋅     ⋅     ⋅     ⋅     ⋅    0.0   1.0  -1.0

No tests for this yet.

[xtalax]

Why isn't this something(findfirst(!iszero, Q.a_r):end ?

[jagot]

Don't know :) I have not thought about this PR since September.

[xtalax]

I just noticed hat I had started a review and not submitted it, probably also since September - better late than never :P