Solve SCCSteadyStateProblem; stop routing SteadyStateProblem through a lowering

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"

docs/src/index.md

@@ -48,6 +48,40 @@ 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`.
+
+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.
+
+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

src/SteadyStateDiffEq.jl

@@ -9,9 +9,13 @@ 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, SymbolicIndexingInterface, finalize_parameters_hook!,
+    parameter_values, set_parameter!, set_state!,
+    state_values, symbolic_type, variable_index, variable_symbols
 
 const infnorm = Base.Fix2(norm, Inf)
 
@@ -20,8 +24,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

src/algorithms.jl

@@ -92,6 +92,11 @@ requires the SCC solver implementation from `SCCNonlinearSolve.jl` to be loaded
 automatically), or an `SCCAlg` configured with separate linear and nonlinear
 block solvers.
 
+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.
+
 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
@@ -148,10 +153,14 @@ 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. 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
 
@@ -251,7 +260,10 @@ 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
+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. 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.