Multiblock Passive Scalar Diffusion (Case 12b)¶
This notebook validates the scalar advection-diffusion module on a multiblock grid against the same exact analytical solution as Case 12a. A single Fourier mode diffuses on a triply periodic 3-D box, but the inner cube \([24, 40]^3\) is refined to level 1, so every step exercises the scalar same-level border comm at every block face and the F2C / C2F pack-and-unpack path at the six refinement-interface faces (plus their edges and corners).
The initial scalar field on the periodic domain \([0, 64)^3\) is
with the analytical decay (eigenfunction of the Laplacian, fluid at rest)
The fluid runs the standard isothermal RR-BGK scheme on D3Q27; the scalar uses D3Q7 RR-BGK with \(D_\phi = 10^{-3}\) in lattice units at level 0. The level-1 region uses the rescaled lattice diffusivity \(D_\phi(\text{lvl}=1) = 2 D_\phi(\text{lvl}=0)\) so the physical diffusivity is the same on both levels.
Pass criteria
Fitted decay rate \(\lambda_{\text{lbm}} = D k^2\) within a few percent of the analytical value (the strongest single multiblock check; insensitive to refinement topology).
The cosine spatial profile remains continuous across the refinement interfaces: no jumps or oscillations near \(x = 24\) or \(x = 40\).
Volume-weighted \(L_2\) scalar error of similar magnitude to the matching \(N = 64\) single-level run from Case 12a, which gives \(\approx 4.3 \times 10^{-4}\).
Setup¶
[1]:
import pathlib
import matplotlib.pyplot as plt
import numpy as np
import nassu.viz as common
from nassu.cfg.model import ConfigScheme
common.use_style()
Load simulation configuration¶
[2]:
filename = "validation/scalar_transport/01_passive_scalar_transport/01.1_passive_scalar_transport_multiblock.nassu.yaml"
sim_cfgs = ConfigScheme.sim_cfgs_from_file_dct(filename)
sim_cfg = sim_cfgs["passiveScalarDiffusion3DMultiblock", 0]
N = sim_cfg.domain.domain_size.x
D_LBM = sim_cfg.models.scalar_transports["scalar"].adv_diff_equation.D
N_STEPS = sim_cfg.n_steps
PHI_AMPLITUDE = 1.0
k = 2.0 * np.pi / N
decay_rate = D_LBM * k**2
print(
f"N = {N}, D = {D_LBM}, k = {k:.4f}, "
f"decay_rate = {decay_rate:.6e}, n_steps = {N_STEPS}, "
f"exp(-D k^2 t_end) = {np.exp(-decay_rate * N_STEPS):.4f}"
)
N = 64, D = 0.001, k = 0.0982, decay_rate = 9.638286e-06, n_steps = 5120, exp(-D k^2 t_end) = 0.9518
Analytical solution¶
The exact field is the IC modulated by the exponential decay \(e^{-D k^2 t}\).
[3]:
def analytical_phi_at_x(x_phys, t):
"""Exact phi at physical x at time t."""
return np.cos(k * x_phys) * np.exp(-D_LBM * k**2 * t)
def phi_l2_norm_analytical(t):
"""Volume-averaged L2 norm of cos(kx) on the periodic box, decayed to t.
The L2 norm of cos(kx) over a full period is 1 / sqrt(2).
"""
return PHI_AMPLITUDE / np.sqrt(2.0) * np.exp(-D_LBM * k**2 * t)
Multiblock-aware HDF5 / XDMF loaders¶
assemble_3d_field from Case 12a’s notebook does not generalise to refinement: blocks at different levels carry different cell spacings, and the volume-averaged \(L_2\) norm therefore needs explicit cell-volume weighting.
A second subtlety is the on-disk axis order. nassu’s volume export stores each block as a C-order ZYX array of shape (nz, ny, nx) and writes the XDMF Origin / Spacing / Topology Dimensions strings slowest-axis-first, i.e. in (z, y, x) order (see vtk_image_cell_arrays_to_zyx and _build_block_grid_xml in nassu/simul/export/macrsExport.py and nassu/utils/xdmf.py). To recover a physically-aligned per-block array we therefore (a) reverse the parsed origin /
spacing / shape strings back to (x, y, z) and
transpose the dataset with
phi_disk.transpose(2, 1, 0); after that steparr_phys[ix, iy, iz]returns the value at physical position(ox + (ix+0.5)*dx, oy + (iy+0.5)*dx, oz + (iz+0.5)*dx).
This matters specifically for the non-cubic level-0 merged blocks: an earlier ravel("C").reshape((nx, ny, nz), order="F") trick only recovered the right mapping for cubic blocks and silently scrambled non-cubic ones. The volume-weighted \(L_2\) norm itself is invariant to any such permutation (it is a sum of squared values), so the decay-rate metric is robust either way; the realignment is required for the per-block error and the spatial profile plot.
[4]:
import h5py
from lxml import etree
from nassu.cfg.schemes.simul import MacrExport
def parse_xdmf_blocks(xdmf_path):
"""Parse the XDMF manifest and return per-timestep block metadata.
nassu emits volume blocks in C-order ZYX: the HDF5 dataset has shape
``(nz, ny, nx)`` and the XDMF ``Origin`` / ``Spacing`` / ``Topology
Dimensions`` strings are written slowest-axis-first, i.e. in
``(z, y, x)`` order (see ``_build_block_grid_xml`` and
``vtk_image_cell_arrays_to_zyx`` in ``nassu/utils/xdmf.py``). We reverse
those strings here so the returned ``shape_xyz`` is ``(nx, ny, nz)`` and
``origin_xyz`` is ``(ox, oy, oz)``.
"""
xdmf_path = pathlib.Path(xdmf_path)
tree = etree.parse(str(xdmf_path))
root = tree.getroot()
domain = root.find("Domain")
time_series = domain.find("Grid")
result = {}
for spatial_grid in time_series.findall("Grid"):
time_elem = spatial_grid.find("Time")
t = float(time_elem.get("Value"))
blocks = []
for block_grid in spatial_grid.findall("Grid"):
block_key = block_grid.get("Name")
geom = block_grid.find("Geometry")
origin_str = geom.find("DataItem[@Name='Origin']").text
spacing_str = geom.find("DataItem[@Name='Spacing']").text
# On-disk strings are ZYX; reverse to XYZ.
origin_xyz = np.array(origin_str.split(), dtype=float)[::-1]
spacing_xyz = np.array(spacing_str.split(), dtype=float)[::-1]
topo = block_grid.find("Topology")
shape_node_zyx = tuple(int(d) for d in topo.get("Dimensions").split())
shape_xyz = tuple(d - 1 for d in shape_node_zyx)[::-1]
attr = block_grid.find("Attribute")
data_item = attr.find("DataItem")
ref = data_item.text.strip()
h5_rel, _dataset_path = ref.split(":", 1)
h5_path = xdmf_path.parent / h5_rel
blocks.append(
{
"block_key": block_key,
"h5_path": h5_path,
"origin_xyz": origin_xyz,
"spacing_xyz": spacing_xyz,
"shape_xyz": shape_xyz,
}
)
result[t] = blocks
return result
def load_block_phi_phys(block_info, t, scalar_name="scalar"):
"""Load one block's phi and re-arrange to physical (x, y, z) order.
The HDF5 dataset is stored C-order ZYX (``phi_disk[iz, iy, ix]``), so a
plain axis transpose recovers the physical layout ``arr[ix, iy, iz]``.
This is correct for both cubic and non-cubic merged blocks (the older
``ravel("C").reshape((nx, ny, nz), order="F")`` trick only worked for
cubic blocks and scrambled non-cubic ones).
"""
phi_key = f"{scalar_name}_phi"
ts_key = f"t{t:.6f}"
with h5py.File(block_info["h5_path"], "r") as f:
phi_disk = f[ts_key][block_info["block_key"]][phi_key][()].astype(np.float64)
return np.ascontiguousarray(phi_disk.transpose(2, 1, 0))
def load_volume_weighted_phi_l2(macr_export: MacrExport, scalar_name: str = "scalar"):
"""Volume-weighted L2 norm of phi at every snapshot.
Each cell contributes phi^2 * dx^3 with dx taken from the block's
own spacing, so refined cells (smaller dx) carry less weight than
coarse cells. Returns (times, phi_l2) with the average taken over
the total physical volume. Uses on-disk arrays directly because
sum-of-squares is invariant to the array reshape permutation.
"""
blocks_by_t = parse_xdmf_blocks(macr_export.xdmf_filename)
phi_key = f"{scalar_name}_phi"
times = []
norms = []
for t in sorted(blocks_by_t.keys()):
sum_sq_vol = 0.0
total_vol = 0.0
for block_info in blocks_by_t[t]:
dx = float(block_info["spacing_xyz"][0])
cell_volume = dx**3
ts_key = f"t{t:.6f}"
with h5py.File(block_info["h5_path"], "r") as f:
if block_info["block_key"] not in f.get(ts_key, {}):
continue
if phi_key not in f[ts_key][block_info["block_key"]]:
continue
phi = f[ts_key][block_info["block_key"]][phi_key][()].astype(np.float64)
sum_sq_vol += float((phi**2).sum()) * cell_volume
total_vol += phi.size * cell_volume
if total_vol > 0:
times.append(t)
norms.append(np.sqrt(sum_sq_vol / total_vol))
return np.array(times), np.array(norms)
Per-block error against the analytical solution¶
[5]:
def per_block_phi_error(macr_export: MacrExport, t_target: float, scalar_name: str = "scalar"):
"""Compute the volume-weighted L2 phi error globally and per block.
Returns `(t_used, l2_global, per_block)` where `per_block` is a
list of dicts keyed by `block_key` with the block's level (inferred
from spacing), cell count, max pointwise error, and per-block RMS
error. The global L2 is normalised by the IC amplitude.
"""
blocks_by_t = parse_xdmf_blocks(macr_export.xdmf_filename)
times = sorted(blocks_by_t.keys())
t_use = min(times, key=lambda t: abs(t - t_target))
sum_sq_vol_total = 0.0
total_vol = 0.0
per_block = []
for block_info in blocks_by_t[t_use]:
dx = float(block_info["spacing_xyz"][0])
ox = block_info["origin_xyz"][0]
nx, ny, nz = block_info["shape_xyz"]
phi_phys = load_block_phi_phys(block_info, t_use, scalar_name=scalar_name)
x_phys = ox + (np.arange(nx) + 0.5) * dx
phi_exact = analytical_phi_at_x(x_phys, t_use)
# Broadcast along axis 0 (x), constant in y and z (axes 1, 2).
phi_exact_3d = np.broadcast_to(phi_exact[:, None, None], phi_phys.shape)
err = phi_phys - phi_exact_3d
cell_volume = dx**3
sum_sq_vol_total += float((err**2).sum()) * cell_volume
total_vol += phi_phys.size * cell_volume
lvl = int(round(np.log2(1.0 / dx))) if dx < 1.0 else 0
per_block.append(
{
"block_key": block_info["block_key"],
"lvl": lvl,
"n_cells": int(phi_phys.size),
"max_err": float(np.abs(err).max()),
"rms_err": float(np.sqrt((err**2).mean())),
"phi_min": float(phi_phys.min()),
"phi_max": float(phi_phys.max()),
}
)
global_l2 = np.sqrt(sum_sq_vol_total / total_vol) / PHI_AMPLITUDE
return t_use, global_l2, per_block
Load simulation output and compute global metrics¶
[6]:
macr_export = sim_cfg.output.exports["default"].volumes["default"].inst
print(f"Reading XDMF: {macr_export.xdmf_filename}")
times, phi_l2_sim = load_volume_weighted_phi_l2(macr_export, scalar_name="scalar")
print(f"Loaded {len(times)} snapshots, t in [{times[0]:.0f}, {times[-1]:.0f}]")
print(
f"L2 norm at t = 0 : {phi_l2_sim[0]:.4f} (analytical: {phi_l2_norm_analytical(0.0):.4f})"
)
print(
f"L2 norm at t = {times[-1]:.0f}: {phi_l2_sim[-1]:.4f} "
f"(analytical: {phi_l2_norm_analytical(times[-1]):.4f})"
)
Reading XDMF: /data/internal/waine/validation/scalar_transport/01_passive_scalar_transport/results/passive_scalar_diffusion_multiblock/passiveScalarDiffusion3DMultiblock__000/outputs/default.volume.default.inst.xdmf
Loaded 21 snapshots, t in [0, 5120]
L2 norm at t = 0 : 0.7071 (analytical: 0.7071)
L2 norm at t = 5120: 0.6726 (analytical: 0.6731)
[7]:
t_used, l2_global_err, per_block = per_block_phi_error(
macr_export, t_target=float(N_STEPS), scalar_name="scalar"
)
print(f"At t = {t_used:.0f}: volume-weighted L2 phi error = {l2_global_err:.4e}")
print()
print(
f"{'block':>10} {'lvl':>4} {'n_cells':>10} {'max_err':>12} {'rms_err':>12} {'phi_range':>22}"
)
for b in per_block:
print(
f"{b['block_key']:>10} {b['lvl']:>4} {b['n_cells']:>10} "
f"{b['max_err']:>12.3e} {b['rms_err']:>12.3e} "
f"[{b['phi_min']:>+.3f}, {b['phi_max']:>+.3f}]"
)
At t = 5120: volume-weighted L2 phi error = 3.3048e-02
block lvl n_cells max_err rms_err phi_range
block0 0 98304 4.683e-02 3.620e-02 [-0.604, +0.951]
block1 0 61440 4.682e-02 3.110e-02 [-0.951, +0.947]
block2 0 38400 4.681e-02 3.109e-02 [-0.951, +0.947]
block3 0 38400 4.682e-02 3.110e-02 [-0.951, +0.947]
block4 0 15360 4.682e-02 3.110e-02 [-0.951, +0.947]
block5 0 6144 4.682e-02 3.669e-02 [-0.673, +0.947]
block6 1 32768 4.817e-02 1.341e-02 [-0.952, -0.673]
Decay rate fit¶
Fit \(\log\|\phi(\cdot, t)\|_{L^2} = a - \lambda t\) on the volume-weighted norm, and compare \(\lambda_{\text{lbm}}\) to the analytical \(\lambda_{\text{exact}} = D k^2\). The match should be insensitive to refinement: spatial truncation does not affect the decay rate to leading order.
[8]:
mask = phi_l2_sim > 0
slope, intercept = np.polyfit(times[mask], np.log(phi_l2_sim[mask]), 1)
lam_lbm = -slope
lam_exact = D_LBM * k**2
rel_err = abs(lam_lbm - lam_exact) / lam_exact
print(f"lambda_exact = {lam_exact:.6e}")
print(f"lambda_lbm = {lam_lbm:.6e}")
print(f"relative err = {rel_err:.2%}")
lambda_exact = 9.638286e-06
lambda_lbm = 9.676997e-06
relative err = 0.40%
[9]:
fig, ax = plt.subplots(figsize=(8, 5))
t_decay = 1.0 / (D_LBM * k**2)
t_norm = times / t_decay
phi_l2_exact = np.array([phi_l2_norm_analytical(t) for t in times])
ax.semilogy(t_norm, phi_l2_exact, **common.markers.exp_line(linestyle="--"), label="Analytical")
ax.semilogy(t_norm, phi_l2_sim, label=f"Multiblock (N = {N}, lvl-1 inner box)", alpha=0.85)
ax.set_xlabel(r"$t / t_D$")
ax.set_ylabel(r"$\|\phi\|_{L_2}(t)$")
ax.set_title(r"Multiblock passive scalar diffusion: $L_2$-norm decay")
ax.legend()
plt.tight_layout()
plt.show()
Spatial profile across the refinement interfaces¶
The cosine should be continuous across the level-0 / level-1 interfaces at \(x = 24\) and \(x = 40\). Plot \(\phi(x)\) along a line through the refined cube (\(y = 32\), \(z = 32\)), assembled from every block the line intersects. Finer cells appear denser, but the profile should read as one smooth cosine with no visible jumps at the interfaces.
[10]:
def line_profile_concat(macr_export, t_target, y_target, z_target, scalar_name="scalar"):
"""Concatenate phi along a (y, z)-fixed line by walking every
block the line intersects. Returns (t_used, x_phys, phi, lvl)
sorted by physical x."""
blocks_by_t = parse_xdmf_blocks(macr_export.xdmf_filename)
times_avail = sorted(blocks_by_t.keys())
t_use = min(times_avail, key=lambda t: abs(t - t_target))
xs, phis, lvls = [], [], []
for block_info in blocks_by_t[t_use]:
ox, oy, oz = block_info["origin_xyz"]
dx = float(block_info["spacing_xyz"][0])
nx, ny, nz = block_info["shape_xyz"]
# Cell-centred bounds.
y_lo, y_hi = oy, oy + ny * dx
z_lo, z_hi = oz, oz + nz * dx
if not (y_lo <= y_target < y_hi and z_lo <= z_target < z_hi):
continue
iy = int((y_target - oy) / dx)
iz = int((z_target - oz) / dx)
iy = max(0, min(ny - 1, iy))
iz = max(0, min(nz - 1, iz))
phi_phys = load_block_phi_phys(block_info, t_use, scalar_name=scalar_name)
line = phi_phys[:, iy, iz]
x_phys = ox + (np.arange(nx) + 0.5) * dx
lvl_block = int(round(np.log2(1.0 / dx))) if dx < 1.0 else 0
xs.extend(x_phys.tolist())
phis.extend(line.tolist())
lvls.extend([lvl_block] * nx)
order = np.argsort(xs)
return t_use, np.array(xs)[order], np.array(phis)[order], np.array(lvls)[order]
[11]:
t_used, x_line, phi_line, lvls_line = line_profile_concat(
macr_export, t_target=float(N_STEPS), y_target=32.0, z_target=32.0
)
print(
f"Line profile at t = {t_used:.0f}: {len(x_line)} points "
f"(lvl-0: {(lvls_line == 0).sum()}, lvl-1: {(lvls_line == 1).sum()})"
)
phi_exact_line = analytical_phi_at_x(x_line, t_used)
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
ax = axes[0]
mask0 = lvls_line == 0
mask1 = lvls_line == 1
ax.plot(x_line[mask0], phi_line[mask0], "o", ms=4, label="Sim, lvl 0", alpha=0.85)
ax.plot(x_line[mask1], phi_line[mask1], "s", ms=3, label="Sim, lvl 1", alpha=0.85)
ax.plot(x_line, phi_exact_line, "k--", lw=1.2, label="Analytical")
for x_face in (24.0, 40.0):
ax.axvline(x_face, color="grey", linestyle=":", alpha=0.6)
ax.set_xlabel("x")
ax.set_ylabel(r"$\phi(x, y=32, z=32, t_{end})$")
ax.set_title("Spatial profile: simulation vs analytical")
ax.legend()
ax.grid(True, alpha=0.3)
err_line = phi_line - phi_exact_line
ax = axes[1]
ax.plot(x_line[mask0], err_line[mask0], "o", ms=4, label="lvl 0", alpha=0.85)
ax.plot(x_line[mask1], err_line[mask1], "s", ms=3, label="lvl 1", alpha=0.85)
for x_face in (24.0, 40.0):
ax.axvline(x_face, color="grey", linestyle=":", alpha=0.6)
ax.axhline(0, color="black", lw=0.5)
ax.set_xlabel("x")
ax.set_ylabel(r"$\phi^{sim} - \phi^{exact}$")
ax.set_title("Pointwise error along the line (refinement interfaces marked)")
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
print(f"Line-profile max error = {np.abs(err_line).max():.3e}")
print(f"Line-profile RMS error = {np.sqrt((err_line**2).mean()):.3e}")
Line profile at t = 5120: 80 points (lvl-0: 48, lvl-1: 32)
Line-profile max error = 4.681e-02
Line-profile RMS error = 2.933e-02
Summary¶
A passing run shows:
Decay rate. \(\lambda_{\text{lbm}}\) matches \(\lambda_{\text{exact}} = D k^2\) to within a few percent, confirming the F2C / C2F scalar pack / unpack kernels do not bias the integrated diffusion across levels.
Continuity at refinement interfaces. The \(\phi(x)\) line plot reads as one cosine with no visible jumps or oscillations at \(x = 24\) and \(x = 40\).
Volume-weighted :math:`L_2` error. Comparable in magnitude to the single-level \(N = 64\) result from Case 12a, indicating no significant error injection at the interfaces.
Version¶
[12]:
sim_info = sim_cfg.output.read_info()
nassu_commit = sim_info["commit"]
nassu_version = sim_info["version"]
print("Version:", nassu_version)
print("Commit hash:", nassu_commit)
Version: 2.0.0a7
Commit hash: 5e47f2762575c2d285254d36bfd354b6af09fda1
Configuration¶
[13]:
from IPython.display import Code
Code(filename=filename)
[13]:
# Passive scalar transport - multiblock 3-D pure-diffusion
#
# Multiblock variant of case 12 (sinusoidal scalar diffusion). Same
# physical setup - single Fourier mode `cos(k x)` diffusing on a triply
# periodic box with the fluid at rest - but the inner slab in x is
# refined to level 1 so the cosine wave passes through two refinement
# interfaces twice per period. This exercises the scalar same-level
# comm at every block face *and* the F2C / C2F pack-and-unpack path at
# the refinement interfaces continuously.
#
# Initial scalar field (cosine mode along x, base level):
# phi(x, y, z, 0) = cos(2*pi*x/N), N = 64
# k = 2*pi/N
#
# Analytical solution (pure diffusion, physical units):
# phi(x, y, z, t) = cos(k x) * exp(-D k^2 t)
#
# The decay rate `lambda = D k^2` is independent of how the domain is
# discretised, so a multiblock run must reproduce the same
# `lambda_lbm` as the matching single-level grid, and the cosine
# profile must remain continuous across the refinement interfaces.
#
# Error metric: fitted decay rate from the L2 norm time series, plus
# a per-block max pointwise error against the analytical solution.
simulations:
- name: passiveScalarDiffusion3DMultiblock
save_path: ./validation/scalar_transport/01_passive_scalar_transport/results/passive_scalar_diffusion_multiblock
# Same scaling as the N=64 single-level run in case 12 so the
# results are directly comparable: D = 1e-3 in lattice units at
# level 0, n_steps = 5120 (one e-folding factor 0.95 of the
# cosine amplitude).
n_steps: 5120
report:
frequency: 1000
domain:
domain_size:
x: 64
y: 64
z: 64
block_size: 8
refinement:
static:
default:
volumes_refine:
# Inner box centred in the domain, refined to level 1.
# With periodic BCs the cosine wave crosses the two x
# refinement interfaces (x = 24 and x = 40) twice per
# period, so F2C / C2F runs at every step. Box-aligned
# start / end values (multiples of `block_size = 8`)
# keep the level-1 region commensurate with the level-0
# grid.
- start: [24, 24, 24]
end: [40, 40, 40]
lvl: 1
is_abs: true
data:
exports:
default:
macrs: [rho, u, scalar_phi, scalar_q_neq]
interval:
frequency: 256
lvl: 0
target:
volumes:
default: {}
outputs:
instantaneous: true
plane_series:
macrs: [scalar_phi]
interval: {frequency: 256, lvl: 0}
target:
planes:
# Mid-domain plane through the refined slab: shows the cosine
# mode crossing both refinement interfaces.
mid_y:
axis: y
axis_pos: 32
dist: 0.5
outputs:
instantaneous: true
models:
precision:
default: single
LBM:
# Fluid stays at the equilibrium fixed point (u = 0, rho = 1);
# tau only damps numerical noise. Matches the case-12 single-
# level setup so multiblock and single-level results are
# directly comparable.
tau: 0.6
vel_set: D3Q27
coll_oper: RRBGK
engine:
name: CUDA
BC:
periodic_dims: [true, true, true]
scalar_transports:
scalar:
velocity_set: D3Q7
collision_operator: RRBGK
adv_diff_equation:
# D in lattice units at level 0. The level-rescaling
# `D_phi(lvl) = D_phi_0 * 2^lvl` keeps the physical
# diffusivity the same on both levels.
D: 1.0e-3
S: "0"
initial_field: "cos(2 * pi * x / 64)"