From d24d72aae719875b22ad9fac9e864829b5c0c3b5 Mon Sep 17 00:00:00 2001 From: ChrisRackauckas Date: Wed, 3 May 2017 08:16:29 -0700 Subject: [PATCH] max. --- src/integrators/explicit_rk_integrator.jl | 2 +- src/integrators/feagin_rk_integrators.jl | 6 +++--- src/integrators/high_order_rk_integrators.jl | 8 ++++---- src/integrators/low_order_rk_integrators.jl | 10 +++++----- src/integrators/rosenbrock_integrators.jl | 4 ++-- src/integrators/verner_rk_integrators.jl | 8 ++++---- 6 files changed, 19 insertions(+), 19 deletions(-) diff --git a/src/integrators/explicit_rk_integrator.jl b/src/integrators/explicit_rk_integrator.jl index bdf86f2985..6663350e5a 100644 --- a/src/integrators/explicit_rk_integrator.jl +++ b/src/integrators/explicit_rk_integrator.jl @@ -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] diff --git a/src/integrators/feagin_rk_integrators.jl b/src/integrators/feagin_rk_integrators.jl index 2bcd6e280f..af9ee5b9d9 100644 --- a/src/integrators/feagin_rk_integrators.jl +++ b/src/integrators/feagin_rk_integrators.jl @@ -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 @@ -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 @@ -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 diff --git a/src/integrators/high_order_rk_integrators.jl b/src/integrators/high_order_rk_integrators.jl index cc098ca2cc..bf5ae1802f 100644 --- a/src/integrators/high_order_rk_integrators.jl +++ b/src/integrators/high_order_rk_integrators.jl @@ -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 @@ -124,8 +124,8 @@ 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)) @@ -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 diff --git a/src/integrators/low_order_rk_integrators.jl b/src/integrators/low_order_rk_integrators.jl index afd2e51023..b6392de925 100644 --- a/src/integrators/low_order_rk_integrators.jl +++ b/src/integrators/low_order_rk_integrators.jl @@ -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 @@ -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 @@ -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 @@ -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 diff --git a/src/integrators/rosenbrock_integrators.jl b/src/integrators/rosenbrock_integrators.jl index 4f9b5b937e..c60b68af43 100644 --- a/src/integrators/rosenbrock_integrators.jl +++ b/src/integrators/rosenbrock_integrators.jl @@ -179,7 +179,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₂ @@ -218,7 +218,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₂ diff --git a/src/integrators/verner_rk_integrators.jl b/src/integrators/verner_rk_integrators.jl index 3858c4b27a..e5036ee5c0 100644 --- a/src/integrators/verner_rk_integrators.jl +++ b/src/integrators/verner_rk_integrators.jl @@ -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; @@ -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; @@ -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; @@ -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;