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
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.

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.

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
bracketed_root	Solve a scalar residual(x, params) = 0 for x in [lo, hi] by bisection.
newton_root	Solve a scalar residual(x, params) = 0 by a damped Newton iteration.
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.
newton_system	Solve a vector root with implicit forward and reverse derivatives.
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.
newton_root_with_info	Solve a scalar root with residual-decreasing steps and a checked report.

bracketed_root

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

Solve a scalar residual(x, params) = 0 for x in [lo, hi] by bisection.

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) routinely exhibit at the bracket ends. The root is differentiated with respect to the parameter pytree params by the implicit function theorem in the custom_jvp rule below; the locators lo/hi carry no gradient.

Parameters:

Name	Type	Description	Default
residual	ResidualFn	Scalar function residual(x, params) -> r with a single sign change on [lo, hi].	required
params	Any	Differentiable parameter pytree forwarded to residual.	required
lo	Array	Lower bracket (residual must straddle zero across [lo, hi]).	required
hi	Array	Upper bracket.	required
tol	float	Absolute width of the final bracket.	1e-12
max_iter	int	Bisection iteration cap.	200

Returns:

Type	Description
Array	The bracketed root x*; differentiable with respect to params.

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). Prefer bracketed_root when a reliable bracket is available; newton_root is for smooth residuals where a good initial guess is cheap (saturation updates, Poynting corrections).

Returns:

Type	Description
Array	The root x*; differentiable with respect to params.

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,
) -> 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

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.

newton_system

newton_system(
    residual: ResidualFn,
    x0: Array,
    theta: Any,
    tol: float = 1e-10,
    max_iter: int = 50,
) -> Array

Solve a vector root with implicit forward and reverse derivatives.

Returns the best iterate for compatibility. Use :func:newton_system_with_info when accepting a result; a finite iterate alone does not prove convergence. Derivatives of an unconverged solution are nonfinite.

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,
) -> SolveResult

Bisect a validated scalar bracket and check the resulting residual.

tol limits bracket width; residual_tol independently limits the function residual in its own units. A sign change across a discontinuity can reduce the width without satisfying the equation and is reported as a failure. Endpoints that already solve the equation take zero iterations.

Raises:

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

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.