Status: planned - not yet validated against reference data. The case is posed on the variable-density low-Mach closure (models.low_mach, Taha et al. 2024): EOS density rho = P/(r T), the exact (rho - rho_inf) g buoyancy and a finite-difference energy field carrying the burner heat-release source. A validated run also needs closure stability at the fire-scale density ratio, a thermal wall BC on the voxelized room shell, an EOS-consistent open doorway, and the 62.9 kW heat-release-to-lattice mapping; see the case README.md.

The analysis below reads the low-Mach temp energy macroscopic and the models.low_mach configuration, and reconstructs the dimensional profiles through the variable-density route (the EOS, the emergent real temperature ratio). The plotted numbers are finalised with the case’s first validated GPU run, when the burner rate and lattice gravity are calibrated to the matched Q*_H and the steady-state averaging window is measured.

Case 03 - Steckler room fire (Steckler, Quintiere & Rinkinen 1982)

Validation of the open, internally-heated compartment flow on the variable-density low-Mach closure (models.low_mach, Taha et al. 2024): a single room with a doorway and a steady floor-centred heat release develops a stratified two-layer flow that spills hot fluid out of the top of the doorway and draws ambient air in through the bottom, with a neutral plane in between.

Scope. This is posed as a real fire, not a reduced-ratio Boussinesq surrogate: the hot layer carries its real T_hot / T_inf ~ 2, the temperature closes the equation of state rho = P / (r T), and the resulting density variation drives the flow through the exact (rho - rho_inf) g buoyancy (see the case README.md). The heat release is a volumetric source on the finite-difference energy field and the comparison is on non-dimensional quantities: the doorway velocity / temperature profile shapes and the neutral-plane height fraction z_n / H_door.

The notebook time-averages the dense doorway plane_series, extracts the door-centreline vertical profiles and the neutral-plane height, and compares them against the digitised Steckler reference when it is available (reference/doorway_profiles.csv).

import pathlib

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

import nassu.viz as common
from nassu.cfg.model import ConfigScheme

common.use_style()

import importlib.util
import os


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}")


project_root = _find_project_root()
# The simulation `save_path` (and therefore every probe/data path read from the
# config below) is repo-root-relative, so run from the project root regardless
# of where the notebook was launched (case dir under Jupyter / nbconvert).
os.chdir(project_root)
case_dir = project_root / "validation/buoyant_flows/03_steckler_room_fire"
case_path = case_dir / "03_steckler_room_fire.nassu.yaml"
reference_csv = case_dir / "reference/doorway_profiles.csv"

# Geometry is in real metres and the grid is driven by one knob, `dx` (level-0
# lattice spacing). Read `dx` from the YAML and map physical (metre) coordinates
# to level-0 lattice with the case's `lattice_mapping` helper, so the analysis
# follows `dx` automatically (no hardcoded lattice numbers).
_spec = importlib.util.spec_from_file_location(
    "lattice_mapping", case_dir / "geometry" / "lattice_mapping.py"
)
lm = importlib.util.module_from_spec(_spec)
_spec.loader.exec_module(lm)

from nassu.cfg.loader import load_cfg_yaml

dx = float(load_cfg_yaml(case_path)["variables"]["dx"])

# Doorway opening in lattice: x in [-0.37, 0.37] m, z in [0, 1.83] m, in the
# front wall y = 0 (door centre at the room centre x = 0).
DOOR_X = (lm.to_lat(-0.37, "x", dx), lm.to_lat(0.37, "x", dx))
H_DOOR = 1.83 / dx  # doorway height (z, lattice units)
DOOR_XC = lm.to_lat(0.0, "x", dx)  # door centreline x (= room centre)

(name, sim_id), cfg = next(iter(ConfigScheme.sim_cfgs_from_file_dct(str(case_path)).items()))
print(f"loaded sim {name} (#{sim_id}); dx = {dx} m, door height H_DOOR = {H_DOOR:g} lattice")

Load the steady doorway plane

The doorway plane_series (plane normal \(y\), sitting in the front-wall doorway) carries the through-door velocity component uy and the low-Mach energy-field temperature temp over the \(x\)-\(z\) opening. We average all snapshots after the spin-up to get the statistically-steady field, exactly as the other thermal plane notebooks do.

