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("\nTesting 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("\nTesting 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("\nModelingToolkit DAE testing completed!")
0 commit comments