[WIP] Correct concretization BCs and ghost derivative operators

Project.toml

@@ -22,6 +22,7 @@ ModelingToolkit = "0.8.0"
 julia = "1"
 
 [extras]
+Compat = "34da2185-b29b-5c13-b0c7-acf172513d20"
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
@@ -30,4 +31,4 @@ SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["FillArrays", "OrdinaryDiffEq", "Random", "SafeTestsets", "SpecialFunctions", "Test"]
+test = ["Compat", "FillArrays", "OrdinaryDiffEq", "Random", "SafeTestsets", "SpecialFunctions", "Test"]

src/derivative_operators/BC_operators.jl

@@ -19,7 +19,7 @@ Qx, Qy, ... = PeriodicBC{T}(size(u)) #When all dimensions are to be extended wit
 
 -------------------------------------------------------------------------------------
 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.
+It is not recommended to concretize this BC type in to a BandedMatrix, since the vast majority of bands will be all 0s. SparseMatrix concretization is recommended.
 """
 struct PeriodicBC{T} <: AtomicBC{T}
     PeriodicBC(T::Type) = new{T}()

src/derivative_operators/concretization.jl

@@ -101,7 +101,10 @@ end
 # Boundary Condition Operator concretizations
 ################################################################################
 
-#Atomic BCs
+# * Atomic BCs
+
+# ** Affine BCs
+
 function LinearAlgebra.Array(Q::AffineBC{T}, N::Int) where {T}
     Q_L = [transpose(Q.a_l) transpose(zeros(T, N-length(Q.a_l))); Diagonal(ones(T,N)); transpose(zeros(T, N-length(Q.a_r))) transpose(Q.a_r)]
     Q_b = [Q.b_l; zeros(T,N); Q.b_r]
@@ -115,45 +118,74 @@ function SparseArrays.SparseMatrixCSC(Q::AffineBC{T}, N::Int) where {T}
 end
 
 function BandedMatrices.BandedMatrix(Q::AffineBC{T}, N::Int) where {T}
-    Q_l = BandedMatrix{T}(Eye(N), (length(Q.a_r)-1, length(Q.a_l)-1))
-    BandedMatrices.inbands_setindex!(Q_l, Q.a_l, 1, 1:length(Q.a_l))
-    BandedMatrices.inbands_setindex!(Q_l, Q.a_r, N, (N-length(Q.a_r)+1):N)
+    # We want the concrete matrix to have as small bandwidths as
+    # possible, and we accomplish this by dropping all trailing
+    # zeros. This way, we do not write outside the bands of the
+    # BandedMatrix.
+    a_r = Q.a_r[1:something(findlast(!iszero, Q.a_r), 0)]
+    a_l = Q.a_l[1:something(findlast(!iszero, Q.a_l), 0)]
+
+    # Compute bandwidths; the BC matrix should have the shape
+    #
+    # [a b c ...
+    #  1
+    #    1
+    #      1
+    #        .
+    #          1
+    #  ... x y z]
+    #
+    # where a,b,c,... and ...,x,y,z are determined by the boundary
+    # conditions. If these coefficients are zero (Dirichlet0BC), then
+    # the proper bandwidths are (l,u) = (1,-1).    
+    l = max(count(!iszero, a_r)+1, 1)
+    u = max(count(!iszero, a_l)-1, -1)
+
+    Q_l = BandedMatrix((-1 => ones(T,N),), (N+2,N), (l, u))
+    for (j,e) ∈ enumerate(a_l)
+        BandedMatrices.inbands_setindex!(Q_l, e, 1, j)
+    end
+    for (j,e) ∈ enumerate(a_r)
+        BandedMatrices.inbands_setindex!(Q_l, e, N+2, N-length(a_r)+j)
+    end
     Q_b = [Q.b_l; zeros(T,N); Q.b_r]
-    return (Q_l, Q_b)
+
+    Q_l, Q_b
 end
 
-function SparseArrays.sparse(Q::AffineBC{T}, N::Int) where {T}
-    SparseMatrixCSC(Q,N)
-end
+"""
+    sparse(Q::AffineBC, N)
 
