JuliaDiffEq Vern9() solver taking too long

General Usage
16

Note that if you use StaticArrays you also get right down to there:

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

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



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

        dX4 = X[1] + 2*X[5] - (1 - mu)*(X[1] + mu)/r1^3 - mu*(X[1] - 1 + mu)/r2^3
        dX5 = X[2] - 2*X[4] - (1 - mu)*X[2]/r1^3 - mu*X[2]/r2^3

        @SVector [X[4],X[5],X[6],dX4,dX5,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

    # Remove compilation run
    sol = FireODE(X0, ODE, tspan, atol, rtol, OdeMethod, prob, DenseOn)

    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 = @SVector [0.994::Float64, 0.0, 0.0]
V0 = @SVector [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

Your old code before:

    DP5: 8.90705e+00, 1.71606e+00
  TSIT5: 9.15936e+00, 7.33951e-02
    DP8: 6.43439e+00, 6.66641e-01
  Vern6: 8.61046e+00, 3.99717e-01
  Vern9: 1.09995e+01, 5.92027e-02

This code:

    DP5: 3.88171e+00, 1.72335e+00
  TSIT5: 5.00398e+00, 6.24401e-02
    DP8: 1.08576e+00, 6.16423e-01
  Vern6: 2.71387e+00, 5.33841e-01
  Vern9: 1.14183e+00, 5.07811e-02
1 Like