rename to CenteredDifference and simplify constructor

Author: ChrisRackauckas

src/DiffEqOperators.jl

@@ -44,6 +44,6 @@ end
 
 export MatrixFreeOperator
 export DiffEqScalar, DiffEqArrayOperator, DiffEqIdentity, JacVecOperator, getops
-export AbstractDerivativeOperator, DerivativeOperator, UpwindOperator, FiniteDifference
-export RobinBC
+export AbstractDerivativeOperator, DerivativeOperator, CenteredDifference, UpwindOperator, FiniteDifference
+export RobinBC, GeneralBC
 end # module

src/derivative_operators/derivative_operator.jl

@@ -9,36 +9,33 @@ struct DerivativeOperator{T<:Real,S1<:SVector,S2<:SVector} <: AbstractDerivative
     boundary_point_count    :: Int
     low_boundary_coefs      :: Vector{S2}
     high_boundary_coefs     :: Vector{S2}
+end
 
-    function DerivativeOperator{T,S1,S2}(derivative_order::Int,
-                                     approximation_order::Int, dx::T,
-                                     dimension::Int) where
-                                     {T<:Real,S1<:SVector,S2<:SVector}
-        dimension               = dimension
-        dx                      = dx
-        stencil_length          = derivative_order + approximation_order - 1 + (derivative_order+approximation_order)%2
-        boundary_stencil_length = derivative_order + approximation_order
-        dummy_x                 = -div(stencil_length,2) : div(stencil_length,2)
-        boundary_x              = -boundary_stencil_length+1:0
-        boundary_point_count    = div(stencil_length,2) - 1 # -1 due to the ghost point
-        # Because it's a N x (N+2) operator, the last stencil on the sides are the [b,0,x,x,x,x] stencils, not the [0,x,x,x,x,x] stencils, since we're never solving for the derivative at the boundary point.
-        deriv_spots             = (-div(stencil_length,2)+1) : -1
-        boundary_deriv_spots    = boundary_x[2:div(stencil_length,2)]
+function CenteredDifference(derivative_order::Int,
+                            approximation_order::Int, dx::T,
+                            dimension::Int) where {T<:Real}
+    dimension               = dimension
+    dx                      = dx
+    stencil_length          = derivative_order + approximation_order - 1 + (derivative_order+approximation_order)%2
+    boundary_stencil_length = derivative_order + approximation_order
+    dummy_x                 = -div(stencil_length,2) : div(stencil_length,2)
+    boundary_x              = -boundary_stencil_length+1:0
+    boundary_point_count    = div(stencil_length,2) - 1 # -1 due to the ghost point
+    # Because it's a N x (N+2) operator, the last stencil on the sides are the [b,0,x,x,x,x] stencils, not the [0,x,x,x,x,x] stencils, since we're never solving for the derivative at the boundary point.
+    deriv_spots             = (-div(stencil_length,2)+1) : -1
+    boundary_deriv_spots    = boundary_x[2:div(stencil_length,2)]
 
-        stencil_coefs           = convert(SVector{stencil_length, T}, calculate_weights(derivative_order, zero(T), dummy_x))
-        low_boundary_coefs      = [convert(SVector{boundary_stencil_length, T}, calculate_weights(derivative_order, oneunit(T)*x0, boundary_x)) for x0 in boundary_deriv_spots]
-        high_boundary_coefs     = reverse!(copy([convert(SVector{boundary_stencil_length, T}, reverse(low_boundary_coefs[i])) for i in 1:boundary_point_count]))
+    stencil_coefs           = convert(SVector{stencil_length, T}, calculate_weights(derivative_order, zero(T), dummy_x))
+    low_boundary_coefs      = SVector{boundary_stencil_length, T}[convert(SVector{boundary_stencil_length, T}, calculate_weights(derivative_order, oneunit(T)*x0, boundary_x)) for x0 in boundary_deriv_spots]
+    high_boundary_coefs     = reverse([reverse(low_boundary_coefs[i]) for i in 1:boundary_point_count])
 
