2-D Taylor-Green Vortex - Spatial Convergence (Case 00a)

This notebook validates second-order spatial accuracy of the RR-BGK LBM against the exact analytical 2-D Taylor-Green vortex solution.

The initial condition on a periodic domain \([0, N)^2\) is:

\[u_x(x,y,0) = U_0 \cos(kx)\sin(ky), \quad u_y(x,y,0) = -U_0 \sin(kx)\cos(ky)\]

where \(k = 2\pi/N\). The nonlinear advection terms vanish for this mode, so the exact time-evolved solution is:

\[u_x(x,y,t) = U_0 \cos(kx)\sin(ky)\,e^{-2\nu k^2 t}, \quad u_y(x,y,t) = -U_0 \sin(kx)\cos(ky)\,e^{-2\nu k^2 t}\]

Protocol: Fix \(Re = U_0 N / \nu = 240\) and \(U_0 = 0.01\) across all grids. Derive \(\nu = U_0 N / Re\) and \(\tau = 0.5 + 3\nu\) per level.

Four grids (\(N = 8, 16, 32, 64\)) are compared. The \(L_2\) velocity error should decrease as \(O(\Delta x^2)\) on a log-log plot.

Setup

import numpy as np
import matplotlib.pyplot as plt
import pathlib
from nassu.cfg.model import ConfigScheme
from tests.validation.notebooks import common

common.use_style()

Load simulation configuration

filename = "tests/validation/cases/00a_taylor_green_vortex_2d.nassu.yaml"
sim_cfgs = ConfigScheme.sim_cfgs_from_file_dct(filename)

# Unrolled simulations: sim_id 0..3 correspond to N = 32, 64, 128, 256
sim_cfgs_list = [sim_cfgs["taylorGreenVortex2D", i] for i in range(4)]

# Protocol constants
U0 = 0.01  # peak velocity amplitude (lattice units)
CS2 = 1.0 / 3.0  # speed of sound squared

# Derive quantities from config
GRID_SIZES = [sim_cfg.domain.domain_size.x for sim_cfg in sim_cfgs_list]
NU = {cfg.domain.domain_size.x: cfg.models.LBM.kinematic_viscosity for cfg in sim_cfgs_list}
TAU = {cfg.domain.domain_size.x: cfg.models.LBM.tau for cfg in sim_cfgs_list}
N_STEPS = {cfg.domain.domain_size.x: cfg.n_steps for cfg in sim_cfgs_list}

for N in GRID_SIZES:
    Ma = U0 * np.sqrt(3)
    Re = U0 * N / NU[N]
    print(
        f"N={N:3d}: nu={NU[N]:.5f}, tau={TAU[N]:.4f}, "
        f"steps={N_STEPS[N]}, Re={Re:.0f}, Ma={Ma:.3f}"
    )

N=  8: nu=0.00033, tau=0.5010, steps=1280, Re=240, Ma=0.017
N= 16: nu=0.00067, tau=0.5020, steps=2560, Re=240, Ma=0.017
N= 32: nu=0.00133, tau=0.5040, steps=5120, Re=240, Ma=0.017
N= 64: nu=0.00267, tau=0.5080, steps=10240, Re=240, Ma=0.017

Analytical solution functions

The exact solution is valid for all \(t > 0\) because the nonlinear advection terms are identically zero for the single-Fourier-mode initial condition.

def analytical_tgv(N, U0, nu, n_steps):
    """Exact 2-D TGV velocity field at step n_steps.

    Positions are at integer lattice sites x = 0, 1, ..., N-1.
    Convention: ux = U0*cos(kx)*sin(ky)*decay, uy = -U0*sin(kx)*cos(ky)*decay.

    Returns:
        ux, uy: arrays of shape (N, N) with indexing='ij' (axis 0 = x, axis 1 = y).
    """
    k = 2.0 * np.pi / (N)
    x, y = np.meshgrid(np.arange(N), np.arange(N), indexing="ij")

    decay = np.exp(-2.0 * nu * k**2 * (n_steps + 1))
    ux = U0 * np.cos(k * x) * np.sin(k * y) * decay
    uy = -U0 * np.sin(k * x) * np.cos(k * y) * decay
    return ux, uy


