Benchmarking all algorithms in DiffEq/Optimisation code

General Usage Performance
, ,
1

If I have some ODE code, say this code from the https://diffeq.sciml.ai/stable/tutorials/advanced_ode_example/:

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) for a given optimisation problem? Or is it again just done by looping over a hand-written list?

1 Like
2

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

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 generally uses handwritten lists so it’s a bit more curated. Here’s an example:

https://benchmarks.sciml.ai/html/StiffODE/Pollution.html

1 Like