LES Intercomparison Study for Neutral Boundary Layers: Horizontal Cross-Sectional Views

Last updated: May 2026

Case Description: This notebook shows horizontal (x-y plane) cross-sections of the velocity components (\(u\), \(v\), \(w\)) from the modified Andren et al. (1994) neutral boundary layer case. Cross-sections are extracted at four height levels near 50 m, 100 m, 250 m, and 500 m. The 3D fields are sampled at the end of the simulation (~80 h). The reference run is 128x128x128_LASDD_SM_SP.

Load the necessary packages

import numpy as np
import matplotlib.pyplot as plt
from pathlib import Path

Input & Output Directories

# Base directory (jaxalfa/)
def find_repo_root(start=None):
    path = Path(start or ('__file__' in globals() and __file__) or Path.cwd()).resolve()
    for candidate in (path, *path.parents):
        if (candidate / 'examples').is_dir() and (candidate / 'docs').is_dir():
            return candidate
    raise FileNotFoundError('Could not locate jaxalfa repository root')

BaseDir = find_repo_root()

# Reference run: 128x128x128 LASDD-SM single precision
RunDir    = BaseDir / 'examples/NBL_A94/runs/128x128x128_LASDD_SM_DP'
OutputDir = RunDir / 'output'

cfg = {}
exec((RunDir / 'Config.py').read_text(), cfg)

Load 3D fields from T = 80 h

T_snapshot = float(cfg['SimTime'])
dt = float(cfg['dt'])
iter_3D = int(T_snapshot / dt)
field_path = OutputDir / f'ALFA_3DFields_Iteration_{iter_3D}.npz'

if field_path.exists():
    File3D = np.load(field_path)
    u3D  = File3D['u']
    v3D  = File3D['v']
    w3D  = File3D['w']
else:
    nx = int(cfg['nx']); ny = int(cfg['ny']); nz = int(cfg['nz'])
    print(f'Missing {field_path}; plotting NaN placeholders for this run.')
    u3D = np.full((nx, ny, nz), np.nan)
    v3D = np.full((nx, ny, nz), np.nan)
    w3D = np.full((nx, ny, nz), np.nan)

Input Information from the Config File

l_x = float(cfg['l_x'])
l_y = float(cfg['l_y'])
l_z = float(cfg['l_z'])
SimTime = float(cfg['SimTime'])

nx = int(cfg['nx'])
ny = int(cfg['ny'])
nz = int(cfg['nz'])
dx = l_x / nx
dy = l_y / ny
dz = l_z / (nz - 1)
x_axis = dx * np.arange(nx)
y_axis = dy * np.arange(ny)

Derived Variables

# Half levels for u, v variables
z_u = np.array([(k + 0.5) * l_z / (nz - 1) for k in range(nz)])

# Full levels for w
z_w = np.array([k * l_z / (nz - 1) for k in range(nz)])

plt.rcParams.update({
    "text.usetex": True,
    "font.size": 14,
    "axes.labelsize": 16,
    "xtick.labelsize": 12,
    "ytick.labelsize": 12
})

Plot horizontal cross-section of longitudinal velocity fields

fig, axes = plt.subplots(2, 2, figsize=(10, 10), constrained_layout=True)
axes = axes.flatten()

# Target heights (m) — nearest grid level selected automatically
target_heights = [50, 100, 250, 500]
k_levels = [int(np.argmin(np.abs(z_u - h))) for h in target_heights]

for i, k in enumerate(k_levels):
    im = axes[i].pcolor(x_axis, y_axis, u3D[:,:,k].T, cmap='inferno')
    axes[i].set_title(f'Longitudinal Velocity at Z = {z_u[k]:.1f} m', fontsize=12)
    axes[i].set_xlabel('X (m)')
    axes[i].set_ylabel('Y (m)')
    axes[i].set_aspect('auto')
    fig.colorbar(im, ax=axes[i], label='u (m/s)')

plt.suptitle('Simulation: 128x128x128 (LASDD-SM, SP)', fontsize=16)
plt.show()

Plot horizontal cross-section of lateral velocity fields

fig, axes = plt.subplots(2, 2, figsize=(10, 10), constrained_layout=True)
axes = axes.flatten()

# Target heights (m) — nearest grid level selected automatically
target_heights = [50, 100, 250, 500]
k_levels = [int(np.argmin(np.abs(z_u - h))) for h in target_heights]

for i, k in enumerate(k_levels):
    im = axes[i].pcolor(x_axis, y_axis, v3D[:,:,k].T, cmap='inferno')
    axes[i].set_title(f'Lateral Velocity at Z = {z_u[k]:.1f} m', fontsize=12)
    axes[i].set_xlabel('X (m)')
    axes[i].set_ylabel('Y (m)')
    axes[i].set_aspect('auto')
    fig.colorbar(im, ax=axes[i], label='v (m/s)')

plt.suptitle('Simulation: 128x128x128 (LASDD-SM, SP)', fontsize=16)
plt.show()

Plot horizontal cross-section of vertical velocity fields

fig, axes = plt.subplots(2, 2, figsize=(10, 10), constrained_layout=True)
axes = axes.flatten()

# Target heights (m) — nearest grid level selected automatically
target_heights = [50, 100, 250, 500]
k_levels = [int(np.argmin(np.abs(z_w - h))) for h in target_heights]

for i, k in enumerate(k_levels):
    im = axes[i].pcolor(x_axis, y_axis, w3D[:,:,k].T, cmap='inferno')
    axes[i].set_title(f'Vertical Velocity at Z = {z_w[k]:.1f} m', fontsize=12)
    axes[i].set_xlabel('X (m)')
    axes[i].set_ylabel('Y (m)')
    axes[i].set_aspect('auto')
    fig.colorbar(im, ax=axes[i], label='w (m/s)')

plt.suptitle('Simulation: 128x128x128 (LASDD-SM, SP)', fontsize=16)
plt.show()