SEM Inlet Reynolds-Stress Tensor and Spectrum

Validates that the Synthetic Eddy Method (SEM) inlet reproduces the prescribed Reynolds-stress tensor R_ij (including the shear components Rxz and Ryz) and a target turbulence spectrum, sampled on a dense inlet-normal plane a few nodes downstream of the inlet.

The flow is statistically homogeneous in the spanwise y and vertical z (both periodic) and the prescribed profile is constant with height, so the sampled second moments are pooled over the whole plane and over time. The pass criteria are:

  • Each sampled covariance <u'_i u'_j> matches the prescribed R_ij within ~10%, with correct sign and magnitude of the shear terms.

  • The 1D streamwise energy spectrum follows the von Karman form in the inertial range.

````{note} Results pending GPU execution. This notebook is wired against the prescribed target and the export schema. The target tensor and all plotting code run as a scaffold now; the measured curves and the verdict populate automatically once the case has been run with

uv run nassu run validation/inlet/01_sem_reynolds_stress/01_sem_reynolds_stress.nassu.yaml

(the RESULTS_AVAILABLE flag below switches on when the plane output exists). ````

Load the configuration

Parse the same YAML the runner ran. All parameters (mean velocity, prescribed R_ij, sampling) are read off the parsed config and the reference profile CSV, never hard-coded.

[ ]:

import pathlib

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

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

common.use_style()

CASE_DIR = pathlib.Path("validation/inlet/01_sem_reynolds_stress")
CONFIG = CASE_DIR / "01_sem_reynolds_stress.nassu.yaml"

sim_cfgs_list = ConfigScheme.sim_cfgs_from_file_dct(str(CONFIG))
sim_cfg = list(sim_cfgs_list.values())[0]
print(sim_cfg.name, "n_steps =", sim_cfg.n_steps, "domain =", sim_cfg.domain.domain_size)

DEV_TIME = 5000  # development steps discarded from the series

# Whether the case has been run: the plane XDMF exists. When False the notebook
# runs as a scaffold (target + plotting code present, measured data absent).
_plane = sim_cfg.output.series["inlet_plane"].planes["measurement"]
RESULTS_AVAILABLE = _plane.xdmf_filename.exists()
print("RESULTS_AVAILABLE =", RESULTS_AVAILABLE)

Prescribed target tensor

