Error when running the example code for GPU-acceleration using diffeqpy

[gonzalovivares]

Hi all, first of all thanks for the development of diffeqpy, since I use it my scientific modelling work has boosted (I'm simply a user of diffeqpy with limited knowledge about software engineering).

I have been planning to upgrade my code in order to use the GPU of my laptop (NVIDIA RTX-5500), by using CUDA to reduce the simulation time of my models. However, when running the example code given in the repository for the use of diffeqpy.cuda, I encountered two run errors. First, the de.jit32() function gave me the following error:

The code I run was:

from diffeqpy import de
import random
from diffeqpy import cuda

def f(u,p,t):
    x, y, z = u
    sigma, rho, beta = p
    return [sigma * (y - x), x * (rho - z) - y, x * y - beta * z]

u0 = [1.0,0.0,0.0]
tspan = (0., 100.)
p = [10.0,28.0,8/3]
prob = de.ODEProblem(f, u0, tspan, p)
fast_prob = de.jit32(prob)
sol = de.solve(fast_prob,saveat=0.01)

def prob_func(prob,i,rep):
  return de.remake(prob,u0=[random.uniform(0, 1)*u0[i] for i in range(0,3)],
            p=[random.uniform(0, 1)*p[i] for i in range(0,3)])

ensembleprob = de.EnsembleProblem(fast_prob, prob_func=prob_func, safetycopy=False)

sol = de.solve(ensembleprob,cuda.GPUTsit5(),cuda.EnsembleGPUKernel(cuda.CUDABackend()),trajectories=10000,saveat=0.01)

And the first error I encountered:

Traceback (most recent call last):
  File "C:\Users\____\Python\block_1_dynamic_adaptation\reporting_results\delete_later.py", line 14, in <module>
    fast_prob = de.jit32(prob)
                ^^^^^^^^^^^^^^
  File "C:\Users\___\.julia\packages\PythonCall\wXfah\src\jlwrap\any.jl", line 208, in __call__
    return self._jl_callmethod($(pyjl_methodnum(pyjlany_call)), args, kwargs)
       ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
juliacall.JuliaError: A completed system is required. Call `complete` or `structural_simplify` on the system before creating an `ODEProblem`
Stacktrace:
  [1] error(s::String)
    @ Base .\error.jl:35
  [2] (SciMLBase.ODEProblem{true, SciMLBase.AutoSpecialize})(sys::ModelingToolkit.ODESystem, u0map::Vector{Float32}, tspan::Tuple{Float32, Float32}, parammap::Vector{Float32}; callback::Nothing, check_length::Bool, warn_initialize_determined::Bool, kwargs::@Kwargs{})
    @ ModelingToolkit C:\Users\___\.julia\packages\ModelingToolkit\Mxj1Q\src\systems\diffeqs\abstractodesystem.jl:1031
  [3] (SciMLBase.ODEProblem{true, SciMLBase.AutoSpecialize})(sys::ModelingToolkit.ODESystem, u0map::Vector{Float32}, tspan::Tuple{Float32, Float32}, parammap::Vector{Float32})
    @ ModelingToolkit C:\Users\___\.julia\packages\ModelingToolkit\Mxj1Q\src\systems\diffeqs\abstractodesystem.jl:1023
  [4] (SciMLBase.ODEProblem{true})(::ModelingToolkit.ODESystem, ::Vector{Float32}, ::Vararg{Any}; kwargs::@Kwargs{})
    @ ModelingToolkit C:\Users\___\.julia\packages\ModelingToolkit\Mxj1Q\src\systems\diffeqs\abstractodesystem.jl:1010
  [5] (SciMLBase.ODEProblem{true})(::ModelingToolkit.ODESystem, ::Vector{Float32}, ::Vararg{Any})
    @ ModelingToolkit C:\Users\___\.julia\packages\ModelingToolkit\Mxj1Q\src\systems\diffeqs\abstractodesystem.jl:1009
  [6] #ODEProblem#728
    @ C:\Users\___\.julia\packages\ModelingToolkit\Mxj1Q\src\systems\diffeqs\abstractodesystem.jl:999 [inlined]
  [7] ODEProblem
    @ C:\Users\___\.julia\packages\ModelingToolkit\Mxj1Q\src\systems\diffeqs\abstractodesystem.jl:998 [inlined]
  [8] jit(x::SciMLBase.ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, false, PyList{Any}, SciMLBase.ODEFunction{false, SciMLBase.AutoSpecialize, ComposedFunction{typeof(SciMLBasePythonCallExt._pyconvert), Py}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing, Dict{Symbol, Union{Expr, Number, Symbol}}}, Nothing, Nothing}, @Kwargs{}, SciMLBase.StandardODEProblem})
    @ Main .\none:1
  [9] pyjlany_call(self::typeof(jit), args_::Py, kwargs_::Py)
    @ PythonCall C:\Users\___\.julia\packages\PythonCall\wXfah\src\jlwrap\any.jl:37
 [10] _pyjl_callmethod(f::Any, self_::Ptr{PythonCall.C.PyObject}, args_::Ptr{PythonCall.C.PyObject}, nargs::Int64)
    @ PythonCall C:\Users\___\.julia\packages\PythonCall\wXfah\src\jlwrap\base.jl:69
 [11] _pyjl_callmethod(o::Ptr{PythonCall.C.PyObject}, args::Ptr{PythonCall.C.PyObject})
    @ PythonCall.C C:\Users\____\.julia\packages\PythonCall\wXfah\src\cpython\jlwrap.jl:47

