OrdinaryDiffEq v7

[ChrisRackauckas]

OrdinaryDiffEq v7 is close. This means that all breaking changes need to make it in by this time. Here's the list right now:

• Split solvers into subpacakges Refactor OrdinaryDiffEq to use Subpackages for Solver Sets #2177 . This is the main change. After this change, we have to now define what we want OrdinaryDiffEq.jl to actually look like. The current plan is to have an OrdinaryDiffEqCore.jl at the bottom which is just the core integrator, then OrdinaryDiffEqDifferentiation.jl on top of that, then OrdinaryDiffEqNonlinearSolve.jl, and then the algorithms depend on 1-3 depending on how much they need. This allows the explicit RK methods to not require ForwardDiff or nonlinear solvers, effectively cutting out most dependencies. Then OrdinaryDiffEqDefault is the default method built from some solver packages, and then OrdinaryDiffEq.jl then becomes a shell with just the most widely used methods like Tsit5, Vern, Rosenbrock, and FBDF re-exported. All other solvers are only available when explicitly importing the specific sublibrary.
• The preconditioner change make preconditioners part of the solver rather than a random extra LinearSolve.jl#514 which @oscardssmith is doing needs to be in because it's a breaking change to how preconditioners are defined in the algorithm
• We will need to really complete the OrdinaryDiffEq documentation because we will need to document algorithms per-module, and we should do this via docstrings, I am asking for help from @ArnoStrouwen
• We want to change autodiff specification to using ADTypes. @gdalle and @avik-pal have already kind of been pushing things down this route. This will make changes to use Enzyme and DifferentiationInterface be non-breaking. Since this has been around for so long, we should have a 2 year deprecation period on the booleans though, but in the constructors we can just change the booleans to AutoForwardDiff and AutoFiniteDiff as a simple path. But with everything ADTypes the integration with NonlinearSolve is simpler as well.
• Default lazy to a Val so we can remove f from non-lazy interpolations, making Leaner serialization for solution types SciMLBase.jl#629 simpler and more robust. I don't think this is necessary breaking?
• Possibly switch default forward autodiff to Enzyme?
• Possibly switch nonlinear solve default to NonlinearSolve?
• Make CheckInit the default initialization algorithm for DAEs (breaking) #2514

Any other major changes to consider @devmotion @YingboMa @ranocha @oscardssmith

For reference, the other big thing is SciMLOperators, but those changes are breaking there and non-breaking here.

Closes #1886

[ChrisRackauckas]

I really want to have Enzyme as default forward mode in v7 as well? But current discussions with @wsmoses and @gdalle suggest it may not make it, as this is the whole caching thing. But if we make a system which is robust and PreallocationTools.jl is not required, then that's a clear want for v7.

[ChrisRackauckas]

I would also like to see #1570 completed to the point that we can change defaults to nonlinearsolve, but we need to make sure that doesn't have a performance regression (needs some jac reuse work in that interface). I don't know if that will be done in time @termi-official @oscardssmith

[oscardssmith]

I don't think #1570 is required for v7. It's not a breaking change.

[gdalle]

