Lid-Driven Cavity - Moment-Based BC Validation (Case 05)

Validation of the moment-based boundary conditions against the Ghia, Ghia & Shin (1982) benchmark for the square lid-driven cavity. The lid is the moment-based moving wall RegularizedVelocityWall; the three stationary walls are the no-slip RegularizedHWBB (the zero-velocity case of the same closure). No population bounce-back is used.

The steady-state centreline velocity profiles are compared against the tabulated Ghia solution at Re = 100, 400, 1000.

Setup

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

filename = "validation/analytical/05_lid_driven_cavity/05_lid_driven_cavity.nassu.yaml"
sim_cfgs = ConfigScheme.sim_cfgs_from_file_dct(filename)

# Unrolled simulations: sim_id 0/1/2 correspond to Re = 100, 400, 1000.
sim_cfgs_list = [sim_cfgs["lidDrivenCavity", i] for i in range(3)]

RE_NUMBERS = [100, 400, 1000]
# Lid speed (lattice units): the `ux` of the moving-wall (N) BC entry.
lid_bc = next(b for b in sim_cfgs_list[0].models.BC.BC_map if b.BC == "RegularizedVelocityWall")
U_LID = float(lid_bc.params["ux"])

for Re, cfg in zip(RE_NUMBERS, sim_cfgs_list):
    L = cfg.domain.domain_size.x
    nu = (cfg.models.LBM.tau - 0.5) / 3.0
    Re_check = U_LID * L / nu
    print(
        f"Re={Re:4d} (sim_id={cfg.sim_id}): L={L}, tau={cfg.models.LBM.tau}, "
        f"nu={nu:.5f}, U={U_LID}, Re=U L/nu={Re_check:.1f}, n_steps={cfg.n_steps}"
    )

Re= 100 (sim_id=0): L=256, tau=0.884, nu=0.12800, U=0.05, Re=U L/nu=100.0, n_steps=300000
Re= 400 (sim_id=1): L=256, tau=0.596, nu=0.03200, U=0.05, Re=U L/nu=400.0, n_steps=500000
Re=1000 (sim_id=2): L=256, tau=0.5384, nu=0.01280, U=0.05, Re=U L/nu=1000.0, n_steps=800000

Reference solution

The accepted benchmark is the multigrid Navier-Stokes solution of Ghia, Ghia & Shin (1982), tabulated as the streamwise velocity u along the vertical centreline (x = L/2) and the vertical velocity v along the horizontal centreline (y = L/2), both normalised by the lid speed. The digitised tables are committed under reference/.

REF_DIR = "validation/analytical/05_lid_driven_cavity/reference"
ghia_u = np.genfromtxt(f"{REF_DIR}/ghia1982_u_vertical_centerline.csv", delimiter=",", names=True)
ghia_v = np.genfromtxt(f"{REF_DIR}/ghia1982_v_horizontal_centerline.csv", delimiter=",", names=True)

# Column name per Reynolds number.
ghia_u_col = {100: "u_Re100", 400: "u_Re400", 1000: "u_Re1000"}
ghia_v_col = {100: "v_Re100", 400: "v_Re400", 1000: "v_Re1000"}
print("Ghia u-table rows:", ghia_u["y"].size, "  v-table rows:", ghia_v["x"].size)

Ghia u-table rows: 17   v-table rows: 17

Centreline extraction

FieldSource returns each macroscopic field in Fortran (x, y) order. The lid sits at y = L (the last y index), so:

• the vertical-centreline profile u(y) at x = L/2 is ux[mid_x, :], indexed by y;

• the horizontal-centreline profile v(x) at y = L/2 is uy[:, mid_y], indexed by x.

The two central columns/rows are averaged so the line sits exactly at 0.5. Coordinates are normalised as i / (N - 1) so node 0 is the wall and node N-1 is the opposite wall (node-coincident regularized walls).

