SciMLOperators.jl linear operator conversion error in SplitODEProblem

@vpuri3 this is the example code using https://docs.sciml.ai/SciMLOperators/dev/tutorials/fftw/ as a starting point for an fft pde (updated)

using SciMLOperators, LinearAlgebra, FFTW, OrdinaryDiffEq

n = 256
l = 2π

dx = l / n
x  = range(-l/2, l/2, n+1)[1:n] |> Array
u  = @. sin(5x)cos(7x)  + im*cos(x) |> complex  
du = @. 5cos(5x)cos(7x) - 7sin(5x)sin(7x)  - im*sin(x) |> complex
ddu = @. -25sin(5x)cos(7x) - 35cos(5x)sin(7x) - 35cos(5x)sin(7x) - 49sin(5x)cos(7x) - im*cos(x) |> complex

k  = fftfreq(n, 2π*n/l) |> Array
m  = length(k)
tr = plan_fft(u)

Lf = FunctionOperator(
     (du,u,p,t) -> mul!(du, tr, u), u,u,
          isinplace=true,
          T=ComplexF64,
          op_adjoint = (du,u,p,t) -> ldiv!(du, tr, u),
          op_inverse = (du,u,p,t) -> ldiv!(du, tr, u),
          op_adjoint_inverse = (du,u,p,t) -> ldiv!(du, tr, u),
)

ik = im * DiagonalOperator(k)
∂x = Lf \ ik * Lf
∂xx= Lf \ -DiagonalOperator(k.^2) * Lf
Dx = cache_operator(∂x, u)
Dxx = cache_operator(∂xx, u)

@show ≈(mul!(copy(u), Dx, u), du; atol=1e-8)
@show ≈(mul!(copy(u), Dxx, u), ddu; atol=1e-8)

## test 
u0 = randn(ComplexF64,n)
F=plan_fft(u0)

## breaks if attempting to use FunctionOperator for linear part
# Error in converting operator
L = -im*0.5*Dxx
N = (u,_,_) -> -im*abs2.(u).*u

prob = SplitODEProblem(L,N,u0,(0,1.))
sol = solve(prob, ETDRK4(), dt=0.05, saveat=0.5);

Throws the error

ERROR: MethodError: Cannot `convert` an object of type 
  var"#3#7" to an object of type 
  AbstractMatrix

Closest candidates are:
  convert(::Type{T}, ::T) where T
   @ Base Base.jl:84
  convert(::Type{T}, ::Factorization) where T<:AbstractArray
   @ LinearAlgebra ~/.julia/juliaup/julia-1.10.0-rc2+0.aarch64.apple.darwin14/share/julia/stdlib/v1.10/LinearAlgebra/src/factorization.jl:108
  convert(::Type{T}, ::LinearAlgebra.AbstractQ) where T<:AbstractArray
   @ LinearAlgebra ~/.julia/juliaup/julia-1.10.0-rc2+0.aarch64.apple.darwin14/share/julia/stdlib/v1.10/LinearAlgebra/src/abstractq.jl:45
  ...

