SciML · ChrisRackauckas-Claude · Jul 29, 2026 · Jul 29, 2026
diff --git a/docs/Project.toml b/docs/Project.toml
@@ -80,7 +80,7 @@ StochasticDiffEqLevyArea = {path = "../lib/StochasticDiffEqLevyArea"}
 
 [compat]
 DiffEqBase = "7.6.1"
-DiffEqDevTools = "3, 4.0"
+DiffEqDevTools = "4"
 Documenter = "0.27, 1"
 GlobalDiffEq = "1.2.1"
 ImplicitDiscreteSolve = "2"

diff --git a/docs/src/devtools/alg_dev/convergence.md b/docs/src/devtools/alg_dev/convergence.md
@@ -20,7 +20,10 @@ A type which holds the data from a convergence simulation.
 
 ### Fields
 
-  - `solutions::Array{<:DESolution}`: Holds all the PdeSolutions.
+  - `solutions::Array{<:DESolution}`: Holds the solutions. For the Monte-Carlo method of
+    `test_convergence` this holds one representative trajectory per convergence-axis
+    point rather than every trajectory, since the error statistics have already been
+    computed from the full ensemble; pass `retain_solutions = true` to keep them all.
 
   - `errors`: Dictionary of the error calculations. Can contain:
 

diff --git a/lib/DelayDiffEq/Project.toml b/lib/DelayDiffEq/Project.toml
@@ -64,7 +64,7 @@ DDEProblemLibrary = "0.1"
 BinaryHeaps = "1"
 DiffEqBase = "7"
 DiffEqCallbacks = "4.17.0"
-DiffEqDevTools = "3"
+DiffEqDevTools = "4"
 DiffEqNoiseProcess = "5.30.0"
 FastBroadcast = "1.3"
 FiniteDiff = "2.27"

diff --git a/lib/DiffEqDevTools/Project.toml b/lib/DiffEqDevTools/Project.toml
@@ -1,7 +1,7 @@
 name = "DiffEqDevTools"
 uuid = "f3b72e0c-5b89-59e1-b016-84e28bfd966d"
 authors = ["Chris Rackauckas <accounts@chrisrackauckas.com>"]
-version = "3.1.5"
+version = "4.0.0"
 
 [deps]
 CommonSolve = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2"