def load_doorway_plane(cfg):
    """Return (xs, zs, T, uy) time-mean grids for the doorway plane.

    The plane normal is y, so the in-plane coordinates are x and z and the
    through-door velocity is the out-of-plane component `uy`. The doorway faces
    -y (the plenum is at y < 0, the room at y > 0), so `uy > 0` is inflow (cold
    ambient drawn into the room through the lower door) and `uy < 0` is outflow
    (hot layer spilling out of the upper door). Averages all snapshots after 60%
    of the run.
    """
    plane = cfg.output.exports["plane_series"].series.planes["doorway"]
    pts = pd.read_csv(plane.points_filename)

    xs = np.sort(pts["x"].unique())
    zs = np.sort(pts["z"].unique())
    ix = {v: i for i, v in enumerate(xs)}
    iz = {v: i for i, v in enumerate(zs)}
    cols = [str(int(i)) for i in pts["idx"]]
    col_xi = np.array([ix[x] for x in pts["x"]])
    col_zi = np.array([iz[z] for z in pts["z"]])

    def _mean_grid(macr):
        df = plane.read_full_data(macr)
        start = int(df["time_step"].max() * 0.6)  # drop the spin-up
        df = df[df["time_step"] >= start]
        arr = df[cols].to_numpy(dtype=float).mean(axis=0)
        grid = np.full((zs.size, xs.size), np.nan)
        grid[col_zi, col_xi] = arr
        return grid

    return xs, zs, _mean_grid("temp"), _mean_grid("uy")


xs, zs, T_grid, uy_grid = load_doorway_plane(cfg)
print(f"doorway plane grid: {xs.size} (x) x {zs.size} (z)")

Door-centreline profiles and the neutral plane

At the door centreline (the room centre \(x = 0\) m) we take the vertical profiles of the through-door velocity and the temperature over the opening height, normalise them, and locate the neutral plane as the height where the through-door velocity changes sign (outflow above, inflow below).

def door_centreline_profiles(xs, zs, T_grid, uy_grid):
    """Vertical profiles at the door centreline over the opening height."""
    xi = int(np.argmin(np.abs(xs - DOOR_XC)))
    mask = (zs >= 0.0) & (zs <= H_DOOR)
    z = zs[mask]
    uy = uy_grid[mask, xi]
    T = T_grid[mask, xi]
    return z, uy, T


def neutral_plane_height(z, uy):
    """Height where the through-door velocity crosses zero (linear interp)."""
    s = np.sign(uy)
    crossings = np.where(np.diff(s) != 0)[0]
    if crossings.size == 0:
        return np.nan
    # Take the dominant crossing (largest velocity swing around it).
    k = crossings[np.argmax(np.abs(uy[crossings + 1] - uy[crossings]))]
    z0, z1 = z[k], z[k + 1]
    u0, u1 = uy[k], uy[k + 1]
    return z0 - u0 * (z1 - z0) / (u1 - u0)


z, uy, T = door_centreline_profiles(xs, zs, T_grid, uy_grid)
z_n = neutral_plane_height(z, uy)

# Non-dimensionalise: height by door height; velocity by its peak magnitude;
# temperature by the centreline span (ambient -> hot layer). This matches the
# normalisation of the digitised reference (`reference/doorway_profiles.csv`).
z_over_H = z / H_DOOR
u_norm = uy / np.nanmax(np.abs(uy))
T_span = np.nanmax(T) - np.nanmin(T)
T_norm = (T - np.nanmin(T)) / T_span if T_span > 0 else T * 0.0

print(f"neutral-plane height: z_n = {z_n:.1f} lattice  ->  z_n / H_door = {z_n / H_DOOR:.3f}")
print("(Steckler test-16 experiment: z_n / H_door = 0.560)")

# Through-door mass balance (uniform spacing). uy > 0 is inflow, uy < 0 is
# outflow; np.trapezoid (NumPy 2.x; np.trapz was removed).
inflow = np.trapezoid(np.clip(uy, 0.0, None), z)
outflow = -np.trapezoid(np.clip(uy, None, 0.0), z)
print(
    f"door flux balance: inflow={inflow:.3g}  outflow={outflow:.3g}  "
    f"imbalance={abs(outflow - inflow) / max(outflow, 1e-30):.1%}"
)

# Reference overlay, read defensively (empty until the profiles are digitised).
ref = None
try:
    _r = pd.read_csv(reference_csv, comment="#")
    if len(_r) > 0:
        ref = _r
