Closed-domain thermodynamic-pressure rise (sealed heated box)

Results pending GPU execution. This notebook is a scaffold: the analysis is complete but it has not yet been run against simulation output. Run the case with uv run nassu run validation/buoyant_flows/05_closed_domain_pressure/05_closed_domain_pressure.nassu.yaml on a GPU, then execute this notebook.

This case validates the closed-domain low-Mach thermodynamic-pressure closure (models.low_mach.domain_closure: closed, Taha et al. 2024). A sealed, rigid, adiabatic box of ideal gas is heated at a known uniform volumetric rate Q. With no inlet/outlet the total mass is conserved, so the closure evolves the spatially-uniform thermodynamic pressure from total-mass conservation P(t) = M r / I(t), with I(t) = integral (1/T) dV.

The analytic target is the constant-volume first-law result for a sealed rigid adiabatic ideal gas heated at total rate Q_total:

\[\frac{dP}{dt} = (\gamma - 1)\,\frac{Q_\text{total}}{V},\]

exact for adiabatic walls (the conduction term integrates to the zero wall heat flux). With whole-box uniform heating the volumetric rate Q gives Q_total = Q V, so the target reduces to dP/dt = (gamma - 1) Q. The notebook recovers P(t) node-wise from the EOS, fits the slope, and asserts it matches the target to within 2%.

[ ]:

import os
import pathlib
import glob

import numpy as np
import matplotlib.pyplot as plt
import h5py

from nassu.cfg.model import ConfigScheme

import nassu.viz as common

common.use_style()

[ ]:

def _find_project_root() -> pathlib.Path:
    here = pathlib.Path.cwd().resolve()
    for cand in [here, *here.parents]:
        if (cand / "pyproject.toml").exists() and (cand / "nassu").is_dir():
            return cand
    raise RuntimeError(f"Could not locate Nassu project root upward from {here}")


# Run from the project root: the simulation `save_path` (and every output path read
# off the parsed config below) is repo-root-relative.
project_root = _find_project_root()
os.chdir(project_root)
case_dir = project_root / "validation/buoyant_flows/05_closed_domain_pressure"
config_path = case_dir / "05_closed_domain_pressure.nassu.yaml"

sim_cfg = ConfigScheme.sim_cfgs_from_file(str(config_path))[0]
n_steps = sim_cfg.n_steps
print(f"loaded {sim_cfg.name}: n_steps={n_steps}, domain={sim_cfg.domain.domain_size}")

[ ]:

# Analytic quantities read straight off the parsed low-Mach config, so the target is
# never hard-coded: gamma = cp / cv with cv = cp - r, Q is the whole-box volumetric
# source rate, and the target slope is (gamma - 1) Q.
lm = sim_cfg.models.low_mach
r_gas = lm.r
cp = lm.cp
cv = cp - r_gas
gamma = cp / cv
P0 = lm.P_thermo
Q = lm.energy.source_region_rate  # volumetric heat-release rate (energy/vol/step)

ds = sim_cfg.domain.domain_size
V = ds.x * ds.y * ds.z
Q_total = Q * V  # whole-box uniform heating

dPdt_analytic = (gamma - 1.0) * Q_total / V  # = (gamma - 1) Q
print(f"r={r_gas}, cp={cp}, cv={cv}, gamma={gamma}")
print(f"V={V}, Q={Q:.3e}, Q_total={Q_total:.3e}")
print(f"analytic dP/dt = (gamma-1) Q_total/V = {dPdt_analytic:.6e} per step")

[ ]:

# Locate the instantaneous volume export off the parsed config and read each snapshot.
# The volume export stores the full domain as a single block per time group:
# `t<step>/block0/<macr>`, each a 3-D array in (z, y, x) layout.
export = sim_cfg.output.instantaneous["default"]
h5_path = sorted(glob.glob(str(export.base_path) + "*.h5"))[0]
print("reading", h5_path)


def read_snapshots(path):
    with h5py.File(path, "r") as h:
        times = sorted({k.split("/")[0] for k in h.keys()}, key=lambda s: float(s[1:]))
        steps = np.array([float(t[1:]) for t in times])
        data = {}
        for t in times:
            g = h[t]["block0"]
            data[t] = {m: g[m][:] for m in ("rho", "ux", "uy", "uz", "temp")}
    return steps, times, data


steps, times, snaps = read_snapshots(h5_path)
print(f"{len(steps)} snapshots, steps {steps[0]:.0f} .. {steps[-1]:.0f}")

[ ]:

# Recover the spatially-uniform thermodynamic pressure node-wise from the EOS
# rho = P/(r T)  =>  P = r * rho * T, exact since P is uniform. Also collect the peak
# speed (must stay quiescent) and the spatial spread of rho/T (must stay uniform).
cs = 1.0 / np.sqrt(3.0)
P_t, umax, rho_spread, T_spread = [], [], [], []
any_nan = False
for t in times:
    s = snaps[t]
    rho, T = s["rho"], s["temp"]
    P_t.append(float(r_gas * np.nanmean(rho * T)))
    speed = np.sqrt(s["ux"] ** 2 + s["uy"] ** 2 + s["uz"] ** 2)
    umax.append(float(np.nanmax(speed)))
    rho_spread.append(float(np.nanstd(rho) / np.nanmean(rho)))
    T_spread.append(float(np.nanstd(T) / np.nanmean(T)))
    any_nan = any_nan or bool(np.isnan(rho).any() or np.isnan(T).any())

