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

Commit 83a9b6c

Browse files
AndiMDAndiMD
committed
Add distinct type for grid vs. function value
Test Robin BC, Dirichlet, Neumann with complex function values
1 parent 76f0105 commit 83a9b6c

File tree

2 files changed

+136
-21
lines changed

2 files changed

+136
-21
lines changed

Diff for: ‎src/derivative_operators/BC_operators.jl

+24-21
Original file line numberDiff line numberDiff line change
@@ -22,28 +22,28 @@ struct RobinBC{T, V<:AbstractVector{T}} <: AffineBC{T}
2222
b_l::T
2323
a_r::V
2424
b_r::T
25-
function RobinBC(l::NTuple{3,T}, r::NTuple{3,T}, dx::T, order = 1) where {T}
25+
function RobinBC(l::NTuple{3,T}, r::NTuple{3,T}, dx::U, order = 1) where {T<:Number,U<:Real}
2626
αl, βl, γl = l
2727
αr, βr, γr = r
2828

29-
s = calculate_weights(1, one(T), Array(one(T):convert(T,order+1))) #generate derivative coefficients about the boundary of required approximation order
29+
s = calculate_weights(1, one(U), Array(one(U):convert(U,order+1))) #generate derivative coefficients about the boundary of required approximation order
3030

31-
a_l = -s[2:end]./(αl*dx/βl + s[1])
32-
a_r = s[end:-1:2]./(αr*dx/βr - s[1]) # for other boundary stencil is flippedlr with *opposite sign*
31+
a_l = -βl*s[2:end]./(αl*dx + βl*s[1])
32+
a_r = βr*s[end:-1:2]./(αr*dx - βr*s[1]) # for other boundary stencil is flippedlr with *opposite sign*
3333

3434
b_l = γl/(αl+βl*s[1]/dx)
3535
b_r = γr/(αr-βr*s[1]/dx)
3636

3737
return new{T, typeof(a_l)}(a_l, b_l, a_r, b_r)
3838
end
39-
function RobinBC(l::Union{NTuple{3,T},AbstractVector{T}}, r::Union{NTuple{3,T},AbstractVector{T}}, dx::AbstractVector{T}, order = 1) where {T}
39+
function RobinBC(l::Union{NTuple{3,T},AbstractVector{T}}, r::Union{NTuple{3,T},AbstractVector{T}}, dx::AbstractVector{U}, order = 1) where {T<:Number,U<:Real}
4040
αl, βl, γl = l
4141
αr, βr, γr = r
4242

43-
s_index = Array(one(T):convert(T,order+1))
43+
s_index = Array(one(U):convert(U,order+1))
4444
dx_l, dx_r = dx[1:length(s_index)], dx[(end-length(s_index)+1):end]
4545

46-
s = calculate_weights(1, one(T), s_index) #generate derivative coefficients about the boundary of required approximation order
46+
s = calculate_weights(1, one(U), s_index) #generate derivative coefficients about the boundary of required approximation order
4747
denom_l = αl+βl*s[1]/dx_l[1]
4848
denom_r = αr-βr*s[1]/dx_r[end]
4949

@@ -76,24 +76,24 @@ struct GeneralBC{T, L<:AbstractVector{T}, R<:AbstractVector{T}} <:AffineBC{T}
7676
b_l::T
7777
a_r::R
7878
b_r::T
79-
function GeneralBC(αl::AbstractVector{T}, αr::AbstractVector{T}, dx::T, order = 1) where {T}
79+
function GeneralBC(αl::AbstractVector{T}, αr::AbstractVector{T}, dx::U, order = 1) where {T<:Number,U<:Real}
8080
nl = length(αl)
8181
nr = length(αr)
8282
S_l = zeros(T, (nl-2, order+nl-2))
8383
S_r = zeros(T, (nr-2, order+nr-2))
8484

8585
for i in 1:(nl-2)
86-
S_l[i,:] = [transpose(calculate_weights(i, one(T), Array(one(T):convert(T, order+i)))) transpose(zeros(T, Int(nl-2-i)))]./(dx^i) #am unsure if the length of the dummy_x is correct here
86+
S_l[i,:] = [transpose(calculate_weights(i, one(U), Array(one(U):convert(U, order+i)))) transpose(zeros(U, Int(nl-2-i)))]./(dx^i) #am unsure if the length of the dummy_x is correct here
8787
end
8888

