LES Intercomparison Study for Neutral Boundary Layers: Sensitivity with respect to Spatial Resolution

Last updated: May 2026

For case setup and physical parameters, see the Description notebook.

Resolutions compared: \(64^3\), \(128^3\), \(256^3\), \(384^3\); SGS: LASDD-SM (DP).

Setup

The next cells load Python packages, locate the simulation outputs, and define the grid and averaging window used throughout the notebook.

import os
import re
import glob
import numpy as np
import matplotlib.pyplot as plt
from pathlib import Path

Output directories

from pathlib import Path

# 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()

def read_config(run_dir):
    cfg = {}
    exec((run_dir / 'Config.py').read_text(), cfg)
    return cfg


optSGS = 1  # LASDD-SM: 1, LASDD-WL: 2, LAD-SM: 3, LAD-WL: 4
optPrecision = 'DP'  # 'DP' or 'SP'

_sgs = {1: 'LASDD_SM', 2: 'LASDD_WL', 3: 'LAD_SM', 4: 'LAD_WL'}
_run = lambda res: f"{res}_{_sgs[optSGS]}_{optPrecision}"

OutputDir1 = BaseDir / f'examples/NBL_A94/runs/{_run("64x64x64")}/output'
OutputDir2 = BaseDir / f'examples/NBL_A94/runs/{_run("128x128x128")}/output'
OutputDir3 = BaseDir / f'examples/NBL_A94/runs/{_run("256x256x256")}/output'
OutputDir4 = BaseDir / f'examples/NBL_A94/runs/{_run("384x384x384")}/output'

Case configuration

cfg_1 = read_config(OutputDir1.parent)
cfg_2 = read_config(OutputDir2.parent)
cfg_3 = read_config(OutputDir3.parent)
cfg_4 = read_config(OutputDir4.parent)

nz_1 = int(cfg_1['nz'])
nz_2 = int(cfg_2['nz'])
nz_3 = int(cfg_3['nz'])
nz_4 = int(cfg_4['nz'])

l_z = float(cfg_1['l_z'])
z_damping = float(cfg_1.get('z_damping', np.nan))
RelaxTime = float(cfg_1.get('RelaxTime', np.nan))
OutputInterval_sec = float(cfg_1.get('OutputInterval_sec', 60.0))

# Averaging window — NBL quasi-steady state (last 10 h of an 83 h run)
T_start = 73 * 3600   # s
T_end   = 83 * 3600   # s

Derived grid and averaging indices

# Half levels — u, v, TH
z_1 = np.array([(k + 0.5) * l_z / (nz_1 - 1) for k in range(nz_1)])
z_2 = np.array([(k + 0.5) * l_z / (nz_2 - 1) for k in range(nz_2)])
z_3 = np.array([(k + 0.5) * l_z / (nz_3 - 1) for k in range(nz_3)])
z_4 = np.array([(k + 0.5) * l_z / (nz_4 - 1) for k in range(nz_4)])

# Full levels — w, uw, vw, wTH, qz
z_w_1 = np.array([k * l_z / (nz_1 - 1) for k in range(nz_1)])
z_w_2 = np.array([k * l_z / (nz_2 - 1) for k in range(nz_2)])
z_w_3 = np.array([k * l_z / (nz_3 - 1) for k in range(nz_3)])
z_w_4 = np.array([k * l_z / (nz_4 - 1) for k in range(nz_4)])

# File indices for the averaging window (same for all runs: same OutputInterval_sec)
T_start_index = int(T_start / OutputInterval_sec) - 1
T_end_index   = int(T_end   / OutputInterval_sec) - 1

print(f'Averaging window: file indices {T_start_index} – {T_end_index}')

Averaging window: file indices 4379 – 4979

Statistics loader

