search

SciML / OrdinaryDiffEq.jl

High performance ordinary differential equation (ODE) and differential-algebraic equation (DAE) solvers, including neural ordinary differential equations (neural ODEs) and scientific machine learning (SciML)
https://diffeq.sciml.ai/latest/
Other
562 stars 211 forks source link

Issues with KLUFactorization() and sparse matrices with SplitODEFunction #2028

Open DanielVandH opened 1 year ago

DanielVandH commented 1 year ago

There are some issues when using a SplitODEFunction with KLUFactorization(). If using a MatrixOperator:

 using OrdinaryDiffEq, LinearSolve, SparseArrays  
 f2 = (du, u, p, t) -> u[2] 
 A = MatrixOperator([1.0 1.0; 1.0 1.0])
 f = SplitFunction(A, f2)
 prob = SplitODEProblem(f, [1.0, 1.0], (0.0, 1.0))
 solve(prob, TRBDF2(linsolve=KLUFactorization()), saveat=0.1)
ERROR: type Nothing has no field f
Stacktrace:
  [1] setproperty!(x::Nothing, f::Symbol, v::Function)
    @ Base .\Base.jl:39
  [2] calc_J!(J::Matrix{…}, integrator::OrdinaryDiffEq.ODEIntegrator{…}, cache::OrdinaryDiffEq.NLNewtonCache{…}, next_step::Bool)
    @ OrdinaryDiffEq C:\Users\User\.julia\packages\OrdinaryDiffEq\NbEDr\src\derivative_utils.jl:138
  [3] calc_W!(W::Matrix{…}, integrator::OrdinaryDiffEq.ODEIntegrator{…}, nlsolver::OrdinaryDiffEq.NLSolver{…}, cache::OrdinaryDiffEq.TRBDF2Cache{…}, dtgamma::Float64, repeat_step::Bool, W_transform::Bool, newJW::Nothing)
    @ OrdinaryDiffEq C:\Users\User\.julia\packages\OrdinaryDiffEq\NbEDr\src\derivative_utils.jl:703
  [4] update_W!(nlsolver::OrdinaryDiffEq.AbstractNLSolver, integrator::SciMLBase.DEIntegrator{…}, cache::Any, dtgamma::Any, repeat_step::Bool, newJW::Any)  
    @ OrdinaryDiffEq C:\Users\User\.julia\packages\OrdinaryDiffEq\NbEDr\src\derivative_utils.jl:813 [inlined]
  [5] update_W!(nlsolver::OrdinaryDiffEq.AbstractNLSolver, integrator::SciMLBase.DEIntegrator{…}, cache::Any, dtgamma::Any, repeat_step::Bool, newJW::Any)  
    @ OrdinaryDiffEq C:\Users\User\.julia\packages\OrdinaryDiffEq\NbEDr\src\derivative_utils.jl:812 [inlined]
  [6] nlsolve!(nlsolver::OrdinaryDiffEq.NLSolver{…}, integrator::OrdinaryDiffEq.ODEIntegrator{…}, cache::OrdinaryDiffEq.TRBDF2Cache{…}, repeat_step::Bool)
    @ OrdinaryDiffEq C:\Users\User\.julia\packages\OrdinaryDiffEq\NbEDr\src\nlsolve\nlsolve.jl:27
  [7] perform_step!(integrator::OrdinaryDiffEq.ODEIntegrator{…}, cache::OrdinaryDiffEq.TRBDF2Cache{…}, repeat_step::Bool)
    @ OrdinaryDiffEq C:\Users\User\.julia\packages\OrdinaryDiffEq\NbEDr\src\perform_step\sdirk_perform_step.jl:490
  [8] perform_step!(integrator::Any, cache::OrdinaryDiffEq.TRBDF2Cache{<:Array}, repeat_step::Any)
    @ OrdinaryDiffEq C:\Users\User\.julia\packages\OrdinaryDiffEq\NbEDr\src\perform_step\sdirk_perform_step.jl:467 [inlined]
  [9] solve!(integrator::OrdinaryDiffEq.ODEIntegrator{…})
    @ OrdinaryDiffEq C:\Users\User\.julia\packages\OrdinaryDiffEq\NbEDr\src\solve.jl:526
 [10] __solve(::ODEProblem{…}, ::TRBDF2{…}; kwargs::@Kwargs{…})
    @ OrdinaryDiffEq C:\Users\User\.julia\packages\OrdinaryDiffEq\NbEDr\src\solve.jl:6
 [11] __solve
    @ OrdinaryDiffEq C:\Users\User\.julia\packages\OrdinaryDiffEq\NbEDr\src\solve.jl:1 [inlined]
 [12] #solve_call#34
    @ DiffEqBase C:\Users\User\.julia\packages\DiffEqBase\d8Zv9\src\solve.jl:539 [inlined]
 [13] solve_up(prob::ODEProblem{…}, sensealg::Nothing, u0::Vector{…}, p::SciMLBase.NullParameters, args::TRBDF2{…}; kwargs::@Kwargs{…})
    @ DiffEqBase C:\Users\User\.julia\packages\DiffEqBase\d8Zv9\src\solve.jl:1000
 [14] solve_up
    @ DiffEqBase C:\Users\User\.julia\packages\DiffEqBase\d8Zv9\src\solve.jl:973 [inlined]
 [15] solve(prob::ODEProblem{…}, args::TRBDF2{…}; sensealg::Nothing, u0::Nothing, p::Nothing, wrap::Val{…}, kwargs::@Kwargs{…})
    @ DiffEqBase C:\Users\User\.julia\packages\DiffEqBase\d8Zv9\src\solve.jl:910
 [16] top-level scope
    @ c:\Users\User\.julia\dev\FiniteVolumeMethod\docs\src\literate_wyos\semilinear_equations.jl:295
