search

SciML / DiffEqOperators.jl

Linear operators for discretizations of differential equations and scientific machine learning (SciML)
https://docs.sciml.ai/DiffEqOperators/stable/
Other
284 stars 74 forks source link

SplitODEProblem broken for MatrixFreeOperator #520

Closed AshtonSBradley closed 2 years ago

AshtonSBradley commented 2 years ago

Does MatrixFreeOperator work with SplitODEProblem? 

using DifferentialEquations, DiffEqOperators
b = randn(ComplexF64,5)

h = (u,p,t) -> begin
    u
end

p = 1.
L = MatrixFreeOperator(h, (p,), size=(5,))

N = (u,p,t) -> begin
    u.^2
end

prob = SplitODEProblem(L,N,b,(0,1.), p)
sol1=solve(prob, Tsit5(),dt=0.05) # success
sol2=solve(prob, ETDRK4(krylov=true),dt=0.05) # error:
ERROR: MethodError: no method matching zero(::Type{Any})
Closest candidates are:
  zero(::Type{Union{Missing, T}}) where T at missing.jl:105
  zero(::Union{Type{P}, P}) where P<:Dates.Period at /Applications/Julia-1.8.app/Contents/Resources/julia/share/julia/stdlib/v1.8/Dates/src/periods.jl:53
  zero(::ForwardDiff.Partials) at ~/.julia/packages/ForwardDiff/wAaVJ/src/partials.jl:39
  ...
