Drop initialization and get tutorial finished

Author: ChrisRackauckas

docs/src/tutorials/modelingtoolkit.md

@@ -7,6 +7,11 @@ the simplest set of equations to solve and exploiting things that normally canno
 by hand. Those exact features are also potentially useful for GPU computing, and thus this
 tutorial showcases how to effectively use MTK with DiffEqGPU.jl.
 
+!!! warn
+    This tutorial currently only works for ODEs defined by ModelingToolkit. More work
+    will be required to support DAEs in full. This is work that is ongoing and expected
+    to be completed by the summer of 2025.
+
 The core aspect to doing this right is two things. First of all, MTK respects the types
 chosen by the user, and thus in order for GPU kernel generation to work the user needs
 to ensure that the problem that is built uses static structures. For example this means
@@ -24,7 +29,6 @@ eqs = [D(D(x)) ~ σ * (y - x),
     D(z) ~ x * y - β * z]
 
 @mtkbuild sys = ODESystem(eqs, t)
-
 u0 = SA[D(x) => 2f0,
     x => 1f0,
     y => 0f0,
@@ -72,6 +76,7 @@ function prob_func2(prob, i, repeat)
     u0, p = sym_setter(prob,@SVector(rand(Float32,3)))
     remake(prob, u0 = u0, p = p)
 end
+
 monteprob = EnsembleProblem(prob, prob_func = prob_func2, safetycopy = false)
 sol = solve(monteprob, GPUTsit5(), EnsembleGPUKernel(CUDA.CUDABackend()),
     trajectories = 10_000)

src/DiffEqGPU.jl

@@ -37,6 +37,7 @@ include("ensemblegpuarray/kernels.jl")
 include("ensemblegpuarray/problem_generation.jl")
 include("ensemblegpuarray/lowerlevel_solve.jl")
 
+include("ensemblegpukernel/problems/ode_problems.jl")
 include("ensemblegpukernel/callbacks.jl")
 include("ensemblegpukernel/lowerlevel_solve.jl")
 include("ensemblegpukernel/gpukernel_algorithms.jl")
@@ -67,7 +68,7 @@ include("ensemblegpukernel/tableaus/verner_tableaus.jl")
 include("ensemblegpukernel/tableaus/rodas_tableaus.jl")
 include("ensemblegpukernel/tableaus/kvaerno_tableaus.jl")
 
-include("ensemblegpukernel/problems/ode_problems.jl")
+
 
 include("utils.jl")
 include("algorithms.jl")

src/ensemblegpukernel/integrators/nonstiff/types.jl

@@ -1,13 +1,24 @@
 ## Fixed TimeStep Integrator
 
-function Adapt.adapt_structure(to, prob::ODEProblem{<:Any, <:Any, iip}) where {iip}
-    ODEProblem{iip, true}(adapt(to, prob.f),
+function Adapt.adapt_structure(to, prob::Union{ODEProblem{<:Any, <:Any, iip}, ImmutableODEProblem{<:Any, <:Any, iip}}) where {iip}
+    ImmutableODEProblem{iip, true}(adapt(to, prob.f),
         adapt(to, prob.u0),
         adapt(to, prob.tspan),
         adapt(to, prob.p);
         adapt(to, prob.kwargs)...)
 end
 
+function Adapt.adapt_structure(to, f::ODEFunction{iip}) where {iip}
+    if f.mass_matrix !== I && f.initialization_data !== nothing
+        error("Adaptation to GPU failed: DAEs of ModelingToolkit currently not supported.")
+    end
+    ODEFunction{iip, SciMLBase.FullSpecialize}(f.f, jac = f.jac, mass_matrix = f.mass_matrix,
+    initializeprobmap = f.initializeprobmap, 
+    initializeprobpmap = f.initializeprobpmap, 
+    update_initializeprob! = f.update_initializeprob!, 
+    initialization_data = nothing, initializeprob = nothing)
+end
+
 mutable struct GPUTsit5Integrator{IIP, S, T, ST, P, F, TS, CB, AlgType} <:
                DiffEqBase.AbstractODEIntegrator{AlgType, IIP, S, T}
     alg::AlgType

src/solve.jl

@@ -273,7 +273,7 @@ function batch_solve_up_kernel(ensembleprob, probs, alg, ensemblealg, I, adaptiv
             _callback.continuous_callbacks)...)
 
     dev = ensemblealg.dev
-    probs = adapt(dev, probs)
+    probs = adapt(dev,adapt.((dev,), probs))
 
     #Adaptive version only works with saveat
     if adaptive