@@ -13,8 +13,8 @@ SUITE["nonstiff"] = BenchmarkGroup()
13
13
function lotka_volterra! (du, u, p, t)
14
14
α, β, γ, δ = p
15
15
x, y = u
16
- du[1 ] = α* x - β* x * y # dx/dt
17
- du[2 ] = δ* x * y - γ* y # dy/dt
16
+ du[1 ] = α * x - β * x * y # dx/dt
17
+ du[2 ] = δ * x * y - γ * y # dy/dt
18
18
nothing
19
19
end
20
20
@@ -33,24 +33,24 @@ function pleiades!(du, u, p, t)
33
33
y = @view u[8 : 14 ]
34
34
vx = @view u[15 : 21 ]
35
35
vy = @view u[22 : 28 ]
36
-
36
+
37
37
# Copy velocities to position derivatives
38
38
du[1 : 7 ] .= vx
39
39
du[8 : 14 ] .= vy
40
-
40
+
41
41
# Calculate accelerations
42
42
fill! (du[15 : 21 ], 0.0 )
43
43
fill! (du[22 : 28 ], 0.0 )
44
-
44
+
45
45
for i in 1 : 7
46
46
for j in 1 : 7
47
47
if i != j
48
48
dx = x[j] - x[i]
49
49
dy = y[j] - y[i]
50
50
r = sqrt (dx^ 2 + dy^ 2 )
51
51
r3 = r^ 3
52
- du[14 + i] += j * dx / r3 # mass j = j
53
- du[21 + i] += j * dy / r3
52
+ du[14 + i] += j * dx / r3 # mass j = j
53
+ du[21 + i] += j * dy / r3
54
54
end
55
55
end
56
56
end
59
59
60
60
function pleiades_prob ()
61
61
# Initial conditions from literature
62
- u0 = [3.0 ,3.0 ,- 1.0 ,- 3.0 ,2.0 ,- 2.0 ,2.0 ,3.0 ,- 3.0 ,2.0 ,0 , 0 , - 4.0 ,4.0 ,
63
- 0 , 0 , 0 , 0 , 0 , 1.75 ,- 1.5 ,0 , 0 , 0 , - 1.25 ,1 , 0 , 0 ]
62
+ u0 = [3.0 , 3.0 , - 1.0 , - 3.0 , 2.0 , - 2.0 , 2.0 , 3.0 , - 3.0 , 2.0 , 0 , 0 , - 4.0 , 4.0 ,
63
+ 0 , 0 , 0 , 0 , 0 , 1.75 , - 1.5 , 0 , 0 , 0 , - 1.25 , 1 , 0 , 0 ]
64
64
tspan = (0.0 , 3.0 )
65
65
ODEProblem (pleiades!, u0, tspan)
66
66
end
69
69
function fitzhugh_nagumo! (du, u, p, t)
70
70
a, b, c = p
71
71
v, w = u
72
- du[1 ] = c * (v - v^ 3 / 3 + w) # dv/dt
73
- du[2 ] = - (v - a + b* w) / c # dw/dt
72
+ du[1 ] = c * (v - v^ 3 / 3 + w) # dv/dt
73
+ du[2 ] = - (v - a + b * w) / c # dw/dt
74
74
nothing
75
75
end
76
76
@@ -91,9 +91,9 @@ SUITE["stiff"] = BenchmarkGroup()
91
91
function rober! (du, u, p, t)
92
92
k1, k2, k3 = p
93
93
y1, y2, y3 = u
94
- du[1 ] = - k1* y1 + k3* y2 * y3 # dy1/dt
95
- du[2 ] = k1* y1 - k2* y2^ 2 - k3* y2 * y3 # dy2/dt
96
- du[3 ] = k2* y2^ 2 # dy3/dt
94
+ du[1 ] = - k1 * y1 + k3 * y2 * y3 # dy1/dt
95
+ du[2 ] = k1 * y1 - k2 * y2^ 2 - k3 * y2 * y3 # dy2/dt
96
+ du[3 ] = k2 * y2^ 2 # dy3/dt
97
97
nothing
98
98
end
99
99
@@ -109,7 +109,7 @@ function van_der_pol!(du, u, p, t)
109
109
μ = p[1 ]
110
110
x, y = u
111
111
du[1 ] = y # dx/dt
112
- du[2 ] = μ * ((1 - x^ 2 )* y - x) # dy/dt
112
+ du[2 ] = μ * ((1 - x^ 2 ) * y - x) # dy/dt
113
113
nothing
114
114
end
115
115
@@ -121,96 +121,95 @@ function van_der_pol_prob()
121
121
end
122
122
123
123
# Pollution Problem (atmospheric chemistry, 20D stiff system)
124
- const k1= .35e0
125
- const k2= .266e2
126
- const k3= .123e5
127
- const k4= .86e-3
128
- const k5= .82e-3
129
- const k6= .15e5
130
- const k7= .13e-3
131
- const k8= .24e5
132
- const k9= .165e5
133
- const k10= .9e4
134
- const k11= .22e-1
135
- const k12= .12e5
136
- const k13= .188e1
137
- const k14= .163e5
138
- const k15= .48e7
139
- const k16= .35e-3
140
- const k17= .175e-1
141
- const k18= .1e9
142
- const k19= .444e12
143
- const k20= .124e4
144
- const k21= .21e1
145
- const k22= .578e1
146
- const k23= .474e-1
147
- const k24= .178e4
148
- const k25= .312e1
124
+ const k1 = .35e0
125
+ const k2 = .266e2
126
+ const k3 = .123e5
127
+ const k4 = .86e-3
128
+ const k5 = .82e-3
129
+ const k6 = .15e5
130
+ const k7 = .13e-3
131
+ const k8 = .24e5
132
+ const k9 = .165e5
133
+ const k10 = .9e4
134
+ const k11 = .22e-1
135
+ const k12 = .12e5
136
+ const k13 = .188e1
137
+ const k14 = .163e5
138
+ const k15 = .48e7
139
+ const k16 = .35e-3
140
+ const k17 = .175e-1
141
+ const k18 = .1e9
142
+ const k19 = .444e12
143
+ const k20 = .124e4
144
+ const k21 = .21e1
145
+ const k22 = .578e1
146
+ const k23 = .474e-1
147
+ const k24 = .178e4
148
+ const k25 = .312e1
149
149
150
150
function pollution! (dy, y, p, t)
151
- r1 = k1 * y[1 ]
152
- r2 = k2 * y[2 ]* y[4 ]
153
- r3 = k3 * y[5 ]* y[2 ]
154
- r4 = k4 * y[7 ]
155
- r5 = k5 * y[7 ]
156
- r6 = k6 * y[7 ]* y[6 ]
157
- r7 = k7 * y[9 ]
158
- r8 = k8 * y[9 ]* y[6 ]
159
- r9 = k9 * y[11 ]* y[2 ]
160
- r10 = k10* y[11 ]* y[1 ]
161
- r11 = k11* y[13 ]
162
- r12 = k12* y[10 ]* y[2 ]
163
- r13 = k13* y[14 ]
164
- r14 = k14* y[1 ]* y[6 ]
165
- r15 = k15* y[3 ]
166
- r16 = k16* y[4 ]
167
- r17 = k17* y[4 ]
168
- r18 = k18* y[16 ]
169
- r19 = k19* y[16 ]
170
- r20 = k20* y[17 ]* y[6 ]
171
- r21 = k21* y[19 ]
172
- r22 = k22* y[19 ]
173
- r23 = k23* y[1 ]* y[4 ]
174
- r24 = k24* y[19 ]* y[1 ]
175
- r25 = k25* y[20 ]
176
-
177
- dy[1 ] = - r1- r10- r14- r23- r24+
178
- r2 + r3 + r9 + r11+ r12+ r22+ r25
179
- dy[2 ] = - r2- r3 - r9 - r12+ r1 + r21
180
- dy[3 ] = - r15+ r1 + r17+ r19+ r22
181
- dy[4 ] = - r2- r16- r17- r23+ r15
182
- dy[5 ] = - r3+ r4 + r4 + r6 + r7 + r13+ r20
183
- dy[6 ] = - r6- r8 - r14- r20+ r3 + r18+ r18
184
- dy[7 ] = - r4- r5 - r6 + r13
185
- dy[8 ] = r4+ r5 + r6 + r7
186
- dy[9 ] = - r7- r8
187
- dy[10 ] = - r12+ r7 + r9
188
- dy[11 ] = - r9- r10+ r8 + r11
151
+ r1 = k1 * y[1 ]
152
+ r2 = k2 * y[2 ] * y[4 ]
153
+ r3 = k3 * y[5 ] * y[2 ]
154
+ r4 = k4 * y[7 ]
155
+ r5 = k5 * y[7 ]
156
+ r6 = k6 * y[7 ] * y[6 ]
157
+ r7 = k7 * y[9 ]
158
+ r8 = k8 * y[9 ] * y[6 ]
159
+ r9 = k9 * y[11 ] * y[2 ]
160
+ r10 = k10 * y[11 ] * y[1 ]
161
+ r11 = k11 * y[13 ]
162
+ r12 = k12 * y[10 ] * y[2 ]
163
+ r13 = k13 * y[14 ]
164
+ r14 = k14 * y[1 ] * y[6 ]
165
+ r15 = k15 * y[3 ]
166
+ r16 = k16 * y[4 ]
167
+ r17 = k17 * y[4 ]
168
+ r18 = k18 * y[16 ]
169
+ r19 = k19 * y[16 ]
170
+ r20 = k20 * y[17 ] * y[6 ]
171
+ r21 = k21 * y[19 ]
172
+ r22 = k22 * y[19 ]
173
+ r23 = k23 * y[1 ] * y[4 ]
174
+ r24 = k24 * y[19 ] * y[1 ]
175
+ r25 = k25 * y[20 ]
176
+
177
+ dy[1 ] = - r1 - r10 - r14 - r23 - r24 +
178
+ r2 + r3 + r9 + r11 + r12 + r22 + r25
179
+ dy[2 ] = - r2 - r3 - r9 - r12 + r1 + r21
180
+ dy[3 ] = - r15 + r1 + r17 + r19 + r22
181
+ dy[4 ] = - r2 - r16 - r17 - r23 + r15
182
+ dy[5 ] = - r3 + r4 + r4 + r6 + r7 + r13 + r20
183
+ dy[6 ] = - r6 - r8 - r14 - r20 + r3 + r18 + r18
184
+ dy[7 ] = - r4 - r5 - r6 + r13
185
+ dy[8 ] = r4 + r5 + r6 + r7
186
+ dy[9 ] = - r7 - r8
187
+ dy[10 ] = - r12 + r7 + r9
188
+ dy[11 ] = - r9 - r10 + r8 + r11
189
189
dy[12 ] = r9
190
- dy[13 ] = - r11+ r10
191
- dy[14 ] = - r13+ r12
190
+ dy[13 ] = - r11 + r10
191
+ dy[14 ] = - r13 + r12
192
192
dy[15 ] = r14
193
- dy[16 ] = - r18- r19+ r16
193
+ dy[16 ] = - r18 - r19 + r16
194
194
dy[17 ] = - r20
195
195
dy[18 ] = r20
196
- dy[19 ] = - r21- r22- r24+ r23+ r25
197
- dy[20 ] = - r25+ r24
196
+ dy[19 ] = - r21 - r22 - r24 + r23 + r25
197
+ dy[20 ] = - r25 + r24
198
198
nothing
199
199
end
200
200
201
201
function pollution_prob ()
202
202
u0 = zeros (20 )
203
- u0[2 ] = 0.2
204
- u0[4 ] = 0.04
205
- u0[7 ] = 0.1
206
- u0[8 ] = 0.3
207
- u0[9 ] = 0.01
203
+ u0[2 ] = 0.2
204
+ u0[4 ] = 0.04
205
+ u0[7 ] = 0.1
206
+ u0[8 ] = 0.3
207
+ u0[9 ] = 0.01
208
208
u0[17 ] = 0.007
209
209
tspan = (0.0 , 60.0 )
210
210
ODEProblem (pollution!, u0, tspan)
211
211
end
212
212
213
-
214
213
# =============================================================================
215
214
# Scaling Problems (different dimensions)
216
215
# =============================================================================
@@ -243,22 +242,25 @@ limit(a, N) = a == N + 1 ? 1 : a == 0 ? N : a
243
242
function brusselator_2d! (du, u, p, t)
244
243
A, B, α, dx, N = p
245
244
α = α / dx^ 2
246
-
245
+
247
246
# Create coordinate arrays for this N
248
247
xyd = range (0 , stop = 1 , length = N)
249
-
248
+
250
249
@inbounds for I in CartesianIndices ((N, N))
251
250
i, j = Tuple (I)
252
251
x, y = xyd[i], xyd[j]
253
- ip1, im1, jp1, jm1 = limit (i + 1 , N), limit (i - 1 , N), limit (j + 1 , N), limit (j - 1 , N)
254
-
252
+ ip1, im1, jp1, jm1 = limit (i + 1 , N), limit (i - 1 , N), limit (j + 1 , N),
253
+ limit (j - 1 , N)
254
+
255
255
# u equation: ∂u/∂t = 1 + u²v - 4.4u + α∇²u + f(x,y,t)
256
- du[i, j, 1 ] = α * (u[im1, j, 1 ] + u[ip1, j, 1 ] + u[i, jp1, 1 ] + u[i, jm1, 1 ] - 4 * u[i, j, 1 ]) +
257
- B + u[i, j, 1 ]^ 2 * u[i, j, 2 ] - (A + 1 ) * u[i, j, 1 ] +
256
+ du[i, j, 1 ] = α * (u[im1, j, 1 ] + u[ip1, j, 1 ] + u[i, jp1, 1 ] + u[i, jm1, 1 ] -
257
+ 4 * u[i, j, 1 ]) +
258
+ B + u[i, j, 1 ]^ 2 * u[i, j, 2 ] - (A + 1 ) * u[i, j, 1 ] +
258
259
brusselator_f (x, y, t)
259
-
260
+
260
261
# v equation: ∂v/∂t = 3.4u - u²v + α∇²v
261
- du[i, j, 2 ] = α * (u[im1, j, 2 ] + u[ip1, j, 2 ] + u[i, jp1, 2 ] + u[i, jm1, 2 ] - 4 * u[i, j, 2 ]) +
262
+ du[i, j, 2 ] = α * (u[im1, j, 2 ] + u[ip1, j, 2 ] + u[i, jp1, 2 ] + u[i, jm1, 2 ] -
263
+ 4 * u[i, j, 2 ]) +
262
264
A * u[i, j, 1 ] - u[i, j, 1 ]^ 2 * u[i, j, 2 ]
263
265
end
264
266
nothing
@@ -270,8 +272,8 @@ function init_brusselator_2d(N::Int)
270
272
for I in CartesianIndices ((N, N))
271
273
x = xyd[I[1 ]]
272
274
y = xyd[I[2 ]]
273
- u[I, 1 ] = 22 * (y * (1 - y))^ (3 / 2 ) # u initial condition
274
- u[I, 2 ] = 27 * (x * (1 - x))^ (3 / 2 ) # v initial condition
275
+ u[I, 1 ] = 22 * (y * (1 - y))^ (3 / 2 ) # u initial condition
276
+ u[I, 2 ] = 27 * (x * (1 - x))^ (3 / 2 ) # v initial condition
275
277
end
276
278
u
277
279
end
@@ -298,14 +300,17 @@ fn_prob = fitzhugh_nagumo_prob()
298
300
explicit_solvers = [Tsit5 (), Vern6 (), Vern7 (), DP5 (), BS3 ()]
299
301
300
302
SUITE[" nonstiff" ][" lotka_volterra" ] = BenchmarkGroup ()
301
- SUITE[" nonstiff" ][" pleiades" ] = BenchmarkGroup ()
303
+ SUITE[" nonstiff" ][" pleiades" ] = BenchmarkGroup ()
302
304
SUITE[" nonstiff" ][" fitzhugh_nagumo" ] = BenchmarkGroup ()
303
305
304
306
for solver in explicit_solvers
305
307
solver_name = string (typeof (solver). name. name)
306
- SUITE[" nonstiff" ][" lotka_volterra" ][solver_name] = @benchmarkable solve ($ lv_prob, $ solver, reltol= 1e-6 , abstol= 1e-8 )
307
- SUITE[" nonstiff" ][" pleiades" ][solver_name] = @benchmarkable solve ($ pl_prob, $ solver, reltol= 1e-6 , abstol= 1e-8 )
308
- SUITE[" nonstiff" ][" fitzhugh_nagumo" ][solver_name] = @benchmarkable solve ($ fn_prob, $ solver, reltol= 1e-6 , abstol= 1e-8 )
308
+ SUITE[" nonstiff" ][" lotka_volterra" ][solver_name] = @benchmarkable solve (
309
+ $ lv_prob, $ solver, reltol = 1e-6 , abstol = 1e-8 )
310
+ SUITE[" nonstiff" ][" pleiades" ][solver_name] = @benchmarkable solve (
311
+ $ pl_prob, $ solver, reltol = 1e-6 , abstol = 1e-8 )
312
+ SUITE[" nonstiff" ][" fitzhugh_nagumo" ][solver_name] = @benchmarkable solve (
313
+ $ fn_prob, $ solver, reltol = 1e-6 , abstol = 1e-8 )
309
314
end
310
315
311
316
# Stiff benchmarks with different solvers
@@ -322,9 +327,12 @@ SUITE["stiff"]["pollution"] = BenchmarkGroup()
322
327
323
328
for solver in stiff_solvers
324
329
solver_name = string (typeof (solver). name. name)
325
- SUITE[" stiff" ][" rober" ][solver_name] = @benchmarkable solve ($ rober_prob_instance, $ solver, reltol= 1e-6 , abstol= 1e-8 )
326
- SUITE[" stiff" ][" van_der_pol" ][solver_name] = @benchmarkable solve ($ vdp_prob, $ solver, reltol= 1e-6 , abstol= 1e-8 )
327
- SUITE[" stiff" ][" pollution" ][solver_name] = @benchmarkable solve ($ pollution_prob_instance, $ solver, reltol= 1e-6 , abstol= 1e-8 )
330
+ SUITE[" stiff" ][" rober" ][solver_name] = @benchmarkable solve (
331
+ $ rober_prob_instance, $ solver, reltol = 1e-6 , abstol = 1e-8 )
332
+ SUITE[" stiff" ][" van_der_pol" ][solver_name] = @benchmarkable solve (
333
+ $ vdp_prob, $ solver, reltol = 1e-6 , abstol = 1e-8 )
334
+ SUITE[" stiff" ][" pollution" ][solver_name] = @benchmarkable solve (
335
+ $ pollution_prob_instance, $ solver, reltol = 1e-6 , abstol = 1e-8 )
328
336
end
329
337
330
338
# Scaling benchmarks
@@ -334,14 +342,15 @@ SUITE["scaling"]["brusselator_2d"] = BenchmarkGroup()
334
342
# Linear ODE scaling (different problem sizes)
335
343
for N in [10 , 50 , 100 ]
336
344
prob = create_linear_prob (N)
337
- SUITE[" scaling" ][" linear" ][" N$N " ] = @benchmarkable solve ($ prob, Tsit5 (), reltol= 1e-6 , abstol= 1e-8 )
345
+ SUITE[" scaling" ][" linear" ][" N$N " ] = @benchmarkable solve (
346
+ $ prob, Tsit5 (), reltol = 1e-6 , abstol = 1e-8 )
338
347
end
339
348
340
349
# Brusselator 2D scaling (different grid sizes)
341
350
for N in [8 , 16 , 32 ]
342
351
prob = create_brusselator_2d_prob (N)
343
- SUITE[" scaling" ][" brusselator_2d" ][" $(N) x$(N) " ] = @benchmarkable solve ($ prob, TRBDF2 (),
344
- reltol= 1e-4 , abstol= 1e-6 , maxiters= 1000 )
352
+ SUITE[" scaling" ][" brusselator_2d" ][" $(N) x$(N) " ] = @benchmarkable solve ($ prob, TRBDF2 (),
353
+ reltol = 1e-4 , abstol = 1e-6 , maxiters = 1000 )
345
354
end
346
355
347
356
# =============================================================================
@@ -353,4 +362,4 @@ SUITE["construction"] = BenchmarkGroup()
353
362
# Test problem construction overhead
354
363
SUITE[" construction" ][" lotka_volterra" ] = @benchmarkable lotka_volterra_prob ()
355
364
SUITE[" construction" ][" rober" ] = @benchmarkable rober_prob ()
356
- SUITE[" construction" ][" linear_N50" ] = @benchmarkable create_linear_prob (50 )
365
+ SUITE[" construction" ][" linear_N50" ] = @benchmarkable create_linear_prob (50 )
0 commit comments