def LoadStatsAverage(stat_files, T_start_index, T_end_index, nz_expected):
    if len(stat_files) == 0:
        print(f'No statistics files available; plotting NaN placeholders for nz={nz_expected}.')
        nan = np.full(nz_expected, np.nan)
        return tuple(nan.copy() for _ in range(15))

    U   = []; V   = []; TH  = []
    u2  = []; v2  = []; w2  = []; TH2 = []
    uv  = []; uw  = []; vw  = []
    txy = []; txz = []; tyz = []
    wTH = []; qz  = []

    for f in stat_files:
        with np.load(f) as d:
            U.append(d['U']);   V.append(d['V']);   TH.append(d['TH'])
            u2.append(d['u2']); v2.append(d['v2']); w2.append(d['w2'])
            TH2.append(d['TH2'])
            uv.append(d['uv']); uw.append(d['uw']); vw.append(d['vw'])
            txy.append(d['txy']); txz.append(d['txz']); tyz.append(d['tyz'])
            wTH.append(d['wTH']); qz.append(d['qz'])

    U  = np.array(U);  V  = np.array(V);  TH  = np.array(TH)
    u2 = np.array(u2); v2 = np.array(v2); w2  = np.array(w2); TH2 = np.array(TH2)
    uv = np.array(uv); uw = np.array(uw); vw  = np.array(vw)
    txy = np.array(txy); txz = np.array(txz); tyz = np.array(tyz)
    wTH = np.array(wTH); qz = np.array(qz)

    sl = slice(T_start_index, min(T_end_index + 1, len(stat_files)))
    if sl.start >= len(stat_files):
        print(f'Averaging window starts after available files; plotting NaN placeholders for nz={nz_expected}.')
        nan = np.full(nz_expected, np.nan)
        return tuple(nan.copy() for _ in range(15))

    return (
        np.mean(U[sl],   axis=0), np.mean(V[sl],   axis=0), np.mean(TH[sl],  axis=0),
        np.mean(u2[sl],  axis=0), np.mean(v2[sl],  axis=0), np.mean(w2[sl],  axis=0),
        np.mean(TH2[sl], axis=0),
        np.mean(uv[sl],  axis=0), np.mean(uw[sl],  axis=0), np.mean(vw[sl],  axis=0),
        np.mean(txy[sl], axis=0), np.mean(txz[sl], axis=0), np.mean(tyz[sl], axis=0),
        np.mean(wTH[sl], axis=0), np.mean(qz[sl],  axis=0)
    )

Available statistics files

def get_stat_files(output_dir):
    files = sorted(
        glob.glob(str(output_dir / 'ALFA_Statistics_Iteration_*.npz')),
        key=lambda x: int(re.search(r'Iteration_(\d+)', x).group(1))
    )
    return files

StatFiles1 = get_stat_files(OutputDir1)
StatFiles2 = get_stat_files(OutputDir2)
StatFiles3 = get_stat_files(OutputDir3)
StatFiles4 = get_stat_files(OutputDir4)

print(f'64^3  : {len(StatFiles1)} files')
print(f'128^3 : {len(StatFiles2)} files')
print(f'256^3 : {len(StatFiles3)} files')
print(f'384^3 : {len(StatFiles4)} files')

64^3  : 5000 files
128^3 : 5000 files
256^3 : 5000 files
384^3 : 0 files

Temporally averaged profiles

(U_avg_1, V_avg_1, TH_avg_1,
 u2_avg_1, v2_avg_1, w2_avg_1, TH2_avg_1,
 uv_avg_1, uw_avg_1, vw_avg_1,
 txy_avg_1, txz_avg_1, tyz_avg_1,
 wTH_avg_1, qz_avg_1) = LoadStatsAverage(StatFiles1, T_start_index, T_end_index, nz_1)

(U_avg_2, V_avg_2, TH_avg_2,
 u2_avg_2, v2_avg_2, w2_avg_2, TH2_avg_2,
 uv_avg_2, uw_avg_2, vw_avg_2,
 txy_avg_2, txz_avg_2, tyz_avg_2,
 wTH_avg_2, qz_avg_2) = LoadStatsAverage(StatFiles2, T_start_index, T_end_index, nz_2)