Stacktrace:
  [1] zero(#unused#::Type{Any})
    @ Base ./missing.jl:106
  [2] zeros(#unused#::Type{Any}, dims::Tuple{Int64, Int64})
    @ Base ./array.jl:589
  [3] arnoldi(A::MatrixFreeOperator{var"#5#6", Tuple{Float64}, Tuple{Int64}, Bool}, b::Vector{ComplexF64}; m::Int64, ishermitian::Bool, kwargs::Base.Pairs{Symbol, Any, Tuple{Symbol, Symbol}, NamedTuple{(:opnorm, :iop), Tuple{typeof(LinearAlgebra.opnorm), Int64}}})
    @ ExponentialUtilities ~/.julia/packages/ExponentialUtilities/gGKOB/src/arnoldi.jl:107
  [4] perform_step!(integrator::OrdinaryDiffEq.ODEIntegrator{ETDRK4{Val{:forward}, true, nothing}, false, Vector{ComplexF64}, Nothing, Float64, Float64, Float64, Float64, Float64, Float64, Vector{Vector{ComplexF64}}, ODESolution{ComplexF64, 2, Vector{Vector{ComplexF64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{ComplexF64}}}, ODEProblem{Vector{ComplexF64}, Tuple{Float64, Float64}, false, Float64, SplitFunction{false, ODEFunction{true, MatrixFreeOperator{var"#5#6", Tuple{Float64}, Tuple{Int64}, Bool}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, ODEFunction{false, var"#7#8", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, LinearAlgebra.UniformScaling{Bool}, Vector{ComplexF64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SplitODEProblem{false}}, ETDRK4{Val{:forward}, true, nothing}, OrdinaryDiffEq.InterpolationData{SplitFunction{false, ODEFunction{true, MatrixFreeOperator{var"#5#6", Tuple{Float64}, Tuple{Int64}, Bool}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, ODEFunction{false, var"#7#8", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, LinearAlgebra.UniformScaling{Bool}, Vector{ComplexF64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{Vector{ComplexF64}}, Vector{Float64}, Vector{Vector{Vector{ComplexF64}}}, OrdinaryDiffEq.ETDRK4ConstantCache{Nothing, Nothing}}, DiffEqBase.DEStats}, SplitFunction{false, ODEFunction{true, MatrixFreeOperator{var"#5#6", Tuple{Float64}, Tuple{Int64}, Bool}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, ODEFunction{false, var"#7#8", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, LinearAlgebra.UniformScaling{Bool}, Vector{ComplexF64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, OrdinaryDiffEq.ETDRK4ConstantCache{Nothing, Nothing}, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, Float64, PIController{Int64}, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(LinearAlgebra.opnorm), Nothing, CallbackSet{Tuple{}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, Nothing, Nothing, Int64, Tuple{}, Tuple{}, Tuple{}}, Vector{ComplexF64}, ComplexF64, Nothing, OrdinaryDiffEq.DefaultInit}, cache::OrdinaryDiffEq.ETDRK4ConstantCache{Nothing, Nothing}, repeat_step::Bool)
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/vGDGq/src/perform_step/exponential_rk_perform_step.jl:356
  [5] perform_step!(integrator::OrdinaryDiffEq.ODEIntegrator{ETDRK4{Val{:forward}, true, nothing}, false, Vector{ComplexF64}, Nothing, Float64, Float64, Float64, Float64, Float64, Float64, Vector{Vector{ComplexF64}}, ODESolution{ComplexF64, 2, Vector{Vector{ComplexF64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{ComplexF64}}}, ODEProblem{Vector{ComplexF64}, Tuple{Float64, Float64}, false, Float64, SplitFunction{false, ODEFunction{true, MatrixFreeOperator{var"#5#6", Tuple{Float64}, Tuple{Int64}, Bool}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, ODEFunction{false, var"#7#8", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, LinearAlgebra.UniformScaling{Bool}, Vector{ComplexF64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SplitODEProblem{false}}, ETDRK4{Val{:forward}, true, nothing}, OrdinaryDiffEq.InterpolationData{SplitFunction{false, ODEFunction{true, MatrixFreeOperator{var"#5#6", Tuple{Float64}, Tuple{Int64}, Bool}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, ODEFunction{false, var"#7#8", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, LinearAlgebra.UniformScaling{Bool}, Vector{ComplexF64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{Vector{ComplexF64}}, Vector{Float64}, Vector{Vector{Vector{ComplexF64}}}, OrdinaryDiffEq.ETDRK4ConstantCache{Nothing, Nothing}}, DiffEqBase.DEStats}, SplitFunction{false, ODEFunction{true, MatrixFreeOperator{var"#5#6", Tuple{Float64}, Tuple{Int64}, Bool}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, ODEFunction{false, var"#7#8", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, LinearAlgebra.UniformScaling{Bool}, Vector{ComplexF64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, OrdinaryDiffEq.ETDRK4ConstantCache{Nothing, Nothing}, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, Float64, PIController{Int64}, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(LinearAlgebra.opnorm), Nothing, CallbackSet{Tuple{}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, Nothing, Nothing, Int64, Tuple{}, Tuple{}, Tuple{}}, Vector{ComplexF64}, ComplexF64, Nothing, OrdinaryDiffEq.DefaultInit}, cache::OrdinaryDiffEq.ETDRK4ConstantCache{Nothing, Nothing})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/vGDGq/src/perform_step/exponential_rk_perform_step.jl:346
  [6] solve!(integrator::OrdinaryDiffEq.ODEIntegrator{ETDRK4{Val{:forward}, true, nothing}, false, Vector{ComplexF64}, Nothing, Float64, Float64, Float64, Float64, Float64, Float64, Vector{Vector{ComplexF64}}, ODESolution{ComplexF64, 2, Vector{Vector{ComplexF64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{ComplexF64}}}, ODEProblem{Vector{ComplexF64}, Tuple{Float64, Float64}, false, Float64, SplitFunction{false, ODEFunction{true, MatrixFreeOperator{var"#5#6", Tuple{Float64}, Tuple{Int64}, Bool}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, ODEFunction{false, var"#7#8", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, LinearAlgebra.UniformScaling{Bool}, Vector{ComplexF64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SplitODEProblem{false}}, ETDRK4{Val{:forward}, true, nothing}, OrdinaryDiffEq.InterpolationData{SplitFunction{false, ODEFunction{true, MatrixFreeOperator{var"#5#6", Tuple{Float64}, Tuple{Int64}, Bool}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, ODEFunction{false, var"#7#8", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, LinearAlgebra.UniformScaling{Bool}, Vector{ComplexF64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{Vector{ComplexF64}}, Vector{Float64}, Vector{Vector{Vector{ComplexF64}}}, OrdinaryDiffEq.ETDRK4ConstantCache{Nothing, Nothing}}, DiffEqBase.DEStats}, SplitFunction{false, ODEFunction{true, MatrixFreeOperator{var"#5#6", Tuple{Float64}, Tuple{Int64}, Bool}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, ODEFunction{false, var"#7#8", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, LinearAlgebra.UniformScaling{Bool}, Vector{ComplexF64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, OrdinaryDiffEq.ETDRK4ConstantCache{Nothing, Nothing}, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, Float64, PIController{Int64}, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(LinearAlgebra.opnorm), Nothing, CallbackSet{Tuple{}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, Nothing, Nothing, Int64, Tuple{}, Tuple{}, Tuple{}}, Vector{ComplexF64}, ComplexF64, Nothing, OrdinaryDiffEq.DefaultInit})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/vGDGq/src/solve.jl:477
  [7] #__solve#562
    @ ~/.julia/packages/OrdinaryDiffEq/vGDGq/src/solve.jl:5 [inlined]
  [8] #solve_call#28
    @ ~/.julia/packages/DiffEqBase/nxrbH/src/solve.jl:287 [inlined]
  [9] #solve_up#30
    @ ~/.julia/packages/DiffEqBase/nxrbH/src/solve.jl:315 [inlined]
 [10] #solve#29
    @ ~/.julia/packages/DiffEqBase/nxrbH/src/solve.jl:299 [inlined]
 [11] top-level scope
    @ ~/Dropbox/Research/Projects/FourierRotatingGPE/rot_gpe_imex.jl:187

Context:

[0c46a032] DifferentialEquations v7.1.0
[9fdde737] DiffEqOperators v4.43.0

julia> versioninfo()
Julia Version 1.8.0-beta3
Commit 3e092a2521 (2022-03-29 15:42 UTC)
Platform Info:
  OS: macOS (arm64-apple-darwin21.3.0)
  CPU: 10 × Apple M1 Max
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-13.0.1 (ORCJIT, apple-m1)
  Threads: 8 on 8 virtual cores
Environment:
  JULIA_PKG_DEVDIR = /Users/abradley/Dropbox/Julia/Dev
  JULIA_NUM_THREADS = 8
  JULIA_EDITOR = code
AshtonSBradley commented 2 years ago

I would be keen to work on a PR for this if someone is willing to some guidance to get started

ChrisRackauckas commented 2 years ago

Try the new FunctionOperator in SciMLOperators.jl. We haven't quite integrated it downstream though, so it might need some hooks in LinearSolve.jl

ChrisRackauckas commented 2 years ago

This operator was deprecated for SciMLOperators.jl's new version. Try there.