Robin tests and some fixes to RobinBC

src/derivative_operators/robin_bc_extended.jl

@@ -32,6 +32,7 @@ end
 
 Base.:*(Q::RobinBC,u) = RobinBCExtended(u, Q.al, Q.bl, Q.cl, Q.dx_l, Q.ar, Q.br, Q.cr, Q.dx_r)
 Base.length(Q::RobinBCExtended) = length(Q.u) + 2
+Base.size(Q::RobinBCExtended) = (length(Q.u)+2,)
 Base.lastindex(Q::RobinBCExtended) = Base.length(Q)
 
 function Base.getindex(Q::RobinBCExtended,i)
@@ -45,8 +46,8 @@ function Base.getindex(Q::RobinBCExtended,i)
 end
 
 function LinearAlgebra.Array(Q::RobinBC, N::Int)
-    Q_L = [(-Q.bl/Q.dx_l)/(Q.al-Q.bl/Q.dx_l) transpose(zeros(N-1)); Diagonal(ones(N)); (Q.br/Q.dx_r)/(Q.ar+Q.br/Q.dx_r) transpose(zeros(N-1))]
-    Q_b = [Q.cl; zeros(N); Q.cr]
+    Q_L = [(-Q.bl/Q.dx_l)/(Q.al-Q.bl/Q.dx_l) transpose(zeros(N-1)); Diagonal(ones(N)); transpose(zeros(N-1)) (Q.br/Q.dx_r)/(Q.ar+Q.br/Q.dx_r)]
+    Q_b = [Q.cl/(Q.al-Q.bl/Q.dx_l); zeros(N); Q.cr/(Q.ar+Q.br/Q.dx_r)]
     return (Q_L, Q_b)
 end
 

test/robin.jl

@@ -0,0 +1,48 @@
+using LinearAlgebra, DiffEqOperators, Random, Test
+
+# Generate random parameters
+al = rand(5)
+bl = rand(5)
+cl = rand(5)
+dx_l = rand(5)
+ar = rand(5)
+br = rand(5)
+cr = rand(5)
+dx_r = rand(5)
+
+# Construct 5 arbitrary RobinBC operators
+for i in 1:5
+    Q = RobinBC(al[i], bl[i], cl[i], dx_l[i], ar[i], br[i], cr[i], dx_r[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_l[i]/bl[i]); zeros(5i-1)]
+    @test Q_L[5i+2,:] ≈ [zeros(5i-1); 1 / (1+ar[i]*dx_r[i]/br[i])]
+
+
+    #Check that Q_b is computed correctly
+    @test Q_b ≈ [cl[i]/(al[i]-bl[i]/dx_l[i]); zeros(5i); cr[i]/(ar[i]+br[i]/dx_r[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(5i)
+
+    Qextended = Q*u
+    CorrectQextended = [(cl[i]-(bl[i]/dx_l[i])*u[1])/(al[i]-bl[i]/dx_l[i]); u; (cr[i]+ (br[i]/dx_r[i])*u[5i])/(ar[i]+br[i]/dx_r[i])]
+
+    @test length(Qextended) ≈ 5i+2
+    @test Qextended ≈ CorrectQextended
+
+    # Check concretization
+    @test Array(Qextended) ≈ CorrectQextended
+
+    # Check that Q_L and Q_b correctly compute RobinBCExtended
+    @test Q_L*u + Q_b ≈ CorrectQextended
+
+    @test [Qextended[1]; Qextended.u; Qextended[5i+2]] ≈ CorrectQextended
+
+end
+
+#TODO: Implement tests for BC's that are contingent on the sign of the coefficient on the operator near the boundary

test/runtests.jl

@@ -2,6 +2,7 @@ using SafeTestsets
 import Base: isapprox
 
 @time @safetestset "Basic Operators Interface" begin include("basic_operators_interface.jl") end
+@time @safetestset "Robin Boundary Condition Operators" begin include("robin.jl") end
 @time @safetestset "JacVec Operators Interface" begin include("jacvec_operators.jl") end
 @time @safetestset "Composite Operators Interface" begin include("composite_operators_interface.jl") end