(U_avg_3, V_avg_3, TH_avg_3,
 u2_avg_3, v2_avg_3, w2_avg_3, TH2_avg_3,
 uv_avg_3, uw_avg_3, vw_avg_3,
 txy_avg_3, txz_avg_3, tyz_avg_3,
 wTH_avg_3, qz_avg_3) = LoadStatsAverage(StatFiles3, T_start_index, T_end_index, nz_3)

(U_avg_4, V_avg_4, TH_avg_4,
 u2_avg_4, v2_avg_4, w2_avg_4, TH2_avg_4,
 uv_avg_4, uw_avg_4, vw_avg_4,
 txy_avg_4, txz_avg_4, tyz_avg_4,
 wTH_avg_4, qz_avg_4) = LoadStatsAverage(StatFiles4, T_start_index, T_end_index, nz_4)

S_avg_1   = np.sqrt(U_avg_1**2 + V_avg_1**2)
uw_tot_1  = uw_avg_1  + txz_avg_1
vw_tot_1  = vw_avg_1  + tyz_avg_1
wTH_tot_1 = wTH_avg_1 + qz_avg_1

S_avg_2   = np.sqrt(U_avg_2**2 + V_avg_2**2)
uw_tot_2  = uw_avg_2  + txz_avg_2
vw_tot_2  = vw_avg_2  + tyz_avg_2
wTH_tot_2 = wTH_avg_2 + qz_avg_2

S_avg_3   = np.sqrt(U_avg_3**2 + V_avg_3**2)
uw_tot_3  = uw_avg_3  + txz_avg_3
vw_tot_3  = vw_avg_3  + tyz_avg_3
wTH_tot_3 = wTH_avg_3 + qz_avg_3

S_avg_4   = np.sqrt(U_avg_4**2 + V_avg_4**2)
uw_tot_4  = uw_avg_4  + txz_avg_4
vw_tot_4  = vw_avg_4  + tyz_avg_4
wTH_tot_4 = wTH_avg_4 + qz_avg_4

print(f'Averaging over {T_end_index - T_start_index + 1} files '
      f'({T_start/3600:.1f}–{T_end/3600:.1f} h)')

No statistics files available; plotting NaN placeholders for nz=384.
Averaging over 601 files (73.0–83.0 h)

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

def plot_profile(x, z, xlabel, ylabel=r"$z$ (m)", linestyle='-k', label=None, ax=None):
    if ax is None:
        fig, ax = plt.subplots(figsize=(5, 6), constrained_layout=True)

    ax.plot(x, z, linestyle, linewidth=2, label=label)
    ax.set_xlabel(xlabel)
    ax.set_ylabel(ylabel)
    ax.grid(False)

Mean Wind and Hodograph

The first comparison shows the mean streamwise and cross-stream wind components, the wind-speed magnitude, and the hodograph over the last 10 h averaging window.

fig, axs = plt.subplots(2, 2, figsize=(10, 10), constrained_layout=True)
axs = axs.ravel()

run_styles = {
    '64x64x64':    {'color': 'red',   'linestyle': '-'},
    '128x128x128': {'color': 'blue',  'linestyle': '-'},
    '256x256x256': {'color': 'green', 'linestyle': '-'},
    '384x384x384': {'color': 'black', 'linestyle': '-'},
}

def plot_run_profile(ax, x, z, xlabel, run_label):
    style = run_styles[run_label]
    ax.plot(x, z, color=style['color'], linestyle=style['linestyle'], linewidth=2, label=run_label)
    ax.set_xlabel(xlabel)
    ax.set_ylabel(r"$z$ (m)")

for lbl, U, V, S, z in [('64x64x64',    U_avg_1, V_avg_1, S_avg_1, z_1),
                          ('128x128x128',    U_avg_2, V_avg_2, S_avg_2, z_2),
                          ('256x256x256', U_avg_3, V_avg_3, S_avg_3, z_3),
                          ('384x384x384', U_avg_4, V_avg_4, S_avg_4, z_4)]:
    plot_run_profile(axs[0], U, z, r"$U$ (m/s)", lbl)
    plot_run_profile(axs[1], V, z, r"$V$ (m/s)", lbl)
    plot_run_profile(axs[2], S, z, r"Wind Speed (m/s)", lbl)

