diff --git a/lib/OrdinaryDiffEqCore/Project.toml b/lib/OrdinaryDiffEqCore/Project.toml index 851e267100..8c039ae6a4 100644 --- a/lib/OrdinaryDiffEqCore/Project.toml +++ b/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" diff --git a/lib/OrdinaryDiffEqCore/src/OrdinaryDiffEqCore.jl b/lib/OrdinaryDiffEqCore/src/OrdinaryDiffEqCore.jl index 2d1d032042..a56cb4ba79 100644 --- a/lib/OrdinaryDiffEqCore/src/OrdinaryDiffEqCore.jl +++ b/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 diff --git a/lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl b/lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl index a17ab531eb..b0a266b13f 100644 --- a/lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl +++ b/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) diff --git a/test/InterfaceIII/exit_diagnostics.jl b/test/InterfaceIII/exit_diagnostics.jl new file mode 100644 index 0000000000..687187624c --- /dev/null +++ b/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 diff --git a/test/runtests.jl b/test/runtests.jl index 00fbfd92ca..293f725d7d 100644 --- a/test/runtests.jl +++ b/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()