SVDSBTL backend — validation against HEOS

This page compares the SVDSBTL surrogate backend against the HEOS Helmholtz-energy reference across four input pairs. For (P, T), (H, P), and (P, S) each panel shows the density error log10 | (rho_SVDSBTL - rho_HEOS) / rho_HEOS |. For (rho, T) density is the input (exact by construction), so the DT-indexed surface is checked on pressure instead: log10 | (p_SVDSBTL - p_HEOS) / p_HEOS | over a (T, rho) rectangle.

The heatmaps are interactive — hover for the local error value, drag to pan, scroll to zoom, double-click to reset. The figures are static HTML (no kernel, no Binder, no JupyterLite); the interactivity is plotly.js embedded by nbsphinx.

Backends: HEOS (reference) vs SVDSBTL&HEOS (rank-truncated SVD surrogate over HEOS).

Inputs: PT_INPUTS, HmassP_INPUTS, PSmass_INPUTS (validate rhomass()) and DmassT_INPUTS (validates p()).

Fluids: Water, Argon, Hydrogen, R1234yf, R245fa, D6, CO2, Helium.

import os
from pathlib import Path

# Route SVDSBTL surface I/O to a docs-dedicated cache directory so the
# build does not collide with developer interactive work in
# ~/.CoolProp/SVDTables/. Setting this *before* importing CoolProp lets
# the C++ config-init read the env var on first access, and joblib's
# loky workers inherit the env so they hit the same cache.
_cache = os.environ.setdefault(
    'COOLPROP_ALTERNATIVE_SVDTABLES_DIRECTORY',
    str(Path.home() / '.CoolProp' / 'SVDTables-docs'),
)
Path(_cache).mkdir(parents=True, exist_ok=True)

import numpy as np
import plotly.io as pio
import plotly.graph_objects as go
import CoolProp.CoolProp as CP
from CoolProp.CoolProp import PT_INPUTS, HmassP_INPUTS, QT_INPUTS, DmassT_INPUTS, PSmass_INPUTS

# Also set via the explicit config setter in the parent process. NOTE:
# this only mutates the parent's Configuration singleton; joblib's loky
# workers spawn fresh interpreters with their own singletons and rely
# entirely on env-var inheritance (which is why os.environ.setdefault
# above MUST stay before the CoolProp import — do not remove it).
CP.set_config_string(CP.ALTERNATIVE_SVDTABLES_DIRECTORY, _cache)

pio.renderers.default = 'notebook_connected'
print(f'SVDSBTL cache directory: {_cache}')

FLUIDS = ['Water', 'Argon', 'Hydrogen', 'R1234yf', 'R245fa', 'D6', 'CO2', 'Helium']
INPUT_PAIRS = ['PT', 'HP', 'DT', 'PS']
N_X, N_Y = 300, 300    # grid resolution per axis (90k cells per panel)
BAND_REL = 0.02            # half-width band around saturation line to mask (PT only)

SVDSBTL cache directory: /__w/CoolProp/CoolProp/.svdtables-cache

def _hp_axes(ref, Tmin, Tmax, pmin, pmax, N_Y):
    """Estimate an enthalpy range covering the single-phase region for
    the HP plot by sampling a coarse (T, p) grid across the rectangle.

    Sampling the interior — not just the four corners — matters for fluids
    like CO2 whose high triple-point pressure and steep melting line put
    most corners outside HEOS validity: 3 of CO2's 4 corners throw (p below
    p_triple, or T below T_melt(p)).  A corners-only range would then
    collapse to a single enthalpy and render a degenerate sliver of a panel.
    """
    hs = []
    for T_q in np.linspace(Tmin, Tmax, 8):
        for p_q in np.geomspace(pmin, pmax, 8):
            try:
                ref.update(PT_INPUTS, p_q, T_q)
                hs.append(ref.hmass())
            except Exception:
                pass
    if len(hs) < 2:
        # Almost nothing valid in the rectangle even after sampling.
        # Return a benign ~1 J/kg-wide axis so the panel still renders;
        # the surface queries there will mostly be NaN.
        h0 = hs[0] if hs else 0.0
        return np.linspace(h0 - 1.0, h0 + 1.0, N_Y)
    return np.linspace(min(hs), max(hs), N_Y)