P_t = np.array(P_t)
umax = np.array(umax)
rho_spread = np.array(rho_spread)
T_spread = np.array(T_spread)

# Least-squares slope of P(t) over the whole run (uniform heating -> linear from step 0).
slope, intercept = np.polyfit(steps, P_t, 1)
rel_err = abs(slope - dPdt_analytic) / abs(dPdt_analytic)
print(f"P(0) recovered    = {P_t[0]:.6f}  (config P0 = {P0})")
print(f"fitted   dP/dt    = {slope:.6e} per step")
print(f"analytic dP/dt    = {dPdt_analytic:.6e} per step")
print(f"relative error    = {rel_err * 100:.3f} %")
print(f"peak |u| over run = {umax.max():.3e}  (Ma_max = {umax.max() / cs:.3e})")

Results

[ ]:

fig, ax = plt.subplots(1, 3, figsize=(16, 4))

# (1) P(t) recovered vs the analytic linear rise.
P_analytic = P0 + dPdt_analytic * steps
ax[0].plot(steps, P_t, "o", ms=4, label=r"$P(t) = r\,\langle\rho T\rangle$ (recovered)")
ax[0].plot(steps, P_analytic, "--", c="k", lw=1.5, label=r"analytic $P_0 + (\gamma-1)\,Q\,t$")
ax[0].set_xlabel("step")
ax[0].set_ylabel("thermodynamic pressure $P$")
ax[0].set_title(f"Closed-domain pressure rise (err {rel_err * 100:.2f}%)")
ax[0].legend()
ax[0].grid(alpha=0.3)

# (2) Residual of the recovered pressure about the analytic rise, with the +/-2%
# acceptance band scaled to the analytic rise at each step (absolute units).
residual = P_t - P_analytic
tol_band = 0.02 * dPdt_analytic * steps
ax[1].plot(steps, residual, "-o", ms=4, c="C3", label=r"$P(t) - P_\mathrm{analytic}$")
ax[1].fill_between(steps, -tol_band, tol_band, color="grey", alpha=0.2, label="2% band")
ax[1].axhline(0.0, ls="--", c="k", lw=1)
ax[1].set_xlabel("step")
ax[1].set_ylabel(r"$P - P_\mathrm{analytic}$")
ax[1].set_title("Pressure residual vs analytic rise")
ax[1].legend()
ax[1].grid(alpha=0.3)

# (3) Quiescence and spatial-uniformity diagnostics.
ax[2].plot(steps, umax, "-o", ms=3, label=r"$\max|u|$")
ax[2].plot(steps, rho_spread, "-s", ms=3, label=r"$\sigma(\rho)/\langle\rho\rangle$")
ax[2].plot(steps, T_spread, "-^", ms=3, label=r"$\sigma(T)/\langle T\rangle$")
ax[2].axhline(0.1 * cs, ls="--", c="grey", lw=1, label=r"$\mathrm{Ma}=0.1$")
ax[2].set_xlabel("step")
ax[2].set_ylabel("value")
ax[2].set_title("Box stays quiescent and spatially uniform")
ax[2].legend()
ax[2].grid(alpha=0.3)
fig.tight_layout()

Summary

Acceptance criteria for the closed-domain pressure closure:

  • Pressure-rise rate - the fitted dP/dt matches the analytic (gamma - 1) Q to within 2%.

  • No divergence - all snapshots are finite (no NaN in rho/T).

  • Quiescent box - max|u| stays far below Ma = 0.1 (no spurious dilatation flow or reduced-pressure checkerboard mode).

  • Spatial uniformity - rho and T stay spatially uniform (small relative spread), confirming the closure keeps the uniformly-heated sealed box homogeneous.

[ ]:

assert not any_nan, "divergence: NaN found in the fields"
assert rel_err < 0.02, f"dP/dt off analytic by {rel_err * 100:.2f}% (> 2% tolerance)"
assert umax.max() < 0.1 * cs, f"velocity not subsonic: max|u|={umax.max():.3e}"
assert rho_spread.max() < 1e-2, f"rho not spatially uniform: max spread {rho_spread.max():.3e}"
assert T_spread.max() < 1e-2, f"T not spatially uniform: max spread {T_spread.max():.3e}"
print("VALIDATION PASSED: closed-domain dP/dt matches (gamma-1) Q within 2%.")

Version

[ ]:

sim_info = sim_cfg.output.read_info()
print(sim_info)

Configuration

[ ]:

from IPython.display import Code

Code(filename=str(config_path), language="yaml")