-        new(derivative_order, approximation_order, dx, dimension, stencil_length,
-            stencil_coefs,
-            boundary_stencil_length,
-            boundary_point_count,
-            low_boundary_coefs,
-            high_boundary_coefs
-            )
-    end
-    DerivativeOperator{T}(dorder::Int,aorder::Int,dx::T,dim::Int) where {T<:Real} =
-        DerivativeOperator{T, SVector{dorder+aorder-1+(dorder+aorder)%2,T}, SVector{dorder+aorder,T}}(dorder, aorder, dx, dim)
+    DerivativeOperator(derivative_order, approximation_order, dx, dimension, stencil_length,
+        stencil_coefs,
+        boundary_stencil_length,
+        boundary_point_count,
+        low_boundary_coefs,
+        high_boundary_coefs
+        )
 end
 
 #=

test/2nd_order_check.jl

@@ -5,8 +5,8 @@ using DiffEqOperators, Test, LinearAlgebra
   Δx = 0.025π
   N = Int(2*(π/Δx)) -2
   x = -π:Δx:π-Δx
-  D1  = DerivativeOperator{Float64}(1,ord,Δx,N,:periodic,:periodic);    # 1nd Derivative
-  D2  = DerivativeOperator{Float64}(2,ord,Δx,N,:periodic,:periodic);    # 2nd Derivative
+  D1  = CenteredDifference(1,ord,Δx,N,:periodic,:periodic);    # 1nd Derivative
+  D2  = CenteredDifference(2,ord,Δx,N,:periodic,:periodic);    # 2nd Derivative
 
   u0 = cos.(x)
   du_true = -cos.(x)

test/KdV.jl

@@ -17,10 +17,10 @@ using DiffEqOperators, OrdinaryDiffEq
     const du3 = zeros(size(x));
     const temp = zeros(size(x));
 
-    # A = DerivativeOperator{Float64}(1,2,Δx,length(x),:Dirichlet0,:Dirichlet0);
+    # A = CenteredDifference(1,2,Δx,length(x),:Dirichlet0,:Dirichlet0);
     A = UpwindOperator{Float64}(1,3,Δx,length(x),true.|BitVector(undef,length(x)),
                                 :Dirichlet0,:Dirichlet0);
-    # C = DerivativeOperator{Float64}(3,2,Δx,length(x),:Dirichlet0,:Dirichlet0);
+    # C = CenteredDifference(3,2,Δx,length(x),:Dirichlet0,:Dirichlet0);
     C = UpwindOperator{Float64}(3,3,Δx,length(x),true.|BitVector(undef,length(x)),
                                 :Dirichlet0,:Dirichlet0);
 

test/derivative_operators_interface.jl

@@ -53,7 +53,7 @@ end
 # Do not modify the following test-set unless you are completely certain of your changes.
 @testset "Correctness of Stencils" begin
     N = 20
-    L = DerivativeOperator{Float64}(4,4, 1.0, N)
+    L = CenteredDifference(4,4, 1.0, N)
     correct = fourth_deriv_approx_stencil(N)
 
     # Check that stencils (according to convert_by_multiplication) agree with correct
@@ -64,7 +64,7 @@ end
     @test sparse(L) ≈ correct
     @test BandedMatrix(L) ≈ correct
 
-    L = DerivativeOperator{Float64}(2,4, 1.0, N)
+    L = CenteredDifference(2,4, 1.0, N)
     correct = second_deriv_fourth_approx_stencil(N)
 
     # Check that stencils (according to convert_by_multiplication) agree with correct
@@ -82,7 +82,7 @@ end
     d_order = 2
     approx_order = 2
     correct = second_derivative_stencil(N)
-    A = DerivativeOperator{Float64}(d_order,approx_order,1.0,N)
+    A = CenteredDifference(d_order,approx_order,1.0,N)
 
     @test convert_by_multiplication(Array,A,N) == correct
     @test Array(A) == second_derivative_stencil(N)
@@ -95,7 +95,7 @@ end
     N = 20
     d_order = 4
     approx_order = 4
-    A = DerivativeOperator{Float64}(d_order,approx_order,1.0,N)
+    A = CenteredDifference(d_order,approx_order,1.0,N)
     correct = convert_by_multiplication(Array,A,N)
 
     @test Array(A) ≈ correct
@@ -105,7 +105,7 @@ end
     N = 26
     d_order = 8
     approx_order = 8
