WIP: Benchmarks for ExpRK and IMEX methods

[MSeeker1340]

nbview for this branch

This is initially intended to resolve #30 but I find it natural to combine both ExpRK and IMEX methods, and also work on tailored semilinear problems instead of appending to exisiting stiff ODE benchmarks. The reason is that these methods aren't really suited to most of the existing benchmark equations: they are usually highly non-linear, has only a few variables and sometimes large positive eigenvalues in their Jacobian. On the other hand, ExpRK/IMEX's target problems are usually semilinear PDE (thus high dimension), with stiffness in the form of large negative eigenvalues in the linear term. (For EPIRK, although it doesn't use the semilinear form, the other criteria still apply). The best we can do about the exisiting stiff ODE benchmarks is to add caching ExpRK methods, which requires rewriting the equations in semilinear form and is unstable if the linear term has large positive eigenvalues (e.g. Van der Pol).

The equations I want to benchmark are (mostly taken from the Tokman EPIRK paper):

• Burgers equation
• Allen-Cahn equation (1D & 2D)
• 2D Brusselator (use existing implementation in DiffEqProblemLibrary?)
• (maybe) nonlinear Schrodinger equation; ExpRK should work in complex but I'm not sure about IMEX

The 2D problems can be a bit tricky as DiffEqOperators support hasn't followed up yet. But I can always hand-write the matrix-free derivative operators for the Krylov ExpRK and IMEX methods (both work well with matrix-free operators). And since the derivative operators are just convolutions it might not be far fetched to use GPU acceleration for that (I presume that it's convenient to wrap a mul! call with cuBLAS matmul in Julia?). And speaking of matrix-free operators, spectral methods and ApproxFun also comes to mind.

BTW, the notebooks in this repo is in dire need of an update...

[ChrisRackauckas]

Great to see this. Yes these are good test equations. However, you are missing all of the interesting comparisons. You need to compare between classes of methods. For example, from the plots it's not clear if one should be using CNAB2 or ETDRK2, or when one should be switching to ETDRK4. It's not clear how these compare to other methods like Rosenbrock23, CVODE_BDF, or Rodas5. Within class benchmarks are used to establish a "champion", and then you show the champion of each class of methods against each other to determine performs best.

Two things to mention about

setups = [Dict(:alg => KenCarp3(linsolve=LinSolveFactorize(lu))),
          Dict(:alg => KenCarp4(linsolve=LinSolveFactorize(lu))),
          Dict(:alg => KenCarp5(linsolve=LinSolveFactorize(lu))),
          Dict(:alg => ARKODE(Sundials.Implicit(), order=3, linear_solver=:GMRES)),
          Dict(:alg => ARKODE(Sundials.Implicit(), order=4, linear_solver=:GMRES)),
          Dict(:alg => ARKODE(Sundials.Implicit(), order=5, linear_solver=:GMRES))]

For one, we should make our default handling automatically do well on sparse matrices and switch to lu. But secondly, you should time ARKODE with and without GMRES, and do one with IterativeSolvers.jl's GMRES on KenCarps. From your benchmarks I do not know anything about the performance differences. Is ARKODE faster, or is using a GMRES linear solver just better on this problem? (Given the previous benchmarks between these methods I'd guess the latter!)

[MSeeker1340]

@ChrisRackauckas Not sure why using GMRES on KenCarp3/4/5 exceeds the maxiter limit. Might be worth looking into?

For the implicit/ExpRK comparison: To gain more insight (say for the presentation next year) we should probably probe into the integrators and count how many Krylov iterations are carried out for each time step. The fact that ExpRK (in its current state) is fixed-dt may also contribute to their under-performance. I'll be curious to see how the adaptive Exprb (and EPIRK methods to some extend) compares with regular implicit methods on a non-semilinear formulated problem.

[ChrisRackauckas]

Not sure why using GMRES on KenCarp3/4/5 exceeds the maxiter limit. Might be worth looking into?

Yes, that's worth looking into. Maybe tweak some of the algorithm options around and see what happens? Take a look at the options on ARKODE and see if changing the KenCarp one closer to the ARKODE defaults leads to better results there. Their defaults are really really bad for highly nonlinear ODE, but I suspected they must do better on semi-stiff PDE. This seems to be the case. The algorithms are very similar except for these defaults so I'd like to see if we can get them inline or even faster. Though we may also want to play with step acceleration on KenCarp. This is a benchmark worth optimizing on.

To gain more insight (say for the presentation next year) we should probably probe into the integrators and count how many Krylov iterations are carried out for each time step.

Yes. SciML/DiffEqBase.jl#54 . This would be a good issue for us to get solved, but it would require a lot of work since every method has a different number of f calls per perform_step!, and some (implicit methods) have variable number of f calls...

[MSeeker1340]

Ran the benchmarks on the Umaryland server. Also changed the GMRES linear solvers to only use 5 Krylov iterations (which is the default of ARKODE). I still get "maxiters reached" warning for KenCarp3/4/5 but the results seemed ok.

One thing to note is that unlike on my local machine, ARKODE and ETDRK4 actually doesn't differ that much. Maybe a result of BLAS optimizations?

[ChrisRackauckas]

One thing to note is that unlike on my local machine, ARKODE and ETDRK4 actually doesn't differ that much. Maybe a result of BLAS optimizations?

Sundials can't use BLAS with 64-bit (I know, right? Right?), so our methods will scale better to more cores. This is seen in other benchmarks as well.

[ChrisRackauckas]

I still get "maxiters reached" warning for KenCarp3/4/5 but the results seemed ok.

They aren't exiting early?

[ChrisRackauckas]

Incorporated into #48 . Thanks!