From 3dc371d22e500deaaba685b5f3cdf45c208c96f4 Mon Sep 17 00:00:00 2001 From: ChrisRackauckas-Claude Date: Wed, 23 Sep 2026 19:51:35 -0400 Subject: [PATCH 1/4] Add __init for DynamicSS/SICNM so init() doesn't fall through to NonlinearSolve DynamicSS and SICNM only had __solve; init(ssprob, DynamicSS(...)) had no matching __init here, so it fell through to NonlinearSolveBase's generic fallback, which treats "no algorithm-specific __init exists" the same as "no algorithm given" and rewrites the call to NonlinearSolve's default-algorithm converter with ::Nothing in place of the algorithm - silently discarding DynamicSS/SICNM. For a plain SteadyStateProblem that converter still finds a NonlinearProblem to solve, so init() "worked" but ignored the requested algorithm; for a SteadyStateProblem whose lowered_problem is an SCCNonlinearProblem, the converter's NonlinearProblem(prob) call hits a FieldError (SCCNonlinearProblem has no u0 field) - the reported CI failure, init(ssprob, DynamicSS(Tsit5())) in SciMLBase.jl's downstream SymbolicIndexingInterface test group. Add SciMLBase.__init(prob::AbstractSteadyStateProblem, alg::DynamicSS/SICNM, ...), mirroring the existing __solve methods (extracted into shared __dynamicss_ode_setup/__sicnm_ode_setup helpers): for a plain problem it returns the live ODE/DAE integrator so init/solve!/step! compose normally; for an SCC lowering there is no single integrator to hand back (blocks solve sequentially), so it eagerly solves via the existing __solve_scc_lowering and returns the finished solution, matching how a stored NonlinearProblem/ LinearProblem lowering is already handled verbatim elsewhere in this file. Being a more specific method than NonlinearSolveBase's generic fallback, this is picked by ordinary dispatch and the buggy conversion path is never reached for these two algorithms - independent of SciML/SciMLBase.jl#1614 (the generic converter itself is fixed separately there and in SciML/NonlinearSolve.jl#1327). Co-Authored-By: Chris Rackauckas Co-Authored-By: Claude Sonnet Agent-Harness: Claude Code 2.1.251 (subagent) Agent-Model: claude-sonnet Agent-Session: https://claude.ai/code/session_01NB3oFTNzGW79UtDt8AyjsM --- src/solve.jl | 167 +++++++++++++++++++++++++++++++++++++++------------ test/scc.jl | 54 +++++++++++++++++ 2 files changed, 182 insertions(+), 39 deletions(-) diff --git a/src/solve.jl b/src/solve.jl index acacbee..d408d4a 100644 --- a/src/solve.jl +++ b/src/solve.jl @@ -161,18 +161,15 @@ function __without_verbose(kwargs) return (; (name => value for (name, value) in pairs(kwargs) if name !== :verbose)...) end -function SciMLBase.__solve( - prob::SciMLBase.AbstractSteadyStateProblem, alg::DynamicSS, - args...; abstol = 1.0e-8, reltol = 1.0e-6, odesolve_kwargs = (;), - save_idxs = nothing, termination_condition = NonlinearSolveBase.NormTerminationMode(infnorm), - alias = SciMLBase.NonlinearAliasSpecifier(), kwargs... - ) - sccsol = __solve_scc_lowering( - prob, alg, args...; abstol, reltol, odesolve_kwargs, - termination_condition, alias, save_idxs, kwargs... +# Shared `DynamicSS` setup for `__solve` and `__init`: builds the ODE problem that +# integrates `prob`'s residual to steady state, along with the termination +# callback that stops the integration (and, for `__solve`, reports convergence). +# Returns `nothing` for an SCC lowering, which has no single ODE trajectory to +# build — callers fall back to `__solve_scc_lowering` for that case instead. +function __dynamicss_ode_setup( + prob::SciMLBase.AbstractSteadyStateProblem, alg::DynamicSS; + abstol, reltol, odesolve_kwargs, termination_condition, alias, kwargs... ) - sccsol !== nothing && return sccsol - tspan = __get_tspan(prob.u0, alg) f = if prob isa SteadyStateProblem @@ -216,19 +213,36 @@ function SciMLBase.__solve( haskey(kwargs, :callback) && (callback = CallbackSet(callback, kwargs[:callback])) haskey(odesolve_kwargs, :callback) && (callback = CallbackSet(callback, odesolve_kwargs[:callback])) - kwargs = pairs(__without_verbose(kwargs)) - # Construct and solve the ODEProblem + run_kwargs = pairs(__without_verbose(kwargs)) odeprob = ODEProblem{isinplace(prob), true}(f, prob.u0, tspan, prob.p) + odealias = SciMLBase.ODEAliasSpecifier(; + alias_p = alias.alias_p, alias_f = alias.alias_f, alias_u0 = alias.alias_u0 + ) + return (; odeprob, tc_cache, abstol, reltol, callback, run_kwargs, odealias) +end + +function SciMLBase.__solve( + prob::SciMLBase.AbstractSteadyStateProblem, alg::DynamicSS, + args...; abstol = 1.0e-8, reltol = 1.0e-6, odesolve_kwargs = (;), + save_idxs = nothing, termination_condition = NonlinearSolveBase.NormTerminationMode(infnorm), + alias = SciMLBase.NonlinearAliasSpecifier(), kwargs... + ) + sccsol = __solve_scc_lowering( + prob, alg, args...; abstol, reltol, odesolve_kwargs, + termination_condition, alias, save_idxs, kwargs... + ) + sccsol !== nothing && return sccsol + + setup = __dynamicss_ode_setup( + prob, alg; abstol, reltol, odesolve_kwargs, termination_condition, alias, kwargs... + ) odesol = solve( - odeprob, alg.alg, args...; abstol, reltol, kwargs..., - odesolve_kwargs..., callback, save_end = true, - alias = SciMLBase.ODEAliasSpecifier(; - alias_p = alias.alias_p, - alias_f = alias.alias_f, alias_u0 = alias.alias_u0 - ) + setup.odeprob, alg.alg, args...; setup.abstol, setup.reltol, + setup.run_kwargs..., odesolve_kwargs..., setup.callback, save_end = true, + alias = setup.odealias ) - resid, u, retcode = __get_result_from_sol(tc_cache, odesol) + resid, u, retcode = __get_result_from_sol(setup.tc_cache, odesol) if save_idxs !== nothing u = u[save_idxs] @@ -241,6 +255,36 @@ function SciMLBase.__solve( ) end +# `init` on a `DynamicSS`-lowered problem hands back the underlying ODE +# integrator (with the steady-state termination callback installed) rather than +# eagerly integrating to steady state, so it composes with `solve!`/`step!` the +# way `init`/`solve!` do for a plain `ODEProblem`. An `SCCNonlinearProblem` +# lowering has no such continuous trajectory — its blocks solve sequentially, not +# through one shared integrator — so there is nothing to "start and pause": it is +# solved eagerly here, exactly as `__solve` does, and the finished solution is +# returned instead of an integrator. +function SciMLBase.__init( + prob::SciMLBase.AbstractSteadyStateProblem, alg::DynamicSS, + args...; abstol = 1.0e-8, reltol = 1.0e-6, odesolve_kwargs = (;), + save_idxs = nothing, termination_condition = NonlinearSolveBase.NormTerminationMode(infnorm), + alias = SciMLBase.NonlinearAliasSpecifier(), kwargs... + ) + sccsol = __solve_scc_lowering( + prob, alg, args...; abstol, reltol, odesolve_kwargs, + termination_condition, alias, save_idxs, kwargs... + ) + sccsol !== nothing && return sccsol + + setup = __dynamicss_ode_setup( + prob, alg; abstol, reltol, odesolve_kwargs, termination_condition, alias, kwargs... + ) + return init( + setup.odeprob, alg.alg, args...; setup.abstol, setup.reltol, + setup.run_kwargs..., odesolve_kwargs..., setup.callback, save_end = true, + alias = setup.odealias + ) +end + # SICNM: Semi-Implicit Continuous Newton Method # Solves 0 = g(y) by integrating the DAE ẏ = z, 0 = J(y)z + g(y) to steady state, # where J is the Jacobian of g. See the SICNM docstring for details and references. @@ -270,19 +314,14 @@ function __sicnm_g_and_jvp!(gval, jvp, g!::G, y, z) where {G} return nothing end -function SciMLBase.__solve( - prob::SciMLBase.AbstractSteadyStateProblem, alg::SICNM, - args...; abstol = 1.0e-8, reltol = 1.0e-6, odesolve_kwargs = (;), - save_idxs = nothing, - termination_condition = NonlinearSolveBase.AbsNormTerminationMode(infnorm), - alias = SciMLBase.NonlinearAliasSpecifier(), kwargs... - ) - sccsol = __solve_scc_lowering( - prob, alg, args...; abstol, reltol, odesolve_kwargs, - termination_condition, alias, save_idxs, kwargs... +# Shared `SICNM` setup for `__solve` and `__init`: builds the extended DAE ODE +# problem whose continuous-Newton flow drives `g(y) = 0`, along with the +# termination callback based on the residual `g`. Mirrors `__dynamicss_ode_setup` +# above; see its docstring for why an SCC lowering is not handled here. +function __sicnm_ode_setup( + prob::SciMLBase.AbstractSteadyStateProblem, alg::SICNM; + abstol, reltol, odesolve_kwargs, termination_condition, kwargs... ) - sccsol !== nothing && return sccsol - prob.u0 isa AbstractVector || throw(ArgumentError("SICNM currently only supports `AbstractVector` initial conditions")) tspan = __get_tspan(prob.u0, alg) @@ -384,17 +423,42 @@ function SciMLBase.__solve( odefun = SciMLBase.ODEFunction{iip, SciMLBase.FullSpecialize}(fext; mass_matrix) odeprob = ODEProblem{iip}(odefun, u0, tspan, p) + run_kwargs = pairs(__without_verbose(kwargs)) + + return (; + odeprob, tc_cache, n, g, gbuf, iip, + ode_abstol, ode_reltol, callback, run_kwargs, + ) +end + +function SciMLBase.__solve( + prob::SciMLBase.AbstractSteadyStateProblem, alg::SICNM, + args...; abstol = 1.0e-8, reltol = 1.0e-6, odesolve_kwargs = (;), + save_idxs = nothing, + termination_condition = NonlinearSolveBase.AbsNormTerminationMode(infnorm), + alias = SciMLBase.NonlinearAliasSpecifier(), kwargs... + ) + sccsol = __solve_scc_lowering( + prob, alg, args...; abstol, reltol, odesolve_kwargs, + termination_condition, alias, save_idxs, kwargs... + ) + sccsol !== nothing && return sccsol + + setup = __sicnm_ode_setup( + prob, alg; abstol, reltol, odesolve_kwargs, termination_condition, kwargs... + ) odesol = solve( - odeprob, alg.alg, args...; abstol = ode_abstol, reltol = ode_reltol, - kwargs..., odesolve_kwargs..., callback, save_end = true + setup.odeprob, alg.alg, args...; abstol = setup.ode_abstol, + reltol = setup.ode_reltol, setup.run_kwargs..., odesolve_kwargs..., + setup.callback, save_end = true ) - u, retcode = __sicnm_result(tc_cache, odesol, n) - resid = if iip - g(gbuf, u) - gbuf + u, retcode = __sicnm_result(setup.tc_cache, odesol, setup.n) + resid = if setup.iip + setup.g(setup.gbuf, u) + setup.gbuf else - g(u) + setup.g(u) end if save_idxs !== nothing @@ -408,6 +472,31 @@ function SciMLBase.__solve( ) end +# See the analogous `DynamicSS` `__init` above for why an SCC lowering is solved +# eagerly here rather than returning a partial integrator. +function SciMLBase.__init( + prob::SciMLBase.AbstractSteadyStateProblem, alg::SICNM, + args...; abstol = 1.0e-8, reltol = 1.0e-6, odesolve_kwargs = (;), + save_idxs = nothing, + termination_condition = NonlinearSolveBase.AbsNormTerminationMode(infnorm), + alias = SciMLBase.NonlinearAliasSpecifier(), kwargs... + ) + sccsol = __solve_scc_lowering( + prob, alg, args...; abstol, reltol, odesolve_kwargs, + termination_condition, alias, save_idxs, kwargs... + ) + sccsol !== nothing && return sccsol + + setup = __sicnm_ode_setup( + prob, alg; abstol, reltol, odesolve_kwargs, termination_condition, kwargs... + ) + return init( + setup.odeprob, alg.alg, args...; abstol = setup.ode_abstol, + reltol = setup.ode_reltol, setup.run_kwargs..., odesolve_kwargs..., + setup.callback, save_end = true + ) +end + function __sicnm_result(tc_cache, odesol, n) u, _, retcode = termination_condition_result( tc_cache, last(odesol.u)[1:n], last(odesol.t), odesol.retcode diff --git a/test/scc.jl b/test/scc.jl index 4780d4f..5622574 100644 --- a/test/scc.jl +++ b/test/scc.jl @@ -435,3 +435,57 @@ end @test sol.prob isa SCCNonlinearProblem @test sol.original isa Tuple{SciMLBase.LinearSolution, NonlinearSolution} end + +# `init` had no `DynamicSS`/`SICNM`-specific dispatch for `AbstractSteadyStateProblem`, +# so it fell through to NonlinearSolve's "no algorithm" default-conversion path, +# which calls `SciMLBase.NonlinearProblem(prob)` and recurses on the result. For a +# plain `SteadyStateProblem` that materializes to an actual `NonlinearProblem`, so +# the recursion terminates; for one whose `lowered_problem` is an +# `SCCNonlinearProblem`, `NonlinearProblem` returns the input unchanged and the +# recursion never terminates (a `FieldError` on unpatched `SciMLBase`, a stack +# overflow once `SciMLBase.NonlinearProblem(::SCCNonlinearProblem)` is an +# identity). `init` on a plain problem returns a live ODE integrator; on an SCC +# lowering there is no single integrator, so it eagerly solves and returns the +# finished solution (see `SciMLBase.__init` in src/solve.jl). +@testset "init on DynamicSS/SICNM does not fall through to NonlinearSolve's default" begin + @testset "plain SteadyStateProblem returns a live integrator" for alg in ( + DynamicSS(Tsit5()), SICNM(Rodas5P()), + ) + prob = SteadyStateProblem((u, p, t) -> 1 .- u, [0.0]) + integ = init(prob, alg; abstol = 1.0e-10, reltol = 1.0e-10) + sol = solve!(integ) + @test successful_retcode(sol) + # `SICNM`'s integrator carries the extended DAE state `[y; z]`, so only + # the first `length(prob.u0)` components are the original residual state. + @test sol.u[end][1:1] ≈ [1.0] atol = 1.0e-6 + end + + @testset "manually-built SCC lowering, alg=$alg" for alg in ( + DynamicSS(Tsit5()), SICNM(Rodas5P()), + ) + sccprob = dynamicss_scc_problem(false, false) + prob = SteadyStateProblem( + (u, p, t) -> 1 .- u, [0.0, 0.0]; lowered_problem = sccprob + ) + sol = init(prob, alg; abstol = 1.0e-10, reltol = 1.0e-10) + @test successful_retcode(sol) + @test sol.u ≈ [1, 2, 1, 2] atol = 1.0e-8 + @test sol.prob === sccprob + end + + @testset "ModelingToolkit SCC decomposition, alg=$(nameof(typeof(alg)))" for alg in ( + DynamicSS(Tsit5()), SICNM(Rodas5P()), + ) + @variables a(t) b(t) x(t) [irreducible = true] + @named model = System( + [D(a) ~ 5 - 3a - b, D(b) ~ 5 - a - 2b, D(x) ~ a + b - x^3], t + ) + sys = mtkcompile(model) + prob = SteadyStateProblem(sys, [a => 0.8, b => 1.8, x => 0.8]) + + sol = init(prob, alg; abstol = 1.0e-10, reltol = 1.0e-10) + @test successful_retcode(sol) + @test sol[[a, b, x]] ≈ [1, 2, cbrt(3)] atol = 1.0e-8 + @test sol.prob isa SCCNonlinearProblem + end +end From c779737153cc8a5f9babafa080b2726b577c0155 Mon Sep 17 00:00:00 2001 From: ChrisRackauckas-Claude Date: Wed, 23 Sep 2026 20:35:10 -0400 Subject: [PATCH 2/4] Update test comment to match SciMLBase's final ArgumentError behavior The comment still described SciMLBase's NonlinearProblem(::SCCNonlinearProblem) conversion as an identity that could recurse forever; it now raises an ArgumentError, and this repo's new __init methods are more specific than the default path anyway, so it is never reached. Co-Authored-By: Chris Rackauckas Co-Authored-By: Claude Sonnet Agent-Harness: Claude Code 2.1.251 (subagent) Agent-Model: claude-sonnet Agent-Session: https://claude.ai/code/session_01NB3oFTNzGW79UtDt8AyjsM --- test/scc.jl | 15 ++++++--------- 1 file changed, 6 insertions(+), 9 deletions(-) diff --git a/test/scc.jl b/test/scc.jl index 5622574..a7f53bc 100644 --- a/test/scc.jl +++ b/test/scc.jl @@ -438,15 +438,12 @@ end # `init` had no `DynamicSS`/`SICNM`-specific dispatch for `AbstractSteadyStateProblem`, # so it fell through to NonlinearSolve's "no algorithm" default-conversion path, -# which calls `SciMLBase.NonlinearProblem(prob)` and recurses on the result. For a -# plain `SteadyStateProblem` that materializes to an actual `NonlinearProblem`, so -# the recursion terminates; for one whose `lowered_problem` is an -# `SCCNonlinearProblem`, `NonlinearProblem` returns the input unchanged and the -# recursion never terminates (a `FieldError` on unpatched `SciMLBase`, a stack -# overflow once `SciMLBase.NonlinearProblem(::SCCNonlinearProblem)` is an -# identity). `init` on a plain problem returns a live ODE integrator; on an SCC -# lowering there is no single integrator, so it eagerly solves and returns the -# finished solution (see `SciMLBase.__init` in src/solve.jl). +# which cannot handle an `SCCNonlinearProblem` lowering (`SciMLBase.NonlinearProblem` +# raises an `ArgumentError` for it). The `__init` methods below are more specific +# than that default path, so it is never reached: on a plain problem `init` returns +# a live ODE integrator, and on an SCC lowering — which has no single integrator to +# hand back — it eagerly solves and returns the finished solution (see +# `SciMLBase.__init` in src/solve.jl). @testset "init on DynamicSS/SICNM does not fall through to NonlinearSolve's default" begin @testset "plain SteadyStateProblem returns a live integrator" for alg in ( DynamicSS(Tsit5()), SICNM(Rodas5P()), From b7e82b4ebdd95abb3f9bde9bbf973564bffc5f91 Mon Sep 17 00:00:00 2001 From: ChrisRackauckas-Claude Date: Wed, 7 Oct 2026 12:55:32 -0400 Subject: [PATCH 3/4] Forward save_idxs in DynamicSS/SICNM init; trim comments Co-Authored-By: Chris Rackauckas Co-Authored-By: Claude Opus 5.5 (1M context) Agent-Harness: Claude Code (subagent of ci-green head) Agent-Model: claude-opus-5-5[1m] Agent-Session: https://claude.ai/code/session_011aK5Xvp9NpdkzKuEwtiMGd --- src/solve.jl | 29 ++++++++++------------------- test/scc.jl | 16 ++++++++-------- 2 files changed, 18 insertions(+), 27 deletions(-) diff --git a/src/solve.jl b/src/solve.jl index d408d4a..e9f91fe 100644 --- a/src/solve.jl +++ b/src/solve.jl @@ -161,11 +161,9 @@ function __without_verbose(kwargs) return (; (name => value for (name, value) in pairs(kwargs) if name !== :verbose)...) end -# Shared `DynamicSS` setup for `__solve` and `__init`: builds the ODE problem that -# integrates `prob`'s residual to steady state, along with the termination -# callback that stops the integration (and, for `__solve`, reports convergence). -# Returns `nothing` for an SCC lowering, which has no single ODE trajectory to -# build — callers fall back to `__solve_scc_lowering` for that case instead. +# Shared `DynamicSS` setup for `__solve` and `__init`: the ODE problem integrating +# `prob`'s residual to steady state and its termination callback. SCC lowerings are +# handled by `__solve_scc_lowering` before this is reached. function __dynamicss_ode_setup( prob::SciMLBase.AbstractSteadyStateProblem, alg::DynamicSS; abstol, reltol, odesolve_kwargs, termination_condition, alias, kwargs... @@ -255,14 +253,9 @@ function SciMLBase.__solve( ) end -# `init` on a `DynamicSS`-lowered problem hands back the underlying ODE -# integrator (with the steady-state termination callback installed) rather than -# eagerly integrating to steady state, so it composes with `solve!`/`step!` the -# way `init`/`solve!` do for a plain `ODEProblem`. An `SCCNonlinearProblem` -# lowering has no such continuous trajectory — its blocks solve sequentially, not -# through one shared integrator — so there is nothing to "start and pause": it is -# solved eagerly here, exactly as `__solve` does, and the finished solution is -# returned instead of an integrator. +# `init` returns the ODE integrator with the steady-state termination callback +# installed. An SCC lowering solves its blocks sequentially with no single +# integrator, so it is solved eagerly and the finished solution is returned. function SciMLBase.__init( prob::SciMLBase.AbstractSteadyStateProblem, alg::DynamicSS, args...; abstol = 1.0e-8, reltol = 1.0e-6, odesolve_kwargs = (;), @@ -281,7 +274,7 @@ function SciMLBase.__init( return init( setup.odeprob, alg.alg, args...; setup.abstol, setup.reltol, setup.run_kwargs..., odesolve_kwargs..., setup.callback, save_end = true, - alias = setup.odealias + alias = setup.odealias, save_idxs ) end @@ -316,8 +309,7 @@ end # Shared `SICNM` setup for `__solve` and `__init`: builds the extended DAE ODE # problem whose continuous-Newton flow drives `g(y) = 0`, along with the -# termination callback based on the residual `g`. Mirrors `__dynamicss_ode_setup` -# above; see its docstring for why an SCC lowering is not handled here. +# termination callback based on the residual `g`. Mirrors `__dynamicss_ode_setup`. function __sicnm_ode_setup( prob::SciMLBase.AbstractSteadyStateProblem, alg::SICNM; abstol, reltol, odesolve_kwargs, termination_condition, kwargs... @@ -472,8 +464,7 @@ function SciMLBase.__solve( ) end -# See the analogous `DynamicSS` `__init` above for why an SCC lowering is solved -# eagerly here rather than returning a partial integrator. +# As for `DynamicSS`, an SCC lowering is solved eagerly. function SciMLBase.__init( prob::SciMLBase.AbstractSteadyStateProblem, alg::SICNM, args...; abstol = 1.0e-8, reltol = 1.0e-6, odesolve_kwargs = (;), @@ -493,7 +484,7 @@ function SciMLBase.__init( return init( setup.odeprob, alg.alg, args...; abstol = setup.ode_abstol, reltol = setup.ode_reltol, setup.run_kwargs..., odesolve_kwargs..., - setup.callback, save_end = true + setup.callback, save_end = true, save_idxs ) end diff --git a/test/scc.jl b/test/scc.jl index a7f53bc..9639045 100644 --- a/test/scc.jl +++ b/test/scc.jl @@ -436,14 +436,9 @@ end @test sol.original isa Tuple{SciMLBase.LinearSolution, NonlinearSolution} end -# `init` had no `DynamicSS`/`SICNM`-specific dispatch for `AbstractSteadyStateProblem`, -# so it fell through to NonlinearSolve's "no algorithm" default-conversion path, -# which cannot handle an `SCCNonlinearProblem` lowering (`SciMLBase.NonlinearProblem` -# raises an `ArgumentError` for it). The `__init` methods below are more specific -# than that default path, so it is never reached: on a plain problem `init` returns -# a live ODE integrator, and on an SCC lowering — which has no single integrator to -# hand back — it eagerly solves and returns the finished solution (see -# `SciMLBase.__init` in src/solve.jl). +# `init` with `DynamicSS`/`SICNM` uses these algorithms rather than NonlinearSolve's +# default: a plain problem gives a live ODE integrator, an SCC lowering is solved +# eagerly and returns the finished solution. @testset "init on DynamicSS/SICNM does not fall through to NonlinearSolve's default" begin @testset "plain SteadyStateProblem returns a live integrator" for alg in ( DynamicSS(Tsit5()), SICNM(Rodas5P()), @@ -455,6 +450,11 @@ end # `SICNM`'s integrator carries the extended DAE state `[y; z]`, so only # the first `length(prob.u0)` components are the original residual state. @test sol.u[end][1:1] ≈ [1.0] atol = 1.0e-6 + + prob2 = SteadyStateProblem((u, p, t) -> [1, 2] .- u, [0.0, 0.0]) + sol2 = solve!(init(prob2, alg; save_idxs = [2], abstol = 1.0e-10, reltol = 1.0e-10)) + @test length(sol2.u[end]) == 1 + @test sol2.u[end] ≈ [2.0] atol = 1.0e-6 end @testset "manually-built SCC lowering, alg=$alg" for alg in ( From aabe095a88153e7a0bfc0a87ad43ce0140f5a139 Mon Sep 17 00:00:00 2001 From: ChrisRackauckas-Claude Date: Thu, 8 Oct 2026 09:57:47 -0400 Subject: [PATCH 4/4] Return steady-state caches from DynamicSS/SICNM init `init(prob, ::DynamicSS/SICNM)` returned the raw ODE integrator, so `solve!` skipped the termination-condition finalization of `solve`: finite-time nonconvergence reported Success and a protective safe termination reported Terminated instead of Unstable. For SCC lowerings `init` returned a finished solution that `solve!` could not accept. `init` now returns a `SteadyStateODECache` (ODE integrator plus the setup holding the termination cache) or a `SteadyStateSCCCache`; `solve!` on either runs the same finalization as `__solve`, and `step!` advances the ODE integrator. Documents the init/solve! contract, including SICNM's extended `[y; z]` integrator state. Co-Authored-By: Chris Rackauckas Co-Authored-By: Claude Opus 5.5 (1M context) Agent-Harness: Claude Code (subagent of ci-green head) Agent-Model: claude-opus-5-5[1m] Agent-Session: https://claude.ai/code/session_011aK5Xvp9NpdkzKuEwtiMGd --- docs/src/index.md | 30 +++++++++++ src/SteadyStateDiffEq.jl | 2 +- src/solve.jl | 104 ++++++++++++++++++++++++++++----------- test/scc.jl | 96 +++++++++++++++++++++++++++--------- 4 files changed, 179 insertions(+), 53 deletions(-) diff --git a/docs/src/index.md b/docs/src/index.md index bae6ad1..b54932c 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -48,6 +48,36 @@ prob = SteadyStateProblem((u, p, t) -> 1 .- u, [0.0]) sol = solve(prob, SICNM(Rodas3d())) ``` +## Initialization and stepping + +`DynamicSS` and `SICNM` also support the `init`/`solve!` interface. `init(prob, alg; +kwargs...)` accepts the same keyword arguments as `solve` and returns a cache, and +`solve!(cache)` returns the same solution as `solve(prob, alg; kwargs...)`: its return +code comes from the termination condition (for example `ReturnCode.Unstable` from a +safe termination mode, or a failure when the time span ends before steady state is +reached), and `save_idxs` selects components of the final state. + +For a problem that is integrated as a whole, `cache.integrator` is the underlying ODE +integrator, and `step!(cache)` advances it. For `SICNM` that integrator solves the +extended continuous-Newton system, so its state is `[y; z]`: the first +`length(prob.u0)` components are the steady-state variables `y` and the rest are the +Newton direction `z`. Its intermediate states are therefore not states of `prob`; +`solve!` returns only `y`. + +When `prob` carries an `SCCNonlinearProblem` lowering, the blocks are solved one after +another and there is no single integrator to step. `solve!` runs that sequential solve. + +```julia +using SciMLBase: SteadyStateProblem, init, solve!, step! +using SteadyStateDiffEq +using OrdinaryDiffEqRosenbrock: Rodas5P + +prob = SteadyStateProblem((u, p, t) -> 1 .- u, [0.0]) +cache = init(prob, SICNM(Rodas5P())) +step!(cache) +sol = solve!(cache) +``` + ## API ```@docs diff --git a/src/SteadyStateDiffEq.jl b/src/SteadyStateDiffEq.jl index e508efc..2d13eb0 100644 --- a/src/SteadyStateDiffEq.jl +++ b/src/SteadyStateDiffEq.jl @@ -10,7 +10,7 @@ using LinearSolve: LinearSolve using SciMLPublic: @public using SciMLBase: SciMLBase, CallbackSet, LinearProblem, NonlinearProblem, ODEProblem, NonlinearSolution, ReturnCode, SteadyStateProblem, SteadyStateSolution, get_du, init, - isinplace, remake, solve, successful_retcode + isinplace, remake, solve, solve!, step!, successful_retcode using SymbolicIndexingInterface: parameter_values const infnorm = Base.Fix2(norm, Inf) diff --git a/src/solve.jl b/src/solve.jl index e9f91fe..461be71 100644 --- a/src/solve.jl +++ b/src/solve.jl @@ -98,16 +98,21 @@ function SciMLBase.solve( return SciMLBase.build_solution(prob, alg, u, resid; retcode, original = sols) end -# A `SteadyStateProblem` that records an `SCCNonlinearProblem` lowering is -# solved block-sequentially in the lowering's ordering instead of one -# monolithic integration. Returns `nothing` when there is no SCC lowering. -function __solve_scc_lowering( - prob, alg, args...; save_idxs = nothing, kwargs... - ) +# The `SCCNonlinearProblem` lowering recorded on a `SteadyStateProblem`, or +# `nothing` when there is none. +function __scc_lowering(prob) lp = prob isa SteadyStateProblem ? prob.lowered_problem : nothing lp === nothing && return nothing lp isa SciMLBase.AbstractSciMLProblem || (lp = lp(prob)) - lp isa SciMLBase.SCCNonlinearProblem || return nothing + return lp isa SciMLBase.SCCNonlinearProblem ? lp : nothing +end + +# A `SteadyStateProblem` that records an `SCCNonlinearProblem` lowering `lp` is +# solved block-sequentially in the lowering's ordering instead of one +# monolithic integration. +function __solve_scc_lowering( + prob, lp, alg, args...; save_idxs = nothing, kwargs... + ) sccsol = solve(lp, alg, args...; kwargs...) save_idxs === nothing && return sccsol return SciMLBase.build_solution( @@ -149,6 +154,39 @@ function SciMLBase.solve(prob::SteadyStateProblem, args...; kwargs...) ) end +# `init(prob, ::DynamicSS/SICNM)` returns one of these caches; `solve!` finishes +# it with the same steady-state finalization as `solve`. +@concrete struct SteadyStateSCCCache + prob + lowered_problem + alg + args + kwargs +end + +function SciMLBase.solve!(cache::SteadyStateSCCCache) + return __solve_scc_lowering( + cache.prob, cache.lowered_problem, cache.alg, cache.args...; cache.kwargs... + ) +end + +@concrete struct SteadyStateODECache + prob + alg + integrator + setup + save_idxs +end + +function SciMLBase.solve!(cache::SteadyStateODECache) + odesol = solve!(cache.integrator) + return __steady_state_solution( + cache.prob, cache.alg, cache.setup, odesol, cache.save_idxs + ) +end + +SciMLBase.step!(cache::SteadyStateODECache, args...) = step!(cache.integrator, args...) + __get_tspan(u0, alg::Union{DynamicSS, SICNM}) = __get_tspan(u0, alg.tspan) __get_tspan(u0, tspan::Tuple) = tspan function __get_tspan(u0, tspan::Number) @@ -225,11 +263,11 @@ function SciMLBase.__solve( save_idxs = nothing, termination_condition = NonlinearSolveBase.NormTerminationMode(infnorm), alias = SciMLBase.NonlinearAliasSpecifier(), kwargs... ) - sccsol = __solve_scc_lowering( - prob, alg, args...; abstol, reltol, odesolve_kwargs, + lp = __scc_lowering(prob) + lp !== nothing && return __solve_scc_lowering( + prob, lp, alg, args...; abstol, reltol, odesolve_kwargs, termination_condition, alias, save_idxs, kwargs... ) - sccsol !== nothing && return sccsol setup = __dynamicss_ode_setup( prob, alg; abstol, reltol, odesolve_kwargs, termination_condition, alias, kwargs... @@ -239,7 +277,10 @@ function SciMLBase.__solve( setup.run_kwargs..., odesolve_kwargs..., setup.callback, save_end = true, alias = setup.odealias ) + return __steady_state_solution(prob, alg, setup, odesol, save_idxs) +end +function __steady_state_solution(prob, alg::DynamicSS, setup, odesol, save_idxs) resid, u, retcode = __get_result_from_sol(setup.tc_cache, odesol) if save_idxs !== nothing @@ -253,29 +294,29 @@ function SciMLBase.__solve( ) end -# `init` returns the ODE integrator with the steady-state termination callback -# installed. An SCC lowering solves its blocks sequentially with no single -# integrator, so it is solved eagerly and the finished solution is returned. function SciMLBase.__init( prob::SciMLBase.AbstractSteadyStateProblem, alg::DynamicSS, args...; abstol = 1.0e-8, reltol = 1.0e-6, odesolve_kwargs = (;), save_idxs = nothing, termination_condition = NonlinearSolveBase.NormTerminationMode(infnorm), alias = SciMLBase.NonlinearAliasSpecifier(), kwargs... ) - sccsol = __solve_scc_lowering( - prob, alg, args...; abstol, reltol, odesolve_kwargs, - termination_condition, alias, save_idxs, kwargs... + lp = __scc_lowering(prob) + lp !== nothing && return SteadyStateSCCCache( + prob, lp, alg, args, (; + abstol, reltol, odesolve_kwargs, termination_condition, alias, + save_idxs, kwargs..., + ) ) - sccsol !== nothing && return sccsol setup = __dynamicss_ode_setup( prob, alg; abstol, reltol, odesolve_kwargs, termination_condition, alias, kwargs... ) - return init( + integrator = init( setup.odeprob, alg.alg, args...; setup.abstol, setup.reltol, setup.run_kwargs..., odesolve_kwargs..., setup.callback, save_end = true, - alias = setup.odealias, save_idxs + alias = setup.odealias ) + return SteadyStateODECache(prob, alg, integrator, setup, save_idxs) end # SICNM: Semi-Implicit Continuous Newton Method @@ -430,11 +471,11 @@ function SciMLBase.__solve( termination_condition = NonlinearSolveBase.AbsNormTerminationMode(infnorm), alias = SciMLBase.NonlinearAliasSpecifier(), kwargs... ) - sccsol = __solve_scc_lowering( - prob, alg, args...; abstol, reltol, odesolve_kwargs, + lp = __scc_lowering(prob) + lp !== nothing && return __solve_scc_lowering( + prob, lp, alg, args...; abstol, reltol, odesolve_kwargs, termination_condition, alias, save_idxs, kwargs... ) - sccsol !== nothing && return sccsol setup = __sicnm_ode_setup( prob, alg; abstol, reltol, odesolve_kwargs, termination_condition, kwargs... @@ -444,7 +485,10 @@ function SciMLBase.__solve( reltol = setup.ode_reltol, setup.run_kwargs..., odesolve_kwargs..., setup.callback, save_end = true ) + return __steady_state_solution(prob, alg, setup, odesol, save_idxs) +end +function __steady_state_solution(prob, alg::SICNM, setup, odesol, save_idxs) u, retcode = __sicnm_result(setup.tc_cache, odesol, setup.n) resid = if setup.iip setup.g(setup.gbuf, u) @@ -464,7 +508,6 @@ function SciMLBase.__solve( ) end -# As for `DynamicSS`, an SCC lowering is solved eagerly. function SciMLBase.__init( prob::SciMLBase.AbstractSteadyStateProblem, alg::SICNM, args...; abstol = 1.0e-8, reltol = 1.0e-6, odesolve_kwargs = (;), @@ -472,20 +515,23 @@ function SciMLBase.__init( termination_condition = NonlinearSolveBase.AbsNormTerminationMode(infnorm), alias = SciMLBase.NonlinearAliasSpecifier(), kwargs... ) - sccsol = __solve_scc_lowering( - prob, alg, args...; abstol, reltol, odesolve_kwargs, - termination_condition, alias, save_idxs, kwargs... + lp = __scc_lowering(prob) + lp !== nothing && return SteadyStateSCCCache( + prob, lp, alg, args, (; + abstol, reltol, odesolve_kwargs, termination_condition, alias, + save_idxs, kwargs..., + ) ) - sccsol !== nothing && return sccsol setup = __sicnm_ode_setup( prob, alg; abstol, reltol, odesolve_kwargs, termination_condition, kwargs... ) - return init( + integrator = init( setup.odeprob, alg.alg, args...; abstol = setup.ode_abstol, reltol = setup.ode_reltol, setup.run_kwargs..., odesolve_kwargs..., - setup.callback, save_end = true, save_idxs + setup.callback, save_end = true ) + return SteadyStateODECache(prob, alg, integrator, setup, save_idxs) end function __sicnm_result(tc_cache, odesol, n) diff --git a/test/scc.jl b/test/scc.jl index 9639045..31ec16b 100644 --- a/test/scc.jl +++ b/test/scc.jl @@ -2,8 +2,9 @@ using SteadyStateDiffEq, NonlinearSolve, OrdinaryDiffEq, Test using ModelingToolkit using ModelingToolkit: t_nounits as t, D_nounits as D using SCCNonlinearSolve: SCCAlg +using NonlinearSolve.NonlinearSolveBase: AbsNormSafeTerminationMode using SciMLBase: HomotopyProblem, LinearProblem, NonlinearProblem, SCCNonlinearProblem, - SteadyStateSolution + SteadyStateSolution, step! function coupled_scc_problem(iip, use_vector) f = if iip @@ -436,38 +437,84 @@ end @test sol.original isa Tuple{SciMLBase.LinearSolution, NonlinearSolution} end -# `init` with `DynamicSS`/`SICNM` uses these algorithms rather than NonlinearSolve's -# default: a plain problem gives a live ODE integrator, an SCC lowering is solved -# eagerly and returns the finished solution. -@testset "init on DynamicSS/SICNM does not fall through to NonlinearSolve's default" begin - @testset "plain SteadyStateProblem returns a live integrator" for alg in ( +# `solve!(init(prob, alg))` goes through the same steady-state finalization as +# `solve(prob, alg)`: termination-condition retcodes, best-state selection, +# `save_idxs`, and the sequential SCC solve of a stored lowering. +@testset "init/solve! on DynamicSS/SICNM matches solve" begin + function test_matches_solve(prob, alg; kwargs...) + sol = solve!(init(prob, alg; kwargs...)) + ref = solve(prob, alg; kwargs...) + @test sol isa NonlinearSolution + @test sol.retcode == ref.retcode + @test sol.u ≈ ref.u + @test sol.resid ≈ ref.resid + return sol + end + + @testset "plain SteadyStateProblem, alg=$(nameof(typeof(alg)))" for alg in ( DynamicSS(Tsit5()), SICNM(Rodas5P()), ) - prob = SteadyStateProblem((u, p, t) -> 1 .- u, [0.0]) - integ = init(prob, alg; abstol = 1.0e-10, reltol = 1.0e-10) - sol = solve!(integ) + prob = SteadyStateProblem((u, p, t) -> [1, 2] .- u, [0.0, 0.0]) + sol = test_matches_solve(prob, alg; abstol = 1.0e-10, reltol = 1.0e-10) @test successful_retcode(sol) - # `SICNM`'s integrator carries the extended DAE state `[y; z]`, so only - # the first `length(prob.u0)` components are the original residual state. - @test sol.u[end][1:1] ≈ [1.0] atol = 1.0e-6 - - prob2 = SteadyStateProblem((u, p, t) -> [1, 2] .- u, [0.0, 0.0]) - sol2 = solve!(init(prob2, alg; save_idxs = [2], abstol = 1.0e-10, reltol = 1.0e-10)) - @test length(sol2.u[end]) == 1 - @test sol2.u[end] ≈ [2.0] atol = 1.0e-6 + @test sol.u ≈ [1.0, 2.0] atol = 1.0e-6 + + sol = test_matches_solve( + prob, alg; save_idxs = [2], abstol = 1.0e-10, reltol = 1.0e-10 + ) + @test sol.u ≈ [2.0] atol = 1.0e-6 + + cache = init(prob, alg; abstol = 1.0e-10, reltol = 1.0e-10) + step!(cache) + @test cache.integrator.t > 0 + @test successful_retcode(solve!(cache)) end - @testset "manually-built SCC lowering, alg=$alg" for alg in ( - DynamicSS(Tsit5()), SICNM(Rodas5P()), + @testset "finite-time nonconvergence, alg=$(nameof(typeof(alg)))" for alg in ( + DynamicSS(Tsit5(); tspan = 1.0e-3), SICNM(Rodas5P(); tspan = 1.0e-3), + ) + prob = SteadyStateProblem((u, p, t) -> [1, 2] .- u, [0.0, 0.0]) + sol = test_matches_solve(prob, alg; abstol = 1.0e-10, reltol = 1.0e-10) + @test !successful_retcode(sol) + end + + @testset "protective termination" begin + # `u' = u` grows past the protective threshold. + prob = SteadyStateProblem((u, p, t) -> u, [1.0]) + tc = AbsNormSafeTerminationMode(u -> maximum(abs, u); protective_threshold = 1.01) + sol = test_matches_solve( + prob, DynamicSS(Tsit5()); abstol = 1.0e-10, reltol = 1.0e-10, + termination_condition = tc + ) + @test sol.retcode == ReturnCode.Unstable + end + + @testset "manually-built SCC lowering, alg=$(nameof(typeof(alg)))" for (alg, shortalg) in ( + (DynamicSS(Tsit5()), DynamicSS(Tsit5(); tspan = 1.0e-3)), + (SICNM(Rodas5P()), SICNM(Rodas5P(); tspan = 1.0e-3)), ) sccprob = dynamicss_scc_problem(false, false) prob = SteadyStateProblem( (u, p, t) -> 1 .- u, [0.0, 0.0]; lowered_problem = sccprob ) - sol = init(prob, alg; abstol = 1.0e-10, reltol = 1.0e-10) + sol = test_matches_solve(prob, alg; abstol = 1.0e-10, reltol = 1.0e-10) @test successful_retcode(sol) @test sol.u ≈ [1, 2, 1, 2] atol = 1.0e-8 - @test sol.prob === sccprob + + sol = test_matches_solve( + prob, alg; save_idxs = [2, 4], abstol = 1.0e-10, reltol = 1.0e-10 + ) + @test sol.u ≈ [2, 2] atol = 1.0e-8 + + # The block cannot reach steady state within the short time span. + failing = SCCNonlinearProblem( + (NonlinearProblem((u, p) -> 1 .- u, [0.0]),), (Returns(nothing),) + ) + prob = SteadyStateProblem( + (u, p, t) -> 1 .- u, [0.0]; lowered_problem = failing + ) + sol = test_matches_solve(prob, shortalg; abstol = 1.0e-10, reltol = 1.0e-10) + @test !successful_retcode(sol) end @testset "ModelingToolkit SCC decomposition, alg=$(nameof(typeof(alg)))" for alg in ( @@ -480,9 +527,12 @@ end sys = mtkcompile(model) prob = SteadyStateProblem(sys, [a => 0.8, b => 1.8, x => 0.8]) - sol = init(prob, alg; abstol = 1.0e-10, reltol = 1.0e-10) + sol = test_matches_solve(prob, alg; abstol = 1.0e-10, reltol = 1.0e-10) @test successful_retcode(sol) @test sol[[a, b, x]] ≈ [1, 2, cbrt(3)] atol = 1.0e-8 - @test sol.prob isa SCCNonlinearProblem + + test_matches_solve( + prob, alg; save_idxs = [1], abstol = 1.0e-10, reltol = 1.0e-10 + ) end end