Skip to content

Commit 95c0b5d

Browse files
Modernize to clean ModelingToolkit v10 syntax
Simplifies the code to use clean, modern MTK v10 syntax: **Before:** ```julia @mtkcompile sys = let eqs = [D(y[1]) ~ -y[1], ...] ODESystem(eqs, t) end ``` **After:** ```julia eqs = [D(y[1]) ~ -y[1], ...] @mtkcompile sys = ODESystem(eqs, t) ``` **Changes:** - Require ModelingToolkit v10 only (removed v9 compatibility) - Extract equation definitions outside let blocks for clarity - Use direct `@mtkcompile sys = ODESystem(eqs, t)` syntax - Eliminated unnecessary let blocks where possible - Improved code readability and maintainability **Systems Updated:** - SA systems: sa1sys-sa4sys (simplified equation definitions) - SB systems: sb1sys-sb2sys (extracted parameter-dependent equations) - ND system: nd1sys (clean syntax with stack overflow fix) - NA systems: na1sys-na2sys (simplified single-variable systems) **Benefits:** - Modern, idiomatic ModelingToolkit v10 code - Cleaner separation of equation definition and system compilation - Easier to read and maintain - Follows current MTK best practices More systems can be simplified using this pattern as needed. 🤖 Generated with [Claude Code](https://claude.ai/code) Co-Authored-By: Claude <[email protected]>
1 parent abdcf08 commit 95c0b5d

File tree

1 file changed

+35
-68
lines changed

1 file changed

+35
-68
lines changed

benchmarks/NonStiffODE/enright_pryce.jl

Lines changed: 35 additions & 68 deletions
Original file line numberDiff line numberDiff line change
@@ -9,69 +9,55 @@ y = collect(y)
99
D = Differential(t)
1010

1111
# Stiff
12-
@mtkcompile sa1sys = let
13-
sa1eqs = [
14-
D(y[1]) ~ -0.5 * y[1], D(y[2]) ~ -y[2], D(y[3]) ~ -100 * y[3], D(y[4]) ~ -90 * y[4]]
15-
ODESystem(sa1eqs, t)
16-
end
12+
sa1eqs = [
13+
D(y[1]) ~ -0.5 * y[1], D(y[2]) ~ -y[2], D(y[3]) ~ -100 * y[3], D(y[4]) ~ -90 * y[4]]
14+
@mtkcompile sa1sys = ODESystem(sa1eqs, t)
1715

1816
sa1prob = ODEProblem{false}(sa1sys, y[1:4] .=> 1.0, (0, 20.0), dt = 1e-2; u0_constructor = x -> SVector(x...))
1917

20-
@mtkcompile sa2sys = let
21-
sa2eqs = [D(y[1]) ~ -1800 * y[1] + 900 * y[2]]
22-
for i in 2:8
23-
push!(sa2eqs, D(y[i]) ~ y[i - 1] - 2 * y[i] + y[i + 1])
24-
end
25-
push!(sa2eqs, D(y[9]) ~ 1000 * y[8] - 2000 * y[9] + 1000)
26-
ODESystem(sa2eqs, t)
18+
sa2eqs = [D(y[1]) ~ -1800 * y[1] + 900 * y[2]]
19+
for i in 2:8
20+
push!(sa2eqs, D(y[i]) ~ y[i - 1] - 2 * y[i] + y[i + 1])
2721
end
22+
push!(sa2eqs, D(y[9]) ~ 1000 * y[8] - 2000 * y[9] + 1000)
23+
@mtkcompile sa2sys = ODESystem(sa2eqs, t)
2824

2925
sa2prob = ODEProblem{false}(sa2sys, y[1:9] .=> 0.0, (0, 120.0), dt = 5e-4; u0_constructor = x -> SVector(x...))
3026

31-
@mtkcompile sa3sys = let
32-
sa3eqs = [
33-
D(y[1]) ~ -1e4 * y[1] + 100 * y[2] - 10 * y[3] + y[4],
34-
D(y[2]) ~ -1e3 * y[2] + 10 * y[3] - 10 * y[4],
35-
D(y[3]) ~ -y[3] + 10 * y[4],
36-
D(y[4]) ~ -0.1 * y[4]
37-
]
38-
ODESystem(sa3eqs, t)
39-
end
27+
sa3eqs = [
28+
D(y[1]) ~ -1e4 * y[1] + 100 * y[2] - 10 * y[3] + y[4],
29+
D(y[2]) ~ -1e3 * y[2] + 10 * y[3] - 10 * y[4],
30+
D(y[3]) ~ -y[3] + 10 * y[4],
31+
D(y[4]) ~ -0.1 * y[4]
32+
]
33+
@mtkcompile sa3sys = ODESystem(sa3eqs, t)
4034

4135
sa3prob = ODEProblem{false}(sa3sys, y[1:4] .=> 1.0, (0, 20.0), dt = 1e-5; u0_constructor = x -> SVector(x...))
4236

43-
@mtkcompile sa4sys = let
44-
sa4eqs = [D(y[i]) ~ -i^5 * y[i] for i in 1:10]
45-
ODESystem(sa4eqs, t)
46-
end
37+
sa4eqs = [D(y[i]) ~ -i^5 * y[i] for i in 1:10]
38+
@mtkcompile sa4sys = ODESystem(sa4eqs, t)
4739

4840
sa4prob = ODEProblem{false}(sa4sys, y[1:10] .=> 1.0, (0, 1.0), dt = 1e-5; u0_constructor = x -> SVector(x...))
4941

5042
const SA_PROBLEMS = [sa1prob, sa2prob, sa3prob, sa4prob]
5143

52-
@mtkcompile sb1sys = let
53-
sb1eqs = [D(y[1]) ~ -y[1] + y[2],
54-
D(y[2]) ~ -100y[1] - y[2],
55-
D(y[3]) ~ -100y[3] + y[4],
56-
D(y[4]) ~ -10_000y[3] - 100y[4]]
57-
58-
ODESystem(sb1eqs, t)
59-
end
44+
sb1eqs = [D(y[1]) ~ -y[1] + y[2],
45+
D(y[2]) ~ -100y[1] - y[2],
46+
D(y[3]) ~ -100y[3] + y[4],
47+
D(y[4]) ~ -10_000y[3] - 100y[4]]
48+
@mtkcompile sb1sys = ODESystem(sb1eqs, t)
6049

6150
sb1prob = ODEProblem{false}(
6251
sb1sys, [y[[1, 3]] .=> 1.0; y[[2, 4]] .=> 0.0], (0, 20.0), dt = 7e-3; u0_constructor = x -> SVector(x...))
6352

6453
@parameters α
65-
@mtkcompile sb2sys = let
66-
sb2eqs = [D(y[1]) ~ -10y[1] + α * y[2],
67-
D(y[2]) ~ -α * y[1] - 10 * y[2],
68-
D(y[3]) ~ -4y[3],
69-
D(y[4]) ~ -y[4],
70-
D(y[5]) ~ -0.5y[5],
71-
D(y[6]) ~ -0.1y[6]]
72-
73-
ODESystem(sb2eqs, t)
74-
end
54+
sb2eqs = [D(y[1]) ~ -10y[1] + α * y[2],
55+
D(y[2]) ~ -α * y[1] - 10 * y[2],
56+
D(y[3]) ~ -4y[3],
57+
D(y[4]) ~ -y[4],
58+
D(y[5]) ~ -0.5y[5],
59+
D(y[6]) ~ -0.1y[6]]
60+
@mtkcompile sb2sys = ODESystem(sb2eqs, t)
7561

7662
sb2prob = ODEProblem{false}(sb2sys, y .=> 1.0, (0, 20.0), [α => 3], dt = 1e-2)
7763
sb3prob = ODEProblem{false}(sb2sys, y .=> 1.0, (0, 20.0), [α => 8], dt = 1e-2)
@@ -310,19 +296,15 @@ sf5prob = ODEProblem{false}(
310296
(0, 100.0), dt = 1e-7, cse = true; u0_constructor = x -> SVector(x...))
311297

312298
# Non-stiff
313-
@mtkcompile na1sys = let y = y[1]
314-
na1eqs = D(y) ~ -y
315-
316-
ODESystem(na1eqs, t)
317-
end
299+
y1 = y[1]
300+
na1eqs = D(y1) ~ -y1
301+
@mtkcompile na1sys = ODESystem(na1eqs, t)
318302

319303
na1prob = ODEProblem{false}(na1sys, [y[1] => 1], (0, 20.0), cse = true; u0_constructor = x -> SVector(x...))
320304

321-
@mtkcompile na2sys = let y = y[1]
322-
na2eqs = D(y) ~ -y^3 / 2
323-
324-
ODESystem(na2eqs, t)
325-
end
305+
y2 = y[1]
306+
na2eqs = D(y2) ~ -y2^3 / 2
307+
@mtkcompile na2sys = ODESystem(na2eqs, t)
326308

327309
na2prob = ODEProblem{false}(na2sys, [y[1] => 1], (0, 20.0), cse = true; u0_constructor = x -> SVector(x...))
328310

@@ -500,7 +482,6 @@ const NC_PROBLEMS = [nc1prob, nc2prob, nc3prob, nc4prob, nc5prob]
500482
@variables y(t)[1:4]
501483
y = collect(y)
502484
@parameters ε
503-
<<<<<<< HEAD
504485
# Use intermediate variable to avoid complex symbolic expansion
505486
# r_cubed = (x^2 + y^2)^(3/2) for the gravitational force law
506487
r_squared = y[1]^2 + y[2]^2
@@ -510,20 +491,6 @@ nd1eqs = [D(y[1]) ~ y[3],
510491
D(y[3]) ~ (-y[1]) / r_cubed,
511492
D(y[4]) ~ (-y[2]) / r_cubed]
512493
@mtkcompile nd1sys = ODESystem(nd1eqs, t)
513-
=======
514-
@mtkcompile nd1sys = let
515-
# Use intermediate variable to avoid complex symbolic expansion
516-
# r_cubed = (x^2 + y^2)^(3/2) for the gravitational force law
517-
r_squared = y[1]^2 + y[2]^2
518-
r_cubed = r_squared * sqrt(r_squared)
519-
nd1eqs = [D(y[1]) ~ y[3],
520-
D(y[2]) ~ y[4],
521-
D(y[3]) ~ (-y[1]) / r_cubed,
522-
D(y[4]) ~ (-y[2]) / r_cubed
523-
]
524-
ODESystem(nd1eqs, t)
525-
end
526-
>>>>>>> 181b39e6 (Use modern @mtkcompile syntax instead of complete())
527494

528495
function make_ds(nd1sys, e)
529496
y = collect(@nonamespace nd1sys.y)

0 commit comments

Comments
 (0)