Add Lobatto IIIa-c and Radau IIa solvers

[axla-io]

Features

• 17 new Fully Implicit Runge-Kutta (FIRK) solvers: Lobatto IIIa-c 2 to 5 stage versions, and Radau IIa 1,2,3,5,7 - stage versions
• Choice of cache with nested nonlinear solve (better performance for large problems) or expanded cache (seems to be more robust)
• Adaptivity and interpolations
• Sparse Jacobians

Issues with the code

Tests pass for the expanded cache, except for in some cases with the convergence tests for the higher order solvers, I think that the reason for this is that the error is low in absolute terms. I've done a check where the tests are marked as broken if the solution error is less than 1e-12

Also there is a runtime dispatch error when calling the JET tests. I have commented the lines in src/solve/firk.jl where the error occurs. I couldn't figure out how to resolve them.

For the nested cache, in addition to the tests that the expanded one doesn't pass, it also seem to have poor performance on NLLNS and Vector of vector initials tests. I haven't found a solver for which the nested ones converges there. For now these tests are commented out.

Also, the nested cache doesn't support ForwardDiff yet.

Checklist

• [x ] Appropriate tests were added
• [x ] Any code changes were done in a way that does not break public API
• [x ] All documentation related to code changes were updated
• The new code follows the
  contributor guidelines, in particular the SciML Style Guide and
  COLPRAC. (I commented out tests for nested cache as a WIP thing)
• Any new documentation only uses public API (I don't understand what this means :( )

[ChrisRackauckas]

Tests pass for the expanded cache, except for in some cases with the convergence tests for the higher order solvers, I think that the reason for this is that the error is low in absolute terms. I've done a check where the tests are marked as broken if the solution error is less than 1e-12

Can you plot what those look like? for higher order you likely need to increase the dts

[avik-pal]

Also let's merge the tests for MIRK and the other methods, except Shooting. The only reason Shooting is separate, is that the other methods currently don't support interpolation in the solution.

[ChrisRackauckas]

That looks like the opposite direction that the convergence plot recipe normally does. Is this not just plot(sim)?

[ChrisRackauckas]

Hit a time limit? We really need to split this into different test groups.

[ErikQQY]

Yeah, I will split them into different test groups

[ChrisRackauckas]

Is the Retest parallelism actually helpful? BLAS will multithread, so we may be slamming the threads on any large enough problem. I don't think that part makes sense in the BVP library, or ML @avik-pal ?

[avik-pal]

Is the Retest parallelism actually helpful? BLAS will multithread, so we may be slamming the threads on any large enough problem. I don't think that part makes sense in the BVP library, or ML @avik-pal ?

Are any of the tests here heavy BLAS? In the ML context, none of the Lux tests are heavy enough to cause issues with oversubscribing, and it helps bring down the build times to about 25 mins (compared to > 1hr) on decent machines (the cuda ones on buildkite). Not the most helpful for github actions.

That said for testing heavy workloads (both compute or memory), I do disable parallelism in testing. I still keep ReTestItems because it lets you isolate the testitem that fails and for testing Enzyme which can segfault/assert crash julia it will just crash the worker, so the test shows up as failed rather than the main process crashing.

[ChrisRackauckas]

The LU/QR will all multithread above 200x200, which I presume we must be hitting?

[avik-pal]

https://github.com/SciML/BoundaryValueDiffEq.jl/actions/runs/10617623160/job/29430931485?pr=178#step:6:397 the AD isn't robust to rank deficient matrices. This probably needs to be handled upstream

[ErikQQY]

https://github.com/SciML/BoundaryValueDiffEq.jl/actions/runs/10617623160/job/29430931485?pr=178#step:6:397 the AD isn't robust to rank deficient matrices. This probably needs to be handled upstream

From local testing, RadauIIa5 just error with GaussNewton, the other nonlinear solvers are just fine, I have just skipped this combination for now

[ChrisRackauckas]

That's the least stable one. Not surprised. It should really default to using a trust region method for this reason.

[ChrisRackauckas]

Oh no I thought it was going to merge 😅

[ErikQQY]

Yeah it was going to, but CI seemed to have a different opinion, let’s see how about now😅

[ErikQQY]

Finally ready to go🎉🎉

[ChrisRackauckas]

Amazing, and thanks @axla-io too.