Solver Algorithm Comparison

Three equilibrium solver algorithms were analysed before deciding on the implementation strategy for Prometheus. This document records the analysis, explains why Gordon-McBride (G-McB) is the production default today, and documents the specific differences between Prometheus’s implementation and the two reference codes: pep-for-zos (Villars-Browne PEP) and cpropep/RocketCEA (Gordon-McBride).

The Three Candidate Algorithms

Gordon-McBride (G-McB)

Source: NASA RP-1311, Gordon & McBride (1994).

Reference implementation: cpropep/libcpropep/src/equilibrium.c, RocketCEA.

The G-McB method solves a modified Newton system where the primary unknowns are the Lagrange multipliers πₖ (one per element), the condensed-species mole corrections Δnc, the correction to total gas moles Δln(n), and — for HP/SP problems — the temperature correction Δln(T). Gas-phase species mole corrections are computed analytically after each Newton step:

Δln(nⱼ) = Σₖ aₖⱼ·πₖ  +  Δln(n)  −  μⱼ   [+ H°ⱼ/RT·Δln(T) for HP/SP]

This analytical substitution keeps the Newton matrix small: (S + nc + 1) for TP problems and (S + nc + 2) for HP/SP, regardless of the total number of product species. For H₂/O₂ with no condensed phases this is a 3×3 system; for a solid-motor problem with one condensed phase (Al₂O₃) it becomes 4×4.

Temperature handling: Δln(T) is an unknown in the same Newton system as the composition. T and composition converge simultaneously. This is elegant but means the iteration can diverge if the starting guess is far from the adiabatic temperature.

Condensed phases: managed by a restart strategy — a condensed species is added or removed and the Newton loop is restarted with a counter that prevents cycling.

Villars-Browne PEP

Source: Cruise, NWC Technical Publication 6037 (1973/1979).

Reference implementation: pep-for-zos/PEP (Fortran 77/IBM z/OS port).

PEP selects an optimised basis of S species (one per element) using Browne’s Gram-Schmidt rank test. All non-basis species are expressed as stoichiometric combinations of the basis via the reaction-coefficient matrix ν = C·B⁻¹. The algorithm then iterates on reaction extents Δζᵢ:

ln Kᵢ  = Σⱼ ν[i,j]·g°ⱼ/RT  −  g°ᵢ/RT        (log equilibrium constant)
ln Qᵢ  = γᵢ·ln(P·nᵢ/n)  −  Σⱼ γⱼ·ν[i,j]·ln(P·nⱼ/n) (log reaction quotient)
Δζᵢ    = (ln Kᵢ − ln Qᵢ) / (γᵢ/nᵢ + Σⱼ γⱼ·ν[i,j]²/nⱼ)

(γⱼ = 1 for gas, 0 for condensed — pure condensed species have unit activity.)

Basis species moles are corrected after each non-basis adjustment to maintain element balance.

Temperature handling: separate outer Newton + interval-halving loop that calls the TP inner loop at each trial temperature. This guarantees convergence as long as the solution lies in [298 K, 6000 K].

Convergence rate: each iteration is a first-order (linear) correction — typically 5–15 inner iterations per temperature step.

pep-for-zos specifics: VLNK is computed using the full chemical potential (H − TS including mixing terms), ensuring VLNK → 0 at equilibrium regardless of sign convention. Basis swaps use a TABLO linear-algebra tableau pivot. Damping is adaptive: VQQ = max(0.05, 0.5 − (JC−1)/20), starting at 50% and tightening with iteration count.

Why PEP is strictly inferior to Hybrid:

  • Same S×S matrix size as Hybrid.

  • Linear rather than quadratic convergence.

  • For H₂/O₂ at 3000 K, PEP typically requires 10–20 iterations; Hybrid converges in 4–6.

  • No simplicity advantage: the reaction-coefficient machinery (ν = C·B⁻¹, basis selection) is needed by Hybrid anyway.

PEP is therefore not implemented as a production solver.

Hybrid Solver

Source: Prometheus (original design).

Concept: combines PEP’s reaction-adjustment architecture with G-McB’s Newton-step quadratic convergence.

The key insight is that any set of S linearly independent species can serve as the primary Newton unknowns — the remaining species follow analytically. Hybrid selects the S “major” species (highest mole amounts, Browne’s test) and forms a compressed Newton system in those S unknowns. Minor species are updated analytically from the element potentials π extracted from the Newton solve, identically to the G-McB gas-correction formula:

