Make CheckInit the default initialization algorithm for DAEs (breaking)

[Ickaser]

Checklist

• Appropriate tests were added
• Any code changes were done in a way that does not break public API
• All documentation related to code changes were updated
• The new code follows the
  contributor guidelines, in particular the SciML Style Guide and
  COLPRAC.
• Any new documentation only uses public API

Additional context

Re: discussion in #2513 , here is a (first-time-contributor, inexperienced) attempt at implementing the solution @ChrisRackauckas suggested.

[Ickaser]

This currently doesn't work--is too eager to trigger the error.

[Ickaser]

Now, in DiffEqBase.initialize_dae!, before calling _initialize_dae!, there is a hasmethod check whether a method exists for BrownFullBasicInit, which is (as far as I can tell) one of the two default algorithms (along with ShampineCollocationInit). If no such method exists (which is used as a proxy for asking if OrdinaryDiffEqNonlinearSolve has been loaded), then an informative error is thrown.

At least as far as the MWE from #2513 is concerned, this would mean that even if you pass a consistent initial condition (so that initialization isn't necessary) and an "algorithm" like CheckInit, it will still throw the error (even if the solver could run just fine).

An alternative I see is to, instead, do the check and throw the error in the methods of _initialize_dae! that handle DefaultInit; that would perhaps be a more robust way to go, since it means confining the error to cases where no initialization algorithm is specified (and not the cases where e.g. NoInit or CheckInit is passed), but means that the error lives in four places instead of one. I don't feel qualified to judge whether or not that would be preferable.

[oscardssmith]

I think the right approach here is probably to make the DefaultInit use CheckInit rather than BrownFullBasicInit, but this would be a breaking change. @ChrisRackauckas does this seem like a viable path forward? I think it is reasonable that we shouldn't mess with the user provided u0 unless the user asked us to.

[ChrisRackauckas]

I'm not sure. Most DAE solvers do just default to BrownBasic. I would be fine with doing this change in OrdinaryDiffEq v7 if the error message in CheckInit is improved to explain this.

[Ickaser]

