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

[WIP] (issue #237) Allow distinct type for grid vs. function value in BC constructor #260

Open
wants to merge 11 commits into
base: master
Choose a base branch
from
Open

[WIP] (issue #237) Allow distinct type for grid vs. function value in BC constructor #260

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
11 commits
Select commit Hold shift + click to select a range
0db40b7
Add distinct type for grid vs. function value
Aug 8, 2020
b6b39f6
Test Robin BC with complex function values
Aug 8, 2020
0d3d564
duplicate comment
Aug 8, 2020
0f63558
Add new test to runtests
Aug 8, 2020
e8a2d96
Restrict types to avoid ambiguity
Aug 8, 2020
b558da8
More tests, type restrictions
Aug 8, 2020
6106632
Fix Dirichlet / Neumann
Aug 8, 2020
d477f23
Add tests for Neuman/Dirichlet w. Complex numbers
Aug 9, 2020
83a9b6c
Add distinct type for grid vs. function value
AndiMD Aug 8, 2020
c6f45a5
Merge remote-tracking branch 'AndiMD/master'
Aug 12, 2020
197c9f2
Express unconjugated dot product by sum(a.*b)
Sep 5, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
45 changes: 24 additions & 21 deletions src/derivative_operators/BC_operators.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,28 +22,28 @@ struct RobinBC{T, V<:AbstractVector{T}} <: AffineBC{T}
b_l::T
a_r::V
b_r::T
function RobinBC(l::NTuple{3,T}, r::NTuple{3,T}, dx::T, order = 1) where {T}
function RobinBC(l::NTuple{3,T}, r::NTuple{3,T}, dx::U, order = 1) where {T<:Number,U<:Real}
αl, βl, γl = l
αr, βr, γr = r

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

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

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

return new{T, typeof(a_l)}(a_l, b_l, a_r, b_r)
end
function RobinBC(l::Union{NTuple{3,T},AbstractVector{T}}, r::Union{NTuple{3,T},AbstractVector{T}}, dx::AbstractVector{T}, order = 1) where {T}
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}
αl, βl, γl = l
αr, βr, γr = r

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

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

Expand Down Expand Up @@ -76,24 +76,24 @@ struct GeneralBC{T, L<:AbstractVector{T}, R<:AbstractVector{T}} <:AffineBC{T}
b_l::T
a_r::R
b_r::T
function GeneralBC(αl::AbstractVector{T}, αr::AbstractVector{T}, dx::T, order = 1) where {T}
function GeneralBC(αl::AbstractVector{T}, αr::AbstractVector{T}, dx::U, order = 1) where {T<:Number,U<:Real}
nl = length(αl)
nr = length(αr)
S_l = zeros(T, (nl-2, order+nl-2))
S_r = zeros(T, (nr-2, order+nr-2))

for i in 1:(nl-2)
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
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
end

for i in 1:(nr-2)
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)
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)
end
s0_l = S_l[:,1] ; Sl = S_l[:,2:end]
s0_r = S_r[:,end] ; Sr = S_r[:,(end-1):-1:1]

denoml = αl[2] .+ αl[3:end] ⋅ s0_l
denomr = αr[2] .+ αr[3:end] ⋅ s0_r
denoml = αl[2] .+ αl[3:end]' ⋅ s0_l
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

why adjoint?

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

⋅(a,b) complex conjugates a. I believe this should not be complex conjugated.

Should probably have been conj. instead of adjoint (same result due to the dot).

denomr = αr[2] .+ αr[3:end]' ⋅ s0_r

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

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

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

for i in 1:(nl-2)
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)
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)
end

for i in 1:(nr-2)
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)
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)
end
s0_l = S_l[:,1] ; Sl = S_l[:,2:end]
s0_r = S_r[:,end] ; Sr = S_r[:,(end-1):-1:1]

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

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


#implement Neumann and Dirichlet as special cases of RobinBC
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)
DirichletBC(αl::T, αr::T) where T = RobinBC((one(T), zero(T), αl), (one(T), zero(T), αr), one(T), 2one(T) )
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)
DirichletBC(αl::T, αr::T) where {T<:Real} = RobinBC((one(T), zero(T), αl), (one(T), zero(T), αr), one(T), 2one(T) )
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) )
#specialized constructors for Neumann0 and Dirichlet0
Dirichlet0BC(T::Type) = DirichletBC(zero(T), zero(T))
Neumann0BC(dx::Union{AbstractVector{T}, T}, order = 1) where T = NeumannBC((zero(T), zero(T)), dx, order)
Neumann0BC(dx::Union{AbstractVector{T}, T}, order = 1) where {T<:Real} = NeumannBC((zero(T), zero(T)), dx, order)
Neumann0BC(::Type{U},dx::Union{AbstractVector{T}, T}, order = 1) where {T<:Real,U<:Number} = NeumannBC((zero(U), zero(U)), dx, order)

# other acceptable argument signatures
#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)

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)
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 )
Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The idea that the dot-product should not complex conjugate came from this line. I figured it should be the same for GeneralBC


Base.:*(Q::PeriodicBC, u::AbstractVector) = BoundaryPaddedVector(u[end], u[1], u)

