prometheus_equilibrium.equilibrium.solver

Equilibrium solvers for Prometheus.

Three concrete solvers are provided, all implementing the same EquilibriumSolver interface. They share the same inputs (EquilibriumProblem) and outputs (EquilibriumSolution) and differ only in the numerical algorithm used internally.

Solver summary

MajorSpeciesSolver (recommended — default)

Compressed Newton on major species combined with an analytical minor- species update and a PEP-style outer temperature search. Combines the quadratic convergence of Newton’s method with PEP’s small matrix size and guaranteed temperature convergence. See MajorSpeciesSolver.

PEPSolver (reference implementation)

Pure Villars-Browne reaction-adjustment method as described in the original NWC Propellant Evaluation Program (Cruise, NWC TP 6037, 1979). Converges linearly; useful as a simpler reference and for validating the more complex Major-Species solver. See PEPSolver.

GordonMcBrideSolver (comparison / validation)

Modified Lagrange multiplier method from NASA RP-1311 (Gordon & McBride, 1994). The algorithm used by CEA. Kept for direct comparison against RocketCEA output. See GordonMcBrideSolver.

Algorithm relationships

All three solvers minimise the same objective — Gibbs free energy subject to element-balance constraints — and find the same solution. They differ in how they parameterise and solve the resulting nonlinear system:

Feature	G-McB	PEP	Major-Species
Core variables	π (Lagrange mult)	nⱼ directly	π for major
Matrix size	(S+cnd+1)²	S×S	S×S
Minor species	log-space Newton	Kᵢ/Qᵢ ratio	analytic from π
Convergence rate	quadratic	linear	quadratic
T convergence	no guarantee	interval halving	interval halving
Condensed phases	restart loop	phase deletion	phase deletion

where S = number of chemical elements, cnd = number of active condensed phases. For H₂/O₂ combustion (S=2, no condensed phases) the Major-Species and PEP matrix is always 2×2 regardless of the number of candidate product species.

Shared infrastructure

PEP and Major-Species solvers both use:

• ElementMatrix.select_basis() — Browne’s optimised basis selection

• ElementMatrix.reaction_coefficients() — ν = C·B⁻¹

• _ReactionAdjustmentBase._temperature_search() — Newton + interval halving

The common base class _ReactionAdjustmentBase holds this shared logic so it is not duplicated.

Classes

EquilibriumSolver([max_iterations, ...])	Abstract base class for all equilibrium solvers.
GordonMcBrideSolver([max_iterations, ...])	Gibbs free energy minimisation via modified Lagrange multipliers.
PEPSolver([max_iterations, tolerance, ...])

class prometheus_equilibrium.equilibrium.solver.EquilibriumSolver(max_iterations: int = 50, tolerance: float = 5e-06, capture_history: bool = True, history_stride: int = 1)

Bases: ABC

Abstract base class for all equilibrium solvers.

Concrete implementations must override solve(). The interface is solver-agnostic: any of the three concrete solvers can be swapped in without changing calling code.

Parameters

max_iterationsint

Maximum inner (composition) iterations before declaring non-convergence.

tolerancefloat

Convergence threshold for species mole-fraction residuals.

abstractmethod solve(problem: EquilibriumProblem, guess: Mixture | None = None, log_failure: bool = True) → EquilibriumSolution

Run the equilibrium calculation and return the converged state.

Parameters

problemEquilibriumProblem

Fully specified problem. Call problem.validate() before passing to catch ill-posed inputs early.

guessMixture, optional

Initial composition guess. If provided, the solver uses these mole amounts as a starting point instead of the default monatomic initialization.

Returns

EquilibriumSolution

The converged state. Always check .converged before using results.

class prometheus_equilibrium.equilibrium.solver.PEPSolver(max_iterations: int = 50, tolerance: float = 5e-06, minor_threshold: float = 0.01, capture_history: bool = True, history_stride: int = 1)

Bases: _ReactionAdjustmentBase

solve(problem: EquilibriumProblem, guess: Mixture | None = None, log_failure: bool = True) → EquilibriumSolution

Run the equilibrium iteration.

Parameters:

• problem – Fully specified equilibrium problem.

• guess – Optional initial composition guess.

Returns:

EquilibriumSolution with the converged composition.

class prometheus_equilibrium.equilibrium.solver.MajorSpeciesSolver(max_iterations: int = 50, tolerance: float = 5e-06, minor_threshold: float = 0.01, capture_history: bool = True, history_stride: int = 1)

Bases: _ReactionAdjustmentBase

Compressed Newton + analytical minor-species update.

Combines the best features of the Gordon-McBride and PEP algorithms:

• PEP basis selection keeps the Newton system small (S×S for major species only, where S = number of elements).

• Newton quadratic convergence on the compressed system instead of PEP’s linear serial corrections.

• Analytical minor-species update from the converged element potentials — exact at any mole-fraction scale, no log-space artefacts.

• PEP outer temperature search with interval-halving fallback.

This is the recommended default solver.

solve(problem: EquilibriumProblem, guess: Mixture | None = None, log_failure: bool = True) → EquilibriumSolution

Run the major-species equilibrium iteration.

Parameters:

• problem – Fully specified equilibrium problem.

• guess – Optional initial composition guess.

Returns:

EquilibriumSolution with the converged composition.

class prometheus_equilibrium.equilibrium.solver.GordonMcBrideSolver(max_iterations: int = 300, tolerance: float = 0.0001, capture_history: bool = True, history_stride: int = 1)

Bases: EquilibriumSolver

Gibbs free energy minimisation via modified Lagrange multipliers.

Implements the full algorithm from NASA RP-1311 (Gordon & McBride, 1994), using cpropep as a line-by-line reference. This is also the algorithm inside RocketCEA, enabling direct numerical validation.

Algorithm overview

The equilibrium condition for gas-phase species j is:

μⱼ/RT = μⱼ°(T)/RT + ln(nⱼ/n_gas) + ln(P/P°) = Σₖ πₖ·A[j,k]

where πₖ are dimensionless element potentials (Lagrange multipliers in units of RT). At each Newton step the system is linearised, yielding a small augmented matrix whose unknowns are [π, Δnc, Δln(n), Δln(T)]:

• π (shape S): solved fresh at each step — not accumulated.

• Δnc (shape nc_active): linear corrections to condensed mole amounts.

• Δln(n): correction to total gas moles.

• Δln(T): temperature correction (HP/SP only).

Matrix size is (S + nc_active + 1) for TP or (S + nc_active + 2) for HP/SP, where S = n_elements and nc_active = currently present condensed species. Gas species are updated analytically from π after each solve.

After each Newton step the condensed phase set is managed:

• Remove: any condensed with nⱼ ≤ 0 is deactivated.

• Include: any inactive condensed with μ°c/RT < Σₖ πₖ·aₖc is a candidate; the most favourable is activated.

Temperature converges in the same Newton loop as composition (HP/SP), not in a separate outer loop.

param max_iterations:

Maximum Newton iterations before declaring non-convergence.

param tolerance:

Convergence threshold: max abs(nⱼ·Δln(nⱼ)) / n_gas (and similar criteria for condensed, total-moles, and temperature).

solve(problem: EquilibriumProblem, guess: Mixture | None = None, log_failure: bool = True) → EquilibriumSolution

Run the Gordon-McBride iteration and return the converged state.

Parameters:

• problem – Fully specified equilibrium problem.

• guess – Optional initial composition guess. If provided, its mole amounts are used as the starting point instead of the default monatomic initialization.

Returns:

EquilibriumSolution with the converged composition.