Speaking as a user, I often work with a system where I often don't make any attempt to have u0 be a consistent initial condition (because I would have to solve a nonlinear equation to do so), and I'm willing to accept an extra dependency for the convenience of not needing to specify an initialization algorithm every time. (But if I've misunderstood, and the Rosenbrock methods only need NonlinearSolve for the initialization, it would be nice to cut down precompilation time by leaving that out.)

But there is definitely a logic to defaulting to CheckInit, so I wouldn't mind as long as the error message that CheckInit throws on inconsistent conditions provides some clear guidance about what has changed and how to set the algorithm to BrownFullBasicInit or NoInit

If we go that route, there will need to be an export for NoInit and CheckInit, because currently I have to reach into OrdinaryDiffEqRosenbrock.SciMLBase.NoInit() and OrdinaryDiffEqRosenbrock.OrdinaryDiffEqCore.CheckInit() to get those two options if I haven't loaded ODENonlinearSolve.

[ChrisRackauckas]

If we go that route, there will need to be an export for NoInit and CheckInit, because currently I have to reach into OrdinaryDiffEqRosenbrock.SciMLBase.NoInit() and OrdinaryDiffEqRosenbrock.OrdinaryDiffEqCore.CheckInit() to get those two options if I haven't loaded ODENonlinearSolve.

Yeah those should get exported. It's kind of a historical accident that they do not.

But there is definitely a logic to defaulting to CheckInit, so I wouldn't mind as long as the error message that CheckInit throws on inconsistent conditions provides some clear guidance about what has changed and how to set the algorithm to BrownFullBasicInit or NoInit

Yes, its error just needs to clearly say:

Initialization failed. Your u0 did not satisfy the initialization requirements. If you wish for the system to automatically change the algebraic variables to satisfy the algebraic constraints, please pass initializealg = BrownBasicInit() to solve (this option will require using OrdinaryDiffEqNonlinearSolve). If you wish to perform an initialization on the complete u0, please pass initializealg = ShampineCollocationInit() to solve. Note that initialization can be a very difficult process for DAEs and in many cases can be numerically intractable without symbolic manipulation of the system. For an automated system that will generate numerically stable initializations, see ModelingToolkit.jl structural simplification for more details.

[Ickaser]

On further reflection, I think I (as a user) would probably prefer to get an error message that suggests using OrdinaryDiffEqNonlinearSolve and keep BrownFullBasicInit as the default. But just to keep this moving, have added a commit where all instances of BrownFullBasicInit are replaced with CheckInit.

1. Using NoInit might be a valid option for some of my use cases, though (where I don't much care about early time behavior of the algebraic variable), so it would still be helpful to have that exported. It looks like exports of NoInit and CheckInit belong in either DiffEqBase or SciMLBase; feat: add implementations of CheckInit and OverrideInit SciMLBase.jl#845 makes it seems like SciMLBase might be the better place, since that is now where they are defined. Should I make a pull request there for that?

2. Also, while I'm looking at this code, I notice that there are two identical methods for _initalize_dae!(..., x::Val{true}) and _initialize_dae!(..., x::Val{false}). Is there a reason those aren't a single method with _initialize_dae!(..., x::Union{Val{true}, Val{false}})?

3. I notice that in refactor: generalize _initialize_dae! to use SciMLBase implementations #2521 @AayushSabharwal adjusts some of this behavior already; should I be waiting until his work is complete to dedicate more time to this?

[ChrisRackauckas]

Using NoInit might be a valid option for some of my use cases, though (where I don't much care about early time behavior of the algebraic variable)

No, that's not how it works. If you don't initialize correctly, it is more than likely that your next step is either unstable or the method is non-convergent. You can have a small error in the starting algebraic variable lead to an O(1) global error in the DAE solve, it's simply unbounded. You are required to start from a consistent state to have any sense of theoretical convergence involved at all, so using NoInit means we give zero guarantee on accuracy. That cannot be a default, it's only an internal for implementing other algorithms that then chain other initializations. It's not really meant to be used by users as a true initialization method, instead they should CheckInit.

Also, while I'm looking at this code, I notice that there are two identical methods for _initalize_dae!(..., x::Val{true}) and _initialize_dae!(..., x::Val{false}). Is there a reason those aren't a single method with _initialize_dae!(..., x::Union{Val{true}, Val{false}})?

Are there identical ones? Usually they are rather different code?

I notice that in #2521 @AayushSabharwal adjusts some of this behavior already; should I be waiting until his work is complete to dedicate more time to this?

If we do a global CheckInit default, then it would make sense to use the results of that.

[Ickaser]

Are there identical ones? Usually they are rather different code?

OrdinaryDiffEq.jl/lib/OrdinaryDiffEqCore/src/initialize_dae.jl

Lines 55 to 105 in ce43d57

	## Default algorithms
	function _initialize_dae!(integrator, prob::ODEProblem,
	alg::DefaultInit, x::Val{true})
	if SciMLBase.has_initializeprob(prob.f)
	_initialize_dae!(integrator, prob,
	OverrideInit(integrator.opts.abstol), x)
	else
	_initialize_dae!(integrator, prob,
	BrownFullBasicInit(integrator.opts.abstol), x)
	end
	end
	function _initialize_dae!(integrator, prob::ODEProblem,
	alg::DefaultInit, x::Val{false})
	if SciMLBase.has_initializeprob(prob.f)
	_initialize_dae!(integrator, prob,
	OverrideInit(integrator.opts.abstol), x)
	else
	_initialize_dae!(integrator, prob,
	BrownFullBasicInit(integrator.opts.abstol), x)
	end
	end
	function _initialize_dae!(integrator, prob::DAEProblem,
	alg::DefaultInit, x::Val{false})
	if SciMLBase.has_initializeprob(prob.f)
	_initialize_dae!(integrator, prob,
	OverrideInit(integrator.opts.abstol), x)
	elseif prob.differential_vars === nothing
	_initialize_dae!(integrator, prob,
	ShampineCollocationInit(), x)
	else
	_initialize_dae!(integrator, prob,
	BrownFullBasicInit(integrator.opts.abstol), x)
	end
	end
	function _initialize_dae!(integrator, prob::DAEProblem,
	alg::DefaultInit, x::Val{true})
	if SciMLBase.has_initializeprob(prob.f)
	_initialize_dae!(integrator, prob,
	OverrideInit(integrator.opts.abstol), x)
	elseif prob.differential_vars === nothing
	_initialize_dae!(integrator, prob,
	ShampineCollocationInit(), x)
	else
	_initialize_dae!(integrator, prob,
	BrownFullBasicInit(integrator.opts.abstol), x)
	end
	end

Strictly speaking, there are 4 methods here; a 2x2 dispatch on ODEProblem vs DAEProblem and isinplace true vs false. There are differences between ODEProblem and DAEProblem, but after replacing SciML.DefaultInit with an appropriate initialization, isinplace just gets passed on as is.

No, that's not how it works. If you don't initialize correctly, it is more than likely that your next step is either unstable or the method is non-convergent.

That makes sense--my intuition was incorrectly informed by the too-simple MWE I have been testing with, which produces correct answers even with an incorrect initial condition. So NoInit probably shouldn't be exported then, just CheckInit?

1. It looks like exports of NoInit and CheckInit belong in either DiffEqBase or SciMLBase; feat: add implementations of CheckInit and OverrideInit SciMLBase.jl#845 makes it seems like SciMLBase might be the better place, since that is now where they are defined. Should I make a pull request there for that?

[ChrisRackauckas]

So NoInit probably shouldn't be exported then, just CheckInit?

Yes.

What is left here? I'm a bit confused as to the current state.

[Ickaser]

What is left here? I'm a bit confused as to the current state.

I'm partly waiting to let #2521 proceed, because that affects where this code should land: if CheckInit lives in SciMLBase then the error message needs to go there, not in OrdinaryDiffEq.

I'm still uncertain of a couple things:

1. Should I go forward with making CheckInit the default, rather than BrownFullBasicInit or ShampineCollocationInit? My feeling is that this seems like an annoying change, if people are going to want BrownFullBasicInit most of the time anyway and just want to know that they need to import ODENonlinearSolve (at least that is how I feel, as a user). But I am happy to move forward with changing the default if desired.
2. To export CheckInit, should that export statement be placed in SciMLBase where CheckInit is defined? Or further down in OrdinaryDiffEq somewhere?

[oscardssmith]

good news: #2521 just proceeded :)

[ChrisRackauckas]

To export CheckInit, should that export statement be placed in SciMLBase where CheckInit is defined? Or further down in OrdinaryDiffEq somewhere?

Yes, SciMLBase should export it. Solvers always reexport SciMLBase, so that would cover it.

[ChrisRackauckas]

Should I go forward with making CheckInit the default, rather than BrownFullBasicInit or ShampineCollocationInit? My feeling is that this seems like an annoying change, if people are going to want BrownFullBasicInit most of the time anyway and just want to know that they need to import ODENonlinearSolve (at least that is how I feel, as a user). But I am happy to move forward with changing the default if desired.

Yeah it's really a feel thing. I think it's the way to go since that's the only way we could ever do the reduced dependencies, and so as long as that error message is really nice I think that would be good. I think down the line it's better, since right now what we do is assume "I know your algebraic variable initial conditions are just guesses", which could be a wrong take (and I've gotten some questions about "why did it change this initial condition? That's not correct!"). The error would explain the routes you can take, and making it explicit "please ignore my initial condition on all algebraic variables and use them as guesses" is really safer. We didn't go that way at the start because we thought it might be too hard for most users to do the initialization stuff, so we fully automated it, but now I think a very informative error message is a better idea instead of making an implicit assumption.

[Ickaser]

Just opened SciML/SciMLBase.jl#886 to export CheckInit and add an informative error as discussed here. So now, what this pull request does is:

• make CheckInit the default DAE initialization algorithm for all ODE and DAE problems
• BrownFullBasicInit and ShampineCollocationInit are exported by OrdinaryDiffEqNonlinearSolve, so that users can pass them to solve

This is now a breaking change. Appropriate tests have not been added yet, and documentation will need to be updated as well.

What this pull request doesn't do any more is provide an error message as originally suggested in #2513. So should there be a separate pull request providing that information for DiffEq <7.0 ?

[ChrisRackauckas]

Alright, so this should go with the OrdinaryDiffEq v7

[ChrisRackauckas]

So should there be a separate pull request providing that information for DiffEq <7.0 ?

It should now already be included by the SciMLBase part?

[ChrisRackauckas]

BrownFullBasicInit and ShampineCollocationInit are exported by OrdinaryDiffEqNonlinearSolve, so that users can pass them to solve

This would be fine pre v7

[Ickaser]

So should there be a separate pull request providing that information for DiffEq <7.0 ?

It should now already be included by the SciMLBase part?

The original error in #2513 that prompted this pull request was that, if ODENonlinearSolve is not imported, _initialize_dae had no method for the default algorithm (BrownFullBasicInit), so what you got was a MethodError (rather than an informative error message about needing using ODENonlinearSolve). This pull request fixes that by changing the default algorithm to something that's in ODECore, but for <7.0 that MethodError will still get hit.

Is it worth splitting this pull request in two, one for the breaking change that fixes things elegantly (7.0) and another for the non-breaking stuff plus an informative error message (<7.0)?

[ChrisRackauckas]

Is it worth splitting this pull request in two, one for the breaking change that fixes things elegantly (7.0) and another for the non-breaking stuff plus an informative error message (<7.0)?

Yes