1+ using ModelingToolkit, DiffEqGPU, OrdinaryDiffEq
2+ using KernelAbstractions: CPU
3+ using ModelingToolkit: t_nounits as t, D_nounits as D
4+
5+ println(" Testing ModelingToolkit DAE support with Cartesian Pendulum..." )
6+
7+ # Define the cartesian pendulum DAE system
8+ @parameters g = 9.81 L = 1.0
9+ @variables x(t) y(t) [state_priority = 10 ] λ(t)
10+
11+ # The cartesian pendulum DAE system:
12+ # m*ddot(x) = (x/L)*λ (simplified with m=1)
13+ # m*ddot(y) = (y/L)*λ - mg (simplified with m=1)
14+ # x^2 + y^2 = L^2 (constraint equation)
15+ eqs = [D(D(x)) ~ λ * x / L
16+ D(D(y)) ~ λ * y / L - g
17+ x^ 2 + y^ 2 ~ L^ 2 ]
18+
19+ @named pendulum = ODESystem(eqs, t, [x, y, λ], [g, L])
20+
21+ # Perform structural simplification with index reduction
22+ pendulum_sys = structural_simplify(dae_index_lowering(pendulum))
23+
24+ # Initial conditions: pendulum starts at bottom right position
25+ u0 = [x => 1.0 , y => 0.0 , λ => - g] # λ initial guess for tension
26+
27+ # Time span
28+ tspan = (0.0f0 , 1.0f0 )
29+
30+ # Create the ODE problem
31+ prob = ODEProblem(pendulum_sys, u0, tspan, Float32[])
32+
33+ println(" ModelingToolkit DAE problem created successfully!" )
34+ println(" Number of equations: " , length(equations(pendulum_sys)))
35+ println(" Variables: " , unknowns(pendulum_sys))
36+
37+ # Test if the problem has initialization data
38+ if SciMLBase. has_initialization_data(prob. f)
39+ println(" Problem has initialization data: ✓" )
40+ else
41+ println(" Problem has initialization data: ✗" )
42+ end
43+
44+ # Test mass matrix presence
45+ if prob. f. mass_matrix != = nothing
46+ println(" Problem has mass matrix: ✓" )
47+ println(" Mass matrix size: " , size(prob. f. mass_matrix))
48+ else
49+ println(" Problem has mass matrix: ✗" )
50+ end
51+
52+ # Create ensemble problem for GPU testing
53+ monteprob = EnsembleProblem(prob, safetycopy = false )
54+
55+ # Test with CPU backend first
56+ println(" \n Testing with GPURosenbrock23 on CPU backend..." )
57+ try
58+ sol = solve(monteprob, GPURosenbrock23(), EnsembleGPUKernel(CPU()),
59+ trajectories = 4 ,
60+ dt = 0.01f0 ,
61+ adaptive = false )
62+
63+ println(" ✓ ModelingToolkit DAE GPU solution successful!" )
64+ println(" Number of solutions: " , length(sol. u))
65+
66+ if length(sol. u) > 0
67+ println(" Final state of first trajectory: " , sol. u[1 ][end ])
68+
69+ # Check constraint satisfaction: x^2 + y^2 ≈ L^2
70+ final_state = sol. u[1 ][end ]
71+ x_final, y_final = final_state[1 ], final_state[2 ]
72+ constraint_error = abs(x_final^ 2 + y_final^ 2 - 1.0f0 )
73+ println(" Constraint error |x² + y² - L²|: " , constraint_error)
74+
75+ if constraint_error < 0.1f0 # Reasonable tolerance for fixed timestep
76+ println(" ✓ Constraint satisfied within tolerance" )
77+ else
78+ println(" ⚠ Constraint violation detected" )
79+ end
80+ end
81+
82+ catch e
83+ println(" ✗ ModelingToolkit DAE GPU solution failed: " , e)
84+ println(" Detailed error: " )
85+ println(sprint(showerror, e, catch_backtrace()))
86+ end
87+
88+ # Test with Rodas4 as well to show mass matrix support
89+ println(" \n Testing with GPURodas4 on CPU backend..." )
90+ try
91+ sol = solve(monteprob, GPURodas4(), EnsembleGPUKernel(CPU()),
92+ trajectories = 4 ,
93+ dt = 0.01f0 ,
94+ adaptive = false )
95+
96+ println(" ✓ ModelingToolkit DAE with GPURodas4 successful!" )
97+ println(" Number of solutions: " , length(sol. u))
98+
99+ catch e
100+ println(" ✗ ModelingToolkit DAE with GPURodas4 failed: " , e)
101+ println(" Error details: " , sprint(showerror, e, catch_backtrace()))
102+ end
103+
104+ println(" \n ModelingToolkit DAE testing completed!" )
0 commit comments