Solvers and diagnostics

Shared numerical reports, checked roots, and implicit differentiation support thermodynamic calculations and process solves. See the reliability guide for convergence checks, scaling, failure handling, and derivative limits.

Solve reports

diagnostics

Numerical solve reports shared by thermodynamics and process simulation.

Reports contain only arrays and are JAX pytrees. Host applications can raise a descriptive exception with :func:require_converged; compiled applications inspect report.converged and retain the report alongside their results.

Classes:

Name	Description
SolveStatus	Stable, machine-readable termination codes.
SolveReport	Termination information for a nonlinear or linear calculation.
SolveResult	A best available solution and its independently checked solve report.
ConvergenceError	A failed calculation, with its machine-readable report attached.

Functions:

Name	Description
nan_unless_converged	Replace every leaf of value with NaN when report did not converge.
with_status	Override a report's status with status (a code or a status array) where failed.
canonical_report	Give every report field a fixed dtype.
require_converged	Raise for a concrete failed report; leave traced reports to compiled callers.
residual_report	Grade a scaled residual independently of an iteration's stopping rule.

SolveStatus

Bases: IntEnum

Stable, machine-readable termination codes.

OUT_OF_DOMAIN marks a well-posed request outside a model's domain (for example, a saturation pressure above the critical temperature). TRIVIAL marks an equilibrium iteration that collapsed onto identical phases where a distinct phase was required.

SolveReport

Bases: NamedTuple

Termination information for a nonlinear or linear calculation.

Attributes:

Name	Type	Description
status	Array	A :class:SolveStatus code.
iterations	Array	Number of attempted iterations.
residual_norm	Array	Maximum absolute scaled residual.
step_norm	Array	Maximum absolute scaled last step.
worst_equation	Array	Index of the largest scaled residual, or -1 if absent.

Methods:

Name	Description
to_dict	Return a strict-JSON-compatible report for a concrete calculation.

converged property

converged: Array

Whether the residual passed its tolerance and is finite.

to_dict

to_dict() -> dict[str, Any]

Return a strict-JSON-compatible report for a concrete calculation.

SolveResult

Bases: NamedTuple

A best available solution and its independently checked solve report.

ConvergenceError

ConvergenceError(
    report: SolveReport,
    context: str = "calculation",
    labels: tuple[str, ...] = (),
)

Bases: RuntimeError

A failed calculation, with its machine-readable report attached.

nan_unless_converged

nan_unless_converged(
    value: Any, report: SolveReport
) -> Any

Replace every leaf of value with NaN when report did not converge.

This is the value-only contract of Fugacio's public equilibrium calls: a failed solve never returns a finite number. The selection is on the primal, so a zero cotangent through a discarded failure stays zero in reverse mode.

with_status

with_status(
    report: SolveReport,
    failed: Array,
    status: SolveStatus | Array,
) -> SolveReport

Override a report's status with status (a code or a status array) where failed.

canonical_report

canonical_report(report: SolveReport) -> SolveReport

Give every report field a fixed dtype.

Reports assembled on different code paths can differ in integer width or JAX weak typing; jax.lax.cond and jax.lax.switch require identical branch outputs, so branches return canonical reports.

require_converged

require_converged(
    report: SolveReport,
    context: str = "calculation",
    labels: tuple[str, ...] = (),
) -> None

Raise for a concrete failed report; leave traced reports to compiled callers.

This function performs no callbacks or side effects inside JAX transforms. Differentiable solver values use nonfinite derivatives on failed solves; compiled callers should also inspect the returned report explicitly.

Raises:

Type	Description
ConvergenceError	If a concrete report indicates failure.

residual_report

residual_report(
    residual: Array,
    tol: float = 1e-08,
    *,
    iterations: Array | int = 0,
    step_norm: Array | float = 0.0,
    failure: SolveStatus = MAX_ITERATIONS,
) -> SolveReport

Grade a scaled residual independently of an iteration's stopping rule.

Implicit solvers

implicit

Implicit differentiation of fixed-point solvers.

Phase-equilibrium calculations are iterative: a flash, a bubble point, or a saturation pressure is the solution of a fixed-point or root-finding loop. Back-propagating through the individual iterations would be wasteful and numerically noisy. Instead Fugacio differentiates the converged solution directly, via the implicit function theorem.

