LES Intercomparison Study for Neutral Boundary Layers: Sensitivity with respect to Numerical Precision

Last updated: May 2026

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

Precisions compared: double (DP), single (SP); SGS: LASDD-SM, grid: \(128^3\).

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


optRes = 2  # 64x64x64: 1, 128x128x128: 2, 256x256x256: 3, 384x384x384: 4

_res = {1: '64x64x64', 2: '128x128x128', 3: '256x256x256', 4: '384x384x384'}
_res_label = {1: r'$64^3$', 2: r'$128^3$', 3: r'$256^3$', 4: r'$384^3$'}

_run = f"{_res[optRes]}_LASDD_SM"

OutputDir_DP = BaseDir / f'examples/NBL_A94/runs/{_run}_DP/output'
OutputDir_SP = BaseDir / f'examples/NBL_A94/runs/{_run}_SP/output'

Case configuration

_cfg_ref = read_config(OutputDir_DP.parent if (OutputDir_DP.parent / 'Config.py').exists() else OutputDir_SP.parent)
nz = int(_cfg_ref['nz'])
l_z = float(_cfg_ref['l_z'])
OutputInterval_sec = float(_cfg_ref.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) for k in range(nz)])

# Full levels — w, uw, vw, wTH, qz
z_w_1 = np.array([k * l_z / (nz - 1) for k in range(nz)])

# File indices for the averaging window
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

StatFiles_DP = get_stat_files(OutputDir_DP)
StatFiles_SP = get_stat_files(OutputDir_SP)

Temporally averaged profiles

(U_avg_DP, V_avg_DP, TH_avg_DP,
 u2_avg_DP, v2_avg_DP, w2_avg_DP, TH2_avg_DP,
 uv_avg_DP, uw_avg_DP, vw_avg_DP,
 txy_avg_DP, txz_avg_DP, tyz_avg_DP,
 wTH_avg_DP, qz_avg_DP) = LoadStatsAverage(StatFiles_DP, T_start_index, T_end_index, nz)

(U_avg_SP, V_avg_SP, TH_avg_SP,
 u2_avg_SP, v2_avg_SP, w2_avg_SP, TH2_avg_SP,
 uv_avg_SP, uw_avg_SP, vw_avg_SP,
 txy_avg_SP, txz_avg_SP, tyz_avg_SP,
 wTH_avg_SP, qz_avg_SP) = LoadStatsAverage(StatFiles_SP, T_start_index, T_end_index, nz)

M_avg_DP   = np.sqrt(U_avg_DP**2 + V_avg_DP**2)
uw_tot_DP  = uw_avg_DP  + txz_avg_DP
vw_tot_DP  = vw_avg_DP  + tyz_avg_DP

M_avg_SP   = np.sqrt(U_avg_SP**2 + V_avg_SP**2)
uw_tot_SP  = uw_avg_SP  + txz_avg_SP
vw_tot_SP  = vw_avg_SP  + tyz_avg_SP

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

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=(11, 10), constrained_layout=True)
axs = axs.ravel()

run_styles = {
    'DP': {'color': 'red',  'linestyle': '-'},
    'SP': {'color': 'blue', '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)")

plot_run_profile(axs[0], U_avg_DP, z_1, r"$U$ (m/s)", 'DP')
plot_run_profile(axs[0], U_avg_SP, z_1, r"$U$ (m/s)", 'SP')

plot_run_profile(axs[1], V_avg_DP, z_1, r"$V$ (m/s)", 'DP')
plot_run_profile(axs[1], V_avg_SP, z_1, r"$V$ (m/s)", 'SP')

plot_run_profile(axs[2], M_avg_DP, z_1, r"Wind Speed (m/s)", 'DP')
plot_run_profile(axs[2], M_avg_SP, z_1, r"Wind Speed (m/s)", 'SP')

axs[3].plot(U_avg_DP, V_avg_DP, color='red',  linestyle='-', marker='o', linewidth=2, markersize=3, label='DP')
axs[3].plot(U_avg_SP, V_avg_SP, color='blue', linestyle='-', marker='o', linewidth=2, markersize=3, label='SP')
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(
    f"Mean Wind Profiles and Hodograph (last 10 h average): "
    f"LASDD-SM, {_res_label[optRes]}",
    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_DP, z_1,   xlabel=r"Resolved $\sigma_u^2$ (m$^2$/s$^2$)", linestyle='-r', ax=axs[0], label='DP')
plot_profile(u2_avg_SP, z_1,   xlabel=r"Resolved $\sigma_u^2$ (m$^2$/s$^2$)", linestyle='-b', ax=axs[0], label='SP')

plot_profile(v2_avg_DP, z_1,   xlabel=r"Resolved $\sigma_v^2$ (m$^2$/s$^2$)", linestyle='-r', ax=axs[1], label='DP')
plot_profile(v2_avg_SP, z_1,   xlabel=r"Resolved $\sigma_v^2$ (m$^2$/s$^2$)", linestyle='-b', ax=axs[1], label='SP')

plot_profile(w2_avg_DP, z_w_1, xlabel=r"Resolved $\sigma_w^2$ (m$^2$/s$^2$)", linestyle='-r', ax=axs[2], label='DP')
plot_profile(w2_avg_SP, z_w_1, xlabel=r"Resolved $\sigma_w^2$ (m$^2$/s$^2$)", linestyle='-b', ax=axs[2], label='SP')

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(
    f"Resolved Velocity Variances (last 10 h average): "
    f"LASDD-SM, {_res_label[optRes]}",
    fontsize=18
)
plt.show()

Total Momentum Fluxes

The total vertical momentum fluxes combine resolved and SGS contributions. Double- and single-precision runs are compared at the selected resolution.

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

plot_profile(uw_tot_DP, z_w_1, xlabel=r"Total $\langle u'w' \rangle$ (m$^2$/s$^2$)", linestyle='-r', label='DP', ax=axs[0])
plot_profile(uw_tot_SP, z_w_1, xlabel=r"Total $\langle u'w' \rangle$ (m$^2$/s$^2$)", linestyle='-b', label='SP', ax=axs[0])
axs[0].grid()
axs[0].legend(frameon=False)

plot_profile(vw_tot_DP, z_w_1, xlabel=r"Total $\langle v'w' \rangle$ (m$^2$/s$^2$)", linestyle='-r', label='DP', ax=axs[1])
plot_profile(vw_tot_SP, z_w_1, xlabel=r"Total $\langle v'w' \rangle$ (m$^2$/s$^2$)", linestyle='-b', label='SP', ax=axs[1])
axs[1].grid()
axs[1].legend(frameon=False)

fig.suptitle(
    f"Total (Resolved + SGS) Momentum Fluxes (last 10 h average): "
    f"LASDD-SM, {_res_label[optRes]}",
    fontsize=18
)
plt.show()