Base.size(Q::AtomicBC) = (Inf, Inf) #Is this nessecary?
Expand Down
112 changes: 112 additions & 0 deletions test/robin.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -104,3 +104,115 @@ ugeneralextended = G*u


#TODO: Implement tests for BC's that are contingent on the sign of the coefficient on the operator near the boundary




# Test complex Robin BC, uniform grid

# Generate random parameters
al = rand(ComplexF64,5)
bl = rand(ComplexF64,5)
cl = rand(ComplexF64,5)
dx = rand(Float64,5)
ar = rand(ComplexF64,5)
br = rand(ComplexF64,5)
cr = rand(ComplexF64,5)

# Construct 5 arbitrary RobinBC operators for each parameter set
for i in 1:5

Q = RobinBC((al[i], bl[i], cl[i]), (ar[i], br[i], cr[i]), dx[i])

Q_L, Q_b = Array(Q,5i)

#Check that Q_L is is correctly computed
@test Q_L[2:5i+1,1:5i] ≈ Array(I, 5i, 5i)
@test Q_L[1,:] ≈ [1 / (1-al[i]*dx[i]/bl[i]); zeros(5i-1)]
@test Q_L[5i+2,:] ≈ [zeros(5i-1); 1 / (1+ar[i]*dx[i]/br[i])]

#Check that Q_b is computed correctly
@test Q_b ≈ [cl[i]/(al[i]-bl[i]/dx[i]); zeros(5i); cr[i]/(ar[i]+br[i]/dx[i])]

# Construct the extended operator and check that it correctly extends u to a (5i+2)
# vector, along with encoding boundary condition information.
u = rand(ComplexF64,5i)

Qextended = Q*u
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])]
@test length(Qextended) ≈ 5i+2

# Check concretization
@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

# Check that Q_L and Q_b correctly compute BoundaryPaddedVector
@test Q_L*u + Q_b ≈ CorrectQextended

@test [Qextended[1]; Qextended.u; Qextended[5i+2]] ≈ CorrectQextended

end

# Test complex Robin BC, w/non-uniform grid
al = rand(ComplexF64,5)
bl = rand(ComplexF64,5)
cl = rand(ComplexF64,5)
dx = rand(Float64,5)
ar = rand(ComplexF64,5)
br = rand(ComplexF64,5)
cr = rand(ComplexF64,5)
for j in 1:2
for i in 1:5
if j == 1
Q = RobinBC((al[i], bl[i], cl[i]), (ar[i], br[i], cr[i]),
dx[i] .* ones(5 * i))
else
Q = RobinBC([al[i], bl[i], cl[i]], [ar[i], br[i], cr[i]],
dx[i] .* ones(5 * i))
end
Q_L, Q_b = Array(Q,5i)

#Check that Q_L is is correctly computed
@test Q_L[2:5i+1,1:5i] ≈ Array(I, 5i, 5i)
@test Q_L[1,:] ≈ [1 / (1-al[i]*dx[i]/bl[i]); zeros(5i-1)]
@test Q_L[5i+2,:] ≈ [zeros(5i-1); 1 / (1+ar[i]*dx[i]/br[i])]

#Check that Q_b is computed correctly
@test Q_b ≈ [cl[i]/(al[i]-bl[i]/dx[i]); zeros(5i); cr[i]/(ar[i]+br[i]/dx[i])]

# Construct the extended operator and check that it correctly extends u to a (5i+2)
# vector, along with encoding boundary condition information.
u = rand(ComplexF64,5i)

Qextended = Q*u
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])]
@test length(Qextended) ≈ 5i+2

# Check concretization
@test Array(Qextended) ≈ CorrectQextended

# Check that Q_L and Q_b correctly compute BoundaryPaddedVector
@test Q_L*u + Q_b ≈ CorrectQextended

@test [Qextended[1]; Qextended.u; Qextended[5i+2]] ≈ CorrectQextended

end
end

# Test NeumannBC, DirichletBC as special cases of RobinBC
let
dx = [0.121, 0.783, 0.317, 0.518, 0.178]
αC = (0.539 + 0.653im, 0.842 + 0.47im)
αR = (0.045, 0.577)
@test NeumannBC(αC, dx).b_l ≈ -0.065219 - 0.079013im
@test DirichletBC(αR...).b_r ≈ 0.577
@test DirichletBC(Float64, αC...).b_l ≈ 0.539 + 0.653im
@test DirichletBC(Float64, αC...).a_r ≈ [-0.0 + 0.0im, 0.0 + 0.0im]

@test Dirichlet0BC(Float64).a_r ≈ [-0.0,0.0]
@test Neumann0BC(dx).a_r ≈ [0.3436293436293436]
@test Neumann0BC(ComplexF64,dx).a_l ≈ [0.15453384418901658 + 0.0im]

@test NeumannBC(αC, first(dx)).b_r ≈ 0.101882 + 0.05687im
@test Neumann0BC(first(dx)).a_r ≈ [1.0 - 0.0im]
@test Neumann0BC(ComplexF64,first(dx)).a_l ≈ [1.0 + 0.0im]
end