For a fixed point x* = g(x*, theta) the sensitivity to the parameters theta satisfies::

(I - dg/dx) dx*/dtheta = dg/dtheta

so a linear solve of the residual Jacobian yields implicit sensitivities regardless of how many iterations the forward solve took. This is the same trick used by the cubic-root fugacio.thermo.eos.compress_factor, generalized to vector unknowns.

Functions:

Name	Description
newton_root	Solve a scalar residual(x, params) = 0 by a damped Newton iteration.
gate_derivative	Return value unchanged, with nonfinite derivatives unless valid.
gate_tree	Apply gate_derivative to every floating leaf of a pytree.
implicit_solution	Attach an implicit derivative to a detached, independently solved root.
newton_system_with_info	Solve a square residual system and report convergence independently.
fixed_point_with_info	Converge a contraction and solve its derivative as a linear system.
fixed_point	Return a contraction's fixed point with implicit forward/reverse derivatives.
bracketed_root_with_info	Bisect a validated scalar bracket and check the resulting residual.
scanned_root_with_info	Locate the first increasing sign change of residual on a grid, then bisect it.
bracketed_root	Return a validated bisection root with implicit forward and reverse derivatives.
newton_root_with_info	Solve a scalar root with residual-decreasing steps and a checked report.

newton_root

newton_root(
    residual: ResidualFn,
    params: Any,
    x0: Array,
    tol: float = 1e-12,
    max_iter: int = 100,
    damping: float = 1.0,
) -> Array

Solve a scalar residual(x, params) = 0 by a damped Newton iteration.

The forward Newton step uses the autodiff slope dr/dx and an optional damping (step multiplier in (0, 1]) for stability; the converged root is differentiated with respect to params by the implicit function theorem (the iteration itself is not traced).

This is the inner-loop primitive for density branches, where a solve is repeated inside every property evaluation and the caller selects among branches. It doesn't verify convergence: callers must check the residual at the returned point (as the SAFT and reference-fluid branch selectors do) before accepting it. Use newton_root_with_info for a checked scalar root.

Returns:

Type	Description
Array	The final iterate; differentiable with respect to params.

gate_derivative

gate_derivative(value: Array, valid: Array) -> Array

Return value unchanged, with nonfinite derivatives unless valid.

Use this when a solve's acceptance depends on checks made after its iteration (a trivial solution, a domain limit, a physical audit). The primal is retained for diagnosis; downstream optimizers can't consume its sensitivity.

gate_tree

gate_tree(tree: Any, valid: Array) -> Any

Apply gate_derivative to every floating leaf of a pytree.

implicit_solution

implicit_solution(
    residual: ResidualFn,
    value: Array,
    theta: Any,
    valid: Array,
    jacobian: JacobianFn
    | Literal["sequential"]
    | None = None,
) -> Array

Attach an implicit derivative to a detached, independently solved root.

The initial guess and iteration history carry no derivative. Both forward and reverse differentiation solve the linearized residual system. A failed primal has a nonfinite sensitivity, preventing optimization from silently consuming derivatives at an unconverged iterate. jacobian optionally supplies an exact structured linearization. The same operator is used for forward derivatives and its transposed adjoint. jacobian="sequential" builds a dense matrix one direction at a time, sharing one residual linearization between state and parameter directions. This avoids repeated nested unit linearizations in a recycle residual.

newton_system_with_info

newton_system_with_info(
    residual: ResidualFn,
    x0: Array,
    theta: Any,
    tol: float = 1e-10,
    max_iter: int = 50,
    *,
    scale: Array | None = None,
    residual_scale: Array | None = None,
    lower: Array | None = None,
    upper: Array | None = None,
    jacobian: JacobianFn | None = None,
    stall_tol: float = 0.0,
    gate: bool = True,
) -> SolveResult

Solve a square residual system and report convergence independently.

Parameters:

