{ "cells": [ { "cell_type": "markdown", "id": "cell-00", "metadata": {}, "source": [ "# SEM Inlet Reynolds-Stress Tensor and Spectrum\n", "\n", "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.\n", "\n", "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:\n", "\n", "- Each sampled covariance `` matches the prescribed `R_ij` within ~10%, with correct sign and magnitude of the shear terms.\n", "- The 1D streamwise energy spectrum follows the von Karman form in the inertial range.\n", "\n", "````{note}\n", "**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\n", "\n", "```bash\n", "uv run nassu run validation/inlet/01_sem_reynolds_stress/01_sem_reynolds_stress.nassu.yaml\n", "```\n", "\n", "(the `RESULTS_AVAILABLE` flag below switches on when the plane output exists).\n", "````" ] }, { "cell_type": "markdown", "id": "cell-01", "metadata": {}, "source": [ "## Load the configuration\n", "\n", "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." ] }, { "cell_type": "code", "execution_count": null, "id": "cell-02", "metadata": {}, "outputs": [], "source": [ "import pathlib\n", "\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "import pandas as pd\n", "\n", "import nassu.viz as common\n", "from nassu.cfg.model import ConfigScheme\n", "\n", "common.use_style()\n", "\n", "CASE_DIR = pathlib.Path(\"validation/inlet/01_sem_reynolds_stress\")\n", "CONFIG = CASE_DIR / \"01_sem_reynolds_stress.nassu.yaml\"\n", "\n", "sim_cfgs_list = ConfigScheme.sim_cfgs_from_file_dct(str(CONFIG))\n", "sim_cfg = list(sim_cfgs_list.values())[0]\n", "print(sim_cfg.name, \"n_steps =\", sim_cfg.n_steps, \"domain =\", sim_cfg.domain.domain_size)\n", "\n", "DEV_TIME = 5000 # development steps discarded from the series\n", "\n", "# Whether the case has been run: the plane XDMF exists. When False the notebook\n", "# runs as a scaffold (target + plotting code present, measured data absent).\n", "_plane = sim_cfg.output.series[\"inlet_plane\"].planes[\"measurement\"]\n", "RESULTS_AVAILABLE = _plane.xdmf_filename.exists()\n", "print(\"RESULTS_AVAILABLE =\", RESULTS_AVAILABLE)" ] }, { "cell_type": "markdown", "id": "cell-03", "metadata": {}, "source": [ "## Prescribed target tensor\n", "\n", "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." ] }, { "cell_type": "code", "execution_count": null, "id": "cell-04", "metadata": {}, "outputs": [], "source": [ "sem = sim_cfg.models.BC.inlet_turbulence\n", "profile_csv = sem.profile.csv_profile_data\n", "df_profile = pd.read_csv(profile_csv)\n", "\n", "# Constant with height: take the prescribed values at any interior row.\n", "row = df_profile.iloc[len(df_profile) // 2]\n", "U_ref = float(row[\"ux\"])\n", "R_target = np.array(\n", " [\n", " [row[\"Rxx\"], row[\"Rxy\"], row[\"Rxz\"]],\n", " [row[\"Rxy\"], row[\"Ryy\"], row[\"Ryz\"]],\n", " [row[\"Rxz\"], row[\"Ryz\"], row[\"Rzz\"]],\n", " ],\n", " dtype=float,\n", ")\n", "\n", "eigs = np.linalg.eigvalsh(R_target)\n", "A_chol = np.linalg.cholesky(R_target) # lower-triangular, A A^T = R\n", "\n", "# Component bookkeeping shared by the table and the plot below.\n", "COMP_NAMES = [\"Rxx\", \"Ryy\", \"Rzz\", \"Rxy\", \"Rxz\", \"Ryz\"]\n", "COMP_IDX = [(0, 0), (1, 1), (2, 2), (0, 1), (0, 2), (1, 2)]\n", "\n", "print(\"U_ref =\", U_ref)\n", "print(\"R_target =\\n\", R_target)\n", "print(\"eigenvalues =\", eigs, \"-> SPD:\", bool(np.all(eigs > 0)))\n", "print(\"Iu, Iv, Iw =\", np.sqrt(np.diag(R_target)) / U_ref)\n", "print(\"rho_uw =\", R_target[0, 2] / np.sqrt(R_target[0, 0] * R_target[2, 2]))\n", "print(\"Cholesky A =\\n\", A_chol)" ] }, { "cell_type": "markdown", "id": "cell-05", "metadata": {}, "source": [ "## Load the sampled inlet plane\n", "\n", "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." ] }, { "cell_type": "code", "execution_count": null, "id": "cell-06", "metadata": {}, "outputs": [], "source": [ "def load_plane_component(sim_cfg, comp: str):\n", " \"\"\"Return (T, P) array of a velocity component over the plane (statistics window).\n", "\n", " T = number of snapshots, P = number of plane points.\n", " \"\"\"\n", " plane = sim_cfg.output.series[\"inlet_plane\"].planes[\"measurement\"]\n", " df = plane.read_full_data(comp)\n", " df = df[df[\"time_step\"] >= DEV_TIME].sort_values(\"time_step\")\n", " point_cols = [c for c in df.columns if c != \"time_step\"]\n", " return df[point_cols].to_numpy(dtype=float)\n", "\n", "\n", "def load_plane_points(sim_cfg):\n", " plane = sim_cfg.output.series[\"inlet_plane\"].planes[\"measurement\"]\n", " return pd.read_csv(plane.points_filename)\n", "\n", "\n", "ux_t = uy_t = uz_t = points = None\n", "if RESULTS_AVAILABLE:\n", " ux_t = load_plane_component(sim_cfg, \"ux\")\n", " uy_t = load_plane_component(sim_cfg, \"uy\")\n", " uz_t = load_plane_component(sim_cfg, \"uz\")\n", " points = load_plane_points(sim_cfg)\n", " print(\"snapshots, points =\", ux_t.shape)\n", "else:\n", " print(\"No results yet - scaffold mode (run the case to populate the plots).\")" ] }, { "cell_type": "markdown", "id": "cell-07", "metadata": {}, "source": [ "## Sampled Reynolds-stress tensor\n", "\n", "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." ] }, { "cell_type": "code", "execution_count": null, "id": "cell-08", "metadata": {}, "outputs": [], "source": [ "def sampled_reynolds_stress(ux_t, uy_t, uz_t):\n", " \"\"\"Pooled covariance tensor over (point, time).\"\"\"\n", " u = np.stack([ux_t.ravel(), uy_t.ravel(), uz_t.ravel()]) # (3, N)\n", " up = u - u.mean(axis=1, keepdims=True)\n", " return (up @ up.T) / up.shape[1]\n", "\n", "\n", "def tensor_report(R_meas, R_target):\n", " rows = []\n", " for n, (i, j) in zip(COMP_NAMES, COMP_IDX):\n", " tgt = R_target[i, j]\n", " scale = abs(tgt) if abs(tgt) > 0 else np.sqrt(R_target[i, i] * R_target[j, j])\n", " mea = np.nan if R_meas is None else R_meas[i, j]\n", " rel = np.nan if R_meas is None else 100.0 * (mea - tgt) / scale\n", " verdict = \"-\" if R_meas is None else (\"PASS\" if abs(rel) <= 10.0 else \"FAIL\")\n", " rows.append((n, tgt, mea, rel, verdict))\n", " return pd.DataFrame(rows, columns=[\"component\", \"target\", \"sampled\", \"rel_err_%\", \"verdict\"])" ] }, { "cell_type": "code", "execution_count": null, "id": "cell-09", "metadata": {}, "outputs": [], "source": [ "R_meas = sampled_reynolds_stress(ux_t, uy_t, uz_t) if RESULTS_AVAILABLE else None\n", "report = tensor_report(R_meas, R_target)\n", "report" ] }, { "cell_type": "code", "execution_count": null, "id": "cell-10", "metadata": {}, "outputs": [], "source": [ "def plot_tensor(R_meas, R_target):\n", " \"\"\"Measured vs prescribed R_ij, one marker pair per component.\n", "\n", " Left: signed component values (target line + measured points), so the\n", " negative shear Rxz and the positive Ryz are read off directly. The shaded\n", " band is the +-10 % acceptance window around each target.\n", " Right: per-component relative error with the +-10 % criterion lines.\n", " \"\"\"\n", " tgt = np.array([R_target[i, j] for (i, j) in COMP_IDX])\n", " scale = np.array(\n", " [abs(t) if abs(t) > 0 else np.sqrt(R_target[i, i] * R_target[j, j])\n", " for t, (i, j) in zip(tgt, COMP_IDX)]\n", " )\n", " x = np.arange(len(COMP_NAMES))\n", "\n", " fig, ax = common.fig_double()\n", "\n", " # Left panel: signed component magnitudes.\n", " ax[0].fill_between(\n", " x, tgt - 0.1 * scale, tgt + 0.1 * scale, color=common.colors.exp, alpha=0.2,\n", " label=r\"target $\\pm$10%\",\n", " )\n", " ax[0].plot(x, tgt, **common.markers.exp_line(linestyle=\"--\"), label=\"prescribed $R_{ij}$\")\n", " if R_meas is not None:\n", " mea = np.array([R_meas[i, j] for (i, j) in COMP_IDX])\n", " ax[0].plot(x, mea, **common.markers.sim(\"o\"), label=r\"sampled $\\langle u'_i u'_j\\rangle$\")\n", " ax[0].axhline(0.0, color=common.colors.grid, linewidth=0.8)\n", " ax[0].set_xticks(x)\n", " ax[0].set_xticklabels(COMP_NAMES)\n", " ax[0].set_ylabel(\"component value (lattice units)\")\n", " ax[0].set_title(\"Reynolds-stress components\")\n", " ax[0].legend()\n", "\n", " # Right panel: relative error vs the +-10% criterion.\n", " ax[1].axhspan(-10.0, 10.0, color=common.colors.exp, alpha=0.15, label=r\"$\\pm$10%\")\n", " ax[1].axhline(0.0, color=common.colors.grid, linewidth=0.8)\n", " if R_meas is not None:\n", " rel = 100.0 * (mea - tgt) / scale\n", " ax[1].plot(x, rel, **common.markers.sim(\"s\"))\n", " else:\n", " ax[1].text(0.5, 0.5, \"results pending\\nGPU execution\", ha=\"center\", va=\"center\",\n", " transform=ax[1].transAxes, color=common.colors.text)\n", " ax[1].set_xticks(x)\n", " ax[1].set_xticklabels(COMP_NAMES)\n", " ax[1].set_ylabel(\"relative error [%]\")\n", " ax[1].set_ylim(-30.0, 30.0)\n", " ax[1].set_title(\"Agreement with target\")\n", " ax[1].legend()\n", " return fig\n", "\n", "\n", "_ = plot_tensor(R_meas, R_target)\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "cell-11", "metadata": {}, "source": [ "## 1D energy spectrum vs von Karman\n", "\n", "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." ] }, { "cell_type": "code", "execution_count": null, "id": "cell-12", "metadata": {}, "outputs": [], "source": [ "def von_karman_reduced(kL_over_U):\n", " \"\"\"Reduced von Karman longitudinal spectrum k S_u(k) / sigma_u^2.\"\"\"\n", " x = np.asarray(kL_over_U, dtype=float)\n", " return 4.0 * x / (1.0 + 70.8 * x**2) ** (5.0 / 6.0)\n", "\n", "\n", "def streamwise_spectrum(ux_t, points):\n", " \"\"\"1D spectrum of u' along the spanwise (y) rows, averaged over z-rows and time.\n", "\n", " Reshapes the flat plane points back to a (ny, nz) grid using the points CSV,\n", " detrends each spanwise row, FFTs along y, and averages |U_hat|^2.\n", " \"\"\"\n", " y = points[\"y\"].to_numpy(dtype=float)\n", " z = points[\"z\"].to_numpy(dtype=float)\n", " yu = np.unique(y)\n", " zu = np.unique(z)\n", " ny, nz = len(yu), len(zu)\n", " dy = np.diff(yu).mean()\n", " iy = np.searchsorted(yu, y)\n", " iz = np.searchsorted(zu, z)\n", " psd = np.zeros(ny // 2 + 1)\n", " count = 0\n", " for t in range(ux_t.shape[0]):\n", " grid = np.full((ny, nz), np.nan)\n", " grid[iy, iz] = ux_t[t]\n", " for j in range(nz):\n", " rowj = grid[:, j]\n", " if np.isnan(rowj).any():\n", " continue\n", " rowj = rowj - rowj.mean()\n", " fhat = np.fft.rfft(rowj)\n", " psd += (np.abs(fhat) ** 2) / ny\n", " count += 1\n", " psd /= max(count, 1)\n", " k = 2.0 * np.pi * np.fft.rfftfreq(ny, d=dy)\n", " return k, psd\n", "\n", "\n", "# Eddy length scale prescribed in the config (lattice units).\n", "L_u = float(sem.eddies.lengthscale.x)\n", "print(\"SEM streamwise length scale L_u =\", L_u)" ] }, { "cell_type": "code", "execution_count": null, "id": "cell-13", "metadata": {}, "outputs": [], "source": [ "def plot_spectrum(ux_t, points):\n", " \"\"\"Reduced streamwise spectrum vs the von Karman curve on log-log axes.\"\"\"\n", " fig, ax = common.fig_single()\n", " xx = np.logspace(-2, 1, 200)\n", " ax.loglog(xx, von_karman_reduced(xx), **common.markers.exp_line(linestyle=\"--\"),\n", " label=\"von Karman\")\n", " # -5/3 inertial-range guide line.\n", " xg = np.logspace(0.0, 1.0, 50)\n", " ax.loglog(xg, 1.2 * xg ** (-2.0 / 3.0), color=common.colors.grid, linewidth=1.0,\n", " label=r\"$-5/3$ slope\")\n", " if ux_t is not None:\n", " k, psd = streamwise_spectrum(ux_t, points)\n", " x = k * L_u / U_ref\n", " reduced = (k * psd) / R_target[0, 0]\n", " ax.loglog(x[1:], reduced[1:], **common.markers.sim(\"o\"), label=\"SEM inlet\")\n", " else:\n", " ax.text(0.5, 0.1, \"results pending GPU execution\", ha=\"center\", va=\"bottom\",\n", " transform=ax.transAxes, color=common.colors.text)\n", " ax.set_xlabel(r\"$k\\, L_u / U$\")\n", " ax.set_ylabel(r\"$k\\, S_u(k) / \\sigma_u^2$\")\n", " ax.set_title(\"Streamwise energy spectrum\")\n", " ax.legend()\n", " return fig\n", "\n", "\n", "_ = plot_spectrum(ux_t, points)\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "cell-14", "metadata": {}, "source": [ "## Summary\n", "\n", "Success criteria for this case:\n", "\n", "- Every component of the sampled covariance tensor `` 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).\n", "- The sampled streamwise turbulence intensity matches `Iu = sqrt(Rxx)/U` (consistency with the existing ABL diagonal-only check).\n", "- 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`).\n", "\n", "**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`)." ] }, { "cell_type": "markdown", "id": "cell-15", "metadata": {}, "source": [ "## Version" ] }, { "cell_type": "code", "execution_count": null, "id": "cell-16", "metadata": {}, "outputs": [], "source": [ "if RESULTS_AVAILABLE:\n", " sim_info = sim_cfg.output.read_info()\n", " print(\"commit:\", sim_info.get(\"commit\"))\n", " print(\"version:\", sim_info.get(\"version\"))\n", "else:\n", " print(\"Version info available after the case has been run.\")" ] }, { "cell_type": "markdown", "id": "cell-17", "metadata": {}, "source": [ "## Configuration" ] }, { "cell_type": "code", "execution_count": null, "id": "cell-18", "metadata": {}, "outputs": [], "source": [ "from IPython.display import Code\n", "\n", "Code(filename=str(CONFIG))" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "name": "python", "version": "3.10" } }, "nbformat": 4, "nbformat_minor": 5 }