def l2_error(ux_lbm, uy_lbm, ux_exact, uy_exact, U0):
    """Normalised L2 velocity error.

    E_L2 = (1/U0) * sqrt( mean( (ux_lbm - ux_exact)^2 + (uy_lbm - uy_exact)^2 ) )
    """
    diff2 = (ux_lbm - ux_exact) ** 2 + (uy_lbm - uy_exact) ** 2
    return np.sqrt(diff2).mean() / U0


def E_analytical(t, U0, N, nu):
    """Exact kinetic energy per unit area at step t.

    E(t) = (U0^2 / 4) * exp(-4 * nu * k^2 * t)
    """
    k = 2.0 * np.pi / N
    return (U0**2 / 4.0) * np.exp(-4.0 * nu * k**2 * t)

HDF5/XDMF output loader

Nassu writes snapshots to rolling HDF5 files indexed by an XDMF manifest. The functions below read that format directly - no VTK dependency required.

import h5py
from lxml import etree
from nassu.cfg.schemes.simul import MacrExport


def load_h5_kinetic_energy(macr_export: MacrExport, dim=2):
    """Compute volume-averaged kinetic energy E(t) from XDMF/HDF5 output.

    Reads all rolling HDF5 files under the macrs/instantaneous directory.
    Returns (times, E) as float arrays sorted by time.
    """
    macrs_dir = macr_export.xdmf_filename.parent
    export_name = macr_export.base_path.name
    h5_files = sorted(macrs_dir.glob(f"{export_name}.*.h5"))

    times = []
    energies = []

    for h5_path in h5_files:
        with h5py.File(h5_path, "r") as f:
            for ts_key in sorted(f.keys()):
                t = float(ts_key[1:])  # strip leading "t"
                total_E = 0.0
                total_cells = 0
                for block_key in f[ts_key].keys():
                    rho = f[ts_key][block_key]["rho"][()]
                    ux = f[ts_key][block_key]["ux"][()]
                    uy = f[ts_key][block_key]["uy"][()]
                    if dim == 3 and "uz" in f[ts_key][block_key]:
                        uz = f[ts_key][block_key]["uz"][()]
                        ke = 0.5 * rho * (ux**2 + uy**2 + uz**2)
                    else:
                        ke = 0.5 * rho * (ux**2 + uy**2)
                    total_E += float(ke.sum())
                    total_cells += ke.size
                times.append(t)
                energies.append(total_E / total_cells)

    idx = np.argsort(times)
    return np.array(times)[idx], np.array(energies)[idx]


def _parse_xdmf_blocks(xdmf_path):
    """Parse XDMF XML and return per-timestep block metadata.

    Returns dict: t (float) -> list of dicts with keys:
      block_key, h5_path, origin_xyz, spacing_xyz, shape_zyx (cell counts).
    """
    xdmf_path = pathlib.Path(xdmf_path)
    tree = etree.parse(str(xdmf_path))
    root = tree.getroot()
    domain = root.find("Domain")
    time_series = domain.find("Grid")  # TimeSeries collection

    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_zyx = np.array(geom.find("DataItem[@Name='Origin']").text.split(), dtype=float)
            spacing_zyx = np.array(
                geom.find("DataItem[@Name='Spacing']").text.split(), dtype=float
            )
            topo = block_grid.find("Topology")
            shape_node_zyx = tuple(int(d) for d in topo.get("Dimensions").split())
            shape_zyx = tuple(d - 1 for d in shape_node_zyx)

            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_zyx[::-1],
                    "spacing_xyz": spacing_zyx[::-1],
                    "shape_zyx": shape_zyx,
                }
            )
        result[t] = blocks
    return result


def assemble_2d_field(macr_export: MacrExport, t_target, macr_names, N=32):
    """Assemble a full 2-D macroscopic field at the snapshot closest to t_target.

    Returns:
        fields: dict {macr_name: ndarray shape (N, N), axes (x, y)}
        t_used: float, actual snapshot time used
        spacing: float, grid spacing dx
    """
    xdmf_path = macr_export.xdmf_filename
    blocks_by_t = _parse_xdmf_blocks(xdmf_path)

    times = sorted(blocks_by_t.keys())
    t_use = min(times, key=lambda t: abs(t - t_target))

    full = {m: np.zeros((N, N), dtype=np.float32) for m in macr_names}
    spacing = None

    for block_info in blocks_by_t[t_use]:
        h5_path = block_info["h5_path"]
        block_key = block_info["block_key"]

        nz, ny, nx = block_info["shape_zyx"]

        ox, oy = block_info["origin_xyz"][:2]
        dx = block_info["spacing_xyz"][0]
        spacing = dx

        ix = int(round(ox / dx))
        iy = int(round(oy / dx))

        ts_key = f"t{t_use:.6f}"
        with h5py.File(h5_path, "r") as f:
            for m in macr_names:
                if m in f[ts_key][block_key]:
                    arr = f[ts_key][block_key][m][()]
                    arr_2d = arr[:, :, 0].T
                    full[m][ix : ix + nx, iy : iy + ny] = arr_2d

    return full, t_use, spacing

