Robust, fast \(H\),\(S\) flash#

This page describes the algorithm CoolProp uses to solve the enthalpy–entropy flash HmolarSmolar_INPUTS for pure fluids: given molar \((h,s)\), recover the thermodynamic state \((T,\rho)\) (single-phase) or \((T,Q)\) (two-phase).

Historically this was the last input pair without a superancillary-accelerated happy path. The legacy solver seeded a Newton iteration from a random \((T,p)\) guess (HS_flash_generate_TP_singlephase_guess), which is the root cause of the failures reported in issue #2022 and relatives. The current implementation replaces that with a deterministic cascade of homotopy solvers that is 100% robust across the fluid library and roughly 1000× faster on total wall time than the legacy path, while never returning a wrong (metastable) root.

The difficulty#

Single-phase \(H\),\(S\) is a genuine 2-D inversion \((h,s)\mapsto(T,\rho)\) with two obstructions:

  1. The two-phase dome. A straight Newton path from a poor guess can wander into the two-phase region, where there is no single-phase solution.

  2. Apparent non-uniqueness. For some states (notably dense, supercritical hydrogen) two distinct \((T,\rho)\) reproduce the same \((h,s)\) and are both mechanically stable (\(\left.\partial p/\partial\rho\right|_T>0\)). Only one is physical.

Both are resolved below.

Strategy: anchor where a matching-entropy state must exist, then march#

Each leg is a homotopy continuation. The idea: instead of solving a hard problem from a blind guess, start from an easy problem whose answer you already know, then deform it continuously into the hard one, nudging the solution along at each step so you never lose it. Concretely, introduce a parameter that sweeps from 0 to 1: at 0 the problem is trivial (an exactly-known anchor state), at 1 it is the target \((h,s)\). Taking small steps in that parameter and re-converging with a few Newton iterations at each one traces a continuous path from the anchor to the answer – and as long as that path never enters the two-phase dome, every point on it is a valid single-phase state, so the march cannot fall into a non-physical region.

Every leg of the cascade is such a march; they differ only in where the easy anchor lives – each a place where a state of the target entropy is guaranteed to exist:

leg

anchor

owns

  1. saturation anchor

the saturated state whose entropy equals \(s\)

sub-critical liquid/vapor (≈99.9%)

  1. supercritical isentrope

the \(T=T_{\max}\) isotherm at the density matching \(s\)

supercritical targets

  1. ideal-gas departure

the ideal-gas limit (\(\lambda=0\), no dome at all)

universal backstop

The cascade tries them in order and accepts the first result that passes the acceptance gate.

Leg 1 — saturation anchor#

Using the fluid’s superancillaries (Chebyshev fits of the saturation curve, evaluated with no EOS calls), find every saturation temperature \(T_0\) where the saturated entropy equals the target: \(s_{\rm sat}(T_0)=s\). This uses the native \(S\)-inversion get_all_intersections('S', s), which returns all roots and so handles the retrograde (non-monotone) saturated entropy of siloxanes. The anchor is that saturated state, \((T_0,\ \rho_0=\rho_{\rm sat}(T_0))\), on the branch (liquid or vapor) whose entropy matches; ties are broken by enthalpy closeness to the target.

Because the anchor entropy already equals \(s\), the \((h,s)\) chord holds entropy fixed — the path is the isentrope \(s=\text{const}\), and only the enthalpy is driven from the saturated value to the target. When the target entropy lies outside the saturation band (no such saturated state exists), this leg defers to leg 2.

Leg 2 — supercritical isentrope#

For a supercritical target there is no saturated state at its entropy, so we anchor on the dome-free hot isotherm \(T=T_{\max}\) at the density \(\rho_0\) where \(s(T_{\max},\rho_0)=s\) (entropy is monotone in density along an isotherm, so this is a clean 1-D bracketed solve). Marching the isentrope from \(T_{\max}\) downward to a supercritical target stays at \(T>T_c\) the whole way and so never touches the dome — structurally guaranteeing a valid path.

Leg 3 — ideal-gas departure homotopy#