Name	Type	Description	Default
residual	ResidualFn	Vector residual F(x, theta) with the same shape as x.	required
x0	Array	Starting vector; may be a previously converged solution.	required
theta	Any	Differentiable parameters. Pass all varying quantities here.	required
tol	float	Maximum scaled residual accepted as converged.	1e-10
max_iter	int	Maximum number of Newton steps.	50
scale	Array | None	Positive characteristic variable magnitudes; defaults to one.	None
residual_scale	Array | None	Positive equation scales; defaults to one.	None
lower	Array | None	Optional lower bounds used only during initialization and search.	None
upper	Array | None	Optional upper bounds used only during initialization and search.	None
jacobian	JacobianFn | None	Optional exact structured Jacobian factory (x, theta). It must describe the original residual, before the supplied scales.	None
stall_tol	float	Residual accepted when no Newton step can reduce it further, i.e. at the floating-point floor of a cancelling residual. Zero (the default) never accepts a stalled iterate.	0.0
gate	bool	Give an unconverged root nonfinite derivatives. Pass False only when the caller selects among roots and gates the result.	True

Returns:

Type	Description
SolveResult	Best iterate and a :class:SolveReport. Bounds and scales are numerical
SolveResult	aids, not additional equations. Sensitivities are defined only when the
SolveResult	original residual converges to a locally nonsingular root.

Raises:

Type	Description
ValueError	If the tolerance or iteration cap is invalid.

fixed_point_with_info

fixed_point_with_info(
    g: Callable[[Array, Any], Array],
    x0: Array,
    theta: Any,
    tol: float = 1e-12,
    max_iter: int = 200,
) -> SolveResult

Converge a contraction and solve its derivative as a linear system.

The derivative does not repeat the forward fixed-point iteration, so a slowly converging adjoint cannot silently exhaust a separate iteration cap.

fixed_point

fixed_point(
    g: Callable[[Array, Any], Array],
    x0: Array,
    theta: Any,
    tol: float = 1e-12,
    max_iter: int = 200,
) -> Array

Return a contraction's fixed point with implicit forward/reverse derivatives.

Use :func:fixed_point_with_info to inspect termination. An unconverged iterate has nonfinite derivatives.

bracketed_root_with_info

bracketed_root_with_info(
    residual: ResidualFn,
    params: Any,
    lo: Array,
    hi: Array,
    tol: float = 1e-12,
    max_iter: int = 200,
    *,
    residual_tol: float = 1e-08,
    relative_tol: float = 1e-06,
    gate: bool = True,
) -> SolveResult

Bisect a validated scalar bracket and check the resulting residual.

The forward pass uses only residual values, so it is robust through the poles and kinks that scalar equilibrium residuals (bubble/dew temperature, Underwood roots, saturation lines) exhibit near bracket ends. tol stops the bisection by bracket width; acceptance is judged independently from the residual at the returned point. A root is converged only when

• both endpoint residuals are finite and differ in sign (a bracket without a sign change is reported as INVALID_INPUT rather than returning an endpoint), and
• |residual(root)| <= max(residual_tol, relative_tol * min(|f(lo)|, |f(hi)|)), so a sign change across a pole or jump, whose residual stays large as the bracket shrinks, is reported as a failure. The residual of the last bisected point is reused for that check: it lies within tol of the returned root, and evaluating it again would compile a second copy of the residual (an exchanger or column solve).

Endpoints that already solve the equation take zero iterations.

Parameters:

Name	Type	Description	Default
residual	ResidualFn	Scalar function residual(x, params) -> r.	required
params	Any	Differentiable parameter pytree forwarded to residual.	required
lo	Array	Lower bracket; carries no derivative.	required
hi	Array	Upper bracket; carries no derivative.	required
tol	float	Bracket-width stopping tolerance.	1e-12
max_iter	int	Bisection iteration cap.	200
residual_tol	float	Absolute residual floor, in the residual's units.	1e-08
relative_tol	float	Residual tolerance relative to the smaller endpoint residual.	1e-06
gate	bool	Give an unconverged root nonfinite derivatives. Pass False only when the caller selects among several roots and gates the selected value itself (a zero cotangent times a NaN gate is still NaN).	True

Returns:

Type	Description
SolveResult	The best root estimate and its report; the value is differentiable with
SolveResult	respect to params by the implicit function theorem.

Raises:

Type	Description
ValueError	If tolerances or the iteration limit are invalid.

scanned_root_with_info