Load simulation output

Kinetic energy is computed directly from the rolling HDF5 files without assembling full fields - this is fast even for large outputs. Full velocity fields are loaded separately for the L2 error analysis.

data = {}
for sim_cfg in sim_cfgs_list:
    N = sim_cfg.domain.domain_size.x
    macr_export = sim_cfg.output.instantaneous["default"]

    if not macr_export.xdmf_filename.exists():
        print(f"N={N}: XDMF output not found ({macr_export.xdmf_filename}) - skipping.")
        continue

    # Load kinetic energy time series (for energy decay plot).
    times, E = load_h5_kinetic_energy(macr_export, dim=2)

    # Load final velocity field for L2 error.
    t_end = float(N_STEPS[N])
    fields, t_used, dx = assemble_2d_field(macr_export, t_end, ["ux", "uy"], N=N)

    data[N] = {
        "times": times,
        "E": E,
        "ux": fields["ux"],
        "uy": fields["uy"],
        "t_used": t_used,
    }
    print(
        f"N={N} (sim_id={sim_cfg.sim_id:03d}): {len(times)} snapshots, "
        f"velocity field loaded at t={t_used:.0f}"
    )

N=8 (sim_id=000): 81 snapshots, velocity field loaded at t=1280
N=16 (sim_id=001): 81 snapshots, velocity field loaded at t=2560
N=32 (sim_id=002): 81 snapshots, velocity field loaded at t=5120
N=64 (sim_id=003): 81 snapshots, velocity field loaded at t=10240

\(L_2\) velocity error at \(t^* = 1\)

The normalised \(L_2\) error is:

\[E_{L2} = \frac{1}{U_0}\sqrt{\frac{1}{N^2}\sum_{i,j}\left[(u_x^{\text{lbm}} - u_x^{\text{exact}})^2 + (u_y^{\text{lbm}} - u_y^{\text{exact}})^2\right]}\]

The analytical field is evaluated at integer lattice positions \(x,y = 0, 1, \ldots, N{-}1\).

# Compute L2 velocity error at t*_end for each grid.
errors = []
N_values = []

for N in GRID_SIZES:
    if N not in data:
        continue
    nu = NU[N]
    n = N_STEPS[N]

    # Exact solution at integer lattice positions.
    ux_exact, uy_exact = analytical_tgv(N, U0, nu, n)

    # Simulation output (shape (N, N), axes (x, y)).
    ux_sim = data[N]["ux"].astype(np.float64)
    uy_sim = data[N]["uy"].astype(np.float64)

    err = l2_error(ux_sim, uy_sim, ux_exact, uy_exact, U0)
    errors.append(err)
    N_values.append(N)
    print(f"N={N:3d}: L2 error = {err:.4e}")

errors = np.array(errors)
N_values = np.array(N_values)

N=  8: L2 error = 3.9538e-02
N= 16: L2 error = 1.0122e-02
N= 32: L2 error = 2.4922e-03
N= 64: L2 error = 7.0099e-04

Diagnostic: \(u_x\) field comparison (simulation vs analytical)

Side-by-side contour plots of the simulated and analytical \(u_x\) fields at \(t^* = 1\), plus the pointwise error field. This helps diagnose axis or assembly issues that would not be visible in the scalar \(L_2\) metric alone.