As far as DifferentiationInterface is concerned, multi-argument support may not be available right away (JuliaDiff/DifferentiationInterface.jl#311). Is it necessary for OrdinaryDiffEq, or would single-argument functions suffice in most cases?

[ChrisRackauckas]

Only single argument here.

[wsmoses]

I really want to have Enzyme as default forward mode in v7 as well? But current discussions with @wsmoses and @gdalle suggest it may not make it, as this is the whole caching thing. But if we make a system which is robust and PreallocationTools.jl is not required, then that's a clear want for v7.

I mean this can just fall back to doing if Enzyme.guaranteed_const(typeof(f)) autodiff(Fwd, Const(f), ...) else autodiff(Fwd, Dup(f, make_zero(f)), ...) in interim. The new default will just be autodiff(Fwd, Const(f) but adding checks to see if f is stored into.

There are separate questions about an optimization-aware preallocation tools [since presently preallocationtools caches appear as existing allocations which could alias with all data]

[gdalle]

Apart from the discussion around constant functions in Enzyme, I would say single-argument DI is pretty much ready for ODE adoption.
The question is whether I make a breaking change right now to make things easier in the future for multi argument, or whether I keep the current API and try to stick to it when multi argument comes around. Thoughts on API design?

Either way, I can help review the DI insertion but I can't launch the PRs from scratch. I have never solved an ODE in Julia, let alone peeked at the code base, so I'll need someone else to give it a first shot before I can be useful.

[ranocha]

Is there some kind of overview of the distribution of individual ODE solvers to subpackages?

[ranocha]

Another aspect to consider might be the adaptivity/controller approach

[wsmoses]

Either way, I can help review the DI insertion but I can't launch the PRs from scratch. I have never solved an ODE in Julia, let alone peeked at the code base, so I'll need someone else to give it a first shot before I can be useful.

Same here with integrating Enzyme FWD mode. If you have any code which does this or issues, I can try to resolve and/or comment best usage practices

[ranocha]

How difficult would it be to enable the usage of arrays of SVectors again for the explicit RK subpackages?

CC @jlchan

[ChrisRackauckas]

Is there some kind of overview of the distribution of individual ODE solvers to subpackages?

Since the split is "complete", you can look at the current libs. I did separate out some very high in-demand methods to be used in relative isolation, like Tsit5 and the Verner methods. We can always subdivide in new ways so I'm happy to entertain any other divisions with a PR. #2317 is currently trying to get documentation up for all of this. Let me know if there is any other splits that would help reduce what's in Trixi.

We did get the load and deps pretty minimal: https://sciml.ai/news/2024/08/10/sciml_small_grants_successes/ which should help the HPC use cases, but there's still a few bits I'm going to be exploring in there like making SparseArrays always an extension down the stack. I'm hoping we can get Tsit5 first run to 0.2 seconds.

[ChrisRackauckas]

Another aspect to consider might be the adaptivity/controller approach

Could you describe what you had in mind?

How difficult would it be to enable the usage of arrays of SVectors again for the explicit RK subpackages?

I'm still a bit weary of this because it's a breakage from the other interface-level tests https://docs.sciml.ai/SciMLBase/stable/interfaces/Array_and_Number/ that are starting to be imposed. I think it makes it quite difficult to enforce interfaces if we have a subset of algorithms that are allowing an interface break. I think the questions I have there are:

1. Is using a wrapper type like VectorOfArray a very large impediment to your workflows? At this point we have had the time to give it a try and iron out the bugs. I assumed that everything was working well after @jlchan 's PRs to RecursiveArrayTools.jl, but if there's still lingering issues then that's something to consider.
2. Is there a way to get this formalized into an interface that we could more broadly setup? For example, I don't see why explicit RK is special here, when for example Broyden methods should in principle be able to support them in the same way. If we were to go this route, it would be good to define a system whereby recursive arrays can have a set of defined behaviors that we could check are true, some defined utilities around them in SciMLBase, some algorithm traits to allow for algorithms to decide to opt into allowing the type, and then treating explicit RK methods as a set of methods which support this interface. And if we do this then we'd want to document it and start rolling that out to some of NonlinearSolve.jl and Optimization.jl as well.

[ranocha]

Another aspect to consider might be the adaptivity/controller approach

Could you describe what you had in mind?

The discussion we had on Slack with Théo Duez

[ranocha]

Is using a wrapper type like VectorOfArray a very large impediment to your workflows? At this point we have had the time to give it a try and iron out the bugs. I assumed that everything was working well after @jlchan 's PRs to RecursiveArrayTools.jl, but if there's still lingering issues then that's something to consider.

@jlchan would be the best to tell more about this. Right now, we still limit DiffEqBase.jl in Trixi.jl because of this. Maybe it's better now since LoopVectorization.jl continues to work in Juli v1.11 - there were several issues at once @jlchan wanted to fix if I remember correctly

[ranocha]

Is there a way to get this formalized into an interface that we could more broadly setup? For example, I don't see why explicit RK is special here, when for example Broyden methods should in principle be able to support them in the same way. If we were to go this route, it would be good to define a system whereby recursive arrays can have a set of defined behaviors that we could check are true, some defined utilities around them in SciMLBase, some algorithm traits to allow for algorithms to decide to opt into allowing the type, and then treating explicit RK methods as a set of methods which support this interface. And if we do this then we'd want to document it and start rolling that out to some of NonlinearSolve.jl and Optimization.jl as well.

I can't guarantee that it's 100% correct, but something along the following lines should be what we need from our site ((abstract) arrays of SVectors) and for explicit RK methods:

• u etc. can be used in broadcasting like @.. u = u + dt * du etc.
• There is some way to determine the underlying scalar type used for dt etc.
• There is a reasonable norm (for error-based step size control)
• similar etc. works already

From a mathematical point of view, we need the operations available in a normed vector space:

• determine the scalar field
• add vectors
• multiply vectors by scalars
• compute the norm

[ranocha]

Is there some kind of overview of the distribution of individual ODE solvers to subpackages?

Since the split is "complete", you can look at the current libs. I did separate out some very high in-demand methods to be used in relative isolation, like Tsit5 and the Verner methods. We can always subdivide in new ways so I'm happy to entertain any other divisions with a PR. #2317 is currently trying to get documentation up for all of this. Let me know if there is any other splits that would help reduce what's in Trixi.

I found some methods that appear to be in the wrong sub-package: #2373 I don't have time for a PR right now but maybe later this week

[ranocha]

Other than that, it looks like a nice improvement in a first test 👍 I'm still wondering why some packages like Test are loaded but it's much better than before where we had to load all of the implicit solvers - that we didn't need

[jlchan]

Is using a wrapper type like VectorOfArray a very large impediment to your workflows? At this point we have had the time to give it a try and iron out the bugs. I assumed that everything was working well after @jlchan 's PRs to RecursiveArrayTools.jl, but if there's still lingering issues then that's something to consider.

@jlchan would be the best to tell more about this. Right now, we still limit DiffEqBase.jl in Trixi.jl because of this. Maybe it's better now since LoopVectorization.jl continues to work in Juli v1.11 - there were several issues at once @jlchan wanted to fix if I remember correctly

I made a few PRs to RecursiveArrayTools.jl, but stopped at SciML/RecursiveArrayTools.jl#378 when things got busy. I also stopped because I wasn't sure how many VectorOfArray methods might have to special case to get things working.

[ChrisRackauckas]

I'm still wondering why some packages like Test are loaded

See my #gripes in the Julia Slack on that. It'll probably take a bit of time to battle against all of those little dependencies, but that's the next phase of the war.