scanned_root_with_info(
    residual: ResidualFn,
    params: Any,
    lo: Array,
    hi: Array,
    tol: float = 1e-12,
    max_iter: int = 200,
    *,
    points: int = 12,
    residual_tol: float = 1e-08,
    relative_tol: float = 1e-06,
) -> SolveResult

Locate the first increasing sign change of residual on a grid, then bisect it.

Saturation residuals are often undefined (NaN) over part of a physically wide bracket: a bubble pressure doesn't exist above a mixture's critical region, and a reference state doesn't exist above a component's critical temperature. A grid of points values finds the first cell whose ends are finite with f <= 0 < f; that cell is bisected by bracketed_root_with_info. No such cell is reported as INVALID_INPUT.

Returns:

Type	Description
SolveResult	The root and its report, differentiable with respect to params.

Raises:

Type	Description
ValueError	If fewer than two grid points are requested.

bracketed_root

bracketed_root(
    residual: ResidualFn,
    params: Any,
    lo: Array,
    hi: Array,
    tol: float = 1e-12,
    max_iter: int = 200,
) -> Array

Return a validated bisection root with implicit forward and reverse derivatives.

The value is the best estimate of :func:bracketed_root_with_info. An invalid bracket or an unconverged residual gives nonfinite derivatives; use the _with_info form to accept or reject the value itself.

newton_root_with_info

newton_root_with_info(
    residual: ResidualFn,
    params: Any,
    x0: Array,
    tol: float = 1e-10,
    max_iter: int = 100,
    *,
    lower: Array | None = None,
    upper: Array | None = None,
) -> SolveResult

Solve a scalar root with residual-decreasing steps and a checked report.

Structured linear algebra

linear

Structured Jacobians and checked, implicitly differentiated linear solves.

The block representation stores a nearest-neighbor chain and a small dense border. Coloring assembles its Jacobian without a stage-sized autodiff batch. Block elimination pivots within each diagonal block; a residual check selects a pivoted dense fallback when that elimination isn't adequate. The fallback changes the linear algorithm, never the nonlinear equations.

Classes:

Name	Description
LinearReport	Independent linear residual and information about the actual algorithm.
LinearResult	Solution and residual evidence; a rejected solution has nonfinite values.
Jacobian	Matrix operations required by Newton and implicit differentiation.
DenseJacobian	A dense Jacobian using JAX's pivoted linear solve.
BorderedBlockJacobian	Block tridiagonal matrix with a small, unrestricted border.
BlockLayout	Static declaration of a chain's exact Jacobian sparsity pattern.

Functions:

Name	Description
dense_jacobian	Linearize a square residual, optionally evaluating directions sequentially.

LinearReport

Bases: NamedTuple

Independent linear residual and information about the actual algorithm.

Methods:

Name	Description
to_dict	Return concrete JSON data, retaining a nonfinite residual as null.

to_dict

to_dict() -> dict[str, Any]

Return concrete JSON data, retaining a nonfinite residual as null.

LinearResult

Bases: NamedTuple

Solution and residual evidence; a rejected solution has nonfinite values.

Jacobian

Bases: Protocol

Matrix operations required by Newton and implicit differentiation.

Methods:

Name	Description
solve	Solve for one or several right-hand sides.
scaled	Return diag(rows) @ self @ diag(columns).
to_dense	Materialize the matrix for a reference solve or diagnosis.

solve

solve(rhs: Array) -> Array

Solve for one or several right-hand sides.

scaled

scaled(rows: Array, columns: Array) -> Jacobian

Return diag(rows) @ self @ diag(columns).

to_dense

to_dense() -> Array

Materialize the matrix for a reference solve or diagnosis.

DenseJacobian

Bases: NamedTuple

A dense Jacobian using JAX's pivoted linear solve.

Methods:

Name	Description
solve	Solve using the differentiable dense reference implementation.
scaled	Apply independent equation and unknown scales.
to_dense	Return the stored matrix.

solve

solve(rhs: Array) -> Array

Solve using the differentiable dense reference implementation.

scaled

scaled(rows: Array, columns: Array) -> DenseJacobian

Apply independent equation and unknown scales.

to_dense

to_dense() -> Array

Return the stored matrix.

