# Solving DAE in Julia using GPUs

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)


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.

• Chris thankyou, do you know of an approach I could use in the meantime, or when the merge would take place? Dec 16, 2021 at 19:46
• Even if that is merged, I'll say it's not useful yet because DFBDF will need to get optimized first. What's more fruitful would be to use ModelingToolkit.jl with structural_simplify to lower the DAEs into ODEs which can more easily be GPU accelerated. In fact, that's a faster way to solve them anyways according to the benchmarks. See for example this benchmark. Dec 16, 2021 at 22:53