-    A = DerivativeOperator{Float64}(d_order,approx_order,1.0,N)
+    A = CenteredDifference(d_order,approx_order,1.0,N)
     correct = convert_by_multiplication(Array,A,N)
 
     @test Array(A) ≈ correct
@@ -119,12 +119,11 @@ end
     y = collect(1:1.0:N+2).^4 - 2*collect(1:1.0:N+2).^3 + collect(1:1.0:N+2).^2;
     y = convert(Array{BigFloat, 1}, y)
 
-    A = DerivativeOperator{BigFloat}(d_order,approx_order,one(BigFloat),N)
+    A = CenteredDifference(d_order,approx_order,one(BigFloat),N)
     correct = convert_by_multiplication(Array,A,N)
     @test Array(A) ≈ correct
     @test sparse(A) ≈ correct
     @test BandedMatrix(A) ≈ correct
-
     @test A*y ≈ Array(A)*y
 end
 
@@ -133,7 +132,7 @@ end
     d_order = 4
     approx_order = 10
 
-    A = DerivativeOperator{Float64}(d_order,approx_order,1.0,N)
+    A = CenteredDifference(d_order,approx_order,1.0,N)
     @test A[1,1] == Array(A)[1,1]
     @test A[10,20] == 0
 
@@ -147,7 +146,7 @@ end
     d_order = 2
     approx_order = 2
 
-    A = DerivativeOperator{Float64}(d_order,approx_order,1.0,N)
+    A = CenteredDifference(d_order,approx_order,1.0,N)
     M = Array(A,1000)
     @test A[1,1] == M[1,1]
     @test A[1:4,1] == M[1:4,1]
@@ -167,8 +166,8 @@ end
     dy = yarr[2]-yarr[1]
     F = [x^2+y for x = xarr, y = yarr]
 
-    A = DerivativeOperator{Float64}(d_order,approx_order,dx,length(xarr)-2)
-    B = DerivativeOperator{Float64}(d_order,approx_order,dy,length(yarr))
+    A = CenteredDifference(d_order,approx_order,dx,length(xarr)-2)
+    B = CenteredDifference(d_order,approx_order,dy,length(yarr))
 
 
     @test A*F ≈ 2*ones(N-2,M) atol=1e-2
@@ -186,7 +185,7 @@ end
     # Only tests the additional functionality defined in "operator_combination.jl"
     N = 10
     Random.seed!(0); LA = DiffEqArrayOperator(rand(N,N))
-    LD = DerivativeOperator{Float64}(2,2,1.0,N)
+    LD = CenteredDifference(2,2,1.0,N)
     @test_broken begin
       L = 1.1*LA - 2.2*LD + 3.3*I
       # Builds convert(L) the brute-force way

test/generic_operator_validation.jl

@@ -11,7 +11,7 @@ y_ = y[2:(end-1)]
 
 @test_broken for dor in 1:6, aor in 1:8
 
-    D1 = DerivativeOperator{Float64}(dor,aor,dx[1],length(x))
+    D1 = CenteredDifference(dor,aor,dx[1],length(x))
     D2 = DiffEqOperators.FiniteDifference{Float64}(dor,aor,dx,length(x))
     D = (D1,D2)
     @test_broken convert(Array, D1) ≈ convert(Array, D2)

test/heat_eqn.jl

@@ -4,7 +4,7 @@ using OrdinaryDiffEq
 @testset "Parabolic Heat Equation with Dirichlet BCs" begin
     x = collect(-pi : 2pi/511 : pi)
     u0 = @. -(x - 0.5).^2 + 1/12
-    A = DerivativeOperator{Float64}(2,2,2π/511,512,:Dirichlet,:Dirichlet;BC=(u0[1],u0[end]))
+    A = CenteredDifference(2,2,2π/511,512,:Dirichlet,:Dirichlet;BC=(u0[1],u0[end]))
     heat_eqn = ODEProblem(A, u0, (0.,10.))
     soln = solve(heat_eqn,Tsit5(),dense=false,tstops=0:0.01:10)
 
@@ -20,10 +20,10 @@ end
     dx = 2π/(N-1)
     x = collect(-pi : dx : pi)
     u0 = @. -(x - 0.5)^2 + 1/12