BorderedBlockJacobian

Bases: NamedTuple

Block tridiagonal matrix with a small, unrestricted border.

lower, diagonal, and upper have shape (n, b, b). The first lower and last upper blocks are zero. Border columns have shape (n, b, k), border rows (k, n, b), and the corner (k, k). Unknowns and residuals are ordered by block, followed by border entries.

Methods:

Name	Description
matvec	Multiply a vector or a matrix of right-hand sides without densifying.
transpose	Transpose both the chain and its border.
scaled	Apply row and column scales while preserving the sparsity pattern.
to_dense	Materialize a reference matrix without differentiating the residual again.
solve_with_info	Solve with implicit gradients and independently checked linear residuals.
solve	Return the checked, differentiable linear solution.

Attributes:

Name	Type	Description
shape	tuple[int, int]	Square matrix shape.

shape property

shape: tuple[int, int]

Square matrix shape.

matvec

matvec(x: Array) -> Array

Multiply a vector or a matrix of right-hand sides without densifying.

transpose

transpose() -> BorderedBlockJacobian

Transpose both the chain and its border.

scaled

scaled(
    rows: Array, columns: Array
) -> BorderedBlockJacobian

Apply row and column scales while preserving the sparsity pattern.

to_dense

to_dense() -> Array

Materialize a reference matrix without differentiating the residual again.

solve_with_info

solve_with_info(rhs: Array) -> LinearResult

Solve with implicit gradients and independently checked linear residuals.

The transposed solve has its own factorization and residual check. Neither the block pivot choices nor dense fallback choices are differentiated. Singular or unresolved equations return NaNs.

solve

solve(rhs: Array) -> Array

Return the checked, differentiable linear solution.

BlockLayout dataclass

BlockLayout(
    blocks: int, block_size: int, border_size: int = 0
)

Static declaration of a chain's exact Jacobian sparsity pattern.

Each core residual block may depend on its own and neighboring unknown blocks and every border unknown. Border equations may depend on every unknown. This is a structural contract, not a pattern inferred from zeros at one operating point. Use check against a dense Jacobian when adding a model adapter.

Methods:

Name	Description
linearize	Assemble the exact block Jacobian using colored JVPs and border VJPs.
check	Compare declared structure with dense autodiff at a concrete point.

Attributes:

Name	Type	Description
size	int	Total number of unknowns and equations.
directions	int	Forward coloring directions, independent of chain length beyond three blocks.

size property

size: int

Total number of unknowns and equations.

directions property

directions: int

Forward coloring directions, independent of chain length beyond three blocks.

linearize

linearize(
    residual: Callable[..., Array], x: Array, theta: Any
) -> BorderedBlockJacobian

Assemble the exact block Jacobian using colored JVPs and border VJPs.

check

check(
    residual: Callable[..., Array],
    x: Array,
    theta: Any,
    *,
    tolerance: float = 1e-09,
) -> dict[str, Any]

Compare declared structure with dense autodiff at a concrete point.

This diagnostic detects missing couplings at the supplied point. It isn't a proof of a pattern's validity throughout the model domain.

dense_jacobian

dense_jacobian(
    residual: Callable[..., Array],
    x: Array,
    theta: Any,
    *,
    vectorize: bool = True,
) -> DenseJacobian

Linearize a square residual, optionally evaluating directions sequentially.

vectorize=False retains one primal linearization and maps over tangent directions with a compiled loop. Nested unit solves then needn't acquire another tangent batch dimension. The resulting matrix and solve remain dense; this option bounds derivative working storage, not matrix storage.

sparsity

Declared residual incidence, deterministic coloring, and structural diagnostics.

Classes:

Name	Description
SparsityPattern	Conservative row dependencies in a fixed scalar unknown order.

SparsityPattern dataclass

SparsityPattern(
    columns: int, rows: tuple[tuple[int, ...], ...]
)

Conservative row dependencies in a fixed scalar unknown order.

Every true derivative entry must appear in rows. Extra entries cost work but don't change the equations. Never infer this contract from zeros at a single operating point. Unknown custom blocks should declare every variable until their dependencies are known.

Methods:

