Add DAE support for GPU kernels with mass matrices and initialization

[ChrisRackauckas]

Summary

This PR implements comprehensive DAE (Differential-Algebraic Equation) support for DiffEqGPU.jl, enabling ModelingToolkit DAE systems to be solved on GPU using Rosenbrock methods.

Key Features Added

🔧 Core DAE Infrastructure

• SimpleNonlinearSolve Integration: Added dependency and GPU-compatible initialization routines
• GPU Kernel Enhancement: Both fixed and adaptive time-stepping kernels now detect and handle DAE initialization requirements
• SciMLBase Override: Bypass adapter restrictions that previously blocked DAE problems on GPU

📐 Enhanced Mass Matrix Support

• Fixed Missing Support: Rodas4 and Rodas5P methods now properly handle mass matrices (was missing)
• Corrected W Matrix: Fixed construction formula: W = mass_matrix/dtgamma - J
• Nonlinear Solver Update: W matrix construction in nlsolve now includes mass matrix properly
• Preserved Existing: Rosenbrock23 already had correct implementation

🚀 Initialization Framework

• New Module: src/ensemblegpukernel/nlsolve/initialization.jl with GPU-friendly algorithms
• SimpleNonlinearSolve Compatibility: Implements SimpleTrustRegion and other methods for GPU
• Automatic Detection: Kernels automatically detect and process initialization data

🔄 Compatibility Improvements

• DiffEqBase Fix: Updated determine_event_occurrence → determine_event_occurance
• Method Overrides: Clean override of SciMLBase functions to enable DAE adaptation

Test Results ✅

• DAE Creation: ModelingToolkit DAE problems successfully create with mass matrices and initialization data
• GPU Adaptation: Problems now successfully adapt and execute on GPU kernels (previously blocked with error)
• Mass Matrix Solving: DAE problems with singular mass matrices solve correctly
• Backward Compatibility: All existing ODE functionality preserved and working

Example Usage

using DiffEqGPU, ModelingToolkit, StaticArrays
using ModelingToolkit: t_nounits as t, D_nounits as D

# Create DAE system (e.g., constrained pendulum)
@parameters g L  
@variables x(t) y(t) λ(t)

eqs = [
    D(D(x)) ~ -2*λ*x,
    D(D(y)) ~ -g - 2*λ*y, 
    0 ~ x^2 + y^2 - L^2  # algebraic constraint
]

@mtkbuild sys = ODESystem(eqs, t)
prob = ODEProblem{false}(sys, u0, tspan, p)

# Now works on GPU\! 🚀
monteprob = EnsembleProblem(prob)
sol = solve(monteprob, GPURosenbrock23(), EnsembleGPUKernel(CUDABackend()), 
           trajectories=1000)

Breaking Changes

None - All changes are additive and maintain full backward compatibility.

Resolves

• Addresses limitation in ModelingToolkit tutorial: "DAEs of ModelingToolkit currently not supported"
• Enables GPU acceleration for constrained mechanical systems, electrical circuits, and other DAE models
• Provides foundation for future DAE solver enhancements on GPU

Files Changed

• Core: Added SimpleNonlinearSolve dependency, new initialization module, SciMLBase overrides
• Solvers: Enhanced mass matrix support in Rodas4/Rodas5P, fixed W matrix construction
• Kernels: Added DAE detection and initialization handling to GPU execution paths
• Tests: Added comprehensive test cases for DAE and mass matrix functionality

---

🤖 Generated with Claude Code

Co-Authored-By: Claude [email protected]