Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Make error check more compatible with SciML interface #896

Draft
wants to merge 2 commits into
base: master
Choose a base branch
from
Draft

Make error check more compatible with SciML interface #896

Changes from 1 commit
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
3611e4a
Make error check more compatible with SciML interface
termi-official Dec 22, 2024
53c19b8
Update src/integrator_interface.jl
ChrisRackauckas Dec 23, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
76 changes: 47 additions & 29 deletions src/integrator_interface.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -168,6 +168,13 @@ function set_proposed_dt!(i::DEIntegrator, dt)
error("set_proposed_dt!: method has not been implemented for the integrator")
end

"""
get_tdir(i::DEIntegrator)

Get the time direction of the integrator. Should return 1 or -1 with the same type as the time type of the integrator.
"""
get_tdir(i::DEIntegrator) = i.tdir
ChrisRackauckas marked this conversation as resolved.
Show resolved Hide resolved

"""
savevalues!(integrator::DEIntegrator,
force_save=false) -> Tuple{Bool, Bool}
Expand Down Expand Up @@ -406,6 +413,15 @@ function get_sol(integrator::DEIntegrator)
return integrator.sol
end

"""
get_retcode(integrator::DEIntegrator)

Get the return code of the integrator.
"""
function get_retcode(integrator::DEIntegrator)
return integrator.sol.retcode
end
Comment on lines +426 to +428
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

To quickly elaborate on this, I essentially have problems where "time sub-integrators" advance parts of the solution and I could not find a sane way to define this distributed solution stuff into the current solution interface. Technically I do not even need every subintegrator to hold a solution. They just need to communicate return codes between each other. Having this interface here would be a compromise.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Most integrators don't set the return code until the end though?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Correct. But for operator splitting methods I have essentially a tree of integrators, each solving subproblems on small time intervals. The easiest example is essentially to have an additivley split right hand side for some ODE:

$d_t u = f(u,p,t) = f_1(u,p,t) + f_2(u,p,t)$

which I want to integrate from [0,T]. Operator splitting method now definine two ODEs

$d_t u^* = f_1(u^*,p,t)$

and

$d_t u^{} = f_2(u^{},p,t)$

and they intrgrate each of these in a specific order on subintervals on $[0,T]$. The simplest rule is to solve them alternatingly on intervals with fixed dt. E.g. the first step would be

  1. Set $u^*(0) = u_0$
  2. Solve $d_t u^* = f_1(u^*,p,t)$ on $[0,dt]$ (with your algorithm of choice)
  3. Set $u^{**}(0) = u^*(dt)$
  4. Solve $d_t u^{} = f_2(u^{},p,t)$ on $[0,dt]$ (with your algorithm of choice)

In my implementation right now I have three separate integrators. One custom time integrator for the coordination of which subintegrator needs to integrate its associated problem next (+ time step controll of the time intervals to solve the subproblems on), and another two integrators for the respective subproblems. To check whether a solve worked or not the outer time integrator essentially checks the inner integrators return code after each solve of the subproblems. Does this explain it?

Probably not the best architecture, but at least it is functional.


### Addat isn't a real thing. Let's make it a real thing Gretchen

function addat!(a::AbstractArray, idxs, val = nothing)
Expand Down Expand Up @@ -564,18 +580,27 @@ end

last_step_failed(integrator::DEIntegrator) = false

# Standard error message for classical PID type controllers
function controller_message_on_dtmin_error(integrator::DEIntegrator)
if isdefined(integrator, :EEst)
return ", and step error estimate = $(integrator.EEst)"
else
return ""
end
end
Comment on lines +589 to +595
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Some voodoo time adaption algorithms actually work without defining an error estimate, so I would like to remove the dependency that the integrator itself needs to carry the error estimate and would like to propose in the future that there is either some kind of ControllerCache or the contoller itself carries the estimate.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yeah... that's going to be rough because there needs to be a tie-in in every stepper for how to calculate and store EEst, if it's needed. I'm not opposed though, because yes some methods like a priori estimates don't need the EEst and so it shouldn't spend the time to calculate it.


"""
check_error(integrator)

