Repository navigation
density_iteration fails for water_3B (na != nb): cancellation in the analytic site fractions #392
Description
After updating to 0.10.2 (thanks for the quick fix in #391!), I ran into a follow-up problem in my water/N2 interfacial tension calculations with the rehner2020 3B water parameters: liquid density iterations and VLEs fail at low temperatures with `density_iteration` did not converge. In 0.10.1 they converged; with the corrected convergence check from #386 the failures become visible.
In my understanding, the cause is the closed form in helmholtz_energy_density_ab_analytic (association/mod.rs, same in association/dft.rs):
let xa = (sqrt + (delta * (rhob - rhoa) + 1.0)).recip() * 2.0;For ρ_A > ρ_B and δ(ρ_A − ρ_B) ≫ 1, sqrt ≈ δ(ρ_A − ρ_B), so the denominator cancels. The relative error of X_A grows like δρ·ε. In liquid water at 20 °C, δρ_B ~ 10³, which gives a noise in p(ρ) of about 1e-5 Pa. That is the size of abstol (1.4e-5 Pa). For na = nb the term vanishes, so 2B and 4C are not affected.
The same happens in mixtures that use the analytic branch, e.g. water 3B + N2 without association sites (tp_flash, Property.bubble_point_pressure, EOS and functional). With an associating second component (Michelsen branch), everything converges.
Reproducer
from decimal import Decimal, getcontext
import feos
import numpy as np
import si_units as si
water = {
"water_2B": feos.PureRecord(feos.Identifier(name="water_2B"), 18.01528, m=1.0, sigma=2.937523956051823,
epsilon_k=272.02757407828676, association_sites=[
{"na": 1.0, "nb": 1.0, "kappa_ab": 0.04448012165716923,
"epsilon_k_ab": 3125.3202766200056}]),
"water_3B": feos.PureRecord(feos.Identifier(name="water_3B"), 18.01528, m=1.6330279759286626,
sigma=2.4570046159550527, epsilon_k=238.31794591948915, association_sites=[
{"na": 2.0, "nb": 1.0, "kappa_ab": 0.037807044304069184,
"epsilon_k_ab": 2749.004567423258}]),
}
# failed density iterations / pure VLEs, 0-125 C, 1-100 bar
for name, record in water.items():
eos = feos.EquationOfState.pcsaft(feos.Parameters.new_pure(record))
n_state = n_vle = 0
for T in np.arange(273.15, 400.0, 2.5):
try:
feos.PhaseEquilibrium.pure(eos, T * si.KELVIN)
except Exception:
n_vle += 1
for p in [1, 5, 20, 40, 80, 100]:
try:
feos.State(eos, T * si.KELVIN, pressure=p * si.BAR, density_initialization="liquid")
except Exception:
n_state += 1
print(f"{name}: State failed {n_state}/306, PhaseEquilibrium.pure failed {n_vle}/51")
# noise of p(rho) in the liquid at 20 C
for name, record in water.items():
eos = feos.EquationOfState.pcsaft(feos.Parameters.new_pure(record))
rho0 = 55.4 * si.KILO * si.MOL / si.METER**3
eps = np.arange(-20, 21) * 2e-16
p = np.array([feos.State(eos, 293.15 * si.KELVIN, density=rho0 * (1 + e)).pressure() / si.PASCAL
for e in eps])
print(f"{name}: noise of p = {np.std(p - np.polyval(np.polyfit(eps, p, 1), eps)):.1e} Pa")
# X_A: current closed form and mass-action variant vs. 50 digits
getcontext().prec = 50
def x_current(d, ra, rb):
s = np.sqrt((d * (ra - rb) + 1.0) ** 2 + d * rb * 4.0)
return 2.0 / (s + (d * (rb - ra) + 1.0)), 2.0 / (s + (d * (ra - rb) + 1.0))
def x_mass_action(d, ra, rb):
xa, xb = x_current(d, ra, rb)
return (1.0 / (1.0 + d * rb * xb), xb) if ra > rb else (xa, 1.0 / (1.0 + d * ra * xa))
def xa_exact(d, ra, rb):
d, ra, rb = Decimal(d), Decimal(ra), Decimal(rb)
return 2 / (((d * (ra - rb) + 1) ** 2 + d * rb * 4).sqrt() + (d * (rb - ra) + 1))
rng = np.random.default_rng(0)
for na, nb in [(1.0, 1.0), (2.0, 1.0)]:
for drho in [1e1, 1e2, 1e3]:
ds = drho * (1 + 1e-3 * rng.random(200))
err = lambda f: max(abs(float((Decimal(f(d, na, nb)[0]) - xa_exact(d, na, nb)) / xa_exact(d, na, nb)))
for d in ds)
print(f"na:nb = {na:.0f}:{nb:.0f}, delta*rho = {drho:4.0f}: "
f"rel. error X_A current {err(x_current):.1e}, mass action {err(x_mass_action):.1e}")Output (0.10.2; with 0.10.1 all failure counts are 0):
water_2B: State failed 0/306, PhaseEquilibrium.pure failed 0/51
water_3B: State failed 5/306, PhaseEquilibrium.pure failed 4/51
water_2B: noise of p = 8.4e-07 Pa
water_3B: noise of p = 1.2e-05 Pa
na:nb = 1:1, delta*rho = 10: rel. error X_A current 1.6e-16, mass action 1.6e-16
na:nb = 1:1, delta*rho = 100: rel. error X_A current 1.5e-16, mass action 1.5e-16
na:nb = 1:1, delta*rho = 1000: rel. error X_A current 2.1e-16, mass action 2.1e-16
na:nb = 2:1, delta*rho = 10: rel. error X_A current 4.6e-16, mass action 1.9e-16
na:nb = 2:1, delta*rho = 100: rel. error X_A current 3.5e-15, mass action 2.4e-16
na:nb = 2:1, delta*rho = 1000: rel. error X_A current 2.2e-14, mass action 1.9e-16
Possible fix
Use the closed form only for the less abundant site and get the other one from mass action, X_A = 1 / (1 + δ ρ_B X_B):
let sqrt = ((delta * (rhoa - rhob) + 1.0).powi(2) + delta * rhob * 4.0).sqrt();
let (xa, xb) = if rhoa.re() > rhob.re() {
let xb = (sqrt + (delta * (rhoa - rhob) + 1.0)).recip() * 2.0;
((delta * rhob * xb + 1.0).recip(), xb)
} else {
let xa = (sqrt + (delta * (rhob - rhoa) + 1.0)).recip() * 2.0;
(xa, (delta * rhoa * xa + 1.0).recip())
};The same change applies elementwise in association/dft.rs.
feos 0.10.2 (PyPI wheel, Windows, Python 3.14)
Activity
prehner commentedon Oct 13, 2026
Good observation, thanks! The solution makes perfect sense, I'll put it in the next patch release.