Type instabilities in DAE solves

[Ickaser]

Describe the bug 🐞

Summarizing what I learned in #2594 and #2698:

• Direct construction of an ODEFunction, as is necessary to specify a mass matrix, defaults to FullSpecialize, where constructing an ODEProblem defaults to AutoSpecialize.
• With AutoSpecialize, autodiff should default to being an inferrable AutoForwardDiff(chunksize=1). In some (but not all) test cases (e.g. OrdinaryDiffEq.jl/test/interface/mass_matrix_tests.jl, and others marked with @test_broken) this does not seem to work as it should; in one test case, the in-place version is not inferred while the out-of-place version is. With FullSpecialize, the chunksize in general needs to be specified for the algorithm and solution type to be inferred.
• In addition to inferring the algorithm (and autodiff) type, successful solution inference requires that get_differential_vars be inferrable. This works if the provided mass matrix is a UniformScaling or is Diagonal; other cases are handled but depend on the values of the mass matrix, so they aren't known at compile time.

Questions to answer

• Should ODEFunction default to AutoSpecialize?
• Is it worth doing linear algebra in the preparation step of get_differential_vars in order to rigorously identify all algebraic variables? (Probably not, but is a call to make.)

Issues to address

• Infer a default chunksize of 1 for all cases with AutoSpecialize, particularly in-place
• Adjust existing mass-matrix code to utilize Diagonal wherever possible
• OrdinaryDiffEqBDF algorithms for DAEs (DFBDF(), etc) are not currently inferred in the test cases, even with specified chunk sizes.
• OrdinaryDiffEqBDF algorithms with mass matrix (FBDF, etc.) are not currently differentiable or inferred in test cases
• Develop a way to make get_differential_vars more inferrable
• Add advice to documentation that assists users to make their DAE solves more inferrable

Personal note
I've committed a lot of time to this in the last few days, and need to start thinking about my regular PhD research again 😅 , so I will probably not be on this immediately

[ChrisRackauckas]

Should ODEFunction default to AutoSpecialize?

I think so. That would decrease latency requirements because AutoSpecialize means that you store the precompilation cache, and "normally" that's the right default.

Is it worth doing linear algebra in the preparation step of get_differential_vars in order to rigorously identify all algebraic variables? (Probably not, but is a call to make.)

Probably not. It's not even clear what the right linear algebra would be 😅, but I'd like to figure it out because I'm curious. I feel like I'm forgetting one piece.

[ChrisRackauckas]

I've committed a lot of time to this in the last few days, and need to start thinking about my regular PhD research again 😅 , so I will probably not be on this immediately

Understood. We got to a decent stopping point there, and there's always more to do. With DAEs, there's still infinite amounts of things to improve 😅