Skip to content

Commit 364d7ae

Browse files
oscardssmithoscardssmith
committed
fixes
1 parent 794c3fa commit 364d7ae

File tree

1 file changed

+24
-38
lines changed

1 file changed

+24
-38
lines changed

‎lib/OrdinaryDiffEqExtrapolation/src/extrapolation_perform_step.jl

Lines changed: 24 additions & 38 deletions
Original file line numberDiff line numberDiff line change
@@ -1173,7 +1173,7 @@ function perform_step!(integrator, cache::ImplicitDeuflhardExtrapolationCache,
11731173
jacobian2W!(W[1], integrator.f.mass_matrix, dt_int, J, true)
11741174
integrator.stats.nw += 1
11751175
@.. broadcast=false u_temp2=uprev
1176-
@.. broadcast=false linsolve_tmps[1]=dt_int * fsalfirst
1176+
@.. broadcast=false linsolve_tmps[1]=fsalfirst
11771177

11781178
linsolve = cache.linsolve[1]
11791179

@@ -1188,13 +1188,12 @@ function perform_step!(integrator, cache::ImplicitDeuflhardExtrapolationCache,
11881188
cache.linsolve[1] = linres.cache
11891189

11901190
integrator.stats.nsolve += 1
1191-
@.. broadcast=false k=-k
1192-
@.. broadcast=false u_temp1=u_temp2 + k # Euler starting step
1191+
@.. broadcast=false u_temp1=u_temp2 - k # Euler starting step
11931192
@.. broadcast=false diff1[1]=u_temp1 - u_temp2
11941193
for j in 2:j_int
11951194
f(k, cache.u_temp1, p, t + (j - 1) * dt_int)
11961195
OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)
1197-
@.. broadcast=false linsolve_tmps[1]=dt_int * k - (u_temp1 - u_temp2)
1196+
@.. broadcast=false linsolve_tmps[1]=k - (u_temp1 - u_temp2)/dt_int
11981197

11991198
linsolve = cache.linsolve[1]
12001199

@@ -1208,8 +1207,7 @@ function perform_step!(integrator, cache::ImplicitDeuflhardExtrapolationCache,
12081207
cache.linsolve[1] = linres.cache
12091208

12101209
integrator.stats.nsolve += 1
1211-
@.. broadcast=false k=-k
1212-
@.. broadcast=false T[i + 1]=2 * u_temp1 - u_temp2 + 2 * k # Explicit Midpoint rule
1210+
@.. broadcast=false T[i + 1]=2 * u_temp1 - u_temp2 - 2 * k # Explicit Midpoint rule
12131211
@.. broadcast=false u_temp2=u_temp1
12141212
@.. broadcast=false u_temp1=T[i + 1]
12151213
if (i <= 1)
@@ -1245,8 +1243,7 @@ function perform_step!(integrator, cache::ImplicitDeuflhardExtrapolationCache,
12451243
jacobian2W!(W[Threads.threadid()], integrator.f.mass_matrix,
12461244
dt_int_temp, J, true)
12471245
@.. broadcast=false u_temp4[Threads.threadid()]=uprev
1248-
@.. broadcast=false linsolve_tmps[Threads.threadid()]=dt_int_temp *
1249-
fsalfirst
1246+
@.. broadcast=false linsolve_tmps[Threads.threadid()]=fsalfirst
12501247

12511248
linsolve = cache.linsolve[Threads.threadid()]
12521249

@@ -1271,10 +1268,9 @@ function perform_step!(integrator, cache::ImplicitDeuflhardExtrapolationCache,
12711268
f(k_tmps[Threads.threadid()],
12721269
cache.u_temp3[Threads.threadid()],
12731270
p, t + (j - 1) * dt_int_temp)
1274-
@.. broadcast=false linsolve_tmps[Threads.threadid()]=dt_int_temp *
1275-
k_tmps[Threads.threadid()] -
1271+
@.. broadcast=false linsolve_tmps[Threads.threadid()]=k_tmps[Threads.threadid()] -
12761272
(u_temp3[Threads.threadid()] -
1277-
u_temp4[Threads.threadid()])
1273+
u_temp4[Threads.threadid()])/dt_int_temp
12781274

12791275
linsolve = cache.linsolve[Threads.threadid()]
12801276

@@ -1290,10 +1286,9 @@ function perform_step!(integrator, cache::ImplicitDeuflhardExtrapolationCache,
12901286
end
12911287
cache.linsolve[Threads.threadid()] = linres.cache
12921288

1293-
@.. broadcast=false k_tmps[Threads.threadid()]=-k_tmps[Threads.threadid()]
12941289
@.. broadcast=false T[index + 1]=2 *
12951290
u_temp3[Threads.threadid()] -
1296-
u_temp4[Threads.threadid()] +
1291+
u_temp4[Threads.threadid()] -
12971292
2 * k_tmps[Threads.threadid()] # Explicit Midpoint rule
12981293
@.. broadcast=false u_temp4[Threads.threadid()]=u_temp3[Threads.threadid()]
12991294
@.. broadcast=false u_temp3[Threads.threadid()]=T[index + 1]
@@ -1332,8 +1327,7 @@ function perform_step!(integrator, cache::ImplicitDeuflhardExtrapolationCache,
13321327
jacobian2W!(W[Threads.threadid()], integrator.f.mass_matrix,
13331328
dt_int_temp, J, true)
13341329
@.. broadcast=false u_temp4[Threads.threadid()]=uprev
1335-
@.. broadcast=false linsolve_tmps[Threads.threadid()]=dt_int_temp *
1336-
fsalfirst
1330+
@.. broadcast=false linsolve_tmps[Threads.threadid()]=fsalfirst
13371331

13381332
linsolve = cache.linsolve[Threads.threadid()]
13391333

@@ -1358,10 +1352,9 @@ function perform_step!(integrator, cache::ImplicitDeuflhardExtrapolationCache,
13581352
f(k_tmps[Threads.threadid()],
13591353
cache.u_temp3[Threads.threadid()],
13601354
p, t + (j - 1) * dt_int_temp)
1361-
@.. broadcast=false linsolve_tmps[Threads.threadid()]=dt_int_temp *
1362-
k_tmps[Threads.threadid()] -
1355+
@.. broadcast=false linsolve_tmps[Threads.threadid()]=k_tmps[Threads.threadid()] -
13631356
(u_temp3[Threads.threadid()] -
1364-
u_temp4[Threads.threadid()])
1357+
u_temp4[Threads.threadid()])/dt_int_temp
13651358

13661359
linsolve = cache.linsolve[Threads.threadid()]
13671360

@@ -1377,10 +1370,9 @@ function perform_step!(integrator, cache::ImplicitDeuflhardExtrapolationCache,
13771370
end
13781371
cache.linsolve[Threads.threadid()] = linres.cache
13791372

1380-
@.. broadcast=false k_tmps[Threads.threadid()]=-k_tmps[Threads.threadid()]
13811373
@.. broadcast=false T[index + 1]=2 *
13821374
u_temp3[Threads.threadid()] -
1383-
u_temp4[Threads.threadid()] +
1375+
u_temp4[Threads.threadid()] -
13841376
2 * k_tmps[Threads.threadid()] # Explicit Midpoint rule
13851377
@.. broadcast=false u_temp4[Threads.threadid()]=u_temp3[Threads.threadid()]
13861378
@.. broadcast=false u_temp3[Threads.threadid()]=T[index + 1]
@@ -1459,7 +1451,7 @@ function perform_step!(integrator, cache::ImplicitDeuflhardExtrapolationCache,
14591451
jacobian2W!(W[1], integrator.f.mass_matrix, dt_int, J, true)
14601452
integrator.stats.nw += 1
14611453
@.. broadcast=false u_temp2=uprev
1462-
@.. broadcast=false linsolve_tmps[1]=dt_int * fsalfirst
1454+
@.. broadcast=false linsolve_tmps[1]=fsalfirst
14631455

14641456
linsolve = cache.linsolve[1]
14651457
linres = dolinsolve(integrator, linsolve; b = _vec(linsolve_tmps[1]),
@@ -2545,7 +2537,7 @@ function perform_step!(integrator, cache::ImplicitHairerWannerExtrapolationCache
25452537
jacobian2W!(W[1], integrator.f.mass_matrix, dt_int, J, true)
25462538
integrator.stats.nw += 1
25472539
@.. broadcast=false u_temp2=uprev
2548-
@.. broadcast=false linsolve_tmps[1]=dt_int * fsalfirst
2540+
@.. broadcast=false linsolve_tmps[1]=fsalfirst
25492541

25502542
linsolve = cache.linsolve[1]
25512543
if !repeat_step
@@ -2558,8 +2550,7 @@ function perform_step!(integrator, cache::ImplicitHairerWannerExtrapolationCache
25582550
cache.linsolve[1] = linres.cache
25592551

25602552
integrator.stats.nsolve += 1
2561-
@.. broadcast=false k=-k
2562-
@.. broadcast=false u_temp1=u_temp2 + k # Euler starting step
2553+
@.. broadcast=false u_temp1=u_temp2 - k # Euler starting step
25632554
@.. broadcast=false diff1[1]=u_temp1 - u_temp2
25642555
for j in 2:(j_int + 1)
25652556
f(k, cache.u_temp1, p, t + (j - 1) * dt_int)
@@ -2619,8 +2610,7 @@ function perform_step!(integrator, cache::ImplicitHairerWannerExtrapolationCache
26192610
jacobian2W!(W[Threads.threadid()], integrator.f.mass_matrix,
26202611
dt_int_temp, J, true)
26212612
@.. broadcast=false u_temp4[Threads.threadid()]=uprev
2622-
@.. broadcast=false linsolve_tmps[Threads.threadid()]=dt_int_temp *
2623-
fsalfirst
2613+
@.. broadcast=false linsolve_tmps[Threads.threadid()]=fsalfirst
26242614

26252615
linsolve = cache.linsolve[Threads.threadid()]
26262616

@@ -2636,8 +2626,7 @@ function perform_step!(integrator, cache::ImplicitHairerWannerExtrapolationCache
26362626
end
26372627
cache.linsolve[Threads.threadid()] = linres.cache
26382628

2639-
@.. broadcast=false k_tmps[Threads.threadid()]=-k_tmps[Threads.threadid()]
2640-
@.. broadcast=false u_temp3[Threads.threadid()]=u_temp4[Threads.threadid()] +
2629+
@.. broadcast=false u_temp3[Threads.threadid()]=u_temp4[Threads.threadid()] -
26412630
k_tmps[Threads.threadid()] # Euler starting step
26422631
@.. broadcast=false diff1[Threads.threadid()]=u_temp3[Threads.threadid()] -
26432632
u_temp4[Threads.threadid()]
@@ -2711,7 +2700,7 @@ function perform_step!(integrator, cache::ImplicitHairerWannerExtrapolationCache
27112700
dt_int_temp = dt / j_int_temp # Stepsize of the ith internal discretisation
27122701
jacobian2W!(W[tid], integrator.f.mass_matrix, dt_int_temp, J, true)
27132702
@.. broadcast=false u_temp4[tid]=uprev
2714-
@.. broadcast=false linsolvetmp=dt_int_temp * fsalfirst
2703+
@.. broadcast=false linsolvetmp=fsalfirst
27152704

27162705
linsolve = cache.linsolve[tid]
27172706
if !repeat_step
@@ -2723,8 +2712,7 @@ function perform_step!(integrator, cache::ImplicitHairerWannerExtrapolationCache
27232712
end
27242713
cache.linsolve[tid] = linres.cache
27252714

2726-
@.. broadcast=false ktmp=-ktmp
2727-
@.. broadcast=false u_temp3[tid]=u_temp4[tid] + ktmp # Euler starting step
2715+
@.. broadcast=false u_temp3[tid]=u_temp4[tid] - ktmp # Euler starting step
27282716
@.. broadcast=false diff1[tid]=u_temp3[tid] - u_temp4[tid]
27292717
for j in 2:(j_int_temp + 1)
27302718
f(ktmp, cache.u_temp3[tid], p, t + (j - 1) * dt_int_temp)
@@ -2826,7 +2814,7 @@ function perform_step!(integrator, cache::ImplicitHairerWannerExtrapolationCache
28262814
jacobian2W!(W[1], integrator.f.mass_matrix, dt_int, J, true)
28272815
integrator.stats.nw += 1
28282816
@.. broadcast=false u_temp2=uprev
2829-
@.. broadcast=false linsolve_tmps[1]=dt_int * fsalfirst
2817+
@.. broadcast=false linsolve_tmps[1]=fsalfirst
28302818

28312819
linsolve = cache.linsolve[1]
28322820

@@ -2840,12 +2828,11 @@ function perform_step!(integrator, cache::ImplicitHairerWannerExtrapolationCache
28402828
cache.linsolve[1] = linres.cache
28412829

28422830
integrator.stats.nsolve += 1
2843-
@.. broadcast=false k=-k
2844-
@.. broadcast=false u_temp1=u_temp2 + k # Euler starting step
2831+
@.. broadcast=false u_temp1=u_temp2 - k # Euler starting step
28452832
for j in 2:(j_int + 1)
28462833
f(k, cache.u_temp1, p, t + (j - 1) * dt_int)
28472834
OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)
2848-
@.. broadcast=false linsolve_tmps[1]=dt_int * k - (u_temp1 - u_temp2)
2835+
@.. broadcast=false linsolve_tmps[1]=k - (u_temp1 - u_temp2)/dt_int
28492836

28502837
linsolve = cache.linsolve[1]
28512838

@@ -2859,8 +2846,7 @@ function perform_step!(integrator, cache::ImplicitHairerWannerExtrapolationCache
28592846
cache.linsolve[1] = linres.cache
28602847

28612848
integrator.stats.nsolve += 1
2862-
@.. broadcast=false k=-k
2863-
@.. broadcast=false T[n_curr + 1]=2 * u_temp1 - u_temp2 + 2 * k # Explicit Midpoint rule
2849+
@.. broadcast=false T[n_curr + 1]=2 * u_temp1 - u_temp2 - 2 * k # Explicit Midpoint rule
28642850
if (j == j_int + 1)
28652851
@.. broadcast=false T[n_curr + 1]=0.5(T[n_curr + 1] + u_temp2)
28662852
end
@@ -3532,7 +3518,7 @@ function perform_step!(integrator, cache::ImplicitEulerBarycentricExtrapolationC
35323518
jacobian2W!(W[1], integrator.f.mass_matrix, dt_int, J, true)
35333519
integrator.stats.nw += 1
35343520
@.. broadcast=false u_temp2=uprev
3535-
@.. broadcast=false linsolve_tmps[1]=dt_int * fsalfirst
3521+
@.. broadcast=false linsolve_tmps[1]=fsalfirst
35363522

35373523
linsolve = cache.linsolve[1]
35383524

0 commit comments

Comments
 (0)