{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "78f85084",
   "metadata": {},
   "source": [
    "# ABL Velocity and Turbulence Intensity from a Velocity Profile\n",
    "\n",
    "A common ABL post-processing task is to extract the mean streamwise velocity\n",
    "profile $U(z)$ and the streamwise turbulence intensity $I_u(z)$ from a line probe\n",
    "placed in the empty domain, and compare them against a target terrain category\n",
    "(e.g., EN 1991-1-4). This notebook walks through the calculation from a probe\n",
    "time series and shows how to overlay the simulated curves on the analytical\n",
    "reference.\n",
    "\n",
    "<div class=\"alert alert-info\">\n",
    "<strong>Note.</strong> This notebook assumes a vertical line probe was configured in the AeroSim\n",
    "setup interface and exported as CSV after the simulation. See the\n",
    "<a href=\"../guided-cases/ABL/setup.md\">ABL guided case setup</a> for the probe configuration and\n",
    "<a href=\"../getting-started/visualizing-results.md\">Visualizing Results</a> for the export workflow.\n",
    "</div>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "33af536d",
   "metadata": {},
   "source": [
    "## Probe Output Format\n",
    "\n",
    "A vertical line probe exports two CSV files: a points file with the coordinates\n",
    "of each probe point, and a velocity file with the time series of `ux`, `uy`,\n",
    "`uz` at each point. For this workflow, only the streamwise component `ux` is\n",
    "needed.\n",
    "\n",
    "| File | Columns | Notes |\n",
    "| --- | --- | --- |\n",
    "| `line.points.csv` | `idx, x, y, z` | One row per probe point, indexed by `idx` |\n",
    "| `line.ux.csv` | `time_step, 0, 1, 2, ...` | One row per time step; remaining columns hold `ux` at each probe point (column header is the point `idx`) |\n",
    "\n",
    "The sampling rate `fs` is implied by the time step column:\n",
    "`dt = time_step[1] - time_step[0]` and `fs = 1 / dt`.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "29b0cbe5",
   "metadata": {},
   "source": [
    "## Step 1: Load the Probe Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fc20cf11",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "points_path = \"line.line_profile line.points.csv\"\n",
    "ux_path = \"line.line_profile line.ux.csv\"\n",
    "\n",
    "df_pts = pd.read_csv(points_path, index_col=\"idx\")\n",
    "df_ux = pd.read_csv(ux_path)\n",
    "\n",
    "time_steps = df_ux[\"time_step\"].to_numpy(dtype=np.float64)\n",
    "df_ux.drop(columns=[\"time_step\"], inplace=True)\n",
    "df_ux.columns = df_ux.columns.astype(int)\n",
    "\n",
    "z = df_pts[\"z\"].to_numpy()\n",
    "dt = time_steps[1] - time_steps[0]\n",
    "print(f\"z range: {z.min():.1f} - {z.max():.1f} m  |  dt = {dt:.4f} s\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6f5d05d1",
   "metadata": {},
   "source": [
    "## Step 2: Compute Mean Velocity and Turbulence Intensity\n",
    "\n",
    "For each probe point, take the time mean and standard deviation of the `ux`\n",
    "time series. Turbulence intensity is the ratio $\\sigma_u / U$.\n",
    "\n",
    "$$\n",
    "U(z) = \\overline{u_x(z, t)}, \\qquad \\sigma_u(z) = \\sqrt{\\overline{u_x'(z, t)^2}}, \\qquad I_u(z) = \\frac{\\sigma_u(z)}{U(z)}\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d323ca47",
   "metadata": {},
   "outputs": [],
   "source": [
    "u_mean = df_ux.mean(axis=0).to_numpy()\n",
    "u_rms  = df_ux.std(axis=0).to_numpy()\n",
    "Iu     = u_rms / u_mean"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d706c14c",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-success\">\n",
    "<strong>Tip.</strong> Discard the initial transient before computing statistics. A common rule of\n",
    "thumb is to drop the first $T \\approx L_x / U_{ref}$ of the time series, where\n",
    "$L_x$ is the streamwise domain length, so that the field has been advected once\n",
    "across the domain before averaging starts.\n",
    "</div>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4051a9b0",
   "metadata": {},
   "source": [
    "## Step 3: EN 1991-1-4 Reference Profiles\n",
    "\n",
    "The Eurocode wind action standard defines five terrain categories. Each is\n",
    "parameterised by a roughness length $z_0$ and a minimum height $z_{min}$ below\n",
    "which the profile is held constant.\n",
    "\n",
    "| Category | $z_0$ (m) | $z_{min}$ (m) | Description |\n",
    "| --- | --- | --- | --- |\n",
    "| 0 | 0.003 | 1.0 | Sea, coastal areas exposed to open sea |\n",
    "| I | 0.01 | 1.0 | Lakes or flat horizontal area with negligible vegetation |\n",
    "| II | 0.05 | 2.0 | Low vegetation, isolated obstacles |\n",
    "| III | 0.30 | 5.0 | Regular cover of vegetation, suburbs, forest |\n",
    "| IV | 1.00 | 10.0 | Urban area, at least 15% covered with buildings |\n",
    "\n",
    "The mean velocity profile is\n",
    "\n",
    "$$\n",
    "\\frac{U(z)}{U_{ref}} = k_r \\, \\ln\\!\\left(\\frac{z}{z_0}\\right), \\qquad k_r = 0.19 \\left(\\frac{z_0}{z_{0,II}}\\right)^{0.07}\n",
    "$$\n",
    "\n",
    "with $z_{0,II} = 0.05$ m. Turbulence intensity is\n",
    "\n",
    "$$\n",
    "I_u(z) = \\frac{1}{\\ln(z / z_0)} \\quad \\text{for } z \\geq z_{min}\n",
    "$$\n",
    "\n",
    "clamped to $z_{min}$ below the minimum height.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dcd426d0",
   "metadata": {},
   "outputs": [],
   "source": [
    "EU_CATEGORIES = {\n",
    "    \"Cat 0\":   {\"z0\": 0.003, \"z_min\": 1.0},\n",
    "    \"Cat I\":   {\"z0\": 0.01,  \"z_min\": 1.0},\n",
    "    \"Cat II\":  {\"z0\": 0.05,  \"z_min\": 2.0},\n",
    "    \"Cat III\": {\"z0\": 0.30,  \"z_min\": 5.0},\n",
    "    \"Cat IV\":  {\"z0\": 1.00,  \"z_min\": 10.0},\n",
    "}\n",
    "Z0_REF = 0.05  # z0_II\n",
    "\n",
    "\n",
    "def eu_velocity_profile(z_arr, z0):\n",
    "    kr = 0.19 * (z0 / Z0_REF) ** 0.07\n",
    "    return kr * np.log(np.maximum(z_arr, 1e-9) / z0)\n",
    "\n",
    "\n",
    "def eu_turbulence_intensity(z_arr, z0, z_min):\n",
    "    z_eff = np.maximum(z_arr, z_min)\n",
    "    return 1.0 / np.log(z_eff / z0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b6684265",
   "metadata": {},
   "source": [
    "## Step 4: Overlay Simulation and Reference\n",
    "\n",
    "Pick the target category, normalise the simulation by the reference velocity\n",
    "used in the setup interface, and plot both curves.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "61129685",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "category = \"Cat II\"\n",
    "cat = EU_CATEGORIES[category]\n",
    "u_ref = 20.0  # reference velocity used in the AeroSim setup [m/s]\n",
    "z_eu = np.linspace(0.5, 200, 500)\n",
    "\n",
    "fig, axes = plt.subplots(1, 2, figsize=(12, 6))\n",
    "\n",
    "# Mean velocity\n",
    "axes[0].plot(u_mean / u_ref, z, label=\"AeroSim LES\", linewidth=2)\n",
    "axes[0].plot(eu_velocity_profile(z_eu, cat[\"z0\"]), z_eu,\n",
    "             \"k--\", label=f\"EN 1991-1-4 {category}\", linewidth=2)\n",
    "axes[0].set_xlabel(\"U / U_ref [-]\")\n",
    "axes[0].set_ylabel(\"z [m]\")\n",
    "axes[0].set_title(\"Mean longitudinal velocity\")\n",
    "axes[0].set_ylim(2, 200)\n",
    "axes[0].grid(True)\n",
    "axes[0].legend(loc=\"lower right\")\n",
    "\n",
    "# Turbulence intensity\n",
    "axes[1].plot(Iu, z, label=\"AeroSim LES\", linewidth=2)\n",
    "axes[1].plot(eu_turbulence_intensity(z_eu, cat[\"z0\"], cat[\"z_min\"]), z_eu,\n",
    "             \"k--\", label=f\"EN 1991-1-4 {category}\", linewidth=2)\n",
    "axes[1].set_xlabel(\"Iu [-]\")\n",
    "axes[1].set_ylabel(\"z [m]\")\n",
    "axes[1].set_title(\"Turbulence intensity\")\n",
    "axes[1].set_ylim(2, 200)\n",
    "axes[1].set_xlim(0, 0.6)\n",
    "axes[1].grid(True)\n",
    "axes[1].legend(loc=\"upper right\")\n",
    "\n",
    "fig.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "db4998d3",
   "metadata": {},
   "source": [
    "## Acceptance Criteria\n",
    "\n",
    "For the empty-domain ABL case, the simulated profiles should match the target\n",
    "across the height range of engineering interest ($z_{ref}$ to $2 z_{ref}$,\n",
    "where $z_{ref}$ is the reference height of the structure under study):\n",
    "\n",
    "| Quantity | Acceptable deviation | Common cause if exceeded |\n",
    "| --- | --- | --- |\n",
    "| $U(z) / U_{ref}$ | within $\\pm 5\\%$ of the target log law | Inlet $z_0$ and ground $z_0$ not aligned; insufficient fetch |\n",
    "| $I_u(z)$ | within $\\pm 10\\%$ of the target | Probe positioned in roughness fin wake; sampling window too short |\n",
    "\n",
    "If the deviation is larger, see the troubleshooting table in\n",
    "[Matching Custom Inlet](matching-custom-inlet.md).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e2c2796a",
   "metadata": {},
   "source": [
    "## Next Steps\n",
    "\n",
    "- [Time Normalization Using the Integral Length Scale](time-normalization-integral-lengthscale.ipynb) - extend the same probe time series to non-dimensional time analysis.\n",
    "- [ABL guided case post-processing](../guided-cases/ABL/post-processing.md) - end-to-end validation of an ABL setup.\n",
    "- [Matching Custom Inlet](matching-custom-inlet.md) - if the target profile comes from a wind tunnel rather than an EN 1991-1-4 category.\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}