for N_diag in GRID_SIZES:
    nu = NU[N_diag]
    n = N_STEPS[N_diag]
    ux_exact, uy_exact = analytical_tgv(N_diag, U0, nu, n)
    ux_sim = data[N_diag]["ux"].astype(np.float64)

    vmin, vmax = ux_exact.min(), ux_exact.max()
    err_field = ux_sim - ux_exact

    fig, axes = plt.subplots(1, 3, figsize=(16, 5))

    # Simulation ux
    im0 = axes[0].imshow(ux_sim.T, origin="lower", cmap="RdBu_r", vmin=vmin, vmax=vmax)
    axes[0].set_title(f"Simulation $u_x$ (N={N_diag})")
    axes[0].set_xlabel("x")
    axes[0].set_ylabel("y")
    axes[0].grid(False)
    fig.colorbar(im0, ax=axes[0], shrink=0.8)

    # Analytical ux
    im1 = axes[1].imshow(ux_exact.T, origin="lower", cmap="RdBu_r", vmin=vmin, vmax=vmax)
    axes[1].set_title(f"Analytical $u_x$ (N={N_diag})")
    axes[1].set_xlabel("x")
    axes[1].set_ylabel("y")
    axes[1].grid(False)
    fig.colorbar(im1, ax=axes[1], shrink=0.8)

    # Pointwise error
    err_max = np.abs(err_field).max()
    im2 = axes[2].imshow(err_field.T, origin="lower", cmap="RdBu_r", vmin=-err_max, vmax=err_max)
    axes[2].set_title("Error $u_x^{\\mathrm{sim}} - u_x^{\\mathrm{exact}}$")
    axes[2].set_xlabel("x")
    axes[2].set_ylabel("y")
    axes[2].grid(False)
    fig.colorbar(im2, ax=axes[2], shrink=0.8)

    plt.tight_layout()
    plt.show()

    # Print key diagnostics
    print(f"Simulation ux range: [{ux_sim.min():.6f}, {ux_sim.max():.6f}]")
    print(f"Analytical ux range: [{ux_exact.min():.6f}, {ux_exact.max():.6f}]")
    print(f"Max pointwise error: {np.abs(err_field).max():.4e}")
    print(f"Decay factor at t={n}: {np.exp(-2.0 * nu * (2*np.pi/N_diag)**2 * n):.6f}")

Simulation ux range: [-0.005293, 0.005293]
Analytical ux range: [-0.005905, 0.005905]
Max pointwise error: 6.1238e-04
Decay factor at t=1280: 0.590740

Simulation ux range: [-0.005757, 0.005757]
Analytical ux range: [-0.005906, 0.005906]
Max pointwise error: 1.4963e-04
Decay factor at t=2560: 0.590740

Simulation ux range: [-0.005871, 0.005870]
Analytical ux range: [-0.005907, 0.005907]
Max pointwise error: 3.7315e-05
Decay factor at t=5120: 0.590740

Simulation ux range: [-0.005899, 0.005900]
Analytical ux range: [-0.005907, 0.005907]
Max pointwise error: 1.1415e-05
Decay factor at t=10240: 0.590740

fig, ax = common.fig_single()

if len(N_values) > 0:
    ax.loglog(
        N_values, errors, **common.markers.sim(shape="o", linestyle="-"), label="AeroSim RR-BGK"
    )

    # O(N^{-2}) reference line (equivalent to O(dx^2)).
    N_ref = np.array([N_values.min() * 0.7, N_values.max() * 1.3])
    scale = errors[0] * N_values[0] ** 2
    ax.loglog(
        N_ref,
        scale * N_ref ** (-2),
        **common.markers.exp_line(linestyle="--"),
        label=r"$O(N^{-2})$",
    )

Re = U0 * GRID_SIZES[0] / NU[GRID_SIZES[0]]
ax.set_xlabel(r"$N$")
ax.set_ylabel(r"$E_{L2}$")
ax.set_xticks([8, 16, 32, 64])
ax.set_xticklabels(["8", "16", "32", "64"])
ax.set_title(f"2-D TGV spatial convergence ($Re = {Re:.0f}$)")
ax.legend()
plt.tight_layout()
plt.show()

Kinetic energy decay (supplementary)

The volume-averaged kinetic energy \(E(t)\) should follow the exact exponential decay \(E_0 \exp(-4\nu k^2 t)\). This plot provides a qualitative check that the viscous dissipation rate is correct at each resolution.

fig, ax = plt.subplots(figsize=(8, 5))

