From 89a643db88e056d7918931d21243894a7ee170d6 Mon Sep 17 00:00:00 2001 From: Daines Date: Sat, 23 Sep 2023 08:00:51 +0100 Subject: [PATCH] Fix handling of small timesteps Sundials solver is allowed to warn and continue if it returns a timestep with no change in t. If the solver doesn't recover, it will fail with eg SciMLBase.ReturnCode.MaxIters Reverts PR https://github.com/SciML/Sundials.jl/pull/416, fixes https://github.com/SciML/Sundials.jl/issues/420 For CVode and ARKode, it looks like this is intended behaviour: the Sundial solver emits a warning message, with an API call CVodeSetMaxHnilWarns ARKStepSetMaxHnilWarns to set the maximum number of warning messages printed, which is exposed to Julia as max_hnil_warns see Sundials driver code https://github.com/LLNL/sundials/blob/e8a3e67e3883bc316c48bc534ee08319a5e8c620/src/cvode/cvode.c#L1339-L1347 For IDA, there is no API like this so it's not so clear what should happen, however the behaviour prior to https://github.com/SciML/Sundials.jl/pull/416 (restored by this PR) for the f_noconverge test added in test/common_interface/ida.jl is to return SciMLBase.ReturnCode.MaxIters which seems reasonable? (NB: @oscardssmith I can't reproduce the issue implied by the title of https://github.com/SciML/Sundials.jl/pull/416 so perhaps what is missing here is the original failure case that demonstrates the issue? The call to solve in the f_noconverge test test/common_interface/ida.jl, returns SciMLBase.ReturnCode.MaxIters, it doesn't return with no error? Also https://github.com/SciML/SciMLBase.jl/issues/458 shows IDA printing [IDAS ERROR] IDASolve At t = 0 and h = 1.49012e-18, the error test failed repeatedly or with |h| = hmin. which suggests that IDA did flag an error in that case ? ) --- src/common_interface/solve.jl | 4 +--- test/common_interface/ida.jl | 2 +- 2 files changed, 2 insertions(+), 4 deletions(-) diff --git a/src/common_interface/solve.jl b/src/common_interface/solve.jl index e3e6f515..a5974d32 100644 --- a/src/common_interface/solve.jl +++ b/src/common_interface/solve.jl @@ -1406,9 +1406,7 @@ function DiffEqBase.solve!(integrator::AbstractSundialsIntegrator; early_free = integrator.userfun.p = integrator.p solver_step(integrator, tstop) integrator.t = first(integrator.tout) - if integrator.t == integrator.tprev - integrator.flag = -3 - end + # NB: CVode, ARKode may warn and then recover if integrator.t == integrator.tprev so don't flag this as an error integrator.flag < 0 && break handle_callbacks!(integrator) # this also updates the interpolation integrator.flag < 0 && break diff --git a/test/common_interface/ida.jl b/test/common_interface/ida.jl index 2c3c3a10..b7a3e8b8 100644 --- a/test/common_interface/ida.jl +++ b/test/common_interface/ida.jl @@ -117,4 +117,4 @@ isapprox(only(sol.u[end]), exp(1), rtol = 1e-3) f_noconverge(out, du, u, p, t) = out .= [du[1] + u[1] / (t - 1)] prob = DAEProblem(f_noconverge, [1.0], [1.0], (0, 2); differential_vars = [true]) sol = solve(prob, IDA()) -@test !(sol.retcode in (ReturnCode.Success, ReturnCode.MaxIters)) +@test !(sol.retcode in (ReturnCode.Success, ))