Skip to content
This repository was archived by the owner on Jul 19, 2023. It is now read-only.

Commit 99e33e2

Browse files
jagotjagot
committed
Correct concretization of affine boundary BCs
1 parent 81a4002 commit 99e33e2

File tree

3 files changed

+181
-30
lines changed

3 files changed

+181
-30
lines changed

‎src/derivative_operators/concretization.jl‎

Lines changed: 48 additions & 23 deletions
Original file line numberDiff line numberDiff line change
@@ -118,49 +118,74 @@ function SparseArrays.SparseMatrixCSC(Q::AffineBC{T}, N::Int) where {T}
118118
end
119119

120120
function BandedMatrices.BandedMatrix(Q::AffineBC{T}, N::Int) where {T}
121-
Q_l = BandedMatrix{T}(Eye(N), (length(Q.a_r)-1, length(Q.a_l)-1))
122-
BandedMatrices.inbands_setindex!(Q_l, Q.a_l, 1, 1:length(Q.a_l))
123-
BandedMatrices.inbands_setindex!(Q_l, Q.a_r, N, (N-length(Q.a_r)+1):N)
121+
# We want the concrete matrix to have as small bandwidths as
122+
# possible, and we accomplish this by dropping all trailing
123+
# zeros. This way, we do not write outside the bands of the
124+
# BandedMatrix.
125+
a_r = Q.a_r[1:something(findlast(!iszero, Q.a_r), 0)]
126+
a_l = Q.a_l[1:something(findlast(!iszero, Q.a_l), 0)]
127+
128+
# Compute bandwidths; the BC matrix should have the shape
129+
#
130+
# [a b c ...
131+
# 1
132+
# 1
133+
# 1
134+
# .
135+
# 1
136+
# ... x y z]
137+
#
138+
# where a,b,c,... and ...,x,y,z are determined by the boundary
139+
# conditions. If these coefficients are zero (Dirichlet0BC), then
140+
# the proper bandwidths are (l,u) = (1,-1).
141+
l = max(count(!iszero, a_r)+1, 1)
142+
u = max(count(!iszero, a_l)-1, -1)
143+
144+
Q_l = BandedMatrix((-1 => ones(T,N),), (N+2,N), (l, u))
145+
for (j,e) ∈ enumerate(a_l)
146+
BandedMatrices.inbands_setindex!(Q_l, e, 1, j)
147+
end
148+
for (j,e) ∈ enumerate(a_r)
149+
BandedMatrices.inbands_setindex!(Q_l, e, N+2, N-length(a_r)+j)
150+
end
124151
Q_b = [Q.b_l; zeros(T,N); Q.b_r]
125-
return (Q_l, Q_b)
152+
153+
Q_l, Q_b
126154
end
127155

128-
function SparseArrays.sparse(Q::AffineBC{T}, N::Int) where {T}
129-
SparseMatrixCSC(Q,N)
130-
end
156+
"""
157+
sparse(Q::AffineBC, N)
158+
159+
Since affine boundary conditions are representable by banded matrices,
160+
that is the default sparse concretization; if you want a
161+
`SparseMatrixCSC`, use `SparseMatrixCSC(Q, N)` instead.
162+
"""
163+
SparseArrays.sparse(Q::AffineBC, N::Int) = BandedMatrix(Q,N)
131164

132165
# ** Periodic BCs
133166

134167
LinearAlgebra.Array(Q::PeriodicBC{T}, N::Int) where T =
135168
([transpose(zeros(T, N-1)) one(T)
136169
Diagonal(ones(T,N))
137-
one(T) transpose(zeros(T, N-1))], zeros(T, N))
170+
one(T) transpose(zeros(T, N-1))],
171+
zeros(T,N+2))
138172

139173
SparseArrays.SparseMatrixCSC(Q::PeriodicBC{T}, N::Int) where T =
140174
(vcat(hcat(zeros(T, 1,N-1), one(T)),
141175
Diagonal(ones(T,N)),
142-
hcat(one(T), zeros(T, 1, N-1))), zeros(T, N))
176+
hcat(one(T), zeros(T, 1, N-1))),
177+
zeros(T,N+2))
143178

144179
SparseArrays.sparse(Q::PeriodicBC{T}, N::Int) where T = SparseMatrixCSC(Q,N)
145180

146181
BandedMatrices.BandedMatrix(::PeriodicBC, ::Int) =
147182
throw(ArgumentError("Periodic boundary conditions should be concretized as sparse matrices"))
148183

