JuliaDiffEq Vern9() solver taking too long

General Usage
14

@ChrisRackauckas The script is same as I had posted originally, just added a few more solvers. Here is it again in its entirety

using DifferentialEquations, ODEInterfaceDiffEq,  BenchmarkTools, Libdl, Statistics, PyCall, LinearAlgebra, ParameterizedFunctions, Plots,Printf,Sundials,LSODA;
gr();

###### @show rtol = 1.0./10.0.^(6:1:13);
@show rtol = 1.0./10.0.^(11);
@show atol = rtol*1e-3;



function cr3bpeom(dX, X, p, t)
    @inbounds begin
        r1 = sqrt((X[1] + mu)^2 + X[2]^2)
        r2 = sqrt((X[1] - 1 + mu)^2 + X[2]^2)

        dX[1:3] .= @view X[4:6]
        dX[4] = X[1] + 2*X[5] - (1 - mu)*(X[1] + mu)/r1^3 - mu*(X[1] - 1 + mu)/r2^3
        dX[5] = X[2] - 2*X[4] - (1 - mu)*X[2]/r1^3 - mu*X[2]/r2^3
        dX[6] = 0.0
    end
end



function JacobiConstant(odesol)
    
    X = odesol[1:3,:]
    Xdot = odesol[4:6,:]

    ~,n = size(X)
    
    # distances
    R1 = copy(X)
    R2 = copy(X)
    R1[1,:] = R1[1,:] .+ mu
    R2[1,:] = R2[1,:] .- 1 .+ mu
    r1 = [norm(R1[:,i]) for i = 1:n]
    r2 = [norm(R2[:,i]) for i = 1:n]
    
    # radius and velocity nagnitudes
    R =  [sqrt(sum((X.^2)[:,i])) for i = 1:n]
    V = [sqrt(sum((Xdot.^2)[:,i])) for i = 1:n]

    # pseudo-potential
    Omega = 1/2*(R.^2) .+ (1-mu)./r1 .+ mu./r2
    
    # Jacobi's Constant
    C = (V.^2)/2 .- Omega
    
    return C
end;




const Runs = 1000;

# The 2nd element of each solver sets the dense=true option
# The 3rd element is the name used in results

struct Solver
    solver
    dense::Bool
    name::String
end

ODESolvers = [
                Solver(dopri5(),false,"DOPRI5"),                              # Hairer's DOPRI5
                Solver(DP5(),true,"DP5"),                                     # Julia's implementation of DOPRI5
                Solver(Tsit5(),true,"TSIT5"),                                 # Tsitouras 5/4 RK
                Solver(dop853(),false,"DOP853"),                              # Hairer's DOP853
                Solver(DP8(),true,"DP8"),                                     # Julia's implementation of DOP853
                #Solver(ARKODE(Sundials.Explicit(), order = 8),true,"ARKODE8"),# Sundial's explicit RK
                Solver(Vern6(),true,"Vern6"),                                 # Verner's 7/6 RK
                #Solver(Vern7(),true,"Vern7"),                                 # Verner's 7/6 RK
                #Solver(Vern8(),true,"Vern8"),                                 # Verner's 7/6 RK    
                Solver(Vern9(),true,"Vern9"),                                 # Verner's 9/8 RK
                Solver(ddeabm(),false,"DDEABM"),                              # Shampine's DDEABM
                #Solver(VCABM(),true,"JABM"),      # forcing dt_min           # Julia's implementation of Shampine's DDEABM
                #Solver(CVODE_Adams(max_order = 12),true,"SunAM12") # Too much error!   # Sundial's Adams-Moulton solver
            ]




# test result Struct
struct ODETestResult
    FinalStateErr
    IOMErr
    MeanTime
    OdeSol
end

# Integrate ODE
function FireODE(X0, ODE, tspan, atol, rtol, method, prob, DenseOn)
 if DenseOn == true
        solve(prob,method,abstol=atol, reltol=rtol,timeseries_errors = false, dense_errors=false, dense=false)  
    else 
        solve(prob,method,abstol=atol, reltol=rtol,timeseries_errors = false, dense_errors=false)  
    end
end
            
# Main ODE Integration Method Test function
function FireODEIntTest(X0, ODE, tspan, OdeMethod,  atol, rtol, FinalState, IOM, IOMFunc::Function, DenseOn)
    
    n = length(FinalState)
    # Define ODE problem
    prob = ODEProblem(ODE, X0, tspan);                
        
    # run solver bunch of times for benchmarking
    sol = 0
    rtime = @elapsed for i in 1:Runs
        sol = FireODE(X0, ODE, tspan, atol, rtol, OdeMethod, prob, DenseOn)
    end
    
    # final state error metric
    serr = norm(sol[1:n,end] - FinalState[1:n])
    
    # compute IOM variations
    IOMsol = IOMFunc(sol)
    ierr = abs.(IOMsol .- IOM)
    
    ODETestResult(serr, ierr, rtime, sol)
end
            
# mass-ratio
const mu = 0.012277471::Float64

# Initial conditions
R0 = [0.994::Float64, 0.0, 0.0]
V0 = [0.0, parse(Float64,"-2.00158510637908252240537862224"), 0.0]

# Time period
Period = parse(Float64,"17.0652165601579625588917206249")

# propagate for bunch of orbits
t0 = 0.0::Float64;
tf = 4.0*Period;

tspan = (t0, tf);


X0 = [R0; V0]

op1 = ODEProblem(cr3bpeom,X0,tspan);

nsol = length(ODESolvers)
ntol = length(rtol)

SolResults3B = Matrix(undef, nsol, ntol)

for itr in 1:ntol
    for ctr in 1:nsol
        print(string("Running: ",ODESolvers[ctr].name, "\n"))
        SolResults3B[ctr,itr] = FireODEIntTest(X0, cr3bpeom, tspan, ODESolvers[ctr].solver,  atol[itr], rtol[itr], [R0; V0], 0, JacobiConstant, ODESolvers[ctr].dense);
    end
    print(string("tol: ", rtol[itr],"\n"))
end

print("Computational time vs errors for the tightest tolerance\n=======================================================\n")
for ctr in 1:nsol
    @Printf.printf("%7s: %.5e, %.5e\n",ODESolvers[ctr].name,SolResults3B[ctr,end].MeanTime,SolResults3B[ctr,end].FinalStateErr)
end