Stacktrace:
  [1] convert(::Type{…}, L::FunctionOperator{…})
    @ SciMLOperators ~/.julia/packages/SciMLOperators/hFLNz/src/func.jl:510
  [2] (::SciMLOperators.var"#118#119")(op::FunctionOperator{…})
    @ SciMLOperators ~/.julia/packages/SciMLOperators/hFLNz/src/basic.jl:538
  [3] MappingRF
    @ ./reduce.jl:100 [inlined]
  [4] afoldl(::Base.MappingRF{…}, ::Base._InitialValue, ::FunctionOperator{…}, ::SciMLOperators.ScaledOperator{…}, ::FunctionOperator{…})
    @ Base ./operators.jl:544
  [5] _foldl_impl(op::Base.MappingRF{…}, init::Base._InitialValue, itr::Tuple{…})
    @ Base ./reduce.jl:68
  [6] foldl_impl(op::Base.MappingRF{…}, nt::Base._InitialValue, itr::Tuple{…})
    @ Base ./reduce.jl:48
  [7] mapfoldl_impl(f::SciMLOperators.var"#118#119", op::typeof(Base.mul_prod), nt::Base._InitialValue, itr::Tuple{…})
    @ Base ./reduce.jl:44
  [8] mapfoldl(f::Function, op::Function, itr::Tuple{…}; init::Base._InitialValue)
    @ Base ./reduce.jl:175
  [9] mapfoldl
    @ Base ./reduce.jl:175 [inlined]
 [10] mapreduce
    @ Base ./reduce.jl:307 [inlined]
 [11] prod(f::Function, a::Tuple{FunctionOperator{…}, SciMLOperators.ScaledOperator{…}, FunctionOperator{…}})
    @ Base ./reduce.jl:591
 [12] convert(::Type{AbstractMatrix}, L::SciMLOperators.ComposedOperator{ComplexF64, Tuple{…}, Tuple{…}})
    @ SciMLOperators ~/.julia/packages/SciMLOperators/hFLNz/src/basic.jl:538
 [13] convert(::Type{…}, L::SciMLOperators.ScaledOperator{…})
    @ SciMLOperators ~/.julia/packages/SciMLOperators/hFLNz/src/basic.jl:223
 [14] alg_cache
    @ ~/.julia/packages/OrdinaryDiffEq/Y6bIX/src/caches/linear_nonlinear_caches.jl:85 [inlined]
 [15] __init(prob::ODEProblem{…}, alg::ETDRK4{…}, timeseries_init::Tuple{}, ts_init::Tuple{}, ks_init::Tuple{}, recompile::Type{…}; saveat::Float64, tstops::Tuple{}, d_discontinuities::Tuple{}, save_idxs::Nothing, save_everystep::Bool, save_on::Bool, save_start::Bool, save_end::Nothing, callback::Nothing, dense::Bool, calck::Bool, dt::Float64, dtmin::Nothing, dtmax::Float64, force_dtmin::Bool, adaptive::Bool, gamma::Int64, abstol::Nothing, reltol::Nothing, qmin::Int64, qmax::Int64, qsteady_min::Int64, qsteady_max::Int64, beta1::Nothing, beta2::Nothing, qoldinit::Int64, controller::Nothing, fullnormalize::Bool, failfactor::Int64, maxiters::Int64, internalnorm::typeof(DiffEqBase.ODE_DEFAULT_NORM), internalopnorm::typeof(opnorm), isoutofdomain::typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), unstable_check::typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), verbose::Bool, timeseries_errors::Bool, dense_errors::Bool, advance_to_tstop::Bool, stop_at_next_tstop::Bool, initialize_save::Bool, progress::Bool, progress_steps::Int64, progress_name::String, progress_message::typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), progress_id::Symbol, userdata::Nothing, allow_extrapolation::Bool, initialize_integrator::Bool, alias_u0::Bool, alias_du0::Bool, initializealg::OrdinaryDiffEq.DefaultInit, kwargs::@Kwargs{})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/Y6bIX/src/solve.jl:334
 [16] __init (repeats 5 times)
    @ ~/.julia/packages/OrdinaryDiffEq/Y6bIX/src/solve.jl:10 [inlined]
 [17] #__solve#746
    @ ~/.julia/packages/OrdinaryDiffEq/Y6bIX/src/solve.jl:5 [inlined]
 [18] __solve
    @ ~/.julia/packages/OrdinaryDiffEq/Y6bIX/src/solve.jl:1 [inlined]
 [19] solve_call(_prob::ODEProblem{…}, args::ETDRK4{…}; merge_callbacks::Bool, kwargshandle::Nothing, kwargs::@Kwargs{…})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/LUZAV/src/solve.jl:561
 [20] solve_up(prob::ODEProblem{…}, sensealg::Nothing, u0::Vector{…}, p::SciMLBase.NullParameters, args::ETDRK4{…}; kwargs::@Kwargs{…})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/LUZAV/src/solve.jl:1010
 [21] solve_up
    @ ~/.julia/packages/DiffEqBase/LUZAV/src/solve.jl:996 [inlined]
 [22] solve(prob::ODEProblem{…}, args::ETDRK4{…}; sensealg::Nothing, u0::Nothing, p::Nothing, wrap::Val{…}, kwargs::@Kwargs{…})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/LUZAV/src/solve.jl:933
 [23] top-level scope
    @ REPL[27]:1
