Skip to content

Commit 78e5c03

Browse files
oscardssmithoscardssmith
committed
wtransform for the extrapolation methods
1 parent e5dd0ff commit 78e5c03

File tree

1 file changed

+32
-43
lines changed

1 file changed

+32
-43
lines changed

‎lib/OrdinaryDiffEqExtrapolation/src/extrapolation_perform_step.jl

Lines changed: 32 additions & 43 deletions
Original file line numberDiff line numberDiff line change
@@ -290,13 +290,13 @@ function perform_step!(integrator, cache::ImplicitEulerExtrapolationCache,
290290
calc_J!(J, integrator, cache) # Store the calculated jac as it won't change in internal discretisation
291291
for index in 1:(n_curr + 1)
292292
dt_temp = dt / sequence[index]
293-
jacobian2W!(W[1], integrator.f.mass_matrix, dt_temp, J, false)
293+
jacobian2W!(W[1], integrator.f.mass_matrix, dt_temp, J, true)
294294
integrator.stats.nw += 1
295295
@.. broadcast=false k_tmps[1]=integrator.fsalfirst
296296
@.. broadcast=false u_tmps[1]=uprev
297297

298298
for j in 1:sequence[index]
299-
@.. broadcast=false linsolve_tmps[1]=dt_temp * k_tmps[1]
299+
@.. broadcast=false linsolve_tmps[1]=k_tmps[1]
300300

301301
linsolve = cache.linsolve[1]
302302

@@ -311,9 +311,8 @@ function perform_step!(integrator, cache::ImplicitEulerExtrapolationCache,
311311
cache.linsolve[1] = linres.cache
312312

313313
integrator.stats.nsolve += 1
314-
@.. broadcast=false k_tmps[1]=-k_tmps[1]
315314
@.. broadcast=false u_tmps2[1]=u_tmps[1]
316-
@.. broadcast=false u_tmps[1]=u_tmps[1] + k_tmps[1]
315+
@.. broadcast=false u_tmps[1]=u_tmps[1] - k_tmps[1]
317316
if index <= 2 && j >= 2
318317
# Deuflhard Stability check for initial two sequences
319318
@.. broadcast=false diff2[1]=u_tmps[1] - u_tmps2[1]
@@ -347,12 +346,11 @@ function perform_step!(integrator, cache::ImplicitEulerExtrapolationCache,
347346
dt_temp = dt / sequence[index]
348347
jacobian2W!(
349348
W[Threads.threadid()], integrator.f.mass_matrix, dt_temp, J,
350-
false)
349+
true)
351350
@.. broadcast=false k_tmps[Threads.threadid()]=integrator.fsalfirst
352351
@.. broadcast=false u_tmps[Threads.threadid()]=uprev
353352
for j in 1:sequence[index]
354-
@.. broadcast=false linsolve_tmps[Threads.threadid()]=dt_temp *
355-
k_tmps[Threads.threadid()]
353+
@.. broadcast=false linsolve_tmps[Threads.threadid()]=k_tmps[Threads.threadid()]
356354

357355
linsolve = cache.linsolve[Threads.threadid()]
358356

@@ -369,9 +367,8 @@ function perform_step!(integrator, cache::ImplicitEulerExtrapolationCache,
369367

370368
cache.linsolve[Threads.threadid()] = linres.cache
371369

372-
@.. broadcast=false k_tmps[Threads.threadid()]=-k_tmps[Threads.threadid()]
373370
@.. broadcast=false u_tmps2[Threads.threadid()]=u_tmps[Threads.threadid()]
374-
@.. broadcast=false u_tmps[Threads.threadid()]=u_tmps[Threads.threadid()] +
371+
@.. broadcast=false u_tmps[Threads.threadid()]=u_tmps[Threads.threadid()] -
375372
k_tmps[Threads.threadid()]
376373
if index <= 2 && j >= 2
377374
# Deuflhard Stability check for initial two sequences
@@ -454,7 +451,7 @@ function perform_step!(integrator, cache::ImplicitEulerExtrapolationCache,
454451
@.. broadcast=false u_tmps[1]=uprev
455452

456453
for j in 1:sequence[n_curr + 1]
457-
@.. broadcast=false linsolve_tmps[1]=dt_temp * k_tmps[1]
454+
@.. broadcast=false linsolve_tmps[1]=k_tmps[1]
458455

459456
linsolve = cache.linsolve[1]
460457

@@ -471,8 +468,7 @@ function perform_step!(integrator, cache::ImplicitEulerExtrapolationCache,
471468
cache.linsolve[1] = linres.cache
472469

473470
integrator.stats.nsolve += 1
474-
@.. broadcast=false k_tmps[1]=-k_tmps[1]
475-
@.. broadcast=false u_tmps[1]=u_tmps[1] + k_tmps[1]
471+
@.. broadcast=false u_tmps[1]=u_tmps[1] - k_tmps[1]
476472
f(k_tmps[1], u_tmps[1], p, t + j * dt_temp)
477473
OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)
478474
end
@@ -1174,7 +1170,7 @@ function perform_step!(integrator, cache::ImplicitDeuflhardExtrapolationCache,
11741170
for i in 0:n_curr
11751171
j_int = 4 * subdividing_sequence[i + 1]
11761172
dt_int = dt / j_int # Stepsize of the ith internal discretisation
1177-
jacobian2W!(W[1], integrator.f.mass_matrix, dt_int, J, false)
1173+
jacobian2W!(W[1], integrator.f.mass_matrix, dt_int, J, true)
11781174
integrator.stats.nw += 1
11791175
@.. broadcast=false u_temp2=uprev
11801176
@.. broadcast=false linsolve_tmps[1]=dt_int * fsalfirst
@@ -1247,7 +1243,7 @@ function perform_step!(integrator, cache::ImplicitDeuflhardExtrapolationCache,
12471243
j_int_temp = 4 * subdividing_sequence[index + 1]
12481244
dt_int_temp = dt / j_int_temp # Stepsize of the ith internal discretisation
12491245
jacobian2W!(W[Threads.threadid()], integrator.f.mass_matrix,
1250-
dt_int_temp, J, false)
1246+
dt_int_temp, J, true)
12511247
@.. broadcast=false u_temp4[Threads.threadid()]=uprev
12521248
@.. broadcast=false linsolve_tmps[Threads.threadid()]=dt_int_temp *
12531249
fsalfirst
@@ -1334,7 +1330,7 @@ function perform_step!(integrator, cache::ImplicitDeuflhardExtrapolationCache,
13341330
j_int_temp = 4 * subdividing_sequence[index + 1]
13351331
dt_int_temp = dt / j_int_temp # Stepsize of the ith internal discretisation
13361332
jacobian2W!(W[Threads.threadid()], integrator.f.mass_matrix,
1337-
dt_int_temp, J, false)
1333+
dt_int_temp, J, true)
13381334
@.. broadcast=false u_temp4[Threads.threadid()]=uprev
13391335
@.. broadcast=false linsolve_tmps[Threads.threadid()]=dt_int_temp *
13401336
fsalfirst
@@ -1460,7 +1456,7 @@ function perform_step!(integrator, cache::ImplicitDeuflhardExtrapolationCache,
14601456
# Update cache.T
14611457
j_int = 4 * subdividing_sequence[n_curr + 1]
14621458
dt_int = dt / j_int # Stepsize of the new internal discretisation
1463-
jacobian2W!(W[1], integrator.f.mass_matrix, dt_int, J, false)
1459+
jacobian2W!(W[1], integrator.f.mass_matrix, dt_int, J, true)
14641460
integrator.stats.nw += 1
14651461
@.. broadcast=false u_temp2=uprev
14661462
@.. broadcast=false linsolve_tmps[1]=dt_int * fsalfirst
@@ -1471,8 +1467,7 @@ function perform_step!(integrator, cache::ImplicitDeuflhardExtrapolationCache,
14711467
cache.linsolve[1] = linres.cache
14721468

14731469
integrator.stats.nsolve += 1
1474-
@.. broadcast=false k=-k
1475-
@.. broadcast=false u_temp1=u_temp2 + k # Euler starting step
1470+
@.. broadcast=false u_temp1=u_temp2 - k # Euler starting step
14761471
for j in 2:j_int
14771472
f(k, cache.u_temp1, p, t + (j - 1) * dt_int)
14781473
OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)
@@ -1484,8 +1479,7 @@ function perform_step!(integrator, cache::ImplicitDeuflhardExtrapolationCache,
14841479
cache.linsolve[1] = linres.cache
14851480

14861481
integrator.stats.nsolve += 1
1487-
@.. broadcast=false k=-k
1488-
@.. broadcast=false T[n_curr + 1]=2 * u_temp1 - u_temp2 + 2 * k # Explicit Midpoint rule
1482+
@.. broadcast=false T[n_curr + 1]=2 * u_temp1 - u_temp2 - 2 * k # Explicit Midpoint rule
14891483
@.. broadcast=false u_temp2=u_temp1
14901484
@.. broadcast=false u_temp1=T[n_curr + 1]
14911485
end
@@ -2548,7 +2542,7 @@ function perform_step!(integrator, cache::ImplicitHairerWannerExtrapolationCache
25482542
for i in 0:n_curr
25492543
j_int = 4 * subdividing_sequence[i + 1]
25502544
dt_int = dt / j_int # Stepsize of the ith internal discretisation
2551-
jacobian2W!(W[1], integrator.f.mass_matrix, dt_int, J, false)
2545+
jacobian2W!(W[1], integrator.f.mass_matrix, dt_int, J, true)
25522546
integrator.stats.nw += 1
25532547
@.. broadcast=false u_temp2=uprev
25542548
@.. broadcast=false linsolve_tmps[1]=dt_int * fsalfirst
@@ -2570,7 +2564,7 @@ function perform_step!(integrator, cache::ImplicitHairerWannerExtrapolationCache
25702564
for j in 2:(j_int + 1)
25712565
f(k, cache.u_temp1, p, t + (j - 1) * dt_int)
25722566
OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)
2573-
@.. broadcast=false linsolve_tmps[1]=dt_int * k - (u_temp1 - u_temp2)
2567+
@.. broadcast=false linsolve_tmps[1]=k - (u_temp1 - u_temp2)/dt_int
25742568

25752569
linsolve = cache.linsolve[1]
25762570

@@ -2584,8 +2578,7 @@ function perform_step!(integrator, cache::ImplicitHairerWannerExtrapolationCache
25842578
cache.linsolve[1] = linres.cache
25852579

25862580
integrator.stats.nsolve += 1
2587-
@.. broadcast=false k=-k
2588-
@.. broadcast=false T[i + 1]=2 * u_temp1 - u_temp2 + 2 * k # Explicit Midpoint rule
2581+
@.. broadcast=false T[i + 1]=2 * u_temp1 - u_temp2 - 2 * k # Explicit Midpoint rule
25892582
if (j == j_int + 1)
25902583
@.. broadcast=false T[i + 1]=0.5(T[i + 1] + u_temp2)
25912584
end
@@ -2624,7 +2617,7 @@ function perform_step!(integrator, cache::ImplicitHairerWannerExtrapolationCache
26242617
j_int_temp = 4 * subdividing_sequence[index + 1]
26252618
dt_int_temp = dt / j_int_temp # Stepsize of the ith internal discretisation
26262619
jacobian2W!(W[Threads.threadid()], integrator.f.mass_matrix,
2627-
dt_int_temp, J, false)
2620+
dt_int_temp, J, true)
26282621
@.. broadcast=false u_temp4[Threads.threadid()]=uprev
26292622
@.. broadcast=false linsolve_tmps[Threads.threadid()]=dt_int_temp *
26302623
fsalfirst
@@ -2652,10 +2645,9 @@ function perform_step!(integrator, cache::ImplicitHairerWannerExtrapolationCache
26522645
f(k_tmps[Threads.threadid()],
26532646
cache.u_temp3[Threads.threadid()],
26542647
p, t + (j - 1) * dt_int_temp)
2655-
@.. broadcast=false linsolve_tmps[Threads.threadid()]=dt_int_temp *
2656-
k_tmps[Threads.threadid()] -
2648+
@.. broadcast=false linsolve_tmps[Threads.threadid()]= k_tmps[Threads.threadid()] -
26572649
(u_temp3[Threads.threadid()] -
2658-
u_temp4[Threads.threadid()])
2650+
u_temp4[Threads.threadid()]) / dt_int_temp
26592651

26602652
linsolve = cache.linsolve[Threads.threadid()]
26612653
if !repeat_step && j == 1
@@ -2670,10 +2662,9 @@ function perform_step!(integrator, cache::ImplicitHairerWannerExtrapolationCache
26702662
end
26712663
cache.linsolve[Threads.threadid()] = linres.cache
26722664

2673-
@.. broadcast=false k_tmps[Threads.threadid()]=-k_tmps[Threads.threadid()]
26742665
@.. broadcast=false T[index + 1]=2 *
26752666
u_temp3[Threads.threadid()] -
2676-
u_temp4[Threads.threadid()] +
2667+
u_temp4[Threads.threadid()] -
26772668
2 * k_tmps[Threads.threadid()] # Explicit Midpoint rule
26782669
if (j == j_int_temp + 1)
26792670
@.. broadcast=false T[index + 1]=0.5(T[index + 1] +
@@ -2718,7 +2709,7 @@ function perform_step!(integrator, cache::ImplicitHairerWannerExtrapolationCache
27182709
index == -1 && continue
27192710
j_int_temp = 4 * subdividing_sequence[index + 1]
27202711
dt_int_temp = dt / j_int_temp # Stepsize of the ith internal discretisation
2721-
jacobian2W!(W[tid], integrator.f.mass_matrix, dt_int_temp, J, false)
2712+
jacobian2W!(W[tid], integrator.f.mass_matrix, dt_int_temp, J, true)
27222713
@.. broadcast=false u_temp4[tid]=uprev
27232714
@.. broadcast=false linsolvetmp=dt_int_temp * fsalfirst
27242715

@@ -2737,8 +2728,8 @@ function perform_step!(integrator, cache::ImplicitHairerWannerExtrapolationCache
27372728
@.. broadcast=false diff1[tid]=u_temp3[tid] - u_temp4[tid]
27382729
for j in 2:(j_int_temp + 1)
27392730
f(ktmp, cache.u_temp3[tid], p, t + (j - 1) * dt_int_temp)
2740-
@.. broadcast=false linsolvetmp=dt_int_temp * ktmp -
2741-
(u_temp3[tid] - u_temp4[tid])
2731+
@.. broadcast=false linsolvetmp=ktmp -
2732+
(u_temp3[tid] - u_temp4[tid])/dt_int_temp
27422733

27432734
linsolve = cache.linsolve[tid]
27442735

@@ -2753,9 +2744,8 @@ function perform_step!(integrator, cache::ImplicitHairerWannerExtrapolationCache
27532744
end
27542745
cache.linsolve[tid] = linres.cache
27552746

2756-
@.. broadcast=false ktmp=-ktmp
27572747
@.. broadcast=false T[index + 1]=2 * u_temp3[tid] -
2758-
u_temp4[tid] + 2 * ktmp # Explicit Midpoint rule
2748+
u_temp4[tid] - 2 * ktmp # Explicit Midpoint rule
27592749
if (j == j_int_temp + 1)
27602750
@.. broadcast=false T[index + 1]=0.5(T[index + 1] +
27612751
u_temp4[tid])
@@ -2833,7 +2823,7 @@ function perform_step!(integrator, cache::ImplicitHairerWannerExtrapolationCache
28332823
# Update cache.T
28342824
j_int = 4 * subdividing_sequence[n_curr + 1]
28352825
dt_int = dt / j_int # Stepsize of the new internal discretisation
2836-
jacobian2W!(W[1], integrator.f.mass_matrix, dt_int, J, false)
2826+
jacobian2W!(W[1], integrator.f.mass_matrix, dt_int, J, true)
28372827
integrator.stats.nw += 1
28382828
@.. broadcast=false u_temp2=uprev
28392829
@.. broadcast=false linsolve_tmps[1]=dt_int * fsalfirst
@@ -3247,7 +3237,7 @@ function perform_step!(integrator, cache::ImplicitEulerBarycentricExtrapolationC
32473237
for i in 0:n_curr
32483238
j_int = sequence_factor * subdividing_sequence[i + 1]
32493239
dt_int = dt / j_int # Stepsize of the ith internal discretisation
3250-
jacobian2W!(W[1], integrator.f.mass_matrix, dt_int, J, false)
3240+
jacobian2W!(W[1], integrator.f.mass_matrix, dt_int, J, true)
32513241
integrator.stats.nw += 1
32523242
@.. broadcast=false u_temp2=uprev
32533243
@.. broadcast=false linsolve_tmps[1]=dt_int * fsalfirst
@@ -3270,7 +3260,7 @@ function perform_step!(integrator, cache::ImplicitEulerBarycentricExtrapolationC
32703260
for j in 2:(j_int + 1)
32713261
f(k, cache.u_temp1, p, t + (j - 1) * dt_int)
32723262
OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1)
3273-
@.. broadcast=false linsolve_tmps[1]=dt_int * k
3263+
@.. broadcast=false linsolve_tmps[1]=k
32743264

32753265
linsolve = cache.linsolve[1]
32763266
if !repeat_step && j == 1
@@ -3283,8 +3273,7 @@ function perform_step!(integrator, cache::ImplicitEulerBarycentricExtrapolationC
32833273
cache.linsolve[1] = linres.cache
32843274

32853275
integrator.stats.nsolve += 1
3286-
@.. broadcast=false k=-k
3287-
@.. broadcast=false T[i + 1]=u_temp1 + k
3276+
@.. broadcast=false T[i + 1]=u_temp1 - k
32883277
if (j == j_int + 1)
32893278
@.. broadcast=false T[i + 1]=0.25(T[i + 1] + 2 * u_temp1 + u_temp2)
32903279
end
@@ -3323,7 +3312,7 @@ function perform_step!(integrator, cache::ImplicitEulerBarycentricExtrapolationC
33233312
j_int_temp = sequence_factor * subdividing_sequence[index + 1]
33243313
dt_int_temp = dt / j_int_temp # Stepsize of the ith internal discretisation
33253314
jacobian2W!(W[Threads.threadid()], integrator.f.mass_matrix,
3326-
dt_int_temp, J, false)
3315+
dt_int_temp, J, true)
33273316
@.. broadcast=false u_temp4[Threads.threadid()]=uprev
33283317
@.. broadcast=false linsolve_tmps[Threads.threadid()]=dt_int_temp *
33293318
fsalfirst
@@ -3414,7 +3403,7 @@ function perform_step!(integrator, cache::ImplicitEulerBarycentricExtrapolationC
34143403
j_int_temp = sequence_factor * subdividing_sequence[index + 1]
34153404
dt_int_temp = dt / j_int_temp # Stepsize of the ith internal discretisation
34163405
jacobian2W!(W[Threads.threadid()], integrator.f.mass_matrix,
3417-
dt_int_temp, J, false)
3406+
dt_int_temp, J, true)
34183407
@.. broadcast=false u_temp4[Threads.threadid()]=uprev
34193408
@.. broadcast=false linsolve_tmps[Threads.threadid()]=dt_int_temp *
34203409
fsalfirst
@@ -3548,7 +3537,7 @@ function perform_step!(integrator, cache::ImplicitEulerBarycentricExtrapolationC
35483537
# Update cache.T
35493538
j_int = sequence_factor * subdividing_sequence[n_curr + 1]
35503539
dt_int = dt / j_int # Stepsize of the new internal discretisation
3551-
jacobian2W!(W[1], integrator.f.mass_matrix, dt_int, J, false)
3540+
jacobian2W!(W[1], integrator.f.mass_matrix, dt_int, J, true)
35523541
integrator.stats.nw += 1
35533542
@.. broadcast=false u_temp2=uprev
35543543
@.. broadcast=false linsolve_tmps[1]=dt_int * fsalfirst

0 commit comments

Comments
 (0)