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.
HybridSolver([max_iterations, tolerance, ...])	GordonMcBrideSolver seeded by composition and a temperature estimate.
PEPSolver([max_iterations, tolerance, ...])

class prometheus_equilibrium.equilibrium.solver.EquilibriumSolver(max_iterations: int = 50, tolerance: float = 5e-06, capture_history: bool = False, 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 = False, 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 = False, 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.

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 = False, 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 and CEA as references.

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.

class prometheus_equilibrium.equilibrium.solver.HybridSolver(max_iterations: int = 300, tolerance: float = 0.0001, seed_max_iterations: int | None = None, seed_tolerance: float = 5e-06, capture_history: bool = False, history_stride: int = 1)

Bases: EquilibriumSolver

GordonMcBrideSolver seeded by composition and a temperature estimate.

Algorithm

1. Seed solve — run MajorSpeciesSolver on a TP sub-problem at problem.t_init. The inner Newton loop (S×S matrix) advances the gas composition away from the flat monatomic initial guess toward the chemical equilibrium at that temperature. Only the inner Newton loop runs — no outer temperature search — so the cost is bounded by seed_max_iterations.

2. Temperature correction — with the partial seed composition, compute one Newton step on the energy constraint:

   ΔT = −f / f′    where  f  = H(T_init) − H₀   (HP)
                           f′ = Cₚ(T_init)
   

   For the T correction to be accurate, the seed composition must be chemically realistic (close to equilibrium), which requires around 20 inner iterations for 3-element systems. For simple 2-element systems 10–12 iterations are sufficient.

3. Main solve — pass both the seed composition and the corrected t_init to GordonMcBrideSolver. Starting from a composition that is already near the equilibrium products eliminates the poorly-damped “build-up” phase of G-McB’s cold start.

Rationale

G-McB’s first iterations from a flat monatomic initialisation suffer heavy damping (λ ≪ 1) because the Newton step size is enormous: every dominant product (H₂O, CO₂, …) has to grow from ~1/N mole fraction to its true value while simultaneously finding the correct temperature. Convergence history for CH₄/O₂ shows T oscillating over a 700 K range for 19 iterations before entering the quadratic basin.

The key threshold is ~20 inner seed iterations for 3-element systems: below this G-McB still oscillates as if cold; at 20 inner iterations the dominant products are established and G-McB converges in ~6 iterations instead of ~24, giving roughly a 2× wall-clock speedup. For simpler 2-element systems (H₂/O₂) 10–12 iterations suffice and the speedup is smaller because G-McB is already fast.

param max_iterations:

Maximum Newton iterations for the G-McB main solve.

param tolerance:

Convergence tolerance for the G-McB main solve.

param seed_max_iterations:

Maximum inner iterations for the seed TP solve. 20 balances seed cost against G-McB iteration reduction across both 2-element (H₂/O₂) and 3-element (CH₄/O₂, N₂H₄/N₂O₄) propellant systems. Formal convergence of the TP solve is not required.

param seed_tolerance:

Convergence tolerance for the seed TP solve.

param capture_history:

Whether to record convergence history (G-McB phase).

param history_stride:

Only every n-th step is stored.

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

Solve by seeding GordonMcBrideSolver with a TP composition from MajorSpeciesSolver.

Parameters:

• problem – Fully specified equilibrium problem.

• guess – Optional initial composition guess (forwarded to the TP seed solve rather than to G-McB directly).

• log_failure – Whether to log convergence failures from G-McB.

Returns:

EquilibriumSolution from the G-McB step. The iterations field is the sum of seed TP iterations and G-McB iterations so that total work is reported accurately.