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
521 stars 198 forks source link

Out of place `linsolve=KrylovJL_GMRES` is broken #2197

Closed oscardssmith closed 2 weeks ago

oscardssmith commented 1 month ago
julia> using LinearSolve, OrdinaryDiffEq
julia> function rober(u, p, t)
           y₁, y₂, y₃ = u
           k₁, k₂, k₃ = p
           [-k₁ * y₁ + k₃ * y₂ * y₃,
             k₁ * y₁ - k₃ * y₂ * y₃ - k₂ * y₂^2,
             k₂ * y₂^2]
       end
rober (generic function with 1 method)

julia> prob = ODEProblem(rober, [1.0,0.0,0.0],(0.0,1e3),(0.04,3e7,1e4));

julia> solve(prob, FBDF(linsolve=KrylovJL_GMRES()))
ERROR: MethodError: Cannot `convert` an object of type LinearAlgebra.LU{Float64, Matrix{Float64}, Vector{Int64}} to an object of type OrdinaryDiffEq.WOperator{false, Float64, LinearAlgebra.UniformScaling{Bool}, Float64, Matrix{Float64}, Vector{Float64}, Matrix{Float64}, FunctionOperator{true, true, false, Float64, SparseDiffTools.FwdModeAutoDiffVecProd{SparseDiffTools.JacFunctionWrapper{false, true, 1, OrdinaryDiffEq.var"#790#793", Vector{Float64}, Tuple{Float64, Float64, Float64}, Float64}, Vector{Float64}, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, typeof(SparseDiffTools.auto_jacvec), typeof(SparseDiffTools.auto_jacvec!)}, Nothing, Nothing, Nothing, @NamedTuple{islinear::Bool, isconvertible::Bool, isconstant::Bool, opnorm::Nothing, issymmetric::Bool, ishermitian::Bool, isposdef::Bool, isinplace::Bool, outofplace::Bool, has_mul5::Bool, ifcache::Bool, T::DataType, batch::Bool, size::Tuple{Int64, Int64}, sizes::Tuple{Tuple{Int64}, Tuple{Int64}}, accepted_kwargs::Tuple{}, kwargs::Dict{Symbol, Any}}, Tuple{Float64, Float64, Float64}, Float64, Tuple{Vector{Float64}, Vector{Float64}}, Float64, Float64}}

Closest candidates are:
  convert(::Type{T}, ::T) where T
   @ Base Base.jl:84

Stacktrace:
  [1] setproperty!(x::OrdinaryDiffEq.NLNewtonConstantCache{…}, f::Symbol, v::LinearAlgebra.LU{…})
    @ Base ./Base.jl:40
  [2] update_W!(nlsolver::OrdinaryDiffEq.NLSolver{…}, integrator::OrdinaryDiffEq.ODEIntegrator{…}, cache::OrdinaryDiffEq.FBDFConstantCache{…}, dtgamma::Float64, repeat_step::Bool, newJW::Nothing)
    @ OrdinaryDiffEq ~/.julia/dev/OrdinaryDiffEq/src/derivative_utils.jl:837
  [3] update_W!(nlsolver::OrdinaryDiffEq.NLSolver{…}, integrator::OrdinaryDiffEq.ODEIntegrator{…}, cache::OrdinaryDiffEq.FBDFConstantCache{…}, dtgamma::Float64, repeat_step::Bool)
    @ OrdinaryDiffEq ~/.julia/dev/OrdinaryDiffEq/src/derivative_utils.jl:827
  [4] nlsolve!(nlsolver::OrdinaryDiffEq.NLSolver{…}, integrator::OrdinaryDiffEq.ODEIntegrator{…}, cache::OrdinaryDiffEq.FBDFConstantCache{…}, repeat_step::Bool)
    @ OrdinaryDiffEq ~/.julia/dev/OrdinaryDiffEq/src/nlsolve/nlsolve.jl:27
  [5] perform_step!(integrator::OrdinaryDiffEq.ODEIntegrator{…}, cache::OrdinaryDiffEq.FBDFConstantCache{…}, repeat_step::Bool)
    @ OrdinaryDiffEq ~/.julia/dev/OrdinaryDiffEq/src/perform_step/bdf_perform_step.jl:1136
  [6] perform_step!
    @ ~/.julia/dev/OrdinaryDiffEq/src/perform_step/bdf_perform_step.jl:1076 [inlined]
  [7] solve!(integrator::OrdinaryDiffEq.ODEIntegrator{…})
    @ OrdinaryDiffEq ~/.julia/dev/OrdinaryDiffEq/src/solve.jl:551
  [8] __solve(::ODEProblem{…}, ::FBDF{…}; kwargs::@Kwargs{})
    @ OrdinaryDiffEq ~/.julia/dev/OrdinaryDiffEq/src/solve.jl:7
  [9] __solve
    @ ~/.julia/dev/OrdinaryDiffEq/src/solve.jl:1 [inlined]
 [10] #solve_call#44
    @ ~/.julia/dev/DiffEqBase/src/solve.jl:612 [inlined]
 [11] solve_call(_prob::ODEProblem{…}, args::FBDF{…})
    @ DiffEqBase ~/.julia/dev/DiffEqBase/src/solve.jl:569
 [12] solve_up(prob::ODEProblem{…}, sensealg::Nothing, u0::Vector{…}, p::Tuple{…}, args::FBDF{…}; kwargs::@Kwargs{})
    @ DiffEqBase ~/.julia/dev/DiffEqBase/src/solve.jl:1080
 [13] solve_up
    @ ~/.julia/dev/DiffEqBase/src/solve.jl:1066 [inlined]
 [14] solve(prob::ODEProblem{…}, args::FBDF{…}; sensealg::Nothing, u0::Nothing, p::Nothing, wrap::Val{…}, kwargs::@Kwargs{})
    @ DiffEqBase ~/.julia/dev/DiffEqBase/src/solve.jl:1003
 [15] solve(prob::ODEProblem{…}, args::FBDF{…})
    @ DiffEqBase ~/.julia/dev/DiffEqBase/src/solve.jl:993
 [16] top-level scope
    @ REPL[23]:1
Some type information was truncated. Use `show(err)` to see complete types.

This problem exists for all KrylongJL methods, but succeeds for in place problems.

prbzrg commented 3 weeks ago

Because of its usage in default algorithm, I'm getting it too https://github.com/SciML/OrdinaryDiffEq.jl/blob/6e9cf4055f8fc53f88af8a97df034d50927dd869/src/composite_algs.jl#L114

ChrisRackauckas commented 3 weeks ago

Yeah now a lot of people might run into this. Try to fix it early next week, or we can remove it from the default alg when out of place if we really have to.