diff --git a/lib/OrdinaryDiffEqExtrapolation/src/extrapolation_perform_step.jl b/lib/OrdinaryDiffEqExtrapolation/src/extrapolation_perform_step.jl index ad26439864..5c33ce5ed3 100644 --- a/lib/OrdinaryDiffEqExtrapolation/src/extrapolation_perform_step.jl +++ b/lib/OrdinaryDiffEqExtrapolation/src/extrapolation_perform_step.jl @@ -290,13 +290,13 @@ function perform_step!(integrator, cache::ImplicitEulerExtrapolationCache, calc_J!(J, integrator, cache) # Store the calculated jac as it won't change in internal discretisation for index in 1:(n_curr + 1) dt_temp = dt / sequence[index] - jacobian2W!(W[1], integrator.f.mass_matrix, dt_temp, J, false) + jacobian2W!(W[1], integrator.f.mass_matrix, dt_temp, J, true) integrator.stats.nw += 1 @.. broadcast=false k_tmps[1]=integrator.fsalfirst @.. broadcast=false u_tmps[1]=uprev for j in 1:sequence[index] - @.. broadcast=false linsolve_tmps[1]=dt_temp * k_tmps[1] + @.. broadcast=false linsolve_tmps[1]=k_tmps[1] linsolve = cache.linsolve[1] @@ -311,9 +311,8 @@ function perform_step!(integrator, cache::ImplicitEulerExtrapolationCache, cache.linsolve[1] = linres.cache integrator.stats.nsolve += 1 - @.. broadcast=false k_tmps[1]=-k_tmps[1] @.. broadcast=false u_tmps2[1]=u_tmps[1] - @.. broadcast=false u_tmps[1]=u_tmps[1] + k_tmps[1] + @.. broadcast=false u_tmps[1]=u_tmps[1] - k_tmps[1] if index <= 2 && j >= 2 # Deuflhard Stability check for initial two sequences @.. broadcast=false diff2[1]=u_tmps[1] - u_tmps2[1] @@ -347,12 +346,11 @@ function perform_step!(integrator, cache::ImplicitEulerExtrapolationCache, dt_temp = dt / sequence[index] jacobian2W!( W[Threads.threadid()], integrator.f.mass_matrix, dt_temp, J, - false) + true) @.. broadcast=false k_tmps[Threads.threadid()]=integrator.fsalfirst @.. broadcast=false u_tmps[Threads.threadid()]=uprev for j in 1:sequence[index] - @.. broadcast=false linsolve_tmps[Threads.threadid()]=dt_temp * - k_tmps[Threads.threadid()] + @.. broadcast=false linsolve_tmps[Threads.threadid()]=k_tmps[Threads.threadid()] linsolve = cache.linsolve[Threads.threadid()] @@ -369,9 +367,8 @@ function perform_step!(integrator, cache::ImplicitEulerExtrapolationCache, cache.linsolve[Threads.threadid()] = linres.cache - @.. broadcast=false k_tmps[Threads.threadid()]=-k_tmps[Threads.threadid()] @.. broadcast=false u_tmps2[Threads.threadid()]=u_tmps[Threads.threadid()] - @.. broadcast=false u_tmps[Threads.threadid()]=u_tmps[Threads.threadid()] + + @.. broadcast=false u_tmps[Threads.threadid()]=u_tmps[Threads.threadid()] - k_tmps[Threads.threadid()] if index <= 2 && j >= 2 # Deuflhard Stability check for initial two sequences @@ -454,7 +451,7 @@ function perform_step!(integrator, cache::ImplicitEulerExtrapolationCache, @.. broadcast=false u_tmps[1]=uprev for j in 1:sequence[n_curr + 1] - @.. broadcast=false linsolve_tmps[1]=dt_temp * k_tmps[1] + @.. broadcast=false linsolve_tmps[1]=k_tmps[1] linsolve = cache.linsolve[1] @@ -471,8 +468,7 @@ function perform_step!(integrator, cache::ImplicitEulerExtrapolationCache, cache.linsolve[1] = linres.cache integrator.stats.nsolve += 1 - @.. broadcast=false k_tmps[1]=-k_tmps[1] - @.. broadcast=false u_tmps[1]=u_tmps[1] + k_tmps[1] + @.. broadcast=false u_tmps[1]=u_tmps[1] - k_tmps[1] f(k_tmps[1], u_tmps[1], p, t + j * dt_temp) OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1) end @@ -1174,10 +1170,10 @@ function perform_step!(integrator, cache::ImplicitDeuflhardExtrapolationCache, for i in 0:n_curr j_int = 4 * subdividing_sequence[i + 1] dt_int = dt / j_int # Stepsize of the ith internal discretisation - jacobian2W!(W[1], integrator.f.mass_matrix, dt_int, J, false) + jacobian2W!(W[1], integrator.f.mass_matrix, dt_int, J, true) integrator.stats.nw += 1 @.. broadcast=false u_temp2=uprev - @.. broadcast=false linsolve_tmps[1]=dt_int * fsalfirst + @.. broadcast=false linsolve_tmps[1]=fsalfirst linsolve = cache.linsolve[1] @@ -1192,13 +1188,12 @@ function perform_step!(integrator, cache::ImplicitDeuflhardExtrapolationCache, cache.linsolve[1] = linres.cache integrator.stats.nsolve += 1 - @.. broadcast=false k=-k - @.. broadcast=false u_temp1=u_temp2 + k # Euler starting step + @.. broadcast=false u_temp1=u_temp2 - k # Euler starting step @.. broadcast=false diff1[1]=u_temp1 - u_temp2 for j in 2:j_int f(k, cache.u_temp1, p, t + (j - 1) * dt_int) OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1) - @.. broadcast=false linsolve_tmps[1]=dt_int * k - (u_temp1 - u_temp2) + @.. broadcast=false linsolve_tmps[1]=k - (u_temp1 - u_temp2)/dt_int linsolve = cache.linsolve[1] @@ -1212,8 +1207,7 @@ function perform_step!(integrator, cache::ImplicitDeuflhardExtrapolationCache, cache.linsolve[1] = linres.cache integrator.stats.nsolve += 1 - @.. broadcast=false k=-k - @.. broadcast=false T[i + 1]=2 * u_temp1 - u_temp2 + 2 * k # Explicit Midpoint rule + @.. broadcast=false T[i + 1]=2 * u_temp1 - u_temp2 - 2 * k # Explicit Midpoint rule @.. broadcast=false u_temp2=u_temp1 @.. broadcast=false u_temp1=T[i + 1] if (i <= 1) @@ -1247,10 +1241,9 @@ function perform_step!(integrator, cache::ImplicitDeuflhardExtrapolationCache, j_int_temp = 4 * subdividing_sequence[index + 1] dt_int_temp = dt / j_int_temp # Stepsize of the ith internal discretisation jacobian2W!(W[Threads.threadid()], integrator.f.mass_matrix, - dt_int_temp, J, false) + dt_int_temp, J, true) @.. broadcast=false u_temp4[Threads.threadid()]=uprev - @.. broadcast=false linsolve_tmps[Threads.threadid()]=dt_int_temp * - fsalfirst + @.. broadcast=false linsolve_tmps[Threads.threadid()]=fsalfirst linsolve = cache.linsolve[Threads.threadid()] @@ -1275,10 +1268,9 @@ function perform_step!(integrator, cache::ImplicitDeuflhardExtrapolationCache, f(k_tmps[Threads.threadid()], cache.u_temp3[Threads.threadid()], p, t + (j - 1) * dt_int_temp) - @.. broadcast=false linsolve_tmps[Threads.threadid()]=dt_int_temp * - k_tmps[Threads.threadid()] - + @.. broadcast=false linsolve_tmps[Threads.threadid()]=k_tmps[Threads.threadid()] - (u_temp3[Threads.threadid()] - - u_temp4[Threads.threadid()]) + u_temp4[Threads.threadid()])/dt_int_temp linsolve = cache.linsolve[Threads.threadid()] @@ -1294,10 +1286,9 @@ function perform_step!(integrator, cache::ImplicitDeuflhardExtrapolationCache, end cache.linsolve[Threads.threadid()] = linres.cache - @.. broadcast=false k_tmps[Threads.threadid()]=-k_tmps[Threads.threadid()] @.. broadcast=false T[index + 1]=2 * u_temp3[Threads.threadid()] - - u_temp4[Threads.threadid()] + + u_temp4[Threads.threadid()] - 2 * k_tmps[Threads.threadid()] # Explicit Midpoint rule @.. broadcast=false u_temp4[Threads.threadid()]=u_temp3[Threads.threadid()] @.. broadcast=false u_temp3[Threads.threadid()]=T[index + 1] @@ -1334,10 +1325,9 @@ function perform_step!(integrator, cache::ImplicitDeuflhardExtrapolationCache, j_int_temp = 4 * subdividing_sequence[index + 1] dt_int_temp = dt / j_int_temp # Stepsize of the ith internal discretisation jacobian2W!(W[Threads.threadid()], integrator.f.mass_matrix, - dt_int_temp, J, false) + dt_int_temp, J, true) @.. broadcast=false u_temp4[Threads.threadid()]=uprev - @.. broadcast=false linsolve_tmps[Threads.threadid()]=dt_int_temp * - fsalfirst + @.. broadcast=false linsolve_tmps[Threads.threadid()]=fsalfirst linsolve = cache.linsolve[Threads.threadid()] @@ -1362,10 +1352,9 @@ function perform_step!(integrator, cache::ImplicitDeuflhardExtrapolationCache, f(k_tmps[Threads.threadid()], cache.u_temp3[Threads.threadid()], p, t + (j - 1) * dt_int_temp) - @.. broadcast=false linsolve_tmps[Threads.threadid()]=dt_int_temp * - k_tmps[Threads.threadid()] - + @.. broadcast=false linsolve_tmps[Threads.threadid()]=k_tmps[Threads.threadid()] - (u_temp3[Threads.threadid()] - - u_temp4[Threads.threadid()]) + u_temp4[Threads.threadid()])/dt_int_temp linsolve = cache.linsolve[Threads.threadid()] @@ -1381,10 +1370,9 @@ function perform_step!(integrator, cache::ImplicitDeuflhardExtrapolationCache, end cache.linsolve[Threads.threadid()] = linres.cache - @.. broadcast=false k_tmps[Threads.threadid()]=-k_tmps[Threads.threadid()] @.. broadcast=false T[index + 1]=2 * u_temp3[Threads.threadid()] - - u_temp4[Threads.threadid()] + + u_temp4[Threads.threadid()] - 2 * k_tmps[Threads.threadid()] # Explicit Midpoint rule @.. broadcast=false u_temp4[Threads.threadid()]=u_temp3[Threads.threadid()] @.. broadcast=false u_temp3[Threads.threadid()]=T[index + 1] @@ -1460,10 +1448,10 @@ function perform_step!(integrator, cache::ImplicitDeuflhardExtrapolationCache, # Update cache.T j_int = 4 * subdividing_sequence[n_curr + 1] dt_int = dt / j_int # Stepsize of the new internal discretisation - jacobian2W!(W[1], integrator.f.mass_matrix, dt_int, J, false) + jacobian2W!(W[1], integrator.f.mass_matrix, dt_int, J, true) integrator.stats.nw += 1 @.. broadcast=false u_temp2=uprev - @.. broadcast=false linsolve_tmps[1]=dt_int * fsalfirst + @.. broadcast=false linsolve_tmps[1]=fsalfirst linsolve = cache.linsolve[1] linres = dolinsolve(integrator, linsolve; b = _vec(linsolve_tmps[1]), @@ -1471,8 +1459,7 @@ function perform_step!(integrator, cache::ImplicitDeuflhardExtrapolationCache, cache.linsolve[1] = linres.cache integrator.stats.nsolve += 1 - @.. broadcast=false k=-k - @.. broadcast=false u_temp1=u_temp2 + k # Euler starting step + @.. broadcast=false u_temp1=u_temp2 - k # Euler starting step for j in 2:j_int f(k, cache.u_temp1, p, t + (j - 1) * dt_int) OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1) @@ -1484,8 +1471,7 @@ function perform_step!(integrator, cache::ImplicitDeuflhardExtrapolationCache, cache.linsolve[1] = linres.cache integrator.stats.nsolve += 1 - @.. broadcast=false k=-k - @.. broadcast=false T[n_curr + 1]=2 * u_temp1 - u_temp2 + 2 * k # Explicit Midpoint rule + @.. broadcast=false T[n_curr + 1]=2 * u_temp1 - u_temp2 - 2 * k # Explicit Midpoint rule @.. broadcast=false u_temp2=u_temp1 @.. broadcast=false u_temp1=T[n_curr + 1] end @@ -2548,10 +2534,10 @@ function perform_step!(integrator, cache::ImplicitHairerWannerExtrapolationCache for i in 0:n_curr j_int = 4 * subdividing_sequence[i + 1] dt_int = dt / j_int # Stepsize of the ith internal discretisation - jacobian2W!(W[1], integrator.f.mass_matrix, dt_int, J, false) + jacobian2W!(W[1], integrator.f.mass_matrix, dt_int, J, true) integrator.stats.nw += 1 @.. broadcast=false u_temp2=uprev - @.. broadcast=false linsolve_tmps[1]=dt_int * fsalfirst + @.. broadcast=false linsolve_tmps[1]=fsalfirst linsolve = cache.linsolve[1] if !repeat_step @@ -2564,13 +2550,12 @@ function perform_step!(integrator, cache::ImplicitHairerWannerExtrapolationCache cache.linsolve[1] = linres.cache integrator.stats.nsolve += 1 - @.. broadcast=false k=-k - @.. broadcast=false u_temp1=u_temp2 + k # Euler starting step + @.. broadcast=false u_temp1=u_temp2 - k # Euler starting step @.. broadcast=false diff1[1]=u_temp1 - u_temp2 for j in 2:(j_int + 1) f(k, cache.u_temp1, p, t + (j - 1) * dt_int) OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1) - @.. broadcast=false linsolve_tmps[1]=dt_int * k - (u_temp1 - u_temp2) + @.. broadcast=false linsolve_tmps[1]=k - (u_temp1 - u_temp2)/dt_int linsolve = cache.linsolve[1] @@ -2584,8 +2569,7 @@ function perform_step!(integrator, cache::ImplicitHairerWannerExtrapolationCache cache.linsolve[1] = linres.cache integrator.stats.nsolve += 1 - @.. broadcast=false k=-k - @.. broadcast=false T[i + 1]=2 * u_temp1 - u_temp2 + 2 * k # Explicit Midpoint rule + @.. broadcast=false T[i + 1]=2 * u_temp1 - u_temp2 - 2 * k # Explicit Midpoint rule if (j == j_int + 1) @.. broadcast=false T[i + 1]=0.5(T[i + 1] + u_temp2) end @@ -2624,10 +2608,9 @@ function perform_step!(integrator, cache::ImplicitHairerWannerExtrapolationCache j_int_temp = 4 * subdividing_sequence[index + 1] dt_int_temp = dt / j_int_temp # Stepsize of the ith internal discretisation jacobian2W!(W[Threads.threadid()], integrator.f.mass_matrix, - dt_int_temp, J, false) + dt_int_temp, J, true) @.. broadcast=false u_temp4[Threads.threadid()]=uprev - @.. broadcast=false linsolve_tmps[Threads.threadid()]=dt_int_temp * - fsalfirst + @.. broadcast=false linsolve_tmps[Threads.threadid()]=fsalfirst linsolve = cache.linsolve[Threads.threadid()] @@ -2643,8 +2626,7 @@ function perform_step!(integrator, cache::ImplicitHairerWannerExtrapolationCache end cache.linsolve[Threads.threadid()] = linres.cache - @.. broadcast=false k_tmps[Threads.threadid()]=-k_tmps[Threads.threadid()] - @.. broadcast=false u_temp3[Threads.threadid()]=u_temp4[Threads.threadid()] + + @.. broadcast=false u_temp3[Threads.threadid()]=u_temp4[Threads.threadid()] - k_tmps[Threads.threadid()] # Euler starting step @.. broadcast=false diff1[Threads.threadid()]=u_temp3[Threads.threadid()] - u_temp4[Threads.threadid()] @@ -2652,10 +2634,9 @@ function perform_step!(integrator, cache::ImplicitHairerWannerExtrapolationCache f(k_tmps[Threads.threadid()], cache.u_temp3[Threads.threadid()], p, t + (j - 1) * dt_int_temp) - @.. broadcast=false linsolve_tmps[Threads.threadid()]=dt_int_temp * - k_tmps[Threads.threadid()] - + @.. broadcast=false linsolve_tmps[Threads.threadid()]= k_tmps[Threads.threadid()] - (u_temp3[Threads.threadid()] - - u_temp4[Threads.threadid()]) + u_temp4[Threads.threadid()]) / dt_int_temp linsolve = cache.linsolve[Threads.threadid()] if !repeat_step && j == 1 @@ -2670,10 +2651,9 @@ function perform_step!(integrator, cache::ImplicitHairerWannerExtrapolationCache end cache.linsolve[Threads.threadid()] = linres.cache - @.. broadcast=false k_tmps[Threads.threadid()]=-k_tmps[Threads.threadid()] @.. broadcast=false T[index + 1]=2 * u_temp3[Threads.threadid()] - - u_temp4[Threads.threadid()] + + u_temp4[Threads.threadid()] - 2 * k_tmps[Threads.threadid()] # Explicit Midpoint rule if (j == j_int_temp + 1) @.. broadcast=false T[index + 1]=0.5(T[index + 1] + @@ -2718,9 +2698,9 @@ function perform_step!(integrator, cache::ImplicitHairerWannerExtrapolationCache index == -1 && continue j_int_temp = 4 * subdividing_sequence[index + 1] dt_int_temp = dt / j_int_temp # Stepsize of the ith internal discretisation - jacobian2W!(W[tid], integrator.f.mass_matrix, dt_int_temp, J, false) + jacobian2W!(W[tid], integrator.f.mass_matrix, dt_int_temp, J, true) @.. broadcast=false u_temp4[tid]=uprev - @.. broadcast=false linsolvetmp=dt_int_temp * fsalfirst + @.. broadcast=false linsolvetmp=fsalfirst linsolve = cache.linsolve[tid] if !repeat_step @@ -2732,13 +2712,12 @@ function perform_step!(integrator, cache::ImplicitHairerWannerExtrapolationCache end cache.linsolve[tid] = linres.cache - @.. broadcast=false ktmp=-ktmp - @.. broadcast=false u_temp3[tid]=u_temp4[tid] + ktmp # Euler starting step + @.. broadcast=false u_temp3[tid]=u_temp4[tid] - ktmp # Euler starting step @.. broadcast=false diff1[tid]=u_temp3[tid] - u_temp4[tid] for j in 2:(j_int_temp + 1) f(ktmp, cache.u_temp3[tid], p, t + (j - 1) * dt_int_temp) - @.. broadcast=false linsolvetmp=dt_int_temp * ktmp - - (u_temp3[tid] - u_temp4[tid]) + @.. broadcast=false linsolvetmp=ktmp - + (u_temp3[tid] - u_temp4[tid])/dt_int_temp linsolve = cache.linsolve[tid] @@ -2753,9 +2732,8 @@ function perform_step!(integrator, cache::ImplicitHairerWannerExtrapolationCache end cache.linsolve[tid] = linres.cache - @.. broadcast=false ktmp=-ktmp @.. broadcast=false T[index + 1]=2 * u_temp3[tid] - - u_temp4[tid] + 2 * ktmp # Explicit Midpoint rule + u_temp4[tid] - 2 * ktmp # Explicit Midpoint rule if (j == j_int_temp + 1) @.. broadcast=false T[index + 1]=0.5(T[index + 1] + u_temp4[tid]) @@ -2833,10 +2811,10 @@ function perform_step!(integrator, cache::ImplicitHairerWannerExtrapolationCache # Update cache.T j_int = 4 * subdividing_sequence[n_curr + 1] dt_int = dt / j_int # Stepsize of the new internal discretisation - jacobian2W!(W[1], integrator.f.mass_matrix, dt_int, J, false) + jacobian2W!(W[1], integrator.f.mass_matrix, dt_int, J, true) integrator.stats.nw += 1 @.. broadcast=false u_temp2=uprev - @.. broadcast=false linsolve_tmps[1]=dt_int * fsalfirst + @.. broadcast=false linsolve_tmps[1]=fsalfirst linsolve = cache.linsolve[1] @@ -2850,12 +2828,11 @@ function perform_step!(integrator, cache::ImplicitHairerWannerExtrapolationCache cache.linsolve[1] = linres.cache integrator.stats.nsolve += 1 - @.. broadcast=false k=-k - @.. broadcast=false u_temp1=u_temp2 + k # Euler starting step + @.. broadcast=false u_temp1=u_temp2 - k # Euler starting step for j in 2:(j_int + 1) f(k, cache.u_temp1, p, t + (j - 1) * dt_int) OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1) - @.. broadcast=false linsolve_tmps[1]=dt_int * k - (u_temp1 - u_temp2) + @.. broadcast=false linsolve_tmps[1]=k - (u_temp1 - u_temp2)/dt_int linsolve = cache.linsolve[1] @@ -2869,8 +2846,7 @@ function perform_step!(integrator, cache::ImplicitHairerWannerExtrapolationCache cache.linsolve[1] = linres.cache integrator.stats.nsolve += 1 - @.. broadcast=false k=-k - @.. broadcast=false T[n_curr + 1]=2 * u_temp1 - u_temp2 + 2 * k # Explicit Midpoint rule + @.. broadcast=false T[n_curr + 1]=2 * u_temp1 - u_temp2 - 2 * k # Explicit Midpoint rule if (j == j_int + 1) @.. broadcast=false T[n_curr + 1]=0.5(T[n_curr + 1] + u_temp2) end @@ -3247,10 +3223,10 @@ function perform_step!(integrator, cache::ImplicitEulerBarycentricExtrapolationC for i in 0:n_curr j_int = sequence_factor * subdividing_sequence[i + 1] dt_int = dt / j_int # Stepsize of the ith internal discretisation - jacobian2W!(W[1], integrator.f.mass_matrix, dt_int, J, false) + jacobian2W!(W[1], integrator.f.mass_matrix, dt_int, J, true) integrator.stats.nw += 1 @.. broadcast=false u_temp2=uprev - @.. broadcast=false linsolve_tmps[1]=dt_int * fsalfirst + @.. broadcast=false linsolve_tmps[1]=fsalfirst linsolve = cache.linsolve[1] @@ -3264,13 +3240,12 @@ function perform_step!(integrator, cache::ImplicitEulerBarycentricExtrapolationC cache.linsolve[1] = linres.cache integrator.stats.nsolve += 1 - @.. broadcast=false k=-k - @.. broadcast=false u_temp1=u_temp2 + k # Euler starting step + @.. broadcast=false u_temp1=u_temp2 - k # Euler starting step @.. broadcast=false diff1[1]=u_temp1 - u_temp2 for j in 2:(j_int + 1) f(k, cache.u_temp1, p, t + (j - 1) * dt_int) OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1) - @.. broadcast=false linsolve_tmps[1]=dt_int * k + @.. broadcast=false linsolve_tmps[1]=k linsolve = cache.linsolve[1] if !repeat_step && j == 1 @@ -3283,8 +3258,7 @@ function perform_step!(integrator, cache::ImplicitEulerBarycentricExtrapolationC cache.linsolve[1] = linres.cache integrator.stats.nsolve += 1 - @.. broadcast=false k=-k - @.. broadcast=false T[i + 1]=u_temp1 + k + @.. broadcast=false T[i + 1]=u_temp1 - k if (j == j_int + 1) @.. broadcast=false T[i + 1]=0.25(T[i + 1] + 2 * u_temp1 + u_temp2) end @@ -3323,10 +3297,9 @@ function perform_step!(integrator, cache::ImplicitEulerBarycentricExtrapolationC j_int_temp = sequence_factor * subdividing_sequence[index + 1] dt_int_temp = dt / j_int_temp # Stepsize of the ith internal discretisation jacobian2W!(W[Threads.threadid()], integrator.f.mass_matrix, - dt_int_temp, J, false) + dt_int_temp, J, true) @.. broadcast=false u_temp4[Threads.threadid()]=uprev - @.. broadcast=false linsolve_tmps[Threads.threadid()]=dt_int_temp * - fsalfirst + @.. broadcast=false linsolve_tmps[Threads.threadid()]=fsalfirst linsolve = cache.linsolve[Threads.threadid()] @@ -3343,7 +3316,7 @@ function perform_step!(integrator, cache::ImplicitEulerBarycentricExtrapolationC cache.linsolve[Threads.threadid()] = linres.cache @.. broadcast=false k_tmps[Threads.threadid()]=-k_tmps[Threads.threadid()] - @.. broadcast=false u_temp3[Threads.threadid()]=u_temp4[Threads.threadid()] + + @.. broadcast=false u_temp3[Threads.threadid()]=u_temp4[Threads.threadid()] - k_tmps[Threads.threadid()] # Euler starting step @.. broadcast=false diff1[Threads.threadid()]=u_temp3[Threads.threadid()] - u_temp4[Threads.threadid()] @@ -3351,8 +3324,7 @@ function perform_step!(integrator, cache::ImplicitEulerBarycentricExtrapolationC f(k_tmps[Threads.threadid()], cache.u_temp3[Threads.threadid()], p, t + (j - 1) * dt_int_temp) - @.. broadcast=false linsolve_tmps[Threads.threadid()]=dt_int_temp * - k_tmps[Threads.threadid()] + @.. broadcast=false linsolve_tmps[Threads.threadid()]=k_tmps[Threads.threadid()] linsolve = cache.linsolve[Threads.threadid()] @@ -3368,8 +3340,7 @@ function perform_step!(integrator, cache::ImplicitEulerBarycentricExtrapolationC end cache.linsolve[Threads.threadid()] = linres.cache - @.. broadcast=false k_tmps[Threads.threadid()]=-k_tmps[Threads.threadid()] - @.. broadcast=false T[index + 1]=u_temp3[Threads.threadid()] + + @.. broadcast=false T[index + 1]=u_temp3[Threads.threadid()] - k_tmps[Threads.threadid()] # Explicit Midpoint rule if (j == j_int_temp + 1) @.. broadcast=false T[index + 1]=0.25(T[index + 1] + @@ -3414,10 +3385,9 @@ function perform_step!(integrator, cache::ImplicitEulerBarycentricExtrapolationC j_int_temp = sequence_factor * subdividing_sequence[index + 1] dt_int_temp = dt / j_int_temp # Stepsize of the ith internal discretisation jacobian2W!(W[Threads.threadid()], integrator.f.mass_matrix, - dt_int_temp, J, false) + dt_int_temp, J, true) @.. broadcast=false u_temp4[Threads.threadid()]=uprev - @.. broadcast=false linsolve_tmps[Threads.threadid()]=dt_int_temp * - fsalfirst + @.. broadcast=false linsolve_tmps[Threads.threadid()]=fsalfirst linsolve = cache.linsolve[Threads.threadid()] @@ -3433,8 +3403,7 @@ function perform_step!(integrator, cache::ImplicitEulerBarycentricExtrapolationC end cache.linsolve[Threads.threadid()] = linres.cache - @.. broadcast=false k_tmps[Threads.threadid()]=-k_tmps[Threads.threadid()] - @.. broadcast=false u_temp3[Threads.threadid()]=u_temp4[Threads.threadid()] + + @.. broadcast=false u_temp3[Threads.threadid()]=u_temp4[Threads.threadid()] - k_tmps[Threads.threadid()] # Euler starting step @.. broadcast=false diff1[Threads.threadid()]=u_temp3[Threads.threadid()] - u_temp4[Threads.threadid()] @@ -3442,8 +3411,7 @@ function perform_step!(integrator, cache::ImplicitEulerBarycentricExtrapolationC f(k_tmps[Threads.threadid()], cache.u_temp3[Threads.threadid()], p, t + (j - 1) * dt_int_temp) - @.. broadcast=false linsolve_tmps[Threads.threadid()]=dt_int_temp * - k_tmps[Threads.threadid()] + @.. broadcast=false linsolve_tmps[Threads.threadid()]=k_tmps[Threads.threadid()] linsolve = cache.linsolve[Threads.threadid()] @@ -3459,8 +3427,7 @@ function perform_step!(integrator, cache::ImplicitEulerBarycentricExtrapolationC end cache.linsolve[Threads.threadid()] = linres.cache - @.. broadcast=false k_tmps[Threads.threadid()]=-k_tmps[Threads.threadid()] - @.. broadcast=false T[index + 1]=u_temp3[Threads.threadid()] + + @.. broadcast=false T[index + 1]=u_temp3[Threads.threadid()] - k_tmps[Threads.threadid()] # Explicit Midpoint rule if (j == j_int_temp + 1) @.. broadcast=false T[index + 1]=0.25(T[index + 1] + @@ -3548,10 +3515,10 @@ function perform_step!(integrator, cache::ImplicitEulerBarycentricExtrapolationC # Update cache.T j_int = sequence_factor * subdividing_sequence[n_curr + 1] dt_int = dt / j_int # Stepsize of the new internal discretisation - jacobian2W!(W[1], integrator.f.mass_matrix, dt_int, J, false) + jacobian2W!(W[1], integrator.f.mass_matrix, dt_int, J, true) integrator.stats.nw += 1 @.. broadcast=false u_temp2=uprev - @.. broadcast=false linsolve_tmps[1]=dt_int * fsalfirst + @.. broadcast=false linsolve_tmps[1]=fsalfirst linsolve = cache.linsolve[1] @@ -3570,7 +3537,7 @@ function perform_step!(integrator, cache::ImplicitEulerBarycentricExtrapolationC for j in 2:(j_int + 1) f(k, cache.u_temp1, p, t + (j - 1) * dt_int) OrdinaryDiffEqCore.increment_nf!(integrator.stats, 1) - @.. broadcast=false linsolve_tmps[1]=dt_int * k + @.. broadcast=false linsolve_tmps[1]=k linsolve = cache.linsolve[1] linres = dolinsolve(integrator, linsolve; b = _vec(linsolve_tmps[1]), @@ -3578,8 +3545,7 @@ function perform_step!(integrator, cache::ImplicitEulerBarycentricExtrapolationC cache.linsolve[1] = linres.cache integrator.stats.nsolve += 1 - @.. broadcast=false k=-k - @.. broadcast=false T[n_curr + 1]=u_temp1 + k # Explicit Midpoint rule + @.. broadcast=false T[n_curr + 1]=u_temp1 - k # Explicit Midpoint rule if (j == j_int + 1) @.. broadcast=false T[n_curr + 1]=0.25(T[n_curr + 1] + 2 * u_temp1 + u_temp2)