As a universal backstop, scale the residual (real-fluid) Helmholtz energy by a parameter \(\lambda\in[0,1]\):

\[\alpha(\tau,\delta;\lambda) = \alpha^0(\tau,\delta) + \lambda\,\alpha^{\mathrm r}(\tau,\delta).\]

At \(\lambda=0\) the fluid is an ideal gas — it has no two-phase dome, and the \((h,s)\) inversion is two trivial monotone 1-D rootfinds (for \(T\) then \(\rho\)). Continuing \(\lambda:0\to1\) grows the real fluid (and its dome) from nothing while the corrector tracks the stable branch. All properties along the way come from the cached \(\alpha^0\) / \(\alpha^{\mathrm r}\) derivatives, so this leg needs no superancillary.

The corrector: Newton in \((T,\ w=\ln\rho)\)#

Working in \(w=\ln\rho\) makes \(\rho>0\) automatic and conditions the entropy direction. At each homotopy step the corrector solves \(h(T,\rho)=h_\lambda,\ s(T,\rho)=s_\lambda\) by Newton, with the Jacobian columns \(\partial/\partial T\) and \(\partial/\partial w=\rho\,\partial/\partial\rho\). Two guards keep it physical:

  • Mechanical-stability damping — each step is halved until \(\left.\partial p/\partial\rho\right|_T>0\), keeping the corrector on the stable branch.

  • Domain clamp with slack — \(T\) is kept in \([\,0.98\,T_{\min},\ 1.5\,T_{\max}]\). The 2% slack below \(T_{\min}\) matters near the triple point: the compressed-liquid isentrope folds a few hundredths of a percent below \(T_{\min}\) (water’s density anomaly) before climbing back to an in-range target, and a hard clamp would stall there.

Failed steps trigger adaptive subdivision (double the number of \(\lambda\) steps).

Acceptance gate and the uniqueness theorem#

A candidate \((T,\rho)\) is accepted only if it

  1. reproduces \((h,s)\) to tolerance,

  2. lies in \([T_{\min},T_{\max}]\),

  3. is fully intrinsically stable: \(\left.\partial p/\partial\rho\right|_T>0\) (mechanical) and \(c_v>0\) (thermal), and

  4. lies outside the two-phase dome.

Conditions 3–4 are what make the answer unique and physical. At fixed entropy, \(\mathrm{d}h = v\,\mathrm{d}p\), and along the adiabat \(\left.\partial p/\partial v\right|_s<0\) for any stable state, so \(h\) is strictly monotone in \(v\) at fixed \(s\): the map \((h,s)\to(T,\rho)\) is injective over the fully stable region. The apparent non-uniqueness (e.g. the dense supercritical-hydrogen pair) is not a real feature — the spurious second root has \(c_v<0\) (thermally unstable, an EOS extrapolation artifact), so the \(c_v>0\) gate removes it. Because the gate certifies the unique physical root, the solver never needs to know the true answer in advance — a matching, stable, outside-dome root is the answer.

For a two-phase \((h,s)\) there is no stable single-phase root: every cascade leg lands on an in-dome metastable state and is rejected by condition 4. The cascade therefore fails cleanly, and the dispatcher falls through to the EOS-exact two-phase solve (the \(Q_h=Q_s\) saturation problem, the \(H\),\(S\) analogue of the \(D\)+\(X\) happy path).

Try it#

The happy path is engaged automatically by PropsSI / the low-level interface for pure fluids when superancillaries are enabled. Below we round-trip a single-phase state and a two-phase state, then time the flash with the low-level interface – a reused AbstractState calling update(HmolarSmolar_INPUTS, ...) – so the measurement reflects the flash cost itself rather than the per-call string-parsing/factory overhead of PropsSI. The superancillary path is toggled with the ENABLE_SUPERANCILLARIES config flag.

[1]:

import time
import CoolProp.CoolProp as CP

fluid = 'Water'
# A single-phase reference state at (p, T); read off its (h, s)
p, T = 50e6, 700.0  # supercritical water
h = CP.PropsSI('Hmolar', 'P', p, 'T', T, fluid)
s = CP.PropsSI('Smolar', 'P', p, 'T', T, fluid)

