Solving DAE in Julia using GPUs

Question

I'm trying to solve a Differential Algebraic Equation (DAE) in Julia which is very computationally expensive using GPUs. I'm brand new to Julia and don't have much experience coding with GPUs. The below problem is just a sample DAE problem I'm trying to solve. The error I'm getting, which occur on the last line of code, is

MethodError: no method matching generate_problem(::DAEProblem{Vector{Float32}, Vector{Float32}, Tuple{Float32, Float32}, true, Vector{Float32}, DAEFunction{true, typeof(f1), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, Vector{Float64}}, ::CuArray{Float32, 2, CUDA.Mem.DeviceBuffer}, ::CuArray{Float32, 2, CUDA.Mem.DeviceBuffer}, ::Nothing, ::Nothing)

The Code that is generating the error is below and based on an examples from

https://github.com/SciML/DiffEqGPU.jl ### This example I can get to run

I've been able to solve DAE's without using GPU's using the code from the below link https://docs.juliahub.com/DifferentialEquations/UQdwS/6.15.0/tutorials/dae_example/

Any help is greatly appreciated

Julia Code:

using Pkg
using Sundials
using DifferentialEquations
using CSV
using DataFrames
using CUDA, LinearAlgebra
using DiffEqGPU



function f1(out,du,u,p,t)
    out[1] = p[1]*(u[2]-u[1]) - du[1]
    out[2] = u[1]*(p[2]-u[3]) - u[2] - du[2]
    out[3] = u[1]*u[2] - p[3]*u[3] - du[3]
end

u₀= Float32[1.0;0.0;0.0]

du₀ = Float32[0.5;0.5;0.5]

tspan = (0.0f0,100.0f0)
p = [10.0f0,28.0f0,8/3f0]

#differential_vars = Float32.(ones(3))

differential_vars = ones(3)



prob = DAEProblem(f1,du₀,u₀,tspan,p,differential_vars = differential_vars)


prob_func = (prob,i,repeat) -> remake(prob,p=rand(Float32,3).*p)
monteprob = EnsembleProblem(prob, prob_func = prob_func, safetycopy=false)

### The below line is where the error occurs
@time sol = solve(monteprob,Tsit5(),EnsembleGPUArray(),trajectories=10_000,saveat=1.0f0)

Answer 1

Note that you need to use a native DAE solver if you're going to use it with DiffEqGPU, and you'll need a release after this merges. When that is all together, DFBDF will be the algorithm that is compatible and the correct code would be:

using OrdinaryDiffEq
using DiffEqGPU

function f1(out,du,u,p,t)
    out[1] = p[1]*(u[2]-u[1]) - du[1]
    out[2] = u[1]*(p[2]-u[3]) - u[2] - du[2]
    out[3] = u[1]*u[2] - p[3]*u[3] - du[3]
end

u₀= Float32[1.0;0.0;0.0]
du₀ = Float32[0.5;0.5;0.5]
tspan = (0.0f0,100.0f0)
p = [10.0f0,28.0f0,8/3f0]
differential_vars = ones(3)

prob = DAEProblem(f1,du₀,u₀,tspan,p,differential_vars = differential_vars)
prob_func = (prob,i,repeat) -> remake(prob,p=rand(Float32,3).*p)
monteprob = EnsembleProblem(prob, prob_func = prob_func, safetycopy=false)

@time sol = solve(monteprob,DFBDF(),EnsembleGPUArray(),trajectories=1_000,saveat=1.0f0)

That said, right now DFBDF is not optimized and so it's unlikely this will be useful until late 2022.