except Exception as exc:  # noqa: BLE001
    print(f"reference not loaded ({exc}); plotting Nassu profiles alone")

fig, axes = plt.subplots(1, 2, figsize=(10, 5), sharey=True)

axes[0].plot(u_norm, z_over_H, "-o", ms=3, label="Nassu")
if ref is not None:
    axes[0].plot(ref["u_norm"], ref["z_over_H"], "ks", ms=4, label="Steckler 1982")
axes[0].axvline(0.0, color="0.6", lw=0.8)
if np.isfinite(z_n):
    axes[0].axhline(z_n / H_DOOR, color="C3", ls="--", lw=1, label="neutral plane")
axes[0].set_xlabel("through-door velocity (normalised)")
axes[0].set_ylabel(r"$z / H_\mathrm{door}$")
axes[0].legend()

axes[1].plot(T_norm, z_over_H, "-o", ms=3, label="Nassu")
if ref is not None:
    axes[1].plot(ref["T_norm"], ref["z_over_H"], "ks", ms=4, label="Steckler 1982")
if np.isfinite(z_n):
    axes[1].axhline(z_n / H_DOOR, color="C3", ls="--", lw=1)
axes[1].set_xlabel("temperature (normalised)")
axes[1].legend()

fig.suptitle("Steckler doorway profiles (door centreline)")
fig.tight_layout()

Dimensional reconstruction by similarity

The run is set up for dynamic similarity with Steckler test 16: the dimensionless heat-release rate \(Q^*_H\), the Prandtl number and the geometry length-ratios all match the experiment (see the case README.md), so the dimensionless fields the solver produces are the experiment’s. That means the dimensional doorway profiles can be reconstructed and overlaid on the raw NIST/FDS measurements in \(^\circ\)C and m/s.

Because the closure is variable-density, the lattice temperature is already the temperature ratio \(\theta = T/T_\mathrm{inf}\) (the ambient is the lattice reference \(T_\mathrm{ref} = 1\)), so the dimensional temperature is simply \(T = \theta\,T_\mathrm{inf}\) - there is no Boussinesq \(\theta_\mathrm{eff} = \beta\,\phi\) expansion-coefficient step. The hot-layer excess ratio \(\theta_\mathrm{hot} - 1 = (T_\mathrm{hot} - T_\mathrm{inf})/T_\mathrm{inf}\) is read emergently off the door-centreline temperature (not fitted), and the density follows the EOS \(\rho = P/(rT) = 1/\theta\) - the warm layer thins to \(\rho_\mathrm{hot} \approx \rho_\mathrm{inf}/2\) at \(\theta_\mathrm{hot} \approx 2\). The reconstruction scales are

\[\Delta T = (\theta_\mathrm{hot}-1)\,T_\mathrm{inf}, \qquad U_b = \sqrt{g\,H_\mathrm{door}\,(\theta_\mathrm{hot}-1)},\]

with \(T_\mathrm{inf} = 301\) K (28 \(^\circ\)C ambient, the experiment) and the physical \(g = 9.81\) m/s\(^2\), \(H_\mathrm{door} = 1.83\) m. With a real temperature ratio this reconstruction carries no Boussinesq linearization error; the residual modelling limit is the reduced Reynolds number. The temperature scale grows like \((\theta_\mathrm{hot}-1)\) but the velocity scale only like its square root.

# Variable-density similarity reconstruction: turn the matched-Q* dimensionless
# profiles into dimensional ones (deg C, m/s) using the run's OWN emergent
# hot-layer temperature RATIO (no Boussinesq beta), and overlay on the raw
# NIST/FDS dimensional data.
T0_C = 28.0  # ambient temperature (deg C), Steckler test 16
T_INF_K = 301.15  # ambient absolute temperature (K); the lattice T_ref = 1 maps here
G_PHYS = 9.81  # gravity (m/s^2)
H_DOOR_M = 1.83  # door height (m)

# In the variable-density closure the lattice temperature IS the ratio
# theta = T / T_inf (ambient T_ref = 1), so the EOS density is rho = 1 / theta
# and the emergent hot-layer excess ratio is read directly off the centreline.
theta_inf = 1.0
theta_hot = np.nanmax(T)
excess = theta_hot - theta_inf  # (T_hot - T_inf) / T_inf, dimensionless

