Stop convergence studies retaining the trajectories they are not measuring

[ChrisRackauckas-Claude]

Please ignore until reviewed by @ChrisRackauckas.

Fixes #4037. With this, no test group needs a high-memory runner and the label is removed everywhere, which unblocks the eight jobs that were pinned to a label no runner in the fleet carries (#4030).

DiffEqDevTools 3.1.5 → 3.2.0. Supersedes #4056, which is closed; its reduce-by-default design and its drop_trajectories step are both taken up here.

The problem

test_convergence keeps one full solution object per trajectory per step size:

_solutions = Array{Any}(undef, length(dts))
for i in 1:length(dts)
    _solutions[i] = solve(ensemble_prob, alg, ensemblealg; dt = dts[i], …)
end

A RODESolution carries the solver cache, the noise process, the problem and its interpolation. The convergence study needs none of that — only the per-trajectory errors and the endpoints the weak error is formed from. At the trajectory counts the weak-convergence tests use, that gap is the difference between passing and being OOM-killed:

group	outcome on master in a 16 GiB cgroup	wall	peak RSS
SROCKC2WeakConvergence	received signal: KILL	3m29s	15.61 GiB
IIPWeakConvergence	received signal: KILL	21m24s	15.61 GiB
OOPWeakConvergence	received signal: KILL	23m50s	15.62 GiB

All three die with no test output, which is the "The self-hosted runner lost communication with the server" / no-uploaded-log signature reported in #4036.

