Skip to content

Commit

Permalink
fix out of place on v0.6
Browse files Browse the repository at this point in the history
  • Loading branch information
ChrisRackauckas committed Apr 4, 2017
1 parent 371331a commit 48a0f68
Show file tree
Hide file tree
Showing 8 changed files with 35 additions and 36 deletions.
4 changes: 2 additions & 2 deletions src/integrators/explicit_rk_integrator.jl
Original file line number Diff line number Diff line change
Expand Up @@ -35,7 +35,7 @@ end
for i = 2:stages
uEEst = @muladd uEEst + αEEst[i]*kk[i]
end
integrator.EEst = integrator.opts.internalnorm( dt*(utilde-uEEst)/@muladd(integrator.opts.abstol+max(abs(uprev),abs(u))*integrator.opts.reltol))
integrator.EEst = integrator.opts.internalnorm( dt*(utilde-uEEst)./@muladd(integrator.opts.abstol+max(abs.(uprev),abs.(u))*integrator.opts.reltol))
end
if isfsal(integrator.alg.tableau)
integrator.fsallast = kk[end]
Expand Down Expand Up @@ -114,7 +114,7 @@ end
end
end
for i in uidx
atmp[i] = (dt*(utilde[i]-uEEst[i])/@muladd(integrator.opts.abstol+max(abs(uprev[i]),abs(u[i]))*integrator.opts.reltol))
atmp[i] = (dt*(utilde[i]-uEEst[i])./@muladd(integrator.opts.abstol+max(abs(uprev[i]),abs(u[i]))*integrator.opts.reltol))
end
integrator.EEst = integrator.opts.internalnorm(atmp)
end
Expand Down
8 changes: 4 additions & 4 deletions src/integrators/feagin_rk_integrators.jl
Original file line number Diff line number Diff line change
Expand Up @@ -27,7 +27,7 @@ end
k17 = f(@muladd(@muladd(t + c16*dt)),@muladd(uprev + dt*(a1600*k1 + a1601*k2 + a1602*k3 + a1604*k5 + a1605*k6 + a1606*k7 + a1607*k8 + a1608*k9 + a1609*k10 + a1610*k11 + a1611*k12 + a1612*k13 + a1613*k14 + a1614*k15 + a1615*k16)))
u = @muladd uprev + dt*(b1*k1 + b2*k2 + b3*k3 + b5*k5 + b7*k7 + b9*k9 + b10*k10 + b11*k11 + b12*k12 + b13*k13 + b14*k14 + b15*k15 + b16*k16 + b17*k17)
if integrator.opts.adaptive
integrator.EEst = integrator.opts.internalnorm((dt*(k2 - k16) * adaptiveConst)./@muladd(integrator.opts.abstol+max(abs(uprev),abs(u))*integrator.opts.reltol))
integrator.EEst = integrator.opts.internalnorm((dt*(k2 - k16) * adaptiveConst)./@muladd(integrator.opts.abstol+max(abs.(uprev),abs.(u))*integrator.opts.reltol))
end
k = f(t+dt,u) # For the interpolation, needs k at the updated point
integrator.fsallast = k
Expand Down Expand Up @@ -170,7 +170,7 @@ end
k = f(t+dt,u)
integrator.fsallast = k
if integrator.opts.adaptive
integrator.EEst = integrator.opts.internalnorm((dt*(k2 - k24) * adaptiveConst)./@muladd(integrator.opts.abstol+max(abs(uprev),abs(u))*integrator.opts.reltol))
integrator.EEst = integrator.opts.internalnorm((dt*(k2 - k24) * adaptiveConst)./@muladd(integrator.opts.abstol+max(abs.(uprev),abs.(u))*integrator.opts.reltol))
end
integrator.k[1] = integrator.fsalfirst
integrator.k[2] = integrator.fsallast
Expand Down Expand Up @@ -295,7 +295,7 @@ end
end
if integrator.opts.adaptive
for i in uidx
atmp[i] = (dt*(k2[i] - k24[i]) * adaptiveConst)/@muladd(integrator.opts.abstol+max(abs(uprev[i]),abs(u[i]))*integrator.opts.reltol)
atmp[i] = (dt*(k2[i] - k24[i]) * adaptiveConst)./@muladd(integrator.opts.abstol+max(abs(uprev[i]),abs(u[i]))*integrator.opts.reltol)
end
integrator.EEst = integrator.opts.internalnorm(atmp)
end
Expand Down Expand Up @@ -350,7 +350,7 @@ end
k35 = f(@muladd(t + c34*dt),@muladd(uprev + dt*(a3400*k1 + a3401*k2 + a3402*k3 + a3404*k5 + a3406*k7 + a3407*k8 + a3409*k10 + a3410*k11 + a3411*k12 + a3412*k13 + a3413*k14 + a3414*k15 + a3415*k16 + a3416*k17 + a3417*k18 + a3418*k19 + a3419*k20 + a3420*k21 + a3421*k22 + a3422*k23 + a3423*k24 + a3424*k25 + a3425*k26 + a3426*k27 + a3427*k28 + a3428*k29 + a3429*k30 + a3430*k31 + a3431*k32 + a3432*k33 + a3433*k34)))
u = @muladd uprev + dt*(b1*k1 + b2*k2 + b3*k3 + b5*k5 + b7*k7 + b8*k8 + b10*k10 + b11*k11 + b12*k12 + b14*k14 + b15*k15 + b16*k16 + b18*k18 + b19*k19 + b20*k20 + b21*k21 + b22*k22 + b23*k23 + b24*k24 + b25*k25 + b26*k26 + b27*k27 + b28*k28 + b29*k29 + b30*k30 + b31*k31 + b32*k32 + b33*k33 + b34*k34 + b35*k35)
if integrator.opts.adaptive
integrator.EEst = integrator.opts.internalnorm((dt*(k2 - k34) * adaptiveConst)./@muladd(integrator.opts.abstol+max(abs(uprev),abs(u))*integrator.opts.reltol))
integrator.EEst = integrator.opts.internalnorm((dt*(k2 - k34) * adaptiveConst)./@muladd(integrator.opts.abstol+max(abs.(uprev),abs.(u))*integrator.opts.reltol))
end
k = f(t+dt,u) # For the interpolation, needs k at the updated point
integrator.fsallast = k
Expand Down
20 changes: 10 additions & 10 deletions src/integrators/high_order_rk_integrators.jl
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,7 @@ end
u = @muladd uprev + dt*(k1*b1+k4*b4+k5*b5+k6*b6+k7*b7+k8*b8+k9*b9)
if integrator.opts.adaptive
utilde = @muladd uprev + dt*(k1*bhat1+k4*bhat4+k5*bhat5+k6*bhat6+k7*bhat7+k8*bhat8+k10*bhat10)
integrator.EEst = integrator.opts.internalnorm( ((utilde-u)/@muladd(integrator.opts.abstol+max(abs(uprev),abs(u))*integrator.opts.reltol)))
integrator.EEst = integrator.opts.internalnorm( ((utilde-u)./@muladd(integrator.opts.abstol+max(abs.(uprev),abs.(u))*integrator.opts.reltol)))
end
k = f(t+dt,u) # For the interpolation, needs k at the updated point
integrator.fsallast = k
Expand Down Expand Up @@ -90,7 +90,7 @@ end
if integrator.opts.adaptive
for i in uidx
utilde[i] = uprev[i] + dt*(k1[i]*bhat1+k4[i]*bhat4+k5[i]*bhat5+k6[i]*bhat6+k7[i]*bhat7+k8[i]*bhat8+k10[i]*bhat10)
atmp[i] = ((utilde[i]-u[i])/(integrator.opts.abstol+max(abs(uprev[i]),abs(u[i]))*integrator.opts.reltol))
atmp[i] = ((utilde[i]-u[i])./(integrator.opts.abstol+max(abs(uprev[i]),abs(u[i]))*integrator.opts.reltol))
end
integrator.EEst = integrator.opts.internalnorm(atmp)
end
Expand Down Expand Up @@ -124,13 +124,13 @@ end
update = dt*kupdate
u = uprev + update
if integrator.opts.adaptive
err5 = integrator.opts.internalnorm(dt*(@muladd(k1*er1 + k6*er6 + k7*er7 + k8*er8 + k9*er9 + k10*er10 + k11*er11 + k12*er12))/@muladd(integrator.opts.abstol+max(abs(uprev),abs(u))*integrator.opts.reltol)) # Order 5
err3 = integrator.opts.internalnorm(@muladd(update - dt*(bhh1*k1 + bhh2*k9 + bhh3*k12))/@muladd(integrator.opts.abstol+max(abs(uprev),abs(u))*integrator.opts.reltol)) # Order 3
err5 = integrator.opts.internalnorm(dt*(@muladd(k1*er1 + k6*er6 + k7*er7 + k8*er8 + k9*er9 + k10*er10 + k11*er11 + k12*er12))./@muladd(integrator.opts.abstol+max(abs.(uprev),abs.(u))*integrator.opts.reltol)) # Order 5
err3 = integrator.opts.internalnorm(@muladd(update - dt*(bhh1*k1 + bhh2*k9 + bhh3*k12))./@muladd(integrator.opts.abstol+max(abs.(uprev),abs.(u))*integrator.opts.reltol)) # Order 3
err52 = err5*err5
if err5 0 && err3 0
integrator.EEst = zero(typeof(integrator.EEst))
else
integrator.EEst = err52/sqrt(err52 + 0.01*err3*err3)
integrator.EEst = err52./sqrt(err52 + 0.01*err3*err3)
end
end
k13 = f(t+dt,u)
Expand Down Expand Up @@ -220,16 +220,16 @@ end
end
if integrator.opts.adaptive
for i in uidx
atmp[i] = (dt*(@muladd(k1[i]*er1 + k6[i]*er6 + k7[i]*er7 + k8[i]*er8 + k9[i]*er9 + k10[i]*er10 + k11[i]*er11 + k12[i]*er12))/@muladd(integrator.opts.abstol+max(abs(uprev[i]),abs(u[i]))*integrator.opts.reltol))
atmp2[i]= (@muladd(update[i] - dt*(bhh1*k1[i] + bhh2*k9[i] + bhh3*k12[i]))/@muladd(integrator.opts.abstol+max(abs(uprev[i]),abs(u[i]))*integrator.opts.reltol))
atmp[i] = (dt*(@muladd(k1[i]*er1 + k6[i]*er6 + k7[i]*er7 + k8[i]*er8 + k9[i]*er9 + k10[i]*er10 + k11[i]*er11 + k12[i]*er12))./@muladd(integrator.opts.abstol+max(abs(uprev[i]),abs(u[i]))*integrator.opts.reltol))
atmp2[i]= (@muladd(update[i] - dt*(bhh1*k1[i] + bhh2*k9[i] + bhh3*k12[i]))./@muladd(integrator.opts.abstol+max(abs(uprev[i]),abs(u[i]))*integrator.opts.reltol))
end
err5 = integrator.opts.internalnorm(atmp) # Order 5
err3 = integrator.opts.internalnorm(atmp2) # Order 3
err52 = err5*err5
if err5 0 && err3 0
integrator.EEst = zero(typeof(integrator.EEst))
else
integrator.EEst = err52/sqrt(err52 + 0.01*err3*err3)
integrator.EEst = err52./sqrt(err52 + 0.01*err3*err3)
end
end
f(t+dt,u,k13)
Expand Down Expand Up @@ -287,7 +287,7 @@ end
update = dt*(@muladd(b1*k1+b6*k6+b7*k7+b8*k8+b9*k9+b10*k10+b11*k11+b12*k12))
u = uprev + update
if integrator.opts.adaptive
integrator.EEst = integrator.opts.internalnorm(@muladd(update - dt*(k1*bhat1 + k6*bhat6 + k7*bhat7 + k8*bhat8 + k9*bhat9 + k10*bhat10 + k13*bhat13))/@muladd(integrator.opts.abstol+max(abs(uprev),abs(u))*integrator.opts.reltol))
integrator.EEst = integrator.opts.internalnorm(@muladd(update - dt*(k1*bhat1 + k6*bhat6 + k7*bhat7 + k8*bhat8 + k9*bhat9 + k10*bhat10 + k13*bhat13))./@muladd(integrator.opts.abstol+max(abs.(uprev),abs.(u))*integrator.opts.reltol))
end
k = f(t+dt,u)
integrator.fsallast = k
Expand Down Expand Up @@ -368,7 +368,7 @@ end
end
if integrator.opts.adaptive
for i in uidx
atmp[i] = (@muladd(update[i] - dt*(k1[i]*bhat1 + k6[i]*bhat6 + k7[i]*bhat7 + k8[i]*bhat8 + k9[i]*bhat9 + k10[i]*bhat10 + k13[i]*bhat13))/@muladd(integrator.opts.abstol+max(abs(uprev[i]),abs(u[i]))*integrator.opts.reltol))
atmp[i] = (@muladd(update[i] - dt*(k1[i]*bhat1 + k6[i]*bhat6 + k7[i]*bhat7 + k8[i]*bhat8 + k9[i]*bhat9 + k10[i]*bhat10 + k13[i]*bhat13))./@muladd(integrator.opts.abstol+max(abs(uprev[i]),abs(u[i]))*integrator.opts.reltol))
end
integrator.EEst = integrator.opts.internalnorm(atmp)
end
Expand Down
14 changes: 7 additions & 7 deletions src/integrators/low_order_rk_integrators.jl
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,7 @@ end
k4 = f(t+dt,u); integrator.fsallast = k4
if integrator.opts.adaptive
utilde = @muladd uprev + dt*(b1*k1 + b2*k2 + b3*k3 + b4*k4)
integrator.EEst = integrator.opts.internalnorm( ((utilde-u)/@muladd(integrator.opts.abstol+max(abs(uprev),abs(u))*integrator.opts.reltol)))
integrator.EEst = integrator.opts.internalnorm( ((utilde-u)./@muladd(integrator.opts.abstol+max(abs.(uprev),abs.(u))*integrator.opts.reltol)))
end
integrator.k[1] = integrator.fsalfirst
integrator.k[2] = integrator.fsallast
Expand Down Expand Up @@ -56,7 +56,7 @@ end
if integrator.opts.adaptive
for i in uidx
utilde[i] = @muladd uprev[i] + dt*(b1*k1[i] + b2*k2[i] + b3*k3[i] + b4*k4[i])
atmp[i] = ((utilde[i]-u[i])/@muladd(integrator.opts.abstol+max(abs(uprev[i]),abs(u[i]))*integrator.opts.reltol))
atmp[i] = ((utilde[i]-u[i])./@muladd(integrator.opts.abstol+max(abs(uprev[i]),abs(u[i]))*integrator.opts.reltol))
end
integrator.EEst = integrator.opts.internalnorm(atmp)
end
Expand Down Expand Up @@ -85,8 +85,8 @@ end
if integrator.opts.adaptive
uhat = dt*@muladd(bhat1*k1 + bhat3*k3 + bhat4*k4 + bhat5*k5 + bhat6*k6)
utilde = @muladd uprev + dt*(btilde1*k1 + btilde2*k2 + btilde3*k3 + btilde4*k4 + btilde5*k5 + btilde6*k6 + btilde7*k7 + btilde8*k8)
EEst1 = integrator.opts.internalnorm( sum(((uhat)./@muladd(integrator.opts.abstol+max(abs(uprev),abs(u))*integrator.opts.reltol))))
EEst2 = integrator.opts.internalnorm( sum(((utilde-u)./@muladd(integrator.opts.abstol+max(abs(uprev),abs(u))*integrator.opts.reltol))))
EEst1 = integrator.opts.internalnorm( sum(((uhat)./@muladd(integrator.opts.abstol+max(abs.(uprev),abs.(u))*integrator.opts.reltol))))
EEst2 = integrator.opts.internalnorm( sum(((utilde-u)./@muladd(integrator.opts.abstol+max(abs.(uprev),abs.(u))*integrator.opts.reltol))))
integrator.EEst = max(EEst1,EEst2)
end
integrator.k[1]=k1; integrator.k[2]=k2; integrator.k[3]=k3;integrator.k[4]=k4;integrator.k[5]=k5;integrator.k[6]=k6;integrator.k[7]=k7;integrator.k[8]=k8
Expand Down Expand Up @@ -173,7 +173,7 @@ end
integrator.fsallast = f(t+dt,u); k7 = integrator.fsallast
if integrator.opts.adaptive
utilde = @muladd uprev + dt*(b1*k1 + b2*k2 + b3*k3 + b4*k4 + b5*k5 + b6*k6 + b7*k7)
integrator.EEst = integrator.opts.internalnorm(((utilde-u)/@muladd(integrator.opts.abstol+max(abs(uprev),abs(u))*integrator.opts.reltol)))
integrator.EEst = integrator.opts.internalnorm(((utilde-u)./@muladd(integrator.opts.abstol+max(abs.(uprev),abs.(u))*integrator.opts.reltol)))
end
integrator.k[1] = k1
integrator.k[2] = k2
Expand Down Expand Up @@ -262,7 +262,7 @@ end
integrator.fsallast = f(t+dt,u); k7 = integrator.fsallast
if integrator.opts.adaptive
utilde = uprev + dt*(b1*k1 + b3*k3 + b4*k4 + b5*k5 + b6*k6 + b7*k7)
integrator.EEst = integrator.opts.internalnorm( ((utilde-u)/@muladd(integrator.opts.abstol+max(abs(uprev),abs(u))*integrator.opts.reltol)))
integrator.EEst = integrator.opts.internalnorm( ((utilde-u)./@muladd(integrator.opts.abstol+max(abs.(uprev),abs.(u))*integrator.opts.reltol)))
end
integrator.k[1] = update
bspl = k1 - update
Expand Down Expand Up @@ -314,7 +314,7 @@ end
if integrator.opts.adaptive
for i in uidx
utilde[i] = @muladd uprev[i] + dt*(b1*k1[i] + b3*k3[i] + b4*k4[i] + b5*k5[i] + b6*k6[i] + b7*k7[i])
atmp[i] = ((utilde[i]-u[i])/@muladd(integrator.opts.abstol+max(abs(uprev[i]),abs(u[i]))*integrator.opts.reltol))
atmp[i] = ((utilde[i]-u[i])./@muladd(integrator.opts.abstol+max(abs(uprev[i]),abs(u[i]))*integrator.opts.reltol))
end
integrator.EEst = integrator.opts.internalnorm(atmp)
end
Expand Down
6 changes: 3 additions & 3 deletions src/integrators/rosenbrock_integrators.jl
Original file line number Diff line number Diff line change
Expand Up @@ -143,7 +143,7 @@ end
end
if integrator.opts.adaptive
for i in uidx
tmp[i] = (dt*(k₁[i] - 2k₂[i] + k₃[i])/6)/@muladd(integrator.opts.abstol+max(abs(uprev[i]),abs(u[i]))*integrator.opts.reltol)
tmp[i] = (dt*(k₁[i] - 2k₂[i] + k₃[i])/6)./@muladd(integrator.opts.abstol+max(abs(uprev[i]),abs(u[i]))*integrator.opts.reltol)
end
integrator.EEst = integrator.opts.internalnorm(tmp)
end
Expand Down Expand Up @@ -176,7 +176,7 @@ end
if integrator.opts.adaptive
integrator.fsallast = f(t+dt,u)
k₃ = W\@muladd(integrator.fsallast - c₃₂*(k₂-f₁)-2(k₁-integrator.fsalfirst)+dt*dT)
integrator.EEst = integrator.opts.internalnorm((dt*(k₁ - 2k₂ + k₃)/6)./@muladd(integrator.opts.abstol+max(abs(uprev),abs(u))*integrator.opts.reltol))
integrator.EEst = integrator.opts.internalnorm((dt*(k₁ - 2k₂ + k₃)/6)./@muladd(integrator.opts.abstol+max(abs.(uprev),abs.(u))*integrator.opts.reltol))
end
integrator.k[1] = k₁
integrator.k[2] = k₂
Expand Down Expand Up @@ -211,7 +211,7 @@ end
k₃ = W\@muladd(integrator.fsallast - c₃₂*(k₂-f₁)-2(k₁-integrator.fsalfirst)+dt*dT)
u = uprev + dt*(k₁ + 4k₂ + k₃)/6
if integrator.opts.adaptive
integrator.EEst = integrator.opts.internalnorm((dt*(k₁ - 2k₂ + k₃)/6)./@muladd(integrator.opts.abstol+max(abs(uprev),abs(u))*integrator.opts.reltol))
integrator.EEst = integrator.opts.internalnorm((dt*(k₁ - 2k₂ + k₃)/6)./@muladd(integrator.opts.abstol+max(abs.(uprev),abs.(u))*integrator.opts.reltol))
end
integrator.k[1] = k₁
integrator.k[2] = k₂
Expand Down
2 changes: 1 addition & 1 deletion src/integrators/threaded_rk_integrators.jl
Original file line number Diff line number Diff line change
Expand Up @@ -80,6 +80,6 @@ end
@noinline function dp5threaded_adaptiveloop(dt,utilde,uprev,b1,k1,b3,k3,b4,k4,b5,k5,b6,k6,b7,k7,atmp,u,abstol,reltol,uidx)
Threads.@threads for i in uidx
utilde[i] = @muladd uprev[i] + dt*(b1*k1[i] + b3*k3[i] + b4*k4[i] + b5*k5[i] + b6*k6[i] + b7*k7[i])
atmp[i] = ((utilde[i]-u[i])/@muladd(abstol+max(abs(uprev[i]),abs(u[i]))*reltol))
atmp[i] = ((utilde[i]-u[i])./@muladd(abstol+max(abs(uprev[i]),abs(u[i]))*reltol))
end
end
16 changes: 8 additions & 8 deletions src/integrators/verner_rk_integrators.jl
Original file line number Diff line number Diff line change
Expand Up @@ -20,7 +20,7 @@ end
integrator.fsallast = f(t+dt,u); k9 = integrator.fsallast
if integrator.opts.adaptive
utilde = @muladd uprev + dt*(b1*k1 + b4*k4 + b5*k5 + b6*k6 + b7*k7 + b8*k8 + b9*k9)
integrator.EEst = integrator.opts.internalnorm( ((utilde-u)/@muladd(integrator.opts.abstol+max(abs(uprev),abs(u))*integrator.opts.reltol)))
integrator.EEst = integrator.opts.internalnorm( ((utilde-u)./@muladd(integrator.opts.abstol+max(abs.(uprev),abs.(u))*integrator.opts.reltol)))
end
integrator.k[1]=k1; integrator.k[2]=k2;
integrator.k[3]=k3; integrator.k[4]=k4;
Expand Down Expand Up @@ -82,7 +82,7 @@ end
if integrator.opts.adaptive
for i in uidx
utilde[i] = @muladd uprev[i] + dt*(b1*k1[i] + b4*k4[i] + b5*k5[i] + b6*k6[i] + b7*k7[i] + b8*k8[i] + b9*k9[i])
atmp[i] = ((utilde[i]-u[i])/@muladd(integrator.opts.abstol+max(abs.(uprev[i]),abs.(u[i]))*integrator.opts.reltol))
atmp[i] = ((utilde[i]-u[i])./@muladd(integrator.opts.abstol+max(abs(uprev[i]),abs(u[i]))*integrator.opts.reltol))
end
integrator.EEst = integrator.opts.internalnorm(atmp)
end
Expand Down Expand Up @@ -112,7 +112,7 @@ end
update = dt*(@muladd(k1*b1 + k4*b4 + k5*b5 + k6*b6 + k7*b7 + k8*b8 + k9*b9))
u = uprev + update
if integrator.opts.adaptive
integrator.EEst = integrator.opts.internalnorm( ((update - dt*(@muladd(bhat1*k1 + bhat4*k4 + bhat5*k5 + bhat6*k6 + bhat7*k7 + bhat10*k10)))/@muladd(integrator.opts.abstol+max(abs.(uprev),abs.(u))*integrator.opts.reltol)))
integrator.EEst = integrator.opts.internalnorm( ((update - dt*(@muladd(bhat1*k1 + bhat4*k4 + bhat5*k5 + bhat6*k6 + bhat7*k7 + bhat10*k10)))./@muladd(integrator.opts.abstol+max(abs.(uprev),abs.(u))*integrator.opts.reltol)))
end
integrator.k[1]=k1; integrator.k[2]=k2;
integrator.k[3]=k3; integrator.k[4]=k4;
Expand Down Expand Up @@ -179,7 +179,7 @@ end
end
if integrator.opts.adaptive
for i in uidx
atmp[i] = (@muladd(update[i] - dt*(bhat1*k1[i] + bhat4*k4[i] + bhat5*k5[i] + bhat6*k6[i] + bhat7*k7[i] + bhat10*k10[i]))/@muladd(integrator.opts.abstol+max(abs.(uprev[i]),abs.(u[i]))*integrator.opts.reltol))
atmp[i] = (@muladd(update[i] - dt*(bhat1*k1[i] + bhat4*k4[i] + bhat5*k5[i] + bhat6*k6[i] + bhat7*k7[i] + bhat10*k10[i]))./@muladd(integrator.opts.abstol+max(abs(uprev[i]),abs(u[i]))*integrator.opts.reltol))
end
integrator.EEst = integrator.opts.internalnorm(atmp)
end
Expand Down Expand Up @@ -212,7 +212,7 @@ end
update = dt*(@muladd(k1*b1 + k6*b6 + k7*b7 + k8*b8 + k9*b9 + k10*b10 + k11*b11 + k12*b12))
u = uprev + update
if integrator.opts.adaptive
integrator.EEst = integrator.opts.internalnorm( (@muladd(update - dt*(bhat1*k1 + bhat6*k6 + bhat7*k7 + bhat8*k8 + bhat9*k9 + bhat10*k10 + bhat13*k13))/@muladd(integrator.opts.abstol+max(abs.(uprev),abs.(u))*integrator.opts.reltol)))
integrator.EEst = integrator.opts.internalnorm( (@muladd(update - dt*(bhat1*k1 + bhat6*k6 + bhat7*k7 + bhat8*k8 + bhat9*k9 + bhat10*k10 + bhat13*k13))./@muladd(integrator.opts.abstol+max(abs.(uprev),abs.(u))*integrator.opts.reltol)))
end
integrator.k[1]=k1; integrator.k[2]=k2;
integrator.k[3]=k3; integrator.k[4]=k4;
Expand Down Expand Up @@ -293,7 +293,7 @@ end
end
if integrator.opts.adaptive
for i in uidx
atmp[i] = (@muladd(update[i] - dt*(bhat1*k1[i] + bhat6*k6[i] + bhat7*k7[i] + bhat8*k8[i] + bhat9*k9[i] + bhat10*k10[i] + bhat13*k13[i]))/@muladd(integrator.opts.abstol+max(abs.(uprev[i]),abs.(u[i]))*integrator.opts.reltol))
atmp[i] = (@muladd(update[i] - dt*(bhat1*k1[i] + bhat6*k6[i] + bhat7*k7[i] + bhat8*k8[i] + bhat9*k9[i] + bhat10*k10[i] + bhat13*k13[i]))./@muladd(integrator.opts.abstol+max(abs(uprev[i]),abs(u[i]))*integrator.opts.reltol))
end
integrator.EEst = integrator.opts.internalnorm(atmp)
end
Expand Down Expand Up @@ -329,7 +329,7 @@ end
update = dt*(@muladd(k1*b1+k8*b8+k9*b9+k10*b10+k11*b11+k12*b12+k13*b13+k14*b14+k15*b15))
u = uprev + update
if integrator.opts.adaptive
integrator.EEst = integrator.opts.internalnorm(@muladd(update - dt*(k1*bhat1 + k8*bhat8 + k9*bhat9 + k10*bhat10 + k11*bhat11 + k12*bhat12 + k13*bhat13 + k16*bhat16))/@muladd(integrator.opts.abstol+max(abs.(uprev),abs.(u))*integrator.opts.reltol))
integrator.EEst = integrator.opts.internalnorm(@muladd(update - dt*(k1*bhat1 + k8*bhat8 + k9*bhat9 + k10*bhat10 + k11*bhat11 + k12*bhat12 + k13*bhat13 + k16*bhat16))./@muladd(integrator.opts.abstol+max(abs.(uprev),abs.(u))*integrator.opts.reltol))
end
integrator.k[1]=k1; integrator.k[2]=k2;
integrator.k[3]=k3; integrator.k[4]=k4;
Expand Down Expand Up @@ -423,7 +423,7 @@ end
end
if integrator.opts.adaptive
for i in uidx
atmp[i] = (@muladd(update[i] - dt*(k1[i]*bhat1 + k8[i]*bhat8 + k9[i]*bhat9 + k10[i]*bhat10 + k11[i]*bhat11 + k12[i]*bhat12 + k13[i]*bhat13 + k16[i]*bhat16))/@muladd(integrator.opts.abstol+max(abs.(uprev[i]),abs.(u[i]))*integrator.opts.reltol))
atmp[i] = (@muladd(update[i] - dt*(k1[i]*bhat1 + k8[i]*bhat8 + k9[i]*bhat9 + k10[i]*bhat10 + k11[i]*bhat11 + k12[i]*bhat12 + k13[i]*bhat13 + k16[i]*bhat16))./@muladd(integrator.opts.abstol+max(abs(uprev[i]),abs(u[i]))*integrator.opts.reltol))
end
integrator.EEst = integrator.opts.internalnorm(atmp)
end
Expand Down
1 change: 0 additions & 1 deletion test/ode/ode_dense_tests.jl
Original file line number Diff line number Diff line change
Expand Up @@ -224,7 +224,6 @@ sol2 =solve(prob,Vern6(),dt=1//2^(7),dense=true,adaptive=false)
#plot!(sol.t,sol[:])
#scatter!(sol.t,sol[:])

f = (t,u,du) -> (du .= linear_bigα.*u)
prob_ode_bigfloatveclinear = ODEProblem(f,[parse(BigFloat,"0.5")],(0.0,1.0))
prob = prob_ode_bigfloatveclinear

Expand Down

0 comments on commit 48a0f68

Please sign in to comment.