Name	Description
coloring	Color the column-intersection graph in a deterministic greedy order.
jacobian	Create an exact colored AD assembler with a dense reference solve.
matching	Find a maximum equation-to-variable matching; unmatched rows contain -1.
diagnose	Describe incidence and unmatched equations/variables without a numerical solve.

coloring

coloring() -> tuple[int, ...]

Color the column-intersection graph in a deterministic greedy order.

jacobian

jacobian(
    residual: Callable[..., Array],
) -> Callable[[Array, Any], DenseJacobian]

Create an exact colored AD assembler with a dense reference solve.

The result still stores a dense global matrix. This adapter reduces differentiation work for sparse flowsheet connections; it doesn't claim a sparse factorization or expand a procedural unit's equations.

matching

matching() -> tuple[int, ...]

Find a maximum equation-to-variable matching; unmatched rows contain -1.

Breadth-first augmenting paths avoid Python recursion limits on long process trains. A complete matching is only a structural upper bound on numerical rank, especially with conservative block dependencies.

diagnose

diagnose(
    *,
    equations: Sequence[str] | None = None,
    variables: Sequence[str] | None = None,
) -> dict[str, Any]

Describe incidence and unmatched equations/variables without a numerical solve.

Reusable sensitivities

sensitivity

Reusable local derivatives of a converged numerical calculation.

One linearization evaluates the primal once and retains the residual data needed by its JVP. Its transpose supplies VJPs without retracing the primal. The object belongs to one operating point; construct a new one when that point, its model parameters, or its initialization changes.

Classes:

Name	Description
Linearization	A value and reusable JVP/VJP at one array-valued operating point.

Functions:

Name	Description
derivative_strategy	Validate derivative options and describe the selected direction count.
linearize	Evaluate once and retain a reusable local derivative, optionally with reports.

Linearization dataclass

Linearization(
    point: Array,
    value: Array,
    auxiliary: Any,
    _push: Callable[[Array], Array],
)

A value and reusable JVP/VJP at one array-valued operating point.

Create this object with linearize. The stored callables are local linear maps, not nonlinear approximations valid at subsequent points. Auxiliary data, such as solve reports, isn't differentiated. Releasing the object releases its retained residuals; no global point cache exists.

Methods:

Name	Description
block_until_ready	Wait for both outputs and the residual buffers retained by the JVP.
jvp	Apply the local derivative to a tangent with the input's shape.
vjp	Apply the transposed derivative to an output cotangent.
jacobian	Build a Jacobian using bounded batches of the cheaper orientation.

block_until_ready

block_until_ready() -> Linearization

Wait for both outputs and the residual buffers retained by the JVP.

JAX's returned linear map is a registered partial pytree whose leaves include those buffers. Waiting for the primal output alone can miss asynchronous work that only a later derivative application consumes. This method also lets jax.block_until_ready synchronize this object.

jvp

jvp(tangent: Array) -> Array

Apply the local derivative to a tangent with the input's shape.

vjp

vjp(cotangent: Array) -> Array

Apply the transposed derivative to an output cotangent.

jacobian

jacobian(
    *, mode: str = "auto", batch_size: int = 1
) -> Array

Build a Jacobian using bounded batches of the cheaper orientation.

The result has shape value.shape + point.shape, matching JAX's Jacobian convention. Forward and reverse modes use the same retained primal and the implicit derivative rules of the underlying solvers.

derivative_strategy

derivative_strategy(
    inputs: int,
    outputs: int,
    *,
    mode: str = "auto",
    batch_size: int = 1,
) -> dict[str, Any]

Validate derivative options and describe the selected direction count.

auto uses the smaller of the input and output dimensions, preferring forward mode on a tie. This is a direction-count heuristic, not a promise that one orientation is faster for every property package. A batch size of one keeps tangent storage bounded for large nested process models.

linearize

linearize(
    function: Callable[..., Any],
    point: Array,
    *,
    has_aux: bool = False,
) -> Linearization

Evaluate once and retain a reusable local derivative, optionally with reports.

Inputs and differentiated outputs must be nonempty real floating arrays. Use an explicit vector adapter for parameter dictionaries or metric trees. The calculation needn't be wrapped in a single JIT; compiled unit kernels can retain their own compilation boundaries in a larger flowsheet.