DelayDiffEq: don't rewrite sol.t[end] in shift_past_discontinuity! (#3437)

PR #3427 made `shift_past_discontinuity!` bump `ode_integrator.sol.t[end]`
by one ULP along with `integrator.t`, on the theory that keeping the
invariant `integrator.t == ode_integrator.sol.t[end]` would let the
equality check in `advance_ode_integrator!` stay strict. The commit
message claimed bit-identical equivalence to the pre-#3427 behavior, but
that verification only checked `sol.u` and `sol(t)` queries, not
`sol.t` equality between the outer DDE solution and the inner ODE
sub-integrator solution.

That invariant is the wrong thing to enforce. The step completed at the
pre-shift time and that is what should appear in the saved history; the
shift is an integration-state nudge to cross a propagated discontinuity,
not a retroactive rewrite of the step that just finished.

Rewriting `ode.sol.t[end]` breaks `lib/DelayDiffEq/test/interface/saveat.jl`
lines 18/46/75 (outer DDE `sol.t == dde_int.integrator.sol.t` for the
default every-step `saveat`). Rewriting the outer DDE `sol.t[end]` as
well (which would restore that equality) breaks the `#190` testset in
`discontinuities.jl` (integer discontinuity times are expected to appear
exactly in user-visible `sol.t`), and is also semantically wrong for
sparse user `saveat` where `sol.t[end]` is an older user-requested time.

Revert to the pre-#3427 approach: shift only `integrator.t` and the ODE
sub-integrator's `t`, leave both `sol.t[end]` values untouched, and
restore the `abs(t - sol.t[end]) > eps(t)` tolerance in
`advance_ode_integrator!`. Both `saveat.jl` and `discontinuities.jl`
pass (verified locally along with Interface, Integrators, Regression
test groups), and the `constant-lag shift past discontinuity (#3426)`
regression test added in #3427 continues to pass since the shift still
crosses the discontinuity — only the saved history is left alone.

Co-authored-by: ChrisRackauckas-Claude <accounts@chrisrackauckas.com>
Co-authored-by: Claude Opus 4.6 (1M context) <noreply@anthropic.com>

Author: 3 people

lib/DelayDiffEq/src/integrators/utils.jl

@@ -51,8 +51,10 @@ function advance_ode_integrator!(integrator::DDEIntegrator, always_calc_begin =
     ode_integrator = integrator.integrator
 
     # algorithm only works if current time of DDE integrator equals final time point
-    # of solution
-    t != ode_integrator.sol.t[end] && error("cannot advance ODE integrator")
+    # of solution (allow one-ULP gap from shift_past_discontinuity!, which nudges
+    # integrator.t past a propagated discontinuity without rewriting sol.t[end])
+    abs(t - ode_integrator.sol.t[end]) > eps(t) &&
+        error("cannot advance ODE integrator")
 
     # complete interpolation data of DDE integrator for time interval [t, t+dt]
     # and copy it to ODE integrator
@@ -231,19 +233,18 @@ function OrdinaryDiffEqCore.handle_discontinuities!(integrator::DDEIntegrator)
     return nothing
 end
 
-# Override shift_past_discontinuity! for DDEIntegrator: shift the DDE's
-# `integrator.t`, the ODE sub-integrator's `t`, and the last saved endpoint
-# `sol.t[end]` together. `u` is continuous across the propagated derivative-
-# discontinuity so only times move. Keeping the invariant
-# `integrator.t == ode_integrator.sol.t[end]` lets the strict equality check
-# in `advance_ode_integrator!` stay strict.
+# Override shift_past_discontinuity! for DDEIntegrator: nudge the DDE's
+# `integrator.t` and the ODE sub-integrator's `t` by one ULP to cross a
+# propagated derivative-discontinuity. `sol.t[end]` on both sides is left alone:
+# the step completed at the pre-shift time and that is what should appear in
+# the saved history (see `#190` testset and `saveat.jl:18/46/75`). The one-ULP
+# gap between `integrator.t` and `ode_integrator.sol.t[end]` is tolerated by
+# the relaxed equality check in `advance_ode_integrator!`.
 function OrdinaryDiffEqCore.shift_past_discontinuity!(integrator::DDEIntegrator)
     integrator.t isa AbstractFloat || return nothing
     newt = integrator.tdir > 0 ? nextfloat(integrator.t) : prevfloat(integrator.t)
     integrator.t = newt
-    ode_integrator = integrator.integrator
-    ode_integrator.t = newt
-    ode_integrator.sol.t[end] = newt
+    integrator.integrator.t = newt
     return nothing
 end