JuliaDiffEq Vern9() solver taking too long

General Usage
15

I profiled it and the view in the differential equation was surprisingly unoptimized, so

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] = X[4]
        dX[2] = X[5]
        dX[3] = X[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

sped it up quite a bit. One thing that was happening was that the compilation was inside of the loop, so the ODE solver needs to be run once before the loop if you want to exclude the compile time. And as I expected Julia’s GC is what you’re hitting in the loop, so I get another reduction by doing:

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 = 0.0

    # Compile before timing
    sol = FireODE(X0, ODE, tspan, atol, rtol, OdeMethod, prob, DenseOn)

    for i in 1:Runs
        GC.gc()
        rtime += @elapsed 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

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

After:

    DP5: 4.15515e+00, 1.71606e+00
  TSIT5: 4.82140e+00, 7.33951e-02
    DP8: 1.15518e+00, 6.66641e-01
  Vern6: 3.15056e+00, 3.99717e-01
  Vern9: 1.02262e+00, 5.92027e-02

That amount of time reduction probably puts it in line with Fortran, so that’s where the difference between the two likely were (I can’t find out how to run Flint on my computer: would you be willing to make a Julia wrapper? It would be nice to benchmark it!)

1 Like