Report the top error-estimate contributors on solver exit

lib/OrdinaryDiffEqCore/Project.toml

@@ -28,6 +28,7 @@ MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
 MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 Preferences = "21216c6a-2e73-6563-6e65-726566657250"
+Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
@@ -82,11 +83,12 @@ Pkg = "1"
 Polyester = "0.7"
 PrecompileTools = "1.2.1, 1.3"
 Preferences = "1.5.0"
+Printf = "1.9"
 Random = "<0.0.1, 1"
 RecursiveArrayTools = "4.2.0"
 Reexport = "1.2.2"
 SafeTestsets = "0.1.0"
-SciMLBase = "3.40"
+SciMLBase = "3.41"
 SciMLLogging = "1.10.1, 2"
 SciMLOperators = "1.24.3"
 SciMLStructures = "1.7"

lib/OrdinaryDiffEqCore/src/OrdinaryDiffEqCore.jl

@@ -50,6 +50,7 @@ import SciMLOperators: MatrixOperator, FunctionOperator,
     isconstant
 
 import Random
+import Printf: @sprintf
 
 import RecursiveArrayTools: recursivecopy!, recursivecopy, recursive_bottom_eltype, recursive_unitless_bottom_eltype, recursive_unitless_eltype, copyat_or_push!, DiffEqArray
 

lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl

@@ -707,6 +707,88 @@ function log_step!(progress_name, progress_id, progress_message, dt, u, p, t, ts
     )
 end
 
+"""
+    error_estimate_residuals(cache)
+
+Return the per-component weighted residuals `atmp` behind the scalar error
+estimate, or `nothing` when the cache does not retain them.
+
+Out-of-place caches build the residuals in a local and drop them, so only
+mutable caches can report a breakdown after the fact.
+"""
+error_estimate_residuals(cache) = hasfield(typeof(cache), :atmp) ? getfield(cache, :atmp) : nothing
+error_estimate_residuals(cache::CompositeCache) =
+    error_estimate_residuals(@inbounds cache.caches[cache.current])
+function error_estimate_residuals(cache::DefaultCache)
+    current = cache.current
+    # Constituent caches are initialized lazily, so an inactive slot can be undefined.
+    for (i, name) in enumerate((:cache1, :cache2, :cache3, :cache4, :cache5, :cache6))
+        if i == current
+            return isdefined(cache, name) ?
+                error_estimate_residuals(getfield(cache, name)) : nothing
+        end
+    end
+    return nothing
+end
+
+# Number of worst-offending state components reported in the exit diagnostic.
+const N_ERROR_CONTRIBUTORS = 3
+
+_diagnostic_value(x) = _format_diagnostic_value(SciMLBase.value(x))
+_format_diagnostic_value(x) = x
+_format_diagnostic_value(x::AbstractFloat) = isfinite(x) ? @sprintf("%.4g", x) : x
+
+# NaN residuals are the most interesting ones, so they sort alongside Inf rather
+# than poisoning the comparisons.
+_contributor_weight(x) = (v = abs(SciMLBase.value(x)); isnan(v) ? typemax(v) : v)
+
+"""
+    top_error_contributors(atmp, n)
+
+Return the indices of the `n` components of `atmp` with the largest magnitude,
+largest first.
+"""
+function top_error_contributors(atmp, n)
+    idxs = collect(eachindex(atmp))
+    n = min(n, length(idxs))
+    partialsort!(idxs, 1:n, by = i -> _contributor_weight(atmp[i]), rev = true)
+    return idxs[1:n]
+end
+
+function SciMLBase.log_error_estimate(integrator::ODEIntegrator)
+    integrator.opts.adaptive || return ""
+    EEst = _diagnostic_value(get_EEst(integrator))
+    stats = integrator.stats
+    lines = [
+        "step error estimate EEst = $EEst (a step is accepted when EEst <= 1) after $(stats.naccept) accepted and $(stats.nreject) rejected steps",
+    ]
+
+    atmp = error_estimate_residuals(integrator.cache)
+    if atmp isa AbstractArray && !isempty(atmp) && eltype(atmp) <: Number
+        nonfinite = count(!isfinite ∘ SciMLBase.value, atmp)
+        nonfinite > 0 && push!(
+            lines,
+            "$nonfinite of $(length(atmp)) weighted residuals are non-finite (NaN/Inf)"
+        )
+        u = integrator.u
+        uprev = integrator.uprev
+        contributors = [
+            "  atmp[$i] = $(_diagnostic_value(atmp[i]))" *
+                (
+                    u isa AbstractArray && eachindex(u) == eachindex(atmp) ?
+                    ", u[$i] = $(_diagnostic_value(u[i])), uprev[$i] = $(_diagnostic_value(uprev[i]))" : ""
+                )
+                for i in top_error_contributors(atmp, N_ERROR_CONTRIBUTORS)
+        ]
+        push!(
+            lines,
+            "largest contributors to EEst = internalnorm(atmp), where atmp is the tolerance-weighted local error per state component:\n" *
+                join(contributors, "\n")
+        )
+    end
+    return "\n\nError estimate diagnostics:\n" * join(lines, "\n\n")
+end
+
 function fixed_t_for_tstop_error!(integrator, ttmp)
     if _get_next_step_tstop(integrator)
         _set_tstop_flag!(integrator, false)

test/InterfaceIII/exit_diagnostics.jl