diff --git a/lib/DiffEqDevTools/src/convergence.jl b/lib/DiffEqDevTools/src/convergence.jl
@@ -54,6 +54,26 @@ function ConvergenceSimulation(
     return (ConvergenceSimulation(solutions, errors, N, auxdata, 𝒪est, convergence_axis))
 end
 
+"""
+    drop_trajectories(sim::EnsembleTestSolution) -> EnsembleTestSolution
+
+Discard every trajectory but the first, keeping the error statistics computed from the
+full ensemble.
+
+`calculate_ensemble_errors` has already reduced the ensemble to the error dictionaries
+that a [`ConvergenceSimulation`](@ref) reads, so for a Monte-Carlo study the retained
+trajectories are dead weight — an SROCK ensemble at 3e6 trajectories holds tens of GiB
+of them. One trajectory is kept rather than none because the `ConvergenceSimulation`
+constructor inspects `solutions[1].u[1]` for its element type.
+"""
+function drop_trajectories(sim::SciMLBase.EnsembleTestSolution{T, N}) where {T, N}
+    u = sim.u[1:min(1, length(sim.u))]
+    return SciMLBase.EnsembleTestSolution{T, N, typeof(u)}(
+        u, sim.errors, sim.weak_errors, sim.error_means, sim.error_medians,
+        sim.elapsedTime, sim.converged
+    )
+end
+
 """
     test_convergence(dts, prob, alg[, ensemblealg]; kwargs...)
     test_convergence(probs, convergence_axis, alg; kwargs...)
@@ -66,6 +86,19 @@ is passed to the solver as `dt`. Remaining keyword arguments are forwarded to `s
 The computed solutions must contain error measurements, usually obtained from an
 analytic solution. Use [`analyticless_test_convergence`](@ref) when only a numerical
 reference solution is available.
+
+`retain_solutions` controls whether the trajectories survive in the returned
+`ConvergenceSimulation`. It defaults to `false`: each step size's ensemble is reduced to
+its error summary as soon as it is solved and then stripped to one representative
+trajectory, so the memory is bounded by a single step size's ensemble rather than by the
+whole study. A Monte-Carlo study at 5e5 trajectories over 8 step sizes is otherwise
+tens of GiB of solutions kept only to compute a mean.
+
+The error statistics are computed from the full ensemble either way, so `errors`,
+`weak_errors`, `error_means` and `𝒪est` are unaffected — the difference is only that
+`sim.solutions[i].u` holds one trajectory rather than `trajectories` of them. Pass
+`retain_solutions = true` when the trajectories themselves are the point, for instance
+to compare two algorithms path by path.
 """
 function test_convergence(
         dts::AbstractArray,
@@ -77,7 +110,7 @@ function test_convergence(
         trajectories, save_start = true, save_everystep = true,
         timeseries_errors = save_everystep, adaptive = false,
         weak_timeseries_errors = false, weak_dense_errors = false,
-        expected_value = nothing, kwargs...
+        expected_value = nothing, retain_solutions = false, kwargs...
     )
     N = length(dts)
 
@@ -87,6 +120,10 @@ function test_convergence(
         ensemble_prob = EnsembleProblem(prob)
     end
 
+    # Reducing each step size's ensemble as soon as it is solved keeps only one of them
+    # alive at a time; deferring to the loop below would hold every step size at once.
+    reduce_eagerly = !retain_solutions && expected_value === nothing
+
     _solutions = Array{Any}(undef, length(dts))
     for i in 1:length(dts)
         sol = solve(
@@ -98,20 +135,30 @@ function test_convergence(
             kwargs...
         )
         @info "dt: $(dts[i]) ($i/$N)"
+        if reduce_eagerly
+            sol = drop_trajectories(
+                SciMLBase.calculate_ensemble_errors(
+                    sol;
+                    weak_timeseries_errors = weak_timeseries_errors,
+                    weak_dense_errors = weak_dense_errors
+                )
+            )
+        end
         _solutions[i] = sol
     end
 
     auxdata = Dict("dts" => dts)
 
     if expected_value == nothing
-        solutions = [
-            SciMLBase.calculate_ensemble_errors(
+        solutions = reduce_eagerly ? _solutions :
+            [
+                SciMLBase.calculate_ensemble_errors(
                     sim;
                     weak_timeseries_errors = weak_timeseries_errors,
                     weak_dense_errors = weak_dense_errors
                 )
                 for sim in _solutions
-        ]
+            ]
         # Now Calculate Weak Errors
         additional_errors = Dict()
         for k in keys(solutions[1].weak_errors)

diff --git a/lib/DiffEqDevTools/test/qa/Project.toml b/lib/DiffEqDevTools/test/qa/Project.toml
@@ -6,6 +6,6 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
 Aqua = "0.8.11"
-DiffEqDevTools = "3"
+DiffEqDevTools = "4"
 SciMLTesting = "2.1"
 julia = "1.10"
diff --git a/lib/DiffEqDevTools/test/retain_solutions_tests.jl b/lib/DiffEqDevTools/test/retain_solutions_tests.jl
@@ -0,0 +1,87 @@
+using DiffEqDevTools, StochasticDiffEq, Random, Test
+
+# The Monte-Carlo `test_convergence` summarises by default so that large studies do not
+# hold every trajectory of every step size at once. That must change only what is kept,
+# never what is measured, so compare the statistics against `retain_solutions = true`
+# exactly. It also pins the default, which is the breaking half of the change.
+
+linear_analytic(u0, p, t, W) = @.(u0 * exp(0.63155t + 0.87W))
+f_linear(du, u, p, t) = @.(du = 1.01 * u)
+σ_linear(du, u, p, t) = @.(du = 0.87 * u)
+prob = SDEProblem(
+    SDEFunction(f_linear, σ_linear, analytic = linear_analytic), [1 / 2], (0.0, 1.0)
+)
+dts = 1 .// 2 .^ (7:-1:4)
+trajectories = 400
+
+function run(retain)
+    Random.seed!(100)
+    return test_convergence(
+        dts, prob, EM(), EnsembleSerial(); save_everystep = false,
+        trajectories = trajectories, weak_timeseries_errors = false,
+        retain_solutions = retain
+    )
+end
+
+kept = run(true)
+dropped = run(false)
+defaulted = let
+    Random.seed!(100)
+    test_convergence(
+        dts, prob, EM(), EnsembleSerial(); save_everystep = false,
+        trajectories = trajectories, weak_timeseries_errors = false
+    )
+end
+
+@testset "retain_solutions preserves the statistics" begin
+    @test keys(kept.𝒪est) == keys(dropped.𝒪est)
+    for k in keys(kept.𝒪est)
+        @test kept.𝒪est[k] == dropped.𝒪est[k]
+    end
+    @test keys(kept.errors) == keys(dropped.errors)
+    for k in keys(kept.errors)
+        @test kept.errors[k] == dropped.errors[k]
+    end
+    for i in eachindex(dts)
+        @test kept.solutions[i].weak_errors == dropped.solutions[i].weak_errors
+        @test kept.solutions[i].error_means == dropped.solutions[i].error_means
+        @test kept.solutions[i].errors == dropped.solutions[i].errors
+    end
+end
+
+@testset "summarising is the default" begin
+    for i in eachindex(dts)
+        @test length(defaulted.solutions[i].u) == 1
+        @test defaulted.solutions[i].errors == kept.solutions[i].errors
+    end
+    for k in keys(kept.𝒪est)
+        @test defaulted.𝒪est[k] == kept.𝒪est[k]
+    end
+end
+
+@testset "retain_solutions drops the trajectories" begin
+    for i in eachindex(dts)
+        @test length(kept.solutions[i].u) == trajectories
+        @test length(dropped.solutions[i].u) == 1
+    end
+    @test Base.summarysize(dropped) < Base.summarysize(kept) / 10
+end
+
+@testset "retain_solutions is ignored when expected_value is given" begin
+    # The expected_value path takes a reducing output_func and never builds an
+    # EnsembleTestSolution, so there is nothing to drop; it must still run and agree.
+    ev = 0.5 * exp(1.01)
+    ensemble_prob = EnsembleProblem(prob; output_func = (sol, ctx) -> (sol.u[end][1], false))
+    Random.seed!(100)
+    a = test_convergence(
+        dts, ensemble_prob, EM(), EnsembleSerial(); save_everystep = false,
+        trajectories = trajectories, expected_value = ev, retain_solutions = false
+    )
+    Random.seed!(100)
+    b = test_convergence(
+        dts, ensemble_prob, EM(), EnsembleSerial(); save_everystep = false,
+        trajectories = trajectories, expected_value = ev
+    )
+    @test a.𝒪est[:weak_final] == b.𝒪est[:weak_final]
+    @test length(a.solutions[1].u) == trajectories
+end
diff --git a/lib/DiffEqDevTools/test/runtests.jl b/lib/DiffEqDevTools/test/runtests.jl
@@ -27,6 +27,9 @@ if TEST_GROUP == "Core" || TEST_GROUP == "ALL"
     @time @testset "Analyticless Convergence Tests" begin
         include("analyticless_convergence_tests.jl")
     end
+    @time @testset "Retain Solutions Tests" begin
+        include("retain_solutions_tests.jl")
+    end
     @time @testset "Analyticless Stochastic WP" begin
         include("analyticless_stochastic_wp.jl")
     end

diff --git a/lib/OrdinaryDiffEqAdamsBashforthMoulton/Project.toml b/lib/OrdinaryDiffEqAdamsBashforthMoulton/Project.toml
@@ -34,7 +34,7 @@ Pkg = "1"
 Test = "<0.0.1, 1"
 FastBroadcast = "1.3"
 Random = "<0.0.1, 1"
-DiffEqDevTools = "3"
+DiffEqDevTools = "4"
 MuladdMacro = "0.2.4"
 SciMLBase = "3.39"
 OrdinaryDiffEqCore = "4"

diff --git a/lib/OrdinaryDiffEqAdamsBashforthMoulton/test/threaded/Project.toml b/lib/OrdinaryDiffEqAdamsBashforthMoulton/test/threaded/Project.toml
@@ -11,7 +11,7 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 OrdinaryDiffEqCore = {path = "../../../OrdinaryDiffEqCore"}
 
 [compat]
-DiffEqDevTools = "3"
+DiffEqDevTools = "4"
 FastBroadcast = "1.3"
 ODEProblemLibrary = "1"
 OrdinaryDiffEqAdamsBashforthMoulton = "2"

diff --git a/lib/OrdinaryDiffEqBDF/Project.toml b/lib/OrdinaryDiffEqBDF/Project.toml
@@ -53,7 +53,7 @@ ForwardDiff = "1.3.3"
 Test = "<0.0.1, 1"
 FastBroadcast = "1.3"
 Random = "<0.0.1, 1"
-DiffEqDevTools = "3"
+DiffEqDevTools = "4"
 DifferentiationInterface = "0.7.18"
 MuladdMacro = "0.2.4"
 LinearSolve = "5.1"

diff --git a/lib/OrdinaryDiffEqCore/Project.toml b/lib/OrdinaryDiffEqCore/Project.toml
@@ -58,7 +58,7 @@ ArrayInterface = "7.28"
 BinaryHeaps = "1"
 CommonSolve = "0.2.6"
 DiffEqBase = "7.8"
-DiffEqDevTools = "3"
+DiffEqDevTools = "4"
 DifferentiationInterface = "0.7.18"
 DocStringExtensions = "0.9.5"
 EnumX = "1.0.4"

diff --git a/lib/OrdinaryDiffEqDefault/Project.toml b/lib/OrdinaryDiffEqDefault/Project.toml
@@ -47,7 +47,7 @@ Pkg = "1"
 OrdinaryDiffEqTsit5 = "2"
 Test = "<0.0.1, 1"
 Random = "<0.0.1, 1"
-DiffEqDevTools = "3"
+DiffEqDevTools = "4"
 OrdinaryDiffEqBDF = "2"
 OrdinaryDiffEqVerner = "2"
 LinearSolve = "5.1"

diff --git a/lib/OrdinaryDiffEqDifferentiation/Project.toml b/lib/OrdinaryDiffEqDifferentiation/Project.toml
@@ -17,7 +17,7 @@ ForwardDiff = "1.3.3"
 FunctionWrappersWrappers = "1"
 FastBroadcast = "1.3"
 Random = "<0.0.1, 1"
-DiffEqDevTools = "3"
+DiffEqDevTools = "4"
 Test = "<0.0.1, 1"
 FiniteDiff = "2.27"
 DifferentiationInterface = "0.7.18"

diff --git a/lib/OrdinaryDiffEqExplicitRK/Project.toml b/lib/OrdinaryDiffEqExplicitRK/Project.toml
@@ -32,7 +32,7 @@ Pkg = "1"
 Test = "<0.0.1, 1"
 FastBroadcast = "1.3"
 Random = "<0.0.1, 1"
-DiffEqDevTools = "3"
+DiffEqDevTools = "4"
 MuladdMacro = "0.2.4"
 LinearAlgebra = "1.10"
 TruncatedStacktraces = "1.4"

diff --git a/lib/OrdinaryDiffEqExplicitTableaus/Project.toml b/lib/OrdinaryDiffEqExplicitTableaus/Project.toml
@@ -27,7 +27,7 @@ OrdinaryDiffEqLowOrderRK = {path = "../OrdinaryDiffEqLowOrderRK"}
 [compat]
 SciMLTesting = "2.1"
 DiffEqBase = "7"
-DiffEqDevTools = "3"
+DiffEqDevTools = "4"
 ODEProblemLibrary = "1"
 OrdinaryDiffEq = "7"
 OrdinaryDiffEqExplicitRK = "2"

diff --git a/lib/OrdinaryDiffEqExponentialRK/Project.toml b/lib/OrdinaryDiffEqExponentialRK/Project.toml
@@ -45,7 +45,7 @@ OrdinaryDiffEqTsit5 = "2"
 Test = "<0.0.1, 1"
 FastBroadcast = "1.3"
 Random = "<0.0.1, 1"
-DiffEqDevTools = "3"
+DiffEqDevTools = "4"
 MuladdMacro = "0.2.4"
 OrdinaryDiffEqVerner = "2"
 LinearSolve = "5.1"

diff --git a/lib/OrdinaryDiffEqExtrapolation/Project.toml b/lib/OrdinaryDiffEqExtrapolation/Project.toml
@@ -38,7 +38,7 @@ Pkg = "1"
 Test = "<0.0.1, 1"
 FastBroadcast = "1.3"
 Random = "<0.0.1, 1"
-DiffEqDevTools = "3"
+DiffEqDevTools = "4"
 MuladdMacro = "0.2.4"
 LinearSolve = "5.1"
 OrdinaryDiffEqDifferentiation = "3"

diff --git a/lib/OrdinaryDiffEqFIRK/Project.toml b/lib/OrdinaryDiffEqFIRK/Project.toml
@@ -41,7 +41,7 @@ Pkg = "1"
 Test = "<0.0.1, 1"
 FastBroadcast = "1.3"
 Random = "<0.0.1, 1"
-DiffEqDevTools = "3"
+DiffEqDevTools = "4"
 FastGaussQuadrature = "1.2"
 MuladdMacro = "0.2.4"
 LinearSolve = "5.1"

diff --git a/lib/OrdinaryDiffEqFeagin/Project.toml b/lib/OrdinaryDiffEqFeagin/Project.toml
@@ -31,7 +31,7 @@ Pkg = "1"
 Test = "<0.0.1, 1"
 FastBroadcast = "1.3"
 Random = "<0.0.1, 1"
-DiffEqDevTools = "3"
+DiffEqDevTools = "4"
 MuladdMacro = "0.2.4"
 SciMLBase = "3.39"
 OrdinaryDiffEqCore = "4"

diff --git a/lib/OrdinaryDiffEqFunctionMap/Project.toml b/lib/OrdinaryDiffEqFunctionMap/Project.toml
@@ -34,7 +34,7 @@ Pkg = "1"
 Test = "<0.0.1, 1"
 FastBroadcast = "1.3"
 Random = "<0.0.1, 1"
-DiffEqDevTools = "3"
+DiffEqDevTools = "4"
 MuladdMacro = "0.2.4"
 SciMLBase = "3.39"
 ODEProblemLibrary = "1"

diff --git a/lib/OrdinaryDiffEqHighOrderRK/Project.toml b/lib/OrdinaryDiffEqHighOrderRK/Project.toml
@@ -31,7 +31,7 @@ Pkg = "1"
 Test = "<0.0.1, 1"
 FastBroadcast = "1.3"
 Random = "<0.0.1, 1"
-DiffEqDevTools = "3"
+DiffEqDevTools = "4"
 MuladdMacro = "0.2.4"
 SciMLBase = "3.39"
 OrdinaryDiffEqCore = "4.10"

diff --git a/lib/OrdinaryDiffEqIMEXMultistep/Project.toml b/lib/OrdinaryDiffEqIMEXMultistep/Project.toml
@@ -31,7 +31,7 @@ Pkg = "1"
 Test = "<0.0.1, 1"
 FastBroadcast = "1.3"
 Random = "<0.0.1, 1"
-DiffEqDevTools = "3"
+DiffEqDevTools = "4"
 OrdinaryDiffEqDifferentiation = "3"
 SciMLBase = "3.39"
 OrdinaryDiffEqCore = "4"

diff --git a/lib/OrdinaryDiffEqImplicitTableaus/Project.toml b/lib/OrdinaryDiffEqImplicitTableaus/Project.toml
@@ -21,7 +21,7 @@ path = "../DiffEqDevTools"
 [compat]
 SciMLTesting = "2.1"
 DiffEqBase = "7"
-DiffEqDevTools = "3"
+DiffEqDevTools = "4"
 Pkg = "1"
 Test = "<0.0.1, 1"
 julia = "1.10"

diff --git a/lib/OrdinaryDiffEqLinear/Project.toml b/lib/OrdinaryDiffEqLinear/Project.toml
@@ -35,7 +35,7 @@ Pkg = "1"
 OrdinaryDiffEqTsit5 = "2"
 Test = "<0.0.1, 1"
 Random = "<0.0.1, 1"
-DiffEqDevTools = "3"
+DiffEqDevTools = "4"
 OrdinaryDiffEqVerner = "2"
 ExponentialUtilities = "1.33"
 LinearAlgebra = "1.10"