for lbl, U, V, style in [('64x64x64',    U_avg_1, V_avg_1, run_styles['64x64x64']),
                           ('128x128x128',    U_avg_2, V_avg_2, run_styles['128x128x128']),
                           ('256x256x256', U_avg_3, V_avg_3, run_styles['256x256x256']),
                           ('384x384x384', U_avg_4, V_avg_4, run_styles['384x384x384'])]:
    axs[3].plot(U, V, color=style['color'], linestyle='-', marker='o',
                linewidth=2, markersize=3, label=lbl)
axs[3].set_xlabel(r"$U$ (m/s)")
axs[3].set_ylabel(r"$V$ (m/s)")
axs[3].set_title('Hodograph')
axs[3].set_aspect('equal')

for ax in axs:
    ax.grid()
    ax.legend(frameon=False)
fig.suptitle("Mean Wind Profiles and Hodograph (last 10 h average)", fontsize=18)
plt.show()

Resolved Velocity Variances

The resolved variance profiles indicate how the resolved turbulent kinetic energy is distributed among the three velocity components.

fig, axs = plt.subplots(1, 3, figsize=(15, 5), constrained_layout=True)

plot_profile(u2_avg_1, z_1,   xlabel=r"Resolved $\sigma_u^2$ (m$^2$/s$^2$)", linestyle='-r', ax=axs[0], label='64x64x64')
plot_profile(v2_avg_1, z_1,   xlabel=r"Resolved $\sigma_v^2$ (m$^2$/s$^2$)", linestyle='-r', ax=axs[1], label='64x64x64')
plot_profile(w2_avg_1, z_w_1, xlabel=r"Resolved $\sigma_w^2$ (m$^2$/s$^2$)", linestyle='-r', ax=axs[2], label='64x64x64')
plot_profile(u2_avg_2, z_2,   xlabel=r"Resolved $\sigma_u^2$ (m$^2$/s$^2$)", linestyle='-b', ax=axs[0], label='128x128x128')
plot_profile(v2_avg_2, z_2,   xlabel=r"Resolved $\sigma_v^2$ (m$^2$/s$^2$)", linestyle='-b', ax=axs[1], label='128x128x128')
plot_profile(w2_avg_2, z_w_2, xlabel=r"Resolved $\sigma_w^2$ (m$^2$/s$^2$)", linestyle='-b', ax=axs[2], label='128x128x128')
plot_profile(u2_avg_3, z_3,   xlabel=r"Resolved $\sigma_u^2$ (m$^2$/s$^2$)", linestyle='-g', ax=axs[0], label='256x256x256')
plot_profile(v2_avg_3, z_3,   xlabel=r"Resolved $\sigma_v^2$ (m$^2$/s$^2$)", linestyle='-g', ax=axs[1], label='256x256x256')
plot_profile(w2_avg_3, z_w_3, xlabel=r"Resolved $\sigma_w^2$ (m$^2$/s$^2$)", linestyle='-g', ax=axs[2], label='256x256x256')
plot_profile(u2_avg_4, z_4,   xlabel=r"Resolved $\sigma_u^2$ (m$^2$/s$^2$)", linestyle='-k', ax=axs[0], label='384x384x384')
plot_profile(v2_avg_4, z_4,   xlabel=r"Resolved $\sigma_v^2$ (m$^2$/s$^2$)", linestyle='-k', ax=axs[1], label='384x384x384')
plot_profile(w2_avg_4, z_w_4, xlabel=r"Resolved $\sigma_w^2$ (m$^2$/s$^2$)", linestyle='-k', ax=axs[2], label='384x384x384')
axs[0].grid(); axs[0].legend(frameon=False)
axs[1].grid(); axs[1].legend(frameon=False)
axs[2].grid(); axs[2].legend(frameon=False)

fig.suptitle("Profiles of Velocity Variances (last 10 h average)", fontsize=18)
plt.show()