Some type information was truncated. Use `show(err)` to see complete types.

If just using normal functions:

 using OrdinaryDiffEq, LinearSolve, SparseArrays  
 f1 = (du, u, p, t) -> u[1]
 f2 = (du, u, p, t) -> u[2] 
 f = SplitFunction(f1, f2, jac_prototype = sparse([1.0 1.0; 1.0 1.0]))
 prob = SplitODEProblem(f, [1.0, 1.0], (0.0, 1.0))
 solve(prob, TRBDF2(linsolve=KLUFactorization()), saveat=0.1)
ERROR: MethodError: no method matching getcolptr(::Matrix{Float64})

Closest candidates are:
  getcolptr(::SubArray{Tv, 2, <:SparseArrays.AbstractSparseMatrixCSC{Tv, Ti}, Tuple{Base.Slice{Base.OneTo{Int64}}, I}} where {Tv, Ti, I<:AbstractUnitRange})
   @ SparseArrays C:\Users\User\.julia\juliaup\julia-1.10.0-beta2+0.x64.w64.mingw32\share\julia\stdlib\v1.10\SparseArrays\src\sparsematrix.jl:186
  getcolptr(::Union{SparseArrays.FixedSparseCSC, SparseMatrixCSC})
   @ SparseArrays C:\Users\User\.julia\juliaup\julia-1.10.0-beta2+0.x64.w64.mingw32\share\julia\stdlib\v1.10\SparseArrays\src\sparsematrix.jl:185

