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