# Convergence and Initialization Lack of convergence is the standard failure mode of a flowsheet simulator. A recycle loop is a fixed point, a fixed point needs a starting guess, and a bad guess turns a correct model into a solver that limit-cycles, overshoots into negative flows, or returns an answer that looks like every other answer and is off by the tear residual. This page is the whole story in one place: where the initial guess comes from, what the three acceleration methods cost, which knobs are tuning and which are correctness, what changes under `jax.grad`, and what to reach for next when the loop still will not close. ## Table of Contents - [When a solve does not converge](#when-a-solve-does-not-converge) - [Where the initial guess comes from](#where-the-initial-guess-comes-from) - [Unit-level `initialize()`](#unit-level-initialize) - [The three acceleration methods](#the-three-acceleration-methods) - [The traced fallback](#the-traced-fallback) - [The equation-oriented route](#the-equation-oriented-route) - [Diagnostics](#diagnostics) - [A worked diagnosis](#a-worked-diagnosis) - [Measured pass rates](#measured-pass-rates) - [Tear selection](#tear-selection) - [What the plugins add](#what-the-plugins-add) --- ## When a solve does not converge `Flowsheet.solve()` warns rather than staying quiet (see [When the Solve Does Not Converge](streams-and-flowsheets.md#when-the-solve-does-not-converge)), so the first thing to do is read the message it already wrote: it names the tear streams, the residual reached, the tolerance, the iteration count and the method that ran. Then work down this list. The snippets on this page share one small loop to run against: a single recycle whose tear is named `tear`, with loop gain 0.97 on a magnitude of about `1e-6` (the corpus case `trace_recycle`). ```python from difflow import Flowsheet, Unit, make_stream from difflow.streams import get_flows def loop(feed, tear): f, t = get_flows(feed), get_flows(tear) return make_stream({"A": f["A"] * 1e-6 + t["A"] * 0.97}, feed["T"], feed["P"]) fs = Flowsheet(["A"], default_flow=0.01) fs.add_feed("feed", make_stream({"A": 1.0}, 300.0, 101325.0)) fs.add_unit(Unit("loop", loop, ["feed", "tear"], ["loop_out"])) fs.add_unit(Unit("out", lambda s: s, ["loop_out"], ["product"])) fs.add_recycle("loop_out", "tear") ``` 1. **Is it actually diverging, or just slow?** Compare `fs.last_solve_residual` against `fs.last_solve_tol`. A residual a factor of two or three above tolerance after `max_iter` iterations is a loop that needs more iterations; a residual of the same order as the flows themselves is a loop that is not converging at all. Raise `max_iter` only for the first. 2. **Give it a starting point.** `tear_initial` is the only guess you control directly, and a guess taken from a converged neighbouring case (a slightly different feed, a slightly different setpoint) is usually worth more than any change of method. See [Explicit guesses](#explicit-guesses) below. 3. **Damp it, or change the method.** The default is Anderson. A loop that *oscillates* — the residual alternating in sign, flows overshooting and coming back — has a tear map with negative eigenvalues, and there are two answers to that. `acceleration="none", damping=0.3` takes a short step toward the substitution value, which makes the iteration contractive where the undamped one is not; it is the one to reach for when the residual *rises*, and the only one that also works under `jax.grad`. `acceleration="wegstein"` gets there by a different route — its acceleration factor is clipped into `[-5, 0]`, which damps toward the previous iterate in exactly that regime — and needs no tuning. A loop that is merely slow and monotone is what Anderson is for. Undamped `acceleration="none"` is the honest baseline: if that diverges, the map itself is not contractive. 4. **Check `clip_negative_flows`.** If any tear flow is legitimately signed, the default clipping is not a safeguard but a bug, and the solve cannot converge. The warning says so when it happened: it counts the iterations on which the projection moved the proposed iterate. See [Clipping negative flows safeguards an extrapolation](#clipping-negative-flows-safeguards-an-extrapolation). 5. **Look for a unit that could not run.** A `TearInitializationWarning` during the solve says the initial pass over the units failed somewhere and that unit's outlets started from the flowsheet default. That is often the real cause: the loop is being started from a stream that is nothing like the answer. Fix it by giving that unit an [`initialize()`](#unit-level-initialize), or by passing `tear_initial` directly. 6. **Scale the default stream.** For a loop carrying trace species, the default tear guess of 0.01 mol/s per species is orders of magnitude too large. `Flowsheet(species_order, default_flow=1e-6)` moves it. 7. **Switch to the equation-oriented solver.** `fs.solve_eo()` assembles every unit's equations into one system and runs Newton on it. Tightly coupled loops that take hundreds of tear iterations often take a handful of Newton steps. By default it initializes *itself* by running the sequential solver first — see [the equation-oriented route](#the-equation-oriented-route). 8. **Make non-convergence loud.** In a script or a test, escalate so a bad solve cannot be mistaken for a good one: ```python streams = fs.solve(on_nonconvergence="raise") ``` --- ## Where the initial guess comes from The tear streams are the **destinations** of the recycles registered with `fs.add_recycle(source, dest)`. A destination is an inlet, so nothing in the flowsheet computes it before the iteration starts, and it has to be guessed. `solve()` tries three sources in order. ### Explicit guesses `tear_initial` takes a dict of `{destination_name: Stream}`. It reaches the solver unmodified and wins over everything below: ```python guess = make_stream({"A": 2.4}, 300.0, 1e5) streams = fs.solve(tear_initial={"tear": guess}) ``` The natural sources of a good guess are a converged solve of a nearby case (continuation by hand), a hand calculation of the loop's material balance, or the result of a looser solve: ```python coarse = fs.solve(tol=1e-4, max_iter=200) fine = fs.solve(tear_initial={"tear": coarse["loop_out"]}, tol=1e-10) ``` ### One pass from the feeds With `use_initialization=True` (the default) and no explicit guess, `solve()` runs every unit once from the feeds with the recycles unknown, and takes the value that pass puts on the *source* end of each recycle. One pass serves every tear, however many recycles there are. Units run in the order they were added, so ordering the flowsheet so that material flows forward through the list makes this pass more useful. A unit whose inlets are not available yet sees the flowsheet default stream. A unit that **raises** on this pass is not necessarily broken: the pass runs before any recycle is known, so a unit inside a loop sees 0.01 mol/s arriving, and a root find or a Rachford-Rice solve can genuinely fail on an almost empty stream. When that happens `solve()` falls back to the unit's [`initialize()`](#unit-level-initialize). When there is no `initialize()` either, the unit's outlets start from the flowsheet default and a `TearInitializationWarning` says so, because a tear guess that is quietly worse than it looks is the thing worth avoiding: ```python from difflow.initialization import TearInitializationWarning ``` (It lives in `difflow.initialization`, not on the top-level package.) ### The flowsheet default stream The last resort, and what is used when `use_initialization=False`: every species at `default_flow`, at `default_T` and `default_P`. Those three are constructor arguments: ```python fs_small = Flowsheet(["A", "B"], default_flow=1e-6, default_T=350.0, default_P=2e5) ``` Under tracing the guessing pass is skipped entirely and the default stream is used, because a best-effort pass that catches exceptions has no business running over tracers — and the initial guess cannot affect a gradient that comes from the converged solution anyway. --- ## Unit-level `initialize()` `difflow.initialization.Initializable` is a runtime-checkable protocol: a unit supports initialization if it has an `initialize(inlet, **kwargs)` method. The flowsheet's only test is `isinstance(operation, Initializable)`, so nothing has to be registered. Three core units implement it, each with an analytic estimate that cannot fail the way the real calculation can: | Unit | Estimate | |---|---| | `CSTR` | residence time from the volumetric flow, then `X = k tau / (1 + k tau)` | | `PFR` | the analytic plug-flow conversion for the same rate constant | | `Flash` | a simplified Rachford-Rice from the thermo package's K-values, with a bubble/dew check for the fully liquid and fully vapor cases | It is consulted only where the real call would have failed: the flowsheet calls `__call__` first on its initialization pass and reaches for `initialize()` when that raises. An initializer is a fallback, not a pre-processing step. The return value is a dict. The flowsheet reads only the outlet guesses out of it, under either of two keys: - `outlets`: the guesses as a **sequence, in the order `__call__` returns them**. That order is the only thing a caller holding `Unit.outlet_names` can match the guesses up by, so a multi-outlet unit has to provide it to be usable by the flowsheet. - `outlet`: the singular spelling, accepted for a unit with one outlet, and meaning the same as a one-element `outlets`. `states` and `info` may carry anything else; the flowsheet ignores them, and they are for a caller that knows which unit it is talking to (`Flash` also returns its guesses under `liquid` and `vapor`). Writing one for a custom unit: ```python from difflow.streams import make_stream class MyReactor: def __call__(self, inlet): ... # the real calculation, which may fail on an empty inlet def initialize(self, inlet, **kwargs): outlet = make_stream({"A": 0.5 * inlet["F_A"]}, inlet["T"], inlet["P"]) return {"outlets": (outlet,), "outlet": outlet} ``` `**kwargs` receives the unit's `params` dict, so an initializer may read the same specifications the call does. Keep it cheap and unconditional: its whole value is that it answers where the real calculation would not. --- ## The three acceleration methods All three iterate on the same packed tear vector, which is, per tear stream and in `species_order`, `[F_1, ..., F_N, T, P]`, with the streams in sorted name order. ### Direct substitution (`"none"`) `x_{k+1} = x_k + \alpha (g(x_k) - x_k)`, run through `optimistix.FixedPointIteration`, with $\alpha$ = `damping` = 1 by default. No history, no extrapolation, and the only path that is traceable — which is why it is also what `jax.grad` falls back to. Undamped it converges linearly at a rate set by the spectral radius of the tear map and diverges when that exceeds one; [damping](#damping) is what extends it to maps that overshoot. Its convergence verdict is judged against the criterion `optimistix` actually stops on — elementwise `|dx| < atol + rtol |x|`, with `tol` serving as both — not against a bare `|dx| < tol`. The two differ by the magnitude of the tear: with flows of order 1 and `tol=1e-8`, a solve `optimistix` calls successful lands near `2e-8`. `last_solve_residual` alongside it stays the plain max-norm, which is the number to compare against `tol` by eye. ### Wegstein Estimates the tear map's slope from the last two iterates and extrapolates: $$x^{(k+1)} = q\, x^{(k)} + (1 - q)\, g(x^{(k)}), \qquad q = \frac{s}{s - 1}, \qquad s = \frac{g(x^{(k)}) - g(x^{(k-1)})}{x^{(k)} - x^{(k-1)}}$$ $q$ is the weight on the **old** iterate and $(1-q)$ the weight on the substitution value. For a linear $g$ this cancels the map's slope exactly and lands on the solution in one step. $q$ is clipped into $[-5, 0]$. That bound is what makes Wegstein the right choice for an oscillating loop: for $s < 0$ — the overshoot case — $q$ lands in $(0, 1)$ and the update is damped toward the previous iterate, while for $0 < s < 1$ it is negative and the update extrapolates past $g$. Cost is one flowsheet pass per iteration plus two vectors of storage. ### Anderson The default. Keeps a rolling history (`anderson_depth`, default 5) and solves a small least-squares problem for the combination of past iterates that minimises the residual — equivalent to GMRES on the fixed-point iteration. It is exact on an affine map within `depth` steps and is the best general-purpose choice for a loop that is slow but not oscillating. Two implementation details matter to a user: - The least-squares problem is solved with `jnp.linalg.lstsq` on the difference matrix directly, **not** through the normal equations. The normal equations square the condition number, which overflows to NaN for a badly scaled packed tear vector — gas transmission networks, whose `[F, T, P]` entries span 1e1 kg/s to 1e6 Pa, have a normal-equation condition number of order `(P/F)^4`. Badly scaled tear vectors are therefore fine. - If an accelerated step comes back non-finite (a degenerate history), the iteration silently falls back to plain substitution for that step rather than propagating NaN. Cost is one flowsheet pass per iteration plus `depth` vectors of storage and a small `lstsq`. On a representative nonlinear one-loop test problem solved to `1e-10`: | `acceleration` | iterations | |---|---| | `"none"` | 19 | | `"wegstein"` | 5 | | `"anderson"` | 8 | Those numbers are problem-specific — Wegstein's scalar slope estimate is at its best on a loop with one dominant mode, and Anderson pulls ahead as the tear vector grows — but the order of magnitude is typical: acceleration is worth a factor of a few, not a factor of a hundred. ### Clipping negative flows safeguards an extrapolation Both accelerated methods extrapolate, and an extrapolated molar flow can land below zero on the way to a solution that is entirely positive. By default the accelerated iterations therefore project the **flow** entries of the tear vector onto `[0, inf)` after each step. Temperature and pressure are never touched. This is right for a chemical flowsheet and wrong wherever a negative tear flow is a legitimate answer. Gas transmission networks are the standing example: flow against an arc's reference direction is routine, the fixed point genuinely has negative components, and clipping means the iteration can never reach it. `difflow_gas.GasNetworkFlowsheet` overrides the default to `False` for exactly this reason. ```python streams = fs.solve(clip_negative_flows=False) # signed tear flows ``` If a solve with signed flows will not converge and nothing else explains it, this is the first thing to check. ### Clipping binds on the accelerated paths only `acceleration="none"` does not clip, and passing the flag to it changes nothing. That is a decision (#263), not an omission, and it follows from what the projection is *for*. | `acceleration` | proposes | clip applies | |---|---|---| | `"none"` | `x + α(g(x) − x)`, a step toward a value the units computed | no | | `"wegstein"` | an extrapolation from the last two iterates | yes | | `"anderson"` | an extrapolation from the last `depth` iterates | yes | The projection guards a *guess no unit computed*. Wegstein and Anderson each propose one, and it can leave the physical orthant on the way to a solution entirely inside it. The unaccelerated step is a convex combination of `x` and `g(x)` (for the usual `α ≤ 1`), so it cannot leave a region both ends are in; there is no guess to guard. Clipping there would clip `g` itself, which is a different act with two costs: - **It invents fixed points.** `clip(g(x)) = x` has solutions where `g(x) = x` has none — the module docstring of `difflow.convergence` names this as one of the ways a solve reports success and lands somewhere that does not audit. - **It breaks gradients.** A traced solve falls back to `"none"`, so that path *is* the one `jax.grad` runs on, and the clip would sit inside the map `optx.fixed_point` implicitly differentiates. Wherever a tear flow converges to exactly zero — an ordinary species absent from a recycle — `jnp.maximum(x, 0)` contributes a derivative of one half at its kink. On the two-species loop in `tests/test_flowsheet_tear_clip.py` the forward answer stays right and the gradient comes back 1.333 against an exact 2.0. What the asymmetry used to cost was silence: an accelerated solve on a signed tear failed with nothing in the output naming the projection as the reason. That is now reported. `fs.last_solve_clip_active` counts the iterations on which the clip actually moved the proposed iterate, and a non-convergence warning with a non-zero count says so and names the remedy: ```text Recycle solve did not converge: tear stream(s) tear reached a residual of 1.324e+00 after 100 of 100 iterations ... The negative-flow clip moved the proposed iterate on 100 of those iterations: if any tear flow here is legitimately signed (a gas network, a power flow, any tear carrying a signed quantity), that projection puts the fixed point out of reach -- pass clip_negative_flows=False. ``` ### Measuring the error, not the step `tol` tests the max-norm of $g(x) - x$, the **step** between successive tear iterates. The quantity a caller cares about is the **error**, $x^* - x$, and those are not the same number: $$x^* - x = (I - M)^{-1}\left(g(x) - x\right), \qquad M = \frac{\partial g}{\partial x}$$ For a loop with one dominant mode of gain $g$ that is a factor of $1/(1-g)$: 33 at a gain of 0.97, 1000 at 0.999. A solve that stops inside `1e-8` on such a loop reports convergence, means it, and is decades outside the tolerance it was given. After a solve converges, `solve` spends `error_probe` extra plain-substitution passes measuring that factor rather than assuming it (#264): $$r_j = g(x_j) - x_j, \quad x_{j+1} = g(x_j), \qquad x^* - x_0 = \sum_j r_j \approx \sum_{j 1$; a map with $s = 3$ has its fixed point half a residual away, not infinitely far. - **Entries already at round-off are dropped.** They say nothing about the error, and a ratio taken from floating-point noise lands near one — which is exactly where $1/(1-s)$ blows up. A tear carrying a pressure of $10^5$ alongside flows of order one has such an entry on almost every solve. #### The warning A converged solve whose measured error is past `tol` raises `TearToleranceWarning` — but only when the error is at least ten times the step, which is a gain of 0.9 or more. Below that the gap is the ordinary slack of a step test (the unaccelerated path deliberately stops on $|g(x)-x| < \mathrm{tol}\,(1+|x|)$), and a warning that fires on every solve is read as noise exactly when it is not. Across the corpus — eleven cases times three accelerations — it fires once, on the one solve that is genuinely wrong. #### Asking for the error instead `tol_basis="error"` makes `tol` mean what it reads as. When the probe says the answer is outside it, `solve` tightens its step test by the factor just measured and solves again, warm-started, up to three times: ```python fs.solve(tol=1e-8, tol_basis="error") # on trace_recycle, gain 0.97 # 53 iterations -> 71, error estimate 3.0e-07 -> 9.7e-09 ``` It is opt-in because it costs iterations: a loop that converges today needs more of them to satisfy it, and changing that for every existing caller is not a cost to impose silently. It is also still an **absolute** test. On a tear whose converged value is $3\times10^{-5}$, an error of `1e-8` is a relative error of $3\times10^{-4}$ — the scale half of the trap, which is the caller's to set with `tol` (or with `Flowsheet(species_order, default_flow=...)` for a trace loop). The gain half is what is measured here. ### Damping `damping` takes a step only part of the way toward the substitution value: $$x^{(k+1)} = x^{(k)} + \alpha\left(g(x^{(k)}) - x^{(k)}\right)$$ The fixed point is unchanged — at $g(x) = x$ the correction vanishes for any $\alpha$ — but the iteration is contractive where plain substitution overshoots. For a tear map whose eigenvalue spectrum lies in $[-m, 0)$ it contracts for $\alpha < 2/(1+m)$, with the optimal scalar damping near $\alpha = 2/(2+m)$. An $\alpha$ of 0.3 is safe up to a spectral radius of about 5, which is why it is the usual first thing to try; when in doubt, measure, since a finite-difference Jacobian of the tear map at a solved state gives $m$ directly. The analysis is general even though it was written down for gas networks, where measured spectra reach $m \approx 3.1$ (GasLib-40). It applies to `acceleration="none"` — and therefore to the traced path, which falls back to it. The default is `1.0`, which is plain undamped substitution: ```python streams = fs.solve(acceleration="none", damping=0.3) ``` On a loop whose tear map has eigenvalue $-2$ — substitution triples the error and flips its sign every step — the difference is not a matter of iteration count: | `damping` | outcome after 200 iterations | |---|---| | `1.0` | residual `4.8e+60`; diverged | | `0.3` | residual `3.0e-12`; converged to the exact fixed point | Three conventions keep this a knob on *how the solve gets there* rather than on what it returns: - **`damping` never moves the answer.** The fixed point is a property of the map. Converged solves at 1.0, 0.5, 0.3 and 0.1 agree to the tolerance. - **`last_solve_residual` is the undamped residual** $|g(x) - x|$, not the $\alpha$-times-smaller step the solver actually took. A heavily damped solve would otherwise look converged by exactly the factor it was damped by. - **`tol` still means `tol`.** `optimistix` stops on the step it takes, so the tolerance handed to it carries the damping factor. Without that, a solve at $\alpha = 0.3$ would stop three times short of the tolerance it was asked for — and the verdict above would then, correctly, call it non-converged. Gradients are unaffected: `optimistix` differentiates the converged solution through the implicit function theorem, so the derivative is exact whatever $\alpha$ and however many iterations ran. That cuts both ways — on the diverging run above, the *gradient* is still 1.0 while the *value* is `-6.4e+60`. A gradient is only as trustworthy as the solve underneath it, and under tracing nothing warns you. `damping` must be positive; `0.0` makes the iteration the identity, under which every point is a fixed point, and is refused. ```{note} `difflow_gas.GasNetworkFlowsheet.solve_differentiable(alpha=...)` predates this and does the same thing with separate `rtol`/`atol` and an optional stats return. It is still the convenient entry point for gas networks, which also need `clip_negative_flows=False` and the builder's tear guesses. ``` --- ## The traced fallback The accelerated solvers are Python loops that branch on the residual to stop early, and a tracer cannot answer `residual < tol`. So when `solve()` sees a tracer in any feed or unit parameter, it **replaces the requested acceleration with the `optimistix` fixed point** and records that it did: ```python import jax def objective(feed_F): fs.feeds["feed"] = make_stream({"A": feed_F}, 300.0, 1e5) return fs.solve(tol=1e-10, max_iter=200)["loop_out"]["F_A"] gradient = jax.grad(objective)(1.0) fs.last_solve_method # 'fixed_point (traced)' fs.last_solve_converged # None -- the residual was a tracer fs.feeds["feed"] = make_stream({"A": 1.0}, 300.0, 101325.0) # undo the traced feed ``` Three consequences worth knowing before differentiating through a flowsheet: - **The gradient is exact, and is not the gradient of the iteration.** `optimistix` carries an implicit-differentiation rule, so the derivative comes from the converged solution through the implicit function theorem, not from an unrolled loop. It does not depend on how many iterations ran. - **You lose the acceleration.** A loop that needs Anderson to converge in 100 iterations eagerly may need many more under `grad`. Raise `max_iter` on the differentiated path. - **You lose the convergence verdict.** There is no concrete residual to judge, so `last_solve_converged` is `None` and the solve stays silent whatever `on_nonconvergence` says — guessing would either raise inside a gradient or warn on a solve that was fine. Check convergence eagerly, at the same parameter values, before trusting a gradient. - **You lose the feed-propagation guess.** The initialization pass catches exceptions, which has no business running over tracers, so a traced solve starts from the flowsheet default stream unless you pass `tear_initial`. For a loop the default start is nowhere near, pass the eagerly converged tear in explicitly. The same rule explains a related gotcha: never record a solve diagnostic with a bare `float()` in code that may be traced. `flowsheet._concrete()` returns `None` under tracing instead of raising. --- ## The equation-oriented route `fs.solve_eo()` assembles every unit's `eo_residuals` into one system $F(x) = 0$ and solves it with Newton via `optimistix` (see [Equation-Oriented (EO) Solver](eo-solver.md)). It converges quadratically near the solution, which is why a tightly coupled loop that takes hundreds of tear iterations often takes a handful of Newton steps — and why it needs a starting point inside the basin. By default it gets one by **running the sequential solver first**: ```python streams = fs.solve_eo(use_sm_init=True) # the default ``` `use_sm_init` runs `fs.solve(tol=1e-6, max_iter=50, acceleration="anderson")` and uses the result — converged or not — as the Newton start, keeping whichever streams the EO state vector needs. A loose sequential solve is cheap and lands close enough for Newton far more often than a feed-propagation guess does. If the sequential solve raises, `solve_eo` falls back to feed propagation rather than failing. Two things follow: - A non-converged sequential solve still makes a fine EO start. That is the point: SM is being used for its robustness, Newton for its accuracy. - Because the sequential pre-solve swallows exceptions, a `solve_eo` that behaves as if it had a bad start may have had one. Run `fs.solve()` on its own to find out. `use_sm_init=False` uses feed propagation instead, and an explicit `initial_guess` (a dict of streams, like `tear_initial` but covering every unknown stream) overrides both. --- ## Diagnostics Every solve records what happened, whichever path ran and whether or not it converged: | Attribute | Meaning | |---|---| | `last_solve_converged` | `True` met the tolerance, `False` ran out of iterations, `None` the residual was a tracer and there is nothing to judge | | `last_solve_residual` | final tear residual, as a plain max-norm of `g(x) - x` | | `last_solve_tol` | the tolerance that was asked for | | `last_solve_iterations` | iterations used; equals `max_iter` when it did not converge | | `last_solve_method` | `"anderson"`, `"wegstein"`, `"none"`, `"fixed_point (traced)"`, or `"direct"` for a recycle-free sequential solve | | `last_solve_tear_streams` | the tear (recycle destination) names | | `last_solve_clip_active` | iterations on which `clip_negative_flows` actually moved the proposed iterate; always `0` on the unaccelerated path, which does not clip | | `last_solve_gain` | measured signed gain of the tear map at the solution, or `None` where there was nothing to measure | | `last_solve_error_estimate` | measured error left in the tear — the number `tol` reads as if it were | The tri-state `last_solve_converged` is the one to test against `is False` rather than falsy: `None` means "not judged", not "failed". `ConvergenceWarning` and `ConvergenceError` are exported from `difflow`, so ```python import warnings from difflow import ConvergenceWarning warnings.simplefilter("error", ConvergenceWarning) ``` turns every non-converged solve in a script into a failure. --- ## A worked diagnosis A one-unit loop whose tear map is $F \mapsto F_{\text{feed}} + k\sqrt{F + 0.1}$, deliberately given too few iterations: ```python import warnings import jax import jax.numpy as jnp from difflow import ConvergenceWarning from difflow.flowsheet import Flowsheet, Unit from difflow.streams import make_stream jax.config.update("jax_enable_x64", True) class Loop: def __init__(self, k): self.k = k def __call__(self, feed, tear): return make_stream( {"A": feed["F_A"] + self.k * jnp.sqrt(tear["F_A"] + 0.1)}, feed["T"], feed["P"], ) def build(): fs = Flowsheet(species_order=["A"], default_flow=0.0) fs.add_feed("feed", make_stream({"A": 1.0}, 300.0, 1e5)) fs.add_unit(Unit("loop", Loop(0.9), ["feed", "tear"], ["loop_out"])) fs.add_recycle("loop_out", "tear") return fs fs = build() with warnings.catch_warnings(): warnings.simplefilter("ignore", ConvergenceWarning) streams = fs.solve(max_iter=3) assert fs.last_solve_converged is False assert fs.last_solve_method == "anderson" assert fs.last_solve_residual > fs.last_solve_tol # 1.9e-2 against 1e-8 ``` The residual is six orders above tolerance, so this is not a loop that needs a few more iterations. A starting point taken from anywhere near the answer buys four orders of magnitude at the *same* three iterations: ```python import warnings import jax import jax.numpy as jnp from difflow import ConvergenceWarning from difflow.flowsheet import Flowsheet, Unit from difflow.streams import make_stream jax.config.update("jax_enable_x64", True) class Loop: def __init__(self, k): self.k = k def __call__(self, feed, tear): return make_stream( {"A": feed["F_A"] + self.k * jnp.sqrt(tear["F_A"] + 0.1)}, feed["T"], feed["P"], ) fs = Flowsheet(species_order=["A"], default_flow=0.0) fs.add_feed("feed", make_stream({"A": 1.0}, 300.0, 1e5)) fs.add_unit(Unit("loop", Loop(0.9), ["feed", "tear"], ["loop_out"])) fs.add_recycle("loop_out", "tear") guess = {"tear": make_stream({"A": 2.4}, 300.0, 1e5)} with warnings.catch_warnings(): warnings.simplefilter("ignore", ConvergenceWarning) fs.solve(tear_initial=guess, max_iter=3) assert fs.last_solve_residual < 1e-4 # 8.3e-6, from 1.9e-2 ``` And the gradient through the converged loop comes from the traced path, matching a central difference on the eager one: ```python import jax import jax.numpy as jnp from difflow.flowsheet import Flowsheet, Unit from difflow.streams import make_stream jax.config.update("jax_enable_x64", True) class Loop: def __init__(self, k): self.k = k def __call__(self, feed, tear): return make_stream( {"A": feed["F_A"] + self.k * jnp.sqrt(tear["F_A"] + 0.1)}, feed["T"], feed["P"], ) fs = Flowsheet(species_order=["A"], default_flow=0.0) fs.add_feed("feed", make_stream({"A": 1.0}, 300.0, 1e5)) fs.add_unit(Unit("loop", Loop(0.9), ["feed", "tear"], ["loop_out"])) fs.add_recycle("loop_out", "tear") def objective(feed_F): fs.feeds["feed"] = make_stream({"A": feed_F}, 300.0, 1e5) return fs.solve(tol=1e-10, max_iter=200)["loop_out"]["F_A"] gradient = jax.grad(objective)(1.0) assert fs.last_solve_method == "fixed_point (traced)" assert fs.last_solve_converged is None h = 1e-6 finite_difference = (objective(1.0 + h) - objective(1.0 - h)) / (2 * h) assert abs(gradient - finite_difference) < 1e-6 ``` --- ## Measured pass rates How often any of this actually works is a measurement, and `difflow.convergence` is where it is made. It holds a corpus of deliberately hard flowsheets and runs each one across every acceleration and every initialization strategy: ```python from difflow.convergence import run_benchmark report = run_benchmark() print(report.as_text()) ``` The full grid is 99 solves and takes several minutes, most of it JAX compilation. Narrow it while iterating: ```python from difflow.convergence import run_benchmark report = run_benchmark(cases=["high_gain_recycle"], accelerations=("anderson",)) ``` ### Passing means converged *and* right A solve reports two separate things, and the benchmark keeps them apart: | | meaning | |---|---| | `Outcome.converged` | the solver's own verdict, `last_solve_converged` | | `Outcome.correct` | the answer closes its material balance, or matches a known one | | `Outcome.passed` | both | They are not the same. `pass_rate` is 46.5% against a `convergence_rate` of 49.5%: three of the 99 solves report success and land somewhere that does not audit. A benchmark that counted only the solver's flag would call those wins. ### What the corpus says today ``` 99 solves (tol=1e-08, max_iter=100) pass rate 46.5% (converged AND the answer audits) convergence rate 49.5% (the solver's own verdict) by acceleration by initialization anderson 28/33 84.8% default 15/33 45.5% wegstein 9/33 27.3% unit 15/33 45.5% none 9/33 27.3% feed 16/33 48.5% mean iterations over passing solves: anderson=4.3, wegstein=18.9, none=41.0 ``` **Every case in the corpus is solved by at least one acceleration.** Eleven cases across the hard families — high loop gain, sharp splits, multi-loop, phase-regime switching, near-pinch columns, trace-magnitude tears, signed tears — and none of them defeats all three methods. That is the most useful thing the benchmark says, and it is a negative result. Four findings worth acting on: **Acceleration dominates; initialization barely registers.** Anderson passes 84.8% where plain substitution passes 27.3%. Against that, the three starting points are separated by a single solve. The first thing to change on a failing recycle is `acceleration`, not the guess. **The `initialize()` path (#247) changed no outcome.** `"unit"` — the propagation pass that #247 wired up — scores exactly what `"default"` scores, the 0.01 mol/s stream it replaced: 15/33 either way, unchanged from when the corpus was six cases. It reliably saves an iteration or two and it flips no verdict anywhere. It made the first step better without making a failing solve succeed. If a recycle is failing, a better guess of the same kind is not the fix. **Anderson is not free, and `clip_negative_flows` is why, once.** `phase_coupled_flash` fails under Anderson and converges under Wegstein in 28 iterations — Anderson walks off somewhere else rather than running out of road. And `signed_tear` inverts the ranking outright: plain substitution solves it and *both* accelerated methods fail. The reason is worth knowing — `clip_negative_flows` defaults to `True` and binds on the Wegstein and Anderson paths but not on the unaccelerated one ([why](#clipping-binds-on-the-accelerated-paths-only)), so on a tear whose answer is genuinely negative the two accelerated methods carry a projection the plain one does not. Passing `clip_negative_flows=False` fixes Anderson on it immediately (100 iterations without convergence → 13 with). This is the same point `difflow_gas` already makes for signed flows; it applies to any tear whose components can go negative. Since #263 the failing solve says this for itself rather than leaving it to be worked out. **A disjunction does not need binaries here.** `regime_switch` is a unit choosing between two linear branches with the answer exactly on the boundary, expanding below the switch and contracting above it. Plain substitution can only chatter across the join, and Wegstein stalls — but Anderson lands on it in two iterations. ### What is *not* hard Worth recording, because these are the families most often assumed difficult: - **A rigorous MESH column in a recycle.** `column_recycle` — twelve stages, 95% of the distillate returned — passes under every setting, at every reflux ratio tried from just above minimum to well above it. - **Nonsmooth tear maps as such.** Both `regime_switch` and the kinked maps tried while building the corpus fall to Anderson quickly. - **Tear dimension.** A tear with twelve components and a spectral radius of 0.985 converges under Anderson as readily as one with three, provided the coupling is non-negative. Where high-dimensional cases did fail, the cause was the negative-flow clipping above, not the dimension. ### The trap the corpus exposes `tol` is a residual on the **step**, and on a high-gain loop that is much weaker than it looks. A step of $d$ on a loop of gain $g$ leaves an error of about $d/(1-g)$, so at $g = 0.97$ the residual understates the remaining error by a factor of 33. `trace_recycle` is built on exactly that: gain 0.97 on a tear whose converged value is about $3 \times 10^{-5}$. Wegstein stops with a residual inside `1e-8`, reports convergence, and is about 1% out. It is right on its own terms — the step really is that small — and wrong on yours. Since #264 the solve measures that instead of leaving you to infer it. See [Measuring the error, not the step](#measuring-the-error-not-the-step): the same solve now reports `last_solve_gain = 0.97` and `last_solve_error_estimate = 3.0e-7`, warns that it is outside the tolerance it was given, and `tol_basis="error"` reaches it. **What was *not* a second instance of this.** `signed_tear` was read as the non-normal version of the same trap: $\rho = 0.75$ suggests an amplification of 4, and the case appeared to leave a relative error of $1.2 \times 10^{-6}$ from a residual of $10^{-8}$ — a factor of about 450, blamed on $\lVert (I-M)^{-1} \rVert$ being the real amplification rather than $1/(1-\rho)$. It was not. The case's reference answer was written out to six decimals, which is itself $4.4 \times 10^{-6}$ from $(I-M)^{-1}\mathbf{1}$, and the audit was measuring that rounding. Plain substitution lands $7.6 \times 10^{-9}$ from the answer — *inside* its own tear residual, an amplification of 0.57. The reference is now at full precision and the case audits at `1e-6` like every other analytic one. The general point stands: for a multicomponent tear the worst-case amplification is $\lVert (I-M)^{-1} \rVert$ rather than $1/(1-\rho)$. The corpus just does not happen to contain an example of it, and the measured estimate reports the amplification along the residual actually present rather than a worst case over directions that are not. ### Adding a case A case is a name, a builder, and a sentence saying what makes it hard: ```python from difflow.convergence import Case, run_case from difflow.flowsheet import Flowsheet, Unit from difflow.streams import make_stream def build_my_flowsheet(): """A FRESH flowsheet every call: solve() records its verdict on the object.""" def mix(feed, tear): return make_stream({"A": feed["F_A"] + tear["F_A"]}, feed["T"], feed["P"]) def split(inlet): # 95% back round: loop gain 0.95 flow = inlet["F_A"] return (make_stream({"A": 0.95 * flow}, inlet["T"], inlet["P"]), make_stream({"A": 0.05 * flow}, inlet["T"], inlet["P"])) fs = Flowsheet(["A"], default_flow=0.01) fs.add_feed("feed", make_stream({"A": 1.0}, 300.0, 101325.0)) fs.add_unit(Unit("mix", mix, ["feed", "tear"], ["mixed"])) fs.add_unit(Unit("split", split, ["mixed"], ["recycle", "product"])) fs.add_recycle("recycle", "tear") return fs case = Case( name="my_hard_loop", build=build_my_flowsheet, # must return a FRESH Flowsheet difficulty="Why this one is hard, in a sentence.", tags=("recycle", "high-gain"), ) run_case(case, acceleration="wegstein") ``` `build` has to return a new flowsheet each call — `solve` records its verdict on the object, so a shared one carries one run's result into the next. The default audit is the overall mole balance; a case whose reaction changes the mole count, or whose answer is known analytically, passes its own `check`. Keep at least one case in the corpus that is *meant* to pass. A benchmark where everything fails cannot distinguish a hard corpus from a broken solver. `two_phase_flash` and `column_recycle` are the controls, and both pass 9/9. And check a new case is not *rigged* — that no initialization strategy starts on the answer it is about to be scored against. `overshoot_loop` originally used $g(x) = 3F - 2x$, whose fixed point is $x = F$, exactly the `"feed"` guess; that strategy collected two passes on a case it had not solved. `TestNoCaseIsRigged` checks every case in the corpus against this. --- ## Tear selection Where you tear a loop is a modelling decision with real consequences: the tear is the fixed point the solve iterates on, so its placement sets the spectral radius of the map, and with it how fast the loop converges and sometimes whether it converges at all. By default difflow does not make that decision for you. The tears are whatever `add_recycle` registered, one per recycle, and a loop closed in the topology but never declared is not torn at all — inferring one silently would change the answer for flowsheets that already run. `fs.tear_analysis()` reports the loops, the declared tears, what each strategy would pick instead, and any inlet the unit order cannot supply. It reads the topology and runs no unit: ```python print(fs.tear_analysis()) ``` `solve(tears="auto")` opts in to the choice, but only when no recycle has been declared; `"minimum"` names the other strategy. The selection is recorded in `last_solve_tear_streams`, and neither it nor the unit sequence is kept, so `fs.recycles` and `fs.units` are what you declared once the solve returns. [Streams and Flowsheets](streams-and-flowsheets.md) covers both in full. The pieces are also usable on their own: ```python from difflow import ( FlowsheetGraph, find_cycles, select_tear_streams, calculation_order, ) graph = FlowsheetGraph.from_flowsheet(fs) cycles = find_cycles(graph) # elementary cycles tears = select_tear_streams(fs, method="heuristic") # or "minimum" calculation_order(graph, tears) # order once those are seeded ``` `"heuristic"` takes one tear per loop, preferring a stream leaving a mixing point, where the stream is the sum of everything entering it so a guess wrong in composition is still right in order of magnitude; `"minimum"` covers every loop with as few tears as it can. The general rule is the one the gas plugin states quantitatively: tear where the map is least sensitive. A chord's tear-map slope there goes as `-sum(beta_e |q_e|) / (beta_c |q_c|)` over its loop, so tearing the *most* resistive element keeps the spectral radius small. --- ## What the plugins add Two plugins solve initialization problems the core does not, and both are worth reading even if you work in neither domain. **Gas networks compute the decomposition from the topology.** `difflow_gas.decompose(net, root=...)` builds a spanning tree with the non-invertible arc kinds forced in-tree and the most resistive pipe of each loop pushed out as the chord, then emits a leaf-to-root balance schedule. The tears are derived rather than declared, and the chord choice is made on the convergence criterion above. See [The computed decomposition](unit-operations-gas.md#the-computed-decomposition) and [Building and solving](unit-operations-gas.md#building-and-solving). **REE mass-action sections initialize by continuation from a correlation.** Mass-action equilibrium is exponentially nonlinear and loses Newton from a poor start, so `difflow_ree` solves the existing L1 correlation first and uses its Kremser profile as the starting point, then walks it into Newton's basin with damped Newton, a trust region and damped Newton again — `n_continuation_steps > 1` additionally ramps the rare-earth feed from dilute to full strength. All of it runs under `stop_gradient`: it moves the starting point and never the answer, which is the pattern to copy for any globalization in a differentiable model. See [Unknowns, equations and how they are solved](unit-operations-ree.md#unknowns-equations-and-how-they-are-solved). --- ## See also - [Streams and Flowsheets](streams-and-flowsheets.md) — building flowsheets, recycles, implicit differentiation - [Equation-Oriented (EO) Solver](eo-solver.md) — the simultaneous alternative - [Solvers and Utilities](solvers-and-utilities.md) — the numerics underneath - [Gas Transmission Network Unit Operations](unit-operations-gas.md) — computed decomposition, damped differentiable tear solve