Correct concretization of periodic BCs

Author: jagot

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]
@@ -126,18 +129,22 @@ function SparseArrays.sparse(Q::AffineBC{T}, N::Int) where {T}
     SparseMatrixCSC(Q,N)
 end
 
-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))
+# ** Periodic BCs
+
+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))
+
+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))
+
 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
+
+BandedMatrices.BandedMatrix(::PeriodicBC, ::Int) =
+    throw(ArgumentError("Periodic boundary conditions should be concretized as sparse matrices"))
 
 function LinearAlgebra.Array(Q::BoundaryPaddedVector)
     return [Q.l; Q.u; Q.r]

test/bc_concretizations.jl

@@ -0,0 +1,31 @@
+using SparseArrays, DiffEqOperators, LinearAlgebra, Random,
+    Test, BandedMatrices, FillArrays, LazyArrays, BlockBandedMatrices
+
+@testset "Concretizations of BCs" begin
+    @testset "Periodic BCs" begin
+        N = 9
+        T = Float64
+        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)
+        end
+
+        @testset "Dense concretization" begin
+            Qm,Qu = Array(Q,N)
+
+            @test Qm == correct
+            @test Qm isa Matrix{T}
+            @test Qu == zeros(T,N)
+        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