Repository navigation
Regressions from the Jan 2026 refactor (01f9ace, 231033e) #399
Description
Hello. These changes from #339 (released in v0.10.0 and later, and on development) appear to have changed some results by accident.
All numbers below use PC-SAFT with the parameter files shipped in this repo: parameters/pcsaft/esper2023.json (pure components) and parameters/pcsaft/rehner2023_binary.json (bip).
1. thermodynamic_factor drops x_i (crates/feos-core/src/state/residual_properties.rs:478, from 01f9ace)
The doc comment and docs/theory/eos/properties.md define Γ_ij = δ_ij + x_i (∂lnφ_i/∂x_j)_{T,p,Σ}. When the code switched to n_dln_phi_dnj(), the factor moles[i] (= x_i·N) was removed, but only N should have gone:
// v0.9 (before 01f9ace)
moles[i] * (dln_phi_dnj[(i, j)] - dln_phi_dnj[(i, n)]) + if i == j { 1.0 } else { 0.0 }
// current
dln_phi_dnj[(i, j)] - dln_phi_dnj[(i, n)] + if i == j { 1.0 } else { 0.0 }Check: for a binary, Γ_11 must not depend on which component is listed first. For ethanol/water, liquid at 300 K, 1 bar, x_ethanol = 0.2, v0.10.3 gives Γ_11 = -3.03 for [ethanol, water] and -0.0067 for [water, ethanol]. With x_i restored, both orders give 0.195.
2. henrys_law_constant uses the liquid composition for the vapor (residual_properties.rs:452, from 01f9ace)
molefracs_vapor is computed but not used:
// current
let vapor = State::new(eos, temperature, vle.vapor().density, molefracs.clone())?;
// fix
let vapor = State::new(eos, temperature, vle.vapor().density, molefracs_vapor)?;This is hidden whenever the solvent is pure (x = y). Check: methane in propane/hexane (1:1) at 300 K. v0.10.3 gives H = 15.73 MPa, while lim_{x1→0} y1·p/x1 from a bubble point at x1 = 1e-8 gives 16.88 MPa (-6.8 %). With the fix the two agree to 2e-8.
3. tp_flash acceptance test (crates/feos-core/src/phase_equilibria/tp_flash.rs:445, from 231033e)
molar_gibbs_energy() returns g^L + g^V instead of (1-β)·g^L + β·g^V. It only decides whether an accelerated step is accepted, so converged results are not affected.
Related, but not a regression (tp_flash.rs:374-378)
With iterations = 5, the condition i >= iterations - 3 stores only k_vec[0..=2]. k_vec[3] stays zero, so the extrapolation starts from K = 1. Together with item 3, this means the acceleration never helps. With both fixed, results are identical and the flash needs fewer iterations (tol 1e-12): methanol/hexane at 280 K, 1 bar, 180 → 36; acetone/water at 340 K, 277 → 99.
I found these with an automated code audit (Claude, run in Claude Code) and checked these four by hand, including the numbers above. The audit produced about 40 more candidates across feos-core, feos-dft, the models and py-feos that I haven't verified yet.
How would you like those reported? I'm happy to send PRs for these four (fixes plus regression tests).
Python: reproduce with the released feos 0.10.3
Run from the root of a feos checkout (for the parameter files):
# /// script
# requires-python = ">=3.10"
# dependencies = ["feos==0.10.3"]
# ///
"""Reproduce the regressions from 01f9ace / 231033e with the released feos 0.10.3 wheel.
Run from the root of a feos checkout (for the parameter files):
uv run repro_v0103.py # or: pip install feos==0.10.3 && python repro_v0103.py
"""
import feos
import numpy as np
from si_units import BAR, KELVIN, MEGA, PASCAL
PURE = "parameters/pcsaft/esper2023.json"
BINARY = "parameters/pcsaft/rehner2023_binary.json"
def pcsaft(components):
return feos.EquationOfState.pcsaft(
feos.Parameters.from_json(components, PURE, BINARY, feos.IdentifierOption.Name)
)
print(f"feos {feos.__version__}\n")
# 1. Thermodynamic factor: for a binary, Gamma_11 must not depend on the component order.
print("1. thermodynamic factor, ethanol/water liquid, 300 K, 1 bar")
ew, we = pcsaft(["ethanol", "water"]), pcsaft(["water", "ethanol"])
for x in [0.2, 0.5, 1e-6]:
g1 = feos.State(ew, 300 * KELVIN, pressure=BAR, composition=np.array([x, 1 - x]), density_initialization="liquid").thermodynamic_factor()[0, 0]
g2 = feos.State(we, 300 * KELVIN, pressure=BAR, composition=np.array([1 - x, x]), density_initialization="liquid").thermodynamic_factor()[0, 0]
print(f" x_ethanol = {x:<6g} [ethanol, water]: {g1:10.6f} [water, ethanol]: {g2:10.6f}")
print(" expected: equal in both orders, and -> 1 for x_ethanol -> 0\n")
# 2. Henry's law constant of methane in a mixed solvent vs the infinite-dilution limit of y1 p / x1.
print("2. Henry's law constant, methane in propane/hexane (1:1), 300 K")
eos = pcsaft(["methane", "propane", "hexane"])
t = 300 * KELVIN
h = feos.State.henrys_law_constant(eos, t, np.array([0.0, 0.5, 0.5]))[0]
x1 = 1e-8
vle = feos.PhaseEquilibrium.bubble_point(eos, t, np.array([x1, 0.5 * (1 - x1), 0.5 * (1 - x1)]))
h_lim = vle.vapor.pressure() * vle.vapor.molefracs[0] / x1
print(f" henrys_law_constant: {h / (MEGA * PASCAL):.6f} MPa")
print(f" lim y1*p/x1: {h_lim / (MEGA * PASCAL):.6f} MPa")
print(f" relative deviation: {h / h_lim - 1:.3e} (expected ~0)\n")
# 3./4. tp_flash: number of iterations (results are unaffected).
print("3./4. tp_flash iterations (tol 1e-12)")
for comps, temp, z in [(["methanol", "hexane"], 280.0, [0.5, 0.5]), (["acetone", "water"], 340.0, [0.2, 0.8])]:
eos = pcsaft(comps)
t = temp * KELVIN
z = np.array(z)
if comps[0] == "methanol":
p = BAR
else: # VLE: midway between bubble and dew pressure
p = 0.5 * (feos.PhaseEquilibrium.bubble_point(eos, t, z).vapor.pressure() + feos.PhaseEquilibrium.dew_point(eos, t, z).vapor.pressure())
print(f" {'/'.join(comps)}, {temp} K: ", end="", flush=True)
feos.PhaseEquilibrium.tp_flash(eos, t, p, z, max_iter=200, tol=1e-12, verbosity=feos.Verbosity.Result)
print(" (with both tp_flash fixes: methanol/hexane 36 steps, acetone/water 99 steps)")Rust: tp_flash iterations (after fix)
Place in crates/feos/examples/ and run from the repo root on 0c889e5: cargo run --release --features pcsaft --example tp_flash_iterations. On v0.10.3 the same example gives 180 and 277 steps.
// Number of tp_flash iterations for two associating mixtures (converged results are unaffected).
// Run from the repository root: cargo run --release --features pcsaft --example tp_flash_iterations
use feos::pcsaft::{PcSaft, PcSaftParameters};
use feos_core::parameter::IdentifierOption;
use feos_core::{Contributions, FeosResult, PhaseEquilibrium, SolverOptions, Verbosity};
use nalgebra::dvector;
use quantity::*;
fn pcsaft(components: Vec<&str>) -> FeosResult<PcSaft> {
Ok(PcSaft::new(PcSaftParameters::from_json(
components,
"parameters/pcsaft/esper2023.json",
Some("parameters/pcsaft/rehner2023_binary.json"),
IdentifierOption::Name,
)?))
}
fn main() -> FeosResult<()> {
let options = SolverOptions::new()
.max_iter(200)
.tol(1e-12)
.verbosity(Verbosity::Result);
// methanol/hexane, liquid-liquid, 280 K, 1 bar
let eos = pcsaft(vec!["methanol", "hexane"])?;
print!("methanol/hexane: ");
PhaseEquilibrium::tp_flash(&&eos, 280.0 * KELVIN, BAR, &dvector![0.5, 0.5], None, options, None)?;
// acetone/water, VLE, 340 K, midway between bubble and dew pressure
let eos = pcsaft(vec!["acetone", "water"])?;
let (t, z) = (340.0 * KELVIN, dvector![0.2, 0.8]);
let p_bubble = PhaseEquilibrium::bubble_point(&&eos, t, &z, None, None, Default::default())?
.vapor()
.pressure(Contributions::Total);
let p_dew = PhaseEquilibrium::dew_point(&&eos, t, &z, None, None, Default::default())?
.vapor()
.pressure(Contributions::Total);
print!("acetone/water: ");
PhaseEquilibrium::tp_flash(&&eos, t, 0.5 * (p_bubble + p_dew), &z, None, options, None)?;
Ok(())
}
Activity
prehner commentedon Oct 21, 2026
That is great, thank you! Keep the issues coming, if you encounter more. I hope those 40 candidates are less severe, but in any case we appreciate the feedback. I fixed this in
g-bauer commentedon Oct 21, 2026
This is great - thanks for taking the time to report those issues!