Skip to content

Commit cd3468b

Browse files
add logging trace to symbolic system
1 parent 053f020 commit cd3468b

1 file changed

Lines changed: 28 additions & 10 deletions

File tree

lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl

Lines changed: 28 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -666,7 +666,7 @@ function SciMLBase.log_instability(integrator::ODEIntegrator)
666666
filter!(t -> !isfinite(t[3]) || abs(t[3]) >= cutoff, entries) #only keep those vals within 1e10 of max or inf/nan
667667
sort!(entries, by = t -> (!isfinite(t[3]), abs(t[3])), rev = true)
668668

669-
# derive rows and columns implicated by the surviving entries
669+
# derive rows and columns from remaining entries
670670
row_set = Set{Int}()
671671
col_set = Set{Int}()
672672
for (i, j, _) in entries
@@ -678,9 +678,15 @@ function SciMLBase.log_instability(integrator::ODEIntegrator)
678678
singularity_cols = sort!(collect(col_set))
679679
end
680680

681+
# trace diagnostics to symbolic system if present
682+
f = integrator.sol.prob.f
683+
sys = (hasproperty(f, :sys) && f.sys !== nothing) ? f.sys : nothing
684+
sym_eqs = (sys !== nothing && hasfield(typeof(sys), :eqs)) ? getfield(sys, :eqs) : nothing
685+
sym_vars = (sys !== nothing && hasfield(typeof(sys), :unknowns)) ? getfield(sys, :unknowns) : nothing
686+
681687
# diagnostic message construction
682688
diagnostic = String[]
683-
if !isempty(nan_inf_idxs) #inf/nan vals in u
689+
if !isempty(nan_inf_idxs) #state vars
684690
if u isa AbstractArray
685691
n_nan = length(nan_inf_idxs)
686692
n_total = length(u)
@@ -693,11 +699,10 @@ function SciMLBase.log_instability(integrator::ODEIntegrator)
693699
push!(diagnostic, "u[$i] = $(u[i]) is non-finite (NaN/Inf)")
694700
end
695701
end
696-
push!(diagnostic, "(Note: reflects the solver's internal state at abort, not sol.u[end] which is the last saved step)")
697702
else
698703
push!(diagnostic, "u = $u is non-finite (NaN/Inf), suggesting blow-up or a NaN in the RHS")
699704
end
700-
elseif !isempty(blown_idxs) #big u values
705+
elseif !isempty(blown_idxs)
701706
if u isa AbstractArray
702707
for i in blown_idxs
703708
push!(diagnostic, "u[$i] = $(@sprintf("%.4g", u[i])) has grown >1e6× its initial value")
@@ -706,15 +711,12 @@ function SciMLBase.log_instability(integrator::ODEIntegrator)
706711
push!(diagnostic, "u = $(@sprintf("%.4g", u)) has grown >1e6× its initial value")
707712
end
708713
end
709-
if bad_entries !== nothing && !isempty(bad_entries) #jacobian analysis
714+
715+
if bad_entries !== nothing && !isempty(bad_entries) #Jacobian analysis
710716
has_nonfinite = false
711717
has_large = false
712718
for (_, _, v) in bad_entries
713-
if isfinite(v)
714-
has_large = true
715-
else
716-
has_nonfinite = true
717-
end
719+
isfinite(v) ? (has_large = true) : (has_nonfinite = true)
718720
end
719721
entry_desc = if has_nonfinite && has_large
720722
"non-finite and large"
@@ -723,13 +725,29 @@ function SciMLBase.log_instability(integrator::ODEIntegrator)
723725
else
724726
"unusually large"
725727
end
728+
726729
example_strs = String[]
727730
for (i, j, v) in first(bad_entries, 5)
728731
push!(example_strs, "J[$i,$j] = $(@sprintf("%.4g", v))")
729732
end
730733
push!(diagnostic, "Jacobian row(s) $singularity_rows have $entry_desc entries (e.g. $(join(example_strs, ", "))), suggesting a singularity in those equation(s)")
734+
if sym_eqs !== nothing
735+
for row in singularity_rows
736+
if row <= length(sym_eqs)
737+
push!(diagnostic, " row $row corresponds to equation: $(sym_eqs[row])") #trace rows back to symbolic eqs
738+
end
739+
end
740+
end
741+
# jac cols
731742
if !isempty(singularity_cols)
732743
push!(diagnostic, "Jacobian column(s) $singularity_cols have $entry_desc entries, suggesting those state component(s) are diverging")
744+
if sym_vars !== nothing
745+
for col in singularity_cols
746+
if col <= length(sym_vars)
747+
push!(diagnostic, " col $col corresponds to variable: $(sym_vars[col])") #trace cols back to symbolic vars
748+
end
749+
end
750+
end
733751
end
734752
end
735753
diagnostic = isempty(diagnostic) ? "" : "\nDiagnostics:\n" * join(diagnostic, "\n") * "."

0 commit comments

Comments
 (0)