From a3b9bd170e55bfbd3050566e421da1f0bcc4363d Mon Sep 17 00:00:00 2001 From: Chris Rackauckas Date: Wed, 14 Oct 2026 19:15:25 +0000 Subject: [PATCH 01/11] Solve SCCSteadyStateProblem; stop routing SteadyStateProblem through a lowering DynamicSS and SICNM solve SCCSteadyStateProblem block by block, SSRootfind and nonlinear algorithms through SCCNonlinearProblem(prob); the solution is built on the SCC problem in its compressed layout, and save_idxs selects in that layout (integers or states via variable_index). SteadyStateProblem solves keep the problem's own layout again, fixing https://github.com/SciML/SteadyStateDiffEq.jl/issues/171. Ports the LV, chain, mixed-SCC, resid, save_idxs and eltype cases of #175/#177. Co-Authored-By: Chris Rackauckas Co-Authored-By: Devin Agent-Harness: Devin CLI 3000.11.3 Agent-Model: fusion-claude-opus-5-5-high-sidekick-swe-2-medium Agent-Session: local headless run, log /Users/chrisrackauckas/sandbox/master/green/scc-ssprob/log.txt --- Project.toml | 6 +- docs/src/index.md | 17 ++ src/SteadyStateDiffEq.jl | 10 +- src/algorithms.jl | 23 ++- src/solve.jl | 165 +++++++++--------- test/qa/qa.jl | 4 +- test/scc.jl | 359 ++++++++++++++++++++++++++------------- 7 files changed, 364 insertions(+), 220 deletions(-) diff --git a/Project.toml b/Project.toml index 9094e20..584fb95 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "SteadyStateDiffEq" uuid = "9672c7b4-1e72-59bd-8a11-6ac3964bc41f" -version = "2.18.0" +version = "2.19.0" [deps] ConcreteStructs = "2569d6c7-a4a2-43d3-a901-331e8e4be471" @@ -24,14 +24,14 @@ DiffEqCallbacks = "4.16.0" ForwardDiff = "1.0.1" LinearAlgebra = "1.10" LinearSolve = "5.4" -ModelingToolkit = "11.44" +ModelingToolkit = "11.46" NonlinearSolve = "4.12" NonlinearSolveBase = "2.40" OrdinaryDiffEq = "7" PrecompileTools = "1.2.1" SafeTestsets = "0.1" SCCNonlinearSolve = "1.13" -SciMLBase = "3.55" +SciMLBase = "3.59" SciMLPublic = "1" SciMLTesting = "2.4" Sundials = "6.4.1" diff --git a/docs/src/index.md b/docs/src/index.md index bae6ad1..2e38d25 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -48,6 +48,23 @@ prob = SteadyStateProblem((u, p, t) -> 1 .- u, [0.0]) sol = solve(prob, SICNM(Rodas3d())) ``` +## SCC-split steady states + +`SCCSteadyStateProblem` is the SCC-split form of a steady-state problem: its +residual decomposes into a block-lower-triangular sequence of `LinearProblem`/ +`NonlinearProblem`/`HomotopyProblem` blocks. Build it from a block list with +`SCCSteadyStateProblem(probs, explicitfuns!, p)`, or from a model with +`ModelingToolkit.SCCSteadyStateProblem(sys, op)`. + +`DynamicSS` and `SICNM` solve the blocks sequentially in SCC order — each +nonlinear block is integrated to steady state on its own residual after +`explicitfuns!` updates its caches from the upstream solutions. `SSRootfind`, a +bare `solve`, or any nonlinear algorithm solve the problem through the SCC +solver (`SCCNonlinearSolve.jl`). The solution's `u` and `resid` are in the +compressed SCC state layout; symbolic indexing (`sol[x]`, `sol.ps[p]`) goes +through `prob.f.sys`, and `save_idxs` is a selection in that layout — integers, +or states resolved through `variable_index`. + ## API ```@docs diff --git a/src/SteadyStateDiffEq.jl b/src/SteadyStateDiffEq.jl index e508efc..c7e5818 100644 --- a/src/SteadyStateDiffEq.jl +++ b/src/SteadyStateDiffEq.jl @@ -9,9 +9,11 @@ using LinearAlgebra: Diagonal, norm using LinearSolve: LinearSolve using SciMLPublic: @public using SciMLBase: SciMLBase, CallbackSet, LinearProblem, NonlinearProblem, ODEProblem, - NonlinearSolution, ReturnCode, SteadyStateProblem, SteadyStateSolution, get_du, init, + NonlinearSolution, ReturnCode, SCCNonlinearProblem, SCCSteadyStateProblem, + SteadyStateProblem, SteadyStateSolution, get_du, init, isinplace, remake, solve, successful_retcode -using SymbolicIndexingInterface: parameter_values +using SymbolicIndexingInterface: + NotSymbolic, parameter_values, state_values, symbolic_type, variable_index const infnorm = Base.Fix2(norm, Inf) @@ -20,8 +22,8 @@ include("solve.jl") include("precompilation.jl") export SSRootfind, DynamicSS, SICNM -export NonlinearProblem, NonlinearSolution, ReturnCode, SteadyStateProblem, - SteadyStateSolution, remake, solve, successful_retcode +export NonlinearProblem, NonlinearSolution, ReturnCode, SCCSteadyStateProblem, + SteadyStateProblem, SteadyStateSolution, remake, solve, successful_retcode @public SteadyStateDiffEqAlgorithm end diff --git a/src/algorithms.jl b/src/algorithms.jl index 9e6e415..c28f2d6 100644 --- a/src/algorithms.jl +++ b/src/algorithms.jl @@ -92,10 +92,10 @@ requires the SCC solver implementation from `SCCNonlinearSolve.jl` to be loaded automatically), or an `SCCAlg` configured with separate linear and nonlinear block solvers. -A `SteadyStateProblem` that records a `lowered_problem` (for example the SCC -decomposition stored by `ModelingToolkit`) is solved through that lowering; -since the lowering defines its own state ordering, the returned solution is -expressed on the lowered problem. +An `SCCSteadyStateProblem` (for example the SCC decomposition built by +`ModelingToolkit.SCCSteadyStateProblem`) is solved through +`SCCNonlinearProblem(prob)`: it is forwarded to `alg` unchanged, and the +solution is expressed on the SCC problem's compressed state ordering. # Arguments @@ -148,10 +148,10 @@ solutions through `explicitfuns!` as in the SCC solve. `LinearProblem` blocks are solved directly, and the remaining blocks are integrated to steady state by `DynamicSS` on the block residual — so convergence only requires each nonlinear block's residual to be attracting under its own pseudo-transient -dynamics, not a globally attracting concatenated field. A `SteadyStateProblem` -that records an `SCCNonlinearProblem` as its `lowered_problem` (for example one -built by `ModelingToolkit`) takes the same sequential solve, and its solution -is expressed on the lowered problem's state ordering. +dynamics, not a globally attracting concatenated field. An +`SCCSteadyStateProblem` (for example one built by +`ModelingToolkit.SCCSteadyStateProblem`) takes the same sequential solve, and +its solution is expressed on the compressed state ordering. # Arguments @@ -251,10 +251,9 @@ solutions through `explicitfuns!` as in the SCC solve. `LinearProblem` blocks are solved directly, and the remaining blocks each run their own SICNM continuous-Newton flow on the block residual — so convergence only requires each nonlinear block to be tractable for SICNM on its own, not the concatenated -system. A `SteadyStateProblem` that records an `SCCNonlinearProblem` as its -`lowered_problem` (for example one built by `ModelingToolkit`) takes the same -sequential solve, and its solution is expressed on the lowered problem's state -ordering. +system. An `SCCSteadyStateProblem` (for example one built by +`ModelingToolkit.SCCSteadyStateProblem`) takes the same sequential solve, and +its solution is expressed on the compressed state ordering. # Arguments diff --git a/src/solve.jl b/src/solve.jl index acacbee..a50d485 100644 --- a/src/solve.jl +++ b/src/solve.jl @@ -10,25 +10,11 @@ function SciMLBase.__solve( args...; kwargs... ) nlprob = NonlinearProblem(prob) - fwd = kwargs - if nlprob isa SciMLBase.SCCNonlinearProblem - # The nonlinear solve pipeline injects `alias`/`verbose` as nonlinear - # specifier types, which `LinearProblem` blocks in the SCC solve do not - # understand. - fwd = (; (n => v for (n, v) in pairs(kwargs) if n !== :alias && n !== :verbose)...) - end - nlsol = solve(nlprob, alg.alg, args...; fwd...) - # A stored `lowered_problem` (e.g. an `SCCNonlinearProblem`) solves in the - # lowering's own state ordering, so its solution is expressed on the - # lowering rather than the steady-state problem. - solprob = if prob isa SteadyStateProblem && prob.lowered_problem !== nothing - nlprob - else - prob - end - return __build_ssrootfind_solution(solprob, nlsol) + nlsol = solve(nlprob, alg.alg, args...; kwargs...) + return __build_ssrootfind_solution(prob, nlsol) end + # An SCCNonlinearProblem has no top-level `u0`/`kwargs` fields, so it cannot go # through the generic AbstractNonlinearProblem solve preprocessing. Forward it # directly to the wrapped algorithm instead, which dispatches to the SCC solver @@ -41,16 +27,16 @@ function SciMLBase.solve( return __build_ssrootfind_solution(prob, nlsol) end -# `DynamicSS`/`SICNM` on an `SCCNonlinearProblem` solve the blocks sequentially +# `DynamicSS`/`SICNM` on an SCC problem solve the blocks sequentially # in SCC order: each block's parameter cache is updated from the upstream # solutions through `explicitfuns!`, `LinearProblem` blocks are solved directly, # and the remaining blocks are integrated to steady state on the block # residual. `explicitfuns!` mutate the block caches, matching the # `SCCNonlinearSolve` protocol; each upstream solution is stripped to a plain # solution so generated explicit functions can index `sols[j][k]`/`sols[j].u`. -function SciMLBase.solve( - prob::SciMLBase.SCCNonlinearProblem, alg::Union{DynamicSS, SICNM}, - args...; kwargs... +function __solve_scc_blocks( + prob::Union{SciMLBase.SCCNonlinearProblem, SciMLBase.SCCSteadyStateProblem}, + alg::Union{DynamicSS, SICNM}, args...; kwargs... ) sols = () for i in eachindex(prob.probs) @@ -58,7 +44,14 @@ function SciMLBase.solve( SciMLBase.invoke_with_despecialized_parameters( prob.explicitfuns![i], (parameter_values(block), sols) ) - _sol = if block isa LinearProblem + _sol = if state_values(block) !== nothing && isempty(state_values(block)) + # A stateless block's residual is already satisfied; there is + # nothing to integrate. + SciMLBase.build_solution( + block, nothing, copy(state_values(block)), + copy(state_values(block)); retcode = ReturnCode.Success + ) + elseif block isa LinearProblem A, b = block.A, block.b # `remake` recomputes `A` and `b` through a `SymbolicLinearInterface` # from the cache that `explicitfun` just updated; for a plain @@ -86,8 +79,8 @@ function SciMLBase.solve( sols = (sols..., _sol) end - u = vcat(map(sol -> sol.u, sols)...) - resid = vcat(map(sol -> sol.resid, sols)...) + u = isempty(sols) ? Float64[] : vcat(map(sol -> sol.u, sols)...) + resid = isempty(sols) ? Float64[] : vcat(map(sol -> sol.resid, sols)...) retcode = SciMLBase.ReturnCode.Success for s in sols if !SciMLBase.successful_retcode(s) @@ -98,54 +91,84 @@ 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... - ) - 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 - sccsol = solve(lp, alg, args...; kwargs...) - save_idxs === nothing && return sccsol - return SciMLBase.build_solution( - lp, sccsol.alg, sccsol.u[save_idxs], sccsol.resid[save_idxs]; - retcode = sccsol.retcode, original = sccsol +function SciMLBase.solve( + prob::SciMLBase.SCCNonlinearProblem, alg::Union{DynamicSS, SICNM}, + args...; kwargs... ) + return __solve_scc_blocks(prob, alg, args...; kwargs...) end -# A `SteadyStateProblem`'s stored SCC lowering materializes through -# `SciMLBase.NonlinearProblem(prob)` during `solve`. When it is an -# `SCCNonlinearProblem`, the generic nonlinear-solve path re-dispatches -# `__solve` on it and the fallbacks try to convert it back via `prob.u0` — a -# field it does not have — so the lowering goes through its own `solve` -# dispatch instead, which routes bare and `AbstractNonlinearAlgorithm` solves -# to the SCC solver (`SCCNonlinearSolve.jl`, loaded downstream). Algorithms -# that do not lower to a nonlinear solve (`DynamicSS`, ODE algorithms) keep -# their own `__solve` dispatch; non-SCC lowerings keep the upstream path -# verbatim via `invoke`. -function SciMLBase.solve(prob::SteadyStateProblem, args...; kwargs...) +function SciMLBase.solve( + prob::SciMLBase.SCCSteadyStateProblem, args...; + save_idxs = nothing, kwargs... + ) alg = isempty(args) ? nothing : first(args) alg === nothing && (alg = get(kwargs, :alg, nothing)) - alg === nothing && prob.kwargs !== nothing && - (alg = get(prob.kwargs, :alg, nothing)) - if alg !== nothing && !(alg isa SciMLBase.AbstractNonlinearAlgorithm) - return invoke( - solve, Tuple{SciMLBase.AbstractNonlinearProblem, Vararg{Any}}, - prob, args...; kwargs... + kwargs = (; (n => v for (n, v) in pairs(kwargs) if n !== :alg)...) + rest = isempty(args) ? () : Base.tail(args) + sol = __solve_scc_steady_state(prob, alg, rest...; kwargs...) + return __scc_save_idxs(prob, sol, save_idxs) +end + +function __solve_scc_steady_state( + prob::SciMLBase.SCCSteadyStateProblem, alg::Union{DynamicSS, SICNM}, + args...; kwargs... + ) + return __solve_scc_blocks(prob, alg, args...; kwargs...) +end + +function __solve_scc_steady_state( + prob::SciMLBase.SCCSteadyStateProblem, alg::SSRootfind, + args...; kwargs... + ) + nlsol = solve(SCCNonlinearProblem(prob), alg.alg, args...; kwargs...) + return SciMLBase.build_solution( + prob, SSRootfind(nlsol.alg), nlsol.u, nlsol.resid; + nlsol.retcode, nlsol.stats, nlsol.left, nlsol.right, original = nlsol + ) +end + +function __solve_scc_steady_state( + prob::SciMLBase.SCCSteadyStateProblem, ::Nothing, + args...; kwargs... + ) + nlsol = solve(SCCNonlinearProblem(prob), args...; kwargs...) + return SciMLBase.build_solution( + prob, nlsol.alg, nlsol.u, nlsol.resid; + nlsol.retcode, nlsol.stats, nlsol.left, nlsol.right, original = nlsol + ) +end + +function __solve_scc_steady_state( + prob::SciMLBase.SCCSteadyStateProblem, alg, + args...; kwargs... + ) + nlsol = solve(SCCNonlinearProblem(prob), alg, args...; kwargs...) + return SciMLBase.build_solution( + prob, nlsol.alg, nlsol.u, nlsol.resid; + nlsol.retcode, nlsol.stats, nlsol.left, nlsol.right, original = nlsol + ) +end + +# `save_idxs` is a selection in the compressed state layout: integers pass +# through, symbolic entries are resolved through `prob`'s index provider and +# must be states of it. +function __scc_save_idxs( + prob::SciMLBase.SCCSteadyStateProblem, sol, save_idxs + ) + save_idxs === nothing && return sol + idxs = map(save_idxs) do s + symbolic_type(s) === NotSymbolic() && return s + idx = variable_index(prob, s) + idx === nothing && throw( + ArgumentError("`save_idxs` entry $s is not a state of the problem") ) + return idx end - # Concretize like `NonlinearSolveBase.solve_up` (evaluates `u0(p, Inf)`-style - # `u0`s, promotes integer `u0`s, applies `u0`/`p` solve kwargs) so a - # callable `lowered_problem` sees the updated operating point. - _prob = NonlinearSolveBase.get_concrete_problem(prob; kwargs...) - nlprob = NonlinearProblem(_prob) - nlprob isa SciMLBase.SCCNonlinearProblem && return solve(nlprob, args...; kwargs...) - return invoke( - solve, Tuple{SciMLBase.AbstractNonlinearProblem, Vararg{Any}}, - _prob, args...; kwargs... + return SciMLBase.build_solution( + prob, sol.alg, sol.u[idxs], sol.resid[idxs]; + retcode = sol.retcode, stats = sol.stats, left = sol.left, + right = sol.right, original = sol ) end @@ -167,12 +190,6 @@ 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, - termination_condition, alias, save_idxs, kwargs... - ) - sccsol !== nothing && return sccsol - tspan = __get_tspan(prob.u0, alg) f = if prob isa SteadyStateProblem @@ -277,12 +294,6 @@ function SciMLBase.__solve( 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 - prob.u0 isa AbstractVector || throw(ArgumentError("SICNM currently only supports `AbstractVector` initial conditions")) tspan = __get_tspan(prob.u0, alg) diff --git a/test/qa/qa.jl b/test/qa/qa.jl index f5224f9..da958a2 100644 --- a/test/qa/qa.jl +++ b/test/qa/qa.jl @@ -2,8 +2,8 @@ using SciMLTesting, SteadyStateDiffEq, Test using JET const REEXPORTS = ( - :NonlinearProblem, :NonlinearSolution, :ReturnCode, :SteadyStateProblem, - :SteadyStateSolution, :remake, :solve, :successful_retcode, + :NonlinearProblem, :NonlinearSolution, :ReturnCode, :SCCSteadyStateProblem, + :SteadyStateProblem, :SteadyStateSolution, :remake, :solve, :successful_retcode, ) const QUALIFIED_PUBLIC_IGNORE = ( diff --git a/test/scc.jl b/test/scc.jl index 4780d4f..c93a8b0 100644 --- a/test/scc.jl +++ b/test/scc.jl @@ -2,6 +2,7 @@ using SteadyStateDiffEq, NonlinearSolve, OrdinaryDiffEq, Test using ModelingToolkit using ModelingToolkit: t_nounits as t, D_nounits as D using SCCNonlinearSolve: SCCAlg +using SymbolicIndexingInterface: state_values, variable_index, variable_symbols using SciMLBase: HomotopyProblem, LinearProblem, NonlinearProblem, SCCNonlinearProblem, SteadyStateSolution @@ -206,78 +207,260 @@ end end end -# A `SteadyStateProblem` recording an `SCCNonlinearProblem` lowering is solved -# through it: `DynamicSS` runs the sequential block solve and `SSRootfind` -# forwards the lowering to the nonlinear solver. The solution is expressed on -# the lowering (whose ordering need not match `prob.u0`), so `sol.prob` is the -# `SCCNonlinearProblem`. -@testset "SteadyStateProblem with an SCC lowering" begin - f_iip(du, u, p, t) = (du .= 1 .- u) - f_oop(u, p, t) = 1 .- u - @testset "iip=$iip builder=$builder alg=$alg" for iip in (false, true), - builder in (false, true), + +# An `SCCSteadyStateProblem` built from the same blocks solves through the same +# sequential/block or SCC-solver paths, with the solution expressed on the SCC +# problem itself. +@testset "SCCSteadyStateProblem from blocks" begin + @testset "iip=$iip vector=$use_vector alg=$alg" for iip in (false, true), + use_vector in (false, true), alg in ( DynamicSS(Tsit5()), SICNM(Rodas5P()), - SSRootfind(NewtonRaphson()), + SSRootfind(SCCAlg()), SSRootfind(SCCAlg(; nlalg = NewtonRaphson())), + NewtonRaphson(), + nothing, ) - sccprob = dynamicss_scc_problem(iip, false) - lowered = builder ? (prob -> sccprob) : sccprob - prob = SteadyStateProblem( - iip ? f_iip : f_oop, [0.0, 0.0]; lowered_problem = lowered + prob = SciMLBase.SCCSteadyStateProblem(dynamicss_scc_problem(iip, use_vector)) + sol = alg === nothing ? solve(prob) : + solve(prob, alg; abstol = 1.0e-10, reltol = 1.0e-10) + @test successful_retcode(sol) + @test sol.prob === prob + @test sol.u ≈ [1, 2, 1, 2] atol = 1.0e-8 + @test sol.u isa Vector{Float64} + @test sol.resid isa Vector{Float64} + @test length(sol.resid) == length(sol.u) + @test sol.resid ≈ zero(sol.resid) atol = 1.0e-7 + if alg isa Union{DynamicSS, SICNM} + @test sol.original isa Tuple + @test sol.original[1] isa SciMLBase.LinearSolution + else + @test sol.original isa NonlinearSolution + @test sol.original.prob isa SCCNonlinearProblem + end + end + + # `save_idxs` selects entries of the compressed SCC state/resid layout and + # keeps the full solution in `original`. + @testset "save_idxs" for alg in (DynamicSS(Tsit5()), NewtonRaphson()) + prob = SciMLBase.SCCSteadyStateProblem(dynamicss_scc_problem(false, false)) + full = solve(prob, alg; abstol = 1.0e-10, reltol = 1.0e-10) + saved = solve( + prob, alg; save_idxs = [2, 4], abstol = 1.0e-10, reltol = 1.0e-10 ) + @test saved.u ≈ full.u[[2, 4]] + @test saved.resid ≈ full.resid[[2, 4]] + @test saved.prob === prob + @test saved.original.u ≈ full.u + end + + # A `tspan` too short to settle gives a failure retcode, like DynamicSS on + # a `SteadyStateProblem`. + @testset "short tspan" for alg in (DynamicSS(Tsit5(); tspan = 1.0e-3), SICNM(Rodas5P(); tspan = 1.0e-3)) + prob = SciMLBase.SCCSteadyStateProblem(dynamicss_scc_problem(false, false)) + sol = solve(prob, alg) + @test !successful_retcode(sol) + @test sol.prob === prob + end +end + +@testset "SSRootfind on a ModelingToolkit SCCSteadyStateProblem" begin + @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 = SciMLBase.SCCSteadyStateProblem(sys, [a => 0.8, b => 1.8, x => 0.8]) + + sol = solve(prob, SSRootfind(NewtonRaphson()); abstol = 1.0e-12, reltol = 1.0e-12) + @test successful_retcode(sol) + @test sol[[a, b, x]] ≈ [1, 2, cbrt(3)] atol = 1.0e-8 + @test sol.prob === prob +end + +@testset "DynamicSS on a ModelingToolkit SCCSteadyStateProblem" begin + # Scalar nonlinear block `x' = 3 - x^3` has a unique real root `∛3`, so the + # per-block pseudo-transient solve and the monolithic integration converge + # to the same steady state. + @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 = SciMLBase.SCCSteadyStateProblem(sys, [a => 0.8, b => 1.8, x => 0.8]) + + sol = solve(prob, DynamicSS(); 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 === prob + @test sol.original isa Tuple{SciMLBase.LinearSolution, NonlinearSolution} +end + +@testset "SICNM on a ModelingToolkit SCCSteadyStateProblem" begin + # Scalar nonlinear block `3 - x^3` has a unique real root `∛3` that is + # attracting under SICNM's continuous-Newton flow. + @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 = SciMLBase.SCCSteadyStateProblem(sys, [a => 0.8, b => 1.8, x => 0.8]) + + sol = solve(prob, SICNM(Rodas5P()); 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 === prob + @test sol.original isa Tuple{SciMLBase.LinearSolution, NonlinearSolution} +end + +# Systems shared by the two solution-layout testsets below: `op` is the +# operating point, `states`/`expected` the steady-state values, `observed` an +# `sym => value` pair of the `osys` system's observed `XY`. +function layout_cases() + @parameters ca = 1.0 cb = 2.0 + @variables x(t) y(t) z(t) + chain = mtkcompile( + System([D(x) ~ ca - x, D(y) ~ x - cb * y, D(z) ~ y - z], t; name = :chain) + ) + + @parameters α = 1.5 β = 1.0 γ = 3.0 δ = 1.0 + @variables u(t) v(t) + lv = mtkcompile( + System([D(u) ~ α * u - β * u * v, D(v) ~ -γ * v + δ * u * v], t; name = :lv) + ) + + @variables a(t) b(t) mx(t) [irreducible = true] + mix = mtkcompile( + System( + [D(a) ~ 5 - 3a - b, D(b) ~ 5 - a - 2b, D(mx) ~ a + b - mx^3], t; + name = :mix + ) + ) + + @parameters kp kd k1 k2 + @variables X(t) Y(t) XY(t) + osys = complete( + System( + [D(X) ~ kp - kd * X - k1 * X + k2 * Y, D(Y) ~ 1 + k1 * X - k2 * Y - Y], t; + observed = [XY ~ X + Y], name = :osys + ) + ) + Xss, Yss = [0.35 -0.5; -0.25 1.5] \ [1.0, 1.0] + + # `DynamicSS(Tsit5())` over `tspan = Inf` does not reach its terminator on + # this system — the integrator exhausts maxiters while the residual is still + # well above the tolerance — so the pseudo-transient runs use a finite span. + allalgs = ( + NewtonRaphson(), SSRootfind(NewtonRaphson()), DynamicSS(Tsit5(); tspan = 100), + SICNM(Rodas5P()), + ) + nlalgs = (NewtonRaphson(), SSRootfind(NewtonRaphson())) + return ( + (; + sys = chain, op = [x => 0.1, y => 0.2, z => 0.3], states = [x, y, z], + expected = [1.0, 0.5, 0.5], algs = allalgs, observed = nothing, + ), + # Lotka–Volterra's nonzero steady state is a center, so the + # pseudo-transient methods do not apply. + (; + sys = lv, op = [u => 2.9, v => 1.6], states = [u, v], + expected = [3.0, 1.5], algs = nlalgs, observed = nothing, + ), + (; + sys = mix, op = [a => 0.8, b => 1.8, mx => 0.8], states = [a, b, mx], + expected = [1.0, 2.0, cbrt(3)], algs = allalgs, observed = nothing, + ), + (; + sys = osys, + op = [X => 4.0, Y => 5.0, kp => 1.0, kd => 0.1, k1 => 0.25, k2 => 0.5], + states = [X, Y], expected = [Xss, Yss], algs = allalgs, + observed = XY => Xss + Yss, + ), + ) +end + +# https://github.com/SciML/SteadyStateDiffEq.jl/issues/171 +# A `SteadyStateProblem` solves in its own layout: `sol.u`/`sol.resid` have the +# length and ordering of `prob.u0`. Nonlinear algorithms return their solution +# on the converted `NonlinearProblem`; the steady-state algorithms on `prob`. +@testset "SteadyStateProblem solutions are in the problem's layout" begin + @testset "case=$(nameof(case.sys)) alg=$(nameof(typeof(alg)))" for + case in layout_cases(), alg in case.algs + + (; sys, op, states, expected, observed) = case + prob = SteadyStateProblem(sys, op) sol = 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 - @test sol.original !== nothing + @test sol.prob isa ( + alg isa SciMLBase.AbstractNonlinearAlgorithm ? + NonlinearProblem : SteadyStateProblem + ) + @test length(sol.u) == length(prob.u0) + @test sol.u[variable_index.((sys,), states)] ≈ expected atol = 1.0e-7 + @test sol[states] ≈ expected atol = 1.0e-7 + observed === nothing || @test sol[first(observed)] ≈ last(observed) atol = 1.0e-7 + @test length(sol.resid) == length(prob.u0) + @test sol.resid ≈ zero(sol.resid) atol = 1.0e-7 + + if alg isa Union{DynamicSS, SICNM} + saved = solve( + prob, alg; save_idxs = [2], abstol = 1.0e-10, reltol = 1.0e-10 + ) + @test saved.u ≈ sol.u[2:2] atol = 1.0e-7 + @test length(saved.resid) == 1 + end end - # Bare and `AbstractNonlinearAlgorithm` solves have no `__solve` dispatch of - # their own here; they land on the generic nonlinear-solve path, which - # converts back through `prob.u0` — absent on `SCCNonlinearProblem`. `solve` - # forwards the lowering to its own `solve` dispatch (the SCC solver) - # instead. - @testset "bare and nonlinear-alg solve" for iip in (false, true), - builder in (false, true), - alg in (nothing, NewtonRaphson()) - - sccprob = dynamicss_scc_problem(iip, false) - lowered = builder ? (prob -> sccprob) : sccprob - prob = SteadyStateProblem( - iip ? f_iip : f_oop, [0.0, 0.0]; lowered_problem = lowered + # The residual solves `x` to an integer literal; the solution must keep + # `prob.u0`'s type. + @variables ex(t) + csys = complete(System([D(ex) ~ 1 - ex], t; name = :csys)) + @testset "eltype, alg=$(nameof(typeof(alg)))" for alg in ( + SSRootfind(Broyden()), DynamicSS(Tsit5()), SICNM(Rodas5P()), ) - sol = alg === nothing ? solve(prob) : solve(prob, alg) + prob = SteadyStateProblem(csys, [ex => 0.2]) + sol = 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-9 + @test sol.u isa typeof(prob.u0) + @test sol.u ≈ [1.0] end +end - # Non-SCC lowerings and problems without a lowering keep the upstream - # nonlinear-solve path. - @testset "non-SCC lowering" begin - nlprob = NonlinearProblem((u, p) -> [u[1]^2 - 4], [1.0]) - prob = SteadyStateProblem(f_oop, [0.0]; lowered_problem = nlprob) - sol = solve(prob) +# The same systems via `SCCSteadyStateProblem`: the solution is expressed on the +# compressed block layout, symbolic indexing resolves eliminated variables as +# observed, and `save_idxs` indexes that layout. +@testset "SCCSteadyStateProblem from ModelingToolkit" begin + @testset "case=$(nameof(case.sys)) alg=$(nameof(typeof(alg)))" for + case in layout_cases(), alg in case.algs + + (; sys, op, states, expected, observed) = case + prob = SciMLBase.SCCSteadyStateProblem(sys, op) + sol = solve(prob, alg; abstol = 1.0e-10, reltol = 1.0e-10) @test successful_retcode(sol) - @test sol.u ≈ [2.0] atol = 1.0e-9 + @test sol.prob === prob + @test length(sol.u) == length(state_values(prob)) + @test sol[states] ≈ expected atol = 1.0e-7 + observed === nothing || @test sol[first(observed)] ≈ last(observed) atol = 1.0e-7 + for s in variable_symbols(prob) + @test sol.u[variable_index(prob, s)] == sol[s] + end + @test length(sol.resid) == length(sol.u) + @test sol.resid ≈ zero(sol.resid) atol = 1.0e-7 - plain = SteadyStateProblem(f_oop, [0.0]) - psol = solve(plain) - @test successful_retcode(psol) - @test psol.u ≈ [1.0] atol = 1.0e-9 + isempty(state_values(prob)) && continue + saved = solve( + prob, alg; save_idxs = [1], abstol = 1.0e-10, reltol = 1.0e-10 + ) + @test saved.u ≈ sol.u[1:1] atol = 1.0e-7 + s1 = first(variable_symbols(prob)) + saved = solve( + prob, alg; save_idxs = [s1], abstol = 1.0e-10, reltol = 1.0e-10 + ) + @test saved.u ≈ [sol[s1]] atol = 1.0e-7 + @test saved.original.u ≈ sol.u atol = 1.0e-7 end - - # `remake`d values reach a callable lowering through the materialized - # problem, so the block solve sees the new operating point. - sccprob = dynamicss_scc_problem(false, false) - prob = SteadyStateProblem( - f_oop, [0.0, 0.0]; lowered_problem = prob -> sccprob - ) - @test NonlinearProblem(prob) === sccprob - prob2 = remake(prob; u0 = [5.0, 6.0]) - @test NonlinearProblem(prob2) === sccprob end @testset "SSRootfind on a ModelingToolkit SCC decomposition" begin @@ -367,71 +550,3 @@ end @test sol.original[2].original isa SteadyStateSolution @test sol.original[2].original.original isa SciMLBase.AbstractODESolution end - -@testset "SteadyStateProblem solves through an SCC lowering" begin - @variables a(t) b(t) x(t) y(t) c(t) d(t) [irreducible = true] - @named model = System( - [ - D(a) ~ 5 - 3a - b, - D(b) ~ 5 - a - 2b, - D(x) ~ a + b - x^2 - y, - D(y) ~ 3a + b - x - y^2, - D(c) ~ x + y + 7 - 2c - d, - D(d) ~ 2x + y + 11 - c - 3d, - ], t - ) - sys = mtkcompile(model) - op = [a => 0.8, b => 1.8, x => 0.8, y => 1.8, c => 2.8, d => 3.8] - prob = SteadyStateProblem(sys, op) - - @test prob.f.initialization_data.initializeprob isa SCCNonlinearProblem - - states = [a, b, x, y, c, d] - expected = [1, 2, 1, 2, 3, 4] - - @test prob.lowered_problem !== nothing - - # The `SSRootfind` solve happens on the `SCCNonlinearProblem` lowering, and - # `sol` is expressed on it. - sol = solve( - prob, SSRootfind(NewtonRaphson()); abstol = 1.0e-12, reltol = 1.0e-12 - ) - @test successful_retcode(sol) - @test sol[states] ≈ expected atol = 1.0e-9 - @test sol.prob isa SCCNonlinearProblem -end - -@testset "DynamicSS on a SteadyStateProblem with an SCC lowering" begin - # Scalar nonlinear block `x' = 3 - x^3` has a unique real root `∛3`, so the - # per-block pseudo-transient solve and the monolithic integration converge - # to the same steady state. - @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 = solve(prob, DynamicSS(); 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 sol.original isa Tuple{SciMLBase.LinearSolution, NonlinearSolution} -end - -@testset "SICNM on a SteadyStateProblem with an SCC lowering" begin - # Scalar nonlinear block `3 - x^3` has a unique real root `∛3` that is - # attracting under SICNM's continuous-Newton flow. - @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 = solve(prob, SICNM(Rodas5P()); 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 sol.original isa Tuple{SciMLBase.LinearSolution, NonlinearSolution} -end From 7d2c9aafd7bbe47f5afa9267ebbec50db4210ae7 Mon Sep 17 00:00:00 2001 From: Chris Rackauckas Date: Wed, 14 Oct 2026 22:51:56 +0000 Subject: [PATCH 02/11] Index save_idxs solutions on the saved layout; handle stateless/empty blocks Solving an SCCSteadyStateProblem with save_idxs now builds the solution on a problem whose index provider (SCCSavedStates) lists the saved states in selection order, so sol[s] resolves saved states in any order, unsaved states and observed expressions throw, parameters index as usual, and sol.original is the full solution. LinearProblem blocks are dispatched before the stateless shortcut so 0x0 blocks solve directly; a stateless nonlinear block's residual is evaluated once and must be empty. Co-Authored-By: Chris Rackauckas Co-Authored-By: Devin Agent-Harness: Devin CLI 3000.11.3 Agent-Model: fusion-claude-opus-5-5-high-sidekick-swe-2-medium Agent-Session: local headless run, log /Users/chrisrackauckas/sandbox/master/green/scc-ssprob/fu1/log.txt --- src/SteadyStateDiffEq.jl | 3 +- src/solve.jl | 98 +++++++++++++++++++++++++++++++++++----- 2 files changed, 89 insertions(+), 12 deletions(-) diff --git a/src/SteadyStateDiffEq.jl b/src/SteadyStateDiffEq.jl index c7e5818..492fae1 100644 --- a/src/SteadyStateDiffEq.jl +++ b/src/SteadyStateDiffEq.jl @@ -13,7 +13,8 @@ using SciMLBase: SciMLBase, CallbackSet, LinearProblem, NonlinearProblem, ODEPro SteadyStateProblem, SteadyStateSolution, get_du, init, isinplace, remake, solve, successful_retcode using SymbolicIndexingInterface: - NotSymbolic, parameter_values, state_values, symbolic_type, variable_index + NotSymbolic, SymbolicIndexingInterface, parameter_values, state_values, + symbolic_type, variable_index, variable_symbols const infnorm = Base.Fix2(norm, Inf) diff --git a/src/solve.jl b/src/solve.jl index a50d485..a83ea64 100644 --- a/src/solve.jl +++ b/src/solve.jl @@ -27,6 +27,25 @@ function SciMLBase.solve( return __build_ssrootfind_solution(prob, nlsol) end +# Evaluate a stateless block's residual once: out-of-place blocks return it, +# in-place blocks write into a `resid_prototype`-sized buffer. `HomotopyProblem` +# residuals take the continuation parameter `λ`, evaluated at its end. +function __scc_stateless_resid( + block::Union{NonlinearProblem, SciMLBase.HomotopyProblem}, u + ) + p = parameter_values(block) + f = block.f + args = block isa SciMLBase.HomotopyProblem ? (u, p, last(block.λspan)) : (u, p) + if isinplace(block) + fbase = f isa SciMLBase.HomotopyNonlinearFunction ? f.f : f + proto = hasproperty(fbase, :resid_prototype) ? fbase.resid_prototype : nothing + resid = proto === nothing ? similar(u) : similar(proto) + f(resid, args...) + return resid + end + return f(args...) +end + # `DynamicSS`/`SICNM` on an SCC problem solve the blocks sequentially # in SCC order: each block's parameter cache is updated from the upstream # solutions through `explicitfuns!`, `LinearProblem` blocks are solved directly, @@ -44,14 +63,7 @@ function __solve_scc_blocks( SciMLBase.invoke_with_despecialized_parameters( prob.explicitfuns![i], (parameter_values(block), sols) ) - _sol = if state_values(block) !== nothing && isempty(state_values(block)) - # A stateless block's residual is already satisfied; there is - # nothing to integrate. - SciMLBase.build_solution( - block, nothing, copy(state_values(block)), - copy(state_values(block)); retcode = ReturnCode.Success - ) - elseif block isa LinearProblem + _sol = if block isa LinearProblem A, b = block.A, block.b # `remake` recomputes `A` and `b` through a `SymbolicLinearInterface` # from the cache that `explicitfun` just updated; for a plain @@ -61,6 +73,19 @@ function __solve_scc_blocks( SciMLBase.build_linear_solution( nothing, linsol.u, resid, nothing; retcode = linsol.retcode ) + elseif state_values(block) !== nothing && isempty(state_values(block)) + # A stateless block carries no integration; its residual must be + # identically empty, which is verified by evaluating it once. + u = copy(state_values(block)) + resid = __scc_stateless_resid(block, u) + isempty(resid) || throw( + DimensionMismatch( + "SCC block $i has no states but a residual of length $(length(resid))" + ) + ) + SciMLBase.build_solution( + block, nothing, u, resid; retcode = ReturnCode.Success + ) else # `HomotopyProblem` has no `solve` preprocessing support # (`get_concrete_problem` is only defined for `NonlinearProblem` @@ -150,14 +175,56 @@ function __solve_scc_steady_state( ) end +# Index provider for the saved-state layout of an `SCCSteadyStateProblem`: +# `idxs` are the positions of the saved states in the full compressed layout. +# Parameters and independent variables forward to the wrapped `sys`; observed +# quantities are unavailable since they need the full state. +struct SCCSavedStates{S} + sys::S + idxs::Vector{Int} +end + +SymbolicIndexingInterface.symbolic_container(w::SCCSavedStates) = w.sys + +function SymbolicIndexingInterface.variable_index(w::SCCSavedStates, sym) + idx = variable_index(w.sys, sym) + idx === nothing && return nothing + if idx isa AbstractArray + out = map(i -> something(findfirst(==(i), w.idxs), 0), idx) + any(iszero, out) && return nothing + return out + end + return findfirst(==(idx), w.idxs) +end + +function SymbolicIndexingInterface.is_variable(w::SCCSavedStates, sym) + return variable_index(w, sym) !== nothing +end + +function SymbolicIndexingInterface.variable_symbols(w::SCCSavedStates) + return variable_symbols(w.sys)[w.idxs] +end + +SymbolicIndexingInterface.all_variable_symbols(w::SCCSavedStates) = + variable_symbols(w) +SymbolicIndexingInterface.is_observed(::SCCSavedStates, sym) = false + +function SymbolicIndexingInterface.all_symbols(w::SCCSavedStates) + return union( + variable_symbols(w), SymbolicIndexingInterface.parameter_symbols(w.sys), + SymbolicIndexingInterface.independent_variable_symbols(w.sys) + ) +end + # `save_idxs` is a selection in the compressed state layout: integers pass # through, symbolic entries are resolved through `prob`'s index provider and -# must be states of it. +# must be states of it. A non-array `save_idxs` is a single entry and keeps the +# scalar shape of `u`/`resid`, matching `DynamicSS` on a `SteadyStateProblem`. function __scc_save_idxs( prob::SciMLBase.SCCSteadyStateProblem, sol, save_idxs ) save_idxs === nothing && return sol - idxs = map(save_idxs) do s + function resolve(s) symbolic_type(s) === NotSymbolic() && return s idx = variable_index(prob, s) idx === nothing && throw( @@ -165,8 +232,17 @@ function __scc_save_idxs( ) return idx end + idxs = save_idxs isa AbstractArray ? map(resolve, save_idxs) : resolve(save_idxs) + saved_prob = if prob.f.sys === nothing + prob + else + SCCSteadyStateProblem( + prob.probs, prob.explicitfuns!, prob.p, prob.parameters_alias; + sys = SCCSavedStates(prob.f.sys, vec(collect(idxs))) + ) + end return SciMLBase.build_solution( - prob, sol.alg, sol.u[idxs], sol.resid[idxs]; + saved_prob, sol.alg, sol.u[idxs], sol.resid[idxs]; retcode = sol.retcode, stats = sol.stats, left = sol.left, right = sol.right, original = sol ) From b8ab9848c975cb34fa361da7f0d4241a146b29dc Mon Sep 17 00:00:00 2001 From: Chris Rackauckas Date: Wed, 14 Oct 2026 22:51:56 +0000 Subject: [PATCH 03/11] Test save_idxs layout and stateless/empty blocks; use Rodas5P in layout cases DynamicSS(Tsit5()) over the layout testset plateaued at a tolerance-sized residual noise floor at the stability-limited step size and flipped between Success and Failure under small perturbations; the L-stable Rodas5P contracts the residual with margin. Documents the save_idxs contract in the SCC docs. Co-Authored-By: Chris Rackauckas Co-Authored-By: Devin Agent-Harness: Devin CLI 3000.11.3 Agent-Model: fusion-claude-opus-5-5-high-sidekick-swe-2-medium Agent-Session: local headless run, log /Users/chrisrackauckas/sandbox/master/green/scc-ssprob/fu1/log.txt --- docs/src/index.md | 7 +++ test/scc.jl | 125 +++++++++++++++++++++++++++++++++++++++++++--- 2 files changed, 126 insertions(+), 6 deletions(-) diff --git a/docs/src/index.md b/docs/src/index.md index 2e38d25..dfef328 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -65,6 +65,13 @@ compressed SCC state layout; symbolic indexing (`sol[x]`, `sol.ps[p]`) goes through `prob.f.sys`, and `save_idxs` is a selection in that layout — integers, or states resolved through `variable_index`. +A solve with `save_idxs` returns a solution whose `prob` shares the original +problem's blocks and parameters but whose index provider lists the saved states +in selection order. `sol[s]` resolves saved states in any order; states that +were not saved and observed expressions (which need the full state) throw. +Parameters and independent variables index as usual, and `sol.original` is the +full solution over all states. + ## API ```@docs diff --git a/test/scc.jl b/test/scc.jl index c93a8b0..6c30ea7 100644 --- a/test/scc.jl +++ b/test/scc.jl @@ -1,8 +1,9 @@ -using SteadyStateDiffEq, NonlinearSolve, OrdinaryDiffEq, Test +using SteadyStateDiffEq, NonlinearSolve, OrdinaryDiffEq, LinearAlgebra, Test using ModelingToolkit using ModelingToolkit: t_nounits as t, D_nounits as D using SCCNonlinearSolve: SCCAlg -using SymbolicIndexingInterface: state_values, variable_index, variable_symbols +import SymbolicIndexingInterface +using SymbolicIndexingInterface: parameter_symbols, state_values, variable_index, variable_symbols using SciMLBase: HomotopyProblem, LinearProblem, NonlinearProblem, SCCNonlinearProblem, SteadyStateSolution @@ -349,11 +350,14 @@ function layout_cases() ) Xss, Yss = [0.35 -0.5; -0.25 1.5] \ [1.0, 1.0] - # `DynamicSS(Tsit5())` over `tspan = Inf` does not reach its terminator on - # this system — the integrator exhausts maxiters while the residual is still - # well above the tolerance — so the pseudo-transient runs use a finite span. + # `DynamicSS` integrates each block's residual to steady state. With the + # explicit `Tsit5` the step size is pinned at the stability limit near the + # fixed point, so the residual plateaus at a tolerance-sized noise floor + # right at the termination threshold and flips between `Success` and + # `Failure` under small perturbations; the L-stable `Rodas5P` grows the + # step size and contracts the residual well below it. allalgs = ( - NewtonRaphson(), SSRootfind(NewtonRaphson()), DynamicSS(Tsit5(); tspan = 100), + NewtonRaphson(), SSRootfind(NewtonRaphson()), DynamicSS(Rodas5P()), SICNM(Rodas5P()), ) nlalgs = (NewtonRaphson(), SSRootfind(NewtonRaphson())) @@ -460,6 +464,115 @@ end ) @test saved.u ≈ [sol[s1]] atol = 1.0e-7 @test saved.original.u ≈ sol.u atol = 1.0e-7 + # States not in `save_idxs` and observed expressions are unresolvable on + # the saved-state layout; parameters still index through it. + for s in Iterators.drop(variable_symbols(prob), 1) + @test_throws ErrorException saved[s] + end + observed === nothing || @test_throws ErrorException saved[first(observed)] + if !isempty(parameter_symbols(prob)) + p1 = first(parameter_symbols(prob)) + @test saved.ps[p1] == sol.ps[p1] + end + + # A symbolic `save_idxs` selection is applied in selection order. + rev = reverse(variable_symbols(prob)) + saved = solve( + prob, alg; save_idxs = rev, abstol = 1.0e-10, reltol = 1.0e-10 + ) + for s in rev + @test saved[s] ≈ sol[s] atol = 1.0e-7 + end + end +end + +@testset "SCCSteadyStateProblem save_idxs symbolic indexing" begin + # A single 2-state linear block with solution x = 1, y = 2; `sys` provides + # symbolic names so `saved[s]` resolves against the saved-state layout. + sysc = SymbolicIndexingInterface.SymbolCache([:x, :y], [], nothing) + linblock() = LinearProblem(Diagonal(ones(2)), [1.0, 2.0]; u0 = zeros(2)) + + @testset "alg=$(nameof(typeof(alg)))" for alg in ( + NewtonRaphson(), SSRootfind(NewtonRaphson()), + DynamicSS(Rodas5P()), SICNM(Rodas5P()), + ) + prob = SCCSteadyStateProblem( + (linblock(),), (Returns(nothing),), nothing, true; sys = sysc + ) + saved = solve(prob, alg; save_idxs = [2, 1]) + @test saved.u == [2, 1] + @test saved[:x] == 1 + @test saved[:y] == 2 + @test saved[[:x, :y]] == [1, 2] + @test saved.original[:x] == 1 + @test saved.original[:y] == 2 + + saved = solve(prob, alg; save_idxs = [:y]) + @test saved.u == [2] + @test saved[:y] == 2 + @test_throws ErrorException saved[:x] + @test saved.original[:x] == 1 + end +end + +@testset "stateless and empty SCC blocks" begin + @testset "alg=$(nameof(typeof(alg)))" for alg in ( + DynamicSS(Rodas5P()), SICNM(Rodas5P()), + ) + # A stateless block with a nonempty residual is rejected — the residual + # is still evaluated once. + for iip in (false, true) + counter = Ref(0) + blk = if iip + NonlinearProblem( + NonlinearFunction( + (res, u, p) -> (counter[] += 1; res .= [1.0]); + resid_prototype = [0.0] + ), Float64[] + ) + else + NonlinearProblem( + NonlinearFunction( + (u, p) -> (counter[] += 1; [1.0]); resid_prototype = [0.0] + ), Float64[] + ) + end + for T in (SCCSteadyStateProblem, SCCNonlinearProblem) + prob = T((blk,), (Returns(nothing),), nothing, true) + @test_throws DimensionMismatch solve(prob, alg) + @test counter[] > 0 + end + end + + # A stateless block whose residual is empty returns Success; the + # residual is evaluated rather than assumed zero. + counter = Ref(0) + blk = NonlinearProblem( + NonlinearFunction( + (u, p) -> (counter[] += 1; Float64[]); resid_prototype = Float64[] + ), Float64[] + ) + for T in (SCCSteadyStateProblem, SCCNonlinearProblem) + prob = T((blk,), (Returns(nothing),), nothing, true) + sol = solve(prob, alg) + @test successful_retcode(sol) + @test counter[] > 0 + end + + # An empty linear block is a no-op, alone or before another block. + lblk() = LinearProblem(zeros(0, 0), Float64[]; u0 = Float64[]) + for T in (SCCSteadyStateProblem, SCCNonlinearProblem) + sol = solve(T((lblk(),), (Returns(nothing),), nothing, true), alg) + @test successful_retcode(sol) + @test sol.u == Float64[] + prob = T( + (lblk(), LinearProblem([2.0;;], [4.0]; u0 = [0.0])), + (Returns(nothing), Returns(nothing)), nothing, true + ) + sol = solve(prob, alg) + @test successful_retcode(sol) + @test sol.u ≈ [2.0] + end end end From ed93e93a67505a8b75b6262c51784cc0f6c62a7a Mon Sep 17 00:00:00 2001 From: Chris Rackauckas Date: Wed, 14 Oct 2026 23:16:28 +0000 Subject: [PATCH 04/11] Reset the call counter per problem type in the SCC block tests Co-Authored-By: Chris Rackauckas Co-Authored-By: Devin Agent-Harness: Devin CLI 3000.11.3 (change, left uncommitted), Claude Code 2.1.285 (commit) Agent-Model: fusion-claude-opus-5-5-high-sidekick-swe-2-medium Agent-Session: local headless run (SCC stack follow-up 1); committed by local Claude session 3c311569-dda6-4687-bc9f-65febb212c33 --- test/scc.jl | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/test/scc.jl b/test/scc.jl index 6c30ea7..f77a942 100644 --- a/test/scc.jl +++ b/test/scc.jl @@ -354,8 +354,9 @@ function layout_cases() # explicit `Tsit5` the step size is pinned at the stability limit near the # fixed point, so the residual plateaus at a tolerance-sized noise floor # right at the termination threshold and flips between `Success` and - # `Failure` under small perturbations; the L-stable `Rodas5P` grows the - # step size and contracts the residual well below it. + # `Failure` under small perturbations. `Rodas5P` is not step-size limited at + # the fixed point, so the residual keeps decaying until the terminator + # fires (here at t ≈ 14). allalgs = ( NewtonRaphson(), SSRootfind(NewtonRaphson()), DynamicSS(Rodas5P()), SICNM(Rodas5P()), @@ -539,6 +540,7 @@ end end for T in (SCCSteadyStateProblem, SCCNonlinearProblem) prob = T((blk,), (Returns(nothing),), nothing, true) + counter[] = 0 @test_throws DimensionMismatch solve(prob, alg) @test counter[] > 0 end @@ -554,6 +556,7 @@ end ) for T in (SCCSteadyStateProblem, SCCNonlinearProblem) prob = T((blk,), (Returns(nothing),), nothing, true) + counter[] = 0 sol = solve(prob, alg) @test successful_retcode(sol) @test counter[] > 0 From 49c272c82b616cb08a7f105646da64af3f21b557 Mon Sep 17 00:00:00 2001 From: Chris Rackauckas Date: Thu, 15 Oct 2026 00:35:33 +0000 Subject: [PATCH 05/11] Give save_idxs solutions a problem in the saved-state layout Co-Authored-By: Chris Rackauckas Co-Authored-By: Devin Agent-Harness: Devin CLI 3000.11.3 Agent-Model: fusion-claude-opus-5-5-high-sidekick-swe-2-medium Agent-Session: local headless run, log /Users/chrisrackauckas/sandbox/master/green/scc-ssprob/fu2/log.txt --- docs/src/index.md | 15 +++++--- src/solve.jl | 96 ++++++++++++++++++++++++++++++++++++++++------- test/scc.jl | 40 +++++++++++++++++++- 3 files changed, 129 insertions(+), 22 deletions(-) diff --git a/docs/src/index.md b/docs/src/index.md index dfef328..887bc14 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -65,12 +65,15 @@ compressed SCC state layout; symbolic indexing (`sol[x]`, `sol.ps[p]`) goes through `prob.f.sys`, and `save_idxs` is a selection in that layout — integers, or states resolved through `variable_index`. -A solve with `save_idxs` returns a solution whose `prob` shares the original -problem's blocks and parameters but whose index provider lists the saved states -in selection order. `sol[s]` resolves saved states in any order; states that -were not saved and observed expressions (which need the full state) throw. -Parameters and independent variables index as usual, and `sol.original` is the -full solution over all states. +A solve with `save_idxs` returns a solution whose `prob` is a saved-state view +of the original problem: it shares the original's blocks and parameters, but +its states are the saved entries in selection order. Solving or remaking that +problem works in the saved layout — `solve(saved.prob, alg)` applies the same +selection and returns a solution over the saved states, and a `save_idxs` +passed there selects among them. `sol[s]` resolves saved states in any order; +states that were not saved and observed expressions (which need the full +state) throw. Parameters and independent variables index as usual, and +`sol.original` is the full solution over all states. ## API diff --git a/src/solve.jl b/src/solve.jl index a83ea64..38209f0 100644 --- a/src/solve.jl +++ b/src/solve.jl @@ -131,8 +131,12 @@ function SciMLBase.solve( alg === nothing && (alg = get(kwargs, :alg, nothing)) kwargs = (; (n => v for (n, v) in pairs(kwargs) if n !== :alg)...) rest = isempty(args) ? () : Base.tail(args) + # The block solves alias `prob`'s `u0`s, so the saved problem is built + # before solving: its coordinates are the saved entries' values as given. + idxs = save_idxs === nothing ? nothing : __resolve_save_idxs(prob, save_idxs) + saved_prob = idxs === nothing ? nothing : SCCSavedProblem(prob, idxs) sol = __solve_scc_steady_state(prob, alg, rest...; kwargs...) - return __scc_save_idxs(prob, sol, save_idxs) + return __scc_save_idxs(sol, idxs, saved_prob) end function __solve_scc_steady_state( @@ -216,14 +220,82 @@ function SymbolicIndexingInterface.all_symbols(w::SCCSavedStates) ) end +# A problem whose state layout is the saved-state selection of an +# `SCCSteadyStateProblem`: `idxs` are the saved positions in `prob`'s +# compressed layout. `state_values`, `solve` and `remake` all act in the saved +# layout, while the wrapped problem keeps the full blocks, parameters and +# index provider, narrowed to the saved states through `SCCSavedStates`. +struct SCCSavedProblem{ + uType, P <: SciMLBase.SCCSteadyStateProblem, I, F, + } <: SciMLBase.AbstractSteadyStateProblem{uType, false} + prob::P + # Saved entries of `prob`'s state values, snapshotted at construction: + # block `u0`s are aliased by the block solves, so a live readback would + # report post-solve values rather than the problem's own coordinates. + u0::uType + idxs::I + f::F +end + +function SCCSavedProblem(prob::SciMLBase.SCCSteadyStateProblem, idxs) + f = if prob.f.sys === nothing + SciMLBase.NonlinearFunction{false}(Returns(nothing)) + else + SciMLBase.NonlinearFunction{false}( + Returns(nothing); sys = SCCSavedStates(prob.f.sys, vec(collect(idxs))) + ) + end + u = state_values(prob)[idxs] + u isa AbstractArray && (u = copy(u)) + return SCCSavedProblem{typeof(u), typeof(prob), typeof(idxs), typeof(f)}( + prob, u, idxs, f + ) +end + +__saved_positions(w::SCCSavedProblem) = w.idxs isa Integer ? (w.idxs,) : w.idxs + +SymbolicIndexingInterface.state_values(w::SCCSavedProblem) = w.u0 +SymbolicIndexingInterface.parameter_values(w::SCCSavedProblem) = + parameter_values(w.prob) +SymbolicIndexingInterface.symbolic_container(w::SCCSavedProblem) = w.f +function SymbolicIndexingInterface.set_state!(w::SCCSavedProblem, val, idx) + set_state!(w.prob, val, __saved_positions(w)[idx]) + w.u0 isa AbstractArray && (w.u0[idx] = val) + return val +end + +function SciMLBase.solve( + w::SCCSavedProblem, args...; save_idxs = nothing, kwargs... + ) + composed = if save_idxs === nothing + w.idxs + else + __saved_positions(w)[__resolve_save_idxs(w, save_idxs)] + end + return solve(w.prob, args...; save_idxs = composed, kwargs...) +end + +function SciMLBase.remake( + w::SCCSavedProblem; u0 = missing, p = missing, kwargs... + ) + if u0 isa Number || (u0 isa AbstractArray && eltype(u0) <: Number) + full = collect(state_values(w.prob)) + positions = collect(__saved_positions(w)) + if u0 isa Number + full[positions] .= u0 + else + full[positions] = u0 + end + u0 = full + end + return SCCSavedProblem(remake(w.prob; u0, p, kwargs...), w.idxs) +end + # `save_idxs` is a selection in the compressed state layout: integers pass # through, symbolic entries are resolved through `prob`'s index provider and # must be states of it. A non-array `save_idxs` is a single entry and keeps the # scalar shape of `u`/`resid`, matching `DynamicSS` on a `SteadyStateProblem`. -function __scc_save_idxs( - prob::SciMLBase.SCCSteadyStateProblem, sol, save_idxs - ) - save_idxs === nothing && return sol +function __resolve_save_idxs(prob, save_idxs) function resolve(s) symbolic_type(s) === NotSymbolic() && return s idx = variable_index(prob, s) @@ -232,15 +304,11 @@ function __scc_save_idxs( ) return idx end - idxs = save_idxs isa AbstractArray ? map(resolve, save_idxs) : resolve(save_idxs) - saved_prob = if prob.f.sys === nothing - prob - else - SCCSteadyStateProblem( - prob.probs, prob.explicitfuns!, prob.p, prob.parameters_alias; - sys = SCCSavedStates(prob.f.sys, vec(collect(idxs))) - ) - end + return save_idxs isa AbstractArray ? map(resolve, save_idxs) : resolve(save_idxs) +end + +function __scc_save_idxs(sol, idxs, saved_prob) + idxs === nothing && return sol return SciMLBase.build_solution( saved_prob, sol.alg, sol.u[idxs], sol.resid[idxs]; retcode = sol.retcode, stats = sol.stats, left = sol.left, diff --git a/test/scc.jl b/test/scc.jl index f77a942..fba25bd 100644 --- a/test/scc.jl +++ b/test/scc.jl @@ -253,7 +253,10 @@ end ) @test saved.u ≈ full.u[[2, 4]] @test saved.resid ≈ full.resid[[2, 4]] - @test saved.prob === prob + @test state_values(saved.prob) == state_values(prob)[[2, 4]] + @test solve( + saved.prob, alg; abstol = 1.0e-10, reltol = 1.0e-10 + ).u ≈ saved.u @test saved.original.u ≈ full.u end @@ -484,6 +487,11 @@ end for s in rev @test saved[s] ≈ sol[s] atol = 1.0e-7 end + re = solve(saved.prob, alg; abstol = 1.0e-10, reltol = 1.0e-10) + for s in rev + @test re[s] ≈ sol[s] atol = 1.0e-7 + end + @test re.u ≈ saved.u end end @@ -491,7 +499,7 @@ end # A single 2-state linear block with solution x = 1, y = 2; `sys` provides # symbolic names so `saved[s]` resolves against the saved-state layout. sysc = SymbolicIndexingInterface.SymbolCache([:x, :y], [], nothing) - linblock() = LinearProblem(Diagonal(ones(2)), [1.0, 2.0]; u0 = zeros(2)) + linblock() = LinearProblem(Diagonal(ones(2)), [1.0, 2.0]; u0 = [10.0, 20.0]) @testset "alg=$(nameof(typeof(alg)))" for alg in ( NewtonRaphson(), SSRootfind(NewtonRaphson()), @@ -508,11 +516,39 @@ end @test saved.original[:x] == 1 @test saved.original[:y] == 2 + # Re-solving the saved problem works in the saved layout: the same + # selection is applied again, and `u0`/`save_idxs` are interpreted in + # saved positions. + re2 = solve(saved.prob, alg) + @test re2.u == [2, 1] + @test re2[:x] == 1 + @test re2[:y] == 2 + @test solve(saved.prob, alg; save_idxs = [:y]).u == [2] + @test state_values(saved.prob) == [20.0, 10.0] + @test SymbolicIndexingInterface.getu(saved.prob, :y)(saved.prob) == 20.0 + + pre_remake = copy(state_values(prob)) + r = remake(saved.prob; u0 = [30.0, 40.0]) + @test state_values(r) == [30.0, 40.0] + @test state_values(prob) == pre_remake + @test solve(r, alg).u == [2, 1] + + # Fresh problem: the saved problem's coordinates are the guesses as + # given, snapshotted before the solve updates the block `u0`s. + prob = SCCSteadyStateProblem( + (linblock(),), (Returns(nothing),), nothing, true; sys = sysc + ) saved = solve(prob, alg; save_idxs = [:y]) @test saved.u == [2] @test saved[:y] == 2 @test_throws ErrorException saved[:x] @test saved.original[:x] == 1 + + re = solve(saved.prob, alg) + @test re.u == [2] + @test re[:y] == 2 + @test_throws ErrorException re[:x] + @test SymbolicIndexingInterface.getu(saved.prob, :y)(saved.prob) == 20.0 end end From 33645eb7cc4948b524cd19a2a5943e50cdc6f556 Mon Sep 17 00:00:00 2001 From: Chris Rackauckas Date: Thu, 15 Oct 2026 00:35:41 +0000 Subject: [PATCH 06/11] Keep solving SteadyStateProblem lowerings until the field is removed Co-Authored-By: Chris Rackauckas Co-Authored-By: Devin Agent-Harness: Devin CLI 3000.11.3 Agent-Model: fusion-claude-opus-5-5-high-sidekick-swe-2-medium Agent-Session: local headless run, log /Users/chrisrackauckas/sandbox/master/green/scc-ssprob/fu2/log.txt --- src/algorithms.jl | 17 +++++++-- src/solve.jl | 85 ++++++++++++++++++++++++++++++++++++++++++-- test/scc.jl | 90 +++++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 188 insertions(+), 4 deletions(-) diff --git a/src/algorithms.jl b/src/algorithms.jl index c28f2d6..0760f69 100644 --- a/src/algorithms.jl +++ b/src/algorithms.jl @@ -97,6 +97,11 @@ An `SCCSteadyStateProblem` (for example the SCC decomposition built by `SCCNonlinearProblem(prob)`: it is forwarded to `alg` unchanged, and the solution is expressed on the SCC problem's compressed state ordering. +A `SteadyStateProblem` that records a `lowered_problem` (for example the SCC +decomposition stored by `ModelingToolkit`) is solved through that lowering; +since the lowering defines its own state ordering, the returned solution is +expressed on the lowered problem. + # Arguments - `alg`: the nonlinear solver algorithm passed to `solve`. When `alg === nothing`, @@ -151,7 +156,11 @@ nonlinear block's residual to be attracting under its own pseudo-transient dynamics, not a globally attracting concatenated field. An `SCCSteadyStateProblem` (for example one built by `ModelingToolkit.SCCSteadyStateProblem`) takes the same sequential solve, and -its solution is expressed on the compressed state ordering. +its solution is expressed on the compressed state ordering. A +`SteadyStateProblem` that records an `SCCNonlinearProblem` as its +`lowered_problem` (for example one built by `ModelingToolkit`) takes the same +sequential solve, and its solution is expressed on the lowered problem's state +ordering. # Arguments @@ -253,7 +262,11 @@ continuous-Newton flow on the block residual — so convergence only requires ea nonlinear block to be tractable for SICNM on its own, not the concatenated system. An `SCCSteadyStateProblem` (for example one built by `ModelingToolkit.SCCSteadyStateProblem`) takes the same sequential solve, and -its solution is expressed on the compressed state ordering. +its solution is expressed on the compressed state ordering. A +`SteadyStateProblem` that records an `SCCNonlinearProblem` as its +`lowered_problem` (for example one built by `ModelingToolkit`) takes the same +sequential solve, and its solution is expressed on the lowered problem's state +ordering. # Arguments diff --git a/src/solve.jl b/src/solve.jl index 38209f0..b6fcd7d 100644 --- a/src/solve.jl +++ b/src/solve.jl @@ -10,8 +10,19 @@ function SciMLBase.__solve( args...; kwargs... ) nlprob = NonlinearProblem(prob) - nlsol = solve(nlprob, alg.alg, args...; kwargs...) - return __build_ssrootfind_solution(prob, nlsol) + fwd = kwargs + if nlprob isa SciMLBase.SCCNonlinearProblem + # The nonlinear solve pipeline injects `alias`/`verbose` as nonlinear + # specifier types, which `LinearProblem` blocks in the SCC solve do not + # understand. + fwd = (; (n => v for (n, v) in pairs(kwargs) if n !== :alias && n !== :verbose)...) + end + nlsol = solve(nlprob, alg.alg, args...; fwd...) + # A stored `lowered_problem` (e.g. an `SCCNonlinearProblem`) solves in the + # lowering's own state ordering, so its solution is expressed on the + # lowering rather than the steady-state problem. + solprob = __lowered_problem(prob) === nothing ? prob : nlprob + return __build_ssrootfind_solution(solprob, nlsol) end @@ -316,6 +327,64 @@ function __scc_save_idxs(sol, idxs, saved_prob) ) end +# `SteadyStateProblem.lowered_problem` is a SciMLBase 3.x field; reads go +# through `hasfield` so this package also loads where the field is gone. +function __lowered_problem(prob) + return prob isa SteadyStateProblem && hasfield(typeof(prob), :lowered_problem) ? + prob.lowered_problem : nothing +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... + ) + lp = __lowered_problem(prob) + lp === nothing && return nothing + lp isa SciMLBase.AbstractSciMLProblem || (lp = lp(prob)) + lp isa SciMLBase.SCCNonlinearProblem || return nothing + sccsol = solve(lp, alg, args...; kwargs...) + save_idxs === nothing && return sccsol + return SciMLBase.build_solution( + lp, sccsol.alg, sccsol.u[save_idxs], sccsol.resid[save_idxs]; + retcode = sccsol.retcode, original = sccsol + ) +end + +# A `SteadyStateProblem`'s stored SCC lowering materializes through +# `SciMLBase.NonlinearProblem(prob)` during `solve`. When it is an +# `SCCNonlinearProblem`, the generic nonlinear-solve path re-dispatches +# `__solve` on it and the fallbacks try to convert it back via `prob.u0` — a +# field it does not have — so the lowering goes through its own `solve` +# dispatch instead, which routes bare and `AbstractNonlinearAlgorithm` solves +# to the SCC solver (`SCCNonlinearSolve.jl`, loaded downstream). Algorithms +# that do not lower to a nonlinear solve (`DynamicSS`, ODE algorithms) keep +# their own `__solve` dispatch; non-SCC lowerings keep the upstream path +# verbatim via `invoke`. +function SciMLBase.solve(prob::SteadyStateProblem, args...; kwargs...) + alg = isempty(args) ? nothing : first(args) + alg === nothing && (alg = get(kwargs, :alg, nothing)) + alg === nothing && prob.kwargs !== nothing && + (alg = get(prob.kwargs, :alg, nothing)) + if alg !== nothing && !(alg isa SciMLBase.AbstractNonlinearAlgorithm) + return invoke( + solve, Tuple{SciMLBase.AbstractNonlinearProblem, Vararg{Any}}, + prob, args...; kwargs... + ) + end + # Concretize like `NonlinearSolveBase.solve_up` (evaluates `u0(p, Inf)`-style + # `u0`s, promotes integer `u0`s, applies `u0`/`p` solve kwargs) so a + # callable `lowered_problem` sees the updated operating point. + _prob = NonlinearSolveBase.get_concrete_problem(prob; kwargs...) + nlprob = NonlinearProblem(_prob) + nlprob isa SciMLBase.SCCNonlinearProblem && return solve(nlprob, args...; kwargs...) + return invoke( + solve, Tuple{SciMLBase.AbstractNonlinearProblem, Vararg{Any}}, + _prob, args...; kwargs... + ) +end + __get_tspan(u0, alg::Union{DynamicSS, SICNM}) = __get_tspan(u0, alg.tspan) __get_tspan(u0, tspan::Tuple) = tspan function __get_tspan(u0, tspan::Number) @@ -334,6 +403,12 @@ 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, + termination_condition, alias, save_idxs, kwargs... + ) + sccsol !== nothing && return sccsol + tspan = __get_tspan(prob.u0, alg) f = if prob isa SteadyStateProblem @@ -438,6 +513,12 @@ function SciMLBase.__solve( 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 + prob.u0 isa AbstractVector || throw(ArgumentError("SICNM currently only supports `AbstractVector` initial conditions")) tspan = __get_tspan(prob.u0, alg) diff --git a/test/scc.jl b/test/scc.jl index fba25bd..41c4526 100644 --- a/test/scc.jl +++ b/test/scc.jl @@ -702,3 +702,93 @@ end @test sol.original[2].original isa SteadyStateSolution @test sol.original[2].original.original isa SciMLBase.AbstractODESolution end +@testset "SteadyStateProblem with an SCC lowering" begin + f_iip(du, u, p, t) = (du .= 1 .- u) + f_oop(u, p, t) = 1 .- u + @testset "iip=$iip builder=$builder alg=$alg" for iip in (false, true), + builder in (false, true), + alg in ( + DynamicSS(Tsit5()), + SICNM(Rodas5P()), + SSRootfind(NewtonRaphson()), + SSRootfind(SCCAlg(; nlalg = NewtonRaphson())), + ) + + sccprob = dynamicss_scc_problem(iip, false) + lowered = builder ? (prob -> sccprob) : sccprob + prob = SteadyStateProblem( + iip ? f_iip : f_oop, [0.0, 0.0]; lowered_problem = lowered + ) + sol = 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 + @test sol.original !== nothing + end + + # Bare and `AbstractNonlinearAlgorithm` solves have no `__solve` dispatch of + # their own here; they land on the generic nonlinear-solve path, which + # converts back through `prob.u0` — absent on `SCCNonlinearProblem`. `solve` + # forwards the lowering to its own `solve` dispatch (the SCC solver) + # instead. + @testset "bare and nonlinear-alg solve" for iip in (false, true), + builder in (false, true), + alg in (nothing, NewtonRaphson()) + + sccprob = dynamicss_scc_problem(iip, false) + lowered = builder ? (prob -> sccprob) : sccprob + prob = SteadyStateProblem( + iip ? f_iip : f_oop, [0.0, 0.0]; lowered_problem = lowered + ) + sol = alg === nothing ? solve(prob) : solve(prob, alg) + @test successful_retcode(sol) + @test sol.u ≈ [1, 2, 1, 2] atol = 1.0e-9 + end + + # Non-SCC lowerings and problems without a lowering keep the upstream + # nonlinear-solve path. + @testset "non-SCC lowering" begin + nlprob = NonlinearProblem((u, p) -> [u[1]^2 - 4], [1.0]) + prob = SteadyStateProblem(f_oop, [0.0]; lowered_problem = nlprob) + sol = solve(prob) + @test successful_retcode(sol) + @test sol.u ≈ [2.0] atol = 1.0e-9 + + plain = SteadyStateProblem(f_oop, [0.0]) + psol = solve(plain) + @test successful_retcode(psol) + @test psol.u ≈ [1.0] atol = 1.0e-9 + end + + # `remake`d values reach a callable lowering through the materialized + # problem, so the block solve sees the new operating point. + sccprob = dynamicss_scc_problem(false, false) + prob = SteadyStateProblem( + f_oop, [0.0, 0.0]; lowered_problem = prob -> sccprob + ) + @test NonlinearProblem(prob) === sccprob + prob2 = remake(prob; u0 = [5.0, 6.0]) + @test NonlinearProblem(prob2) === sccprob +end + +# A `SteadyStateProblem` may record a non-SCC `lowered_problem`; it solves +# through the lowering and the solution is expressed on the lowered problem's +# own coordinates. +@testset "SteadyStateProblem with a generic lowering" begin + ssys = SymbolicIndexingInterface.SymbolCache([:x, :y, :z], [], nothing) + lowsys = SymbolicIndexingInterface.SymbolCache([:y], [], nothing) + low = NonlinearProblem( + SciMLBase.NonlinearFunction((u, p) -> [2 - u[1]]; sys = lowsys), [0.0] + ) + for lp in (low, prob -> low), alg in (SSRootfind(NewtonRaphson()), SSRootfind()) + prob = SteadyStateProblem( + SciMLBase.ODEFunction((u, p, t) -> [1 - u[1], 2 - u[2], 3 - u[3]]; sys = ssys), + zeros(3); lowered_problem = lp + ) + sol = solve(prob, alg) + @test sol.u ≈ [2] + @test variable_symbols(sol) == [:y] + @test sol[:y] ≈ 2 + @test sol.prob === low + end +end From 69c9c1d53d4cebed71c8c4e3933e132692ed47cd Mon Sep 17 00:00:00 2001 From: Chris Rackauckas Date: Thu, 15 Oct 2026 00:39:20 +0000 Subject: [PATCH 07/11] Import set_state! for the saved-state problem Co-Authored-By: Chris Rackauckas Co-Authored-By: Devin Agent-Harness: Devin CLI 3000.11.3 Agent-Model: fusion-claude-opus-5-5-high-sidekick-swe-2-medium Agent-Session: local headless run, log /Users/chrisrackauckas/sandbox/master/green/scc-ssprob/fu2/log.txt --- src/SteadyStateDiffEq.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/SteadyStateDiffEq.jl b/src/SteadyStateDiffEq.jl index 492fae1..9a645b4 100644 --- a/src/SteadyStateDiffEq.jl +++ b/src/SteadyStateDiffEq.jl @@ -13,8 +13,8 @@ using SciMLBase: SciMLBase, CallbackSet, LinearProblem, NonlinearProblem, ODEPro SteadyStateProblem, SteadyStateSolution, get_du, init, isinplace, remake, solve, successful_retcode using SymbolicIndexingInterface: - NotSymbolic, SymbolicIndexingInterface, parameter_values, state_values, - symbolic_type, variable_index, variable_symbols + NotSymbolic, SymbolicIndexingInterface, parameter_values, set_state!, + state_values, symbolic_type, variable_index, variable_symbols const infnorm = Base.Fix2(norm, Inf) From 5edd7c30600f479aad39dbdaabdf624af3c73075 Mon Sep 17 00:00:00 2001 From: Chris Rackauckas Date: Thu, 15 Oct 2026 00:52:05 +0000 Subject: [PATCH 08/11] Make the saved-state problem a live view of the full problem Co-Authored-By: Chris Rackauckas Co-Authored-By: Devin Agent-Harness: Devin CLI 3000.11.3 Agent-Model: fusion-claude-opus-5-5-high-sidekick-swe-2-medium Agent-Session: local headless run, log /Users/chrisrackauckas/sandbox/master/green/scc-ssprob/fu2/log.txt --- src/solve.jl | 44 ++++++++++++++++++-------------------------- test/scc.jl | 15 +++++++-------- 2 files changed, 25 insertions(+), 34 deletions(-) diff --git a/src/solve.jl b/src/solve.jl index b6fcd7d..abaffe8 100644 --- a/src/solve.jl +++ b/src/solve.jl @@ -142,12 +142,8 @@ function SciMLBase.solve( alg === nothing && (alg = get(kwargs, :alg, nothing)) kwargs = (; (n => v for (n, v) in pairs(kwargs) if n !== :alg)...) rest = isempty(args) ? () : Base.tail(args) - # The block solves alias `prob`'s `u0`s, so the saved problem is built - # before solving: its coordinates are the saved entries' values as given. - idxs = save_idxs === nothing ? nothing : __resolve_save_idxs(prob, save_idxs) - saved_prob = idxs === nothing ? nothing : SCCSavedProblem(prob, idxs) sol = __solve_scc_steady_state(prob, alg, rest...; kwargs...) - return __scc_save_idxs(sol, idxs, saved_prob) + return __scc_save_idxs(prob, sol, save_idxs) end function __solve_scc_steady_state( @@ -231,19 +227,15 @@ function SymbolicIndexingInterface.all_symbols(w::SCCSavedStates) ) end -# A problem whose state layout is the saved-state selection of an -# `SCCSteadyStateProblem`: `idxs` are the saved positions in `prob`'s -# compressed layout. `state_values`, `solve` and `remake` all act in the saved -# layout, while the wrapped problem keeps the full blocks, parameters and -# index provider, narrowed to the saved states through `SCCSavedStates`. +# A live view of an `SCCSteadyStateProblem` whose state layout is the +# saved-state selection: `idxs` are the saved positions in `prob`'s compressed +# layout. `state_values`, `solve` and `remake` all act in the saved layout; +# the view shares the wrapped problem's blocks, parameters and index provider, +# narrowed to the saved states through `SCCSavedStates`. struct SCCSavedProblem{ uType, P <: SciMLBase.SCCSteadyStateProblem, I, F, } <: SciMLBase.AbstractSteadyStateProblem{uType, false} prob::P - # Saved entries of `prob`'s state values, snapshotted at construction: - # block `u0`s are aliased by the block solves, so a live readback would - # report post-solve values rather than the problem's own coordinates. - u0::uType idxs::I f::F end @@ -256,23 +248,20 @@ function SCCSavedProblem(prob::SciMLBase.SCCSteadyStateProblem, idxs) Returns(nothing); sys = SCCSavedStates(prob.f.sys, vec(collect(idxs))) ) end - u = state_values(prob)[idxs] - u isa AbstractArray && (u = copy(u)) - return SCCSavedProblem{typeof(u), typeof(prob), typeof(idxs), typeof(f)}( - prob, u, idxs, f - ) + return SCCSavedProblem{ + typeof(state_values(prob)[idxs]), typeof(prob), typeof(idxs), typeof(f), + }(prob, idxs, f) end __saved_positions(w::SCCSavedProblem) = w.idxs isa Integer ? (w.idxs,) : w.idxs -SymbolicIndexingInterface.state_values(w::SCCSavedProblem) = w.u0 +SymbolicIndexingInterface.state_values(w::SCCSavedProblem) = + state_values(w.prob)[w.idxs] SymbolicIndexingInterface.parameter_values(w::SCCSavedProblem) = parameter_values(w.prob) SymbolicIndexingInterface.symbolic_container(w::SCCSavedProblem) = w.f function SymbolicIndexingInterface.set_state!(w::SCCSavedProblem, val, idx) - set_state!(w.prob, val, __saved_positions(w)[idx]) - w.u0 isa AbstractArray && (w.u0[idx] = val) - return val + return set_state!(w.prob, val, __saved_positions(w)[idx]) end function SciMLBase.solve( @@ -318,10 +307,13 @@ function __resolve_save_idxs(prob, save_idxs) return save_idxs isa AbstractArray ? map(resolve, save_idxs) : resolve(save_idxs) end -function __scc_save_idxs(sol, idxs, saved_prob) - idxs === nothing && return sol +function __scc_save_idxs( + prob::SciMLBase.SCCSteadyStateProblem, sol, save_idxs + ) + save_idxs === nothing && return sol + idxs = __resolve_save_idxs(prob, save_idxs) return SciMLBase.build_solution( - saved_prob, sol.alg, sol.u[idxs], sol.resid[idxs]; + SCCSavedProblem(prob, idxs), sol.alg, sol.u[idxs], sol.resid[idxs]; retcode = sol.retcode, stats = sol.stats, left = sol.left, right = sol.right, original = sol ) diff --git a/test/scc.jl b/test/scc.jl index 41c4526..169de2b 100644 --- a/test/scc.jl +++ b/test/scc.jl @@ -524,8 +524,11 @@ end @test re2[:x] == 1 @test re2[:y] == 2 @test solve(saved.prob, alg; save_idxs = [:y]).u == [2] - @test state_values(saved.prob) == [20.0, 10.0] - @test SymbolicIndexingInterface.getu(saved.prob, :y)(saved.prob) == 20.0 + @test state_values(saved.prob) == state_values(prob)[[2, 1]] + @test SymbolicIndexingInterface.getu(saved.prob, :y)(saved.prob) == + SymbolicIndexingInterface.getu(prob, :y)(prob) + @test SymbolicIndexingInterface.getu(saved.prob, :x)(saved.prob) == + SymbolicIndexingInterface.getu(prob, :x)(prob) pre_remake = copy(state_values(prob)) r = remake(saved.prob; u0 = [30.0, 40.0]) @@ -533,11 +536,6 @@ end @test state_values(prob) == pre_remake @test solve(r, alg).u == [2, 1] - # Fresh problem: the saved problem's coordinates are the guesses as - # given, snapshotted before the solve updates the block `u0`s. - prob = SCCSteadyStateProblem( - (linblock(),), (Returns(nothing),), nothing, true; sys = sysc - ) saved = solve(prob, alg; save_idxs = [:y]) @test saved.u == [2] @test saved[:y] == 2 @@ -548,7 +546,8 @@ end @test re.u == [2] @test re[:y] == 2 @test_throws ErrorException re[:x] - @test SymbolicIndexingInterface.getu(saved.prob, :y)(saved.prob) == 20.0 + @test SymbolicIndexingInterface.getu(saved.prob, :y)(saved.prob) == + SymbolicIndexingInterface.getu(prob, :y)(prob) end end From bd2e6fe6f12bc9b56c6cbdc9236e8a51c99c843c Mon Sep 17 00:00:00 2001 From: Chris Rackauckas Date: Thu, 15 Oct 2026 00:52:44 +0000 Subject: [PATCH 09/11] Separate the SCC lowering testset with a blank line Co-Authored-By: Chris Rackauckas Co-Authored-By: Devin Agent-Harness: Devin CLI 3000.11.3 Agent-Model: fusion-claude-opus-5-5-high-sidekick-swe-2-medium Agent-Session: local headless run, log /Users/chrisrackauckas/sandbox/master/green/scc-ssprob/fu2/log.txt --- test/scc.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/test/scc.jl b/test/scc.jl index 169de2b..7611062 100644 --- a/test/scc.jl +++ b/test/scc.jl @@ -701,6 +701,7 @@ end @test sol.original[2].original isa SteadyStateSolution @test sol.original[2].original.original isa SciMLBase.AbstractODESolution end + @testset "SteadyStateProblem with an SCC lowering" begin f_iip(du, u, p, t) = (du .= 1 .- u) f_oop(u, p, t) = 1 .- u From 86e034a5b5d911c51a03a481c544458dc764ca2f Mon Sep 17 00:00:00 2001 From: Chris Rackauckas Date: Thu, 15 Oct 2026 09:42:40 +0000 Subject: [PATCH 10/11] Forward parameters, guess-less states and promoted remakes through saved problems The saved-state view of an SCCSteadyStateProblem forwarded parameter_values but not set_parameter!, so setp on it updated only the top-level parameter object and left unaliased block caches stale: re-solving returned y = 2 while sol.ps[:k] == 3. set_parameter! and finalize_parameters_hook! now forward to the wrapped problem. Building the view indexed state_values(prob) unconditionally, so save_idxs threw MethodError getindex(::Nothing, ...) on all-LinearProblem problems without initial guesses. Such views have uType Nothing and state_values nothing; remaking one with a saved-layout u0 fills the unsaved positions with zeros. remake copied the old state vector and assigned into it, truncating BigFloat replacements to Float64 and throwing on ForwardDiff duals. The expanded state now uses the promoted element type. Co-Authored-By: Chris Rackauckas Co-Authored-By: Devin Agent-Harness: Devin CLI 3000.11.92 Agent-Model: fusion-claude-opus-5-5-high-sidekick-swe-2-medium Agent-Session: local headless run, log /Users/chrisrackauckas/sandbox/master/green/fu/SciML_SteadyStateDiffEq.jl-178-0909/log.txt --- src/SteadyStateDiffEq.jl | 3 +- src/solve.jl | 32 ++++++++++--- test/scc.jl | 98 +++++++++++++++++++++++++++++++++++++++- 3 files changed, 125 insertions(+), 8 deletions(-) diff --git a/src/SteadyStateDiffEq.jl b/src/SteadyStateDiffEq.jl index 9a645b4..0059871 100644 --- a/src/SteadyStateDiffEq.jl +++ b/src/SteadyStateDiffEq.jl @@ -13,7 +13,8 @@ using SciMLBase: SciMLBase, CallbackSet, LinearProblem, NonlinearProblem, ODEPro SteadyStateProblem, SteadyStateSolution, get_du, init, isinplace, remake, solve, successful_retcode using SymbolicIndexingInterface: - NotSymbolic, SymbolicIndexingInterface, parameter_values, set_state!, + NotSymbolic, SymbolicIndexingInterface, finalize_parameters_hook!, + parameter_values, set_parameter!, set_state!, state_values, symbolic_type, variable_index, variable_symbols const infnorm = Base.Fix2(norm, Inf) diff --git a/src/solve.jl b/src/solve.jl index abaffe8..50a6a22 100644 --- a/src/solve.jl +++ b/src/solve.jl @@ -248,21 +248,30 @@ function SCCSavedProblem(prob::SciMLBase.SCCSteadyStateProblem, idxs) Returns(nothing); sys = SCCSavedStates(prob.f.sys, vec(collect(idxs))) ) end - return SCCSavedProblem{ - typeof(state_values(prob)[idxs]), typeof(prob), typeof(idxs), typeof(f), - }(prob, idxs, f) + uType = state_values(prob) === nothing ? Nothing : typeof(state_values(prob)[idxs]) + return SCCSavedProblem{uType, typeof(prob), typeof(idxs), typeof(f)}( + prob, idxs, f + ) end __saved_positions(w::SCCSavedProblem) = w.idxs isa Integer ? (w.idxs,) : w.idxs -SymbolicIndexingInterface.state_values(w::SCCSavedProblem) = - state_values(w.prob)[w.idxs] +function SymbolicIndexingInterface.state_values(w::SCCSavedProblem) + states = state_values(w.prob) + return states === nothing ? nothing : states[w.idxs] +end SymbolicIndexingInterface.parameter_values(w::SCCSavedProblem) = parameter_values(w.prob) SymbolicIndexingInterface.symbolic_container(w::SCCSavedProblem) = w.f function SymbolicIndexingInterface.set_state!(w::SCCSavedProblem, val, idx) return set_state!(w.prob, val, __saved_positions(w)[idx]) end +function SymbolicIndexingInterface.set_parameter!(w::SCCSavedProblem, val, idx) + return set_parameter!(w.prob, val, idx) +end +function SymbolicIndexingInterface.finalize_parameters_hook!(w::SCCSavedProblem, sym) + return finalize_parameters_hook!(w.prob, sym) +end function SciMLBase.solve( w::SCCSavedProblem, args...; save_idxs = nothing, kwargs... @@ -279,7 +288,18 @@ function SciMLBase.remake( w::SCCSavedProblem; u0 = missing, p = missing, kwargs... ) if u0 isa Number || (u0 isa AbstractArray && eltype(u0) <: Number) - full = collect(state_values(w.prob)) + old = state_values(w.prob) + T = promote_type( + old === nothing ? Union{} : eltype(old), + u0 isa Number ? typeof(u0) : eltype(u0) + ) + full = if old === nothing + # A guess-less problem has no values at the unsaved positions; + # zeros match the guess `LinearSolve` uses when none is given. + zeros(T, sum(size(b.A, 2) for b in w.prob.probs)) + else + copyto!(Vector{T}(undef, length(old)), old) + end positions = collect(__saved_positions(w)) if u0 isa Number full[positions] .= u0 diff --git a/test/scc.jl b/test/scc.jl index 7611062..a4ab9c9 100644 --- a/test/scc.jl +++ b/test/scc.jl @@ -1,9 +1,12 @@ using SteadyStateDiffEq, NonlinearSolve, OrdinaryDiffEq, LinearAlgebra, Test +using ForwardDiff using ModelingToolkit using ModelingToolkit: t_nounits as t, D_nounits as D using SCCNonlinearSolve: SCCAlg import SymbolicIndexingInterface -using SymbolicIndexingInterface: parameter_symbols, state_values, variable_index, variable_symbols +using SymbolicIndexingInterface: + parameter_symbols, parameter_values, setp, state_values, variable_index, + variable_symbols using SciMLBase: HomotopyProblem, LinearProblem, NonlinearProblem, SCCNonlinearProblem, SteadyStateSolution @@ -551,6 +554,99 @@ end end end +@testset "saved-problem parameter updates reach the blocks" begin + blk(c, s) = NonlinearProblem( + NonlinearFunction((u, p) -> [c * p[1] - u[1]]; sys = SymbolicIndexingInterface.SymbolCache([s], [:k])), + [0.0], [1.0] + ) + mkprob() = SCCSteadyStateProblem( + (blk(1, :x), blk(2, :y)), (Returns(nothing), Returns(nothing)), + [1.0], false; sys = SymbolicIndexingInterface.SymbolCache([:x, :y], [:k]) + ) + @testset "alg=$(nameof(typeof(alg)))" for alg in ( + NewtonRaphson(), SSRootfind(NewtonRaphson()), + DynamicSS(Rodas5P()), SICNM(Rodas5P()), + ) + prob = mkprob() + setp(prob, :k)(prob, 3.0) + @test solve(prob, alg; abstol = 1.0e-10, reltol = 1.0e-10).u ≈ [3, 6] atol = 1.0e-6 + + prob = mkprob() + saved = solve(prob, alg; save_idxs = [:y], abstol = 1.0e-10, reltol = 1.0e-10) + @test saved.u ≈ [2] atol = 1.0e-6 + setp(saved.prob, :k)(saved.prob, 3.0) + @test all(b -> parameter_values(b) == [3.0], saved.prob.prob.probs) + re = solve(saved.prob, alg; abstol = 1.0e-10, reltol = 1.0e-10) + @test re.ps[:k] == 3 + @test re.u ≈ [6.0] atol = 1.0e-6 + @test re[:y] ≈ 6 atol = 1.0e-6 + end +end + +@testset "guess-less blocks support save_idxs" begin + @testset "alg=$(nameof(typeof(alg)))" for alg in ( + NewtonRaphson(), SSRootfind(NewtonRaphson()), + DynamicSS(Rodas5P()), SICNM(Rodas5P()), + ) + prob = SCCSteadyStateProblem( + (LinearProblem(Diagonal(ones(2)), [1.0, 2.0]),), + (Returns(nothing),); sys = SymbolicIndexingInterface.SymbolCache([:x, :y], []) + ) + @test state_values(prob) === nothing + @test solve(prob, alg).u == [1, 2] + + for (sel, expected) in ( + ([2], [2]), ([:y], [2]), ([2, 1], [2, 1]), + ([:y, :x], [2, 1]), (2, 2), (:y, 2), + ) + saved = solve(prob, alg; save_idxs = sel) + @test saved.u == expected + @test saved[:y] == 2 + @test state_values(saved.prob) === nothing + @test solve(saved.prob, alg).u == expected + if sel isa AbstractArray + @test solve(saved.prob, alg; save_idxs = [:y]).u == [2] + end + end + + saved = solve(prob, alg; save_idxs = [2, 1]) + r = remake(saved.prob; u0 = [5.0, 6.0]) + @test state_values(r) == [5.0, 6.0] + @test solve(r, alg).u == [2, 1] + end +end + +@testset "saved remake preserves replacement numeric type" begin + linblock() = LinearProblem(Diagonal(ones(2)), [1.0, 2.0]; u0 = [10.0, 20.0]) + @testset "alg=$(nameof(typeof(alg)))" for alg in ( + NewtonRaphson(), SSRootfind(NewtonRaphson()), + DynamicSS(Rodas5P()), SICNM(Rodas5P()), + ) + prob = SCCSteadyStateProblem( + (linblock(),), (Returns(nothing),), nothing, true; + sys = SymbolicIndexingInterface.SymbolCache([:x, :y], [], nothing) + ) + saved = solve(prob, alg; save_idxs = [2]) + unsaved = state_values(prob)[1] + b = big"1.2345678901234567890123456789" + r = remake(saved.prob; u0 = [b]) + @test state_values(r) == [b] + @test eltype(state_values(r)) == BigFloat + @test state_values(r.prob)[1] == unsaved + @test state_values(remake(prob; u0 = [10.0, b]))[2] == b + + d = ForwardDiff.Dual(1.5, 1.0) + r = remake(saved.prob; u0 = [d]) + @test state_values(r) == [d] + @test ForwardDiff.partials(only(state_values(r))) == ForwardDiff.partials(d) + @test state_values(r.prob) == [unsaved, d] + + saved = solve(prob, alg; save_idxs = 2) + r = remake(saved.prob; u0 = b) + @test state_values(r) == b + end +end + @testset "stateless and empty SCC blocks" begin @testset "alg=$(nameof(typeof(alg)))" for alg in ( DynamicSS(Rodas5P()), SICNM(Rodas5P()), From f278d45d1ae25a23528931062fbdc71a34c5747f Mon Sep 17 00:00:00 2001 From: Chris Rackauckas Date: Thu, 15 Oct 2026 11:13:38 +0000 Subject: [PATCH 11/11] Apply solve-time u0 and p overrides on SCC solve entrypoints Solve-time `u0`/`p` keywords on `SCCSteadyStateProblem`, `SCCNonlinearProblem` and the saved-state view were forwarded into the block solves instead of remaking the problem: nonlinear algorithms silently ignored them, and `DynamicSS`/`SICNM` block solves received full-layout `u0`/`p`, leaving the solution's `prob` reporting stale state and parameter values. Route the keywords through `remake` before dispatch so `u0` is chunked into the blocks, `parameters_alias` is enforced, and the returned solution is expressed on the updated problem. Co-Authored-By: Chris Rackauckas Co-Authored-By: Devin Agent-Harness: Devin CLI 3000.11.92 Agent-Model: fusion-claude-opus-5-5-high-sidekick-swe-2-medium Agent-Session: local headless run, log /Users/chrisrackauckas/sandbox/master/green/fu/SciML_SteadyStateDiffEq.jl-178-1100/log.txt --- docs/src/index.md | 7 ++++ src/solve.jl | 18 +++++++++- test/scc.jl | 90 +++++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 114 insertions(+), 1 deletion(-) diff --git a/docs/src/index.md b/docs/src/index.md index 887bc14..0d79063 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -75,6 +75,13 @@ states that were not saved and observed expressions (which need the full state) throw. Parameters and independent variables index as usual, and `sol.original` is the full solution over all states. +Solve-time `u0` and `p` keywords are applied to the problem with `remake` +before dispatch: `u0` is given in the compressed state layout — in the saved +layout when solving a saved-state view — and `p` is subject to +`parameters_alias`, so an unaliased problem throws since its blocks hold +separate parameter objects. The returned solution's `prob` is the updated +problem. + ## API ```@docs diff --git a/src/solve.jl b/src/solve.jl index 50a6a22..d6940eb 100644 --- a/src/solve.jl +++ b/src/solve.jl @@ -26,6 +26,15 @@ function SciMLBase.__solve( end +# Solve-time `u0`/`p` go through `remake`, which chunks `u0` into the blocks and +# enforces `parameters_alias`; the solve then runs on the updated problem. +function __scc_apply_overrides(prob; u0 = nothing, p = nothing, kwargs...) + if u0 !== nothing || p !== nothing + prob = remake(prob; u0 = u0 === nothing ? missing : u0, p = p === nothing ? missing : p) + end + return prob, kwargs +end + # An SCCNonlinearProblem has no top-level `u0`/`kwargs` fields, so it cannot go # through the generic AbstractNonlinearProblem solve preprocessing. Forward it # directly to the wrapped algorithm instead, which dispatches to the SCC solver @@ -34,6 +43,7 @@ function SciMLBase.solve( prob::SciMLBase.SCCNonlinearProblem, alg::SSRootfind, args...; kwargs... ) + prob, kwargs = __scc_apply_overrides(prob; kwargs...) nlsol = solve(prob, alg.alg, args...; kwargs...) return __build_ssrootfind_solution(prob, nlsol) end @@ -131,6 +141,7 @@ function SciMLBase.solve( prob::SciMLBase.SCCNonlinearProblem, alg::Union{DynamicSS, SICNM}, args...; kwargs... ) + prob, kwargs = __scc_apply_overrides(prob; kwargs...) return __solve_scc_blocks(prob, alg, args...; kwargs...) end @@ -138,6 +149,7 @@ function SciMLBase.solve( prob::SciMLBase.SCCSteadyStateProblem, args...; save_idxs = nothing, kwargs... ) + prob, kwargs = __scc_apply_overrides(prob; kwargs...) alg = isempty(args) ? nothing : first(args) alg === nothing && (alg = get(kwargs, :alg, nothing)) kwargs = (; (n => v for (n, v) in pairs(kwargs) if n !== :alg)...) @@ -274,8 +286,12 @@ function SymbolicIndexingInterface.finalize_parameters_hook!(w::SCCSavedProblem, end function SciMLBase.solve( - w::SCCSavedProblem, args...; save_idxs = nothing, kwargs... + w::SCCSavedProblem, args...; save_idxs = nothing, u0 = nothing, p = nothing, + kwargs... ) + if u0 !== nothing || p !== nothing + w = remake(w; u0 = u0 === nothing ? missing : u0, p = p === nothing ? missing : p) + end composed = if save_idxs === nothing w.idxs else diff --git a/test/scc.jl b/test/scc.jl index a4ab9c9..4270405 100644 --- a/test/scc.jl +++ b/test/scc.jl @@ -647,6 +647,96 @@ end end end +@testset "solve-time u0 and p overrides" begin + tol = (; abstol = 1.0e-10, reltol = 1.0e-10) + kblock(p) = NonlinearProblem( + NonlinearFunction( + (u, p) -> [p[1] - u[1]]; + sys = SymbolicIndexingInterface.SymbolCache([:x], [:k]) + ), [0.0], p + ) + function xblock(u0) + return NonlinearProblem((u, p) -> [u[1] - u[1]^3], [u0]) + end + function yblock(u0) + return NonlinearProblem((u, p) -> [4u[1] - u[1]^3], [u0]) + end + cubic() = NonlinearProblem( + (u, p) -> [u[1] - u[1]^3, 4u[2] - u[2]^3], [0.8, 1.5] + ) + function mkparam() + p = [1.0] + return SCCSteadyStateProblem( + (kblock(p),), (Returns(nothing),), p, true; + sys = SymbolicIndexingInterface.SymbolCache([:x], [:k]) + ) + end + mkblocks() = ( + SCCSteadyStateProblem( + (xblock(0.8), yblock(1.5)), (Returns(nothing), Returns(nothing)), + nothing, true; sys = SymbolicIndexingInterface.SymbolCache([:x, :y], []) + ), + SCCSteadyStateProblem( + (cubic(),), (Returns(nothing),), + nothing, true; sys = SymbolicIndexingInterface.SymbolCache([:x, :y], []) + ), + ) + @testset "alg=$(nameof(typeof(alg)))" for alg in ( + NewtonRaphson(), SSRootfind(NewtonRaphson()), + DynamicSS(Rodas5P()), SICNM(Rodas5P()), + ) + prob = mkparam() + sol = solve(prob, alg; p = [3.0], tol...) + @test sol.u ≈ [3] atol = 1.0e-6 + @test sol.ps[:k] == 3 + @test parameter_values(sol.prob) == [3.0] + @test sol.resid ≈ [parameter_values(sol.prob)[1] - sol.u[1]] atol = 1.0e-8 + @test prob.ps[:k] == 1 + + saved = solve(prob, alg; save_idxs = [:x], tol...) + re = solve(saved.prob, alg; p = [3.0], tol...) + @test re.u ≈ [3] atol = 1.0e-6 + @test re.ps[:k] == 3 + + unaliased = SCCSteadyStateProblem( + (kblock([1.0]),), (Returns(nothing),), [1.0], false; + sys = SymbolicIndexingInterface.SymbolCache([:x], [:k]) + ) + @test_throws ArgumentError solve(unaliased, alg; p = [3.0]) + + if !(alg isa typeof(NewtonRaphson())) + p = [1.0] + nlprob = SCCNonlinearProblem( + (kblock(p),), (Returns(nothing),), p, true; + sys = SymbolicIndexingInterface.SymbolCache([:x], [:k]) + ) + nlsol = solve(nlprob, alg; p = [3.0], tol...) + @test nlsol.u ≈ [3] atol = 1.0e-6 + @test parameter_values(nlsol.prob) == [3.0] + end + + for prob in mkblocks() + sol = solve(prob, alg; u0 = [-0.8, -1.5], tol...) + @test sol.u ≈ [-1, -2] atol = 1.0e-6 + @test state_values(prob) == [0.8, 1.5] + + saved = solve(prob, alg; save_idxs = [:y, :x], tol...) + @test saved.u ≈ [2, 1] atol = 1.0e-6 + re = solve(saved.prob, alg; u0 = [-1.5, 0.8], tol...) + @test re.u ≈ [-2, 1] atol = 1.0e-6 + @test re.u ≈ solve(remake(saved.prob; u0 = [-1.5, 0.8]), alg; tol...).u atol = 1.0e-6 + + saved = solve(prob, alg; save_idxs = [:y], tol...) + r = solve(saved.prob, alg; u0 = [-1.5], tol...) + @test r.u ≈ [-2] atol = 1.0e-6 + @test r.original.u ≈ [1, -2] atol = 1.0e-6 + + saved = solve(prob, alg; save_idxs = :y, tol...) + @test solve(saved.prob, alg; u0 = -1.5, tol...).u ≈ -2 atol = 1.0e-6 + end + end +end + @testset "stateless and empty SCC blocks" begin @testset "alg=$(nameof(typeof(alg)))" for alg in ( DynamicSS(Rodas5P()), SICNM(Rodas5P()),