8989
for i in 1:(nr-2)
90-
S_r[i,:] = [transpose(calculate_weights(i, convert(T, order+i), Array(one(T):convert(T, order+i)))) transpose(zeros(T, Int(nr-2-i)))]./(dx^i)
90+
S_r[i,:] = [transpose(calculate_weights(i, convert(U, order+i), Array(one(U):convert(U, order+i)))) transpose(zeros(U, Int(nr-2-i)))]./(dx^i)
9191
end
9292
s0_l = S_l[:,1] ; Sl = S_l[:,2:end]
9393
s0_r = S_r[:,end] ; Sr = S_r[:,(end-1):-1:1]
9494

95-
denoml = αl[2] .+ αl[3:end] ⋅ s0_l
96-
denomr = αr[2] .+ αr[3:end] ⋅ s0_r
95+
denoml = αl[2] .+ αl[3:end]' ⋅ s0_l
96+
denomr = αr[2] .+ αr[3:end]' ⋅ s0_r
9797

9898
a_l = -transpose(transpose(αl[3:end]) * Sl) ./denoml
9999
a_r = reverse(-transpose(transpose(αr[3:end]) * Sr) ./denomr)
@@ -103,7 +103,7 @@ struct GeneralBC{T, L<:AbstractVector{T}, R<:AbstractVector{T}} <:AffineBC{T}
103103
new{T, typeof(a_l), typeof(a_r)}(a_l,b_l,a_r,b_r)
104104
end
105105

106-
function GeneralBC(αl::AbstractVector{T}, αr::AbstractVector{T}, dx::AbstractVector{T}, order = 1) where {T}
106+
function GeneralBC(αl::AbstractVector{T}, αr::AbstractVector{T}, dx::AbstractVector{U}, order = 1) where {T<:Number,U<:Real}
107107

108108
nl = length(αl)
109109
nr = length(αr)
@@ -112,17 +112,17 @@ struct GeneralBC{T, L<:AbstractVector{T}, R<:AbstractVector{T}} <:AffineBC{T}
112112
S_r = zeros(T, (nr-2, order+nr-2))
113113

114114
for i in 1:(nl-2)
115-
S_l[i,:] = [transpose(calculate_weights(i, one(T), Array(one(T):convert(T, order+i)))) transpose(zeros(T, Int(nl-2-i)))]./(dx_l.^i)
115+
S_l[i,:] = [transpose(calculate_weights(i, one(U), Array(one(U):convert(U, order+i)))) transpose(zeros(U, Int(nl-2-i)))]./(dx_l.^i)
116116
end
117117

118118
for i in 1:(nr-2)
119-
S_r[i,:] = [transpose(calculate_weights(i, convert(T, order+i), Array(one(T):convert(T, order+i)))) transpose(zeros(T, Int(nr-2-i)))]./(dx_r.^i)
119+
S_r[i,:] = [transpose(calculate_weights(i, convert(U, order+i), Array(one(U):convert(U, order+i)))) transpose(zeros(U, Int(nr-2-i)))]./(dx_r.^i)
120120
end
121121
s0_l = S_l[:,1] ; Sl = S_l[:,2:end]
122122
s0_r = S_r[:,end] ; Sr = S_r[:,(end-1):-1:1]
123123

124-
denoml = αl[2] .+ αl[3:end] ⋅ s0_l
125-
denomr = αr[2] .+ αr[3:end] ⋅ s0_r
124+
denoml = αl[2] .+ αl[3:end]' ⋅ s0_l
125+
denomr = αr[2] .+ αr[3:end]' ⋅ s0_r
126126

127127
a_l = -transpose(transpose(αl[3:end]) * Sl) ./denoml
128128
a_r = reverse(-transpose(transpose(αr[3:end]) * Sr) ./denomr)
@@ -136,16 +136,19 @@ end
136136

137137