# Recover (T, rho) from (h, s)
T_rt   = CP.PropsSI('T',      'Hmolar', h, 'Smolar', s, fluid)
rho_rt = CP.PropsSI('Dmolar', 'Hmolar', h, 'Smolar', s, fluid)
rho_ref = CP.PropsSI('Dmolar', 'P', p, 'T', T, fluid)
print(f'single-phase: T {T:.4f} -> {T_rt:.4f} K, rho {rho_ref:.3f} -> {rho_rt:.3f} mol/m^3')

single-phase: T 700.0000 -> 700.0000 K, rho 27256.484 -> 27256.484 mol/m^3

[2]:

# Two-phase round-trip: take a (T, Q) state, read (h, s), recover (T, Q)
Tsat, Q = 450.0, 0.3
h2 = CP.PropsSI('Hmolar', 'T', Tsat, 'Q', Q, fluid)
s2 = CP.PropsSI('Smolar', 'T', Tsat, 'Q', Q, fluid)
T2 = CP.PropsSI('T', 'Hmolar', h2, 'Smolar', s2, fluid)
Q2 = CP.PropsSI('Q', 'Hmolar', h2, 'Smolar', s2, fluid)
print(f'two-phase:    T {Tsat:.4f} -> {T2:.4f} K, Q {Q:.4f} -> {Q2:.4f}')

two-phase:    T 450.0000 -> 450.0000 K, Q 0.3000 -> 0.3000

[3]:

# Low-level timing: reuse one AbstractState and call update() directly, so the
# per-call PropsSI overhead (string parsing, backend lookup) is excluded.
state = CP.AbstractState('HEOS', fluid)
HS = CP.HmolarSmolar_INPUTS
state.update(HS, h, s)  # warm up

def time_hs(n=20000):
    t0 = time.perf_counter()
    for _ in range(n):
        state.update(HS, h, s)
    return 1e6 * (time.perf_counter() - t0) / n  # us per flash

CP.set_config_bool(CP.ENABLE_SUPERANCILLARIES, True)
us_fast = time_hs()
CP.set_config_bool(CP.ENABLE_SUPERANCILLARIES, False)
us_legacy = time_hs()
CP.set_config_bool(CP.ENABLE_SUPERANCILLARIES, True)  # restore
print(f'happy path: {us_fast:6.1f} us/flash   legacy: {us_legacy:8.1f} us/flash   speedup: {us_legacy/us_fast:5.1f}x')

happy path:   86.1 us/flash   legacy:   1868.0 us/flash   speedup:  21.7x

Performance map#

The figure below (generated by CoolProp.Plots.HSFlashTimingMap from a dense \((\log p, T)\) sweep) shows the cost of the cascade over the entire single-phase domain of water, in \(p\)–\(T\) and \(h\)–\(s\) coordinates: EOS-evaluation count, solve time, and which leg solved each point. The body of the domain is the saturation-anchor leg at ~6–20 evaluations; a thin band just above the critical point is the supercritical-isentrope leg. Across 130 pure fluids and ~190k single-phase points the cascade solves 100% with 0 regressions against the legacy path, at a mean of ~16 EOS evaluations and a worst case under 40 (vs the legacy heavy tail of ~0.5 s on the hardest near-critical points).

./build_rel/CatchTestRunner "[HS_watermap]"          # writes water_hs_map.csv
python -m CoolProp.Plots.HSFlashTimingMap water_hs_map.csv   # --coords pt,hs,both

References & implementation#

  • Implementation: FlashRoutines::HS_flash in src/Backends/Helmholtz/FlashRoutines.cpp (the hs_corrector / hs_leg_* / hs_cascade helpers), mirroring the HSU_D_flash hybrid.

  • Characterization & stress tests: src/Tests/CoolProp-Tests-HS.cpp ([HS_grid], [HS_casc], [HS_stress], [HS_prod2ph]).

  • Superancillaries: see the Superancillaries page.

  • Original issue: #2022.