Add distinct type for grid vs. function value

Test Robin BC, Dirichlet, Neumann with complex function values

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] .+ α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)
@@ -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] .+ α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)
@@ -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 )
+
 Base.:*(Q::PeriodicBC, u::AbstractVector) = BoundaryPaddedVector(u[end], u[1], u)
 
 Base.size(Q::AtomicBC) = (Inf, Inf) #Is this nessecary?

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