diff --git a/src/derivative_operators/BC_operators.jl b/src/derivative_operators/BC_operators.jl index 0f05711e0..debcb6f74 100644 --- a/src/derivative_operators/BC_operators.jl +++ b/src/derivative_operators/BC_operators.jl @@ -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] @@ -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] .+ sum(αl[3:end] .* s0_l) # dot product without complex conjugation + denomr = αr[2] .+ sum(αr[3:end] .* s0_r) a_l = -transpose(transpose(αl[3:end]) * Sl) ./denoml 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} 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) @@ -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] .+ sum(αl[3:end] .* s0_l) + denomr = αr[2] .+ sum(αr[3:end] .* s0_r) a_l = -transpose(transpose(αl[3:end]) * Sl) ./denoml a_r = reverse(-transpose(transpose(αr[3:end]) * Sr) ./denomr) @@ -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( sum(Q.a_l .* u[1:length(Q.a_l)]) + Q.b_l, sum(Q.a_r .* u[(end-length(Q.a_r)+1):end]) + Q.b_r, u ) + Base.:*(Q::PeriodicBC, u::AbstractVector) = BoundaryPaddedVector(u[end], u[1], u) Base.size(Q::AtomicBC) = (Inf, Inf) #Is this nessecary? diff --git a/test/robin.jl b/test/robin.jl index d7c8b3882..aaf115036 100644 --- a/test/robin.jl +++ b/test/robin.jl @@ -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 diff --git a/test/robin_complex.jl b/test/robin_complex.jl new file mode 100644 index 000000000..e0a737f04 --- /dev/null +++ b/test/robin_complex.jl @@ -0,0 +1,45 @@ +include("src/DiffEqOperators.jl") +using .DiffEqOperators +using LinearAlgebra, Random, Test + +# 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