Solving linear system with summed DiffEqOperators

General Usage
1

I’m trying to solve an elliptic differential equation that has multiple derivative terms and wanted to use DiffEqOperators.jl to build the matrices. Is it possible to sum DiffEqOperators and then solve using backslash? When I try this, I get an error. Here’s an example, based on https://discourse.julialang.org/t/how-to-discretize-and-solve-a-3d-poisson-equation-with-dirichlet-boundary-conditions-owing-to-diffeqoperators/40796:

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
1 Like
2

u = SparseMatrixCSC((0.5*Δ+0.5*Δ)*bc) \ f Hmm that’s missing a concretization. Can you open an issue?

1 Like
3

I’m not sure what the solution is, but FWIW, this happens because the ghost operators have their own backslash defined here, which I think is necessary since the ghost operator is not strictly a linear operator but an affine one AFAIU.

Maybe you need to add a *(operator combination, bc), similar to this for operator composition?

4

Yeah exactly. It seems like it’s just a missing method.

5

Submitted as
https://github.com/SciML/DiffEqOperators.jl/issues/391

1 Like
6

This does work:

A1, b1 = sparse(0.5*Δ*bc)
A2, b2 = sparse(0.5*Δ*bc)
As = A1 + A2
u = As \ f

so I can combine linear operators this way, though it’s less elegant than adding a method for DiffEqOperatorCombination as @briochemc suggested.

Separately but related, calling sparse on a DiffEqOperator is much slower and more memory intensive than building the same SparseMatrixCSC from three vectors. Here’s some example text, but I can open a separate issue/thread for this if it’s too off-topic:

N = 10001
Rᵢ = 20.
r = collect(range(0.,Rᵢ,length=N))
Δr = r[2]

function by_hand()

    R = zeros(3N-2)
    C = zeros(3N-2)
    V = zeros(3N-2)
    for i in 2:(N-1)

        # diagonal - 1
        j = 3*(i-1)
        R[j], C[j], V[j] = i, i-1, 1.0/Δr^2

        # diagonal
        j = 3*(i-1) + 1
        R[j], C[j], V[j] = i, i, -2.0/Δr^2

        # diagonal - 1
        j = 3*(i-1) + 2
        R[j], C[j], V[j] = i, i+1, 1.0/Δr^2

    end
    # Zero derivative at boundaries (Neumann)
    j= 1
    R[j], C[j], V[j] = 1, 1, 1.
    j = 2
    R[j], C[j], V[j] = 1, 2, -1.

    j = 3N - 2
    R[j], C[j], V[j] = N, N, 1.
    j = 3N - 3
    R[j], C[j], V[j] = N, N-1, -1.

    A = sparse(R,C,V)
end
function by_DEO()
    Δ = CenteredDifference(2,2,Δr,N)
    bc = Neumann0BC(Δr,1)
    A,b = sparse(Δ*bc)
    return A
end

println("Create by hand")
@btime Abh = by_hand();
println("Create by DiffEqOperators")
@btime Abd = by_DEO();

Output:

Create by hand
  6.089 ms (185603 allocations: 5.27 MiB)
Create by DiffEqOperators
  5.218 s (257528 allocations: 6.01 GiB)

I would have thought sparse on a DerivativeOperator or GhostDerivativeOperator would have done something like my by_hand function. Any thoughts?

7

At it to the issue. It should be doing it like you did by hand, it may be falling back to an iteration on the indices method.

8

I opened a separate issue on git

https://github.com/SciML/DiffEqOperators.jl/issues/392