SVDSBTL backend — validation against HEOS#
This page compares the SVDSBTL surrogate backend against the HEOS Helmholtz-energy reference for rho as a function of two different input pairs: (P, T) and (H, P). Each fluid below shows log10 | (rho_SVDSBTL - rho_HEOS) / rho_HEOS | over single-phase rectangles in the two coordinate systems.
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 and HmassP_INPUTS. Output: rhomass() [kg/m^3].Fluids:
Water, Argon, Hydrogen, R1234yf, R245fa, D6.[1]:
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
# 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']
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)
---------------------------------------------------------------------------
ModuleNotFoundError Traceback (most recent call last)
Cell In[1], line 16
13 Path(_cache).mkdir(parents=True, exist_ok=True)
15 import numpy as np
---> 16 import plotly.io as pio
17 import plotly.graph_objects as go
18 import CoolProp.CoolProp as CP
ModuleNotFoundError: No module named 'plotly'
[2]:
def _hp_axes(ref, Tmin, Tmax, pmin, pmax, N_Y):
"""Probe the four corners of the (T, p) rectangle to estimate an
enthalpy range that covers the single-phase region for the HP plot.
"""
hs = []
for T_q, p_q in [(Tmin, pmin), (Tmin, pmax), (Tmax, pmin), (Tmax, pmax)]:
try:
ref.update(PT_INPUTS, p_q, T_q)
hs.append(ref.hmass())
except Exception:
pass
if not hs:
return np.linspace(0.0, 1e6, N_Y)
h_min, h_max = min(hs), max(hs)
if h_max == h_min:
# All surviving probes returned the same enthalpy. Expand by a
# small relative delta so np.linspace produces a non-degenerate
# axis; the panel will be uninformative but at least renders.
delta = max(abs(h_min) * 1e-6, 1e-6)
h_min -= delta
h_max += delta
return np.linspace(h_min, h_max, N_Y)
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' or 'HP'.
'PT' -> x = p (log, in Pa), y = T (linear, in K).
Mask a +/- BAND_REL band around the saturation line.
'HP' -> x = p (log, in Pa), y = h (linear, in J/kg).
No explicit saturation mask; the SVDSBTL surface is undefined
inside the dome and queries there are caught and set to NaN.
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
p_axis = np.geomspace(pmin, pmax, N_X)
if input_pair == 'PT':
y_axis = np.linspace(Tmin, Tmax, N_Y)
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':
y_axis = _hp_axes(ref, Tmin, Tmax, pmin, pmax, N_Y)
y_label = 'specific enthalpy h [kJ/kg]'
psat = None
else:
raise ValueError(f'unknown input_pair: {input_pair!r}')
x_label = 'reduced pressure p / p_c'
err = np.full((N_Y, N_X), np.nan)
try:
svd = CP.AbstractState('SVDSBTL&HEOS', fluid)
except Exception as e:
return (p_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, p in enumerate(p_axis):
if psat is not None and np.isfinite(psat[i]) and abs(p - psat[i]) / psat[i] < BAND_REL:
continue
try:
if input_pair == 'PT':
ref.update(PT_INPUTS, p, y)
svd.update(PT_INPUTS, p, y)
else: # HP
ref.update(HmassP_INPUTS, y, p)
svd.update(HmassP_INPUTS, y, p)
rho_ref = ref.rhomass()
rho_svd = svd.rhomass()
if rho_ref > 0:
err[i, j] = abs(rho_svd - rho_ref) / rho_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 being shown as NaN/white — which would falsely suggest a failed
# query. SVDSBTL inside the two-phase dome returns the same density as
# HEOS (both go through the same saturation lookup), so err is 0 there.
with np.errstate(divide='ignore', invalid='ignore'):
log_err = np.log10(np.maximum(err, 1e-16))
return p_axis, y_axis, log_err, Tc, pc, x_label, y_label, None
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[2], line 24
21 h_max += delta
22 return np.linspace(h_min, h_max, N_Y)
---> 24 def compute_error_grid(fluid, input_pair, N_X=N_X, N_Y=N_Y):
25 """Return (x_axis, y_axis, log10_err, Tc, pc, x_label, y_label, err_msg).
26
27 input_pair is 'PT' or 'HP'.
(...) 35 the returned grid is all-NaN in that case.
36 """
37 ref = CP.AbstractState('HEOS', fluid)
NameError: name 'N_X' is not defined
[3]:
def plot_from_grid(fluid, input_pair, grid):
p_axis, y_axis, log_err, Tc, pc, x_label, y_label, err_msg = grid
pr = p_axis / pc
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 (log), y = T (linear). z is stored as [y_index, x_index].
Tr = y_axis / Tc
heat = go.Heatmap(
x=pr, y=Tr, 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)
else: # HP — Mollier convention: x = h (linear), y = p (log).
h_disp = y_axis / 1e3 # J/kg -> kJ/kg for display
heat = go.Heatmap(
x=h_disp, y=pr, 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')
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
Parallel grid computation#
Each (fluid, input_pair) pair is an independent unit of work, so we compute all 12 of them concurrently via joblib (process-based workers; n_jobs=min(len(FLUIDS)*len(INPUT_PAIRS), cpu_count)). Wall-clock and per-panel timings are printed below.
[4]:
import time, os
from joblib import Parallel, delayed
def compute_timed(fluid, input_pair):
t0 = time.perf_counter()
grid = compute_error_grid(fluid, input_pair)
return (fluid, input_pair), grid, time.perf_counter() - t0
tasks = [(f, ip) for f in FLUIDS for ip in INPUT_PAIRS]
n_jobs = min(len(tasks), os.cpu_count() or 1)
t0 = time.perf_counter()
results = Parallel(n_jobs=n_jobs, backend='loky')(
delayed(compute_timed)(f, ip) for f, ip in tasks
)
wall = time.perf_counter() - t0
GRIDS = {key: grid for key, grid, _ in results}
print(f'wall-clock: {wall:.2f} s (n_jobs={n_jobs}, '
f'cores={os.cpu_count()}, panels={len(tasks)})')
for key, _, dt in results:
print(f' {key[0]:9s} {key[1]:>3s} {dt:6.2f} s')
print(f' sum {sum(dt for _, _, dt in results):6.2f} s '
f'speedup vs serial: {sum(dt for _, _, dt in results) / wall:.2f}x')
---------------------------------------------------------------------------
ModuleNotFoundError Traceback (most recent call last)
Cell In[4], line 2
1 import time, os
----> 2 from joblib import Parallel, delayed
4 def compute_timed(fluid, input_pair):
5 t0 = time.perf_counter()
ModuleNotFoundError: No module named 'joblib'
Water#
PT#
[5]:
plot_from_grid('Water', 'PT', GRIDS[('Water', 'PT')])
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[5], line 1
----> 1 plot_from_grid('Water', 'PT', GRIDS[('Water', 'PT')])
NameError: name 'GRIDS' is not defined
HP#
[6]:
plot_from_grid('Water', 'HP', GRIDS[('Water', 'HP')])
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[6], line 1
----> 1 plot_from_grid('Water', 'HP', GRIDS[('Water', 'HP')])
NameError: name 'GRIDS' is not defined
Argon#
PT#
[7]:
plot_from_grid('Argon', 'PT', GRIDS[('Argon', 'PT')])
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[7], line 1
----> 1 plot_from_grid('Argon', 'PT', GRIDS[('Argon', 'PT')])
NameError: name 'GRIDS' is not defined
HP#
[8]:
plot_from_grid('Argon', 'HP', GRIDS[('Argon', 'HP')])
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[8], line 1
----> 1 plot_from_grid('Argon', 'HP', GRIDS[('Argon', 'HP')])
NameError: name 'GRIDS' is not defined
Hydrogen#
PT#
[9]:
plot_from_grid('Hydrogen', 'PT', GRIDS[('Hydrogen', 'PT')])
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[9], line 1
----> 1 plot_from_grid('Hydrogen', 'PT', GRIDS[('Hydrogen', 'PT')])
NameError: name 'GRIDS' is not defined
HP#
[10]:
plot_from_grid('Hydrogen', 'HP', GRIDS[('Hydrogen', 'HP')])
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[10], line 1
----> 1 plot_from_grid('Hydrogen', 'HP', GRIDS[('Hydrogen', 'HP')])
NameError: name 'GRIDS' is not defined
R1234yf#
PT#
[11]:
plot_from_grid('R1234yf', 'PT', GRIDS[('R1234yf', 'PT')])
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[11], line 1
----> 1 plot_from_grid('R1234yf', 'PT', GRIDS[('R1234yf', 'PT')])
NameError: name 'GRIDS' is not defined
HP#
[12]:
plot_from_grid('R1234yf', 'HP', GRIDS[('R1234yf', 'HP')])
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[12], line 1
----> 1 plot_from_grid('R1234yf', 'HP', GRIDS[('R1234yf', 'HP')])
NameError: name 'GRIDS' is not defined
R245fa#
PT#
[13]:
plot_from_grid('R245fa', 'PT', GRIDS[('R245fa', 'PT')])
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[13], line 1
----> 1 plot_from_grid('R245fa', 'PT', GRIDS[('R245fa', 'PT')])
NameError: name 'GRIDS' is not defined
HP#
[14]:
plot_from_grid('R245fa', 'HP', GRIDS[('R245fa', 'HP')])
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[14], line 1
----> 1 plot_from_grid('R245fa', 'HP', GRIDS[('R245fa', 'HP')])
NameError: name 'GRIDS' is not defined
D6#
PT#
[15]:
plot_from_grid('D6', 'PT', GRIDS[('D6', 'PT')])
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[15], line 1
----> 1 plot_from_grid('D6', 'PT', GRIDS[('D6', 'PT')])
NameError: name 'GRIDS' is not defined
HP#
[16]:
plot_from_grid('D6', 'HP', GRIDS[('D6', 'HP')])
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[16], line 1
----> 1 plot_from_grid('D6', 'HP', GRIDS[('D6', 'HP')])
NameError: name 'GRIDS' is not defined
CO2#
PT#
[17]:
plot_from_grid('CO2', 'PT', GRIDS[('CO2', 'PT')])
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[17], line 1
----> 1 plot_from_grid('CO2', 'PT', GRIDS[('CO2', 'PT')])
NameError: name 'GRIDS' is not defined
HP#
[18]:
plot_from_grid('CO2', 'HP', GRIDS[('CO2', 'HP')])
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[18], line 1
----> 1 plot_from_grid('CO2', 'HP', GRIDS[('CO2', 'HP')])
NameError: name 'GRIDS' is not defined
Helium#
PT#
[19]:
plot_from_grid('Helium', 'PT', GRIDS[('Helium', 'PT')])
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[19], line 1
----> 1 plot_from_grid('Helium', 'PT', GRIDS[('Helium', 'PT')])
NameError: name 'GRIDS' is not defined
HP#
[20]:
plot_from_grid('Helium', 'HP', GRIDS[('Helium', 'HP')])
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[20], line 1
----> 1 plot_from_grid('Helium', 'HP', GRIDS[('Helium', 'HP')])
NameError: name 'GRIDS' is not defined
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.
R1234yf,R245fa, andD6have narrow nativeTmin/Tmaxenvelopes that may clip the plot domain below the nominal0.6 Tclower bound.The 12 panels are computed concurrently via
joblib(process workers;n_jobs=min(panels, cpu_count)). Wall-clock and per-panel breakdown are in the parallel-compute cell above.This notebook is re-executed on every
make htmlbuild (Web/conf.pypre-executes all.ipynbviajupyter nbconvert --execute). Surfaces are cached at the directory printed in the first code cell (default~/.CoolProp/SVDTables-docs/, override viaCOOLPROP_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.