LLNL / sundials

Official development repository for SUNDIALS - a SUite of Nonlinear and DIfferential/ALgebraic equation Solvers. Pull requests are welcome for bug fixes and minor changes.
https://computing.llnl.gov/projects/sundials
BSD 3-Clause "New" or "Revised" License
517 stars 126 forks source link

The purpose of jcur in PrecSetupFn #107

Closed ChrisRackauckas closed 2 years ago

ChrisRackauckas commented 2 years ago

The jcur in the PrecSetupFn has always confused me https://sundials.readthedocs.io/en/latest/cvode/Usage/index.html#preconditioner-setup-iterative-linear-solvers

jcur – a pointer to a flag which should be set to SUNTRUE if Jacobian data was recomputed, or set to SUNFALSE if Jacobian data was not recomputed, but saved data was still reused.

What exactly is the use case for this, and how does it effect the solver algorithm after the user triggers it? It cannot be for matrix-free Newton-Krylov because there's no Jacobian to recompute. But even if there is a Jacobian, where would that information even be used? Usually jcur is for knowing whether to perform a new factorization, but this is in preconditioners so there is no factorization. Is this data even used or is it simply for updating the counters for number of Jacobians calculated?

drreynolds commented 2 years ago

Hi Chris, you may have already figured this out since you partially answered this question in your explanation. It is frequently the case that a preconditioner is formed (by the user) based on an incomplete factorization of the Jacobian (e.g., ILU for the full Jacobian, or an exact factorization of a banded or block-diagonal Jacobian approximation). Thus although our "matrix-free Newton-Krylov" algorithms may indeed be matrix-free, your preconditioner need not be.

ChrisRackauckas commented 2 years ago

But why would the preconditioner need to signal back to the algorithm "I recomputed a Jacobian inside of this preconditioner call"? What you're describing seems to be what jok is for:

jok – an input flag indicating whether the Jacobian-related data needs to be updated. The jok argument provides for the reuse of Jacobian data in the preconditioner solve function. jok = SUNFALSE means that the Jacobian-related data must be recomputed from scratch. jok = SUNTRUE means that the Jacobian data, if saved from the previous call to this function, can be reused (with the current value of ). A call with jok = SUNTRUE can only occur after a call with jok = SUNFALSE.

jok is clear: if someone is doing an ILU then that should be tied into the Jacobian reuse schema so that way the ILU is only recomputed when a new Jacobian is calculated, and they could provide analytical Jacobian expressions for some approximation. That would seemingly cover that whole use case without the existence of jcur, unless I'm mistaken somewhere in that and Sundials does not allow tying the approximate Jacobian into the standard Jacobian machinery.

So then what's the case where the user is choosing to override the Jacobian update functionality inside of a preconditioner, and then how exactly would the time stepper want to make use of it? And where is that done in the code? Because for the most part it almost looks like the code is ignoring this signal back.

drreynolds commented 2 years ago

Thank you for clarifying your question -- I believe I now understand what you're asking.

We need to know how "current" the preconditioner in case the Newton iteration fails to converge. Specifically, if the Newton iteration diverges with a "stale" preconditioner, then we will first call psetup with jok = SUNFALSE so that you can recompute things before we try again. However, if the Newton iteration fails to converge with a "current" preconditioner, then we have no recourse but to reduce the time step size and retry the time step.

All of this is discussed in the CVODE documentation: https://sundials.readthedocs.io/en/latest/cvode/Mathematics_link.html#ivp-solution. in the paragraphs following equation (4.8). Similar discussion is included in the other user guides.

ChrisRackauckas commented 2 years ago

I do not see where in that documentation this detail is described, but the description here makes sense. Thanks!