Variable-density low-Mach closure - warm-bubble sanity check

This is a sanity check, not a benchmark. There is no closed-form reference solution. The case validates two things about the variable-density low-Mach thermal closure (Taha et al. 2024), run with the in-collision bulk-viscosity stabilization (models.LBM.bulk_viscosity) enabled:

  1. Stability - the closure runs the full 6000 steps with no divergence; the peak speed stays bounded and subsonic, i.e. there is no exponentially growing reduced-pressure checkerboard / acoustic mode (with the bulk viscosity off the same setup diverges within a few hundred steps).

  2. Physical buoyancy - a warm Gaussian bubble seeded low in the box, made lighter by the equation of state rho = P/(r T), rises (mean uz > 0) while diffusing toward ambient.

[ ]:

import glob
import os
import pathlib

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

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

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/01_lowmach_bubble"
config_path = case_dir / "01_lowmach_bubble.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}")

[ ]:

# 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.exports["default"].volumes["default"].inst
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}")

[ ]:

# Per-snapshot diagnostics. All are reductions over the whole field, so the
# (z, y, x) axis layout is immaterial here.
umax, uz_mean, t_max, rho_min, rho_max = [], [], [], [], []
any_nan = False
for t in times:
    s = snaps[t]
    speed = np.sqrt(s["ux"] ** 2 + s["uy"] ** 2 + s["uz"] ** 2)
    umax.append(float(np.nanmax(speed)))
    uz_mean.append(float(np.nanmean(s["uz"])))
    t_max.append(float(np.nanmax(s["temp"])))
    rho_min.append(float(np.nanmin(s["rho"])))
    rho_max.append(float(np.nanmax(s["rho"])))
    any_nan = any_nan or bool(np.isnan(s["rho"]).any() or np.isnan(s["uz"]).any())

umax = np.array(umax)
uz_mean = np.array(uz_mean)
t_max = np.array(t_max)
rho_min = np.array(rho_min)
rho_max = np.array(rho_max)
cs = 1.0 / np.sqrt(3.0)
print(f"peak |u| over run = {umax.max():.3e}  (Ma_max = {umax.max() / cs:.3e})")
print(f"final mean uz     = {uz_mean[-1]:.3e}  (buoyant rise, > 0)")
print(f"NaN encountered   = {any_nan}")

Results

[ ]:

fig, ax = plt.subplots(1, 2, figsize=(11, 4))

ax[0].plot(steps, umax, "-o", ms=3, label=r"$\max|u|$")
ax[0].plot(steps, uz_mean, "-s", ms=3, label=r"mean $u_z$ (buoyant rise)")
ax[0].axhline(0.1 * cs, ls="--", c="grey", lw=1, label=r"$\mathrm{Ma}=0.1$")
ax[0].set_xlabel("step")
ax[0].set_ylabel("lattice velocity")
ax[0].set_title("Velocity stays bounded, subsonic; mean $u_z>0$")
ax[0].legend()
ax[0].grid(alpha=0.3)

ax[1].plot(steps, t_max, "-o", ms=3, label=r"$T_{\max}$")
ax[1].plot(steps, rho_min, "-s", ms=3, label=r"$\rho_{\min}$ (warm core)")
ax[1].plot(steps, rho_max, "-^", ms=3, label=r"$\rho_{\max}$")
ax[1].axhline(1.0, ls="--", c="grey", lw=1, label="ambient")
ax[1].set_xlabel("step")
ax[1].set_ylabel("value")
ax[1].set_title("Bubble diffuses toward ambient; warm core stays lighter")
ax[1].legend()
ax[1].grid(alpha=0.3)
fig.tight_layout()

[ ]:

# Vertical mid-plane (y = Ny/2) temperature slices: the warm bubble diffuses
# and gently rises along -gravity (+z). Array layout is (z, y, x).
sel = [times[0], times[len(times) // 3], times[2 * len(times) // 3], times[-1]]
fig, axes = plt.subplots(1, len(sel), figsize=(13, 4.2), sharey=True)
ny = snaps[times[0]]["temp"].shape[1]
vmin, vmax = 1.0, t_max[0]
for ax, t in zip(axes, sel):
    sl = snaps[t]["temp"][:, ny // 2, :]  # (z, x)
    im = ax.imshow(sl, origin="lower", aspect="auto", cmap="inferno", vmin=vmin, vmax=vmax)
    ax.set_title(f"step {int(float(t[1:]))}")
    ax.set_xlabel("x")
axes[0].set_ylabel("z (gravity is -z)")
fig.colorbar(im, ax=axes, label="T", fraction=0.025)
fig.suptitle("Mid-plane temperature: warm bubble diffuses and rises")

Summary

Success criteria for this sanity check:

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

  • Bounded, subsonic velocity - max|u| stays well below Ma = 0.1, with no exponential growth (the checkerboard / acoustic mode is suppressed by the in-collision bulk viscosity).

  • Physical buoyant rise - the mean vertical velocity is positive.

  • Diffusive relaxation - the peak temperature decays toward ambient and the warm core stays EOS-consistently lighter than ambient (rho_min < 1).

[ ]:

assert not any_nan, "divergence: NaN found in the fields"
assert umax.max() < 0.1 * cs, f"velocity not subsonic: max|u|={umax.max():.3e}"
assert umax.max() < 5e-3, f"velocity grew unexpectedly large: {umax.max():.3e}"
assert uz_mean[-1] > 0, f"no buoyant rise: final mean uz={uz_mean[-1]:.3e}"
assert t_max[-1] < t_max[0], "peak temperature did not relax toward ambient"
assert rho_min.min() < 1.0, "warm core is not lighter than ambient"
print("SANITY CHECK PASSED: stable, subsonic, buoyant, diffusing.")

Version

[ ]:

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

Configuration

[ ]:

from IPython.display import Code

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