Resolved Momentum Fluxes

These profiles show the resolved turbulent momentum fluxes before adding the SGS contribution.

fig, axs = plt.subplots(1, 3, figsize=(15, 5), constrained_layout=True)

plot_profile(uv_avg_1, z_1,   xlabel=r"$\langle u'v' \rangle$ (m$^2$/s$^2$)", linestyle='-r', ax=axs[0], label='64x64x64')
plot_profile(uw_avg_1, z_w_1, xlabel=r"$\langle u'w' \rangle$ (m$^2$/s$^2$)", linestyle='-r', ax=axs[1], label='64x64x64')
plot_profile(vw_avg_1, z_w_1, xlabel=r"$\langle v'w' \rangle$ (m$^2$/s$^2$)", linestyle='-r', ax=axs[2], label='64x64x64')
plot_profile(uv_avg_2, z_2,   xlabel=r"$\langle u'v' \rangle$ (m$^2$/s$^2$)", linestyle='-b', ax=axs[0], label='128x128x128')
plot_profile(uw_avg_2, z_w_2, xlabel=r"$\langle u'w' \rangle$ (m$^2$/s$^2$)", linestyle='-b', ax=axs[1], label='128x128x128')
plot_profile(vw_avg_2, z_w_2, xlabel=r"$\langle v'w' \rangle$ (m$^2$/s$^2$)", linestyle='-b', ax=axs[2], label='128x128x128')
plot_profile(uv_avg_3, z_3,   xlabel=r"$\langle u'v' \rangle$ (m$^2$/s$^2$)", linestyle='-g', ax=axs[0], label='256x256x256')
plot_profile(uw_avg_3, z_w_3, xlabel=r"$\langle u'w' \rangle$ (m$^2$/s$^2$)", linestyle='-g', ax=axs[1], label='256x256x256')
plot_profile(vw_avg_3, z_w_3, xlabel=r"$\langle v'w' \rangle$ (m$^2$/s$^2$)", linestyle='-g', ax=axs[2], label='256x256x256')
plot_profile(uv_avg_4, z_4,   xlabel=r"$\langle u'v' \rangle$ (m$^2$/s$^2$)", linestyle='-k', ax=axs[0], label='384x384x384')
plot_profile(uw_avg_4, z_w_4, xlabel=r"$\langle u'w' \rangle$ (m$^2$/s$^2$)", linestyle='-k', ax=axs[1], label='384x384x384')
plot_profile(vw_avg_4, z_w_4, xlabel=r"$\langle v'w' \rangle$ (m$^2$/s$^2$)", linestyle='-k', ax=axs[2], label='384x384x384')
axs[0].grid(); axs[0].legend(frameon=False)
axs[1].grid(); axs[1].legend(frameon=False)
axs[2].grid(); axs[2].legend(frameon=False)

fig.suptitle("Profiles of Resolved Momentum Fluxes (last 10 h average)", fontsize=18)
plt.show()

SGS Stresses

The SGS stress profiles show the modeled stress contribution retained by each grid resolution.

fig, axs = plt.subplots(1, 3, figsize=(15, 5), constrained_layout=True)