def _ps_axes(ref, Tmin, Tmax, pmin, pmax, N_Y):
    """Estimate a specific-entropy range covering the single-phase region
    for the PS plot by sampling a coarse (T, p) grid across the rectangle.
    Entropy twin of _hp_axes — entropy has no natural reduced form, so the
    axis stays in absolute J/(kg.K).  Interior sampling matters for the
    same reason as _hp_axes (CO2's corners fall outside HEOS validity).
    """
    ss = []
    for T_q in np.linspace(Tmin, Tmax, 8):
        for p_q in np.geomspace(pmin, pmax, 8):
            try:
                ref.update(PT_INPUTS, p_q, T_q)
                ss.append(ref.smass())
            except Exception:
                pass
    if len(ss) < 2:
        s0 = ss[0] if ss else 0.0
        return np.linspace(s0 - 1.0, s0 + 1.0, N_Y)
    return np.linspace(min(ss), max(ss), N_Y)

def _dt_axes(ref, Tmin, Tmax, pmin, pmax, N_X):
    """Density axis for the DT panel.  rho is the *input* for DmassT, so
    we need a rho range, not a p range: sample the densities reached across
    the (T, p) rectangle and span their min..max on a log axis.  Sampling
    the interior (not just corners) keeps fluids like CO2 — whose corners
    fall below p_triple / the melting line — from collapsing it.  The
    spanned range also covers the two-phase dome, whose densities lie
    between the vapor and liquid branches sampled here.
    """
    rhos = []
    for T_q in np.linspace(Tmin, Tmax, 8):
        for p_q in np.geomspace(pmin, pmax, 8):
            try:
                ref.update(PT_INPUTS, p_q, T_q)
                rhos.append(ref.rhomass())
            except Exception:
                pass
    if len(rhos) < 2:
        r0 = rhos[0] if rhos else 1.0
        return np.geomspace(max(r0 * 0.5, 1e-6), r0 * 2.0, N_X)
    return np.geomspace(max(min(rhos), 1e-6), max(rhos), N_X)

