Hello all!

I’m using DifferentialEquations.jl to solve a stiff DAE. I’ve tried both mass matrix and in-place methods with Rodas4 and IDA respectively at tolerances of 1e-8, and the result is the same (albeit IDA tends to be faster).

This model has a small set of tuning parameters that I would like to optimize. I thought parallel ensemble simulations was the best fit, especially since I have a GPU from NVIDIA. As I’ve read, GPUs work better with Float32 numbers. Additionally, as I understand it the problem now must be in mass matrix form because ensemble simulations require solvers written in Julia (so farewell Sundials.jl).

However, once I change the inputs to Float32 and attempt to solve the system, no solver is able to converge, i.e. I get the InitialFailure retcode everytime. All I did was change the type of float. I’ve tried a few things:

- Different solvers - Rosenbrock23, Rodas4, Rodas4P, Rodas4P2, etc.
- Relaxing the tolerances (1e-4, which is still reasonable)
- Modelingtoolkitize for analytical jacobians and time gradients (works with Float64 and yields the same expected results)

Same results. So I thought I’d stick to Float64 and try the ensemble simulations all the same:

```
using DifferentialEquations, ModelingToolkit, CUDA, DiffEqGPU
f = ODEFunction(model, mass_matrix=Diagonal([1,1,1,1,1,1,0,0]))
prob = ODEProblem(f, x₀, tspan, p)
sys = modelingtoolkitize(prob)
f = ODEFunction(sys, jac=true, tgrad=true, mass_matrix=Diagonal([1,1,1,1,1,1,0,0]))
prob_j = ODEProblem(f, x₀, tspan, p)
ens = EnsembleProblem(prob_j)
sol = solve(ens, Rodas4(), EnsembleGPUArray(CUDA.CUDABackend()), trajectories=10, saveat=1.0f0)
```

Which yields the following error:

```
ERROR: GPU compilation of MethodInstance for DiffEqGPU.gpu_gpu_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}}}}}, ::ModelingToolkit.var"#f#531"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x77563588, 0x499ec9b7, 0x3b31503d, 0xdff47b82, 0xf6f5cef8), Expr}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x41231819, 0x0205e81c, 0xb065ed56, 0x20ebc83e, 0xf4e4193f), Expr}}, ::CUDA.CuDeviceMatrix{Float64, 1}, ::CUDA.CuDeviceMatrix{Float64, 1}, ::CUDA.CuDeviceMatrix{NTuple{11, Float64}, 1}, ::Float64) failed
KernelError: passing and using non-bitstype argument
Argument 3 to your kernel function is of type ModelingToolkit.var"#f#531"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x77563588, 0x499ec9b7, 0x3b31503d, 0xdff47b82, 0xf6f5cef8), Expr}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x41231819, 0x0205e81c, 0xb065ed56, 0x20ebc83e, 0xf4e4193f), Expr}}, which is not isbits:
.f_oop is of type RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x77563588, 0x499ec9b7, 0x3b31503d, 0xdff47b82, 0xf6f5cef8), Expr} which is not isbits.
.body is of type Expr which is not isbits.
.head is of type Symbol which is not isbits.
.args is of type Vector{Any} which is not isbits.
.f_iip is of type RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x41231819, 0x0205e81c, 0xb065ed56, 0x20ebc83e, 0xf4e4193f), Expr} which is not isbits.
.body is of type Expr which is not isbits.
.head is of type Symbol which is not isbits.
.args is of type Vector{Any} which is not isbits.
```

I’m at a loss. Any help is greatly appreciated! However, I cannot provide a MWE. These are the versions of the packages:

CUDA v4.3.0

DiffEqGPU v2.2.1

DifferentialEquations v7.7.0

ModelingToolkit v8.57.0