I read a previous threat about a similar error, and the conclusion was that de.jit() might not be useful in some cases, but not sure if it would also apply to this example, or if there is something I might need to do on my computer (I'm missing any installation in Julia, I need a conda interpreter rather than a virtualenv in Pycharm...?).

Anyway, I skipped this error by using the normal de.ODEProblem to see if the solver using GPU could work, and then I got this other error, related to the cuda sub-package:

Traceback (most recent call last):
  File "C:\Users\___\Python\block_1_dynamic_adaptation\reporting_results\delete_later.py", line 23, in <module>
    sol = de.solve(ensembleprob,cuda.GPUTsit5(),cuda.EnsembleGPUKernel(cuda.CUDABackend()),trajectories=10000,saveat=0.01)
          ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "C:\Users\___\.julia\packages\PythonCall\wXfah\src\jlwrap\any.jl", line 208, in __call__
    return self._jl_callmethod($(pyjl_methodnum(pyjlany_call)), args, kwargs)
       ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
juliacall.JuliaError: InvalidIRError: compiling MethodInstance for DiffEqGPU.gpu_ode_asolve_kernel(::KernelAbstractions.CompilerMetadata{KernelAbstractions.NDIteration.DynamicSize, KernelAbstractions.NDIteration.DynamicCheck, Nothing, CartesianIndices{1, Tuple{Base.OneTo{Int64}}}, KernelAbstractions.NDIteration.NDRange{1, KernelAbstractions.NDIteration.DynamicSize, KernelAbstractions.NDIteration.DynamicSize, CartesianIndices{1, Tuple{Base.OneTo{Int64}}}, CartesianIndices{1, Tuple{Base.OneTo{Int64}}}}}, ::CUDA.CuDeviceVector{DiffEqGPU.ImmutableODEProblem{StaticArraysCore.SVector{3, Float64}, Tuple{Float64, Float64}, false, PyList{Any}, SciMLBase.ODEFunction{false, SciMLBase.AutoSpecialize, ComposedFunction{typeof(SciMLBasePythonCallExt._pyconvert), Py}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing, Dict{Symbol, Union{Expr, Number, Symbol}}}, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, SciMLBase.StandardODEProblem}, 1}, ::DiffEqGPU.GPUTsit5, ::CUDA.CuDeviceMatrix{StaticArraysCore.SVector{3, Float64}, 1}, ::CUDA.CuDeviceMatrix{Float64, 1}, ::Float64, ::SciMLBase.CallbackSet{Tuple{}, Tuple{}}, ::Nothing, ::Float64, ::Float64, ::StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, ::Val{false}) resulted in invalid LLVM IR
Reason: unsupported call to an unknown function (call to julia.new_gc_frame)
Reason: unsupported call to an unknown function (call to julia.push_gc_frame)
Reason: unsupported call to an unknown function (call to julia.get_gc_frame_slot)
Reason: unsupported dynamic function invocation (call to var"#pycall#59"(kwargs::Base.Pairs{Symbol, V, Tuple{Vararg{Symbol, N}}, NamedTuple{names, T}} where {V, N, names, T<:Tuple{Vararg{Any, N}}}, ::typeof(pycall), f, args...) @ PythonCall C:\Users\___\.julia\packages\PythonCall\wXfah\src\abstract\object.jl:218)

And a long list of Stacktraces that I could share if you need it.

I am wondering if there is a solution to this problem. I have been able to succesfully run the examples of the diffeqGPU.jl in Julia, so the GPU seems to work so far, but if there is no solution for me using Python, I might need to rewrite my models from Python to Julia (and learn Julia btw) to get advantage of my GPU, what do you think?

Thanks a lot.

[jaswin90]

I am also facing the same problem!

[ChrisRackauckas]

Fixed in the new release. Let me know if there are any remaining issues.

[gonzalovivares]

Hi Chris, thanks for your fix! Now the function de.fit32( ) works perfectly, but when using the solver via GPU I still get the following error:

  File "C:\Users\_____.py", line 31, in <module>
    sol = de.solve(ensembleprob,cuda.GPUTsit5(),cuda.EnsembleGPUKernel(cuda.CUDABackend()),trajectories=10000,saveat=0.01)
          ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "C:\Users\vivar002\.julia\packages\PythonCall\vW5OI\src\JlWrap\any.jl", line 223, in __call__
    return self._jl_callmethod($(pyjl_methodnum(pyjlany_call)), args, kwargs)
       ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
juliacall.JuliaError: InvalidIRError: compiling MethodInstance for DiffEqGPU.gpu_ode_asolve_kernel(::KernelAbstractions.CompilerMetadata{KernelAbstractions.NDIteration.DynamicSize, KernelAbstractions.NDIteration.DynamicCheck, Nothing, CartesianIndices{1, Tuple{Base.OneTo{Int64}}}, KernelAbstractions.NDIteration.NDRange{1, KernelAbstractions.NDIteration.DynamicSize, KernelAbstractions.NDIteration.DynamicSize, CartesianIndices{1, Tuple{Base.OneTo{Int64}}}, CartesianIndices{1, Tuple{Base.OneTo{Int64}}}}}, ::CUDA.CuDeviceVector{DiffEqGPU.ImmutableODEProblem{StaticArraysCore.SVector{3, Float64}, Tuple{Float32, Float32}, true, PyList{Any}, SciMLBase.ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#f#684"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x87082b8a, 0x9b0188e9, 0x33d32870, 0x3c85d5f8, 0x086ff446), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xa674921c, 0xfd1aa460, 0xb5b2c07d, 0xcacabc44, 0x7cc4471f), Nothing}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, ModelingToolkit.var"#834#generated_observed#693"{Bool, ModelingToolkit.ODESystem, Dict{Any, Any}, Vector{Any}}, Nothing, ModelingToolkit.ODESystem, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, SciMLBase.StandardODEProblem}, 1}, ::DiffEqGPU.GPUTsit5, ::CUDA.CuDeviceMatrix{StaticArraysCore.SVector{3, Float32}, 1}, ::CUDA.CuDeviceMatrix{Float32, 1}, ::Float32, ::SciMLBase.CallbackSet{Tuple{}, Tuple{}}, ::Nothing, ::Float32, ::Float32, ::StepRangeLen{Float32, Float64, Float64, Int64}, ::Val{false}) resulted in invalid LLVM IR
Reason: unsupported call to an unknown function (call to julia.new_gc_frame)
Stacktrace:
  [1] macro expansion
    @ C:\Users\vivar002\.julia\packages\PythonCall\vW5OI\src\Convert\pyconvert.jl:355
  [2] macro expansion
    @ C:\Users\vivar002\.julia\packages\PythonCall\vW5OI\src\Core\Py.jl:132
  [3] pyconvert
    @ C:\Users\vivar002\.julia\packages\PythonCall\vW5OI\src\Convert\pyconvert.jl:372
  [4] getindex
    @ C:\Users\vivar002\.julia\packages\PythonCall\vW5OI\src\Wrap\PyList.jl:29
  [5] macro expansion
    @ C:\Users\vivar002\.julia\packages\SymbolicUtils\c0xQb\src\code.jl:468
  [6] macro expansion
    @ C:\Users\vivar002\.julia\packages\RuntimeGeneratedFunctions\Yo8zx\src\RuntimeGeneratedFunctions.jl:163
  [7] macro expansion
    @ .\none:0
  [8] generated_callfunc
    @ .\none:0
  [9] RuntimeGeneratedFunction
    @ C:\Users\vivar002\.julia\packages\RuntimeGeneratedFunctions\Yo8zx\src\RuntimeGeneratedFunctions.jl:150
 [10] f
    @ C:\Users\vivar002\.julia\packages\ModelingToolkit\Mxj1Q\src\systems\diffeqs\abstractodesystem.jl:336
 [11] ODEFunction
    @ C:\Users\vivar002\.julia\packages\SciMLBase\Dwomw\src\scimlfunctions.jl:2168
 [12] step!
    @ C:\Users\vivar002\.julia\packages\DiffEqGPU\I999k\src\ensemblegpukernel\perform_step\gpu_tsit5_perform_step.jl:90
 [13] macro expansion
    @ C:\Users\vivar002\.julia\packages\DiffEqGPU\I999k\src\ensemblegpukernel\kernels.jl:97
 [14] gpu_ode_asolve_kernel
    @ C:\Users\vivar002\.julia\packages\KernelAbstractions\zPAn3\src\macros.jl:95
 [15] gpu_ode_asolve_kernel
    @ .\none:0
Reason: unsupported call to an unknown function (call to julia.push_gc_frame)

Followed by many more stacktrace messages that I can share if needed...

Not sure if @jaswin90 can reproduce the same issue as well?

[ChrisRackauckas]

@utkarsh530 are you able to reproduce this with just Julia?

[utkarsh530]

https://github.com/SciML/diffeqpy?tab=readme-ov-file#benchmark The Julia example here works for me.

[utkarsh530]

@ChrisRackauckas I meant here the code works fine, and I cannot reproduce the error 😅

[ChrisRackauckas]

@gonzalovivares does the Julia code work for you? Maybe it's an installation issue or just using an older version.

[gonzalovivares]

@ChrisRackauckas Yes, running the examples of DiffEqGPU using Julia (via VScode) works fine for me as well. I created a new project and Python environment to reinstall diffeqpy and discard an installation issue (diffEqPy 2.4.1 and Julia 1.10.2), but the error in Python when calling the cuda.EnsembleGPUKernel still shows up... Not sure if any other users are also getting this error in Python...

I have found discussions from other packages reporting the same error Reason: unsupported call to an unknown function (call to julia.new_gc_frame), maybe their solutions might be useful for you?

https://discourse.julialang.org/t/using-mapreduce-on-gpu-with-cuda-jl/108023

CliMA/Oceananigans.jl#3420