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.

[1]:

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()

[2]:

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

loaded lowmach_bubble: n_steps=6000, domain=x=48 y=48 z=96

[3]:

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

reading /data/internal/waine/validation/buoyant_flows/01_lowmach_bubble/results/lowmach_bubble/lowmach_bubble__000/outputs/default.volume.default.inst.000.h5
21 snapshots, steps 0 .. 6000

[4]:

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

peak |u| over run = 1.078e-03  (Ma_max = 1.868e-03)
final mean uz     = 9.073e-04  (buoyant rise, > 0)
NaN encountered   = False

Results

[5]:

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()
../../../_images/validation_buoyant_flows_01_lowmach_bubble_01_lowmach_bubble_6_0.png 
[6]:

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

[6]:

Text(0.5, 0.98, 'Mid-plane temperature: warm bubble diffuses and rises')
../../../_images/validation_buoyant_flows_01_lowmach_bubble_01_lowmach_bubble_7_1.png

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).

[7]:

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

SANITY CHECK PASSED: stable, subsonic, buoyant, diffusing.

Version

[8]:

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

{'commit': '5e47f2762575c2d285254d36bfd354b6af09fda1', 'sim_aux_info': {'n_blocks': 432, 'n_blocks_per_level': {0: 432}, 'n_nodes': 221184, 'n_nodes_per_block': 512, 'n_nodes_per_level': {0: 221184}, 'n_nodes_updated_per_step': 221184}, 'sim_balancer': [{'blocks_per_lvl': {0: 432}, 'memory_weight': 432, 'time_weight': 432}], 'sim_configs': {'checkpoint': {'export': {'finish_save': False, 'interval': {'end_step': 0, 'frequency': 0, 'start_step': 0}, 'keep_only_last_checkpoint': False}, 'load': {'checkpoint_start': False, 'folderpath': None, 'reset_time_step': False}}, 'data': {'body_nodes': {}, 'derived_macrs': {}, 'export_IBM_nodes': {}, 'export_rescales': {}, 'exports': {'default': {'interval': {'constant_dt': True, 'end_step': 0, 'frequency': 300.0, 'lvl': 0, 'start_step': 0}, 'interval_flush': 250, 'interval_group': 1000, 'macrs': ['rho', 'u', 'temp'], 'outputs': {'instantaneous': True, 'stats': None}, 'rescale': None, 'system': None, 'target': {'bodies': {}, 'csvs': {}, 'lines': {}, 'planes': {}, 'points': {}, 'volumes': {'default': {'max_h5_size_gb': 4.0, 'max_lvl': -1, 'volume': {'end': [1.0, 1.0, 1.0], 'is_abs': False, 'start': [0.0, 0.0, 0.0], 'system': None, 'transformation': None}}}}}}, 'guards': {}, 'monitors': {'fields': {}}}, 'debug': {'IBM': {'no_force_spread': False, 'no_nodes_export': False}, 'LBM': {'collision_only': False, 'no_interblock_comm': False, 'no_macrs_export': False, 'streaming_only': False}, 'code_generation': {'load_generated_files': False, 'save_generated_files': False}, 'isolate': None, 'multiblock': {'export_comm_vtk': False, 'export_used_nodes': False, 'run_communication': True}, 'output_IBM_nodes': False, 'output_only': False, 'profile': False, 'profile_kernels': False}, 'domain': {'block_forest_load': None, 'block_forest_save': None, 'block_size': 8, 'bodies': {}, 'bodies_domain_limits': {'end': [1.0, 1.0, 1.0], 'is_abs': False, 'start': [0.0, 0.0, 0.0], 'system': None, 'transformation': None}, 'default_system': 'lbm', 'domain_extent': None, 'domain_extent_system': None, 'domain_size': {'x': 48, 'y': 48, 'z': 96}, 'global_transformations': [], 'point_clouds': {}, 'refinement': {'static': {}}, 'systems': {}}, 'models': {'BC': {'BC_map': [], 'WM_cfg': {'NeqWM_pres_grad_mult': 1.0, 'NeqWM_u_friction_floor': 0.0001, 'TDMA_max_div': 51, 'TDMA_max_error': 5e-06, 'TDMA_max_iters': 20, 'TDMA_min_div': 21, 'TDMA_stretch_beta': 2.0, 'TDMA_yp_target': 0.2}, 'inlet_turbulence': None, 'periodic_dims': [True, True, True], 'rho_normalization': []}, 'IBM': {'body_cfgs': {}, 'dirac_delta': '3_points', 'forces_accomodate_time': 0, 'forces_spread_limit': 0.1, 'min_dirac_sum': 0.99, 'reset_forces': True}, 'LBM': {'F': {'x': 0.0, 'y': 0.0, 'z': 0.0}, 'bulk_viscosity': 0.16667, 'coll_oper': 'HRRBGK', 'coll_oper_params': {'mode_hrrbgk': 'dynamic', 'sigma_hrrbgk': 0.99}, 'tau': 0.6, 'thermal_model': False, 'vel_set': 'D3Q27'}, 'LES': {'model': 'Smagorinsky', 'sgs_cte': 0.1}, 'energy': None, 'engine': {'devices_numbers': [3], 'n_devices': 1, 'name': 'CUDA'}, 'initialization': {'equations': {'rho': '1', 'ux': '0', 'uy': '0', 'uz': '0'}, 'inlet_field': False, 'macrs_filename': None}, 'low_mach': {'P_thermo': 1.0, 'Pr': 0.71, 'T_ref': 1.0, 'cp': 1.0, 'domain_closure': 'open', 'energy': {'initial_field': '1.0 + 0.05*exp(-(((x-24)**2)+((y-24)**2)+((z-24)**2))/72.0)', 'source_regions': [], 'wall_bcs': []}, 'gravity': [0.0, 0.0, -0.0002], 'mu_exponent': 0.0, 'mu_ref': None, 'r': 1.0, 'rho_inf': None}, 'multiblock': {'custom_overlap_F2C': {}, 'mark_nodes_as_unused': True, 'overlap_F2C': 2}, 'precision': {'calculations': 'double', 'default': 'single', 'macroscopics': 'default', 'populations': 'default'}, 'rheology': None, 'scalar_transports': {}, 'volumetric_regions': []}, 'n_steps': 6000, 'name': 'lowmach_bubble', 'parent': None, 'report': {'end_step': 0, 'frequency': 1000, 'start_step': 0, 'warmup_steps': 0}, 'run_simul': True, 'save_path': 'validation/buoyant_flows/01_lowmach_bubble/results/lowmach_bubble', 'sim_id': 0}, 'sim_memory': {'border_idxs': {'bytes_memory': 1184.0, 'num_buffers': 1, 'total_memory': '0.00 Mb'}, 'bulk_idxs': {'bytes_memory': 1728.0, 'num_buffers': 1, 'total_memory': '0.00 Mb'}, 'comm_same_lvl': {'bytes_memory': 6925824.0, 'num_buffers': 3, 'total_memory': '6.60 Mb'}, 'comm_same_lvl_scalar_energy': {'bytes_memory': 6483456.0, 'num_buffers': 1, 'total_memory': '6.18 Mb'}, 'indexes': {'bytes_memory': 1728.0, 'num_buffers': 1, 'total_memory': '0.00 Mb'}, 'lt_map': {'bytes_memory': 221184.0, 'num_buffers': 1, 'total_memory': '0.21 Mb'}, 'macrs_full_domain': {'bytes_memory': 17694720.0, 'num_buffers': 1, 'total_memory': '16.88 Mb'}, 'macrs_scalar_energy': {'bytes_memory': 6193152.0, 'num_buffers': 1, 'total_memory': '5.91 Mb'}, 'structs': {'bytes_memory': 128504.0, 'num_buffers': 3, 'total_memory': '0.12 Mb'}}, 'sim_runtime': {'devices': [{'AsyncEngineCount': 2, 'CanMapHostMemory': 1, 'ClockRate': 1665000, 'ComputeMode': 0, 'ConcurrentKernels': 1, 'DirectManagedMemAccessFromHost': 0, 'EccEnabled': 0, 'GlobalL1CacheSupported': 1, 'GlobalMemoryBusWidth': 384, 'GpuOverlap': 1, 'IsMultiGpuBoard': 0, 'L2CacheSize': 6291456, 'LocalL1CacheSupported': 1, 'ManagedMemory': 1, 'MaxBlockDimX': 1024, 'MaxBlockDimY': 1024, 'MaxBlockDimZ': 64, 'MaxGridDimX': 2147483647, 'MaxGridDimY': 65535, 'MaxGridDimZ': 65535, 'MaxRegistersPerBlock': 65536, 'MaxRegistersPerMultiprocessor': 65536, 'MaxSharedMemoryPerBlock': 49152, 'MaxSharedMemoryPerBlockOptin': 101376, 'MaxSharedMemoryPerMultiprocessor': 102400, 'MaxThreadsPerBlock': 1024, 'MaxThreadsPerMultiProcessor': 1536, 'MemoryClockRate': 8001000, 'MultiProcessorCount': 80, 'Name': 'NVIDIA RTX A5500', 'PciBusId': 196, 'PciDeviceId': 0, 'SingleToDoublePrecisionPerfRatio': 64, 'TotalConstantMemory': 65536, 'WarpSize': 32}], 'diverged_scalars': {}, 'events': [{'device': None, 'name': 'comm.unpack_interp_fork_events[0]', 'recorded_by': [], 'waited_by': []}, {'device': None, 'name': 'comm.unpack_interp_join_events[0][0]', 'recorded_by': [], 'waited_by': []}, {'device': None, 'name': 'comm.unpack_interp_join_events[0][1]', 'recorded_by': [], 'waited_by': []}, {'device': None, 'name': 'comm.unpack_interp_join_events[0][2]', 'recorded_by': [], 'waited_by': []}, {'device': None, 'name': 'comm.unpack_interp_join_events[0][3]', 'recorded_by': [], 'waited_by': []}, {'device': 0, 'name': 'lbm.lbm_main_bulk_events[0]', 'recorded_by': ['lbm.lbm_main_bulk_queues[0]'], 'waited_by': []}, {'device': 0, 'name': 'scalar.main_bulk_events[energy][0]', 'recorded_by': ['scalar.main_bulk_queues[energy][0]'], 'waited_by': []}], 'host_overhead': {}, 'info': {'MLUPS': 59.479637145996094, 'curr_step': 6000, 'sim_aux_info': {'n_blocks': 432, 'n_blocks_per_level': {0: 432}, 'n_nodes': 221184, 'n_nodes_per_block': 512, 'n_nodes_per_level': {0: 221184}, 'n_nodes_updated_per_step': 221184}, 'total_runtime': 22.311904191970825}, 'kernel_timings': {}, 'kernels': {'nassu_LBM_init_S_kernel': {'call': {'block_size': [256, 1, 1], 'blocks_per_dim': [2, 432, 1], 'grid_size': [512, 432, 1], 'kernel_name': 'nassu_LBM_init_S_kernel', 'only_1D': False, 'perc_wasted_threads': 0.0, 'total_threads': 221184, 'total_wasted_threads': 0, 'wasted_threads': [0, 0, 0]}, 'info': {'constSizeBytes': 204, 'localSizeBytes': 0, 'maxDynamicSharedSizeBytes': 49152, 'maxThreadsPerBlock': 1024, 'numRegs': 28, 'preferredShmemCarveout': -1, 'sharedSizeBytes': 0}, 'n_calls': 1, 'n_grids': 1, 'occupancy': {'active_blocks_per_sm': 6, 'active_warps_per_sm': 48.0, 'occupancy_fraction': 1.0}}, 'nassu_LBM_init_borders_kernel': {'call': {'block_size': [256, 1, 1], 'blocks_per_dim': [2, 432, 1], 'grid_size': [512, 432, 1], 'kernel_name': 'nassu_LBM_init_borders_kernel', 'only_1D': False, 'perc_wasted_threads': 0.0, 'total_threads': 221184, 'total_wasted_threads': 0, 'wasted_threads': [0, 0, 0]}, 'info': {'constSizeBytes': 204, 'localSizeBytes': 0, 'maxDynamicSharedSizeBytes': 49152, 'maxThreadsPerBlock': 384, 'numRegs': 164, 'preferredShmemCarveout': -1, 'sharedSizeBytes': 0}, 'n_calls': 1, 'n_grids': 1, 'occupancy': {'active_blocks_per_sm': 1, 'active_warps_per_sm': 8.0, 'occupancy_fraction': 0.1667}}, 'nassu_LBM_init_macrs_kernel': {'call': {'block_size': [256, 1, 1], 'blocks_per_dim': [2, 432, 1], 'grid_size': [512, 432, 1], 'kernel_name': 'nassu_LBM_init_macrs_kernel', 'only_1D': False, 'perc_wasted_threads': 0.0, 'total_threads': 221184, 'total_wasted_threads': 0, 'wasted_threads': [0, 0, 0]}, 'info': {'constSizeBytes': 204, 'localSizeBytes': 0, 'maxDynamicSharedSizeBytes': 49152, 'maxThreadsPerBlock': 1024, 'numRegs': 36, 'preferredShmemCarveout': -1, 'sharedSizeBytes': 0}, 'n_calls': 1, 'n_grids': 1, 'occupancy': {'active_blocks_per_sm': 6, 'active_warps_per_sm': 48.0, 'occupancy_fraction': 1.0}}, 'nassu_LBM_main_bulk_kernel|lvls=(0,)': {'call': {'block_size': [512, 1, 1], 'blocks_per_dim': [1, 432, 1], 'grid_size': [512, 432, 1], 'kernel_name': 'nassu_LBM_main_bulk_kernel', 'only_1D': False, 'perc_wasted_threads': 0.0, 'total_threads': 221184, 'total_wasted_threads': 0, 'wasted_threads': [0, 0, 0]}, 'info': {'constSizeBytes': 204, 'localSizeBytes': 144, 'maxDynamicSharedSizeBytes': 53248, 'maxThreadsPerBlock': 512, 'numRegs': 128, 'preferredShmemCarveout': -1, 'sharedSizeBytes': 0}, 'n_calls': 6000, 'n_grids': 1, 'occupancy': {'active_blocks_per_sm': 1, 'active_warps_per_sm': 16.0, 'occupancy_fraction': 0.3333}}, 'nassu_LBM_update_borders_marked_kernel|lvls=(0,)': {'call': {'block_size': [296, 1, 1], 'blocks_per_dim': [1, 432, 1], 'grid_size': [296, 432, 1], 'kernel_name': 'nassu_LBM_update_borders_marked_kernel', 'only_1D': False, 'perc_wasted_threads': 0.0, 'total_threads': 127872, 'total_wasted_threads': 0, 'wasted_threads': [0, 0, 0]}, 'info': {'constSizeBytes': 204, 'localSizeBytes': 24, 'maxDynamicSharedSizeBytes': 49152, 'maxThreadsPerBlock': 296, 'numRegs': 168, 'preferredShmemCarveout': -1, 'sharedSizeBytes': 0}, 'n_calls': 6000, 'n_grids': 1, 'occupancy': {'active_blocks_per_sm': 1, 'active_warps_per_sm': 9.25, 'occupancy_fraction': 0.1927}}, 'nassu_scalar_energy_energy_correction_kernel|energy=energy|lvls=(0,)': {'call': {'block_size': [256, 1, 1], 'blocks_per_dim': [2, 432, 1], 'grid_size': [512, 432, 1], 'kernel_name': 'nassu_scalar_energy_energy_correction_kernel', 'only_1D': False, 'perc_wasted_threads': 0.0, 'total_threads': 221184, 'total_wasted_threads': 0, 'wasted_threads': [0, 0, 0]}, 'info': {'constSizeBytes': 204, 'localSizeBytes': 0, 'maxDynamicSharedSizeBytes': 49152, 'maxThreadsPerBlock': 896, 'numRegs': 70, 'preferredShmemCarveout': -1, 'sharedSizeBytes': 0}, 'n_calls': 6000, 'n_grids': 1, 'occupancy': {'active_blocks_per_sm': 3, 'active_warps_per_sm': 24.0, 'occupancy_fraction': 0.5}}, 'nassu_scalar_energy_init_kernel': {'call': {'block_size': [256, 1, 1], 'blocks_per_dim': [2, 432, 1], 'grid_size': [512, 432, 1], 'kernel_name': 'nassu_scalar_energy_init_kernel', 'only_1D': False, 'perc_wasted_threads': 0.0, 'total_threads': 221184, 'total_wasted_threads': 0, 'wasted_threads': [0, 0, 0]}, 'info': {'constSizeBytes': 204, 'localSizeBytes': 0, 'maxDynamicSharedSizeBytes': 49152, 'maxThreadsPerBlock': 512, 'numRegs': 108, 'preferredShmemCarveout': -1, 'sharedSizeBytes': 0}, 'n_calls': 1, 'n_grids': 1, 'occupancy': {'active_blocks_per_sm': 2, 'active_warps_per_sm': 16.0, 'occupancy_fraction': 0.3333}}, 'nassu_scalar_energy_main_bulk_kernel|scalar=energy|lvls=(0,)': {'call': {'block_size': [512, 1, 1], 'blocks_per_dim': [1, 432, 1], 'grid_size': [512, 432, 1], 'kernel_name': 'nassu_scalar_energy_main_bulk_kernel', 'only_1D': False, 'perc_wasted_threads': 0.0, 'total_threads': 221184, 'total_wasted_threads': 0, 'wasted_threads': [0, 0, 0]}, 'info': {'constSizeBytes': 204, 'localSizeBytes': 8, 'maxDynamicSharedSizeBytes': 53248, 'maxThreadsPerBlock': 512, 'numRegs': 128, 'preferredShmemCarveout': -1, 'sharedSizeBytes': 0}, 'n_calls': 6000, 'n_grids': 1, 'occupancy': {'active_blocks_per_sm': 1, 'active_warps_per_sm': 16.0, 'occupancy_fraction': 0.3333}}, 'nassu_update_finite_differences_kernel|lvls=(0,)': {'call': {'block_size': [512, 1, 1], 'blocks_per_dim': [1, 432, 1], 'grid_size': [512, 432, 1], 'kernel_name': 'nassu_update_finite_differences_kernel', 'only_1D': False, 'perc_wasted_threads': 0.0, 'total_threads': 221184, 'total_wasted_threads': 0, 'wasted_threads': [0, 0, 0]}, 'info': {'constSizeBytes': 204, 'localSizeBytes': 0, 'maxDynamicSharedSizeBytes': 49152, 'maxThreadsPerBlock': 1024, 'numRegs': 24, 'preferredShmemCarveout': -1, 'sharedSizeBytes': 0}, 'n_calls': 6000, 'n_grids': 1, 'occupancy': {'active_blocks_per_sm': 3, 'active_warps_per_sm': 48.0, 'occupancy_fraction': 1.0}}}, 'streams': [{'device': 0, 'name': 'comm.pack_queues[0]'}, {'device': 0, 'name': 'comm.unpack_queues[0]'}, {'device': 0, 'name': 'comm.unpack_interp_queues[0][0]'}, {'device': 0, 'name': 'comm.unpack_interp_queues[0][1]'}, {'device': 0, 'name': 'comm.unpack_interp_queues[0][2]'}, {'device': 0, 'name': 'comm.unpack_interp_queues[0][3]'}, {'device': 0, 'name': 'lbm.lbm_init_queues[0]'}, {'device': 0, 'name': 'lbm.lbm_main_bulk_queues[0]'}, {'device': 0, 'name': 'bc.bc_queues[0][0]'}, {'device': 0, 'name': 'scalar.init_queues[energy][0]'}, {'device': 0, 'name': 'scalar.main_bulk_queues[energy][0]'}, {'counts_per_device': {0: 11}}]}, 'sim_setup_timings': {'balancer': {'seconds': 0.001574, 'tag': 'python'}, 'block_forest_build': {'seconds': 0.122195, 'tag': 'python'}, 'codegen': {'seconds': 36.529412, 'tag': 'python'}, 'config_parse': {'seconds': 0.055664, 'tag': 'python'}, 'engine_init': {'seconds': 0.270602, 'tag': 'cuda'}, 'export_setup': {'seconds': 0.022029, 'tag': 'python'}, 'ibm_build': {'seconds': 0.00034, 'tag': 'python'}, 'init_kernels': {'seconds': 0.001912, 'tag': 'cuda'}, 'init_macrs': {'seconds': 0.165855, 'tag': 'python'}, 'kernel_calls_build': {'seconds': 0.007943, 'tag': 'python'}, 'lattice_type_map': {'seconds': 0.041839, 'tag': 'python'}, 'lbm_handler': {'seconds': 1.138989, 'tag': 'python'}, 'mem_alloc': {'seconds': 0.498316, 'tag': 'cuda'}, 'nvrtc_compile': {'seconds': 11.056222, 'tag': 'cuda'}, 'probes_setup': {'seconds': 4.2e-05, 'tag': 'python'}, 'sem_seed': {'seconds': 1.2e-05, 'tag': 'python'}, 'templates_init': {'seconds': 0.002449, 'tag': 'python'}}, 'sim_timedeltas': {'build': 49.75455403327942, 'full': 72.79738020896912, 'initialization': 0.16783452033996582, 'runtime': 22.410943031311035}, 'sim_timestamps': {'build_end': datetime.datetime(2026, 7, 23, 0, 50, 1, 20192, tzinfo=datetime.timezone(datetime.timedelta(0), '+00:00')), 'build_start': datetime.datetime(2026, 7, 23, 0, 49, 11, 265633, tzinfo=datetime.timezone(datetime.timedelta(0), '+00:00')), 'full_end': datetime.datetime(2026, 7, 23, 0, 50, 24, 63002, tzinfo=datetime.timezone(datetime.timedelta(0), '+00:00')), 'full_start': datetime.datetime(2026, 7, 23, 0, 49, 11, 265628, tzinfo=datetime.timezone(datetime.timedelta(0), '+00:00')), 'initialization_end': datetime.datetime(2026, 7, 23, 0, 50, 1, 516155, tzinfo=datetime.timezone(datetime.timedelta(0), '+00:00')), 'initialization_start': datetime.datetime(2026, 7, 23, 0, 50, 1, 348328, tzinfo=datetime.timezone(datetime.timedelta(0), '+00:00')), 'runtime_end': datetime.datetime(2026, 7, 23, 0, 50, 24, 62996, tzinfo=datetime.timezone(datetime.timedelta(0), '+00:00')), 'runtime_start': datetime.datetime(2026, 7, 23, 0, 50, 1, 652060, tzinfo=datetime.timezone(datetime.timedelta(0), '+00:00'))}, 'timestamp': datetime.datetime(2026, 7, 23, 0, 50, 24, 63005, tzinfo=datetime.timezone(datetime.timedelta(0), '+00:00')), 'version': '2.0.0a7'}