Stacktrace:
  [1] solve!(cache::LinearSolve.LinearCache{…}, alg::KLUFactorization; kwargs::@Kwargs{…})
    @ LinearSolve C:\Users\User\.julia\packages\LinearSolve\19ohw\src\factorization.jl:829
  [2] solve!(::LinearSolve.LinearCache{…}; kwargs::@Kwargs{…})
    @ LinearSolve C:\Users\User\.julia\packages\LinearSolve\19ohw\src\common.jl:197
  [3] dolinsolve(integrator::OrdinaryDiffEq.ODEIntegrator{…}, linsolve::LinearSolve.LinearCache{…}; A::Matrix{…}, linu::Vector{…}, b::Vector{…}, du::Nothing, u::Nothing, p::Nothing, t::Nothing, weight::Nothing, solverdata::Nothing, reltol::Float64)
    @ OrdinaryDiffEq C:\Users\User\.julia\packages\OrdinaryDiffEq\NbEDr\src\misc_utils.jl:107
  [4] compute_step!(nlsolver::OrdinaryDiffEq.NLSolver{…}, integrator::OrdinaryDiffEq.ODEIntegrator{…}, γW::Float64)
    @ OrdinaryDiffEq C:\Users\User\.julia\packages\OrdinaryDiffEq\NbEDr\src\nlsolve\newton.jl:193
  [5] nlsolve!(nlsolver::OrdinaryDiffEq.NLSolver{…}, integrator::OrdinaryDiffEq.ODEIntegrator{…}, cache::OrdinaryDiffEq.TRBDF2Cache{…}, repeat_step::Bool)
    @ OrdinaryDiffEq C:\Users\User\.julia\packages\OrdinaryDiffEq\NbEDr\src\nlsolve\nlsolve.jl:52
  [6] perform_step!(integrator::OrdinaryDiffEq.ODEIntegrator{…}, cache::OrdinaryDiffEq.TRBDF2Cache{…}, repeat_step::Bool)
    @ OrdinaryDiffEq C:\Users\User\.julia\packages\OrdinaryDiffEq\NbEDr\src\perform_step\sdirk_perform_step.jl:490
  [7] perform_step!(integrator::Any, cache::OrdinaryDiffEq.TRBDF2Cache{<:Array}, repeat_step::Any)
    @ OrdinaryDiffEq C:\Users\User\.julia\packages\OrdinaryDiffEq\NbEDr\src\perform_step\sdirk_perform_step.jl:467 [inlined]
  [8] solve!(integrator::OrdinaryDiffEq.ODEIntegrator{…})
    @ OrdinaryDiffEq C:\Users\User\.julia\packages\OrdinaryDiffEq\NbEDr\src\solve.jl:526
  [9] __solve(::ODEProblem{…}, ::TRBDF2{…}; kwargs::@Kwargs{…})
    @ OrdinaryDiffEq C:\Users\User\.julia\packages\OrdinaryDiffEq\NbEDr\src\solve.jl:6
 [10] __solve
    @ OrdinaryDiffEq C:\Users\User\.julia\packages\OrdinaryDiffEq\NbEDr\src\solve.jl:1 [inlined]
 [11] #solve_call#34
    @ DiffEqBase C:\Users\User\.julia\packages\DiffEqBase\d8Zv9\src\solve.jl:539 [inlined]
 [12] solve_up(prob::ODEProblem{…}, sensealg::Nothing, u0::Vector{…}, p::SciMLBase.NullParameters, args::TRBDF2{…}; kwargs::@Kwargs{…})
    @ DiffEqBase C:\Users\User\.julia\packages\DiffEqBase\d8Zv9\src\solve.jl:1000
 [13] solve_up
    @ DiffEqBase C:\Users\User\.julia\packages\DiffEqBase\d8Zv9\src\solve.jl:973 [inlined]
 [14] solve(prob::ODEProblem{…}, args::TRBDF2{…}; sensealg::Nothing, u0::Nothing, p::Nothing, wrap::Val{…}, kwargs::@Kwargs{…})
    @ DiffEqBase C:\Users\User\.julia\packages\DiffEqBase\d8Zv9\src\solve.jl:910
 [15] top-level scope
    @ c:\Users\User\.julia\dev\FiniteVolumeMethod\docs\src\literate_wyos\semilinear_equations.jl:295
Some type information was truncated. Use `show(err)` to see complete types.

I think the issue in this case is that jac_prototype doesn't get correctly passed into prob?

julia> f.jac_prototype
2×2 SparseMatrixCSC{Float64, Int64} with 4 stored entries:
 1.0  1.0
 1.0  1.0

julia> prob.f.jac_prototype

julia> prob.f.f1.jac_prototype

julia> prob.f.f2.jac_prototype
ChrisRackauckas commented 1 year ago

Indeed, solving with a non-split method throws a warning if the solver isn't set:

┌ Warning: Split ODE problem do not work with factorization linear solvers. Bug detailed in https://github.com/SciML/OrdinaryDiffEq.jl/pull/1643. Defaulting to linsolve=KrylovJL()
└ @ OrdinaryDiffEq C:\Users\accou\.julia\packages\OrdinaryDiffEq\NbEDr\src\alg_utils.jl:248

referencing https://github.com/SciML/OrdinaryDiffEq.jl/pull/1643. Some of the PDE benchmarks need this fixed.

CC: @vpuri3