v4.19.4 CVODE recoverable warnings now cause failures

[sjdaines]

Looks like with the Sundials.jl v4.19.4 release what was previously a CVODE warning followed by recovery has now become an error and then failure

The CVODE warning is

  Internal t = -1e+09 and h = 6.17217e-13 are such that t + h = t on the next step. The solver will continue anyway.

which prior to v4.19.4 did just what it says (continued), but now errors.

Plausibly the change is PR #416, https://github.com/oscardssmith/Sundials.jl/blob/35abe0140c6b70ff962e2971e024a75b0af04552/src/common_interface/solve.jl#L1403C12-L1405C16

[sjdaines]

A failure is visible in CI log from https://github.com/PALEOtoolkit/PALEOcopse.jl/actions/runs/6207071757/job/16851946212 :

┌ Info: 
│ ================================================================================
│ PALEOmodel.ODE.integrate:
│     tspan: (-1.0e9, 0)
│     algorithm: Sundials.CVODE_BDF{:Newton, :Dense, Nothing, Nothing}(0, 0, 0, false, 10, 5, 7, 3, 10, nothing, nothing, 0)
│     Jacobian: NoJacobian
│     steadystate: false
│     using 1 BLAS threads
└ ================================================================================
[ Info: ODEfunction: using Jacobian NoJacobian
[ Info: jac_config_ode: jac_ad=NoJacobian

[CVODES WARNING]  CVode
  Internal t = -1e+09 and h = 6.17217e-13 are such that t + h = t on the next step. The solver will continue anyway.

  3.526897 seconds (1.44 M allocations: 85.717 MiB, 0.80% gc time)
┌ Warning: `sol.destats` is deprecated. Use `sol.stats` instead.
│   caller = ip:0x0
└ @ Core :-1
┌ Info: 
│ ================================================================================
│ print_sol_stats:
│     retcode=Unstable
│     nsteps 1
│ 
│     alg=Sundials.CVODE_BDF{:Newton, :Dense, Nothing, Nothing}(0, 0, 0, false, 10, 5, 7, 3, 10, nothing, nothing, 0)
│     DiffEqBase.Stats
│ Number of function 1 evaluations:                  2
│ Number of function 2 evaluations:                  0
│ Number of W matrix evaluations:                    1
│ Number of linear solves:                           0
│ Number of Jacobians created:                       1
│ Number of nonlinear solver iterations:             1
│ Number of nonlinear solver convergence failures:   0
│ Number of rootfind condition calls:                0
│ Number of accepted steps:                          1
│ Number of rejected steps:                          0
│     length(sol.t) 1
│     size(sol) (19, 1)
└ ================================================================================

[ChrisRackauckas]

@oscardssmith can you look into that?

[oscardssmith]

this was an intentional change to match the rest of the ode ecosystem. also, how is the integration not getting stuck in an infinite loop if this happens? is sundials storing a double double time internally?

[sjdaines]

Or perhaps sundials has some additional internal state (numerical Jacobian? ??) which is modified/reset?

Would it be reasonable just to keep calling sundials CVode until it either returns an error code or the maximum number of timesteps is reached?

[ivborissov]

The same problem. After Sundials update to >=v4.19.4 the model returns Unstable retcode with CVODE_BDF(). Previouly CVODE_BDF() handled this model with Success retcode, but printing dt warnings as mentioned above:

Internal t = -1e+09 and h = 6.17217e-13 are such that t + h = t on the next step. The solver will continue anyway.

Setting dtmin=0. doesn't help Sundials here.
FBDF() solves it, but it is ~7 times slower on this problem. Unfortunaly, can't share the code of the model, but will be glad to contribute to this as well as to making FBDF() faster, as CVODE_BDF() is indeed far from being accurate.