# HIRES Work-Precision Diagrams

##### Chris Rackauckas
using OrdinaryDiffEq, ParameterizedFunctions, Plots, ODE, ODEInterfaceDiffEq, LSODA, DiffEqDevTools, Sundials
using LinearAlgebra

gr() #gr(fmt=:png)

f = @ode_def Hires begin
dy1 = -1.71*y1 + 0.43*y2 + 8.32*y3 + 0.0007
dy2 = 1.71*y1 - 8.75*y2
dy3 = -10.03*y3 + 0.43*y4 + 0.035*y5
dy4 = 8.32*y2 + 1.71*y3 - 1.12*y4
dy5 = -1.745*y5 + 0.43*y6 + 0.43*y7
dy6 = -280.0*y6*y8 + 0.69*y4 + 1.71*y5 -
0.43*y6 + 0.69*y7
dy7 = 280.0*y6*y8 - 1.81*y7
dy8 = -280.0*y6*y8 + 1.81*y7
end

u0 = zeros(8)
u0[1] = 1
u0[8] = 0.0057
prob = ODEProblem(f,u0,(0.0,321.8122))

sol = solve(prob,Rodas5(),abstol=1/10^14,reltol=1/10^14)
test_sol = TestSolution(sol)
abstols = 1.0 ./ 10.0 .^ (4:11)
reltols = 1.0 ./ 10.0 .^ (1:8);

plot(sol)

plot(sol,tspan=(0.0,5.0))


## Omissions

The following were omitted from the tests due to convergence failures. ODE.jl's adaptivity is not able to stabilize its algorithms, while GeometricIntegratorsDiffEq has not upgraded to Julia 1.0. GeometricIntegrators.jl's methods used to be either fail to converge at comparable dts (or on some computers errors due to type conversions).

#sol = solve(prob,ode23s()); println("Total ODE.jl steps: \$(length(sol))")
#using GeometricIntegratorsDiffEq
#try
#catch e
#    println(e)
#end


The stabilized explicit methods are not stable enough to handle this problem well. While they don't diverge, they are really slow.

setups = [
#Dict(:alg=>ROCK2()),
#Dict(:alg=>ROCK4())
#Dict(:alg=>ESERK5())
]

0-element Array{Any,1}


## High Tolerances

This is the speed when you just want the answer.

abstols = 1.0 ./ 10.0 .^ (5:8)
reltols = 1.0 ./ 10.0 .^ (1:4);
setups = [Dict(:alg=>Rosenbrock23()),
Dict(:alg=>Rodas3()),
Dict(:alg=>TRBDF2()),
Dict(:alg=>CVODE_BDF()),
Dict(:alg=>rodas()),
Dict(:alg=>ROS34PW1a()),
Dict(:alg=>lsoda()),
]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;
save_everystep=false,appxsol=test_sol,maxiters=Int(1e5),numruns=10)
plot(wp)

wp = WorkPrecisionSet(prob,abstols,reltols,setups;dense = false,verbose=false,
appxsol=test_sol,maxiters=Int(1e5),error_estimate=:l2)
plot(wp)

wp = WorkPrecisionSet(prob,abstols,reltols,setups;
appxsol=test_sol,maxiters=Int(1e5),error_estimate=:L2)
plot(wp)

