@@ -3240,7 +3240,7 @@ function perform_step!(integrator, cache::ImplicitEulerBarycentricExtrapolationC
3240
3240
jacobian2W! (W[1 ], integrator. f. mass_matrix, dt_int, J, true )
3241
3241
integrator. stats. nw += 1
3242
3242
@. . broadcast= false u_temp2= uprev
3243
- @. . broadcast= false linsolve_tmps[1 ]= dt_int * fsalfirst
3243
+ @. . broadcast= false linsolve_tmps[1 ]= fsalfirst
3244
3244
3245
3245
linsolve = cache. linsolve[1 ]
3246
3246
@@ -3254,8 +3254,7 @@ function perform_step!(integrator, cache::ImplicitEulerBarycentricExtrapolationC
3254
3254
cache. linsolve[1 ] = linres. cache
3255
3255
3256
3256
integrator. stats. nsolve += 1
3257
- @. . broadcast= false k= - k
3258
- @. . broadcast= false u_temp1= u_temp2 + k # Euler starting step
3257
+ @. . broadcast= false u_temp1= u_temp2 - k # Euler starting step
3259
3258
@. . broadcast= false diff1[1 ]= u_temp1 - u_temp2
3260
3259
for j in 2 : (j_int + 1 )
3261
3260
f (k, cache. u_temp1, p, t + (j - 1 ) * dt_int)
@@ -3314,8 +3313,7 @@ function perform_step!(integrator, cache::ImplicitEulerBarycentricExtrapolationC
3314
3313
jacobian2W! (W[Threads. threadid ()], integrator. f. mass_matrix,
3315
3314
dt_int_temp, J, true )
3316
3315
@. . broadcast= false u_temp4[Threads. threadid ()]= uprev
3317
- @. . broadcast= false linsolve_tmps[Threads. threadid ()]= dt_int_temp *
3318
- fsalfirst
3316
+ @. . broadcast= false linsolve_tmps[Threads. threadid ()]= fsalfirst
3319
3317
3320
3318
linsolve = cache. linsolve[Threads. threadid ()]
3321
3319
@@ -3332,16 +3330,15 @@ function perform_step!(integrator, cache::ImplicitEulerBarycentricExtrapolationC
3332
3330
cache. linsolve[Threads. threadid ()] = linres. cache
3333
3331
3334
3332
@. . broadcast= false k_tmps[Threads. threadid ()]= - k_tmps[Threads. threadid ()]
3335
- @. . broadcast= false u_temp3[Threads. threadid ()]= u_temp4[Threads. threadid ()] +
3333
+ @. . broadcast= false u_temp3[Threads. threadid ()]= u_temp4[Threads. threadid ()] -
3336
3334
k_tmps[Threads. threadid ()] # Euler starting step
3337
3335
@. . broadcast= false diff1[Threads. threadid ()]= u_temp3[Threads. threadid ()] -
3338
3336
u_temp4[Threads. threadid ()]
3339
3337
for j in 2 : (j_int_temp + 1 )
3340
3338
f (k_tmps[Threads. threadid ()],
3341
3339
cache. u_temp3[Threads. threadid ()],
3342
3340
p, t + (j - 1 ) * dt_int_temp)
3343
- @. . broadcast= false linsolve_tmps[Threads. threadid ()]= dt_int_temp *
3344
- k_tmps[Threads. threadid ()]
3341
+ @. . broadcast= false linsolve_tmps[Threads. threadid ()]= k_tmps[Threads. threadid ()]
3345
3342
3346
3343
linsolve = cache. linsolve[Threads. threadid ()]
3347
3344
@@ -3357,8 +3354,7 @@ function perform_step!(integrator, cache::ImplicitEulerBarycentricExtrapolationC
3357
3354
end
3358
3355
cache. linsolve[Threads. threadid ()] = linres. cache
3359
3356
3360
- @. . broadcast= false k_tmps[Threads. threadid ()]= - k_tmps[Threads. threadid ()]
3361
- @. . broadcast= false T[index + 1 ]= u_temp3[Threads. threadid ()] +
3357
+ @. . broadcast= false T[index + 1 ]= u_temp3[Threads. threadid ()] -
3362
3358
k_tmps[Threads. threadid ()] # Explicit Midpoint rule
3363
3359
if (j == j_int_temp + 1 )
3364
3360
@. . broadcast= false T[index + 1 ]= 0.25 (T[index + 1 ] +
@@ -3405,8 +3401,7 @@ function perform_step!(integrator, cache::ImplicitEulerBarycentricExtrapolationC
3405
3401
jacobian2W! (W[Threads. threadid ()], integrator. f. mass_matrix,
3406
3402
dt_int_temp, J, true )
3407
3403
@. . broadcast= false u_temp4[Threads. threadid ()]= uprev
3408
- @. . broadcast= false linsolve_tmps[Threads. threadid ()]= dt_int_temp *
3409
- fsalfirst
3404
+ @. . broadcast= false linsolve_tmps[Threads. threadid ()]= fsalfirst
3410
3405
3411
3406
linsolve = cache. linsolve[Threads. threadid ()]
3412
3407
@@ -3422,17 +3417,15 @@ function perform_step!(integrator, cache::ImplicitEulerBarycentricExtrapolationC
3422
3417
end
3423
3418
cache. linsolve[Threads. threadid ()] = linres. cache
3424
3419
3425
- @. . broadcast= false k_tmps[Threads. threadid ()]= - k_tmps[Threads. threadid ()]
3426
- @. . broadcast= false u_temp3[Threads. threadid ()]= u_temp4[Threads. threadid ()] +
3420
+ @. . broadcast= false u_temp3[Threads. threadid ()]= u_temp4[Threads. threadid ()] -
3427
3421
k_tmps[Threads. threadid ()] # Euler starting step
3428
3422
@. . broadcast= false diff1[Threads. threadid ()]= u_temp3[Threads. threadid ()] -
3429
3423
u_temp4[Threads. threadid ()]
3430
3424
for j in 2 : (j_int_temp + 1 )
3431
3425
f (k_tmps[Threads. threadid ()],
3432
3426
cache. u_temp3[Threads. threadid ()],
3433
3427
p, t + (j - 1 ) * dt_int_temp)
3434
- @. . broadcast= false linsolve_tmps[Threads. threadid ()]= dt_int_temp *
3435
- k_tmps[Threads. threadid ()]
3428
+ @. . broadcast= false linsolve_tmps[Threads. threadid ()]= k_tmps[Threads. threadid ()]
3436
3429
3437
3430
linsolve = cache. linsolve[Threads. threadid ()]
3438
3431
@@ -3448,8 +3441,7 @@ function perform_step!(integrator, cache::ImplicitEulerBarycentricExtrapolationC
3448
3441
end
3449
3442
cache. linsolve[Threads. threadid ()] = linres. cache
3450
3443
3451
- @. . broadcast= false k_tmps[Threads. threadid ()]= - k_tmps[Threads. threadid ()]
3452
- @. . broadcast= false T[index + 1 ]= u_temp3[Threads. threadid ()] +
3444
+ @. . broadcast= false T[index + 1 ]= u_temp3[Threads. threadid ()] -
3453
3445
k_tmps[Threads. threadid ()] # Explicit Midpoint rule
3454
3446
if (j == j_int_temp + 1 )
3455
3447
@. . broadcast= false T[index + 1 ]= 0.25 (T[index + 1 ] +
@@ -3559,16 +3551,15 @@ function perform_step!(integrator, cache::ImplicitEulerBarycentricExtrapolationC
3559
3551
for j in 2 : (j_int + 1 )
3560
3552
f (k, cache. u_temp1, p, t + (j - 1 ) * dt_int)
3561
3553
OrdinaryDiffEqCore. increment_nf! (integrator. stats, 1 )
3562
- @. . broadcast= false linsolve_tmps[1 ]= dt_int * k
3554
+ @. . broadcast= false linsolve_tmps[1 ]= k
3563
3555
3564
3556
linsolve = cache. linsolve[1 ]
3565
3557
linres = dolinsolve (integrator, linsolve; b = _vec (linsolve_tmps[1 ]),
3566
3558
linu = _vec (k))
3567
3559
cache. linsolve[1 ] = linres. cache
3568
3560
3569
3561
integrator. stats. nsolve += 1
3570
- @. . broadcast= false k= - k
3571
- @. . broadcast= false T[n_curr + 1 ]= u_temp1 + k # Explicit Midpoint rule
3562
+ @. . broadcast= false T[n_curr + 1 ]= u_temp1 - k # Explicit Midpoint rule
3572
3563
if (j == j_int + 1 )
3573
3564
@. . broadcast= false T[n_curr + 1 ]= 0.25 (T[n_curr + 1 ] + 2 * u_temp1 +
3574
3565
u_temp2)
0 commit comments