{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Robust, fast $H$,$S$ flash\n", "\n", "This page describes the algorithm CoolProp uses to solve the **enthalpy–entropy** flash\n", "`HmolarSmolar_INPUTS` for pure fluids: given molar $(h,s)$, recover the thermodynamic state\n", "$(T,\\rho)$ (single-phase) or $(T,Q)$ (two-phase).\n", "\n", "Historically this was the last input pair without a superancillary-accelerated *happy path*.\n", "The legacy solver seeded a Newton iteration from a **random** $(T,p)$ guess\n", "(`HS_flash_generate_TP_singlephase_guess`), which is the root cause of the failures reported in\n", "[issue #2022](https://github.com/CoolProp/CoolProp/issues/2022) and relatives. The current\n", "implementation replaces that with a deterministic *cascade* of homotopy solvers that is\n", "**100% robust** across the fluid library and roughly **1000× faster on total wall time** than the\n", "legacy path, while never returning a wrong (metastable) root." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## The difficulty\n", "\n", "Single-phase $H$,$S$ is a genuine 2-D inversion $(h,s)\\mapsto(T,\\rho)$ with two obstructions:\n", "\n", "1. **The two-phase dome.** A straight Newton path from a poor guess can wander into the\n", " two-phase region, where there is no single-phase solution.\n", "2. **Apparent non-uniqueness.** For some states (notably dense, supercritical hydrogen) *two*\n", " distinct $(T,\\rho)$ reproduce the same $(h,s)$ and are both mechanically stable\n", " ($\\left.\\partial p/\\partial\\rho\\right|_T>0$). Only one is physical.\n", "\n", "Both are resolved below." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Strategy: anchor where a matching-entropy state must exist, then march\n", "\n", "Each leg is a **homotopy continuation**. The idea: instead of solving a hard problem from a\n", "blind guess, start from an *easy* problem whose answer you already know, then deform it\n", "continuously into the hard one, nudging the solution along at each step so you never lose it.\n", "Concretely, introduce a parameter that sweeps from 0 to 1: at 0 the problem is trivial (an\n", "exactly-known anchor state), at 1 it is the target $(h,s)$. Taking small steps in that\n", "parameter and re-converging with a few Newton iterations at each one traces a continuous path\n", "from the anchor to the answer -- and as long as that path never enters the two-phase dome, every\n", "point on it is a valid single-phase state, so the march cannot fall into a non-physical region.\n", "\n", "Every leg of the cascade is such a march; they differ only in *where* the easy anchor lives --\n", "each a place where a state of the **target entropy** is guaranteed to exist:\n", "\n", "| leg | anchor | owns |\n", "|----|--------|------|\n", "| 1. saturation anchor | the saturated state whose entropy equals $s$ | sub-critical liquid/vapor (≈99.9%) |\n", "| 2. supercritical isentrope | the $T=T_{\\max}$ isotherm at the density matching $s$ | supercritical targets |\n", "| 3. ideal-gas departure | the ideal-gas limit ($\\lambda=0$, no dome at all) | universal backstop |\n", "\n", "The cascade tries them in order and accepts the first result that passes the acceptance gate.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Leg 1 — saturation anchor\n", "\n", "Using the fluid's **superancillaries** (Chebyshev fits of the saturation curve, evaluated with no\n", "EOS calls), find every saturation temperature $T_0$ where the saturated entropy equals the target:\n", "$s_{\\rm sat}(T_0)=s$. This uses the native $S$-inversion `get_all_intersections('S', s)`, which\n", "returns *all* roots and so handles the retrograde (non-monotone) saturated entropy of siloxanes.\n", "The anchor is that saturated state, $(T_0,\\ \\rho_0=\\rho_{\\rm sat}(T_0))$, on the branch (liquid or\n", "vapor) whose entropy matches; ties are broken by enthalpy closeness to the target.\n", "\n", "Because the anchor entropy already equals $s$, the $(h,s)$ chord holds entropy fixed — the path is\n", "the **isentrope** $s=\\text{const}$, and only the enthalpy is driven from the saturated value to the\n", "target. When the target entropy lies outside the saturation band (no such saturated state exists),\n", "this leg defers to leg 2." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Leg 2 — supercritical isentrope\n", "\n", "For a supercritical target there is no saturated state at its entropy, so we anchor on the\n", "**dome-free** hot isotherm $T=T_{\\max}$ at the density $\\rho_0$ where $s(T_{\\max},\\rho_0)=s$\n", "(entropy is monotone in density along an isotherm, so this is a clean 1-D bracketed solve).\n", "Marching the isentrope from $T_{\\max}$ downward to a supercritical target stays at $T>T_c$ the\n", "whole way and so **never touches the dome** — structurally guaranteeing a valid path." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Leg 3 — ideal-gas departure homotopy\n", "\n", "As a universal backstop, scale the *residual* (real-fluid) Helmholtz energy by a parameter\n", "$\\lambda\\in[0,1]$:\n", "\n", "$$\\alpha(\\tau,\\delta;\\lambda) = \\alpha^0(\\tau,\\delta) + \\lambda\\,\\alpha^{\\mathrm r}(\\tau,\\delta).$$\n", "\n", "At $\\lambda=0$ the fluid is an **ideal gas** — it has *no* two-phase dome, and the $(h,s)$ inversion\n", "is two trivial monotone 1-D rootfinds (for $T$ then $\\rho$). Continuing $\\lambda:0\\to1$ grows the\n", "real fluid (and its dome) from nothing while the corrector tracks the stable branch. All\n", "properties along the way come from the cached $\\alpha^0$ / $\\alpha^{\\mathrm r}$ derivatives, so this\n", "leg needs no superancillary." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## The corrector: Newton in $(T,\\ w=\\ln\\rho)$\n", "\n", "Working in $w=\\ln\\rho$ makes $\\rho>0$ automatic and conditions the entropy direction. At each\n", "homotopy step the corrector solves $h(T,\\rho)=h_\\lambda,\\ s(T,\\rho)=s_\\lambda$ by Newton, with the\n", "Jacobian columns $\\partial/\\partial T$ and $\\partial/\\partial w=\\rho\\,\\partial/\\partial\\rho$.\n", "Two guards keep it physical:\n", "\n", "* **Mechanical-stability damping** — each step is halved until $\\left.\\partial p/\\partial\\rho\\right|_T>0$,\n", " keeping the corrector on the stable branch.\n", "* **Domain clamp with slack** — $T$ is kept in $[\\,0.98\\,T_{\\min},\\ 1.5\\,T_{\\max}]$. The 2% slack\n", " below $T_{\\min}$ matters near the triple point: the compressed-liquid isentrope *folds* a few\n", " hundredths of a percent below $T_{\\min}$ (water's density anomaly) before climbing back to an\n", " in-range target, and a hard clamp would stall there.\n", "\n", "Failed steps trigger adaptive subdivision (double the number of $\\lambda$ steps)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Acceptance gate and the uniqueness theorem\n", "\n", "A candidate $(T,\\rho)$ is accepted only if it\n", "\n", "1. reproduces $(h,s)$ to tolerance,\n", "2. lies in $[T_{\\min},T_{\\max}]$,\n", "3. is **fully intrinsically stable**: $\\left.\\partial p/\\partial\\rho\\right|_T>0$ (mechanical) **and** $c_v>0$ (thermal), and\n", "4. lies **outside** the two-phase dome.\n", "\n", "Conditions 3–4 are what make the answer unique and physical. At fixed entropy,\n", "$\\mathrm{d}h = v\\,\\mathrm{d}p$, and along the *adiabat* $\\left.\\partial p/\\partial v\\right|_s<0$ for any\n", "stable state, so $h$ is strictly monotone in $v$ at fixed $s$: the map $(h,s)\\to(T,\\rho)$ is\n", "**injective over the fully stable region**. The apparent non-uniqueness (e.g. the dense\n", "supercritical-hydrogen pair) is *not* a real feature — the spurious second root has $c_v<0$\n", "(thermally unstable, an EOS extrapolation artifact), so the $c_v>0$ gate removes it. Because the\n", "gate certifies the unique physical root, **the solver never needs to know the true answer in\n", "advance** — a matching, stable, outside-dome root *is* the answer.\n", "\n", "For a **two-phase** $(h,s)$ there is no stable single-phase root: every cascade leg lands on an\n", "in-dome metastable state and is rejected by condition 4. The cascade therefore *fails cleanly*,\n", "and the dispatcher falls through to the EOS-exact two-phase solve (the $Q_h=Q_s$ saturation\n", "problem, the $H$,$S$ analogue of the $D$+$X$ happy path)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Try it\n", "\n", "The happy path is engaged automatically by `PropsSI` / the low-level interface for pure\n", "fluids when superancillaries are enabled. Below we round-trip a single-phase state and a\n", "two-phase state, then time the flash with the **low-level interface** -- a reused\n", "`AbstractState` calling `update(HmolarSmolar_INPUTS, ...)` -- so the measurement reflects the\n", "flash cost itself rather than the per-call string-parsing/factory overhead of `PropsSI`.\n", "The superancillary path is toggled with the `ENABLE_SUPERANCILLARIES` config flag.\n" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "execution": { "iopub.execute_input": "2026-05-29T01:48:37.750090Z", "iopub.status.busy": "2026-05-29T01:48:37.749587Z", "iopub.status.idle": "2026-05-29T01:48:38.793742Z", "shell.execute_reply": "2026-05-29T01:48:38.793297Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "single-phase: T 700.0000 -> 700.0000 K, rho 27256.484 -> 27256.484 mol/m^3\n" ] } ], "source": [ "import time\n", "import CoolProp.CoolProp as CP\n", "\n", "fluid = 'Water'\n", "# A single-phase reference state at (p, T); read off its (h, s)\n", "p, T = 50e6, 700.0 # supercritical water\n", "h = CP.PropsSI('Hmolar', 'P', p, 'T', T, fluid)\n", "s = CP.PropsSI('Smolar', 'P', p, 'T', T, fluid)\n", "\n", "# Recover (T, rho) from (h, s)\n", "T_rt = CP.PropsSI('T', 'Hmolar', h, 'Smolar', s, fluid)\n", "rho_rt = CP.PropsSI('Dmolar', 'Hmolar', h, 'Smolar', s, fluid)\n", "rho_ref = CP.PropsSI('Dmolar', 'P', p, 'T', T, fluid)\n", "print(f'single-phase: T {T:.4f} -> {T_rt:.4f} K, rho {rho_ref:.3f} -> {rho_rt:.3f} mol/m^3')" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "execution": { "iopub.execute_input": "2026-05-29T01:48:38.821876Z", "iopub.status.busy": "2026-05-29T01:48:38.821640Z", "iopub.status.idle": "2026-05-29T01:48:38.826078Z", "shell.execute_reply": "2026-05-29T01:48:38.825644Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "two-phase: T 450.0000 -> 450.0000 K, Q 0.3000 -> 0.3000\n" ] } ], "source": [ "# Two-phase round-trip: take a (T, Q) state, read (h, s), recover (T, Q)\n", "Tsat, Q = 450.0, 0.3\n", "h2 = CP.PropsSI('Hmolar', 'T', Tsat, 'Q', Q, fluid)\n", "s2 = CP.PropsSI('Smolar', 'T', Tsat, 'Q', Q, fluid)\n", "T2 = CP.PropsSI('T', 'Hmolar', h2, 'Smolar', s2, fluid)\n", "Q2 = CP.PropsSI('Q', 'Hmolar', h2, 'Smolar', s2, fluid)\n", "print(f'two-phase: T {Tsat:.4f} -> {T2:.4f} K, Q {Q:.4f} -> {Q2:.4f}')" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "execution": { "iopub.execute_input": "2026-05-29T01:48:38.827606Z", "iopub.status.busy": "2026-05-29T01:48:38.827445Z", "iopub.status.idle": "2026-05-29T01:49:17.914922Z", "shell.execute_reply": "2026-05-29T01:49:17.914468Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "happy path: 86.1 us/flash legacy: 1868.0 us/flash speedup: 21.7x\n" ] } ], "source": [ "# Low-level timing: reuse one AbstractState and call update() directly, so the\n", "# per-call PropsSI overhead (string parsing, backend lookup) is excluded.\n", "state = CP.AbstractState('HEOS', fluid)\n", "HS = CP.HmolarSmolar_INPUTS\n", "state.update(HS, h, s) # warm up\n", "\n", "def time_hs(n=20000):\n", " t0 = time.perf_counter()\n", " for _ in range(n):\n", " state.update(HS, h, s)\n", " return 1e6 * (time.perf_counter() - t0) / n # us per flash\n", "\n", "CP.set_config_bool(CP.ENABLE_SUPERANCILLARIES, True)\n", "us_fast = time_hs()\n", "CP.set_config_bool(CP.ENABLE_SUPERANCILLARIES, False)\n", "us_legacy = time_hs()\n", "CP.set_config_bool(CP.ENABLE_SUPERANCILLARIES, True) # restore\n", "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')\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Performance map\n", "\n", "The figure below (generated by `CoolProp.Plots.HSFlashTimingMap` from a dense $(\\log p, T)$ sweep)\n", "shows the cost of the cascade over the entire single-phase domain of water, in $p$–$T$ and $h$–$s$\n", "coordinates: EOS-evaluation count, solve time, and which leg solved each point. The body of the\n", "domain is the saturation-anchor leg at ~6–20 evaluations; a thin band just above the critical\n", "point is the supercritical-isentrope leg. Across 130 pure fluids and ~190k single-phase\n", "points the cascade solves **100%** with **0 regressions** against the legacy path, at a mean of\n", "~16 EOS evaluations and a worst case under 40 (vs the legacy heavy tail of ~0.5 s on the\n", "hardest near-critical points).\n", "\n", "```\n", "./build_rel/CatchTestRunner \"[HS_watermap]\" # writes water_hs_map.csv\n", "python -m CoolProp.Plots.HSFlashTimingMap water_hs_map.csv # --coords pt,hs,both\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## References & implementation\n", "\n", "* Implementation: `FlashRoutines::HS_flash` in `src/Backends/Helmholtz/FlashRoutines.cpp`\n", " (the `hs_corrector` / `hs_leg_*` / `hs_cascade` helpers), mirroring the `HSU_D_flash` hybrid.\n", "* Characterization & stress tests: `src/Tests/CoolProp-Tests-HS.cpp` (`[HS_grid]`, `[HS_casc]`,\n", " `[HS_stress]`, `[HS_prod2ph]`).\n", "* Superancillaries: see the [Superancillaries](SuperAncillary.ipynb) page.\n", "* Original issue: [#2022](https://github.com/CoolProp/CoolProp/issues/2022)." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.13.3" } }, "nbformat": 4, "nbformat_minor": 4 }