Read the prescribed R_ij from the SEM profile CSV (constant with height) and verify it is symmetric positive-definite, so its Cholesky factor A (u'_a = A_ab u~_b) is real. This is the target the sampled covariance must reproduce.

[ ]:

sem = sim_cfg.models.BC.inlet_turbulence
profile_csv = sem.profile.csv_profile_data
df_profile = pd.read_csv(profile_csv)

# Constant with height: take the prescribed values at any interior row.
row = df_profile.iloc[len(df_profile) // 2]
U_ref = float(row["ux"])
R_target = np.array(
    [
        [row["Rxx"], row["Rxy"], row["Rxz"]],
        [row["Rxy"], row["Ryy"], row["Ryz"]],
        [row["Rxz"], row["Ryz"], row["Rzz"]],
    ],
    dtype=float,
)

eigs = np.linalg.eigvalsh(R_target)
A_chol = np.linalg.cholesky(R_target)  # lower-triangular, A A^T = R

# Component bookkeeping shared by the table and the plot below.
COMP_NAMES = ["Rxx", "Ryy", "Rzz", "Rxy", "Rxz", "Ryz"]
COMP_IDX = [(0, 0), (1, 1), (2, 2), (0, 1), (0, 2), (1, 2)]

print("U_ref =", U_ref)
print("R_target =\n", R_target)
print("eigenvalues =", eigs, "-> SPD:", bool(np.all(eigs > 0)))
print("Iu, Iv, Iw =", np.sqrt(np.diag(R_target)) / U_ref)
print("rho_uw =", R_target[0, 2] / np.sqrt(R_target[0, 0] * R_target[2, 2]))
print("Cholesky A =\n", A_chol)

Load the sampled inlet plane

Read the inlet_plane.measurement plane series. read_full_data(macr) returns one column per sample point (indexed by global point idx) plus a time_step column; the companion .points.csv gives the (x, y, z) of each point. We keep only the statistics window (time_step >= DEV_TIME). When the run has not happened yet (RESULTS_AVAILABLE == False) the velocity arrays stay None and every downstream cell falls back to its scaffold branch.

[ ]:

def load_plane_component(sim_cfg, comp: str):
    """Return (T, P) array of a velocity component over the plane (statistics window).

    T = number of snapshots, P = number of plane points.
    """
    plane = sim_cfg.output.series["inlet_plane"].planes["measurement"]
    df = plane.read_full_data(comp)
    df = df[df["time_step"] >= DEV_TIME].sort_values("time_step")
    point_cols = [c for c in df.columns if c != "time_step"]
    return df[point_cols].to_numpy(dtype=float)


def load_plane_points(sim_cfg):
    plane = sim_cfg.output.series["inlet_plane"].planes["measurement"]
    return pd.read_csv(plane.points_filename)


ux_t = uy_t = uz_t = points = None
if RESULTS_AVAILABLE:
    ux_t = load_plane_component(sim_cfg, "ux")
    uy_t = load_plane_component(sim_cfg, "uy")
    uz_t = load_plane_component(sim_cfg, "uz")
    points = load_plane_points(sim_cfg)
    print("snapshots, points =", ux_t.shape)
else:
    print("No results yet - scaffold mode (run the case to populate the plots).")

Sampled Reynolds-stress tensor

Remove the mean from each component (homogeneous in (y, z), so the per-component mean over all plane points and time is the inlet mean), then form the six independent covariances over the pooled (point, time) samples and compare against R_target component by component. The comparison is reported both as a table and as a measured-vs-target plot (project convention: a tabular result is always accompanied by a figure), with the shear terms Rxz/Ryz highlighted.

[ ]:

def sampled_reynolds_stress(ux_t, uy_t, uz_t):
    """Pooled covariance tensor <u'_i u'_j> over (point, time)."""
    u = np.stack([ux_t.ravel(), uy_t.ravel(), uz_t.ravel()])  # (3, N)
    up = u - u.mean(axis=1, keepdims=True)
    return (up @ up.T) / up.shape[1]


def tensor_report(R_meas, R_target):
    rows = []
    for n, (i, j) in zip(COMP_NAMES, COMP_IDX):
        tgt = R_target[i, j]
        scale = abs(tgt) if abs(tgt) > 0 else np.sqrt(R_target[i, i] * R_target[j, j])
        mea = np.nan if R_meas is None else R_meas[i, j]
        rel = np.nan if R_meas is None else 100.0 * (mea - tgt) / scale
        verdict = "-" if R_meas is None else ("PASS" if abs(rel) <= 10.0 else "FAIL")
        rows.append((n, tgt, mea, rel, verdict))
    return pd.DataFrame(rows, columns=["component", "target", "sampled", "rel_err_%", "verdict"])

[ ]:

R_meas = sampled_reynolds_stress(ux_t, uy_t, uz_t) if RESULTS_AVAILABLE else None
report = tensor_report(R_meas, R_target)
report

[ ]:

def plot_tensor(R_meas, R_target):
    """Measured vs prescribed R_ij, one marker pair per component.

    Left: signed component values (target line + measured points), so the
    negative shear Rxz and the positive Ryz are read off directly. The shaded
    band is the +-10 % acceptance window around each target.
    Right: per-component relative error with the +-10 % criterion lines.
    """
    tgt = np.array([R_target[i, j] for (i, j) in COMP_IDX])
    scale = np.array(
        [abs(t) if abs(t) > 0 else np.sqrt(R_target[i, i] * R_target[j, j])
         for t, (i, j) in zip(tgt, COMP_IDX)]
    )
    x = np.arange(len(COMP_NAMES))

    fig, ax = common.fig_double()

    # Left panel: signed component magnitudes.
    ax[0].fill_between(
        x, tgt - 0.1 * scale, tgt + 0.1 * scale, color=common.colors.exp, alpha=0.2,
        label=r"target $\pm$10%",
    )
    ax[0].plot(x, tgt, **common.markers.exp_line(linestyle="--"), label="prescribed $R_{ij}$")
    if R_meas is not None:
        mea = np.array([R_meas[i, j] for (i, j) in COMP_IDX])
        ax[0].plot(x, mea, **common.markers.sim("o"), label=r"sampled $\langle u'_i u'_j\rangle$")
    ax[0].axhline(0.0, color=common.colors.grid, linewidth=0.8)
    ax[0].set_xticks(x)
    ax[0].set_xticklabels(COMP_NAMES)
    ax[0].set_ylabel("component value (lattice units)")
    ax[0].set_title("Reynolds-stress components")
    ax[0].legend()

    # Right panel: relative error vs the +-10% criterion.
    ax[1].axhspan(-10.0, 10.0, color=common.colors.exp, alpha=0.15, label=r"$\pm$10%")
    ax[1].axhline(0.0, color=common.colors.grid, linewidth=0.8)
    if R_meas is not None:
        rel = 100.0 * (mea - tgt) / scale
        ax[1].plot(x, rel, **common.markers.sim("s"))
    else:
        ax[1].text(0.5, 0.5, "results pending\nGPU execution", ha="center", va="center",
                   transform=ax[1].transAxes, color=common.colors.text)
    ax[1].set_xticks(x)
    ax[1].set_xticklabels(COMP_NAMES)
    ax[1].set_ylabel("relative error [%]")
    ax[1].set_ylim(-30.0, 30.0)
    ax[1].set_title("Agreement with target")
    ax[1].legend()
    return fig


_ = plot_tensor(R_meas, R_target)
plt.show()

1D energy spectrum vs von Karman

Compute the one-sided streamwise velocity spectrum along the homogeneous spanwise rows of the plane, averaged over rows and snapshots, and overlay the von Karman target with the prescribed variance sigma_u^2 = Rxx and an integral length scale set by the SEM eddy size.

[ ]:

def von_karman_reduced(kL_over_U):
    """Reduced von Karman longitudinal spectrum k S_u(k) / sigma_u^2."""
    x = np.asarray(kL_over_U, dtype=float)
    return 4.0 * x / (1.0 + 70.8 * x**2) ** (5.0 / 6.0)


def streamwise_spectrum(ux_t, points):
    """1D spectrum of u' along the spanwise (y) rows, averaged over z-rows and time.

    Reshapes the flat plane points back to a (ny, nz) grid using the points CSV,
    detrends each spanwise row, FFTs along y, and averages |U_hat|^2.
    """
    y = points["y"].to_numpy(dtype=float)
    z = points["z"].to_numpy(dtype=float)
    yu = np.unique(y)
    zu = np.unique(z)
    ny, nz = len(yu), len(zu)
    dy = np.diff(yu).mean()
    iy = np.searchsorted(yu, y)
    iz = np.searchsorted(zu, z)
    psd = np.zeros(ny // 2 + 1)
    count = 0
    for t in range(ux_t.shape[0]):
        grid = np.full((ny, nz), np.nan)
        grid[iy, iz] = ux_t[t]
        for j in range(nz):
            rowj = grid[:, j]
            if np.isnan(rowj).any():
                continue
            rowj = rowj - rowj.mean()
            fhat = np.fft.rfft(rowj)
            psd += (np.abs(fhat) ** 2) / ny
            count += 1
    psd /= max(count, 1)
    k = 2.0 * np.pi * np.fft.rfftfreq(ny, d=dy)
    return k, psd


# Eddy length scale prescribed in the config (lattice units).
L_u = float(sem.eddies.lengthscale.x)
print("SEM streamwise length scale L_u =", L_u)

[ ]:

def plot_spectrum(ux_t, points):
    """Reduced streamwise spectrum vs the von Karman curve on log-log axes."""
    fig, ax = common.fig_single()
    xx = np.logspace(-2, 1, 200)
    ax.loglog(xx, von_karman_reduced(xx), **common.markers.exp_line(linestyle="--"),
              label="von Karman")
    # -5/3 inertial-range guide line.
    xg = np.logspace(0.0, 1.0, 50)
    ax.loglog(xg, 1.2 * xg ** (-2.0 / 3.0), color=common.colors.grid, linewidth=1.0,
              label=r"$-5/3$ slope")
    if ux_t is not None:
        k, psd = streamwise_spectrum(ux_t, points)
        x = k * L_u / U_ref
        reduced = (k * psd) / R_target[0, 0]
        ax.loglog(x[1:], reduced[1:], **common.markers.sim("o"), label="SEM inlet")
    else:
        ax.text(0.5, 0.1, "results pending GPU execution", ha="center", va="bottom",
                transform=ax.transAxes, color=common.colors.text)
    ax.set_xlabel(r"$k\, L_u / U$")
    ax.set_ylabel(r"$k\, S_u(k) / \sigma_u^2$")
    ax.set_title("Streamwise energy spectrum")
    ax.legend()
    return fig


_ = plot_spectrum(ux_t, points)
plt.show()

Summary

Success criteria for this case:

  • Every component of the sampled covariance tensor <u'_i u'_j> agrees with the prescribed R_ij within ~10% (all rows of the report show PASS and all markers sit inside the shaded band), in particular the shear terms Rxz (negative) and Ryz (positive).

  • The sampled streamwise turbulence intensity matches Iu = sqrt(Rxx)/U (consistency with the existing ABL diagonal-only check).

  • The 1D streamwise energy spectrum collapses onto the von Karman curve in the inertial range (a -5/3 slope rolling off into the energetic large scales near k L_u / U ~ 1).

Results pending GPU execution - the plots render in scaffold mode now and fill in automatically once the case has produced its results/ outputs (RESULTS_AVAILABLE flips to True).

Version

[ ]:

if RESULTS_AVAILABLE:
    sim_info = sim_cfg.output.read_info()
    print("commit:", sim_info.get("commit"))
    print("version:", sim_info.get("version"))
else:
    print("Version info available after the case has been run.")

Configuration

[ ]:

from IPython.display import Code

Code(filename=str(CONFIG))