Some type information was truncated. Use `show(err)` to see complete types.

Using

  [7a1cc6ca] FFTW v1.7.2
  [1dea7af3] OrdinaryDiffEq v6.60.0
  [c0aeaf25] SciMLOperators v0.3.7
  [37e2e46d] LinearAlgebra

julia> versioninfo()
Julia Version 1.10.0-rc2
Commit dbb9c46795b (2023-12-03 15:25 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: macOS (arm64-apple-darwin22.4.0)
  CPU: 10 × Apple M1 Max
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-15.0.7 (ORCJIT, apple-m1)
  Threads: 11 on 8 virtual cores
Environment:
  JULIA_PKG_DEVDIR = /Users/abradley/Dropbox/Julia/Dev
  JULIA_NUM_THREADS = 8
  JULIA_PKG_SERVER = us-west.pkg.julialang.org
  JULIA_EDITOR = code

May 27 '23 18:05 AshtonSBradley

It looks like your time-stepper, ETDRK2, and most other exponential methods require full matrices. They depend on ExponentialUtilities.jl that, upon a cursory look, requires full matrices.

I tried your problem with Tsit5 which doesn't require a full matrix, and it worked fine. In general, the fastest way to solve such a problem is with an IMEX timestepper and a krylov linear solver.

F = SplitFunction(L, N; jac_prototype = L)
prob = ODEProblem(F, u0, (0,1.))
sol = solve(prob, TRBDF2(autodiff = false, linsolve = KrylovJL_GMRES()))

However, that functionality is broken right now. Split operators won't work with implicit/semiimplicit solvers, but we're looking to fix that in the next 1-2 weeks. Till then I'd recommend sticking with any of the stiff explicit solvers.

May 27 '23 19:05 vpuri3

Thanks for your reply, I'll follow your advice. Being able to write to this interface is already great progress. Thanks a lot for your work!

May 27 '23 22:05 AshtonSBradley

of course. I'll let you know here when the split function + imex is done. also, if you are using fourier stuff extensively, I've implemented almost all relevant operators with scimloperators here:

https://github.com/CalculustJL/FourierSpaces.jl https://github.com/CalculustJL/FourierSpaces.jl/blob/master/src/trans_operators.jl https://github.com/CalculustJL/FourierSpaces.jl/blob/master/src/phys_operators.jl

here's an example of usage

https://github.com/CalculustJL/FourierSpaces.jl/blob/master/examples/fourier1d/heat.jl

May 28 '23 01:05 vpuri3

Thanks a lot for the links, I am reading the code and finding it quite enlightening. I think I could contribute some additional examples such as vortex lattice formation in a Bose-Einstein condensate. Any updates regarding split + imex?

Jun 28 '23 17:06 AshtonSBradley

In FourierSpaces.jl, or SciMLOperators.jl will it be possibly to evaluate mixed spectral operators of the form, e.g. x*d/dy, or f(x,d/dy)+g(d/dx,y)

Jun 30 '23 22:06 AshtonSBradley

@AshtonSBradley sorry for the delayed response. Contributions are welcome. Feel free to make an issue in FourierSpaces.jl and we can take it from there.

Jul 02 '23 14:07 vpuri3

No updates on the Split ODEs prob + implicit methods. That is going to be a substantial change in OrdinaryDiffEq and would require a week or so of dedicated effort. Note this is one of the two remaining pieces of https://github.com/SciML/SciMLOperators.jl/issues/142 - a six-month long effort.

I won't hit that use case in my school work till at least the fall. Best case scenario I get to it during vacation in August, unless @gaurav-arya has time to hit that before me.

Jul 02 '23 14:07 vpuri3

In FourierSpaces.jl, or SciMLOperators.jl will it be possibly to evaluate mixed spectral operators of the form, e.g. x*d/dy, or f(x,d/dy)+g(d/dx,y)

In FourierSpaces,

V = FourierSpace(Nx, Ny)
x, y = points(V)
Dx, Dy = gradientOp(V)

L = DiagonalOperator(x) * Dy

@AshtonSBradley is this what you were looking for?

Jul 02 '23 14:07 vpuri3

No updates on the Split ODEs prob + implicit methods. That is going to be a substantial change in OrdinaryDiffEq and would require a week or so of dedicated effort. Note this is one of the two remaining pieces of #142 - a six-month long effort.

I won't hit that use case in my school work till at least the fall. Best case scenario I get to it during vacation in August, unless @gaurav-arya has time to hit that before me.

I am curious if there is any progress along this line of Split ODE prob + implicit methods?

Dec 07 '23 00:12 AshtonSBradley

I thought I handled that a few weeks ago and closed the issue in OrdinaryDiffEq?

Dec 07 '23 01:12 ChrisRackauckas

the error has changed

sol = solve(prob, ETDRK4(), dt=0.05, saveat=0.5);
ERROR: MethodError: Cannot `convert` an object of type 
  var"#11#15" to an object of type 
  AbstractMatrix

Closest candidates are:
  convert(::Type{T}, ::T) where T
   @ Base Base.jl:84
  convert(::Type{T}, ::Factorization) where T<:AbstractArray
   @ LinearAlgebra ~/.julia/juliaup/julia-1.10.0-rc1+0.aarch64.apple.darwin14/share/julia/stdlib/v1.10/LinearAlgebra/src/factorization.jl:108
  convert(::Type{T}, ::T) where T<:AbstractArray
   @ Base abstractarray.jl:16
  ...

Stacktrace:
  [1] convert(::Type{…}, L::FunctionOperator{…})
    @ SciMLOperators ~/.julia/packages/SciMLOperators/hFLNz/src/func.jl:510
  [2] (::SciMLOperators.var"#118#119")(op::FunctionOperator{…})
    @ SciMLOperators ~/.julia/packages/SciMLOperators/hFLNz/src/basic.jl:538
  [3] MappingRF
    @ ./reduce.jl:100 [inlined]
  [4] afoldl(::Base.MappingRF{…}, ::Base._InitialValue, ::FunctionOperator{…}, ::SciMLOperators.ScaledOperator{…}, ::FunctionOperator{…})
    @ Base ./operators.jl:544
  [5] _foldl_impl(op::Base.MappingRF{…}, init::Base._InitialValue, itr::Tuple{…})
    @ Base ./reduce.jl:68
  [6] foldl_impl(op::Base.MappingRF{…}, nt::Base._InitialValue, itr::Tuple{…})
    @ Base ./reduce.jl:48
  [7] mapfoldl_impl(f::SciMLOperators.var"#118#119", op::typeof(Base.mul_prod), nt::Base._InitialValue, itr::Tuple{…})
    @ Base ./reduce.jl:44
  [8] mapfoldl(f::Function, op::Function, itr::Tuple{…}; init::Base._InitialValue)
    @ Base ./reduce.jl:175
  [9] mapfoldl
    @ Base ./reduce.jl:175 [inlined]
 [10] mapreduce
    @ Base ./reduce.jl:307 [inlined]
 [11] prod(f::Function, a::Tuple{FunctionOperator{…}, SciMLOperators.ScaledOperator{…}, FunctionOperator{…}})
    @ Base ./reduce.jl:591
 [12] convert(::Type{AbstractMatrix}, L::SciMLOperators.ComposedOperator{ComplexF64, Tuple{…}, Tuple{…}})
    @ SciMLOperators ~/.julia/packages/SciMLOperators/hFLNz/src/basic.jl:538
 [13] convert(::Type{…}, L::SciMLOperators.ScaledOperator{…})
    @ SciMLOperators ~/.julia/packages/SciMLOperators/hFLNz/src/basic.jl:223
 [14] alg_cache
    @ ~/.julia/packages/OrdinaryDiffEq/Y6bIX/src/caches/linear_nonlinear_caches.jl:85 [inlined]
 [15] __init(prob::ODEProblem{…}, alg::ETDRK4{…}, timeseries_init::Tuple{}, ts_init::Tuple{}, ks_init::Tuple{}, recompile::Type{…}; saveat::Float64, tstops::Tuple{}, d_discontinuities::Tuple{}, save_idxs::Nothing, save_everystep::Bool, save_on::Bool, save_start::Bool, save_end::Nothing, callback::Nothing, dense::Bool, calck::Bool, dt::Float64, dtmin::Nothing, dtmax::Float64, force_dtmin::Bool, adaptive::Bool, gamma::Int64, abstol::Nothing, reltol::Nothing, qmin::Int64, qmax::Int64, qsteady_min::Int64, qsteady_max::Int64, beta1::Nothing, beta2::Nothing, qoldinit::Int64, controller::Nothing, fullnormalize::Bool, failfactor::Int64, maxiters::Int64, internalnorm::typeof(DiffEqBase.ODE_DEFAULT_NORM), internalopnorm::typeof(opnorm), isoutofdomain::typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), unstable_check::typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), verbose::Bool, timeseries_errors::Bool, dense_errors::Bool, advance_to_tstop::Bool, stop_at_next_tstop::Bool, initialize_save::Bool, progress::Bool, progress_steps::Int64, progress_name::String, progress_message::typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), progress_id::Symbol, userdata::Nothing, allow_extrapolation::Bool, initialize_integrator::Bool, alias_u0::Bool, alias_du0::Bool, initializealg::OrdinaryDiffEq.DefaultInit, kwargs::@Kwargs{})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/Y6bIX/src/solve.jl:334
 [16] __init (repeats 5 times)
    @ ~/.julia/packages/OrdinaryDiffEq/Y6bIX/src/solve.jl:10 [inlined]
 [17] #__solve#746
    @ ~/.julia/packages/OrdinaryDiffEq/Y6bIX/src/solve.jl:5 [inlined]
 [18] __solve
    @ ~/.julia/packages/OrdinaryDiffEq/Y6bIX/src/solve.jl:1 [inlined]
 [19] solve_call(_prob::ODEProblem{…}, args::ETDRK4{…}; merge_callbacks::Bool, kwargshandle::Nothing, kwargs::@Kwargs{…})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/8lmsn/src/solve.jl:561
 [20] solve_up(prob::ODEProblem{…}, sensealg::Nothing, u0::Vector{…}, p::SciMLBase.NullParameters, args::ETDRK4{…}; kwargs::@Kwargs{…})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/8lmsn/src/solve.jl:1010
 [21] solve_up
    @ ~/.julia/packages/DiffEqBase/8lmsn/src/solve.jl:996 [inlined]
 [22] solve(prob::ODEProblem{…}, args::ETDRK4{…}; sensealg::Nothing, u0::Nothing, p::Nothing, wrap::Val{…}, kwargs::@Kwargs{…})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/8lmsn/src/solve.jl:933
 [23] top-level scope
    @ REPL[57]:1
Some type information was truncated. Use `show(err)` to see complete types.

Dec 07 '23 08:12 AshtonSBradley

I thought I handled that a few weeks ago and closed the issue in OrdinaryDiffEq?

@ChrisRackauckas do you mean https://github.com/SciML/OrdinaryDiffEq.jl/issues/2009 ? It doesn't seem to have entirely fixed things, as above

Dec 07 '23 20:12 AshtonSBradley

Can you open it up on OrdinaryDiffEq? I don't know if I'll get to this soon but it would be good to track.

Dec 10 '23 10:12 ChrisRackauckas

Can you open it up on OrdinaryDiffEq? I don't know if I'll get to this soon but it would be good to track.

#2078

Dec 10 '23 20:12 AshtonSBradley

SciMLOperators.jl SciMLOperators.jl copied to clipboard

linear operator conversion error in SplitODEProblem

SciMLOperators.jl
SciMLOperators.jl copied to clipboard