def compute_error_grid(fluid, input_pair, N_X=N_X, N_Y=N_Y):
    """Return (x_axis, y_axis, log10_err, Tc, pc, x_label, y_label, err_msg).

    input_pair is 'PT', 'HP', 'DT', or 'PS'.
      'PT' -> x = p (log, Pa),        y = T (linear, K).   Validates rho.
             Mask a +/- BAND_REL band around the saturation line.
      'HP' -> x = p (log, Pa),        y = h (linear, J/kg). Validates rho.
             No mask; the SVDSBTL surface is undefined inside the dome and
             queries there are caught and set to NaN.
      'PS' -> x = p (log, Pa),        y = s (linear, J/kg/K). Validates rho.
             No mask (same dome handling as HP — the two-phase wedge is a
             true 2D region in (s, p) and surface queries there are NaN).
      'DT' -> x = rho (log, kg/m^3),  y = T (linear, K).   Validates p.
             rho is the input (exact), so the DT-indexed surface is checked
             on p(rho, T).  No mask: the two-phase dome is a valid DT region
             (SVDSBTL routes it via the dome lever-rule, matching p_sat(T)).

    SVDSBTL surface build failures are caught and reported via err_msg;
    the returned grid is all-NaN in that case.
    """
    ref = CP.AbstractState('HEOS', fluid)
    Tc = ref.T_critical()
    pc = ref.p_critical()
    Tmin = max(ref.Tmin(), 0.6 * Tc)
    Tmax = min(ref.Tmax(), 1.4 * Tc)
    pmin = max(1e3, 0.01 * pc)
    pmax = 3.0 * pc

    psat = None
    if input_pair == 'PT':
        x_axis = np.geomspace(pmin, pmax, N_X)
        y_axis = np.linspace(Tmin, Tmax, N_Y)
        x_label = 'reduced pressure p / p_c'
        y_label = 'reduced temperature T / T_c'
        # Saturation line p_sat(T) for masking.
        psat = np.full_like(y_axis, np.nan)
        for i, T in enumerate(y_axis):
            if T < Tc:
                try:
                    ref.update(QT_INPUTS, 0.0, T)
                    psat[i] = ref.p()
                except Exception:
                    pass
    elif input_pair == 'HP':
        x_axis = np.geomspace(pmin, pmax, N_X)
        y_axis = _hp_axes(ref, Tmin, Tmax, pmin, pmax, N_Y)
        x_label = 'reduced pressure p / p_c'
        y_label = 'specific enthalpy h [kJ/kg]'
    elif input_pair == 'DT':
        x_axis = _dt_axes(ref, Tmin, Tmax, pmin, pmax, N_X)
        y_axis = np.linspace(Tmin, Tmax, N_Y)
        x_label = 'density rho [kg/m^3]'
        y_label = 'reduced temperature T / T_c'
    elif input_pair == 'PS':
        x_axis = np.geomspace(pmin, pmax, N_X)
        y_axis = _ps_axes(ref, Tmin, Tmax, pmin, pmax, N_Y)
        x_label = 'reduced pressure p / p_c'
        y_label = 'specific entropy s [J/(kg K)]'
    else:
        raise ValueError(f'unknown input_pair: {input_pair!r}')

    err = np.full((N_Y, N_X), np.nan)

    try:
        svd = CP.AbstractState('SVDSBTL&HEOS', fluid)
    except Exception as e:
        return (x_axis, y_axis, err, Tc, pc, x_label, y_label,
                f'SVDSBTL surface build failed: {e}')

    for i, y in enumerate(y_axis):
        for j, xv in enumerate(x_axis):
            if psat is not None and np.isfinite(psat[i]) and abs(xv - psat[i]) / psat[i] < BAND_REL:
                continue
            try:
                if input_pair == 'PT':
                    ref.update(PT_INPUTS, xv, y)
                    svd.update(PT_INPUTS, xv, y)
                    q_ref, q_svd = ref.rhomass(), svd.rhomass()
                elif input_pair == 'HP':
                    ref.update(HmassP_INPUTS, y, xv)
                    svd.update(HmassP_INPUTS, y, xv)
                    q_ref, q_svd = ref.rhomass(), svd.rhomass()
                elif input_pair == 'PS':
                    # PSmass takes (p, s) = (xv, y); validate density.
                    ref.update(PSmass_INPUTS, xv, y)
                    svd.update(PSmass_INPUTS, xv, y)
                    q_ref, q_svd = ref.rhomass(), svd.rhomass()
                else:  # DT — rho is the input; validate pressure
                    ref.update(DmassT_INPUTS, xv, y)
                    svd.update(DmassT_INPUTS, xv, y)
                    q_ref, q_svd = ref.p(), svd.p()
                if q_ref > 0:
                    err[i, j] = abs(q_svd - q_ref) / q_ref
            except Exception:
                pass
    # Floor zero / near-zero errors at 1e-16 so exact-match cells render at
    # the colorbar floor (best-possible accuracy) instead of going to -inf
    # and showing as NaN/white. Inside the two-phase dome SVDSBTL returns the
    # same value as HEOS (shared saturation lookup), so err is ~0 there.
    with np.errstate(divide='ignore', invalid='ignore'):
        log_err = np.log10(np.maximum(err, 1e-16))
    return x_axis, y_axis, log_err, Tc, pc, x_label, y_label, None

