GPU Kernel Compilation Error

I’m porting some code from my CPU to GPU. I’m obtaining the error that I’ve copied below. The only thing I can think is the problem that my array of initial conditions are scoped incorrectly for problem func. But this is only a guess – I’m struggling to interpret this error. My code is below.

ERROR: InvalidIRError: compiling kernel #gpu_gpu_kernel(KernelAbstractions.CompilerMetadata{KernelAbstractions.NDIt$
ration.DynamicSize, KernelAbstractions.NDIteration.DynamicCheck, Nothing, CartesianIndices{1, Tuple{Base.OneTo{Int6$
}}}, KernelAbstractions.NDIteration.NDRange{1, KernelAbstractions.NDIteration.DynamicSize, KernelAbstractions.NDIte$
ation.DynamicSize, CartesianIndices{1, Tuple{Base.OneTo{Int64}}}, CartesianIndices{1, Tuple{Base.OneTo{Int64}}}}}, t
ypeof(f), CuDeviceMatrix{Float64, 1}, CuDeviceMatrix{Float64, 1}, CuDeviceMatrix{Float64, 1}, Float32) resulted in i
nvalid LLVM IR
Reason: unsupported dynamic function invocation (call to -)
Stacktrace:
 [1] Ik
   @ ./REPL[25]:1
 [2] f
   @ ./REPL[27]:5
 [3] macro expansion
   @ ~/.julia/packages/DiffEqGPU/hWrYl/src/DiffEqGPU.jl:24
 [4] gpu_gpu_kernel
   @ ~/.julia/packages/KernelAbstractions/DqITC/src/macros.jl:81
 [5] gpu_gpu_kernel
   @ ./none:0

using DiffEqGPU
using DifferentialEquations
using CUDA
using GPUInspector

#################################################
##### ODE params

# initialize parameter vector
p = fill(NaN,8);

# Time Constants for Homeostasis -- (g_Ca, g_K, nHalf)
p[1] = 500; p[2] = 500; p[3] = 500;

# C_T & Delta
p[4] = 20; p[5] = 5; 

# Sigmoid Ranges -- (g_Ca, g_K, nHalf)
p[6] = 3; p[7] = 6; p[8] = 20;

# Calcium Handling
A = -1; k = -1/.1;

I_app = 0;
C = 1; g_L = .5;
E_L = -50; E_Ca = 100; E_K = -70;

#################################################
## Initial Condition Testing -- Tests for Multistability

IC_V = -75.0:15.0:75.0;
IC_n = 0:.2:1;

# IC_gCa = 0:p[6]/5:p[6];
# IC_gK = 0:p[7]/5:p[7];
IC_gCa = 0:.3:1.5;
IC_gK = 0:.6:3;

IC_Ca = 0.0:50.0:300.0;
IC_nHalf = 0:p[8]/5:p[8];

AllICs = collect.(Iterators.product(IC_V,IC_n,IC_gCa,IC_gK,IC_Ca,IC_nHalf));

#################################################

# Define functions
sigmoid(x) = 1/(1+exp(-x));
tau(V, halfN) = 3/cosh(((V-halfN)-10)/29);
Ninf(V, halfN) = sigmoid(((V-halfN)-10)/7.25);

# Define Currents
Ica(g_Ca,V) = g_Ca*(sigmoid((V+1)/7.5)+.1)*(V-E_Ca);
Ik(g_K,n,V) = g_K*n*(V-E_K);
Il(V) = g_L*(V-E_L);

# u[1] - V cell 1, u[2] - slow var cell, u[3] - g_Ca, u[4] - g_K, u[5] - [Ca], u[6] - halfN
# p[1] - tau_L.gCa , p[2] - tau_L.gK, p[3] - tau_L.nHalf, p[4] - C_T, p[5] - Delta, p[6] - g_Ca, p[7] - g_K, p[8] - nHalf
function f(du,u,p,t)
    du[1] = (-1*Il(u[1]) + -1*Ica(u[3],u[1]) + -1*Ik(u[4],u[2],u[1]) + I_app)/C
    du[2] = (Ninf(u[1],u[6])-u[2])/tau(u[1],u[6])
    du[3] = (p[6]*sigmoid((p[4]-u[5])/p[5]) - u[3])/p[1]
    du[4] = (p[7]*sigmoid((u[5]-p[4])/p[5]) - u[4])/p[2]
    du[5] = -1*k*(A*Ica(u[3],u[1]) - u[5])
    du[6] = (p[8]*(sigmoid((p[4]-u[5])/p[5])) - u[6])/p[3]
end

# Simulation Parameters
#Trajs = 99792;
Trajs = 997;
SaveTimes = [7000:.1:7150;7150.1];
ICs = [0.0f0 0.0f0 0.0f0 0.0f0 0.0f0 0.0f0];
tspan = (0.0f0,8000.0f0);
prob = ODEProblem(f,ICs,tspan,p);
function prob_func(prob,i,repeat)
  remake(prob,u0=collect(AllICs[i]))
end
ensemble_prob = EnsembleProblem(prob, prob_func = prob_func, safetycopy=false)

monitoring_start()    
@time sim = solve(ensemble_prob,Tsit5(),EnsembleGPUArray(),trajectories=Trajs,saveat=SaveTimes)
results = monitoring_stop();
plot_monitoring_results(results)