setups = [Dict(:alg=>Rosenbrock23()),
Dict(:alg=>Kvaerno3()),
Dict(:alg=>CVODE_BDF()),
Dict(:alg=>KenCarp4()),
Dict(:alg=>TRBDF2()),
Dict(:alg=>KenCarp3()),
# Dict(:alg=>SDIRK2()), # Removed because it's bad
wp = WorkPrecisionSet(prob,abstols,reltols,setups;
save_everystep=false,appxsol=test_sol,maxiters=Int(1e5))
plot(wp)

wp = WorkPrecisionSet(prob,abstols,reltols,setups;dense = false,verbose=false,
appxsol=test_sol,maxiters=Int(1e5),error_estimate=:l2)
plot(wp)

wp = WorkPrecisionSet(prob,abstols,reltols,setups;
appxsol=test_sol,maxiters=Int(1e5),error_estimate=:L2)
plot(wp)

setups = [Dict(:alg=>Rosenbrock23()),
Dict(:alg=>KenCarp5()),
Dict(:alg=>KenCarp4()),
Dict(:alg=>KenCarp3()),
Dict(:alg=>ARKODE(order=5)),
Dict(:alg=>ARKODE()),
Dict(:alg=>ARKODE(order=3))]
names = ["Rosenbrock23" "KenCarp5" "KenCarp4" "KenCarp3" "ARKODE5" "ARKODE4" "ARKODE3"]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;
names=names,save_everystep=false,appxsol=test_sol,maxiters=Int(1e5))
plot(wp)

wp = WorkPrecisionSet(prob,abstols,reltols,setups;dense = false,verbose=false,
appxsol=test_sol,maxiters=Int(1e5),error_estimate=:l2)
plot(wp)

setups = [Dict(:alg=>Rosenbrock23()),
Dict(:alg=>TRBDF2()),
Dict(:alg=>ImplicitEulerExtrapolation()),
#Dict(:alg=>ImplicitDeuflhardExtrapolation()), # Diverges
#Dict(:alg=>ImplicitHairerWannerExtrapolation()), # Diverges
Dict(:alg=>ABDF2()),
Dict(:alg=>QNDF()),
Dict(:alg=>Exprb43()),
Dict(:alg=>Exprb32()),
]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;
save_everystep=false,appxsol=test_sol,maxiters=Int(1e5))
plot(wp)


### Low Tolerances

This is the speed at lower tolerances, measuring what's good when accuracy is needed.

abstols = 1.0 ./ 10.0 .^ (7:13)
reltols = 1.0 ./ 10.0 .^ (4:10)

setups = [Dict(:alg=>GRK4A()),
Dict(:alg=>Rodas4P()),
Dict(:alg=>CVODE_BDF()),
Dict(:alg=>ddebdf()),
Dict(:alg=>Rodas5()),
Dict(:alg=>rodas()),
Dict(:alg=>lsoda()),
]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;
save_everystep=false,appxsol=test_sol,maxiters=Int(1e5))
plot(wp)

wp = WorkPrecisionSet(prob,abstols,reltols,setups;verbose=false,
dense=false,appxsol=test_sol,maxiters=Int(1e5),error_estimate=:l2)
plot(wp)

wp = WorkPrecisionSet(prob,abstols,reltols,setups;
appxsol=test_sol,maxiters=Int(1e5),error_estimate=:L2)
plot(wp)

setups = [
Dict(:alg=>Rodas5()),
Dict(:alg=>Kvaerno5()),
Dict(:alg=>CVODE_BDF()),
Dict(:alg=>KenCarp4()),
Dict(:alg=>Rodas4()),
wp = WorkPrecisionSet(prob,abstols,reltols,setups;
save_everystep=false,appxsol=test_sol,maxiters=Int(1e5))
plot(wp)

wp = WorkPrecisionSet(prob,abstols,reltols,setups;verbose=false,
dense=false,appxsol=test_sol,maxiters=Int(1e5),error_estimate=:l2)
plot(wp)

wp = WorkPrecisionSet(prob,abstols,reltols,setups;
appxsol=test_sol,maxiters=Int(1e5),error_estimate=:L2)
plot(wp)

setups = [Dict(:alg=>Rosenbrock23()),
Dict(:alg=>KenCarp5()),
Dict(:alg=>KenCarp4()),
Dict(:alg=>KenCarp3()),
Dict(:alg=>ARKODE(order=5)),
Dict(:alg=>ARKODE()),
Dict(:alg=>ARKODE(order=3))]
names = ["Rosenbrock23" "KenCarp5" "KenCarp4" "KenCarp3" "ARKODE5" "ARKODE4" "ARKODE3"]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;
names=names,save_everystep=false,appxsol=test_sol,maxiters=Int(1e5))
plot(wp)