Merge pull request #511 from xtalax/completedifferences

Complete Differences

Author: ChrisRackauckas

src/DiffEqOperators.jl

@@ -73,7 +73,7 @@ end
 export MatrixFreeOperator
 export AnalyticalJacVecOperator, JacVecOperator, getops
 export AbstractDerivativeOperator, DerivativeOperator,
-       CenteredDifference, UpwindDifference, nonlinear_diffusion, nonlinear_diffusion!,
+       CenteredDifference, UpwindDifference, CompleteCenteredDifference, CompleteUpwindDifference, CompleteHalfCenteredDifference, nonlinear_diffusion, nonlinear_diffusion!,
        GradientOperator, Gradient, CurlOperator, Curl, DivergenceOperator, Divergence
 export DirichletBC, Dirichlet0BC, NeumannBC, Neumann0BC, RobinBC, GeneralBC, MultiDimBC, PeriodicBC,
        MultiDimDirectionalBC, ComposedMultiDimBC

src/composite_operators.jl

@@ -29,19 +29,19 @@ struct DiffEqOperatorCombination{T,O<:Tuple{Vararg{AbstractDiffEqLinearOperator{
         new{T,typeof(ops),typeof(cache)}(ops, cache)
     end
 end
-+(ops::AbstractDiffEqLinearOperator...) = DiffEqOperatorCombination(ops)
-+(op::AbstractDiffEqLinearOperator, A::AbstractMatrix) = op + DiffEqArrayOperator(A)
-+(op::AbstractDiffEqLinearOperator{T}, α::UniformScaling) where T = op + UniformScaling(T(α.λ))(size(op,1))
-+(A::AbstractMatrix, op::AbstractDiffEqLinearOperator) = op + A
-+(α::UniformScaling, op::AbstractDiffEqLinearOperator) = op + α
-+(L1::DiffEqOperatorCombination, L2::AbstractDiffEqLinearOperator) = DiffEqOperatorCombination((L1.ops..., L2))
-+(L1::AbstractDiffEqLinearOperator, L2::DiffEqOperatorCombination) = DiffEqOperatorCombination((L1, L2.ops...))
-+(L1::DiffEqOperatorCombination, L2::DiffEqOperatorCombination) = DiffEqOperatorCombination((L1.ops..., L2.ops...))
--(L1::AbstractDiffEqLinearOperator, L2::AbstractDiffEqLinearOperator) = L1 + (-L2)
--(L::AbstractDiffEqLinearOperator, A::AbstractMatrix) = L + (-A)
--(A::AbstractMatrix, L::AbstractDiffEqLinearOperator) = A + (-L)
--(L::AbstractDiffEqLinearOperator, α::UniformScaling) = L + (-α)
--(α::UniformScaling, L::AbstractDiffEqLinearOperator) = α + (-L)
+Base.:+(ops::AbstractDiffEqLinearOperator...) = DiffEqOperatorCombination(ops)
+Base.:+(op::AbstractDiffEqLinearOperator, A::AbstractMatrix) = op + DiffEqArrayOperator(A)
+Base.:+(op::AbstractDiffEqLinearOperator{T}, α::UniformScaling) where T = op + UniformScaling(T(α.λ))(size(op,1))
+Base.:+(A::AbstractMatrix, op::AbstractDiffEqLinearOperator) = op + A
+Base.:+(α::UniformScaling, op::AbstractDiffEqLinearOperator) = op + α
+Base.:+(L1::DiffEqOperatorCombination, L2::AbstractDiffEqLinearOperator) = DiffEqOperatorCombination((L1.ops..., L2))
+Base.:+(L1::AbstractDiffEqLinearOperator, L2::DiffEqOperatorCombination) = DiffEqOperatorCombination((L1, L2.ops...))
+Base.:+(L1::DiffEqOperatorCombination, L2::DiffEqOperatorCombination) = DiffEqOperatorCombination((L1.ops..., L2.ops...))
+Base.:-(L1::AbstractDiffEqLinearOperator, L2::AbstractDiffEqLinearOperator) = L1 + (-L2)
+Base.:-(L::AbstractDiffEqLinearOperator, A::AbstractMatrix) = L + (-A)
+Base.:-(A::AbstractMatrix, L::AbstractDiffEqLinearOperator) = A + (-L)
+Base.:-(L::AbstractDiffEqLinearOperator, α::UniformScaling) = L + (-α)
+Base.:-(α::UniformScaling, L::AbstractDiffEqLinearOperator) = α + (-L)
 getops(L::DiffEqOperatorCombination) = L.ops
 Matrix(L::DiffEqOperatorCombination) = sum(Matrix, L.ops)
 convert(::Type{AbstractMatrix}, L::DiffEqOperatorCombination) =
@@ -51,8 +51,8 @@ SparseArrays.sparse(L::DiffEqOperatorCombination) = sum(sparse1, L.ops)
 size(L::DiffEqOperatorCombination, args...) = size(L.ops[1], args...)
 getindex(L::DiffEqOperatorCombination, i::Int) = sum(op -> op[i], L.ops)
 getindex(L::DiffEqOperatorCombination, I::Vararg{Int, N}) where {N} = sum(op -> op[I...], L.ops)
-*(L::DiffEqOperatorCombination, x::AbstractArray) = sum(op -> op * x, L.ops)
-*(x::AbstractArray, L::DiffEqOperatorCombination) = sum(op -> x * op, L.ops)
+Base.:*(L::DiffEqOperatorCombination, x::AbstractArray) = sum(op -> op * x, L.ops)
+Base.:*(x::AbstractArray, L::DiffEqOperatorCombination) = sum(op -> x * op, L.ops)
 /(L::DiffEqOperatorCombination, x::AbstractArray) = sum(op -> op / x, L.ops)
 \(x::AbstractArray, L::DiffEqOperatorCombination) = sum(op -> x \ op, L.ops)
 function mul!(y::AbstractVector, L::DiffEqOperatorCombination, b::AbstractVector)
@@ -88,13 +88,20 @@ struct DiffEqOperatorComposition{T,O<:Tuple{Vararg{AbstractDiffEqLinearOperator{
     new{T,typeof(ops),typeof(caches)}(ops, caches)
   end
 end
-*(ops::AbstractDiffEqLinearOperator...) = DiffEqOperatorComposition(reverse(ops))
+# this is needed to not break dispatch in MethodOfLines
+function Base.:*(ops::AbstractDiffEqLinearOperator...) 
+  try 
+    return DiffEqOperatorComposition(reverse(ops))
+  catch e
+    return 1
+  end
+end
 ∘(L1::AbstractDiffEqLinearOperator, L2::AbstractDiffEqLinearOperator) = DiffEqOperatorComposition((L2, L1))
-*(L1::DiffEqOperatorComposition, L2::AbstractDiffEqLinearOperator) = DiffEqOperatorComposition((L2, L1.ops...))
+Base.:*(L1::DiffEqOperatorComposition, L2::AbstractDiffEqLinearOperator) = DiffEqOperatorComposition((L2, L1.ops...))
 ∘(L1::DiffEqOperatorComposition, L2::AbstractDiffEqLinearOperator) = DiffEqOperatorComposition((L2, L1.ops...))
-*(L1::AbstractDiffEqLinearOperator, L2::DiffEqOperatorComposition) = DiffEqOperatorComposition((L2.ops..., L1))
+Base.:*(L1::AbstractDiffEqLinearOperator, L2::DiffEqOperatorComposition) = DiffEqOperatorComposition((L2.ops..., L1))
 ∘(L1::AbstractDiffEqLinearOperator, L2::DiffEqOperatorComposition) = DiffEqOperatorComposition((L2.ops..., L1))
-*(L1::DiffEqOperatorComposition, L2::DiffEqOperatorComposition) = DiffEqOperatorComposition((L2.ops..., L1.ops...))
+Base.:*(L1::DiffEqOperatorComposition, L2::DiffEqOperatorComposition) = DiffEqOperatorComposition((L2.ops..., L1.ops...))
 ∘(L1::DiffEqOperatorComposition, L2::DiffEqOperatorComposition) = DiffEqOperatorComposition((L2.ops..., L1.ops...))
 getops(L::DiffEqOperatorComposition) = L.ops
 Matrix(L::DiffEqOperatorComposition) = prod(Matrix, reverse(L.ops))
@@ -105,8 +112,8 @@ SparseArrays.sparse(L::DiffEqOperatorComposition) = prod(sparse1, reverse(L.ops)
 size(L::DiffEqOperatorComposition) = (size(L.ops[end], 1), size(L.ops[1], 2))
 size(L::DiffEqOperatorComposition, m::Integer) = size(L)[m]
 opnorm(L::DiffEqOperatorComposition) = prod(opnorm, L.ops)
-*(L::DiffEqOperatorComposition, x::AbstractArray) = foldl((acc, op) -> op*acc, L.ops; init=x)
-*(x::AbstractArray, L::DiffEqOperatorComposition) = foldl((acc, op) -> acc*op, reverse(L.ops); init=x)
+Base.:*(L::DiffEqOperatorComposition, x::AbstractArray) = foldl((acc, op) -> op*acc, L.ops; init=x)
+Base.:*(x::AbstractArray, L::DiffEqOperatorComposition) = foldl((acc, op) -> acc*op, reverse(L.ops); init=x)
 /(L::DiffEqOperatorComposition, x::AbstractArray) = foldl((acc, op) -> op/acc, L.ops; init=x)
 /(x::AbstractArray, L::DiffEqOperatorComposition) = foldl((acc, op) -> acc/op, L.ops; init=x)
 \(L::DiffEqOperatorComposition, x::AbstractArray) = foldl((acc, op) -> op\acc, reverse(L.ops); init=x)

src/derivative_operators/derivative_operator.jl

@@ -105,6 +105,8 @@ function CenteredDifference{N}(derivative_order::Int,
         )
 end
 
+
+
 function generate_coordinates(i::Int, stencil_x, dummy_x, dx::AbstractVector{T}) where T <: Real
     len = length(stencil_x)
     stencil_x .= stencil_x.*zero(T)
@@ -161,6 +163,98 @@ function CenteredDifference{N}(derivative_order::Int,
         )
 end
 
+
+# TODO: Set weight to zero if its below a certain threshold
+"""
+A helper function to compute the coefficients of a derivative operator including the boundary coefficients in the centered scheme.
+"""
+function CompleteCenteredDifference(derivative_order::Int,
+    approximation_order::Int, dx::T) where {T<:Real,N}
+    @assert approximation_order>1 "approximation_order must be greater than 1."
+    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)
+    left_boundary_x         = 0:(boundary_stencil_length-1)
+    right_boundary_x        = reverse(-boundary_stencil_length+1:0)
+
+    boundary_point_count    = div(stencil_length,2) # -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  # unused
+    L_boundary_deriv_spots  = left_boundary_x[1:div(stencil_length,2)]
+
+    stencil_coefs           = convert(SVector{stencil_length, T}, (1/dx^derivative_order) * calculate_weights(derivative_order, zero(T), dummy_x))
+    _low_boundary_coefs     = SVector{boundary_stencil_length, T}[convert(SVector{boundary_stencil_length, T}, (1/dx^derivative_order) * calculate_weights(derivative_order, oneunit(T)*x0, left_boundary_x)) for x0 in L_boundary_deriv_spots]
+    low_boundary_coefs      = convert(SVector{boundary_point_count},vcat(_low_boundary_coefs))
+
+    # _high_boundary_coefs    = SVector{boundary_stencil_length, T}[convert(SVector{boundary_stencil_length, T}, (1/dx^derivative_order) * calculate_weights(derivative_order, oneunit(T)*x0, reverse(right_boundary_x))) for x0 in R_boundary_deriv_spots]
+    high_boundary_coefs      = convert(SVector{boundary_point_count},[reverse(v)*(-1)^derivative_order for v in _low_boundary_coefs])
+
+    offside = 0
+
+    coefficients            = nothing
+
+
+    DerivativeOperator{T,Nothing,false,T,typeof(stencil_coefs),
+    typeof(low_boundary_coefs),typeof(high_boundary_coefs),typeof(coefficients),
+    Nothing}(
+    derivative_order, approximation_order, dx, 1, stencil_length,
+    stencil_coefs,
+    boundary_stencil_length,
+    boundary_point_count,
+    low_boundary_coefs,
+    high_boundary_coefs,offside,coefficients,nothing
+    )
+end
+
+
+"""
+A helper function to compute the coefficients of a derivative operator including the boundary coefficients in the half offset centered scheme. See table 2 in https://web.njit.edu/~jiang/math712/fornberg.pdf
+"""
+function CompleteHalfCenteredDifference(derivative_order::Int,
+    approximation_order::Int, dx::T) where {T<:Real,N}
+    @assert approximation_order>1 "approximation_order must be greater than 1."
+    centered_stencil_length          = approximation_order + 2*Int(floor(derivative_order/2)) + (approximation_order%2)
+    boundary_stencil_length = derivative_order + approximation_order
+    endpoint = div(centered_stencil_length, 2)
+    dummy_x                 = 1-endpoint : endpoint
+    left_boundary_x         = 1:(boundary_stencil_length)
+    right_boundary_x        = reverse(-boundary_stencil_length:-1)
+
+    boundary_point_count    = div(centered_stencil_length,2) # -1 due to the ghost point
+    # ? Is fornberg valid when taking an x0 outside of the stencil i.e at the boundary?
+    xoffset = range(1.5, length=boundary_point_count, step = 1.0)
+
+    # 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  # unused
+    L_boundary_deriv_spots  = xoffset[1:div(centered_stencil_length,2)]
+
+    stencil_coefs           = convert(SVector{centered_stencil_length, T}, (1/dx^derivative_order) * calculate_weights(derivative_order, convert(T, 0.5), dummy_x))
+    # For each boundary point, for each tappoint in the half offset central difference stencil, we need to calculate the coefficients for the stencil.
+
+
+    _low_boundary_coefs     = [convert(SVector{boundary_stencil_length, T}, (1/dx^derivative_order) * calculate_weights(derivative_order, offset, left_boundary_x)) for offset in L_boundary_deriv_spots]
+    low_boundary_coefs     = convert(SVector{boundary_point_count}, _low_boundary_coefs)
+
+    _high_boundary_coefs = [reverse(stencil)*(-1)^derivative_order for stencil in low_boundary_coefs]
+    high_boundary_coefs = convert(SVector{boundary_point_count}, _high_boundary_coefs)
+    # _high_boundary_coefs    = SVector{boundary_stencil_length, T}[convert(SVector{boundary_stencil_length, T}, (1/dx^derivative_order) * calculate_weights(derivative_order, oneunit(T)*x0, reverse(right_boundary_x))) for x0 in R_boundary_deriv_spots]
+
+    offside = 0
+
+    coefficients            = nothing
+
+
+    DerivativeOperator{T,Nothing,false,T,typeof(stencil_coefs),
+    typeof(low_boundary_coefs),typeof(high_boundary_coefs),typeof(coefficients),Nothing}(
+    derivative_order, approximation_order, dx, 1, centered_stencil_length,
+    stencil_coefs,
+    boundary_stencil_length,
+    boundary_point_count,
+    low_boundary_coefs,
+    high_boundary_coefs,offside,coefficients,nothing
+    )
+end
+
 struct UpwindDifference{N} end
 
 """
@@ -287,6 +381,50 @@ function UpwindDifference{N}(derivative_order::Int,
         )
 end
 
+"""
+A helper function to compute the coefficients of a derivative operator including the boundary coefficients in the upwind scheme.
+"""
+function CompleteUpwindDifference(derivative_order::Int,
+    approximation_order::Int, dx::T,
+    offside::Int=0) where {T<:Real,N}
+
+    @assert offside > -1 "Number of offside points should be non-negative"
+    @assert offside <= div(derivative_order + approximation_order - 1,2) "Number of offside points should not exceed the primary wind points"
+
+    stencil_length          = derivative_order + approximation_order
+    boundary_stencil_length = derivative_order + approximation_order
+    boundary_point_count    = boundary_stencil_length - 1 - offside
+
+    # TODO: Clean up the implementation here so that it is more readable and easier to extend in the future
+    dummy_x = (0.0 - offside) : stencil_length - 1.0 - offside
+    stencil_coefs = convert(SVector{stencil_length, T}, (1/dx^derivative_order) * calculate_weights(derivative_order, 0.0, dummy_x))
+
+    low_boundary_x         = 0.0:(boundary_stencil_length-1)
+    L_boundary_deriv_spots = 0.0:boundary_stencil_length - 2.0 - offside
+    _low_boundary_coefs     = SVector{boundary_stencil_length, T}[convert(SVector{boundary_stencil_length, T}, (1/dx^derivative_order) * calculate_weights(derivative_order, oneunit(T)*x0, low_boundary_x)) for x0 in L_boundary_deriv_spots]
+    low_boundary_coefs      = convert(SVector{boundary_point_count},_low_boundary_coefs)
+
+    high_boundary_x         = 0.0:-1.0:-(boundary_stencil_length-1.0)
+    R_boundary_deriv_spots = 0.0:-1.0:-(boundary_stencil_length-2.0)
+    _high_boundary_coefs     = SVector{boundary_stencil_length, T}[convert(SVector{boundary_stencil_length, T}, ((-1/dx)^derivative_order) * calculate_weights(derivative_order, oneunit(T)*x0, high_boundary_x)) for x0 in R_boundary_deriv_spots]
+    high_boundary_coefs = convert(SVector{boundary_point_count + offside},_high_boundary_coefs)
+
+    coefficients = fill!(Vector{T}(undef,len),0)
+
+
+    DerivativeOperator{T,N,true,T,typeof(stencil_coefs),
+    typeof(low_boundary_coefs),typeof(high_boundary_coefs),Vector{T},
+    typeof(coeff_func)}(
+    derivative_order, approximation_order, dx, len, stencil_length,
+    stencil_coefs,
+    boundary_stencil_length,
+    boundary_point_count,
+    low_boundary_coefs,
+    high_boundary_coefs,offside,coefficients,coeff_func
+    )
+end
+
+
 CenteredDifference(args...) = CenteredDifference{1}(args...)
 UpwindDifference(args...;kwargs...) = UpwindDifference{1}(args...;kwargs...)
 nonlinear_diffusion(args...) = nonlinear_diffusion{1}(args...)

src/derivative_operators/fornberg.jl

@@ -55,6 +55,8 @@ function calculate_weights(order::Int, x0::T, x::AbstractVector) where T<:Real
         http://epubs.siam.org/doi/pdf/10.1137/S0036144596322507 - Modified Fornberg Algorithm
     =#
     _C = C[:,end]
-    _C[div(N,2)+1] -= sum(_C)
+    if order != 0
+        _C[div(N,2)+1] -= sum(_C)
+    end
     return _C
 end

test/DerivativeOperators/derivative_operators_interface.jl

@@ -144,6 +144,34 @@ end
     @test BandedMatrix(L) ≈ correct
 end
 
+@testset "Correctness of Uniform Stencils, Complete" begin
+    weights = []
+    
+    push!(weights, ([-0.5,0,0.5], [1.,-2.,1.], [-1/2,1.,0.,-1.,1/2]))
+    push!(weights, ([1/12, -2/3,0,2/3,-1/12], [-1/12,4/3,-5/2,4/3,-1/12], [1/8,-1.,13/8,0.,-13/8,1.,-1/8]))
+
+    for d in 1:3
+        for (i,a) in enumerate([2,4])
+            D = CompleteCenteredDifference(d,a,1.0)
+
+            @test all(isapprox.(D.stencil_coefs, weights[i][d], atol=1e-10))
+        end
+    end
+end
+
+@testset "Correctness of Uniform Stencils, Complete Half" begin
+    weights = (([.5, .5], [-1/16, 9/16, 9/16, -1/16]),
+               ([-1., 1.], [1/24, -9/8, 9/8, -1/24]))
+    for (i,a) in enumerate([2,4])
+        for d in 0:1
+            D = CompleteHalfCenteredDifference(d,a,1.0)
+            @test all(D.stencil_coefs .≈ weights[d+1][i]) 
+        end
+    end
+
+end
+
+
 @testset "Correctness of Non-Uniform Stencils" begin
     x = [0., 0.08, 0.1, 0.15, 0.19, 0.26, 0.29]
     nx = length(x)