Cannot perform linear solves with DiffEqOperatorCombination

[bclyons12]

I'm trying to solve an elliptic differential equation that has multiple derivative terms and wanted to use DiffEqOperators.jl to build the matrices. It does not seem possible to sum DiffEqOperators and then solve using backslash? When I try this, I get an error:

using DiffEqOperators
Δx = 0.1
x = Δx:Δx:1-Δx # Solve only for the interior: the endpoints are known to be zero!
N = length(x)
f = sin.(2π*x)

Δ = CenteredDifference(2, 2, Δx, N)
bc = Dirichlet0BC(Float64)

# Works
u = (Δ*bc) \ f
println("This solve is successful")

#Fails
u = ((0.5*Δ+0.5*Δ)*bc) \ f
println("But this solve fails")

with output

This solve is successful
ERROR: LoadError: MethodError: Cannot `convert` an object of type 
  GhostDerivativeOperator{Float64, DerivativeOperator{Float64, 1, false, Float64, StaticArrays.SVector{3, Float64}, StaticArrays.SVector{0, StaticArrays.SVector{4, Float64}}, StaticArrays.SVector{0, StaticArrays.SVector{4, Float64}}, Float64, Nothing}, RobinBC{Float64, Vector{Float64}}} to an object of type 
  AbstractMatrix{T} where T
Closest candidates are:
  convert(::Type{AbstractMatrix{T} where T}, ::DerivativeOperator) at /Users/Lyons/.julia/packages/DiffEqOperators/DMNmH/src/derivative_operators/concretization.jl:63
  convert(::Type{AbstractMatrix{T} where T}, ::SciMLBase.DiffEqScaledOperator) at /Users/Lyons/.julia/packages/SciMLBase/XuLdB/src/operators/diffeq_operator.jl:97
  convert(::Type{AbstractMatrix{T} where T}, ::SciMLBase.DiffEqArrayOperator) at /Users/Lyons/.julia/packages/SciMLBase/XuLdB/src/operators/basic_operators.jl:102
  ...
