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