149-
function LinearAlgebra.Array(Q::BoundaryPaddedVector)
150-
return [Q.l; Q.u; Q.r]
151-
end
152-
153-
function Base.convert(::Type{Array},A::AbstractBC{T}) where T
154-
Array(A)
155-
end
184+
LinearAlgebra.Array(Q::BoundaryPaddedVector) = [Q.l; Q.u; Q.r]
156185

157-
function Base.convert(::Type{SparseMatrixCSC},A::AbstractBC{T}) where T
158-
SparseMatrixCSC(A)
159-
end
160-
161-
function Base.convert(::Type{AbstractMatrix},A::AbstractBC{T}) where T
162-
SparseMatrixCSC(A)
163-
end
186+
Base.convert(::Type{Array},A::AbstractBC) = Array(A)
187+
Base.convert(::Type{SparseMatrixCSC},A::AbstractBC) = SparseMatrixCSC(A)
188+
Base.convert(::Type{AbstractMatrix},A::AbstractBC) = SparseMatrixCSC(A)
164189

165190
# Multi dimensional BC operators
166191

‎test/bc_coeff_compositions.jl‎

Lines changed: 2 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -93,9 +93,8 @@ end
9393
@test SparseMatrixCSC(A)[2] ≈ (SparseMatrixCSC(L)*SparseMatrixCSC(Q,N)[1], SparseMatrixCSC(L)*SparseMatrixCSC(Q,N)[2])[2]
9494
@test sparse(A)[1] ≈ (sparse(L)*sparse(Q,N)[1], sparse(L)*sparse(Q,N)[2])[1]
9595
@test sparse(A)[2] ≈ (sparse(L)*sparse(Q,N)[1], sparse(L)*sparse(Q,N)[2])[2]
96-
# BandedMatrix not implemeted for boundary operator
97-
@test_broken BandedMatrix(A)[1] ≈ (BandedMatrix(L)*BandedMatrix(Q,N)[1], BandedMatrix(L)*BandedMatrix(Q,N)[2])[1]
98-
@test_broken BandedMatrix(A)[2] ≈ (BandedMatrix(L)*BandedMatrix(Q,N)[1], BandedMatrix(L)*BandedMatrix(Q,N)[2])[2]
96+
@test BandedMatrix(A)[1] ≈ (BandedMatrix(L)*BandedMatrix(Q,N)[1], BandedMatrix(L)*BandedMatrix(Q,N)[2])[1]
97+
@test BandedMatrix(A)[2] ≈ (BandedMatrix(L)*BandedMatrix(Q,N)[1], BandedMatrix(L)*BandedMatrix(Q,N)[2])[2]
9998

10099
# Test that concretization works with multiplication
101100
u = rand(20)

‎test/bc_concretizations.jl‎

Lines changed: 131 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -2,9 +2,132 @@ using SparseArrays, DiffEqOperators, LinearAlgebra, Random,
22
Test, BandedMatrices, FillArrays, LazyArrays, BlockBandedMatrices
33