def plot_from_grid(fluid, input_pair, grid):
    x_axis, y_axis, log_err, Tc, pc, x_label, y_label, err_msg = grid
    title = f'{fluid} [{input_pair}] \u2014 Tc={Tc:.2f} K, pc={pc/1e5:.3f} bar'
    if err_msg is not None:
        title += '  [SVDSBTL surface unavailable]'
    if input_pair == 'PT':
        # x = p/p_c (log), y = T/T_c (linear). z is stored as [y_index, x_index].
        heat = go.Heatmap(
            x=x_axis / pc, y=y_axis / Tc, z=log_err,
            colorscale='Viridis', zmin=-12, zmax=-2,
            colorbar=dict(title='log10 |\u0394\u03c1/\u03c1|'),
            hovertemplate=('p/p_c=%{x:.3f}<br>'
                           'T/T_c=%{y:.3f}<br>'
                           'log10 err=%{z:.2f}<extra></extra>'),
        )
        xaxis = dict(title=x_label, type='log')
        yaxis = dict(title=y_label)
    elif input_pair == 'HP':  # Mollier convention: x = h (linear), y = p (log).
        heat = go.Heatmap(
            x=y_axis / 1e3, y=x_axis / pc, z=log_err.T,
            colorscale='Viridis', zmin=-12, zmax=-2,
            colorbar=dict(title='log10 |\u0394\u03c1/\u03c1|'),
            hovertemplate=('h=%{x:.1f} kJ/kg<br>'
                           'p/p_c=%{y:.3f}<br>'
                           'log10 err=%{z:.2f}<extra></extra>'),
        )
        xaxis = dict(title=y_label)                     # 'specific enthalpy h [kJ/kg]'
        yaxis = dict(title='reduced pressure p / p_c', type='log')
    elif input_pair == 'PS':  # Mollier-like: x = s (linear), y = p (log).
        heat = go.Heatmap(
            x=y_axis, y=x_axis / pc, z=log_err.T,
            colorscale='Viridis', zmin=-12, zmax=-2,
            colorbar=dict(title='log10 |\u0394\u03c1/\u03c1|'),
            hovertemplate=('s=%{x:.1f} J/(kg K)<br>'
                           'p/p_c=%{y:.3f}<br>'
                           'log10 err=%{z:.2f}<extra></extra>'),
        )
        xaxis = dict(title=y_label)                     # 'specific entropy s [J/(kg K)]'
        yaxis = dict(title='reduced pressure p / p_c', type='log')
    else:  # DT — x = rho (log), y = T/T_c (linear). Validates pressure.
        heat = go.Heatmap(
            x=x_axis, y=y_axis / Tc, z=log_err,
            colorscale='Viridis', zmin=-12, zmax=-2,
            colorbar=dict(title='log10 |\u0394p/p|'),
            hovertemplate=('rho=%{x:.4g} kg/m^3<br>'
                           'T/T_c=%{y:.3f}<br>'
                           'log10 err=%{z:.2f}<extra></extra>'),
        )
        xaxis = dict(title=x_label, type='log')
        yaxis = dict(title=y_label)
    fig = go.Figure(heat)
    fig.update_layout(
        title=title,
        xaxis=xaxis,
        yaxis=yaxis,
        width=720, height=480,
        margin=dict(l=70, r=20, t=50, b=60),
    )
    if err_msg is not None:
        print(err_msg)
    return fig

Environment preflight

Before the parallel grid below, eagerly materialize all four input-pair surfaces for a single fluid via the SVDSBTL prebuild option. This fails loudly here if the CoolProp build / kernel environment can’t produce them (e.g. a stale wheel without PSmass support) — rather than letting every PS / DT panel downstream render as a silent blank. prebuild strips itself from the cache key, so these surfaces are shared with the plain SVDSBTL&HEOS constructions the panels use.

# prebuild=true builds PT + HmassP + DmassT + PSmass for Water in a single
# construction; a build / env failure raises here and fails the docs build,
# which is the intended loud signal (vs. a quietly blank panel).
CP.AbstractState('SVDSBTL&HEOS', 'Water?{"prebuild": true}')
print('preflight OK: all SVDSBTL input-pair surfaces build for Water')

preflight OK: all SVDSBTL input-pair surfaces build for Water

Grid computation

Each (fluid, input_pair) pair is an independent unit of work, so we compute all of them concurrently via joblib (process-based workers; n_jobs=min(len(FLUIDS)*len(INPUT_PAIRS), cpu_count)). Most of the wall-clock here is the one-time SVDSBTL surface build (a dense SVD at production resolution), not the per-point lookups; surfaces are cached and reused on subsequent builds.

import os
from joblib import Parallel, delayed

def _compute(fluid, input_pair):
    return (fluid, input_pair), compute_error_grid(fluid, input_pair)

tasks = [(f, ip) for f in FLUIDS for ip in INPUT_PAIRS]
n_jobs = min(len(tasks), os.cpu_count() or 1)
results = Parallel(n_jobs=n_jobs, backend='loky')(
    delayed(_compute)(f, ip) for f, ip in tasks
)
GRIDS = {key: grid for key, grid in results}
print(f'computed {len(GRIDS)} panels  (n_jobs={n_jobs}, cores={os.cpu_count()})')

computed 32 panels  (n_jobs=4, cores=4)

Water

PT

plot_from_grid('Water', 'PT', GRIDS[('Water', 'PT')])

HP

plot_from_grid('Water', 'HP', GRIDS[('Water', 'HP')])

DT

plot_from_grid('Water', 'DT', GRIDS[('Water', 'DT')])

PS

plot_from_grid('Water', 'PS', GRIDS[('Water', 'PS')])

Argon

PT

plot_from_grid('Argon', 'PT', GRIDS[('Argon', 'PT')])

HP

plot_from_grid('Argon', 'HP', GRIDS[('Argon', 'HP')])

