linear operator conversion error in SplitODEProblem #191

AshtonSBradley · 2023-05-27T18:40:34Z

@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

vpuri3 · 2023-05-27T19:09:17Z

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.

AshtonSBradley · 2023-05-27T22:47:08Z

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!

vpuri3 · 2023-05-28T01:47:04Z

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

AshtonSBradley · 2023-06-28T17:56:42Z

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?

AshtonSBradley · 2023-06-30T22:26:30Z

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)

vpuri3 · 2023-07-02T14:36:03Z

@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.

vpuri3 · 2023-07-02T14:45:51Z

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.

vpuri3 · 2023-07-02T14:48:11Z

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?

AshtonSBradley · 2023-12-07T00:58:00Z

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?

ChrisRackauckas · 2023-12-07T01:06:57Z

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

AshtonSBradley · 2023-12-07T08:59:13Z

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.

AshtonSBradley · 2023-12-07T20:15:08Z

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

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

ChrisRackauckas · 2023-12-10T10:09:54Z

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.

AshtonSBradley · 2023-12-10T20:38:07Z

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

vpuri3 mentioned this issue Jul 24, 2023

Support user-specified W_prototype <: AbstractSciMLOperator SciML/OrdinaryDiffEq.jl#1996

Merged

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

linear operator conversion error in SplitODEProblem #191

linear operator conversion error in SplitODEProblem #191

AshtonSBradley commented May 27, 2023 •

edited

Loading

vpuri3 commented May 27, 2023

AshtonSBradley commented May 27, 2023 •

edited

Loading

vpuri3 commented May 28, 2023

AshtonSBradley commented Jun 28, 2023

AshtonSBradley commented Jun 30, 2023

vpuri3 commented Jul 2, 2023

vpuri3 commented Jul 2, 2023 •

edited

Loading

vpuri3 commented Jul 2, 2023 •

edited

Loading

AshtonSBradley commented Dec 7, 2023

ChrisRackauckas commented Dec 7, 2023

AshtonSBradley commented Dec 7, 2023

AshtonSBradley commented Dec 7, 2023

ChrisRackauckas commented Dec 10, 2023

AshtonSBradley commented Dec 10, 2023 •

edited

Loading

linear operator conversion error in SplitODEProblem #191

linear operator conversion error in SplitODEProblem #191

Comments

AshtonSBradley commented May 27, 2023 • edited Loading

vpuri3 commented May 27, 2023

AshtonSBradley commented May 27, 2023 • edited Loading

vpuri3 commented May 28, 2023

AshtonSBradley commented Jun 28, 2023

AshtonSBradley commented Jun 30, 2023

vpuri3 commented Jul 2, 2023

vpuri3 commented Jul 2, 2023 • edited Loading

vpuri3 commented Jul 2, 2023 • edited Loading

AshtonSBradley commented Dec 7, 2023

ChrisRackauckas commented Dec 7, 2023

AshtonSBradley commented Dec 7, 2023

AshtonSBradley commented Dec 7, 2023

ChrisRackauckas commented Dec 10, 2023

AshtonSBradley commented Dec 10, 2023 • edited Loading

AshtonSBradley commented May 27, 2023 •

edited

Loading

AshtonSBradley commented May 27, 2023 •

edited

Loading

vpuri3 commented Jul 2, 2023 •

edited

Loading

vpuri3 commented Jul 2, 2023 •

edited

Loading

AshtonSBradley commented Dec 10, 2023 •

edited

Loading