44
@testset "Concretizations of BCs" begin
5+
T = Float64
6+
L = 10one(T)
7+
N = 9
8+
δx = L/(N+1)
9+
10+
@testset "Affine BCs" begin
11+
@testset "Dirichlet0BC" begin
12+
Q = Dirichlet0BC(T)
13+
14+
correct = vcat(zeros(T,1,N),
15+
Diagonal(ones(T,N)),
16+
zeros(T,1,N))
17+
18+
@testset "$mode concretization" for (mode,Mat,Expected,ExpectedBandwidths) in [
19+
("sparse -> Banded", sparse, BandedMatrix{T}, (1,-1)),
20+
("Banded", BandedMatrix, BandedMatrix{T}, (1,-1)),
21+
("Sparse", SparseMatrixCSC, SparseMatrixCSC{T}, nothing),
22+
("Dense", Array, Matrix{T}, nothing)
23+
]
24+
Qm,Qu = Mat(Q,N)
25+
26+
@test Qm == correct
27+
@test Qm isa Expected
28+
@test Qu == zeros(T,N+2)
29+
30+
!isnothing(ExpectedBandwidths) &&
31+
@test bandwidths(Qm) == ExpectedBandwidths
32+
end
33+
end
34+
35+
@testset "Neumann0BC" begin
36+
Q = Neumann0BC(δx)
37+
38+
correct = vcat(hcat(one(T),zeros(T,1,N-1)),
39+
Diagonal(ones(T,N)),
40+
hcat(zeros(T,1,N-1),one(T)))
41+
42+
@testset "$mode concretization" for (mode,Mat,Expected,ExpectedBandwidths) in [
43+
("sparse -> Banded", sparse, BandedMatrix{T}, (2,0)),
44+
("Banded", BandedMatrix, BandedMatrix{T}, (2,0)),
45+
("Sparse", SparseMatrixCSC, SparseMatrixCSC{T}, nothing),
46+
("Dense", Array, Matrix{T}, nothing)
47+
]
48+
Qm,Qu = Mat(Q,N)
49+
50+
@test Qm == correct
51+
@test Qm isa Expected
52+
@test Qu == zeros(T,N+2)
53+
54+
!isnothing(ExpectedBandwidths) &&
55+
@test bandwidths(Qm) == ExpectedBandwidths
56+
end
57+
58+
@testset "Banded concretization, extra zeros" begin
59+
@testset "lz = $lz" for lz = 0:3
60+
@testset "rz = $rz" for rz = 0:3
61+
Q′ = Neumann0BC(δx)
62+
# Artificially add some zero coefficients, which should
63+
# not increase the bandwidth of the concretized BC.
64+
append!(Q′.a_l, zeros(lz))
65+
append!(Q′.a_r, zeros(rz))
66+
67+
Q′m,Q′u = sparse(Q′,N)
68+
@test bandwidths(Q′m) == (2,0)
69+
70+
@test Q′m == correct
71+
@test Q′u == zeros(T,N+2)
72+
end
73+
end
74+
end
75+
end
76+
77+
@testset "General BCs" begin
78+
@testset "Left BC order = $ld" for ld = 2:5
79+
@testset "Right BC order = $rd" for rd = 2:5
80+
αl = 0.0:ld-1
81+
αr = 0.0:rd-1
82+
83+
Q = GeneralBC(αl, αr, δx)
84+
85+
correct = vcat(hcat(Q.a_l',zeros(T,1,N-(ld-2))),
86+
Diagonal(ones(T,N)),
87+
hcat(zeros(T,1,N-(rd-2)),Q.a_r'))
88+
89+
Qm,Qu = sparse(Q,N)
90+
91+
@test Qm == correct
92+
@test Qm isa BandedMatrix{T}
93+
@test bandwidths(Qm) == (rd-1,ld-3)
94+
95+
@test Qu == vcat(Q.b_l,zeros(T,N),Q.b_r)
96+
end
97+
end
98+
end
99+
100+
@testset "Dirichlet0BC" begin
101+
# This is equivalent to a Dirichlet0BC; the trailing zeros
102+
# should be dropped and the bandwidths optimal.
103+
Q = GeneralBC([0.0, 1.0, 0.0, 0.0], [0.0, 1.0, 0.0, 0.0, 0.0], δx)
104+
105+
correct = vcat(zeros(T,1,N),
106+
Diagonal(ones(T,N)),
107+
zeros(T,1,N))
108+
109+
Qm = first(sparse(Q,N))
110+
@test Qm == correct
111+
@test bandwidths(Qm) == (1,-1)
112+
end
113+
114+
@testset "Almost DirichletBC" begin
115+
Q = GeneralBC([1.0, 1.0, 0.0, 0.0, eps(Float64)],
116+
[1.0, 1.0, 0.0, 0.0, 0.0], δx)
117+
118+
correct = vcat(zeros(T,1,N),
119+
Diagonal(ones(T,N)),
120+
zeros(T,1,N))
121+
122+
Qm,Qu = sparse(Q,N)
123+
124+
@test Qm ≈ correct
125+
@test bandwidths(Qm) == (1,2)
126+
@test Qu ≈ vcat(-one(T),zeros(T,N),-one(T))
127+
end
128+
end
129+
5130
@testset "Periodic BCs" begin
6-
N = 9
7-
T = Float64
8131
Q = PeriodicBC(T)
9132
@test_throws ArgumentError BandedMatrix(Q,N)
10133

@@ -17,15 +140,19 @@ using SparseArrays, DiffEqOperators, LinearAlgebra, Random,
17140

18141
@test Qm == correct
19142
@test Qm isa SparseMatrixCSC{T}
20-
@test Qu == zeros(T,N)
143+
@test Qu == zeros(T,N+2)
144+
145+
Qm′ = first(sparse(Q, N))
146+
@test Qm′ == correct
147+
@test Qm′ isa SparseMatrixCSC{T}
21148
end
22149

23150
@testset "Dense concretization" begin
24151
Qm,Qu = Array(Q,N)
25152

26153
@test Qm == correct
27154
@test Qm isa Matrix{T}
28-
@test Qu == zeros(T,N)
155+
@test Qu == zeros(T,N+2)
29156
end
30157
end
31158
end

0 commit comments

Comments
 (0)