Check state of `integrator` and return one of the
[Return Codes](https://docs.sciml.ai/DiffEqDocs/stable/basics/solution/#retcodes)
"""
function check_error(integrator::DEIntegrator)
if integrator.sol.retcode ∉ (ReturnCode.Success, ReturnCode.Default)
return integrator.sol.retcode
if get_retcode(integrator) ∉ (ReturnCode.Success, ReturnCode.Default)
return get_retcode(integrator)
end
opts = integrator.opts
verbose = opts.verbose
verbose = hasproperty(opts, :verbose) && opts.verbose
# This implementation is intended to be used for ODEIntegrator and
# SDEIntegrator.
if isnan(integrator.dt)
Expand All @@ -584,7 +609,7 @@ function check_error(integrator::DEIntegrator)
end
return ReturnCode.DtNaN
end
if integrator.iter > opts.maxiters
if hasproperty(opts, :maxiters) && 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
Expand All @@ -597,42 +622,35 @@ function check_error(integrator::DEIntegrator)
# 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
if !opts.force_dtmin && opts.adaptive
if abs(integrator.dt) <= abs(opts.dtmin) &&
(!step_accepted || (hasproperty(opts, :tstops) ?
integrator.t + integrator.dt < integrator.tdir * first(opts.tstops) :
true))
step_rejected = last_step_failed(integrator)
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Right now this does not work, as last_step_failed is always false if the integrator is adaptive. If everyone agrees that this is closer to the semantics that we really want, then I will go through all subpackages and make the corresponding change.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think this could be reasonable.

step_accepted = !step_rejected # For better readability
force_dtmin = hasproperty(integrator, :force_dtmin) && integrator.force_dtmin
if !force_dtmin && isadaptive(integrator)
dt_below_min = abs(integrator.dt) ≤ abs(opts.dtmin)
before_next_tstop = has_tstop(integrator) ? integrator.t + integrator.dt < get_tdir(integrator) * first_tstop(integrator) : true
if dt_below_min && (step_rejected || before_next_tstop)
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.")
controller_string = controller_message_on_dtmin_error(integrator)
@warn("dt($(integrator.dt)) <= dtmin($(opts.dtmin)) at t=$(integrator.t)$(controller_string). Aborting. There is either an error in your model specification or the true solution is unstable.")
end
return ReturnCode.DtLessThanMin
elseif !step_accepted && integrator.t isa AbstractFloat &&
abs(integrator.dt) <= abs(eps(integrator.t))
elseif step_rejected && integrator.t isa AbstractFloat &&
abs(integrator.dt) <= abs(eps(integrator.t)) # = DiffEqBase.timedepentdtmin(integrator)
if verbose
if isdefined(integrator, :EEst)
EEst = ", and step error estimate = $(integrator.EEst)"
else
EEst = ""
end
@warn("At t=$(integrator.t), dt was forced below floating point epsilon $(integrator.dt)$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))).")
controller_string = controller_message_on_dtmin_error(integrator)
@warn("At t=$(integrator.t), dt was forced below floating point epsilon $(integrator.dt)$(controller_string). 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
return ReturnCode.Unstable
end
end
if step_accepted &&
if step_accepted && hasproperty(opts, :unstable_check) &&
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 last_step_failed(integrator) && !isadaptive(integrator)
if verbose
@warn("Newton steps could not converge and algorithm is not adaptive. Use a lower dt.")
end
Expand Down Expand Up @@ -873,7 +891,7 @@ Base.length(iter::TimeChoiceIterator) = length(iter.ts)
end
end

xflip --> integrator.tdir < 0
xflip --> get_tdir(integrator) < 0

if denseplot
seriestype --> :path
Expand Down Expand Up @@ -904,11 +922,11 @@ Base.length(iter::TimeChoiceIterator) = length(iter.ts)
end

function step!(integ::DEIntegrator, dt, stop_at_tdt = false)
(dt * integ.tdir) < 0 * oneunit(dt) && error("Cannot step backward.")
(dt * get_tdir(integ)) < 0 * oneunit(dt) && error("Cannot step backward.")
t = integ.t
next_t = t + dt
stop_at_tdt && add_tstop!(integ, next_t)
while integ.t * integ.tdir < next_t * integ.tdir
while integ.t * get_tdir(integ) < next_t * get_tdir(integ)
step!(integ)
integ.sol.retcode in (ReturnCode.Default, ReturnCode.Success) || break
end
Expand Down