138138
#implement Neumann and Dirichlet as special cases of RobinBC
139-
NeumannBC(α::NTuple{2,T}, dx::Union{AbstractVector{T}, T}, order = 1) where T = RobinBC((zero(T), one(T), α[1]), (zero(T), one(T), α[2]), dx, order)
140-
DirichletBC(αl::T, αr::T) where T = RobinBC((one(T), zero(T), αl), (one(T), zero(T), αr), one(T), 2one(T) )
139+
NeumannBC(α::NTuple{2,T}, dx::Union{AbstractVector{U}, U}, order = 1) where {T<:Number,U<:Real} = RobinBC((zero(T), one(T), α[1]), (zero(T), one(T), α[2]), dx, order)
140+
DirichletBC(αl::T, αr::T) where {T<:Real} = RobinBC((one(T), zero(T), αl), (one(T), zero(T), αr), one(T), 2one(T) )
141+
DirichletBC(::Type{U},αl::T, αr::T) where {T<:Number,U<:Real} = RobinBC((one(T), zero(T), αl), (one(T), zero(T), αr), one(U), 2one(U) )
141142
#specialized constructors for Neumann0 and Dirichlet0
142143
Dirichlet0BC(T::Type) = DirichletBC(zero(T), zero(T))
143-
Neumann0BC(dx::Union{AbstractVector{T}, T}, order = 1) where T = NeumannBC((zero(T), zero(T)), dx, order)
144+
Neumann0BC(dx::Union{AbstractVector{T}, T}, order = 1) where {T<:Real} = NeumannBC((zero(T), zero(T)), dx, order)
145+
Neumann0BC(::Type{U},dx::Union{AbstractVector{T}, T}, order = 1) where {T<:Real,U<:Number} = NeumannBC((zero(U), zero(U)), dx, order)
144146

145147
# other acceptable argument signatures
146148
#RobinBC(al::T, bl::T, cl::T, dx_l::T, ar::T, br::T, cr::T, dx_r::T, order = 1) where T = RobinBC([al,bl, cl], [ar, br, cr], dx_l, order)
147149

148-
Base.:*(Q::AffineBC, u::AbstractVector) = BoundaryPaddedVector(Q.a_l ⋅ u[1:length(Q.a_l)] + Q.b_l, Q.a_r ⋅ u[(end-length(Q.a_r)+1):end] + Q.b_r, u)
150+
Base.:*(Q::AffineBC, u::AbstractVector) = BoundaryPaddedVector( Q.a_l'⋅ u[1:length(Q.a_l)] + Q.b_l, Q.a_r' ⋅ u[(end-length(Q.a_r)+1):end] + Q.b_r, u )
151+
149152
Base.:*(Q::PeriodicBC, u::AbstractVector) = BoundaryPaddedVector(u[end], u[1], u)
150153

151154
Base.size(Q::AtomicBC) = (Inf, Inf) #Is this nessecary?

Diff for: ‎test/robin.jl

+112
Original file line numberDiff line numberDiff line change
@@ -104,3 +104,115 @@ ugeneralextended = G*u
104104

105105