-    B = DerivativeOperator{Float64}(1,2,dx,N,:None,:None)
+    B = CenteredDifference(1,2,dx,N,:None,:None)
     deriv_start, deriv_end = (B*u0)[1], (B*u0)[end]
 
-    A = DerivativeOperator{Float64}(2,2,dx,N,:Neumann,:Neumann;BC=(deriv_start,deriv_end))
+    A = CenteredDifference(2,2,dx,N,:Neumann,:Neumann;BC=(deriv_start,deriv_end))
 
     heat_eqn = ODEProblem(A, u0, (0.,10.))
     soln = solve(heat_eqn,Tsit5(),dense=false,tstops=0:0.01:10)
@@ -42,14 +42,14 @@ end
     dx = 2π/(N-1)
     x = collect(-pi : dx : pi)
     u0 = @. -(x - 0.5).^2 + 1/12
-    B = DerivativeOperator{Float64}(1,2,dx,N,:None,:None)
+    B = CenteredDifference(1,2,dx,N,:None,:None)
     deriv_start, deriv_end = (B*u0)[1], (B*u0)[end]
     params = [1.0,0.5]
 
     left_RBC = params[1]*u0[1] - params[2]*deriv_start
     right_RBC = params[1]*u0[end] + params[2]*deriv_end
 
-    A = DerivativeOperator{Float64}(2,2,dx,N,:Robin,:Dirichlet;BC=((params[1],params[2],left_RBC),u0[end]));
+    A = CenteredDifference(2,2,dx,N,:Robin,:Dirichlet;BC=((params[1],params[2],left_RBC),u0[end]));
 
     heat_eqn = ODEProblem(A, u0, (0.,10.));
     soln = solve(heat_eqn,Tsit5(),dense=false,tstops=0:0.01:10);

test/runtests.jl

@@ -11,31 +11,3 @@ import Base: isapprox
 #@time @safetestset "KdV" begin include("KdV.jl") end # KdV times out and all fails
 #@time @safetestset "Heat Equation" begin include("heat_eqn.jl") end
 @time @safetestset "Matrix-Free Operators" begin include("matrixfree.jl") end
-
-#=
-using StaticArrays, DiffEqOperators
-function isapprox(x::DerivativeOperator{T,S,LBC,RBC},y::FiniteDifference{T,S,LBC,RBC}; kwargs...) where {T<:Real,S<:StaticArrays.SVector,LBC,RBC}
-    der_order           = (x,y) -> x.derivative_order == y.derivative_order
-    approximation_order = (x,y) -> x.approximation_order == y.approximation_order
-    dx                  = (x,y) -> !any(x.dx .!= y.dx)
-    dimension           = (x,y) -> x.dimension == y.dimension
-    stencil_length      = (x,y) -> x.stencil_length == y.stencil_length
-    stencil_coefs       = (x,y) -> !any(!isapprox.(x.stencil_coefs,y.stencil_coefs; kwargs...))
-    boundary_point_count= (x,y) -> x.boundary_point_count == y.boundary_point_count
-    boundary_length     = (x,y) -> x.boundary_length == y.boundary_length
-    low_boundary_coefs  = (x,y) -> !any(!isapprox.(x.low_boundary_coefs,y.low_boundary_coefs; kwargs...))
-    high_boundary_coefs = (x,y) -> !any(!isapprox.(x.high_boundary_coefs,y.high_boundary_coefs; kwargs...))
-    boundary_condition  = (x,y) -> !any(!isapprox.(x.high_boundary_coefs,y.high_boundary_coefs; kwargs...))
-    t                   = (x,y) -> x.t ≈ y.t
-
-    #in order of fast to slow comparison
-    for f in (der_order,approximation_order,dimension,t,boundary_point_count,boundary_length,dx,boundary_condition,low_boundary_coefs,high_boundary_coefs,stencil_coefs)
-        if !f(x,y)
-            return false
-        end
-    end
-    return true
-end
-
-isapprox(x::FiniteDifference{T,S,LBC,RBC},y::DerivativeOperator{T,S,LBC,RBC}; kwargs...) where {T<:Real,S<:StaticArrays.SVector,LBC,RBC} = isapprox(y,x; kwargs...)
-=#