Fix ExplicitTaylor fsallast aliasing to state buffer - #3975
Open
singhharsh1708 wants to merge 1 commit into
Open
Fix ExplicitTaylor fsallast aliasing to state buffer#3975singhharsh1708 wants to merge 1 commit into
singhharsh1708 wants to merge 1 commit into
Conversation
singhharsh1708
force-pushed
the
fix-explicittaylor-fsallast-aliasing
branch
3 times, most recently
from
August 1, 2026 21:26
72e7246 to
eb143a7
Compare
get_fsalfirstlast returned (cache.u, cache.u) for ExplicitTaylorCache and ExplicitTaylorAdaptiveOrderCache, aliasing the FSAL buffer to the state buffer itself. ExplicitTaylor does not override isfsal, so the auto-dt heuristic takes its FSAL branch and writes an extra f evaluation into fsallast; with the aliasing that call is f(u, u, p, t), which corrupts the state mid-evaluation for any in-place RHS that reads u after writing du. The corruption lands before the first step: integrator.u is already wrong after init, and the resulting NaN in the initial-step estimate collapses dt0 to dtmin, so the controller burns hundreds of steps recovering. On Pleiades at 1e-10 that is 905 steps versus 587 after the fix. Return independent buffers of the du type instead, matching the core OrdinaryDiffEqConstantCache default.
singhharsh1708
force-pushed
the
fix-explicittaylor-fsallast-aliasing
branch
from
August 1, 2026 22:10
eb143a7 to
651ad20
Compare
Comment on lines
+84
to
+91
| # FSAL currently not used. `cache.u` must NOT be reused here: the generic | ||
| # auto-dt heuristic (`ode_determine_initdt`) writes an extra `f` evaluation | ||
| # into whichever buffer `fsallast` points to, so aliasing it to the state | ||
| # itself makes that call `f(u, u, p, t)` and can corrupt `u` mid-evaluation | ||
| # for RHS functions that read from `u` after writing to `du`, which corrupts the | ||
| # state before the first step. Use independent buffers of the `du` type, matching | ||
| # the core `OrdinaryDiffEqConstantCache` default. | ||
| get_fsalfirstlast(cache::ExplicitTaylorCache, u) = (zero(u), zero(u)) |
Member
There was a problem hiding this comment.
where is this actually called? This function shouldn't allocate.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
get_fsalfirstlastforExplicitTaylorCacheandExplicitTaylorAdaptiveOrderCachereturned(cache.u, cache.u), aliasing the FSAL buffer to the state buffer itself (cache.u === integrator.u).ExplicitTaylordoes not overrideisfsal, so it inheritsisfsal(alg) == trueand the generic auto-dt heuristic (ode_determine_initdt) takes its FSAL branch, writing an extrafevaluation intointegrator.fsallast. With the aliasing, that call is literallyf(u, u, p, t), which corruptsumid-evaluation for any in-place RHS that reads fromuafter writing todu— a normal, valid pattern.Effects on current master:
integrator.uis silently wrong immediately afterinit(...), before a single step is taken. Withu0 = [1.0, 2.0]and an RHS that writesduthen readsu,integrator.ucomes back as[-16.0, 16.0].NaN, emits a spurious "First function call produced NaNs" warning, and falls back totdir * dtmin. Onprob_ode_pleiadesatabstol = reltol = 1e-10the first step size is1.0e-323, and the controller spends hundreds of steps climbing out.Measured on
prob_ode_pleiadeswithExplicitTaylor(order = Val(6)):abstol = reltol = 1e-10Solutions agree to about six significant figures on this chaotic problem.
Fix: return independent buffers of the
dutype rather than the state, matching the core defaultget_fsalfirstlast(cache::OrdinaryDiffEqConstantCache, u) = (zero(u), zero(u)). These were the only two caches in the repo with this aliasing; the correspondingConstantCaches already used the core default.Known follow-up, not fixed here
ExplicitTaylorAdaptiveOrder()had the same aliasing. Before this fix its auto-dt collapsed todtminand the solve died loudly withretcode = MaxItersafter 1,000,001 steps. After the fix that failure is gone, but a separate pre-existing bug in its per-step order/error controller becomes visible: on Pleiades it returnsretcode = Successafter 3 steps withsol.u[end]bit-identical tou0, i.e. the state never moved. That is a silent wrong answer where there used to be a loud failure. It is filed as #3976 rather than fixed here, and it is worth weighing that trade-off when reviewing.Tests
fsallast must not alias the state bufferasserts directly thatintegrator.fsallast !== integrator.uand thatinit(...)leavesuunmodified, for both affected caches, using a two-line RHS that readsuafter writingdu. It fails four times on master and passes here, in about two seconds.Auto-dt on Pleiades does not collapse to dtminkeeps the end-to-end check with a bound that actually discriminates (length(sol.t) < 700; master takes 905 steps, this branch 587).Existing
OrdinaryDiffEqTaylorSeriesCore group passes unchanged. Runic--checkclean on both changed files.AI Disclosure
Claude assisted with this work.