for i, N in enumerate(GRID_SIZES[::-1]):
    if N not in data:
        continue
    nu = NU[N]
    t_arr = data[N]["times"]
    E_sim = data[N]["E"]
    # Normalise time by viscous decay scale t_nu = 1 / (4 * nu * k^2).
    k = 2.0 * np.pi / N
    t_decay = 1.0 / (4.0 * nu * k**2)
    t_norm = t_arr / t_decay
    E_exact = np.array([E_analytical(t, U0, N, nu) for t in t_arr])
    if i == 0:
        ax.semilogy(t_norm, E_exact, **common.markers.exp_line(linestyle="--"), label="Analytical")
    ax.semilogy(t_norm, E_sim, label=f"N = {N}", alpha=0.7)

ax.set_xlabel(r"$t / t_\nu$")
ax.set_ylabel(r"$E(t)$")
ax.set_title("2-D TGV: kinetic energy decay")
ax.legend()
plt.tight_layout()
plt.show()

Summary

A passing result shows:

1. The \(L_2\) velocity error decreases with grid refinement, with all errors well below 1%.

2. The convergence order between consecutive levels is close to 2, confirming second-order spatial accuracy of the RR-BGK LBM.

3. The kinetic energy \(E(t)\) tracks the exact analytical decay curve for all grid sizes.

The \(L_2\) error is normalised by \(U_0\) so that it represents the relative velocity error independent of the Mach number.

Version

sim_info = sim_cfgs_list[0].output.read_info()

nassu_commit = sim_info["commit"]
nassu_version = sim_info["version"]
print("Version:", nassu_version)
print("Commit hash:", nassu_commit)

Version: 1.6.45
Commit hash: 0262a0dd797d4c1ad8d58db0226f9706a87f359e

Configuration

from IPython.display import Code

Code(filename=filename)

# Taylor-Green Vortex - 2D spatial convergence
#
# Validates second-order spatial accuracy of the RR-BGK LBM against the exact
# analytical 2D TGV solution at Re = 240. Four grid resolutions (N = 8, 16,
# 32, 64) are compared at the same dimensionless end time.
#
# Protocol (constant Re and U_lbm across grids):
#   Re = U0 * N / nu = 240
#   U0 = 0.01, tau in [0.501, 0.508]
#   Ma = U0 * sqrt(3) = 0.017
#
# Initial condition (cos*sin convention):
#   ux =  U0 * cos(kx) * sin(ky)
#   uy = -U0 * sin(kx) * cos(ky)
#   k  = 2*pi/N
#
# Error metric: L2 velocity error at t_end normalised by U0.
# Expected convergence: O(dx^2) on a log-log plot.
#
# Note: requires models.initialization.equations (PR #443 / feat/init-equation-sem-field).

variables:
  U0: 0.01

simulations:
  - name: taylorGreenVortex2D
    save_path: ./tests/validation/results/00a_taylor_green_vortex_2d

    n_steps: !unroll [1280, 2560, 5120, 10240]

    report:
      frequency: 1000

    domain:
      domain_size:
        x: !unroll [8, 16, 32, 64]
        y: !unroll [8, 16, 32, 64]
      block_size: 8

    data:
      divergence: { frequency: 100 }
      instantaneous:
        default:
          # ~20 snapshots per case (frequency = N)
          interval: { frequency: !unroll [16, 32, 64, 128] }
          macrs: [rho, u]
      statistics:
        interval: { frequency: 0 }

    models:
      precision:
        default: single

      LBM:
        tau: !unroll [0.501, 0.502, 0.504, 0.508]
        vel_set: D2Q9
        coll_oper: RRBGK

      engine:
        name: CUDA

      BC:
        periodic_dims: [true, true]

      initialization:
        equations:
          rho: "1"
          ux: !unroll
            - "0.01 * cos(2 * pi * x / 8) * sin(2 * pi * y / 8)"
            - "0.01 * cos(2 * pi * x / 16) * sin(2 * pi * y / 16)"
            - "0.01 * cos(2 * pi * x / 32) * sin(2 * pi * y / 32)"
            - "0.01 * cos(2 * pi * x / 64) * sin(2 * pi * y / 64)"
          uy: !unroll
            - "-0.01 * sin(2 * pi * x / 8) * cos(2 * pi * y / 8)"
            - "-0.01 * sin(2 * pi * x / 16) * cos(2 * pi * y / 16)"
            - "-0.01 * sin(2 * pi * x / 32) * cos(2 * pi * y / 32)"
            - "-0.01 * sin(2 * pi * x / 64) * cos(2 * pi * y / 64)"
          uz: "0"