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 matplotlib.pyplot as plt
import matplotlib.ticker as mticker
import numpy as np

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

common.use_style()

Load simulation configuration

filename = "validation/analytical/04_taylor_green_vortex_2d/04_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}, 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)

Field loaders

Fields are read through nassu.viz.FieldSource, the shared reader for Nassu’s XDMF/HDF5 macroscopic exports. Its read_arrays method assembles a snapshot into numpy arrays indexed (x, y, z), so no manual HDF5 axis handling is needed here - the helpers below just reduce those arrays to the kinetic energy time series and the final velocity field.

def load_kinetic_energy(src):
    """Volume-averaged kinetic energy E(t) over every exported step.

    E(t) = mean( 0.5 * rho * (ux^2 + uy^2) ) over all cells.
    Returns (times, E) as float arrays sorted by time.
    """
    times = []
    energies = []
    for step in src.steps:
        arrays, _ = src.read_arrays(step, ["rho", "ux", "uy"])
        ke = 0.5 * arrays["rho"] * (arrays["ux"] ** 2 + arrays["uy"] ** 2)
        times.append(src.steps_to_time[step])
        energies.append(float(np.nanmean(ke)))

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


def load_velocity_field(src, step):
    """Return (ux, uy) as (N, N) arrays indexed (x, y) at ``step``."""
    arrays, _ = src.read_arrays(step, ["ux", "uy"])
    return arrays["ux"][:, :, 0], arrays["uy"][:, :, 0]

Load simulation output

Each simulation is opened with FieldSource. The kinetic energy time series drives the energy-decay plot, and the final-step velocity field feeds the \(L_2\) error analysis.

PROJECT_ROOT = common.find_project_root()

data = {}
for sim_cfg in sim_cfgs_list:
    N = sim_cfg.domain.domain_size.x

    try:
        src = common.FieldSource.from_cfg(sim_cfg, project_root=PROJECT_ROOT)
    except FileNotFoundError as exc:
        print(f"N={N}: {exc} - skipping.")
        continue

    # Kinetic energy time series (for the energy-decay plot).
    times, E = load_kinetic_energy(src)

    # Final velocity field (for the L2 error).
    step_end = src.steps[-1]
    ux, uy = load_velocity_field(src, step_end)
    t_used = src.steps_to_time[step_end]

    data[N] = {
        "times": times,
        "E": E,
        "ux": ux,
        "uy": 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.9536e-02
N= 16: L2 error = 1.0126e-02
N= 32: L2 error = 2.4894e-03
N= 64: L2 error = 6.9372e-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.1245e-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.4951e-04
Decay factor at t=2560: 0.590740

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

Simulation ux range: [-0.005899, 0.005899]
Analytical ux range: [-0.005907, 0.005907]
Max pointwise error: 1.1408e-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}$")
# Label only the grid sizes on the log x-axis. The default LogFormatter also
# labels the minor ticks (6x10^0, 2x10^1, ...), which overlap these and produce
# a garbled axis - suppress them.
ax.set_xticks([8, 16, 32, 64])
ax.set_xticklabels(["8", "16", "32", "64"])
ax.xaxis.set_minor_formatter(mticker.NullFormatter())
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: 2.0.0a7
Commit hash: 5e47f2762575c2d285254d36bfd354b6af09fda1

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: ./validation/analytical/04_taylor_green_vortex_2d/results

    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:
      exports:
        default:
          macrs: [rho, u]
          interval:
            frequency: !unroll [16, 32, 64, 128]
            lvl: 0
          target:
            volumes:
              default: {}
          outputs:
            instantaneous: true

    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"