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

\[\phi(x, y, z, 0) = \cos(k x), \quad k = 2\pi / 64,\]

with the analytical decay (eigenfunction of the Laplacian, fluid at rest)

\[\phi(x, y, z, t) = \cos(k x) \, e^{-D k^2 t}.\]

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

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

  2. The cosine spatial profile remains continuous across the refinement interfaces: no jumps or oscillations near \(x = 24\) or \(x = 40\).

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

  1. transpose the dataset with phi_disk.transpose(2, 1, 0); after that step arr_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()
../../../_images/validation_scalar_transport_01_passive_scalar_transport_01b_passive_scalar_diffusion_multiblock_16_0.png

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)
../../../_images/validation_scalar_transport_01_passive_scalar_transport_01b_passive_scalar_diffusion_multiblock_19_1.png
Line-profile max error  = 4.681e-02
Line-profile RMS error  = 2.933e-02

Summary

A passing run shows:

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

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

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