plot_profile(txy_avg_1, z_1,   xlabel=r"$\tau_{xy}$ (m$^2$/s$^2$)", linestyle='-r', ax=axs[0], label='64x64x64')
plot_profile(txz_avg_1, z_w_1, xlabel=r"$\tau_{xz}$ (m$^2$/s$^2$)", linestyle='-r', ax=axs[1], label='64x64x64')
plot_profile(tyz_avg_1, z_w_1, xlabel=r"$\tau_{yz}$ (m$^2$/s$^2$)", linestyle='-r', ax=axs[2], label='64x64x64')
plot_profile(txy_avg_2, z_2,   xlabel=r"$\tau_{xy}$ (m$^2$/s$^2$)", linestyle='-b', ax=axs[0], label='128x128x128')
plot_profile(txz_avg_2, z_w_2, xlabel=r"$\tau_{xz}$ (m$^2$/s$^2$)", linestyle='-b', ax=axs[1], label='128x128x128')
plot_profile(tyz_avg_2, z_w_2, xlabel=r"$\tau_{yz}$ (m$^2$/s$^2$)", linestyle='-b', ax=axs[2], label='128x128x128')
plot_profile(txy_avg_3, z_3,   xlabel=r"$\tau_{xy}$ (m$^2$/s$^2$)", linestyle='-g', ax=axs[0], label='256x256x256')
plot_profile(txz_avg_3, z_w_3, xlabel=r"$\tau_{xz}$ (m$^2$/s$^2$)", linestyle='-g', ax=axs[1], label='256x256x256')
plot_profile(tyz_avg_3, z_w_3, xlabel=r"$\tau_{yz}$ (m$^2$/s$^2$)", linestyle='-g', ax=axs[2], label='256x256x256')
plot_profile(txy_avg_4, z_4,   xlabel=r"$\tau_{xy}$ (m$^2$/s$^2$)", linestyle='-k', ax=axs[0], label='384x384x384')
plot_profile(txz_avg_4, z_w_4, xlabel=r"$\tau_{xz}$ (m$^2$/s$^2$)", linestyle='-k', ax=axs[1], label='384x384x384')
plot_profile(tyz_avg_4, z_w_4, xlabel=r"$\tau_{yz}$ (m$^2$/s$^2$)", linestyle='-k', ax=axs[2], label='384x384x384')
axs[0].grid(); axs[0].legend(frameon=False)
axs[1].grid(); axs[1].legend(frameon=False)
axs[2].grid(); axs[2].legend(frameon=False)

fig.suptitle("Profiles of SGS Stresses (last 10 h average)", fontsize=18)
plt.show()

Total Momentum Fluxes

The total vertical momentum flux combines resolved and SGS contributions.

fig, axs = plt.subplots(1, 2, figsize=(10, 5), constrained_layout=True)

plot_profile(uw_tot_1, z_w_1, xlabel=r"$\langle u'w' \rangle + \tau_{xz}$ (m$^2$/s$^2$)", linestyle='-r', label='64x64x64',    ax=axs[0])
plot_profile(uw_tot_2, z_w_2, xlabel=r"$\langle u'w' \rangle + \tau_{xz}$ (m$^2$/s$^2$)", linestyle='-b', label='128x128x128',    ax=axs[0])
plot_profile(uw_tot_3, z_w_3, xlabel=r"$\langle u'w' \rangle + \tau_{xz}$ (m$^2$/s$^2$)", linestyle='-g', label='256x256x256', ax=axs[0])
plot_profile(uw_tot_4, z_w_4, xlabel=r"$\langle u'w' \rangle + \tau_{xz}$ (m$^2$/s$^2$)", linestyle='-k', label='384x384x384', ax=axs[0])
axs[0].grid()
axs[0].legend(frameon=False)

plot_profile(vw_tot_1, z_w_1, xlabel=r"$\langle v'w' \rangle + \tau_{yz}$ (m$^2$/s$^2$)", linestyle='-r', label='64x64x64',    ax=axs[1])
plot_profile(vw_tot_2, z_w_2, xlabel=r"$\langle v'w' \rangle + \tau_{yz}$ (m$^2$/s$^2$)", linestyle='-b', label='128x128x128',    ax=axs[1])
plot_profile(vw_tot_3, z_w_3, xlabel=r"$\langle v'w' \rangle + \tau_{yz}$ (m$^2$/s$^2$)", linestyle='-g', label='256x256x256', ax=axs[1])
plot_profile(vw_tot_4, z_w_4, xlabel=r"$\langle v'w' \rangle + \tau_{yz}$ (m$^2$/s$^2$)", linestyle='-k', label='384x384x384', ax=axs[1])
axs[1].grid()
axs[1].legend(frameon=False)

fig.suptitle("Total Momentum Fluxes: Resolved + SGS (last 10 h average)", fontsize=18)
plt.show()