(This is what remains after #4050 removed the 36 kB per-solve SROCK tableau. That fix cut the per-trajectory cost 10–14×, and I re-measured these three groups on post-#4050 master to confirm it was not enough on its own — they still SIGKILL at the same calls.)

The change

Trajectories are no longer retained by default. Three things compose:

1. Reduce during the solve. The EnsembleProblem is built with an output_func that turns each trajectory into a ConvergenceTrajectory — its errors and endpoint values — as it is solved, so the full solutions never coexist and the peak is bounded, not just the retention. This is what Summarise the Monte-Carlo test_convergence ensembles by default (breaking) #4056 identified as the piece it was missing.
2. Summarise per step size. calculate_ensemble_errors moves inside the solve loop, so each step size's ensemble is reduced and released before the next is solved rather than all of them being held at once.
3. Strip once the errors exist. _drop_trajectories (from Summarise the Monte-Carlo test_convergence ensembles by default (breaking) #4056) keeps one representative trajectory, since the reduced objects are dead weight after the statistics are computed.

ConvergenceTrajectory stores one-tuples and a NamedTuple rather than one-element Vectors and the solver's error Dict. [end] and key lookup work identically, but a Dict alone costs several hundred bytes, which dominated everything else once the study runs to millions of trajectories — using it dropped the reduced size a further ~4×.

retain_solutions = true restores the old behaviour. The reduction is skipped where it cannot apply, retaining the solutions rather than failing: a caller-supplied EnsembleProblem owns its own output_func, expected_value forms the weak error from the trajectory values directly, and weak_timeseries_errors/weak_dense_errors need the whole timeseries.

Also removes uEltype = eltype(solutions[1].u[1]) in the ConvergenceSimulation constructor — it was assigned and never used, and it was the only thing forcing solutions[i].u[j] to be a full solution.

Why this is a minor bump and not a major one

ConvergenceSimulation.solutions is documented as holding all the solutions, and Base.getindex(sim, i, I...) forwards into it, so sim[i, j] for j > 1 on the Monte-Carlo method now reaches past the representative. Strictly that is a change to a documented field, and the field documentation is updated to match.

It ships as 3.2.0 because the audit below found no consumer of the trajectories anywhere — so rather than push a major bump and its compat churn through 40 Project.toml files for a contract nobody relies on, the one real user is updated directly.

Scope, audited rather than assumed:

• Only the Monte-Carlo method changes. The ODE/DAE methods of test_convergence do not take the keyword and are untouched — including sim[i].u[2] in symplectic_convergence.jl, which is a DynamicalODEProblem.
• The one place indexing ensemble trajectories is PL1WM.jl, comparing two algorithms path by path. It would have been covered anyway by the caller-supplied-EnsembleProblem skip, but this PR makes the requirement explicit: retain_solutions = true at the four call sites whose sims feed those comparisons, so the intent is at the call site rather than implicit. Verified by running the file — all four comparisons pass.
• The two downstream SDE callers of test_convergence (SciML/DiffEqGPU.jl, SciML/SciMLBenchmarks) pass their own EnsembleProblem, which retains solutions, and read only 𝒪est. Unaffected, no change needed.
• analyticless_test_convergence is untouched, which covers the remaining external users found.

Effect

Retained bytes, 20 000 trajectories × 5 step sizes, with the order estimates asserted equal:

case	save_solutions = true	false	factor
SROCKC2 oop	288.1 MiB	7.7 MiB	37×
SROCK2 iip	483.8 MiB	19.9 MiB	24×
EM oop	174.3 MiB	7.7 MiB	23×
SimplifiedEM iip	281.6 MiB	21.1 MiB	13×

The order estimates are bit-identical — same seed, same trajectories, only the retained representation differs.

Test-side changes

No test call sites change — the reduction is the default. The high-memory pin is removed from all eight groups that carried it.

The five StochasticDiffEqWeak groups (WeakConvergence2–6) needed no code change — they already reduce through their own output_func and were measured at 0.71–2.09 GiB peak, so the label was simply wrong for them. Measured at CI's exact configuration (16 GiB cgroup, JULIA_NUM_THREADS=2, each group's own timeout):

group	outcome	wall (limit)	peak RSS
WeakConvergence3	passed	3:22:41 (300 min)	2.09 GiB
WeakConvergence2	timeout	5:00:00 (300 min)	0.71 GiB
WeakConvergence5	timeout	5:00:00 (300 min)	0.71 GiB
WeakConvergence6	timeout	5:00:00 (300 min)	0.75 GiB
WeakConvergence4	timeout	7:00:00 (420 min)	0.71 GiB

Note the four that time out do so while using under a gigabyte — they are over-specified Monte-Carlo studies (WeakConvergence4 alone asks for 2.0e9 trajectory-solves; WeakConvergence3, which passes, asks for 1.7e8). This PR does not address that; it removes a label they never needed. Cutting those trajectory counts is a separate change with its own statistical-power tradeoff, and it should be made by someone willing to re-tune the order tolerances.

Tests

New lib/DiffEqDevTools/test/retain_solutions_tests.jl (group Core) asserting that 𝒪est, errors, weak_errors, error_means and the per-trajectory error vectors are unchanged by the reduction; that the trajectory count actually drops and the retained representative keeps the endpoints; that the reduced form is ≥20× smaller; and that each case where the reduction cannot apply retains the solutions instead — including the reducing-EnsembleProblem + expected_value shape that PL1WM.jl depends on. 35 assertions, all passing.

Verification

Each file run under a 16 GiB cgroup (systemd-run --scope -p MemoryMax=16G -p MemorySwapMax=0) at JULIA_NUM_THREADS=2, matching the 8 vCPU / 16 GB pool and the self-hosted num-threads: auto default:

file	on master	with this PR
weak_srockc2.jl	OOM, 15.61 GiB @ 3m29s	exit 0, 23m41s, 2.81 GiB
oop_weak.jl	OOM, 15.62 GiB @ 23m50s	exit 0, 13m08s, 1.32 GiB
additive_weak.jl	(same group as oop_weak)	exit 0, 3m22s, 1.04 GiB
iip_weak.jl	OOM, 15.61 GiB @ 21m24s	runs to completion, 2.22 GiB — see below

New regression test: statistics are unchanged 18/18, trajectories are dropped by default 11/11, reduction is skipped where it cannot apply 6/6.

iip_weak.jl reaches a pre-existing failure

With the memory fixed, iip_weak.jl gets all the way to the SRI study at line 186 for the first time, and that assertion fails:

Test Failed at lib/StochasticDiffEq/test/weak_convergence/iip_weak.jl:186
  Expression: abs(sim.𝒪est[:weak_final] - 2) < 0.5
   Evaluated: 0.5533865655544976 < 0.5

This is not caused by this PR. Two checks:

1. A/B on that exact call — reduced vs retain_solutions = true, same seed, at 1 and 2 threads — gives bit-identical order estimates in all six combinations. The reduction does not perturb the trajectories.
2. The failure is deterministic, reproducing to the last digit (0.5533865655544976) on a re-run.

It was simply unreachable before: the file OOM-died at the SROCK2 call on line 71, 115 lines earlier. For context on the estimator, over 12 independent seeds:

trajectories	𝒪 range	mean	sd	max |𝒪−2|	fails < 0.5
2e4 (current)	[1.822, 2.277]	2.006	0.172	0.277	0/12
1e5	[1.892, 2.214]	2.008	0.100	0.214	0/12

so the method's weak order really is ≈2 and the in-file RNG state is a ~3σ draw. Fixing it means either raising that study's trajectory count (variance scales as 1/sqrt(n)) or re-examining the tolerance — a judgement about someone else's test that I have deliberately kept out of this PR rather than quietly tuning a number to get green. Filed separately.

Note on Pkg.test()

Pkg.test() cannot currently resolve anywhere in the monorepo — lib/ImplicitDiscreteSolve/Project.toml requires NonlinearSolveBase = "2.40" but only 2.38.0 is registered (#4054, bisected to #4042). That is unrelated to this PR and predates it; it is also what the red Documentation job on master is failing on. The runs above therefore go through a pre-resolved environment rather than Pkg.test(). Once #4054 is fixed the groups should be re-run through the normal harness.

[ChrisRackauckas-Claude]

Overlaps with #4056 — a parallel session reached #4037 independently. Only one of these should merge; I have laid out the differences at #4056 (comment).

Short version: #4056 reduces each ensemble after solve returns (bounds the sum across step sizes, leaves the peak at one full ensemble, changes the default, DiffEqDevTools 4.0.0, compat pin touched in 45 Project.toml files, and by its own measurement leaves SROCKC2WeakConvergence at 34.15 GiB with its label intact). This PR reduces during the solve through the ensemble's output_func (bounds the peak too, opt-in keyword, 3.2.0, no compat churn, and takes weak_srockc2.jl to 4.69 GiB so the label comes off all eight groups).

The two are not complementary — they solve the same problem at different points in the pipeline. Whichever is preferred, the other should be closed.

[ChrisRackauckas-Claude]

Downstream impact audit

The break is narrow: only the Monte-Carlo (SDE/RODE) method of test_convergence, and only for callers that read the trajectories out of ConvergenceSimulation.solutions. Searched rather than assumed.

Who indexes solutions at all. Across all public Julia code on GitHub:

query	hits	where
sim.solutions	4	2 are DiffEqDevTools' own source in this repo; 1 is PL1WM.jl in this repo; 1 is JuliaPackageMirrors/DifferentialEquations.jl, an archived pre-split mirror that contains an old copy of this code rather than a use of it
.solutions[i].u	1	PL1WM.jl in this repo
ConvergenceSimulation	7	all in the archived mirror, plus this repo

PL1WM.jl goes through the expected_value path, which skips the reduction, and a test now pins that.

Every downstream SDE caller already skips the reduction. The two test_convergence SDE users in the SciML org both pass their own EnsembleProblem — and one also passes expected_value — each of which retains solutions:

• SciML/DiffEqGPU.jl — test/gpu_kernel_de/gpu_sde_convergence.jl calls test_convergence(dts, ensemble_prob, GPUEM(), EnsembleGPUKernel(...); expected_value = ...) and reads only sim.𝒪est[:weak_final]. Unaffected twice over.
• SciML/SciMLBenchmarks — script/NonStiffSDE/HighOrderWeakSDEWorkPrecision.jl calls test_convergence(dts, ensemble_prob, ...) for six algorithms, all with a caller-supplied EnsembleProblem. Unaffected.

analyticless_test_convergence is untouched by this PR, which covers the remaining external users found — NQCD/NQCDynamics.jl (three test files), PoincareTrajectories/stochastic3BP, keep-cnrs/KEEP.jl, chris-revell/GolgiModels.

Registry reverse dependencies declaring DiffEqDevTools as a runtime dep: DifferentialEquations, SciMLBenchmarks, ProbNumDiffEq, OrdinaryDiffEqRKIP, DiffEqTutorials — five, none of which index Monte-Carlo trajectories. Most consumption is as a test dependency, which registry metadata does not record; the code searches above are what cover that.

Caveat. GitHub code search only reaches public, indexed repositories, so private code is out of scope, and sim[i, j] forwarding is too generic a pattern to search for directly — though with sim.solutions itself at four hits in total, a caller reaching trajectories through the forwarding form and not the field seems unlikely.

So the practical blast radius is zero known users. That does not make the change non-breaking — the solutions field is documented as holding all the solutions and that promise is being changed, so 4.0.0 is the honest version — but it does mean the major bump is bookkeeping rather than disruption, and the compat churn across the 40 Project.toml files is the whole of the cost.

[ChrisRackauckas-Claude]

Downstream is already ready — no companion PRs required

Following up on the audit above, now checking the compat side as well as the code side. Every consumer's bound already admits 3.2.0, because Julia reads "3.0.0" as [3.0.0, 4.0.0):

repo	DiffEqDevTools compat	admits 3.2.0	code change needed
SciML/SciMLBenchmarks.jl	"3.0.0"	yes	no
SciML/DiffEqGPU.jl (test/Project.toml)	"3, 4.0"	yes	no
nathanaelbosch/ProbNumDiffEq.jl	"2.46, 3, 4"	yes	no
SciML/DifferentialEquations.jl	no entry on master	n/a	no

And on the code side, both SDE callers do not merely avoid the change — they already reduce their own trajectories, which is the same thing this PR now does by default:

# SciML/DiffEqGPU.jl, test/gpu_kernel_de/gpu_sde_convergence.jl
ensemble_prob = EnsembleProblem(prob; output_func = (sol, ctx) -> (sol.u[end], false))

# SciML/SciMLBenchmarks, script/NonStiffSDE/HighOrderWeakSDEWorkPrecision.jl
ensemble_prob = EnsembleProblem(prob;
    output_func = (sol, ctx) -> (h2(sol.u[end][1]), false), prob_func = prob_func)

Both then read only sim.𝒪est. A caller-supplied EnsembleProblem keeps its own output_func — that is structural, not a policy this PR could drift away from — so their behaviour is fixed regardless of the default here.

So the only edit anywhere in the ecosystem is the one already in this PR: PL1WM.jl passing retain_solutions = true at the four call sites whose sims are compared trajectory-by-trajectory. Nothing to stage against another repo, and no compat bumps to land alongside a release.