# 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


## High Tolerances

This is the speed when you just want the answer.

solve(prob, ddebdf())
solve(prob, rodas())
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=>lsoda()),
Dict(:alg=>ROCK2()),
Dict(:alg=>ROCK4())
]
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)


### 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()),
#Dict(:alg=>ROCK2()),   #Needs more iterations
Dict(:alg=>ROCK4())
]
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)

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