From 429c570f7078147f2234a000d85c08b7463d1a16 Mon Sep 17 00:00:00 2001 From: oscarddssmith Date: Thu, 16 May 2024 11:42:14 -0400 Subject: [PATCH 1/7] improve dtmin and unstable detection --- src/integrator_interface.jl | 67 +++++++++++++++++++++++-------------- 1 file changed, 41 insertions(+), 26 deletions(-) diff --git a/src/integrator_interface.jl b/src/integrator_interface.jl index 9593baeb4..f41cae8f0 100644 --- a/src/integrator_interface.jl +++ b/src/integrator_interface.jl @@ -584,53 +584,68 @@ function check_error(integrator::DEIntegrator) if integrator.sol.retcode ∉ (ReturnCode.Success, ReturnCode.Default) return integrator.sol.retcode end + opts = integrator.opts + verbose = opts.verbose # This implementation is intended to be used for ODEIntegrator and # SDEIntegrator. if isnan(integrator.dt) - if integrator.opts.verbose + if verbose @warn("NaN dt detected. Likely a NaN value in the state, parameters, or derivative value caused this outcome.") end return ReturnCode.DtNaN end - if integrator.iter > integrator.opts.maxiters - if integrator.opts.verbose + if integrator.iter > opts.maxiters + if verbose @warn("Interrupted. Larger maxiters is needed. If you are using an integrator for non-stiff ODEs or an automatic switching algorithm (the default), you may want to consider using a method for stiff equations. See the solver pages for more details (e.g. https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/#Stiff-Problems).") end return ReturnCode.MaxIters end # The last part: - # If you are close to the end, don't exit: let the user hit the end! - # However, if we try that and the step fails, exit instead of infinite loop - if !integrator.opts.force_dtmin && integrator.opts.adaptive && - abs(integrator.dt) <= abs(integrator.opts.dtmin) && - (((hasproperty(integrator, :opts) && hasproperty(integrator.opts, :tstops)) ? - integrator.t + integrator.dt < integrator.tdir * first(integrator.opts.tstops) : - true) || (hasproperty(integrator, :accept_step) && !integrator.accept_step)) - if integrator.opts.verbose - if isdefined(integrator, :EEst) - EEst = ", and step error estimate = $(integrator.EEst)" - else - EEst = "" + # Bail out if we take a step with dt less than the minimum value (which may be time dependent) + # except if we are sucessfully taking such a small timestep is to hit a tstop exactly + # We also exit if the ODE is unstable (by default this is the same as nonfinite u) + # but only consider the ODE unstable if the error is somewhat controlled + # (to prevent from bailing out as unstable when we just took way too big a step) + if !opts.force_dtmin && opts.adaptive + if abs(integrator.dt) <= abs(opts.dtmin) && + (((hasproperty(integrator, :opts) && hasproperty(opts, :tstops)) ? + integrator.t + integrator.dt < integrator.tdir * first(opts.tstops) : + true) || (hasproperty(integrator, :accept_step) && !integrator.accept_step)) + if verbose + if isdefined(integrator, :EEst) + EEst = ", and step error estimate = $(integrator.EEst)" + else + EEst = "" + end + @warn("dt($(integrator.dt)) <= dtmin($(opts.dtmin)) at t=$(integrator.t)$EEst. Aborting. There is either an error in your model specification or the true solution is unstable.") + end + return ReturnCode.DtLessThanMin + elseif abs(integrator.dt) <= abs(eps(integrator.t)) && + hasproperty(integrator, :accept_step) && !integrator.accept_step + if verbose + if isdefined(integrator, :EEst) + EEst = ", and step error estimate = $(integrator.EEst)" + else + EEst = "" + end + @warn("dt($(integrator.dt)) <= eps(t)($(integrator.t)) $EEst. Aborting. There is either an error in your model specification or the true solution is unstable (or the true solution can not be represented in the precision of $(eltype(integrator.u))).") end - @warn("dt($(integrator.dt)) <= dtmin($(integrator.opts.dtmin)) at t=$(integrator.t)$EEst. Aborting. There is either an error in your model specification or the true solution is unstable.") + return ReturnCode.Unstable end - return ReturnCode.DtLessThanMin - end - if integrator.opts.unstable_check(integrator.dt, integrator.u, integrator.p, - integrator.t) - bigtol = max(integrator.opts.reltol, integrator.opts.abstol) - # only declare instability if the dt is very small - # or we have at least one digit of accuracy in the solution - if integrator.dt Date: Thu, 16 May 2024 11:52:08 -0400 Subject: [PATCH 2/7] typo --- src/integrator_interface.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/integrator_interface.jl b/src/integrator_interface.jl index f41cae8f0..2dad8bfb3 100644 --- a/src/integrator_interface.jl +++ b/src/integrator_interface.jl @@ -603,7 +603,7 @@ function check_error(integrator::DEIntegrator) # The last part: # Bail out if we take a step with dt less than the minimum value (which may be time dependent) - # except if we are sucessfully taking such a small timestep is to hit a tstop exactly + # except if we are successfully taking such a small timestep is to hit a tstop exactly # We also exit if the ODE is unstable (by default this is the same as nonfinite u) # but only consider the ODE unstable if the error is somewhat controlled # (to prevent from bailing out as unstable when we just took way too big a step) From 8c09f210f99f78c425b7cdddb98d9b2770e94ae2 Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Thu, 16 May 2024 13:18:15 -0400 Subject: [PATCH 3/7] Update src/integrator_interface.jl Co-authored-by: Christopher Rackauckas --- src/integrator_interface.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/integrator_interface.jl b/src/integrator_interface.jl index 2dad8bfb3..c9361682b 100644 --- a/src/integrator_interface.jl +++ b/src/integrator_interface.jl @@ -634,7 +634,7 @@ function check_error(integrator::DEIntegrator) return ReturnCode.Unstable end end - bigtol = max(opts.reltol, opts.abstol) + bigtol = max(max(opts.reltol), max(opts.abstol)) if isdefined(integrator, :EEst) && integrator.EEst * bigtol < .1 if opts.unstable_check(integrator.dt, integrator.u, integrator.p, integrator.t) From d2961351a8980ecc5dd455217207af11e7314fde Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Mon, 20 May 2024 17:45:59 -0400 Subject: [PATCH 4/7] bugfix --- src/integrator_interface.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/integrator_interface.jl b/src/integrator_interface.jl index c9361682b..c52ff0bf9 100644 --- a/src/integrator_interface.jl +++ b/src/integrator_interface.jl @@ -634,7 +634,7 @@ function check_error(integrator::DEIntegrator) return ReturnCode.Unstable end end - bigtol = max(max(opts.reltol), max(opts.abstol)) + bigtol = max(maximum(opts.reltol), maximum(opts.abstol)) if isdefined(integrator, :EEst) && integrator.EEst * bigtol < .1 if opts.unstable_check(integrator.dt, integrator.u, integrator.p, integrator.t) From 1326803048543d1c09f3bd7290ea4f1503ab0771 Mon Sep 17 00:00:00 2001 From: oscarddssmith Date: Tue, 21 May 2024 14:48:53 -0400 Subject: [PATCH 5/7] I think this now works as intended --- src/integrator_interface.jl | 26 +++++++++++--------------- 1 file changed, 11 insertions(+), 15 deletions(-) diff --git a/src/integrator_interface.jl b/src/integrator_interface.jl index c52ff0bf9..d32f8661e 100644 --- a/src/integrator_interface.jl +++ b/src/integrator_interface.jl @@ -604,14 +604,15 @@ function check_error(integrator::DEIntegrator) # The last part: # Bail out if we take a step with dt less than the minimum value (which may be time dependent) # except if we are successfully taking such a small timestep is to hit a tstop exactly - # We also exit if the ODE is unstable (by default this is the same as nonfinite u) - # but only consider the ODE unstable if the error is somewhat controlled - # (to prevent from bailing out as unstable when we just took way too big a step) + # We also exit if the ODE is unstable acording to a user chosen callbakc + # but only if we accepted the step to prevent from bailing out as unstable + # when we just took way too big a step) + step_accepted = !hasproperty(integrator, :accept_step) || integrator.accept_step if !opts.force_dtmin && opts.adaptive if abs(integrator.dt) <= abs(opts.dtmin) && - (((hasproperty(integrator, :opts) && hasproperty(opts, :tstops)) ? + (!step_accepted || ((hasproperty(integrator, :opts) && hasproperty(opts, :tstops)) ? integrator.t + integrator.dt < integrator.tdir * first(opts.tstops) : - true) || (hasproperty(integrator, :accept_step) && !integrator.accept_step)) + true)) if verbose if isdefined(integrator, :EEst) EEst = ", and step error estimate = $(integrator.EEst)" @@ -621,8 +622,7 @@ function check_error(integrator::DEIntegrator) @warn("dt($(integrator.dt)) <= dtmin($(opts.dtmin)) at t=$(integrator.t)$EEst. Aborting. There is either an error in your model specification or the true solution is unstable.") end return ReturnCode.DtLessThanMin - elseif abs(integrator.dt) <= abs(eps(integrator.t)) && - hasproperty(integrator, :accept_step) && !integrator.accept_step + elseif !step_accepted && integrator.t isa AbstractFloat && abs(integrator.dt) <= abs(eps(integrator.t)) if verbose if isdefined(integrator, :EEst) EEst = ", and step error estimate = $(integrator.EEst)" @@ -634,15 +634,11 @@ function check_error(integrator::DEIntegrator) return ReturnCode.Unstable end end - bigtol = max(maximum(opts.reltol), maximum(opts.abstol)) - if isdefined(integrator, :EEst) && integrator.EEst * bigtol < .1 - if opts.unstable_check(integrator.dt, integrator.u, integrator.p, - integrator.t) - if verbose - @warn("Instability detected. Aborting") - end - return ReturnCode.Unstable + if step_accepted && opts.unstable_check(integrator.dt, integrator.u, integrator.p, integrator.t) + if verbose + @warn("Instability detected. Aborting") end + return ReturnCode.Unstable end if last_step_failed(integrator) if verbose From bb3a749fbc06946867fa96fabd07bf1439b6deeb Mon Sep 17 00:00:00 2001 From: oscarddssmith Date: Tue, 21 May 2024 14:50:04 -0400 Subject: [PATCH 6/7] typo --- src/integrator_interface.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/integrator_interface.jl b/src/integrator_interface.jl index d32f8661e..b0281bdda 100644 --- a/src/integrator_interface.jl +++ b/src/integrator_interface.jl @@ -604,7 +604,7 @@ function check_error(integrator::DEIntegrator) # The last part: # Bail out if we take a step with dt less than the minimum value (which may be time dependent) # except if we are successfully taking such a small timestep is to hit a tstop exactly - # We also exit if the ODE is unstable acording to a user chosen callbakc + # We also exit if the ODE is unstable according to a user chosen callback # but only if we accepted the step to prevent from bailing out as unstable # when we just took way too big a step) step_accepted = !hasproperty(integrator, :accept_step) || integrator.accept_step From 761a8c7b72da029ac643e62eb3e4b23120d6b2bb Mon Sep 17 00:00:00 2001 From: oscarddssmith Date: Wed, 22 May 2024 10:47:38 -0400 Subject: [PATCH 7/7] remove unnecessary check --- src/integrator_interface.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/integrator_interface.jl b/src/integrator_interface.jl index b0281bdda..9ca7c5ddf 100644 --- a/src/integrator_interface.jl +++ b/src/integrator_interface.jl @@ -610,7 +610,7 @@ function check_error(integrator::DEIntegrator) step_accepted = !hasproperty(integrator, :accept_step) || integrator.accept_step if !opts.force_dtmin && opts.adaptive if abs(integrator.dt) <= abs(opts.dtmin) && - (!step_accepted || ((hasproperty(integrator, :opts) && hasproperty(opts, :tstops)) ? + (!step_accepted || (hasproperty(opts, :tstops) ? integrator.t + integrator.dt < integrator.tdir * first(opts.tstops) : true)) if verbose