DT

plot_from_grid('Argon', 'DT', GRIDS[('Argon', 'DT')])

PS

plot_from_grid('Argon', 'PS', GRIDS[('Argon', 'PS')])

Hydrogen

PT

plot_from_grid('Hydrogen', 'PT', GRIDS[('Hydrogen', 'PT')])

HP

plot_from_grid('Hydrogen', 'HP', GRIDS[('Hydrogen', 'HP')])

DT

plot_from_grid('Hydrogen', 'DT', GRIDS[('Hydrogen', 'DT')])

PS

plot_from_grid('Hydrogen', 'PS', GRIDS[('Hydrogen', 'PS')])

R1234yf

PT

plot_from_grid('R1234yf', 'PT', GRIDS[('R1234yf', 'PT')])

HP

plot_from_grid('R1234yf', 'HP', GRIDS[('R1234yf', 'HP')])

DT

plot_from_grid('R1234yf', 'DT', GRIDS[('R1234yf', 'DT')])

PS

plot_from_grid('R1234yf', 'PS', GRIDS[('R1234yf', 'PS')])

R245fa

PT

plot_from_grid('R245fa', 'PT', GRIDS[('R245fa', 'PT')])

HP

plot_from_grid('R245fa', 'HP', GRIDS[('R245fa', 'HP')])

DT

plot_from_grid('R245fa', 'DT', GRIDS[('R245fa', 'DT')])

PS

plot_from_grid('R245fa', 'PS', GRIDS[('R245fa', 'PS')])

D6

PT

plot_from_grid('D6', 'PT', GRIDS[('D6', 'PT')])

HP

plot_from_grid('D6', 'HP', GRIDS[('D6', 'HP')])

DT

plot_from_grid('D6', 'DT', GRIDS[('D6', 'DT')])

PS

plot_from_grid('D6', 'PS', GRIDS[('D6', 'PS')])

CO2

PT

plot_from_grid('CO2', 'PT', GRIDS[('CO2', 'PT')])

HP

plot_from_grid('CO2', 'HP', GRIDS[('CO2', 'HP')])

DT

plot_from_grid('CO2', 'DT', GRIDS[('CO2', 'DT')])

PS

plot_from_grid('CO2', 'PS', GRIDS[('CO2', 'PS')])

Helium

PT

plot_from_grid('Helium', 'PT', GRIDS[('Helium', 'PT')])

HP

plot_from_grid('Helium', 'HP', GRIDS[('Helium', 'HP')])

DT

plot_from_grid('Helium', 'DT', GRIDS[('Helium', 'DT')])

PS

plot_from_grid('Helium', 'PS', GRIDS[('Helium', 'PS')])

Notes

• PT panels mask a 2-percent band around the saturation line, since (T, p) inside that band straddles the dome at ULP scale and SVDSBTL is single-phase only.

• HP panels have no explicit mask: the two-phase dome is a true 2D region in (h, p), and SVDSBTL surface queries inside it throw — caught and reported as NaN. The white wedge in each HP heatmap is the dome.

• PS panels mirror HP: the surface tabulates (rho, T, h, u, w) over (p, s) with s as the input axis. No mask — the two-phase dome is a 2D wedge in (s, p), shown as the white region where surface queries are NaN.

• DT panels validate pressure (rho is the input, hence exact): log10 |\u0394p/p| over (T, rho). No mask — the two-phase dome is a valid DT region (SVDSBTL routes it through the dome lever-rule, matching HEOS p_sat(T)), so it shows near-zero error rather than a NaN wedge.

• R1234yf, R245fa, and D6 have narrow native Tmin/Tmax envelopes that may clip the plot domain below the nominal 0.6 Tc lower bound.

• All panels are computed concurrently via joblib (process workers; n_jobs=min(panels, cpu_count)). Most of the run time is the one-time SVDSBTL surface build, not the per-point lookups.

• This notebook is re-executed on every make html build (Web/conf.py pre-executes all .ipynb via jupyter nbconvert --execute). Surfaces are cached at the directory printed in the first code cell (default ~/.CoolProp/SVDTables-docs/, override via COOLPROP_ALTERNATIVE_SVDTABLES_DIRECTORY). First build of a fresh (fluid, input_pair) takes 20-80 s; subsequent builds reuse the serialized surface and run in ~1 s.

• Source generator: Web/coolprop/_gen/gen_SVDSBTLValidation.py.