# Predicted physical scales (not fitted to the experiment).
dT_pred = excess * T_INF_K  # hot-layer temperature rise (K)
Ub_pred = np.sqrt(G_PHYS * H_DOOR_M * excess)  # door buoyant velocity (m/s)

print(f"emergent theta_hot = T_hot/T_inf = {theta_hot:.3f}   (experiment ~ 2)")
print(f"EOS hot-layer density rho_hot = 1/theta_hot = {1.0 / theta_hot:.3f}   (~ 0.5)")
print(f"-> predicted dT = {dT_pred:.0f} K     (experiment ~ 90 K)")
print(f"-> predicted U_b = {Ub_pred:.2f} m/s  (experiment ~ 2.3 m/s, measured peak 2.12)")

# Reconstructed Nassu dimensional profiles. Temperature is EOS-direct
# (T = theta * T_inf); the through-door velocity uses the buoyant scale (nassu
# convention: u > 0 = inflow).
u_phys = u_norm * Ub_pred  # m/s
T_phys = T * T_INF_K - 273.15  # deg C, straight from the temperature ratio

# Raw experimental dimensional profiles (door-centreline column *_C).
fds = pd.read_csv(case_dir / "reference" / "steckler_test16_fds_exp.csv")
fds = fds[(fds["Height"] <= H_DOOR_M) & fds["V_C"].notna() & fds["T_C"].notna()].sort_values(
    "Height"
)
z_over_H_exp = fds["Height"].to_numpy() / H_DOOR_M
u_exp = -fds["V_C"].to_numpy()  # flip FDS (outflow +) to nassu (inflow +)
T_exp = fds["T_C"].to_numpy()  # deg C

fig, axes = plt.subplots(1, 2, figsize=(10, 5), sharey=True)

axes[0].plot(u_phys, z_over_H, "-o", ms=3, label="Nassu (reconstructed)")
axes[0].plot(u_exp, z_over_H_exp, "ks", ms=4, label="Steckler 1982 (NIST/FDS)")
axes[0].axvline(0.0, color="0.6", lw=0.8)
if np.isfinite(z_n):
    axes[0].axhline(z_n / H_DOOR, color="C3", ls="--", lw=1, label="neutral plane")
axes[0].set_xlabel("through-door velocity (m/s)")
axes[0].set_ylabel(r"$z / H_\mathrm{door}$")
axes[0].legend()

axes[1].plot(T_phys, z_over_H, "-o", ms=3, label="Nassu (reconstructed)")
axes[1].plot(T_exp, z_over_H_exp, "ks", ms=4, label="Steckler 1982 (NIST/FDS)")
if np.isfinite(z_n):
    axes[1].axhline(z_n / H_DOOR, color="C3", ls="--", lw=1)
axes[1].set_xlabel(r"temperature ($^\circ$C)")
axes[1].legend()

fig.suptitle("Steckler doorway profiles - dimensional reconstruction by similarity")
fig.tight_layout()

Centreline plane - plume and stratification

The room mid-plane through the burner and the door centre (the room centre \(x = 0\) m, plane normal \(x\)) shows the buoyant plume, the hot upper layer and the exchange through the doorway into the ambient plenum.

def load_centreline_plane(cfg):
    plane = cfg.output.exports["plane_series"].series.planes["centreline"]
    pts = pd.read_csv(plane.points_filename)
    ys = np.sort(pts["y"].unique())
    zs = np.sort(pts["z"].unique())
    iy = {v: i for i, v in enumerate(ys)}
    iz = {v: i for i, v in enumerate(zs)}
    cols = [str(int(i)) for i in pts["idx"]]
    col_yi = np.array([iy[y] for y in pts["y"]])
    col_zi = np.array([iz[z] for z in pts["z"]])

    df = plane.read_full_data("temp")
    start = int(df["time_step"].max() * 0.6)
    df = df[df["time_step"] >= start]
    arr = df[cols].to_numpy(dtype=float).mean(axis=0)
    grid = np.full((zs.size, ys.size), np.nan)
    grid[col_zi, col_yi] = arr
    return ys, zs, grid


ys, zc, Tc = load_centreline_plane(cfg)

