# Benchmarking all algorithms in DiffEq/Optimisation code

**URL:** <https://discourse.julialang.org/t/benchmarking-all-algorithms-in-diffeq-optimisation-code/82358>\
**Category:** Performance\
**Tags:** diffeq, benchmark, optimization\
**Created:** [June 7, 2022, 6:02am UTC](https://discourse.julialang.org/t/benchmarking-all-algorithms-in-diffeq-optimisation-code/82358 "2022-06-07T06:02:48Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![DanielVandH](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/danielvandh/32/31134_2.png) [@DanielVandH](https://discourse.julialang.org/u/DanielVandH)\
**Post date:** [June 7, 2022, 6:02am UTC](https://discourse.julialang.org/t/benchmarking-all-algorithms-in-diffeq-optimisation-code/82358/1 "2022-06-07T06:02:48Z")

</div>

If I have some ODE code, say this code from the [https://diffeq.sciml.ai/stable/tutorials/advanced\_ode\_example/:](https://diffeq.sciml.ai/stable/tutorials/advanced_ode_example/:)

```julia
using DifferentialEquations, LinearAlgebra, SparseArrays

const N = 32
const xyd_brusselator = range(0,stop=1,length=N)
brusselator_f(x, y, t) = (((x-0.3)^2 + (y-0.6)^2) <= 0.1^2) * (t >= 1.1) * 5.
limit(a, N) = a == N+1 ? 1 : a == 0 ? N : a
function brusselator_2d_loop(du, u, p, t)
  A, B, alpha, dx = p
  alpha = alpha/dx^2
  @inbounds for I in CartesianIndices((N, N))
    i, j = Tuple(I)
    x, y = xyd_brusselator[I[1]], xyd_brusselator[I[2]]
    ip1, im1, jp1, jm1 = limit(i+1, N), limit(i-1, N), limit(j+1, N), limit(j-1, N)
    du[i,j,1] = alpha*(u[im1,j,1] + u[ip1,j,1] + u[i,jp1,1] + u[i,jm1,1] - 4u[i,j,1]) +
                B + u[i,j,1]^2*u[i,j,2] - (A + 1)*u[i,j,1] + brusselator_f(x, y, t)
    du[i,j,2] = alpha*(u[im1,j,2] + u[ip1,j,2] + u[i,jp1,2] + u[i,jm1,2] - 4u[i,j,2]) +
                A*u[i,j,1] - u[i,j,1]^2*u[i,j,2]
    end
end
p = (3.4, 1., 10., step(xyd_brusselator))

function init_brusselator_2d(xyd)
  N = length(xyd)
  u = zeros(N, N, 2)
  for I in CartesianIndices((N, N))
    x = xyd[I[1]]
    y = xyd[I[2]]
    u[I,1] = 22*(y*(1-y))^(3/2)
    u[I,2] = 27*(x*(1-x))^(3/2)
  end
  u
end
u0 = init_brusselator_2d(xyd_brusselator)
prob_ode_brusselator_2d = ODEProblem(brusselator_2d_loop,u0,(0.,11.5),p)

```

is there an easy way to efficiently benchmark solving this problem with `solve` for all (or most) available algorithms, and say all linear solvers and all preconditioners from `DifferentialEquations.jl` and `LinearSolve.jl`? Or is the way that this is done just simply looping over a hand-written list of algorithms? The problems I’m interested in arise from a discretisation of a 2D finite volume code over a general domain.

Similarly: Is there a similar way to benchmark optimisation algorithms from [GitHub - SciML/Optimization.jl: Local, global, and beyond optimization for scientific machine learning (SciML)](https://github.com/SciML/Optimization.jl) for a given optimisation problem? Or is it again just done by looping over a hand-written list?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [June 7, 2022, 11:05am UTC](https://discourse.julialang.org/t/benchmarking-all-algorithms-in-diffeq-optimisation-code/82358/2 "2022-06-07T11:05:41Z")

</div>

Generally it’s done by looping over a handwritten list. You can do something fancier using `subtypes`:

```julia
julia> println(subtypes(OrdinaryDiffEq.OrdinaryDiffEqAlgorithm))
Any[AB3, AB4, AB5, ABM32, ABM43, ABM54, Anas5, CFRLDDRK64, CalvoSanz4, CandyRoz4, CarpenterKennedy2N54, CayleyEuler, DGLDDRK73_C, DGLDDRK84_C, DGLDDRK84_F, Euler, FunctionMap, HSLDDRK64, IRKN3, IRKN4, KYK2014DGSSPRK_3S2, KYKSSPRK42, KahanLi6, KahanLi8, KuttaPRK2p5, MSRK5, McAte2, McAte3, McAte4, McAte42, McAte5, McAte8, NDBLSRK124, NDBLSRK134, NDBLSRK144, Nystrom4, Nystrom4VelocityIndependent, Nystrom5VelocityIndependent, ORK256, OrdinaryDiffEq.OrdinaryDiffEqAdaptiveAlgorithm, OrdinaryDiffEq.OrdinaryDiffEqCompositeAlgorithm, OrdinaryDiffEq.OrdinaryDiffEqExponentialAlgorithm, OrdinaryDiffEq.OrdinaryDiffEqImplicitAlgorithm, ParsaniKetchesonDeconinck3S105, ParsaniKetchesonDeconinck3S173, ParsaniKetchesonDeconinck3S184, ParsaniKetchesonDeconinck3S205, ParsaniKetchesonDeconinck3S32, ParsaniKetchesonDeconinck3S53, ParsaniKetchesonDeconinck3S82, ParsaniKetchesonDeconinck3S94, PseudoVerletLeapfrog, RK46NL, RKM, RKO65, Ruth3, SHLDDRK52, SHLDDRK64, SHLDDRK_2N, SSPRK104, SSPRK22, SSPRK33, SSPRK53, SSPRK53_2N1, SSPRK53_2N2, SSPRK53_H, SSPRK54, SSPRK63, SSPRK73, SSPRK83, SofSpa10, SymplecticEuler, TSLDDRK74, VelocityVerlet, VerletLeapfrog, Yoshida6]

julia> println(subtypes(OrdinaryDiffEq.OrdinaryDiffEqAdaptiveAlgorithm))
Any[AN5, BS3, BS5, CKLLSRK43_2, CKLLSRK54_3C, CKLLSRK54_3C_3R, CKLLSRK54_3M_3R, CKLLSRK54_3M_4R, CKLLSRK54_3N_3R, CKLLSRK54_3N_4R, CKLLSRK65_4M_4R, CKLLSRK75_4M_5R, CKLLSRK85_4C_3R, CKLLSRK85_4FM_4R, CKLLSRK85_4M_3R, CKLLSRK85_4P_3R, CKLLSRK95_4C, CKLLSRK95_4M, CKLLSRK95_4S, DP5, DP8, DPRKN12, DPRKN6, DPRKN8, ERKN4, ERKN5, ESERK4, ESERK5, ExplicitRK, FRK65, Feagin10, Feagin12, Feagin14, Heun, MagnusAdapt4, Midpoint, OrdinaryDiffEq.OrdinaryDiffEqAdamsVarOrderVarStepAlgorithm, OrdinaryDiffEq.OrdinaryDiffEqAdaptiveExponentialAlgorithm, OrdinaryDiffEq.OrdinaryDiffEqAdaptiveImplicitAlgorithm, OrdinaryDiffEq.OrdinaryDiffEqExtrapolationVarOrderVarStepAlgorithm, OwrenZen3, OwrenZen4, OwrenZen5, PFRK87, RDPK3Sp35, RDPK3Sp49, RDPK3Sp510, RDPK3SpFSAL35, RDPK3SpFSAL49, RDPK3SpFSAL510, RK4, RKC, ROCK2, ROCK4, Ralston, SERK2, SSPRK43, SSPRK432, SSPRK932, SSPRKMSVS32, SSPRKMSVS43, TanYam7, Tsit5, TsitPap8, VCAB3, VCAB4, VCAB5, VCABM3, VCABM4, VCABM5, Vern6, Vern7, Vern8, Vern9]

julia> println(subtypes(OrdinaryDiffEq.OrdinaryDiffEqAdaptiveImplicitAlgorithm))
Any[OrdinaryDiffEq.OrdinaryDiffEqImplicitExtrapolationAlgorithm, OrdinaryDiffEq.OrdinaryDiffEqNewtonAdaptiveAlgorithm, OrdinaryDiffEq.OrdinaryDiffEqRosenbrockAdaptiveAlgorithm]

julia> println(subtypes(OrdinaryDiffEq.OrdinaryDiffEqNewtonAdaptiveAlgorithm))
Any[ABDF2, Cash4, ESDIRK54I8L2SA, FBDF, Hairer4, Hairer42, IRKC, ImplicitEuler, KenCarp3, KenCarp4, KenCarp47, KenCarp5, KenCarp58, Kvaerno3, Kvaerno4, Kvaerno5, QNDF, QNDF1, QNDF2, RadauIIA3, RadauIIA5, SDIRK2, SDIRK22, TRBDF2, Trapezoid]

```

though you’d probably have to weed out a lot that you don’t actually want to be trying all of the time. Similarly you can grab all linear solvers.

[SciMLBenchmarks](https://github.com/SciML/SciMLBenchmarks.jl) generally uses handwritten lists so it’s a bit more curated. Here’s an example:

[https://benchmarks.sciml.ai/html/StiffODE/Pollution.html](https://benchmarks.sciml.ai/html/StiffODE/Pollution.html)