Stacktrace:
  [1] (::DiffEqOperators.var"#117#118")(op::GhostDerivativeOperator{Float64, DerivativeOperator{Float64, 1, false, Float64, StaticArrays.SVector{3, Float64}, StaticArrays.SVector{0, StaticArrays.SVector{4, Float64}}, StaticArrays.SVector{0, StaticArrays.SVector{4, Float64}}, Float64, Nothing}, RobinBC{Float64, Vector{Float64}}})
    @ DiffEqOperators ~/.julia/packages/DiffEqOperators/DMNmH/src/composite_operators.jl:38
  [2] MappingRF
    @ ./reduce.jl:93 [inlined]
  [3] afoldl(::Base.MappingRF{DiffEqOperators.var"#117#118", Base.BottomRF{typeof(Base.add_sum)}}, ::Base._InitialValue, ::GhostDerivativeOperator{Float64, DerivativeOperator{Float64, 1, false, Float64, StaticArrays.SVector{3, Float64}, StaticArrays.SVector{0, StaticArrays.SVector{4, Float64}}, StaticArrays.SVector{0, StaticArrays.SVector{4, Float64}}, Float64, Nothing}, RobinBC{Float64, Vector{Float64}}}, ::GhostDerivativeOperator{Float64, DerivativeOperator{Float64, 1, false, Float64, StaticArrays.SVector{3, Float64}, StaticArrays.SVector{0, StaticArrays.SVector{4, Float64}}, StaticArrays.SVector{0, StaticArrays.SVector{4, Float64}}, Float64, Nothing}, RobinBC{Float64, Vector{Float64}}})
    @ Base ./operators.jl:533
  [4] _foldl_impl(op::Base.MappingRF{DiffEqOperators.var"#117#118", Base.BottomRF{typeof(Base.add_sum)}}, init::Base._InitialValue, itr::Tuple{GhostDerivativeOperator{Float64, DerivativeOperator{Float64, 1, false, Float64, StaticArrays.SVector{3, Float64}, StaticArrays.SVector{0, StaticArrays.SVector{4, Float64}}, StaticArrays.SVector{0, StaticArrays.SVector{4, Float64}}, Float64, Nothing}, RobinBC{Float64, Vector{Float64}}}, GhostDerivativeOperator{Float64, DerivativeOperator{Float64, 1, false, Float64, StaticArrays.SVector{3, Float64}, StaticArrays.SVector{0, StaticArrays.SVector{4, Float64}}, StaticArrays.SVector{0, StaticArrays.SVector{4, Float64}}, Float64, Nothing}, RobinBC{Float64, Vector{Float64}}}})
    @ Base ./tuple.jl:263
  [5] foldl_impl(op::Base.MappingRF{DiffEqOperators.var"#117#118", Base.BottomRF{typeof(Base.add_sum)}}, nt::Base._InitialValue, itr::Tuple{GhostDerivativeOperator{Float64, DerivativeOperator{Float64, 1, false, Float64, StaticArrays.SVector{3, Float64}, StaticArrays.SVector{0, StaticArrays.SVector{4, Float64}}, StaticArrays.SVector{0, StaticArrays.SVector{4, Float64}}, Float64, Nothing}, RobinBC{Float64, Vector{Float64}}}, GhostDerivativeOperator{Float64, DerivativeOperator{Float64, 1, false, Float64, StaticArrays.SVector{3, Float64}, StaticArrays.SVector{0, StaticArrays.SVector{4, Float64}}, StaticArrays.SVector{0, StaticArrays.SVector{4, Float64}}, Float64, Nothing}, RobinBC{Float64, Vector{Float64}}}})
    @ Base ./reduce.jl:48
  [6] mapfoldl_impl(f::DiffEqOperators.var"#117#118", op::typeof(Base.add_sum), nt::Base._InitialValue, itr::Tuple{GhostDerivativeOperator{Float64, DerivativeOperator{Float64, 1, false, Float64, StaticArrays.SVector{3, Float64}, StaticArrays.SVector{0, StaticArrays.SVector{4, Float64}}, StaticArrays.SVector{0, StaticArrays.SVector{4, Float64}}, Float64, Nothing}, RobinBC{Float64, Vector{Float64}}}, GhostDerivativeOperator{Float64, DerivativeOperator{Float64, 1, false, Float64, StaticArrays.SVector{3, Float64}, StaticArrays.SVector{0, StaticArrays.SVector{4, Float64}}, StaticArrays.SVector{0, StaticArrays.SVector{4, Float64}}, Float64, Nothing}, RobinBC{Float64, Vector{Float64}}}})
    @ Base ./reduce.jl:44
  [7] mapfoldl(f::Function, op::Function, itr::Tuple{GhostDerivativeOperator{Float64, DerivativeOperator{Float64, 1, false, Float64, StaticArrays.SVector{3, Float64}, StaticArrays.SVector{0, StaticArrays.SVector{4, Float64}}, StaticArrays.SVector{0, StaticArrays.SVector{4, Float64}}, Float64, Nothing}, RobinBC{Float64, Vector{Float64}}}, GhostDerivativeOperator{Float64, DerivativeOperator{Float64, 1, false, Float64, StaticArrays.SVector{3, Float64}, StaticArrays.SVector{0, StaticArrays.SVector{4, Float64}}, StaticArrays.SVector{0, StaticArrays.SVector{4, Float64}}, Float64, Nothing}, RobinBC{Float64, Vector{Float64}}}}; init::Base._InitialValue)
    @ Base ./reduce.jl:160
  [8] mapfoldl
    @ ./reduce.jl:160 [inlined]
  [9] #mapreduce#218
    @ ./reduce.jl:287 [inlined]
 [10] mapreduce
    @ ./reduce.jl:287 [inlined]
 [11] #sum#221
    @ ./reduce.jl:501 [inlined]
 [12] sum(f::Function, a::Tuple{GhostDerivativeOperator{Float64, DerivativeOperator{Float64, 1, false, Float64, StaticArrays.SVector{3, Float64}, StaticArrays.SVector{0, StaticArrays.SVector{4, Float64}}, StaticArrays.SVector{0, StaticArrays.SVector{4, Float64}}, Float64, Nothing}, RobinBC{Float64, Vector{Float64}}}, GhostDerivativeOperator{Float64, DerivativeOperator{Float64, 1, false, Float64, StaticArrays.SVector{3, Float64}, StaticArrays.SVector{0, StaticArrays.SVector{4, Float64}}, StaticArrays.SVector{0, StaticArrays.SVector{4, Float64}}, Float64, Nothing}, RobinBC{Float64, Vector{Float64}}}})
    @ Base ./reduce.jl:501
 [13] convert(#unused#::Type{AbstractMatrix{T} where T}, L::DiffEqOperators.DiffEqOperatorCombination{Float64, Tuple{GhostDerivativeOperator{Float64, DerivativeOperator{Float64, 1, false, Float64, StaticArrays.SVector{3, Float64}, StaticArrays.SVector{0, StaticArrays.SVector{4, Float64}}, StaticArrays.SVector{0, StaticArrays.SVector{4, Float64}}, Float64, Nothing}, RobinBC{Float64, Vector{Float64}}}, GhostDerivativeOperator{Float64, DerivativeOperator{Float64, 1, false, Float64, StaticArrays.SVector{3, Float64}, StaticArrays.SVector{0, StaticArrays.SVector{4, Float64}}, StaticArrays.SVector{0, StaticArrays.SVector{4, Float64}}, Float64, Nothing}, RobinBC{Float64, Vector{Float64}}}}, Vector{Float64}})
    @ DiffEqOperators ~/.julia/packages/DiffEqOperators/DMNmH/src/composite_operators.jl:37
 [14] \(L::DiffEqOperators.DiffEqOperatorCombination{Float64, Tuple{GhostDerivativeOperator{Float64, DerivativeOperator{Float64, 1, false, Float64, StaticArrays.SVector{3, Float64}, StaticArrays.SVector{0, StaticArrays.SVector{4, Float64}}, StaticArrays.SVector{0, StaticArrays.SVector{4, Float64}}, Float64, Nothing}, RobinBC{Float64, Vector{Float64}}}, GhostDerivativeOperator{Float64, DerivativeOperator{Float64, 1, false, Float64, StaticArrays.SVector{3, Float64}, StaticArrays.SVector{0, StaticArrays.SVector{4, Float64}}, StaticArrays.SVector{0, StaticArrays.SVector{4, Float64}}, Float64, Nothing}, RobinBC{Float64, Vector{Float64}}}}, Vector{Float64}}, x::Vector{Float64})
    @ SciMLBase ~/.julia/packages/SciMLBase/XuLdB/src/operators/common_defaults.jl:14

For some futher discussion, please see https://discourse.julialang.org/t/solving-linear-system-with-summed-diffeqoperators/60998

[bclyons12]

@ChrisRackauckas Here's one of the issues I mentioned in our call today

[ChrisRackauckas]

This is probably best replaced with MethodOfLines.jl

[sjkelly]

There is a strange dispatch bug here. Changing:

u = ((0.5*Δ+0.5*Δ)*bc) \ f

to:

u = (0.5*Δ*bc) \f + (0.5*Δ*bc) \f

fixes this.

[bclyons12]

@sjkelly I don't believe the \ operator can be distributed like that. (A + B) * x = S gives x = (A + B)^-1 * S, which is not x = A^-1 * S + B^-1 * S.

[sjkelly]

Woops, yes, you are correct. Something didn't seem right about that!

[xtalax]

Yes this is a known issue, nobody was able to fix this back when ghost derivs were being implemented and I only noticed recently that they just deleted the test. Again, best practice these days is to use MethodOfLines

[ChrisRackauckas]

Yes, this is rather difficult to solve in the DiffEqOperators.jl form and is one of the reasons why the library was deprecated in favor of MethodOfLines.jl.