@@ -0,0 +1,107 @@
+using OrdinaryDiffEqCore, OrdinaryDiffEqTsit5, OrdinaryDiffEqRosenbrock, OrdinaryDiffEqBDF
+using OrdinaryDiffEqCore: top_error_contributors, error_estimate_residuals
+using SciMLBase, StaticArrays, Logging, Test
+import OrdinaryDiffEqCore.SciMLLogging as SciMLLogging
+
+# Component 7 is the only stiff one, so it alone drives the explicit method's
+# step size down and must be named as the dominant error contributor.
+function stiff_component!(du, u, p, t)
+    @. du = -0.1 * u
+    du[7] = -1.0e5 * u[7]
+    return nothing
+end
+stiff_prob = ODEProblem(stiff_component!, ones(10), (0.0, 10.0))
+
+# Component 4 blows up in finite time while the rest stay bounded.
+function blowup_component!(du, u, p, t)
+    @. du = u^3
+    du[4] = u[4]^5
+    return nothing
+end
+blowup_prob = ODEProblem(blowup_component!, ones(6), (0.0, 10.0))
+
+function solve_logs(args...; kwargs...)
+    io = IOBuffer()
+    sol = with_logger(SimpleLogger(io)) do
+        solve(args...; kwargs...)
+    end
+    return sol, String(take!(io))
+end
+
+@testset "top_error_contributors" begin
+    @test top_error_contributors([0.1, -5.0, 3.0, 0.0], 2) == [2, 3]
+    @test top_error_contributors([1.0, 2.0], 5) == [2, 1]
+    @test top_error_contributors(Float64[], 3) == Int[]
+    # NaN outranks every finite residual rather than poisoning the comparisons.
+    @test first(top_error_contributors([1.0e9, NaN, 2.0], 3)) == 2
+end
+
+@testset "maxiters names the dominant error contributor" begin
+    sol, logs = solve_logs(stiff_prob, Tsit5(); maxiters = 500)
+    @test sol.retcode == ReturnCode.MaxIters
+    @test occursin("Error estimate diagnostics", logs)
+    @test occursin("step error estimate EEst", logs)
+    @test occursin("atmp[7]", logs)
+    @test occursin("u[7]", logs)
+    # The exit point is reported so the log alone locates the failure in time.
+    @test occursin("Reached t=", logs)
+end
+
+@testset "instability exits report the error estimate too" begin
+    sol, logs = solve_logs(blowup_prob, Tsit5())
+    @test sol.retcode in (ReturnCode.Unstable, ReturnCode.DtLessThanMin)
+    @test occursin("atmp[4]", logs)
+
+    sol, logs = solve_logs(blowup_prob, FBDF())
+    @test sol.retcode in (ReturnCode.Unstable, ReturnCode.DtLessThanMin)
+    @test occursin("atmp[4]", logs)
+end
+
+@testset "non-finite residuals are called out" begin
+    function nan_component!(du, u, p, t)
+        @. du = -u
+        du[2] = t > 0.5 ? NaN : -u[2]
+        return nothing
+    end
+    sol, logs = solve_logs(ODEProblem(nan_component!, ones(4), (0.0, 1.0)), Tsit5())
+    @test sol.retcode != ReturnCode.Success
+    @test occursin("of 4 weighted residuals are non-finite", logs)
+    @test occursin("atmp[2] = NaN", logs)
+end
+
+@testset "composite and default caches report the active cache" begin
+    for alg in (AutoTsit5(Rosenbrock23()),)
+        sol, logs = solve_logs(stiff_prob, alg; maxiters = 30)
+        @test sol.retcode == ReturnCode.MaxIters
+        @test occursin("atmp[7]", logs)
+    end
+end
+
+@testset "silent verbosity emits nothing" begin
+    sol, logs = solve_logs(
+        stiff_prob, Tsit5(); maxiters = 500, verbose = SciMLLogging.None()
+    )
+    @test sol.retcode == ReturnCode.MaxIters
+    @test isempty(logs)
+end
+
+@testset "degrades gracefully without a stored residual" begin
+    # Out-of-place caches build the residuals in a local, so only the scalar
+    # estimate is available. The diagnostic must still be emitted.
+    oop = ODEProblem(
+        (u, p, t) -> SA[-0.1u[1], -0.1u[2], -1.0e5 * u[3]], SA[1.0, 1.0, 1.0],
+        (0.0, 10.0)
+    )
+    sol, logs = solve_logs(oop, Tsit5(); maxiters = 500)
+    @test sol.retcode == ReturnCode.MaxIters
+    @test occursin("step error estimate EEst", logs)
+    @test !occursin("atmp[", logs)
+
+    integ = init(oop, Tsit5())
+    @test error_estimate_residuals(integ.cache) === nothing
+end
+
+@testset "non-adaptive integration has no error estimate to report" begin
+    integ = init(stiff_prob, Tsit5(); adaptive = false, dt = 1.0e-4)
+    @test SciMLBase.log_error_estimate(integ) == ""
+end

test/runtests.jl

@@ -66,7 +66,8 @@ function interface_iii()
     @time @safetestset "No Index Tests" include("InterfaceIII/noindex_tests.jl")
     @time @safetestset "Events + DAE addsteps Tests" include("InterfaceIII/event_dae_addsteps.jl")
     @time @safetestset "Units Tests" include("InterfaceIII/units_tests.jl")
-    return @time @safetestset "DEVerbosity Tests" include("InterfaceIII/verbosity.jl")
+    @time @safetestset "DEVerbosity Tests" include("InterfaceIII/verbosity.jl")
+    return @time @safetestset "Solver Exit Diagnostics Tests" include("InterfaceIII/exit_diagnostics.jl")
 end
 
 function interface_iv()