This repository was archived by the owner on Jul 19, 2023. It is now read-only.
This repository was archived by the owner on Jul 19, 2023. It is now read-only.
Banded concretization of derivatives with boundary conditions gives dense matrices #182
Description
Consider the following example:
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> BandedMatrix(Δ)
9×11 BandedMatrix{Float64,Array{Float64,2},Base.OneTo{Int64}}:
1.0 -2.0 1.0 0.0 0.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
0.0 1.0 -2.0 1.0 0.0 0.0 ⋅ ⋅ ⋅ ⋅ ⋅
0.0 0.0 1.0 -2.0 1.0 0.0 0.0 ⋅ ⋅ ⋅ ⋅
0.0 0.0 0.0 1.0 -2.0 1.0 0.0 0.0 ⋅ ⋅ ⋅
0.0 0.0 0.0 0.0 1.0 -2.0 1.0 0.0 0.0 ⋅ ⋅
⋅ 0.0 0.0 0.0 0.0 1.0 -2.0 1.0 0.0 0.0 ⋅
⋅ ⋅ 0.0 0.0 0.0 0.0 1.0 -2.0 1.0 0.0 0.0
⋅ ⋅ ⋅ 0.0 0.0 0.0 0.0 1.0 -2.0 1.0 0.0
⋅ ⋅ ⋅ ⋅ 0.0 0.0 0.0 0.0 1.0 -2.0 1.0(Should not the bandwidths be smaller?)
If I now add BCs, I get a dense matrix upon concretization:
julia> Q = Dirichlet0BC(typeof(δx))
RobinBC{Float64,Array{Float64,1}}([-0.0, 0.0], 0.0, [-0.0, 0.0], 0.0)
julia> M,u = BandedMatrix(Δ*Q);
julia> M
9×9 Array{Float64,2}:
-2.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
1.0 -2.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 1.0 -2.0 1.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 1.0 -2.0 1.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 1.0 -2.0 1.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 1.0 -2.0 1.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 1.0 -2.0 1.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 1.0 -2.0 1.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 -2.0This is because the BCs are concretized as dense matrices. Arguably, they should be BandedMatrices as well, however that function does not work either:
inbands_setindex!does not support ranges anymore, it seems,- some typos (
Q_LvsQ_l)
-
BandedMatrixof a boundary should return a banded matrix. It is fine ifPeriodicBCthrows an error here. -
sparseof a BC should return aBandedMatrixunless it's periodic, in which case it should beSparseMatrixCSC -
sparseon aGhostDerivativeOperatorshould preferBandedMatrixunlessPeriodicBC - The bandwidths are too large