# Using DiffEqGPU for very large systems: size app 1e7-1e8

**URL:** <https://discourse.julialang.org/t/using-diffeqgpu-for-very-large-systems-size-app-1e7-1e8/134718>\
**Category:** Modelling & Simulations\
**Tags:** diffeq\
**Created:** [December 24, 2025, 6:22am UTC](https://discourse.julialang.org/t/using-diffeqgpu-for-very-large-systems-size-app-1e7-1e8/134718 "2025-12-24T06:22:22Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![bpeters](https://avatars.discourse-cdn.com/v4/letter/b/ec9cab/32.png) [@bpeters](https://discourse.julialang.org/u/bpeters)\
**Post date:** [December 24, 2025, 6:22am UTC](https://discourse.julialang.org/t/using-diffeqgpu-for-very-large-systems-size-app-1e7-1e8/134718/1 "2025-12-24T06:22:22Z")

</div>

@ChrisRackauckas

Dear Chris,

I see a huge benefit using DiffEqGPU for my applications (luxdem.uni.lu), in particular to use the power of GPUs of different architectures. My application comes from thermodynamics and requires to solve the one-dimensional and transient heat equation (diffusion type differential equation). I have already seen MoL, however, I am not able to use it due to conservation issues and working in different frames of references. However, that is not the issue. Thus, my ode system comes from the discretisation of the one-dimensional heat equation. My application requires to solve app. n=1e6 small and similar ode-systems coming from individual particles. Each small ode system representing a particle consisting of app 10 - 100 unknowns with individual boundary conditions. I am not clear if I should use  
EnsembleGPUKernel for one big ODE system or EnsembleGPUArray to handle each small ODE system. From my understanding, EnsembleGPUKernel should be the method to be applied because

- the function f to evaluate the system i.e. u and would be a loop over all small ode system for which transport properties are dependent on u

- the boundary conditions for each small ode system come from the parameters defined in advance before solving

- is able to account for interactions between the small ode systems

That would then require to run EnsembleGPUKernel for one trajectory. Trying this with the following example taken from the homepage and modified the trajectory number:

````julia-auto
```julia
using DiffEqGPU, OrdinaryDiffEq, StaticArrays, CUDA

function lorenz(u, p, t)
σ = p[1]
ρ = p[2]
β = p[3]
du1 = σ * (u[2] - u[1])
du2 = u[1] * (ρ - u[3]) - u[2]
du3 = u[1] * u[2] - β * u[3]
return SVector{3}(du1, du2, du3)
end

u0 = _at_SVector [1.0f0; 0.0f0; 0.0f0]
tspan = (0.0f0, 10.0f0)
p = _at_SVector [10.0f0, 28.0f0, 8 / 3.0f0]
prob = ODEProblem{false}(lorenz, u0, tspan, p)
prob_func = (prob, i, repeat) → remake(prob, p = (_at_SVector rand(Float32, 3)) .* p)
monteprob = EnsembleProblem(prob, prob_func = prob_func, safetycopy = false)
sol = solve(monteprob, GPUTsit5(), EnsembleGPUKernel(CUDA.CUDABackend()), trajectories = 1)
'''

````

gives me (very surprising):

julia\> include(“runODETestsFloat32GPUEnsembleGPUKernel.jl”)  
ERROR: LoadError: MethodError: Cannot `convert` an object of type GPUTsit5 to an object of type Vector{SVector{3, Float32}}  
The function `convert` exists, but no method is defined for this combination of argument types.

Closest candidates are:  
convert(::Type{Array{T, N}}, ::SizedArray{S, T, N, N, Array{T, N}}) where {S, T, N}  
\_at\_ StaticArrays ~/.julia/packages/StaticArrays/DsPgf/src/SizedArray.jl:88  
convert(::Type{Array{T, N}}, ::SizedArray{S, T, N, M, TData} where {M, TData\<:AbstractArray{T, M}}) where {T, S, N}  
\_at\_ StaticArrays ~/.julia/packages/StaticArrays/DsPgf/src/SizedArray.jl:82  
convert(::Type{Array{S, N}}, ::PooledArrays.PooledArray{T, R, N}) where {S, T, R, N}  
\_at\_ PooledArrays ~/.julia/packages/PooledArrays/Vy2X0/src/PooledArrays.jl:499  
…

Stacktrace:  
[1] \_\_init(prob::ODEProblem{…}, alg::CompositeAlgorithm{…}, timeseries\_init::GPUTsit5, ts\_init::Tuple{}, ks\_init::Tuple{}; saveat::Tuple{}, tstops::Tuple{}, d\_discontinuities::Tuple{}, save\_idxs::Nothing, save\_everystep::Bool, save\_on::Bool, save\_discretes::Bool, save\_start::Bool, save\_end::Nothing, callback::Nothing, dense::Bool, calck::Bool, dt::Float32, dtmin::Float32, dtmax::Float32, force\_dtmin::Bool, adaptive::Bool, gamma::Rational{…}, abstol::Nothing, reltol::Nothing, qmin::Rational{…}, qmax::Int64, qsteady\_min::Int64, qsteady\_max::Int64, beta1::Nothing, beta2::Nothing, qoldinit::Rational{…}, controller::Nothing, fullnormalize::Bool, failfactor::Int64, maxiters::Int64, internalnorm::typeof(DiffEqBase.ODE\_DEFAULT\_NORM), internalopnorm::typeof(LinearAlgebra.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::ODEAliasSpecifier, initializealg::DiffEqBase.DefaultInit, kwargs::\_at\_KwargsKwargsKwargsKwargsKwargsKwargsKwargsKwargs{…})  
\_at\_ OrdinaryDiffEqCore ~/.julia/packages/OrdinaryDiffEqCore/GMkz9/src/solve.jl:334  
[2] \_\_solve(prob::ODEProblem{…}, alg::CompositeAlgorithm{…}, args::GPUTsi\_at\_Kwargs5; kw\_at\_KwargsKwargsrgs::\_at\_Kwargs{…})  
\_at\_ OrdinaryDiffEqCore ~/.julia/packages/OrdinaryDiffEqCore/GMkz9/src/solve.jl:6  
[3] \_\_solve(prob::ODEProblem{…}, ::Nothing, args::\_at\_KwargsPUTsi\_at\_KwargsKwargs5; kwargs::\_at\_Kwargs{…})  
\_at\_ OrdinaryDiffEqDefault ~/.julia/packages/OrdinaryDiffEqDefault/MdlB6/src/default\_alg.jl:48  
[4] \_\_solve(prob::ODEProblem{…}, args::GPUTsit5; default\_set::Bool, sec\_at\_Kwargsnd\_ti\_at\_KwargsKwargse::Bool, kwargs::\_at\_Kwargs{})  
\_at\_ DiffEqBase ~/.julia/packages/DiffEqBase/6sydR/src/solve.jl:758  
[5] \_\_solve(prob::ODEProblem{…}, args::GPUTsit5)  
\_at\_ DiffEqBase ~/.julia/packages/DiffEqBase/6sydR/src/solve.jl:749  
[6] solve\_call(\_prob::ODEProblem{…}, args::GPUTsit5; merge\_callbacks::Bool, k\_at\_Kwargsargsh\_at\_KwargsKwargsndle::Nothing, kwargs::\_at\_Kwargs{})  
\_at\_ DiffEqBase ~/.julia/packages/DiffEqBase/6sydR/src/solve.jl:142  
[7] solve\_call(\_prob::ODEProblem{…}, args::GPUTsit5)  
\_at\_ DiffEqBase ~/.julia/packages/DiffEqBase/6sydR/src/solve.jl:109  
[8] solve\_up(prob::ODEProblem{…}, sensealg::Nothing, u0::SVector{…}, p::SVector{…}, ar\_at\_KwargshainRulesOri\_at\_KwargshainRulesOriginatorinators::G\_at\_KwargsUTsit\_at\_Kwargs; originator::SciMLBase.\_at\_KwargshainRulesOriginator, kwargs::\_at\_Kwargs{})  
\_at\_ DiffEqBase ~/.julia/packages/DiffEqBase/6sydR/src/solve.jl:578  
[9] solve(prob::ODEProblem{…}, args::GPUTsit5;\_at\_Kwargssense\_at\_Kwargslg::Nothing, u0::Nothing\_at\_Kwargs p::Nothing, wrap::Val{…}, kwargs::\_at\_Kwargs{})  
\_at\_ DiffEqBase ~/.julia/packages/DiffEqBase/6sydR/src/solve.jl:545  
[10] solve(prob::ODEProblem{…}, args::GPUTsit5)  
\_at\_ DiffEqBase ~/.julia/packages/DiffEqBase/6sydR/src/solve.jl:535  
[1\_at\_Kwargs] bat\_at\_Kwargsh\_func(i::Int64, prob::E\_at\_KwargssembleProblem{…}, alg::GPUTsit5; kwargs::\_at\_Kwargs{})  
\_at\_ SciMLBase ~/.julia/packages/SciMLBase/hHrRi/src/ensemble/basic\_ensemble\_solve.jl:235  
[12] batch\_func  
\_at\_ ~/.julia/packages/SciMLBase/hHrRi/src/ensemble/basic\_ensemble\_solve.jl:222 [inlined]  
[13] #811  
\_at\_ ~/.julia/packages/SciMLBase/hHrRi/src/ensemble/basic\_ensemble\_solve.jl:294 [inlined]  
[14] responsible\_map  
\_at\_ ~/.julia/packages/SciMLBase/hHrRi/src/ensemble/basic\_ensemble\_solve.jl:287 [inlined]  
[15] solve\_batch(prob::EnsembleProblem{…}, alg::GPUTsit5, ::Ensem\_at\_KwargsKwargsleSerial, II::Unit\_at\_Kwargsange{…}, pmap\_batch\_size::Int64; kwarg\_at\_Kwarg\_at\_Kwarg\_at\_Kwarg\_at\_Kwargs::\_at\_Kwargs{})  
\_at\_ SciMLBase ~/.julia/packages/SciMLBase/hHrRi/src/ensemble/basic\_ensemble\_solve.jl:293  
[16] solve\_batch  
\_at\_ ~/.julia/packages/SciMLBase/hHrRi/src/ensemble/basic\_e\_at\_Kwarg\_at\_Kwargssemble\_solve.jl:292 [inlined]  
[17] macro expansion  
\_at\_ ./timi\_at\_Kwarg\_at\_Kwargsg.jl:461 [inlined]  
[1\_at\_Kwarg\_at\_Kwargs] \_at\_Kwargs::SciMLBase.var"#800#801"{Int64, Int64, I\_at\_Kwarg\_at\_Kwarg\_at\_Kwarg\_at\_Kwargst64, \_at\_Kwargs{}, EnsembleProblem{…}, GPUTsit5, EnsembleSerial})()  
\_at\_ SciMLBase ~/.julia/packages/SciMLBase/hHrRi/src/ensemble/basic\_ensemble\_solve.jl:188  
[19] with\_logstate(f::SciMLBase.var"#800#801"{…}, logstate::Base.CoreLogging.LogState)  
\_at\_ Base.CoreLogging ./logging/logging.jl:540  
[20] with\_logger  
\_at\_ ./logging/logging.jl:651 [inlined]  
[21] \_\_solve(prob::EnsemblePro\_at\_KwargsInt64lem{…}, alg::GPUTsit5, ens\_at\_KwargsInt64mblealg::EnsembleSerial; traje\_at\_Kwargsnt64t\_at\_Kwargskwargsri\_at\_Kwargss::Int64, batch\_si\_at\_Kw\_at\_Kwargsnt64r\_at\_Kwargskwargs\_at\_K\_at\_KwargsargsInt6\_at\_KwargsInt64kwargse:\_at\_KwargsInt64, progress\_aggregate::Bool, pmap\_batch\_siz\_at\_Kwarg\_at\_Kwargsnt64k\_at\_Kwargskwargsar\_at\_\_at\_Kwargsnt64w\_at\_Kwargskwargsrg\_at\_Kwargss::\_at\_Kwargsnt64,\_at\_Kwargskwargs::\_at\_Kwargs{})  
\_at\_ SciMLBase ~/.julia/packages/SciMLBase/hHrRi/src/ensemble/basic\_ensemble\_solve.jl:172  
[22] \_\_solve  
\_at\_ ~/.julia/packages/SciMLBase/hHrRi/src/ensemble/basic\_ensemble\_solve.jl:164 [inlined]  
[23] \_\_solve(ensembleprob::EnsembleProblem{…}, alg::GPUTsit5, ense\_at\_KwargsBoolb\_at\_Kwargskwargsea\_at\_Kwargsg::Ensembl\_at\_KwargsBoolG\_at\_KwargskwargsUK\_at\_Kwargsrnel{…}; tra\_at\_Kwargskwargsec\_at\_Kwargsories::Int64, batch\_size::Int64, unstable\_check::Function, adapti\_at\_Kwarg\_at\_KwargsBoolk\_at\_Kwargskwargsar\_at\_\_at\_KwargsBoolw\_at\_Kwargskwargsrg\_at\_Kwargsse:\_at\_KwargsBool,\_at\_Kwargskwargs::\_at\_Kwargs{})  
\_at\_ DiffEqGPU ~/.julia/packages/DiffEqGPU/JsBHq/src/solve.jl:10  
[24] \_\_solve  
\_at\_ ~/.julia/packages/DiffEqGPU/JsBHq/src/solve.jl:1 [inlined]  
[25] #solve#825  
\_at\_ ~/.julia/packages/SciMLBase/hHrRi/src/ensemble/basic\_ensemble\_solve.jl:359 [inlined]  
[26] top-level scope  
\_at\_ ~/xdemExaScale/xdem/odeSolver/test/runODETestsFloat32GPUEnsembleGPUKernel.jl:19  
[27] include(mapexpr::Function, mod::Module, \_path::String)  
\_at\_ Base ./Base.jl:307  
[28] top-level scope  
\_at\_ REPL[3]:1  
in expression starting at /home/bernhard/xdemExaScale/xdem/odeSolver/test/runODETestsFloat32GPUEnsembleGPUKernel.jl:19  
Some type information was truncated. Use `show(err)` to see complete types.

I hope I made myself clear and appreciate very much your advice.

Best

Bernhard

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [January 21, 2026, 11:34am UTC](https://discourse.julialang.org/t/using-diffeqgpu-for-very-large-systems-size-app-1e7-1e8/134718/2 "2026-01-21T11:34:14Z")

</div>

For this kind of problem it would be much more effective to write out the PDE and use the array form rather than the ensemble form. See the tutorial that is PDE focused

> **[GPU-Accelerated Stochastic Partial Differential Equations · Overview of...](https://docs.sciml.ai/Overview/stable/showcase/gpu_spde/)**
>
> Documentation for Overview of Julia's SciML.

---

<div class="post-metadata">

**Author:** ![bpeters](https://avatars.discourse-cdn.com/v4/letter/b/ec9cab/32.png) [@bpeters](https://discourse.julialang.org/u/bpeters)\
**Post date:** [January 23, 2026, 2:10pm UTC](https://discourse.julialang.org/t/using-diffeqgpu-for-very-large-systems-size-app-1e7-1e8/134718/3 "2026-01-23T14:10:10Z")

</div>

Thanks for the confirmation. Without diving deeply into the package, I used that more or less intuitively.

Best

BP