Configuration

[9]:

from IPython.display import Code

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

[9]:

# Variable-density low-Mach closure - warm-bubble sanity check
# =============================================================
#
# SANITY CHECK (not a benchmark): there is no closed-form reference solution.
# This case validates that the variable-density low-Mach thermal closure
# (Taha et al. 2024) runs STABLY over a sustained run and produces
# the PHYSICALLY EXPECTED buoyant response, with the in-collision bulk-viscosity
# stabilization enabled (`models.LBM.bulk_viscosity`).
#
# Setup: a smooth warm Gaussian temperature bubble (+5% over T_ref, sigma ~ 6
# lattice units) seeded low in a periodic box; the equation of state
# rho = P / (r T) makes the warm core lighter, and gravity along -z lifts it.
# No continuous heat source - the bubble simply diffuses and rises.
#
# Expected (physical) behaviour, asserted by the notebook:
#   - the run completes all steps with NO divergence (no NaN in rho / u);
#   - the peak speed stays bounded and subsonic (|u| ~ 1e-3, Ma << 0.1) -
#     i.e. NO exponentially growing checkerboard / acoustic mode;
#   - the mean vertical velocity is positive and grows smoothly (a buoyant
#     rise), the asymmetric mean-upward signal that is the genuine physics;
#   - the peak temperature relaxes toward ambient (heat diffusion) and the
#     warm core stays EOS-consistently lighter than ambient.
#
# With `bulk_viscosity` OFF the same setup diverges within a few hundred steps
# (the reduced-pressure checkerboard instability); see the thermodynamics
# stability theory page and the reproducer examples/debug_lowmach/.