106106
#TODO: Implement tests for BC's that are contingent on the sign of the coefficient on the operator near the boundary
107+
108+
109+
110+
111+
# Test complex Robin BC, uniform grid
112+
113+
# Generate random parameters
114+
al = rand(ComplexF64,5)
115+
bl = rand(ComplexF64,5)
116+
cl = rand(ComplexF64,5)
117+
dx = rand(Float64,5)
118+
ar = rand(ComplexF64,5)
119+
br = rand(ComplexF64,5)
120+
cr = rand(ComplexF64,5)
121+
122+
# Construct 5 arbitrary RobinBC operators for each parameter set
123+
for i in 1:5
124+
125+
Q = RobinBC((al[i], bl[i], cl[i]), (ar[i], br[i], cr[i]), dx[i])
126+
127+
Q_L, Q_b = Array(Q,5i)
128+
129+
#Check that Q_L is is correctly computed
130+
@test Q_L[2:5i+1,1:5i] ≈ Array(I, 5i, 5i)
131+
@test Q_L[1,:] ≈ [1 / (1-al[i]*dx[i]/bl[i]); zeros(5i-1)]
132+
@test Q_L[5i+2,:] ≈ [zeros(5i-1); 1 / (1+ar[i]*dx[i]/br[i])]
133+
134+
#Check that Q_b is computed correctly
135+
@test Q_b ≈ [cl[i]/(al[i]-bl[i]/dx[i]); zeros(5i); cr[i]/(ar[i]+br[i]/dx[i])]
136+
137+
# Construct the extended operator and check that it correctly extends u to a (5i+2)
138+
# vector, along with encoding boundary condition information.
139+
u = rand(ComplexF64,5i)
140+
141+
Qextended = Q*u
142+
CorrectQextended = [(cl[i]-(bl[i]/dx[i])*u[1])/(al[i]-bl[i]/dx[i]); u; (cr[i]+ (br[i]/dx[i])*u[5i])/(ar[i]+br[i]/dx[i])]
143+
@test length(Qextended) ≈ 5i+2
144+
145+
# Check concretization
146+
@test Array(Qextended) ≈ CorrectQextended # Q.a_l ⋅ u[1:length(Q.a_l)] + Q.b_l, Q.a_r ⋅ u[(end-length(Q.a_r)+1):end] + Q.b_r
147+
148+
# Check that Q_L and Q_b correctly compute BoundaryPaddedVector
149+
@test Q_L*u + Q_b ≈ CorrectQextended
150+
151+
@test [Qextended[1]; Qextended.u; Qextended[5i+2]] ≈ CorrectQextended
152+
153+
end
154+
155+
# Test complex Robin BC, w/non-uniform grid
156+
al = rand(ComplexF64,5)
157+
bl = rand(ComplexF64,5)
158+
cl = rand(ComplexF64,5)
159+
dx = rand(Float64,5)
160+
ar = rand(ComplexF64,5)
161+
br = rand(ComplexF64,5)
162+
cr = rand(ComplexF64,5)
163+
for j in 1:2
164+
for i in 1:5
165+
if j == 1
166+
Q = RobinBC((al[i], bl[i], cl[i]), (ar[i], br[i], cr[i]),
167+
dx[i] .* ones(5 * i))
168+
else
169+
Q = RobinBC([al[i], bl[i], cl[i]], [ar[i], br[i], cr[i]],
170+
dx[i] .* ones(5 * i))
171+
end
172+
Q_L, Q_b = Array(Q,5i)
173+
174+
#Check that Q_L is is correctly computed
175+
@test Q_L[2:5i+1,1:5i] ≈ Array(I, 5i, 5i)
176+
@test Q_L[1,:] ≈ [1 / (1-al[i]*dx[i]/bl[i]); zeros(5i-1)]
177+
@test Q_L[5i+2,:] ≈ [zeros(5i-1); 1 / (1+ar[i]*dx[i]/br[i])]
178+
179+
#Check that Q_b is computed correctly
180+
@test Q_b ≈ [cl[i]/(al[i]-bl[i]/dx[i]); zeros(5i); cr[i]/(ar[i]+br[i]/dx[i])]
181+
182+
# Construct the extended operator and check that it correctly extends u to a (5i+2)
183+
# vector, along with encoding boundary condition information.
184+
u = rand(ComplexF64,5i)
185+
186+
Qextended = Q*u
187+
CorrectQextended = [(cl[i]-(bl[i]/dx[i])*u[1])/(al[i]-bl[i]/dx[i]); u; (cr[i]+ (br[i]/dx[i])*u[5i])/(ar[i]+br[i]/dx[i])]
188+
@test length(Qextended) ≈ 5i+2
189+
190+
# Check concretization
191+
@test Array(Qextended) ≈ CorrectQextended
192+
193+
# Check that Q_L and Q_b correctly compute BoundaryPaddedVector
194+
@test Q_L*u + Q_b ≈ CorrectQextended
195+
196+
@test [Qextended[1]; Qextended.u; Qextended[5i+2]] ≈ CorrectQextended
197+
198+
end
199+
end
200+
201+
# Test NeumannBC, DirichletBC as special cases of RobinBC
202+
let
203+
dx = [0.121, 0.783, 0.317, 0.518, 0.178]
204+
αC = (0.539 + 0.653im, 0.842 + 0.47im)
205+
αR = (0.045, 0.577)
206+
@test NeumannBC(αC, dx).b_l ≈ -0.065219 - 0.079013im
207+
@test DirichletBC(αR...).b_r ≈ 0.577
208+
@test DirichletBC(Float64, αC...).b_l ≈ 0.539 + 0.653im
209+
@test DirichletBC(Float64, αC...).a_r ≈ [-0.0 + 0.0im, 0.0 + 0.0im]
210+
211+
@test Dirichlet0BC(Float64).a_r ≈ [-0.0,0.0]
212+
@test Neumann0BC(dx).a_r ≈ [0.3436293436293436]
213+
@test Neumann0BC(ComplexF64,dx).a_l ≈ [0.15453384418901658 + 0.0im]
214+
215+
@test NeumannBC(αC, first(dx)).b_r ≈ 0.101882 + 0.05687im
216+
@test Neumann0BC(first(dx)).a_r ≈ [1.0 - 0.0im]
217+
@test Neumann0BC(ComplexF64,first(dx)).a_l ≈ [1.0 + 0.0im]
218+
end

0 commit comments

Comments
 (0)