SciML / OrdinaryDiffEq.jl

High performance ordinary differential equation (ODE) and differential-algebraic equation (DAE) solvers, including neural ordinary differential equations (neural ODEs) and scientific machine learning (SciML)
https://diffeq.sciml.ai/latest/
Other
533 stars 200 forks source link

Successive Newton iteration GMRES starting points #724

Open ChrisRackauckas opened 5 years ago

ChrisRackauckas commented 5 years ago

GMRES can use an initial condition. In the coming defaults, the x part is zero'd out in the start because the Rosenbrock methods don't necessarily have a good guess for the k vectors (or do they? We should test this out). But in Newton-Krylov iterations, the previous Jv=b calculation might give a nice guess for v, so we'd want to override the zeroing in the linsolve and let it know to use the initial condition. Maybe this can even be preserved between steps?

ChrisRackauckas commented 5 years ago

This seems to be what Sundials does as well. https://github.com/LLNL/sundials/blob/master/src/sunlinsol/spgmr/sunlinsol_spgmr.c#L486-L487

It initializes the correction vector at 0. It also resets the Krylov dimension at 0 each time:

https://github.com/LLNL/sundials/blob/master/src/sunlinsol/spgmr/sunlinsol_spgmr.c#L405-L407

Such things don't seem to be mentioned in https://computation.llnl.gov/casc/nsde/pubs/221215.pdf

isaacsas commented 5 years ago

I think most of what Sundials does is documented (though it takes a fair amount of reading) within the "User Documentation for CVODE" pdf. Unfortunately it is over 300 pages, so a fair amount of work to parse through. My reading of it is consistent that they reset the guess to zero and reinitialize the Krylov space each time, but also that they solve a (diagonal matrix) rescaled system where the rescaling is recalculated at each Newton iteration step (and based on the solution-dependent weights they use in their error norm). I think they also fix their GMRES solves to only five iterations with no restarts.

There are also some differences in that they stop GMRES based on the absolute residual I think, with a stopping tolerance that is related to the Newton method tolerance and using their weighted error norm.