simulations:
  - name: lowmach_bubble
    save_path: ./validation/buoyant_flows/01_lowmach_bubble/results/lowmach_bubble

    n_steps: 6000

    report:
      frequency: 1000

    domain:
      domain_size:
        x: 48
        y: 48
        z: 96
      block_size: 8

    data:
      exports:
        default:
          macrs: [rho, u, temp]
          interval:
            frequency: 300
            lvl: 0
          target:
            volumes:
              default: {}
          outputs:
            instantaneous: true

    models:
      precision:
        default: single
        calculations: double

      LBM:
        F: {x: 0, y: 0, z: 0}
        tau: 0.6
        vel_set: D3Q27
        coll_oper: HRRBGK
        # In-collision bulk viscosity (omega_bulk = 1.0): damps the reduced-pressure
        # checkerboard through the trace of the non-equilibrium stress.
        bulk_viscosity: 0.16667

      LES:
        model: Smagorinsky
        sgs_cte: 0.1

      engine:
        name: CUDA

      # Variable-density low-Mach closure. The reduced-pressure checkerboard is
      # damped by the in-collision bulk viscosity set above (models.LBM.bulk_viscosity).
      low_mach:
        r: 1.0
        P_thermo: 1.0
        T_ref: 1.0
        Pr: 0.71
        gravity: [0, 0, -2.0e-4]
        energy:
          initial_field: "1.0 + 0.05*exp(-(((x-24)**2)+((y-24)**2)+((z-24)**2))/72.0)"