prometheus_equilibrium.equilibrium.element_matrix
ElementMatrix — stoichiometric matrix relating species to elements.
The element matrix A is the core linear-algebra object shared by all three
equilibrium solvers. Entry A[i, k] gives the number of atoms of element k
in species i (from the species’ elements dict).
Mathematical role
The equilibrium constraints are:
A^T · n = b₀
where n is the vector of species mole amounts and b₀ is the vector of total element abundances (fixed by the reactant composition). The element residual
Δb = b₀ − A^T · n
appears in the right-hand side of the element-balance equations in all three solvers.
Basis selection (Browne’s method — used by PEP and Hybrid solvers)
A basis is a subset of S species (S = number of elements) whose S×S composition sub-matrix B is non-singular. The optimised basis (Browne, 1960) picks the S species with the largest mole amounts that satisfy this condition, found by sorting species descending by mole amount and accepting candidates via the Gram-Schmidt orthogonality test.
The reaction-coefficient matrix ν = C · B⁻¹ expresses every non-basis species as a stoichiometric combination of basis species. Together B and ν are the foundation of the Villars reaction-adjustment iteration and the Hybrid solver’s compressed Newton system.
Charge balance
The electron pseudo-element "e-" is included as an ordinary column so
that ionic species are balanced automatically.
Singular-system handling (Gordon-McBride solver)
If the reactants do not supply every element represented in the product
species list, the corresponding column of A has no matching b₀ entry and
the Jacobian will be singular. independent_elements() identifies
the linearly independent subset of elements via QR decomposition, allowing
the G-McB solver to drop redundant constraints (RP-1311 §3.2).
Classes
|
Stoichiometric matrix A mapping species to element abundances. |
- class prometheus_equilibrium.equilibrium.element_matrix.ElementMatrix(species: List[Species], elements: List[str])
Bases:
objectStoichiometric matrix A mapping species to element abundances.
Parameters
- specieslist of Species
All species in the mixture (gas first, condensed second).
- elementslist of str
Ordered list of element symbols that define the columns of A. Typically derived from the union of all elements present in species.
- classmethod from_mixture(mixture: Mixture) ElementMatrix
Build from a
Mixture.The element list is the union of all elements across all species, sorted alphabetically (with
"e-"placed last if present).
- element_abundances(moles: ndarray) ndarray
Compute element abundances b = A^T · n [mol of element k].
Parameters
- molesarray-like, shape (n_species,)
Species mole amounts nⱼ.
Returns
- np.ndarray, shape (n_elements,)
Total moles of each element in the mixture.
- element_residuals(moles: ndarray, b0: ndarray) ndarray
Compute the element-balance residual Δb = b₀ − A^T · n.
This vector must be zero at equilibrium. It forms the right-hand side of the element-balance equations in all three solvers.
- Parameters:
moles – Current species mole amounts, shape (n_species,).
b0 – Target element abundances fixed by the reactant composition, shape (n_elements,).
- Returns:
Residual vector b₀ − A^T·n, shape (n_elements,).
- select_basis(moles: ndarray) Tuple[List[int], List[int]]
Select the optimised basis using Browne’s method (Gram-Schmidt).
The optimised basis is the set of S = n_elements species with the largest mole amounts whose S×S composition sub-matrix B is non-singular. It is found by sorting species in descending mole- amount order and accepting each candidate via the Gram-Schmidt orthogonality test (Cruise, NWC TP 6037, eq. 5).
At the start of each outer iteration the basis is re-optimised (since the composition has changed). This automatic re-selection serves two purposes:
It keeps the basis well-conditioned (dominant species are always in the basis, so B⁻¹ is numerically stable).
It relieves the user from choosing a basis manually.
Parameters
- molesarray-like, shape (n_species,)
Current mole amounts nⱼ. Species with nⱼ ≤ 0 are skipped.
Returns
- basis_indiceslist of int
Indices (into
self.species) of the S basis species, in the order they were accepted.- nonbasis_indiceslist of int
All remaining species indices, in original species order.
Raises
- RuntimeError
If fewer than S linearly independent species are found (the system is under-determined — reactant elements are not all representable).
Notes
The Gram-Schmidt test replaces the more expensive determinant test. A candidate row cₘ is orthogonalised against the current incomplete basis rows; if the result is the zero vector, the candidate is linearly dependent and rejected (Cruise eq. 5). The PEP paper also describes updating the basis via a linear-programming tableau pivot (Smith & Missen, 1968) to avoid recomputing B⁻¹ from scratch each iteration; this optimisation can be added once the basic loop works.
- basis_matrix(basis_indices: List[int]) ndarray
Return B = A[basis_indices, :], the S×S basis composition matrix.
B[j, k] = atoms of element k in the j-th basis species.
Parameters
- basis_indiceslist of int, length S
Row indices of the basis species in A.
Returns
np.ndarray, shape (S, S) where S = n_elements
- reaction_coefficients(basis_indices: List[int]) ndarray
Compute ν = C · B⁻¹, the reaction-coefficient matrix.
ν[i, j] is the stoichiometric coefficient of basis species j consumed when one mole of (non-basis) species i is formed from the basis. The chemical reaction for species i is:
Σⱼ ν[i,j] · basis_species[j] → species[i]
Used to:
Compute equilibrium constants (PEP solver):
ln Kᵢ = −Σⱼ ν[i,j] · gⱼ°/RTCorrect basis mole amounts after each stoichiometric adjustment (PEP eq. 10).
Parameters
basis_indices : list of int, length S
Returns
- np.ndarray, shape (n_species, n_elements)
Full ν matrix (rows for all species, including basis species themselves — their rows will be identity-like by construction).
- rank() int
Numerical rank of A (number of independent element constraints).
Computed via singular-value decomposition. A rank less than
n_elementsindicates redundant elements that will make the G-McB Jacobian singular if included naively.
- independent_elements() List[str]
Return the subset of element names forming a full-rank system.
Uses QR decomposition with column pivoting on Aᵀ to identify a maximal linearly independent set of element columns. Used by the G-McB solver to drop redundant element rows from the Jacobian (RP-1311 §3.2).
Returns
- list of str
Element symbols in a linearly independent ordering (most influential elements first, as determined by QR pivoting).
- reduced(active_elements: List[str] | None = None) ElementMatrix
Return a new ElementMatrix restricted to active_elements.
If active_elements is None,
independent_elements()is called automatically. Used by the G-McB solver to remove redundant constraints before building the Newton system.Parameters
- active_elementslist of str, optional
Element symbols to keep. Must be a subset of
self.elements.
Returns
- ElementMatrix
A new instance with the same species but only the specified element columns.