CUDA memory error while using DifferentialEquations.jl

Not sure if I could get an answer to this kind of question here, but hopefully someone has an insight.

I am solving a differential equation whose solution is in the form of f(z,ω). For each z, f(z,ω) is some sort of distribution in ω, and this distribution evolves as z increases. The solution is obtained by integrating ∂f/∂z in z for a given initial condition f(0,ω) using DifferentialEquations.jl.

The number of sampling points in ω is very large, so I store f(z,ω) for each z as a CuVector. The integration in z is done for 0 ≤ z ≤ L. The integration works well in general, but for a sufficiently large L that requires a large number of integration steps (typically more than 400,000 steps, but the number changes with problem specs), I get the following “illegal memory access” error:

ERROR: LoadError: CUDA error: an illegal memory access was encountered (code 700, ERROR_ILLEGAL_ADDRESS)
Stacktrace:
  [1] throw_api_error(res::CUDA.cudaError_enum)
    @ CUDA ~/.julia/packages/CUDA/lwSps/lib/cudadrv/error.jl:105
  [2] query
    @ ~/.julia/packages/CUDA/lwSps/lib/cudadrv/stream.jl:102 [inlined]
  [3] synchronize(stream::CuStream; blocking::Bool)
    @ CUDA ~/.julia/packages/CUDA/lwSps/lib/cudadrv/stream.jl:130
  [4] synchronize (repeats 2 times)
    @ ~/.julia/packages/CUDA/lwSps/lib/cudadrv/stream.jl:117 [inlined]
  [5] unsafe_copyto!(dest::Vector{Float64}, doffs::Int64, src::CuArray{Float64, 1}, soffs::Int64, n::Int64)
    @ CUDA ~/.julia/packages/CUDA/lwSps/src/array.jl:356
  [6] copyto!
    @ ~/.julia/packages/CUDA/lwSps/src/array.jl:316 [inlined]
  [7] getindex
    @ ~/.julia/packages/GPUArrays/8dzSJ/src/host/indexing.jl:89 [inlined]
  [8] #24
    @ ~/.julia/packages/GPUArrays/8dzSJ/src/host/indexing.jl:75 [inlined]
  [9] task_local_storage(body::GPUArrays.var"#24#27"{CuArray{Float64, 1}}, key::Symbol, val::Bool)
    @ Base ./task.jl:281
 [10] macro expansion
    @ ~/.julia/packages/GPUArrays/8dzSJ/src/host/indexing.jl:74 [inlined]
 [11] #_mapreduce#21
    @ ~/.julia/packages/GPUArrays/8dzSJ/src/host/mapreduce.jl:65 [inlined]
 [12] #mapreduce#19
    @ ~/.julia/packages/GPUArrays/8dzSJ/src/host/mapreduce.jl:28 [inlined]
 [13] mapreduce
    @ ~/.julia/packages/GPUArrays/8dzSJ/src/host/mapreduce.jl:28 [inlined]
 [14] #_sum#682
    @ ./reducedim.jl:878 [inlined]
 [15] _sum
    @ ./reducedim.jl:878 [inlined]
 [16] #sum#680
    @ ./reducedim.jl:874 [inlined]
 [17] sum
    @ ./reducedim.jl:874 [inlined]
 [18] ODE_DEFAULT_NORM(u::CuArray{ComplexF64, 1}, t::Float64)
    @ DiffEqBase ~/.julia/packages/DiffEqBase/niZxn/src/init.jl:210
 [19] perform_step!(integrator::OrdinaryDiffEq.ODEIntegrator{OrdinaryDiffEq.Tsit5, true, CuArray{ComplexF64, 1}, Nothing, Float64, PulseFD.DiffEqIPDriverParams{CuArray{ComplexF64, 1}, Tuple{DispersionEffect{CuArray{ComplexF64, 1}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{DispersionEffect{CuArray{ComplexF64, 1}}, Nothing, typeof(identity), Float64, FourierManager{Float64, CuArray{Float64, 1}, CuArray{ComplexF64, 1}, CUDA.CUFFT.cCuFFTPlan{ComplexF64, -1, true, 1}}}}, Nothing}, Tuple{Tuple{KerrEffect{Float64, CuArray{Float64, 1}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{KerrEffect{Float64, CuArray{Float64, 1}}, Nothing, typeof(identity), Float64, FourierManager{Float64, CuArray{Float64, 1}, CuArray{ComplexF64, 1}, CUDA.CUFFT.cCuFFTPlan{ComplexF64, -1, true, 1}}}}, Nothing}, Tuple{RamanEffect{Float64, CuArray{Float64, 1}, CuArray{ComplexF64, 1}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{RamanEffect{Float64, CuArray{Float64, 1}, CuArray{ComplexF64, 1}}, Nothing, typeof(identity), Float64, FourierManager{Float64, CuArray{Float64, 1}, CuArray{ComplexF64, 1}, CUDA.CUFFT.cCuFFTPlan{ComplexF64, -1, true, 1}}}}, Nothing}}, typeof(identity), FourierManager{Float64, CuArray{Float64, 1}, CuArray{ComplexF64, 1}, CUDA.CUFFT.cCuFFTPlan{ComplexF64, -1, true, 1}}, PulseFD.MonitorVars{Float64}}, Float64, Float64, Float64, Float64, Vector{CuArray{ComplexF64, 1}}, SciMLBase.ODESolution{ComplexF64, 2, Vector{CuArray{ComplexF64, 1}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{CuArray{ComplexF64, 1}}}, SciMLBase.ODEProblem{CuArray{ComplexF64, 1}, Tuple{Float64, Float64}, true, PulseFD.DiffEqIPDriverParams{CuArray{ComplexF64, 1}, Tuple{DispersionEffect{CuArray{ComplexF64, 1}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{DispersionEffect{CuArray{ComplexF64, 1}}, Nothing, typeof(identity), Float64, FourierManager{Float64, CuArray{Float64, 1}, CuArray{ComplexF64, 1}, CUDA.CUFFT.cCuFFTPlan{ComplexF64, -1, true, 1}}}}, Nothing}, Tuple{Tuple{KerrEffect{Float64, CuArray{Float64, 1}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{KerrEffect{Float64, CuArray{Float64, 1}}, Nothing, typeof(identity), Float64, FourierManager{Float64, CuArray{Float64, 1}, CuArray{ComplexF64, 1}, CUDA.CUFFT.cCuFFTPlan{ComplexF64, -1, true, 1}}}}, Nothing}, Tuple{RamanEffect{Float64, CuArray{Float64, 1}, CuArray{ComplexF64, 1}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{RamanEffect{Float64, CuArray{Float64, 1}, CuArray{ComplexF64, 1}}, Nothing, typeof(identity), Float64, FourierManager{Float64, CuArray{Float64, 1}, CuArray{ComplexF64, 1}, CUDA.CUFFT.cCuFFTPlan{ComplexF64, -1, true, 1}}}}, Nothing}}, typeof(identity), FourierManager{Float64, CuArray{Float64, 1}, CuArray{ComplexF64, 1}, CUDA.CUFFT.cCuFFTPlan{ComplexF64, -1, true, 1}}, PulseFD.MonitorVars{Float64}}, SciMLBase.ODEFunction{true, typeof(PulseFD.calc_driver!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, OrdinaryDiffEq.Tsit5, OrdinaryDiffEq.InterpolationData{SciMLBase.ODEFunction{true, typeof(PulseFD.calc_driver!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{CuArray{ComplexF64, 1}}, Vector{Float64}, Vector{Vector{CuArray{ComplexF64, 1}}}, OrdinaryDiffEq.Tsit5Cache{CuArray{ComplexF64, 1}, CuArray{ComplexF64, 1}, CuArray{ComplexF64, 1}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}}}, DiffEqBase.DEStats}, SciMLBase.ODEFunction{true, typeof(PulseFD.calc_driver!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, OrdinaryDiffEq.Tsit5Cache{CuArray{ComplexF64, 1}, CuArray{ComplexF64, 1}, CuArray{ComplexF64, 1}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}}, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, Float64, OrdinaryDiffEq.PIController{Rational{Int64}}, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(LinearAlgebra.opnorm), Nothing, DiffEqBase.CallbackSet{Tuple{}, Tuple{DiffEqBase.DiscreteCallback{typeof(PulseFD.cond_always), typeof(PulseFD.update_mv!), typeof(DiffEqBase.INITIALIZE_DEFAULT), typeof(DiffEqBase.FINALIZE_DEFAULT)}, DiffEqBase.DiscreteCallback{typeof(PulseFD.cond_always), typeof(PulseFD.print_mv), typeof(DiffEqBase.INITIALIZE_DEFAULT), typeof(DiffEqBase.FINALIZE_DEFAULT)}}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryMinHeap{Float64}, DataStructures.BinaryMinHeap{Float64}, Nothing, Nothing, Int64, Tuple{}, Vector{Float64}, Tuple{}}, CuArray{ComplexF64, 1}, ComplexF64, Nothing, OrdinaryDiffEq.DefaultInit}, cache::OrdinaryDiffEq.Tsit5Cache{CuArray{ComplexF64, 1}, CuArray{ComplexF64, 1}, CuArray{ComplexF64, 1}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}}, repeat_step::Bool)
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/9mxZY/src/perform_step/low_order_rk_perform_step.jl:658
 [20] perform_step!
    @ ~/.julia/packages/OrdinaryDiffEq/9mxZY/src/perform_step/low_order_rk_perform_step.jl:628 [inlined]
 [21] solve!(integrator::OrdinaryDiffEq.ODEIntegrator{OrdinaryDiffEq.Tsit5, true, CuArray{ComplexF64, 1}, Nothing, Float64, PulseFD.DiffEqIPDriverParams{CuArray{ComplexF64, 1}, Tuple{DispersionEffect{CuArray{ComplexF64, 1}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{DispersionEffect{CuArray{ComplexF64, 1}}, Nothing, typeof(identity), Float64, FourierManager{Float64, CuArray{Float64, 1}, CuArray{ComplexF64, 1}, CUDA.CUFFT.cCuFFTPlan{ComplexF64, -1, true, 1}}}}, Nothing}, Tuple{Tuple{KerrEffect{Float64, CuArray{Float64, 1}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{KerrEffect{Float64, CuArray{Float64, 1}}, Nothing, typeof(identity), Float64, FourierManager{Float64, CuArray{Float64, 1}, CuArray{ComplexF64, 1}, CUDA.CUFFT.cCuFFTPlan{ComplexF64, -1, true, 1}}}}, Nothing}, Tuple{RamanEffect{Float64, CuArray{Float64, 1}, CuArray{ComplexF64, 1}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{RamanEffect{Float64, CuArray{Float64, 1}, CuArray{ComplexF64, 1}}, Nothing, typeof(identity), Float64, FourierManager{Float64, CuArray{Float64, 1}, CuArray{ComplexF64, 1}, CUDA.CUFFT.cCuFFTPlan{ComplexF64, -1, true, 1}}}}, Nothing}}, typeof(identity), FourierManager{Float64, CuArray{Float64, 1}, CuArray{ComplexF64, 1}, CUDA.CUFFT.cCuFFTPlan{ComplexF64, -1, true, 1}}, PulseFD.MonitorVars{Float64}}, Float64, Float64, Float64, Float64, Vector{CuArray{ComplexF64, 1}}, SciMLBase.ODESolution{ComplexF64, 2, Vector{CuArray{ComplexF64, 1}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{CuArray{ComplexF64, 1}}}, SciMLBase.ODEProblem{CuArray{ComplexF64, 1}, Tuple{Float64, Float64}, true, PulseFD.DiffEqIPDriverParams{CuArray{ComplexF64, 1}, Tuple{DispersionEffect{CuArray{ComplexF64, 1}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{DispersionEffect{CuArray{ComplexF64, 1}}, Nothing, typeof(identity), Float64, FourierManager{Float64, CuArray{Float64, 1}, CuArray{ComplexF64, 1}, CUDA.CUFFT.cCuFFTPlan{ComplexF64, -1, true, 1}}}}, Nothing}, Tuple{Tuple{KerrEffect{Float64, CuArray{Float64, 1}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{KerrEffect{Float64, CuArray{Float64, 1}}, Nothing, typeof(identity), Float64, FourierManager{Float64, CuArray{Float64, 1}, CuArray{ComplexF64, 1}, CUDA.CUFFT.cCuFFTPlan{ComplexF64, -1, true, 1}}}}, Nothing}, Tuple{RamanEffect{Float64, CuArray{Float64, 1}, CuArray{ComplexF64, 1}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{RamanEffect{Float64, CuArray{Float64, 1}, CuArray{ComplexF64, 1}}, Nothing, typeof(identity), Float64, FourierManager{Float64, CuArray{Float64, 1}, CuArray{ComplexF64, 1}, CUDA.CUFFT.cCuFFTPlan{ComplexF64, -1, true, 1}}}}, Nothing}}, typeof(identity), FourierManager{Float64, CuArray{Float64, 1}, CuArray{ComplexF64, 1}, CUDA.CUFFT.cCuFFTPlan{ComplexF64, -1, true, 1}}, PulseFD.MonitorVars{Float64}}, SciMLBase.ODEFunction{true, typeof(PulseFD.calc_driver!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, OrdinaryDiffEq.Tsit5, OrdinaryDiffEq.InterpolationData{SciMLBase.ODEFunction{true, typeof(PulseFD.calc_driver!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{CuArray{ComplexF64, 1}}, Vector{Float64}, Vector{Vector{CuArray{ComplexF64, 1}}}, OrdinaryDiffEq.Tsit5Cache{CuArray{ComplexF64, 1}, CuArray{ComplexF64, 1}, CuArray{ComplexF64, 1}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}}}, DiffEqBase.DEStats}, SciMLBase.ODEFunction{true, typeof(PulseFD.calc_driver!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, OrdinaryDiffEq.Tsit5Cache{CuArray{ComplexF64, 1}, CuArray{ComplexF64, 1}, CuArray{ComplexF64, 1}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}}, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, Float64, OrdinaryDiffEq.PIController{Rational{Int64}}, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(LinearAlgebra.opnorm), Nothing, DiffEqBase.CallbackSet{Tuple{}, Tuple{DiffEqBase.DiscreteCallback{typeof(PulseFD.cond_always), typeof(PulseFD.update_mv!), typeof(DiffEqBase.INITIALIZE_DEFAULT), typeof(DiffEqBase.FINALIZE_DEFAULT)}, DiffEqBase.DiscreteCallback{typeof(PulseFD.cond_always), typeof(PulseFD.print_mv), typeof(DiffEqBase.INITIALIZE_DEFAULT), typeof(DiffEqBase.FINALIZE_DEFAULT)}}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryMinHeap{Float64}, DataStructures.BinaryMinHeap{Float64}, Nothing, Nothing, Int64, Tuple{}, Vector{Float64}, Tuple{}}, CuArray{ComplexF64, 1}, ComplexF64, Nothing, OrdinaryDiffEq.DefaultInit})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/9mxZY/src/solve.jl:478
 [22] #__solve#467
    @ ~/.julia/packages/OrdinaryDiffEq/9mxZY/src/solve.jl:5 [inlined]
 [23] __solve(::SciMLBase.ODEProblem{CuArray{ComplexF64, 1}, Tuple{Float64, Float64}, true, PulseFD.DiffEqIPDriverParams{CuArray{ComplexF64, 1}, Tuple{DispersionEffect{CuArray{ComplexF64, 1}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{DispersionEffect{CuArray{ComplexF64, 1}}, Nothing, typeof(identity), Float64, FourierManager{Float64, CuArray{Float64, 1}, CuArray{ComplexF64, 1}, CUDA.CUFFT.cCuFFTPlan{ComplexF64, -1, true, 1}}}}, Nothing}, Tuple{Tuple{KerrEffect{Float64, CuArray{Float64, 1}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{KerrEffect{Float64, CuArray{Float64, 1}}, Nothing, typeof(identity), Float64, FourierManager{Float64, CuArray{Float64, 1}, CuArray{ComplexF64, 1}, CUDA.CUFFT.cCuFFTPlan{ComplexF64, -1, true, 1}}}}, Nothing}, Tuple{RamanEffect{Float64, CuArray{Float64, 1}, CuArray{ComplexF64, 1}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{RamanEffect{Float64, CuArray{Float64, 1}, CuArray{ComplexF64, 1}}, Nothing, typeof(identity), Float64, FourierManager{Float64, CuArray{Float64, 1}, CuArray{ComplexF64, 1}, CUDA.CUFFT.cCuFFTPlan{ComplexF64, -1, true, 1}}}}, Nothing}}, typeof(identity), FourierManager{Float64, CuArray{Float64, 1}, CuArray{ComplexF64, 1}, CUDA.CUFFT.cCuFFTPlan{ComplexF64, -1, true, 1}}, PulseFD.MonitorVars{Float64}}, SciMLBase.ODEFunction{true, typeof(PulseFD.calc_driver!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, ::Nothing; default_set::Bool, kwargs::Base.Iterators.Pairs{Symbol, Any, NTuple{5, Symbol}, NamedTuple{(:second_time, :save_everystep, :saveat, :abstol, :callback), Tuple{Bool, Bool, Vector{Float64}, Float64, DiffEqBase.CallbackSet{Tuple{}, Tuple{DiffEqBase.DiscreteCallback{typeof(PulseFD.cond_always), typeof(PulseFD.update_mv!), typeof(DiffEqBase.INITIALIZE_DEFAULT), typeof(DiffEqBase.FINALIZE_DEFAULT)}, DiffEqBase.DiscreteCallback{typeof(PulseFD.cond_always), typeof(PulseFD.print_mv), typeof(DiffEqBase.INITIALIZE_DEFAULT), typeof(DiffEqBase.FINALIZE_DEFAULT)}}}}}})
    @ DifferentialEquations ~/.julia/packages/DifferentialEquations/el96s/src/default_solve.jl:7
 [24] #__solve#73
    @ ~/.julia/packages/DiffEqBase/niZxn/src/solve.jl:285 [inlined]
 [25] #solve_call#58
    @ ~/.julia/packages/DiffEqBase/niZxn/src/solve.jl:61 [inlined]
 [26] #solve_up#60
    @ ~/.julia/packages/DiffEqBase/niZxn/src/solve.jl:88 [inlined]
 [27] solve(::SciMLBase.ODEProblem{CuArray{ComplexF64, 1}, Tuple{Float64, Float64}, true, PulseFD.DiffEqIPDriverParams{CuArray{ComplexF64, 1}, Tuple{DispersionEffect{CuArray{ComplexF64, 1}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{DispersionEffect{CuArray{ComplexF64, 1}}, Nothing, typeof(identity), Float64, FourierManager{Float64, CuArray{Float64, 1}, CuArray{ComplexF64, 1}, CUDA.CUFFT.cCuFFTPlan{ComplexF64, -1, true, 1}}}}, Nothing}, Tuple{Tuple{KerrEffect{Float64, CuArray{Float64, 1}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{KerrEffect{Float64, CuArray{Float64, 1}}, Nothing, typeof(identity), Float64, FourierManager{Float64, CuArray{Float64, 1}, CuArray{ComplexF64, 1}, CUDA.CUFFT.cCuFFTPlan{ComplexF64, -1, true, 1}}}}, Nothing}, Tuple{RamanEffect{Float64, CuArray{Float64, 1}, CuArray{ComplexF64, 1}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{RamanEffect{Float64, CuArray{Float64, 1}, CuArray{ComplexF64, 1}}, Nothing, typeof(identity), Float64, FourierManager{Float64, CuArray{Float64, 1}, CuArray{ComplexF64, 1}, CUDA.CUFFT.cCuFFTPlan{ComplexF64, -1, true, 1}}}}, Nothing}}, typeof(identity), FourierManager{Float64, CuArray{Float64, 1}, CuArray{ComplexF64, 1}, CUDA.CUFFT.cCuFFTPlan{ComplexF64, -1, true, 1}}, PulseFD.MonitorVars{Float64}}, SciMLBase.ODEFunction{true, typeof(PulseFD.calc_driver!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}; sensealg::Nothing, u0::Nothing, p::Nothing, kwargs::Base.Iterators.Pairs{Symbol, Any, NTuple{4, Symbol}, NamedTuple{(:save_everystep, :saveat, :abstol, :callback), Tuple{Bool, Vector{Float64}, Float64, DiffEqBase.CallbackSet{Tuple{}, Tuple{DiffEqBase.DiscreteCallback{typeof(PulseFD.cond_always), typeof(PulseFD.update_mv!), typeof(DiffEqBase.INITIALIZE_DEFAULT), typeof(DiffEqBase.FINALIZE_DEFAULT)}, DiffEqBase.DiscreteCallback{typeof(PulseFD.cond_always), typeof(PulseFD.print_mv), typeof(DiffEqBase.INITIALIZE_DEFAULT), typeof(DiffEqBase.FINALIZE_DEFAULT)}}}}}})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/niZxn/src/solve.jl:73

Because the code works well for a fairly large number of integration steps, I suspect this might be a bug in either CUDA.jl or DifferentialEquations.jl rather than in my own code. Has anyone had a similar error, or does anyone have any insights on this problem?

I’ll need an example for this. This seems, deep :sweat_smile:.

@ChrisRackauckas, while trying to create MWE, I found out that it was a CUDA error. I filed an issue at CUDA.jl: Illegal memory access during complex exponential with large imaginary part as exponent · Issue #1085 · JuliaGPU/CUDA.jl · GitHub.