-LinearAlgebra.Array(Q::PeriodicBC{T}, N::Int) where T = (Array([transpose(zeros(T, N-1)) one(T); Diagonal(ones(T,N)); one(T) transpose(zeros(T, N-1))]), zeros(T, N))
-SparseArrays.SparseMatrixCSC(Q::PeriodicBC{T}, N::Int) where T = ([transpose(zeros(T, N-1)) one(T); Diagonal(ones(T,N)); one(T) transpose(zeros(T, N-1))], zeros(T, N))
-SparseArrays.sparse(Q::PeriodicBC{T}, N::Int) where T = SparseMatrixCSC(Q,N)
-function BandedMatrices.BandedMatrix(Q::PeriodicBC{T}, N::Int) where T #Not reccomended!
-    Q_array = BandedMatrix{T}(Eye(N), (N-1, N-1))
-    Q_array[1, end] = one(T)
-    Q_array[1, 1] = zero(T)
-    Q_array[end, 1] = one(T)
-    Q_array[end, end] = zero(T)
-
-    return (Q_array, zeros(T, N))
-end
+Since affine boundary conditions are representable by banded matrices,
+that is the default sparse concretization; if you want a
+`SparseMatrixCSC`, use `SparseMatrixCSC(Q, N)` instead.
+"""
+SparseArrays.sparse(Q::AffineBC, N::Int) = BandedMatrix(Q,N)
 
-function LinearAlgebra.Array(Q::BoundaryPaddedVector)
-    return [Q.l; Q.u; Q.r]
-end
+# ** Periodic BCs
 
-function Base.convert(::Type{Array},A::AbstractBC{T}) where T
-    Array(A)
-end
+LinearAlgebra.Array(Q::PeriodicBC{T}, N::Int) where T =
+    ([transpose(zeros(T, N-1)) one(T)
+      Diagonal(ones(T,N))
+      one(T) transpose(zeros(T, N-1))],
+     zeros(T,N+2))
 
-function Base.convert(::Type{SparseMatrixCSC},A::AbstractBC{T}) where T
-    SparseMatrixCSC(A)
-end
+SparseArrays.SparseMatrixCSC(Q::PeriodicBC{T}, N::Int) where T =
+    (vcat(hcat(zeros(T, 1,N-1), one(T)),
+          Diagonal(ones(T,N)),
+          hcat(one(T), zeros(T, 1, N-1))),
+     zeros(T,N+2))
 
-function Base.convert(::Type{AbstractMatrix},A::AbstractBC{T}) where T
-    SparseMatrixCSC(A)
-end
+SparseArrays.sparse(Q::PeriodicBC{T}, N::Int) where T = SparseMatrixCSC(Q,N)
+
+BandedMatrices.BandedMatrix(::PeriodicBC, ::Int) =
+    throw(ArgumentError("Periodic boundary conditions should be concretized as sparse matrices"))
+
+LinearAlgebra.Array(Q::BoundaryPaddedVector) = [Q.l; Q.u; Q.r]
+
+Base.convert(::Type{Array},A::AbstractBC) = Array(A)
+Base.convert(::Type{SparseMatrixCSC},A::AbstractBC) = SparseMatrixCSC(A)
+Base.convert(::Type{AbstractMatrix},A::AbstractBC) = SparseMatrixCSC(A)
 
 # Multi dimensional BC operators
 

test/bc_coeff_compositions.jl

@@ -86,21 +86,29 @@ end
     u = rand(22)
     @test (L + L2) * u ≈ convert(AbstractMatrix,L + L2) * u ≈ (BandedMatrix(L) + BandedMatrix(L2)) * u
 
-    # Test concretization
-    @test Array(A)[1] ≈ (Array(L)*Array(Q,N)[1], Array(L)*Array(Q,N)[2])[1]
-    @test Array(A)[2] ≈ (Array(L)*Array(Q,N)[1], Array(L)*Array(Q,N)[2])[2]
-    @test SparseMatrixCSC(A)[1] ≈ (SparseMatrixCSC(L)*SparseMatrixCSC(Q,N)[1], SparseMatrixCSC(L)*SparseMatrixCSC(Q,N)[2])[1]
-    @test SparseMatrixCSC(A)[2] ≈ (SparseMatrixCSC(L)*SparseMatrixCSC(Q,N)[1], SparseMatrixCSC(L)*SparseMatrixCSC(Q,N)[2])[2]
-    @test sparse(A)[1] ≈ (sparse(L)*sparse(Q,N)[1], sparse(L)*sparse(Q,N)[2])[1]
-    @test sparse(A)[2] ≈ (sparse(L)*sparse(Q,N)[1], sparse(L)*sparse(Q,N)[2])[2]
-    # BandedMatrix not implemeted for boundary operator
-    @test_broken BandedMatrix(A)[1] ≈ (BandedMatrix(L)*BandedMatrix(Q,N)[1], BandedMatrix(L)*BandedMatrix(Q,N)[2])[1]
-    @test_broken BandedMatrix(A)[2] ≈ (BandedMatrix(L)*BandedMatrix(Q,N)[1], BandedMatrix(L)*BandedMatrix(Q,N)[2])[2]
-
-    # Test that concretization works with multiplication
-    u = rand(20)
-    @test Array(A)[1]*u + Array(A)[2] ≈ L*(Q*u) ≈ A*u
-    @test sparse(A)[1]*u + sparse(A)[2] ≈ L*(Q*u) ≈ A*u
+    @testset "$mode concretization" for (mode,Mat) in [("Dense", Array),
+                                                       ("Sparse", SparseMatrixCSC),
+                                                       ("Best sparse", sparse),
+                                                       ("BandedMatrix", BandedMatrix)]
+        Am,Au = Mat(A)
+        Lm = Mat(L)
+        Qm,Qu = Mat(Q,N)
+
+        @test Am ≈ Lm*Qm
+        @test Au ≈ Lm*Qu
+    end
+
+    @testset "$mode concrete multiplication" for (mode,Mat) in [("Dense", Array),
+                                                                ("Sparse", SparseMatrixCSC),
+                                                                ("Best sparse", sparse),
+                                                                ("BandedMatrix", BandedMatrix)]
+        u = rand(20)
+        Am,Au = Mat(A)
+        Lm = Mat(L)
+        Qm,Qu = Mat(Q,N)
+
+        @test Am*u + Au ≈ L*(Q*u) ≈ A*u
+    end
 end
 
 @testset "Test Left Division L2 (second order)" begin

test/bc_concretizations.jl

@@ -0,0 +1,158 @@
+using SparseArrays, DiffEqOperators, LinearAlgebra, Random,
+    Test, BandedMatrices, FillArrays, LazyArrays, BlockBandedMatrices, Compat
+
+@testset "Concretizations of BCs" begin
+    T = Float64
+    L = 10one(T)
+    N = 9
+    δx = L/(N+1)
+
+    @testset "Affine BCs" begin
+        @testset "Dirichlet0BC" begin
+            Q = Dirichlet0BC(T)
+
+            correct = vcat(zeros(T,1,N),
+                           Diagonal(ones(T,N)),
+                           zeros(T,1,N))
+
+            @testset "$mode concretization" for (mode,Mat,Expected,ExpectedBandwidths) in [
+                ("sparse -> Banded", sparse, BandedMatrix{T}, (1,-1)),
+                ("Banded", BandedMatrix, BandedMatrix{T}, (1,-1)),
+                ("Sparse", SparseMatrixCSC, SparseMatrixCSC{T}, nothing),
+                ("Dense", Array, Matrix{T}, nothing)
+            ]
+                Qm,Qu = Mat(Q,N)
+
+                @test Qm == correct
+                @test Qm isa Expected
+                @test Qu == zeros(T,N+2)
+
+                !isnothing(ExpectedBandwidths) &&
+                    @test bandwidths(Qm) == ExpectedBandwidths
+            end
+        end
+
+        @testset "Neumann0BC" begin
+            Q = Neumann0BC(δx)
+
+            correct = vcat(hcat(one(T),zeros(T,1,N-1)),
+                           Diagonal(ones(T,N)),
+                           hcat(zeros(T,1,N-1),one(T)))
+
+            @testset "$mode concretization" for (mode,Mat,Expected,ExpectedBandwidths) in [
+                ("sparse -> Banded", sparse, BandedMatrix{T}, (2,0)),
+                ("Banded", BandedMatrix, BandedMatrix{T}, (2,0)),
+                ("Sparse", SparseMatrixCSC, SparseMatrixCSC{T}, nothing),
+                ("Dense", Array, Matrix{T}, nothing)
+            ]
+                Qm,Qu = Mat(Q,N)
+
+                @test Qm == correct
+                @test Qm isa Expected
+                @test Qu == zeros(T,N+2)
+
+                !isnothing(ExpectedBandwidths) &&
+                    @test bandwidths(Qm) == ExpectedBandwidths
+            end
+
+            @testset "Banded concretization, extra zeros" begin
+                @testset "lz = $lz" for lz = 0:3
+                    @testset "rz = $rz" for rz = 0:3
+                        Q′ = Neumann0BC(δx)
+                        # Artificially add some zero coefficients, which should
+                        # not increase the bandwidth of the concretized BC.
+                        append!(Q′.a_l, zeros(lz))
+                        append!(Q′.a_r, zeros(rz))
+
+                        Q′m,Q′u = sparse(Q′,N)
+                        @test bandwidths(Q′m) == (2,0)
+
+                        @test Q′m == correct
+                        @test Q′u == zeros(T,N+2)
+                    end
+                end
+            end
+        end
+
+        @testset "General BCs" begin
+            @testset "Left BC order = $ld" for ld = 2:5
+                @testset "Right BC order = $rd" for rd = 2:5
+                    αl = 0.0:ld-1
+                    αr = 0.0:rd-1
+
+                    Q = GeneralBC(αl, αr, δx)
+
+                    correct = vcat(hcat(Q.a_l',zeros(T,1,N-(ld-2))),
+                                   Diagonal(ones(T,N)),
+                                   hcat(zeros(T,1,N-(rd-2)),Q.a_r'))
+
+                    Qm,Qu = sparse(Q,N)
+
+                    @test Qm == correct
+                    @test Qm isa BandedMatrix{T}
+                    @test bandwidths(Qm) == (rd-1,ld-3)
+
+                    @test Qu == vcat(Q.b_l,zeros(T,N),Q.b_r)
+                end
+            end
+        end
+
+        @testset "Dirichlet0BC" begin
+            # This is equivalent to a Dirichlet0BC; the trailing zeros
+            # should be dropped and the bandwidths optimal.
+            Q = GeneralBC([0.0, 1.0, 0.0, 0.0], [0.0, 1.0, 0.0, 0.0, 0.0], δx)
+
+            correct = vcat(zeros(T,1,N),
+                           Diagonal(ones(T,N)),
+                           zeros(T,1,N))
+
+            Qm = first(sparse(Q,N))
+            @test Qm == correct
+            @test bandwidths(Qm) == (1,-1)
+        end
+
+        @testset "Almost DirichletBC" begin
+            Q = GeneralBC([1.0, 1.0, 0.0, 0.0, eps(Float64)],
+                          [1.0, 1.0, 0.0, 0.0, 0.0], δx)
+
+            correct = vcat(zeros(T,1,N),
+                           Diagonal(ones(T,N)),
+                           zeros(T,1,N))
+
+            Qm,Qu = sparse(Q,N)
+
+            @test Qm ≈ correct
+            @test bandwidths(Qm) == (1,2)
+            @test Qu ≈ vcat(-one(T),zeros(T,N),-one(T))
+        end
+    end
+
+    @testset "Periodic BCs" begin
+        Q = PeriodicBC(T)
+        @test_throws ArgumentError BandedMatrix(Q,N)
+
+        correct = vcat(hcat(zeros(T,1,N-1),one(T)),
+                       Diagonal(ones(T,N)),
+                       hcat(one(T),zeros(T,1,N-1)))
+
+        @testset "Sparse concretization" begin
+            Qm,Qu = SparseMatrixCSC(Q,N)
+
+            @test Qm == correct
+            @test Qm isa SparseMatrixCSC{T}
+            @test Qu == zeros(T,N+2)
+
+            Qm′ = first(sparse(Q, N))
+            @test Qm′ == correct
+            @test Qm′ isa SparseMatrixCSC{T}
+        end
+
+        @testset "Dense concretization" begin
+            Qm,Qu = Array(Q,N)
+
+            @test Qm == correct
+            @test Qm isa Matrix{T}
+            @test Qu == zeros(T,N+2)
+        end
+    end
+end

test/runtests.jl

@@ -19,3 +19,4 @@ import Base: isapprox
 @time @safetestset "Higher Dimensional Concretization" begin include("concretization.jl") end
 @time @safetestset "Upwind Operator Interface" begin include("upwind_operators_interface.jl") end
 @time @safetestset "MOLFiniteDifference Interface" begin include("MOLtest.jl") end
+@time @safetestset "BC concretizations" begin include("bc_concretizations.jl") end