fig, ax = plt.subplots(figsize=(8, 4))
pc = ax.pcolormesh(ys, zc, Tc, shading="auto", cmap="inferno")
fig.colorbar(pc, ax=ax, label=r"$T$ (temperature)")
ax.set_xlabel("y (lattice) - plenum -> doorway -> room interior -> back wall")
ax.set_ylabel("z (lattice)")
ax.set_title("Time-mean temperature, room centreline (x = 0 m)")
fig.tight_layout()

Centreline plane - instantaneous velocity

The same room mid-plane (\(x = 0\) m), but the last instantaneous velocity snapshot rather than the time mean: the buoyant plume rising off the burner, the ceiling jet, and the two-way exchange through the doorway into the plenum. The colour is the velocity magnitude; the white unit arrows show the in-plane direction.

def load_centreline_velocity_last(cfg):
    """Last-snapshot velocity components on the room centreline plane (x = 0).

    Unlike the temperature map above (time-mean), this returns the single most
    recent instantaneous frame so the plume and doorway jet show their live
    structure rather than a smoothed average.
    """
    plane = cfg.output.exports["plane_series"].series.planes["centreline"]
    pts = pd.read_csv(plane.points_filename)
    ys = np.sort(pts["y"].unique())
    zs = np.sort(pts["z"].unique())
    iy = {v: i for i, v in enumerate(ys)}
    iz = {v: i for i, v in enumerate(zs)}
    cols = [str(int(i)) for i in pts["idx"]]
    col_yi = np.array([iy[y] for y in pts["y"]])
    col_zi = np.array([iz[z] for z in pts["z"]])

    comp = {}
    step = 0
    for macr in ("ux", "uy", "uz"):
        df = plane.read_full_data(macr)
        step = int(df["time_step"].max())
        row = df.loc[df["time_step"].idxmax(), cols].to_numpy(dtype=float)
        grid = np.full((zs.size, ys.size), np.nan)
        grid[col_zi, col_yi] = row
        comp[macr] = grid
    return ys, zs, comp["ux"], comp["uy"], comp["uz"], step


ys, zs, ux, uy, uz, step = load_centreline_velocity_last(cfg)
speed = np.sqrt(ux**2 + uy**2 + uz**2)

fig, ax = plt.subplots(figsize=(8, 4))
pc = ax.pcolormesh(ys, zs, speed, shading="auto", cmap="viridis")
fig.colorbar(pc, ax=ax, label=r"$|u|$ (lattice)")

# In-plane velocity (uy horizontal, uz vertical), coarsely subsampled and scaled
# by magnitude so the plume and doorway jet stand out and the quiescent ambient
# stays clean. Masked solid/exterior nodes are zeroed (no arrow).
sy = max(1, ys.size // 32)
sz = max(1, zs.size // 24)
Yq, Zq = np.meshgrid(ys[::sy], zs[::sz])
ax.quiver(
    Yq, Zq, np.nan_to_num(uy[::sz, ::sy]), np.nan_to_num(uz[::sz, ::sy]), color="w", width=0.003
)

ax.set_xlabel("y (lattice) - plenum -> doorway -> room interior -> back wall")
ax.set_ylabel("z (lattice)")
ax.set_title(f"Instantaneous velocity magnitude, room centreline (x = 0), step {step}")
fig.tight_layout()

Version and configuration

# Version / configuration, read from the run's info.yaml (same as the other
# validation notebooks).
sim_info = cfg.output.read_info()
print("nassu version:", sim_info["version"])
print("commit hash:", sim_info["commit"])
print("case:", case_path.relative_to(project_root))
print("domain:", cfg.domain.domain_size)

lm = cfg.models.low_mach
nu = (cfg.models.LBM.tau - 0.5) / 3.0
rho_inf = lm.rho_inf if lm.rho_inf is not None else lm.P_thermo / (lm.r * lm.T_ref)
print("operator:", cfg.models.LBM.coll_oper, " bulk_viscosity:", cfg.models.LBM.bulk_viscosity)
print("tau (fluid):", cfg.models.LBM.tau, " nu:", nu)
print("EOS: r =", lm.r, " P_thermo =", lm.P_thermo, " T_ref =", lm.T_ref, " rho_inf =", rho_inf)
print("Pr:", lm.Pr, " thermal diffusivity alpha = nu/Pr:", nu / lm.Pr)
print("gravity:", lm.gravity)
print("energy source rate:", lm.energy.source_regions[0].rate)