Minor-species update (no iteration needed):
  ln(nⱼ) = Σₖ πₖ·A[j,k]  −  g°ⱼ/RT  −  ln(P/P°)  +  ln(n_gas)

  • Matrix size: S×S (identical to PEP).

  • Convergence rate: quadratic (identical to G-McB).

  • Temperature handling: outer Newton + interval-halving (identical to PEP, more robust than G-McB’s embedded-T approach).

  • Condensed phases: phase-deletion rule (same as PEP, no restart required).

Comparison Table

Feature

G-McB

PEP

Hybrid

Primary unknowns

π (Lagrange mult.)

reaction extents Δζ

π for major species

Matrix size (TP)

(S + nc + 1)²

S×S

S×S

Matrix size (HP/SP)

(S + nc + 2)²

S×S + outer loop

S×S + outer loop

Convergence rate

quadratic

linear

quadratic

Temperature

embedded in Newton

outer interval-halving

outer interval-halving

T convergence guarantee

none

yes

yes

Condensed handling

restart Newton

phase deletion

phase deletion

Reference code

cpropep, RocketCEA

pep-for-zos

—

Implementation status

complete

not implemented

complete

For H₂/O₂ combustion (S = 2, nc = 0):

  • G-McB: 3×3 matrix, no outer loop.

  • Hybrid: 2×2 matrix, outer loop (fast).

  • PEP: 2×2 matrix, outer loop (slow — linear convergence).

For a solid-motor problem (S = 3 elements: Al, H, O; nc = 1: Al₂O₃(s)):

  • G-McB: 5×5 matrix.

  • Hybrid: 3×3 matrix.

  • PEP: 3×3 matrix (linear convergence).

Why G-McB Is The Preferred Default

G-McB is both the preferred and fastest solver in the current Prometheus implementation for three practical reasons:

  1. Reference code available. cpropep/libcpropep/src/equilibrium.c implements G-McB with the exact same RP-1311 formulas, providing a function-by-function reference for debugging.

  2. Used by RocketCEA. The torture-test benchmark compares Prometheus against RocketCEA output. Implementing G-McB first ensures we can validate composition, temperature, and performance numbers against the same underlying algorithm.

  3. Current runtime profile. In the present codebase, G-McB delivers the best end-to-end runtime for the standard benchmark workloads while preserving robust convergence behaviour.

Differences: Prometheus G-McB vs cpropep

Aspect

cpropep

Prometheus G-McB

Language

C

Python / NumPy

Rank analysis

simple rank() check

QR with column pivoting

Independent elements

drop elements with zero b₀

QR pivot gives maximally independent subset

Condensed inclusion

most-negative μ°c/RT − Σπa

same criterion, same π source

Convergence criterion

5×10⁻⁶ (matches RP-1311)

5×10⁻⁶ (same)

Temperature handling

embedded in Newton (HP/SP)

same: Δln(T) as Newton unknown

Initial moles

equal among gas species

equal among gas species

Damping

λ₁ caps Δln to ≤ 2; λ₂ trace floor

same two-step strategy

Note

fill_matrix in cpropep builds the composition-search system (what Prometheus calls _assemble_jacobian). fill_equilibrium_matrix is a separate function called only for HP/SP that adds the Δln(T) row and column. Prometheus merges both into a single _assemble_jacobian call that is aware of problem_type.

Differences: Prometheus Hybrid vs pep-for-zos PEP

Aspect

pep-for-zos

Prometheus Hybrid

Inner loop

reaction-extent Δζ (linear)

compressed Newton on major species (quadratic)

Minor species update

K/Q ratio, every 4th pass

analytical from π, every pass

Basis swap

TABLO tableau pivot

Gram-Schmidt re-selection

Damping

adaptive VQQ

RP-1311-style λ₁/λ₂

Mole space

linear

log-space for stability

Temperature

outer HBAL/SBAL loop

outer _temperature_search Newton + halving

VLNK convention

full μ (includes mixing)

not applicable

Summary

Solver

Role

Status

GordonMcBrideSolver

Production default / RocketCEA comparison

Complete

HybridSolver

Alternative solver

Complete

PEPSolver

Not implemented (strictly inferior to Hybrid)

Deferred indefinitely

GordonMcBrideSolver is the recommended solver for production use and is the current default throughout Prometheus. It is the best-performing solver in the current implementation and matches the CEA algorithm used by RocketCEA, enabling direct comparison of composition, temperature, and performance metrics. HybridSolver remains available as an alternative implementation for algorithmic comparison and future optimisation work.