def _mid_indices(N):
    return [N // 2 - 1, N // 2]


def centerlines(src, step, U):
    """Return (y_norm, u_over_U, x_norm, v_over_U) at a given exported step."""
    arrays, _ = src.read_arrays(step, ["ux", "uy"])
    ux = np.asarray(arrays["ux"]).squeeze()  # (x, y)
    uy = np.asarray(arrays["uy"]).squeeze()
    N = ux.shape[0]
    mid = _mid_indices(N)
    coord = np.arange(N) / (N - 1)
    u_y = ux[mid, :].mean(axis=0) / U   # u(y) at x = L/2
    v_x = uy[:, mid].mean(axis=1) / U   # v(x) at y = L/2
    return coord, u_y, coord, v_x

Load simulation output

PROJECT_ROOT = common.find_project_root()

data = {}
for Re, cfg in zip(RE_NUMBERS, sim_cfgs_list):
    src = common.FieldSource.from_cfg(cfg, project_root=PROJECT_ROOT)
    step_end = src.steps[-1]
    y, u_y, x, v_x = centerlines(src, step_end, U_LID)

    # Steady-state check: how much the centreline moved over the last interval.
    _, u_prev, _, _ = centerlines(src, src.steps[-2], U_LID)
    steady_resid = float(np.max(np.abs(u_y - u_prev)))

    data[Re] = dict(y=y, u_y=u_y, x=x, v_x=v_x, steady_resid=steady_resid)
    print(
        f"Re={Re:4d}: {len(src.steps)} snapshots, last step {step_end}, "
        f"steady residual max|du/U| over last interval = {steady_resid:.2e}"
    )

Re= 100: 21 snapshots, last step 300000, steady residual max|du/U| over last interval = 6.41e-07
Re= 400: 21 snapshots, last step 500000, steady residual max|du/U| over last interval = 1.51e-06
Re=1000: 21 snapshots, last step 800000, steady residual max|du/U| over last interval = 1.91e-06

\(u\) along the vertical centreline

Streamwise velocity u(y) at x = L/2, simulation (lines) vs Ghia (markers).

fig, ax = common.fig_single()
colors = plt.rcParams["axes.prop_cycle"].by_key()["color"]
for c, Re in zip(colors, RE_NUMBERS):
    d = data[Re]
    ax.plot(d["u_y"], d["y"], **{**common.markers.sim_line(), "color": c}, label=f"Nassu Re={Re}")
    ax.plot(
        ghia_u[ghia_u_col[Re]], ghia_u["y"],
        **{**common.markers.exp(), "color": c}, label=f"Ghia Re={Re}",
    )
ax.set_xlabel(r"$u / U$")
ax.set_ylabel(r"$y / L$")
ax.set_title("Vertical-centreline streamwise velocity")
ax.legend(fontsize=8, ncol=2)
fig.tight_layout()

\(v\) along the horizontal centreline

Vertical velocity v(x) at y = L/2, simulation (lines) vs Ghia (markers).

fig, ax = common.fig_single()
for c, Re in zip(colors, RE_NUMBERS):
    d = data[Re]
    ax.plot(d["x"], d["v_x"], **{**common.markers.sim_line(), "color": c}, label=f"Nassu Re={Re}")
    ax.plot(
        ghia_v["x"], ghia_v[ghia_v_col[Re]],
        **{**common.markers.exp(), "color": c}, label=f"Ghia Re={Re}",
    )
ax.set_xlabel(r"$x / L$")
ax.set_ylabel(r"$v / U$")
ax.set_title("Horizontal-centreline vertical velocity")
ax.legend(fontsize=8, ncol=2)
fig.tight_layout()

Quantitative error vs Ghia

Mean absolute deviation (normalised by the lid speed) between the simulated centreline profiles, interpolated to the Ghia node positions, and the tabulated values.

print(f"{'Re':>5}  {'MAE(u)/U':>10}  {'MAE(v)/U':>10}")
mae_u, mae_v = [], []
for Re in RE_NUMBERS:
    d = data[Re]
    u_at = np.interp(ghia_u["y"], d["y"], d["u_y"])
    v_at = np.interp(ghia_v["x"], d["x"], d["v_x"])
    eu = float(np.mean(np.abs(u_at - ghia_u[ghia_u_col[Re]])))
    ev = float(np.mean(np.abs(v_at - ghia_v[ghia_v_col[Re]])))
    mae_u.append(eu)
    mae_v.append(ev)
    print(f"{Re:>5}  {eu:>10.4f}  {ev:>10.4f}")

fig, ax = common.fig_single()
xpos = np.arange(len(RE_NUMBERS))
ax.bar(xpos - 0.2, mae_u, 0.4, label=r"MAE$(u)/U$")
ax.bar(xpos + 0.2, mae_v, 0.4, label=r"MAE$(v)/U$")
ax.set_xticks(xpos)
ax.set_xticklabels([f"Re={Re}" for Re in RE_NUMBERS])
ax.set_ylabel("mean absolute deviation / U")
ax.set_title("Centreline deviation from Ghia (1982)")
ax.legend()
fig.tight_layout()

   Re    MAE(u)/U    MAE(v)/U
  100      0.0017      0.0020
  400      0.0027      0.0120
 1000      0.0044      0.0049

Flow field: primary and secondary vortices

Streamlines of the steady Re = 1000 solution, showing the primary vortex and the secondary corner vortices.

cfg = sim_cfgs_list[-1]  # Re = 1000
src = common.FieldSource.from_cfg(cfg, project_root=PROJECT_ROOT)
arrays, _ = src.read_arrays(src.steps[-1], ["ux", "uy"])
ux = np.asarray(arrays["ux"]).squeeze()
uy = np.asarray(arrays["uy"]).squeeze()
N = ux.shape[0]
xs = np.arange(N) / (N - 1)
X, Y = np.meshgrid(xs, xs, indexing="ij")
speed = np.sqrt(ux**2 + uy**2) / U_LID

fig, ax = plt.subplots(figsize=(5.2, 5))
strm = ax.streamplot(
    X.T, Y.T, ux.T, uy.T, density=1.6, color=speed.T, cmap="viridis", linewidth=0.7,
)
fig.colorbar(strm.lines, ax=ax, label=r"$|u| / U$")
ax.set_xlabel(r"$x / L$")
ax.set_ylabel(r"$y / L$")
ax.set_title("Lid-driven cavity streamlines, Re = 1000")
ax.set_aspect("equal")
fig.tight_layout()

Summary

Success criteria:

• The simulated centreline profiles overlay the Ghia (1982) benchmark at every Reynolds number, reproducing the primary-vortex sign reversal and the near-wall velocity extrema.

• The mean absolute deviation from the tabulated points stays within a few percent of the lid speed, decreasing with Reynolds number resolution.

• The steady-state residual confirms the flow has converged (the centreline stops moving).

• The moment-based moving wall imposes u = U exactly at the lid and u = 0 at the no-slip walls, with the density staying close to 1 (weakly compressible).

Version

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

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

Version: 2.0.0a7
Commit hash: 5e47f2762575c2d285254d36bfd354b6af09fda1

Configuration

from IPython.display import Code

Code(filename=filename)

# Lid-driven cavity flow - moment-based boundary-condition validation (issue #921).
#
# A square cavity (L x L, 2-D D2Q9) with three stationary no-slip walls and a
# top lid translating at a constant velocity U. This is the canonical benchmark
# for validating moment-based boundary conditions in the LBM
# (Mohammed & Reis 2017): the lid is the moment-based moving wall
# `RegularizedVelocityWall` (imposes u = (U, 0), takes the density
# zero-normal-gradient from the interior, derives the rate-of-strain from the
# one-sided wall-normal finite difference), and the three static walls are the
# no-slip `RegularizedHWBB` (the u_wall = 0 case of the same closure). All walls
# use the default extrapolated (zero-normal-gradient) density - no `rho_wall` is
# set, which is the physically correct closure for a recirculating cavity where
# the wall pressure varies. No population bounce-back is used anywhere.
#
# Reference solution
# ------------------
# There is no closed form; the accepted benchmark is the multigrid Navier-Stokes
# solution of Ghia, Ghia & Shin (1982), tabulated as:
#   - u along the vertical centreline x = L/2, as a function of y, and
#   - v along the horizontal centreline y = L/2, as a function of x,
# for Re = 100, 400, 1000 (digitised tables in reference/). The simulation must
# reproduce these centreline profiles, including the sign reversals of the
# primary vortex and the near-wall extrema.
#
# Reynolds-number sweep
# ---------------------
# Re = U L / nu is swept over {100, 400, 1000} via `!unroll`, holding the grid
# (L = 256) and the lid speed (U = 0.05, so Ma = U*sqrt(3) ~= 0.087 < 0.1) fixed
# and varying the viscosity through tau = 3 nu + 1/2:
#   Re = 100  -> nu = 0.128   -> tau = 0.8840
#   Re = 400  -> nu = 0.032   -> tau = 0.5960
#   Re = 1000 -> nu = 0.0128  -> tau = 0.5384
# The 256-grid resolves the near-wall velocity extrema that the Ghia tables
# probe densely. Higher Re converges to steady state more slowly, so n_steps
# grows with Re. sim_id 0 = Re 100, 1 = Re 400, 2 = Re 1000 (the notebook
# indexes by this order).

simulations:
  - name: lidDrivenCavity
    save_path: !unroll
      - ./validation/analytical/05_lid_driven_cavity/results/re_100
      - ./validation/analytical/05_lid_driven_cavity/results/re_400
      - ./validation/analytical/05_lid_driven_cavity/results/re_1000

    n_steps: !unroll [300000, 500000, 800000]

    report:
      frequency: 40000

    domain:
      domain_size:
        x: 256
        y: 256
      block_size: 8

    data:
      # Default guard (rho non-finite -> stop, checked every 10 steps).
      monitors:
        fields:
          # Watch the flow settle to steady state (the velocity extrema stop
          # moving) and catch any blow-up early.
          macrs_stats:
            macrs: [rho, u]
            stats: [min, max, mean]
            interval: {frequency: 1000}
      exports:
        # ~20 volume snapshots per case so the notebook can confirm the
        # centreline profiles have reached steady state and read the last one.
        default:
          macrs: [rho, u]
          interval: !unroll
            - {frequency: 15000, lvl: 0}
            - {frequency: 25000, lvl: 0}
            - {frequency: 40000, lvl: 0}
          target:
            volumes:
              default: {}
          outputs:
            instantaneous: true

    models:
      precision:
        default: single
        calculations: double

      LBM:
        tau: !unroll [0.8840, 0.5960, 0.5384]
        vel_set: D2Q9
        coll_oper: RRBGK

      engine:
        name: CUDA

      BC:
        periodic_dims: [false, false]
        BC_map:
          # Moving lid: moment-based velocity wall translating at +x.
          - pos: N
            BC: RegularizedVelocityWall
            wall_normal: N
            order: 1
            params:
              ux: 0.05
              uy: 0
          # Three stationary no-slip walls (u_wall = 0 case of the same closure).
          - pos: S
            BC: RegularizedHWBB
            wall_normal: S
            order: 1
          - pos: W
            BC: RegularizedHWBB
            wall_normal: W
            order: 0
          - pos: E
            BC: RegularizedHWBB
            wall_normal: E
            order: 0