zacheryasc/airfRANS-model-exploration
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions
airfRANS-model-exploration/notebooks/explore_raw_subset.ipynb

1655 lines
772 KiB
Text
Raw Normal View History

init
2026-07-21 08:32:30 +00:00
{
"cells": [
{
"cell_type": "markdown",
"id": "3df2f77e",
"metadata": {},
"source": [
"# Explore the local AirfRANS raw subset\n",
"\n",
"This notebook is a guided tour through one raw OpenFOAM-style AirfRANS simulation and the 50-simulation local subset. It imports the package utilities in `src/airfrans_frontier`; it does not train models, download data, or mutate the dataset.\n",
"\n",
"AirfRANS cases are CFD simulations around 2D aerofoils. The raw case is a solver directory, not a tidy ML table. Each file is tied to a physical object:\n",
"\n",
"- `system/`: solver controls, freestream direction, reference values, and output functions.\n",
"- `constant/polyMesh/`: mesh topology; points, faces, and named boundary patches.\n",
"- `0/`: initial field values and boundary conditions.\n",
"- `40000/`: final solved fields after the steady solver reached its last iteration.\n",
"- `postProcessing/`: histories such as drag/lift coefficients over solver iterations.\n",
"\n",
"The charts below answer three questions: which simulations are in the subset, what one selected simulation represents physically, and what the final solver fields look like on the aerofoil surface and throughout the domain.\n"
]
},
{
"cell_type": "markdown",
"id": "745e9b81",
"metadata": {},
"source": [
"## Setup\n",
"\n",
"Run from the repository root or the `notebooks/` directory. Select the repository `.venv` as the notebook kernel. If imports fail, run `uv sync --dev` from the repository root, restart VS Code's notebook kernel picker, and choose the `.venv` interpreter.\n"
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "255840f4",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"python: /home/aaron/data/airfrans/.venv/bin/python\n",
"repo: /home/aaron/data/airfrans\n",
"data: /home/aaron/data/airfrans/data/raw/OF_dataset\n"
]
}
],
"source": [
"from pathlib import Path\n",
"import gzip\n",
"import math\n",
"import re\n",
"import sys\n",
"print(f\"python: {sys.executable}\")\n",
"\n",
"\n",
"try:\n",
" import matplotlib.pyplot as plt\n",
" import numpy as np\n",
"except ModuleNotFoundError as exc:\n",
" raise ModuleNotFoundError(\n",
" f\"{exc.name!r} is not installed in this notebook kernel: {sys.executable}. \"\n",
" \"From the repository root, run `uv sync --dev`, then restart the VS Code kernel.\"\n",
" ) from exc\n",
"\n",
"\n",
"\n",
"def find_repo_root(start: Path) -> Path:\n",
" for candidate in (start, *start.parents):\n",
" if (candidate / \"pyproject.toml\").exists() and (candidate / \"src\" / \"airfrans_frontier\").exists():\n",
" return candidate\n",
" raise RuntimeError(f\"Could not find repository root from {start}\")\n",
"\n",
"\n",
"REPO_ROOT = find_repo_root(Path.cwd())\n",
"SRC_DIR = REPO_ROOT / \"src\"\n",
"if str(SRC_DIR) not in sys.path:\n",
" sys.path.insert(0, str(SRC_DIR))\n",
"\n",
"from airfrans_frontier.paths import DEFAULT_RAW_DATA_DIR, DEFAULT_RAW_MANIFEST_PATH\n",
"from airfrans_frontier.raw.inspect import format_raw_inspection, inspect_raw_subset\n",
"from airfrans_frontier.raw.manifest import load_raw_subset_manifest\n",
"\n",
"DATA_DIR = REPO_ROOT / DEFAULT_RAW_DATA_DIR\n",
"MANIFEST_PATH = REPO_ROOT / DEFAULT_RAW_MANIFEST_PATH\n",
"print(f\"repo: {REPO_ROOT}\")\n",
"print(f\"data: {DATA_DIR}\")\n"
]
},
{
"cell_type": "markdown",
"id": "a920e5bc",
"metadata": {},
"source": [
"## Verify the local subset\n",
"\n",
"This is the same package-backed check as the CLI. If this fails, fix local data before interpreting plots.\n"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "792f49bd",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"AirfRANS raw subset\n",
"status: ok\n",
"data_dir: /home/aaron/data/airfrans/data/raw/OF_dataset\n",
"manifest: /home/aaron/data/airfrans/data/raw/OF_dataset_subset_manifest.json\n",
"simulations: 50 / 50\n",
"files: 5820 / 5820\n",
"bytes: 7565631222 / 7565631222 (7.57 GB)\n",
"source_zip_bytes: 71310335730\n",
"sample_seed: 20260719\n",
"sample_method: uniform random sample without replacement from the 1000 top-level OF_dataset simulation directories, using sorted names as the sampling population\n",
"sample_simulations:\n",
"- airFoil2D_SST_32.137_12.122_4.854_5.202_9.247\n",
"- airFoil2D_SST_33.238_-0.071_0.667_2.065_6.479\n",
"- airFoil2D_SST_33.816_-1.984_6.59_0.0_7.983\n",
"- airFoil2D_SST_33.846_6.261_0.922_3.023_1.0_12.36\n",
"- airFoil2D_SST_36.155_8.69_3.296_7.636_1.0_8.571\n",
"- airFoil2D_SST_38.797_3.751_1.748_7.04_0.0_19.453\n",
"- airFoil2D_SST_40.175_10.442_1.962_3.054_0.0_6.878\n",
"- airFoil2D_SST_40.585_14.545_0.315_2.716_9.25\n",
"... 42 more\n"
]
}
],
"source": [
"report = inspect_raw_subset(DATA_DIR, MANIFEST_PATH)\n",
"print(format_raw_inspection(report, sample_limit=8))\n"
]
},
{
"cell_type": "markdown",
"id": "1aad48db",
"metadata": {},
"source": [
"## Subset-level simulation conditions\n",
"\n",
"Simulation names encode the sampled run conditions and NACA shape parameters. This section turns those names into columns so we can see what local cases are available before opening a single case.\n",
"\n",
"Name pattern used here:\n",
"\n",
"`airFoil2D_<turbulence>_<U_inf>_<alpha>_<naca_a>_<naca_b>_<naca_c>`\n",
"\n",
"Meaning of the parsed values:\n",
"\n",
"- `turbulence`: turbulence closure family used to generate the case.\n",
"- `U_inf`: freestream speed in m/s. Higher values increase the Reynolds number when viscosity is fixed.\n",
"- `alpha`: angle of attack in degrees. Positive/negative values rotate the incoming flow relative to the aerofoil and usually change lift sign/magnitude.\n",
"- `naca_a`, `naca_b`, `naca_c`: NACA shape parameters encoded by the dataset. Treat them as geometry descriptors: they identify aerofoil shape variation, not solver outputs.\n",
"\n",
"The histograms show how many local cases occupy each range of `U_inf` and `alpha`; the scatter plot shows whether speed and angle are sampled independently or clustered in this subset.\n"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "30001151",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"simulations: 50\n",
"U_inf range: 32.137 .. 93.213 m/s\n",
"alpha range: -2.718 .. 14.794 deg\n",
"First parsed rows show one case per simulation directory:\n"
]
},
{
"data": {
"text/plain": [
"[{'name': 'airFoil2D_SST_32.137_12.122_4.854_5.202_9.247',\n",
" 'turbulence': 'SST',\n",
" 'u_inf': 32.137,\n",
" 'alpha': 12.122,\n",
" 'naca_a': 4.854,\n",
" 'naca_b': 5.202,\n",
" 'naca_c': 9.247},\n",
" {'name': 'airFoil2D_SST_33.238_-0.071_0.667_2.065_6.479',\n",
" 'turbulence': 'SST',\n",
" 'u_inf': 33.238,\n",
" 'alpha': -0.071,\n",
" 'naca_a': 0.667,\n",
" 'naca_b': 2.065,\n",
" 'naca_c': 6.479},\n",
" {'name': 'airFoil2D_SST_33.816_-1.984_6.59_0.0_7.983',\n",
" 'turbulence': 'SST',\n",
" 'u_inf': 33.816,\n",
" 'alpha': -1.984,\n",
" 'naca_a': 6.59,\n",
" 'naca_b': 0.0,\n",
" 'naca_c': 7.983},\n",
" {'name': 'airFoil2D_SST_33.846_6.261_0.922_3.023_1.0_12.36'},\n",
" {'name': 'airFoil2D_SST_36.155_8.69_3.296_7.636_1.0_8.571'}]"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"manifest = load_raw_subset_manifest(MANIFEST_PATH)\n",
"sim_names = list(manifest.simulation_names)\n",
"\n",
"# Simulation directory names carry both operating conditions and geometry parameters.\n",
"# The named capture groups below become explicit variables instead of leaving the\n",
"# notebook reader to decode a long string by eye.\n",
"SIM_RE = re.compile(\n",
" r\"^airFoil2D_(?P<turbulence>[^_]+)_\" # solver/turbulence family label\n",
" r\"(?P<u_inf>-?\\d+(?:\\.\\d+)?)_\" # freestream speed U_inf [m/s]\n",
" r\"(?P<alpha>-?\\d+(?:\\.\\d+)?)_\" # angle of attack alpha [degrees]\n",
" r\"(?P<naca_a>-?\\d+(?:\\.\\d+)?)_\" # first encoded NACA geometry parameter\n",
" r\"(?P<naca_b>-?\\d+(?:\\.\\d+)?)_\" # second encoded NACA geometry parameter\n",
" r\"(?P<naca_c>-?\\d+(?:\\.\\d+)?)$\" # third encoded NACA geometry parameter\n",
")\n",
"\n",
"\n",
"def parse_sim_name(name: str) -> dict[str, float | str]:\n",
" \"\"\"Parse one AirfRANS simulation directory name into readable columns.\"\"\"\n",
" match = SIM_RE.match(name)\n",
" if not match:\n",
" # Keep unparsed names visible rather than dropping them silently.\n",
" return {\"name\": name}\n",
"\n",
" row: dict[str, float | str] = {\"name\": name, \"turbulence\": match.group(\"turbulence\")}\n",
"\n",
" # Convert numeric groups to floats so range checks, histograms, and scatter\n",
" # plots operate on real values rather than lexicographic strings.\n",
" for key in [\"u_inf\", \"alpha\", \"naca_a\", \"naca_b\", \"naca_c\"]:\n",
" row[key] = float(match.group(key))\n",
" return row\n",
"\n",
"\n",
"sim_meta = [parse_sim_name(name) for name in sim_names]\n",
"u_inf = np.array([row[\"u_inf\"] for row in sim_meta if \"u_inf\" in row], dtype=float)\n",
"alpha = np.array([row[\"alpha\"] for row in sim_meta if \"alpha\" in row], dtype=float)\n",
"\n",
"print(f\"simulations: {len(sim_names)}\")\n",
"print(f\"U_inf range: {u_inf.min():.3f} .. {u_inf.max():.3f} m/s\")\n",
"print(f\"alpha range: {alpha.min():.3f} .. {alpha.max():.3f} deg\")\n",
"print(\"First parsed rows show one case per simulation directory:\")\n",
"sim_meta[:5]\n"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "15776396",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAABW4AAAGMCAYAAABK9zuEAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjExLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvctoD+AAAAAlwSFlzAAAPYQAAD2EBqD+naQAA1c9JREFUeJzs3XdYU+f7P/B3DoEAsgUVBCeiVal7VRH3HnXPOqrW2WrVOro+1fqt2vqptbVV66qtraPOWke1dVftcFStouJmKIJsWeE8vz/8kY8xAZIQSEjer+viusiT55xz3+cc8iQ3J89RCCEEiIiIiIiIiIiIiMhqSJYOgIiIiIiIiIiIiIi0sXBLREREREREREREZGVYuCUiIiIiIiIiIiKyMizcEhEREREREREREVkZFm6JiIiIiIiIiIiIrAwLt0RERERERERERERWhoVbIiIiIiIiIiIiIivDwi0RERERERERERGRlWHhloiIiIokOjoaK1euRHR0tKVDKZKi5mHM8o8ePcLKlStx584dk7ZVUusk01nD8SiuGKwht+dlZGRg5cqVuHTpkqbNGuMkIiIiMgYLt0RERDYiJSUFK1euzPfnl19+KZbtXr16FRMnTsTVq1eLZf0lxZQ8Tp06hU2bNhm9/N27dzFx4kRcuHDB1HBLZJ1kOms4HkWJ4cGDB1i5ciXu3btn1vUWl+TkZEycOBG//fabps2ScR47dgwbNmwo8e0SERGRbVFaOgAiIiIyj7i4OEycOBH169dHs2bNLB2OXZg+fTpq1KiBIUOGGLVcuXLlMH78eFStWtVssRTHOsl0pf14REZGYuLEidizZw8qVaqk9Vxpyc2ScX733XfYtWsXRo4cWeLbJiIiItvBwi0REZGN6dy5MxYtWmTpMGzegwcP8Oeff2LmzJlGL1upUiWsXLnSrPEUxzrJdLZ8PEpLbqUlTiIiIqL8sHBLRERkh9LT0/Hdd99pHqtUKlSrVg0tW7aEUqn79iArKwu///477t27h4CAALz00ktwc3PT6Zebm4tff/0VsbGxqF+/PurXr19oLNHR0dizZw969uyJcuXK4dChQ0hMTESbNm1QsWJFAEBOTg5+/fVXxMXFoWXLlggODjYpn2e3Vb58eRw+fBh3795F//79843vypUrOH78OOrXr4/mzZtr2vfs2QMnJyd07txZZ5mcnBwcPHgQ8fHxaNSoEerWrav1/KNHj7B9+3Z06dIFVapU0YmtQoUKRu9HfevMi+XEiROIioqCn58fGjduDD8/v0LXBwCPHz/G8ePHkZiYiMDAQISFhcHZ2VnzvCkxZ2ZmauLx8vJCmzZt4O3tXWAcBR03JyenQo+9Wq3G2rVr0bRpUzRo0EDT98cff0RCQgJGjBgBV1dXAE+vXN+xYwc6d+5c4JWahe3Xwo6xqef6tWvXcOTIEQwePBheXl6a9ry/gbCwMNSpUyffuA35W7lx4wZ2794NANi/fz+ioqIAAC1atEC9evXyPdeA4jln9Ll79y5OnDgBFxcXdOrUSW+fwo7B8+dS3nkYHR2NU6dOIS0tDcHBwWjZsiUkSXeGudjYWJw+fRppaWmoU6cOGjVqBAD4+eefERERgczMTK3C8ahRo7T2xb///ovz589DCIEGDRrovE4cPnwYMTExGD58OB4+fIijR4/C0dERvr6+iIiIwJgxY+Dg4KC1zK1bt3Dw4EF0794dQUFBBu9PIiIisk4s3BIREdmhnJwcrXkf84otrq6uOHToEKpXr655bu/evRg7dixUKhWaNWuGlJQUTJ48GStXrkT79u01/RISEtCxY0d4eXkhLS0No0ePxty5c/HRRx8VGEve3LDu7u5YtWoVypcvj6ioKIwbNw4///wzatasiQEDBiAgIACxsbEYO3YsNm7ciEGDBhmdz7PbWrFiBfz8/HDz5k20aNFCb2z79u3D4MGD0bFjR4wYMULrud27d6Ndu3Zwd3fXan/w4AFat26NgIAApKSk4NVXX8Vrr72Gr776CgqFAsD/5t7cuXOnpqCUF5u3tzdWrVpl9H7Ut84rV66gU6dOKFOmDJo1a4bk5GRcunQJ06ZNwxtvvFHg+latWoXp06cjJCQEL7zwAk6dOgW1Wo0tW7agZcuWJsW8f/9+jB49Gu7u7mjatClu376NkSNHYu3atRgwYEC+sRR03AIDAws99kqlEv/9739Rr149/PjjjwCAJ0+e4JVXXkFWVhYqVaqEbt26AQB27tyJiRMn4tatW/nGY8h+LegYF+VcP336NCZOnIg2bdpoFW4TExMxceJEfPHFFwUWbg35W0lMTMSNGzcAALdv30Zubi4AaIrI+nIDiuec0WfJkiWYO3cuGjVqhCpVqmD+/PlYvHixTr/CjsHz55KHhwemT5+OFStWoFWrVqhQoQLee+89lC9fHj/99JOmuK5WqzX9mjZtimrVqmH58uUoU6YM9uzZg5s3byIhIQG5ubla+zpvP6alpeGVV17BgQMH0L59eygUCowfPx4dOnTA999/r3lNWbduHY4ePQpXV1fMmzcPtWrVQnx8PF5//XWMHz8e5cuXR+/evbVynj9/PrZt24Zhw4YVuh+JiIioFBBERERkE27cuCEAiM6dO4sVK1bo/Fy8eLHA5Z88eSKaNGkiOnbsqGm7cOGCcHJyEsOGDROZmZma9ri4OHH69GkhhBCHDh0SAETjxo3FzZs3NX3efvtt4eDgICIjIwvcbt7yjRo1Enfv3hVCCCHLsujatasICQkRQ4cOFbdv39b079mzpwgKChLZ2dlG55O3rRdffFFcv35dCCFERkaGiI+P1zx36NAhIYQQn3/+uXBwcBBvvfWWkGVZa91paWnC2dlZrFy5Umfd9erVEzdu3NC0b9q0SQAQK1as0LT99ddfAoDYuXOnzvKm7kd96+zdu7eoV6+eyMnJ0bRlZGSI/fv3F7iu48ePC4VCIaZOnaqVc+vWrYWPj4+Ij483OuZLly4JlUolRo8erRXPBx98IJycnERERES+8RR03PTRd+wnTZokvL29RW5urhBCiP379wtJkkTNmjW18uzXr5+oXr16gfvHkP1a0DEuyrm+fv16AUBcvXpVK6b79+8LAOKLL74oMAZ99O2vEydOCABiz549Ov31rbc4zhl9jhw5IgCI999/X9MWHx8vOnXqJACIpUuXFhhnQedSXgy//PKLpn9SUpKoV6+eCA8P17TNmjVLSJIkdu/erRXbsWPHRFJSkhBCiDFjxoiyZcvqzWHMmDHC0dFRnDlzRitWlUolRowYoWkbNmyY8PDwEKNGjdKcA/fv3xdqtVpUqlRJdOrUSWu9CQkJwtnZWYwePbqgXUhERESliO53foiIiKhUe/jwIS5cuKDz8+jRI52+//77LzZt2oSvv/4aGzZsQLly5XDy5EnIsgwA+PTTT+Hg4IDPP/8cKpVKs5yfn5/WtAEA0KNHD1SrVk3z+NVXX0Vubi6OHTtmUNy9evXS3ARJoVBgyJAhuH79OqpXr651Vd+QIUNw//59XL9+3eh88nTt2hU1atQAADg7O6Ns2bKa53JzczFlyhTMmDEDq1atwscff6y5UjbPL7/8gqysLPTs2VMnhu7du2t9vX3w4MFo3LgxvvzyS4P2Q1H347MSEhKgUCgghNC0OTs7o0uXLgUu9+WXX8LV1RUffvihpq1MmTJYuHAhHj9+jI0bNxod89KlSyGEwLJly7Smr5g7dy5UKhXWrVtXaD4FHbfCjn3Hjh2RmJiIv/76CwBw6NAhNGrUCAMGDMChQ4cAALIs4/Dhw+jYsWOBcZi6X/OY41wvKkP/VgxVHOeMPitXroSPjw/efvttTVvZsmV1rogvzPPnUpkyZfDZZ59h4MCBWlMveHp6Ys6cOTh27BiuXr2KjIwMfP755xgwYAB69eqltc7WrVvD09OzwO2mpqZiw4YNGDJkiNZNJBs3bowRI0bg+++/R2JioqY9JSUFs2fPhqOjIwAgMDAQDg4OeO2113Do0CFERkZq+q5fvx6ZmZkYM2aMUfuCiIiIrBenSiAiIrIxhtycLDExEX369MFff/2l+Vq
"text/plain": [
"<Figure size 1400x380 with 3 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fig, axes = plt.subplots(1, 3, figsize=(14, 3.8))\n",
"axes[0].hist(u_inf, bins=12, edgecolor=\"white\")\n",
"axes[0].set_title(\"Freestream speed distribution\")\n",
"axes[0].set_xlabel(\"U_inf [m/s]; imposed incoming speed\")\n",
"axes[0].set_ylabel(\"number of simulation cases\")\n",
"axes[0].grid(alpha=0.25)\n",
"\n",
"axes[1].hist(alpha, bins=12, edgecolor=\"white\")\n",
"axes[1].set_title(\"Angle-of-attack distribution\")\n",
"axes[1].set_xlabel(\"alpha [deg]; inflow angle relative to aerofoil\")\n",
"axes[1].set_ylabel(\"number of simulation cases\")\n",
"axes[1].grid(alpha=0.25)\n",
"\n",
"axes[2].scatter(u_inf, alpha, s=28, alpha=0.8)\n",
"axes[2].set_title(\"Local subset coverage\")\n",
"axes[2].set_xlabel(\"U_inf [m/s]\")\n",
"axes[2].set_ylabel(\"alpha [deg]\")\n",
"axes[2].grid(alpha=0.25)\n",
"fig.suptitle(\"Each mark/bin is one raw simulation directory\", y=1.03)\n",
"fig.tight_layout()\n"
]
},
{
"cell_type": "markdown",
"id": "619e4a24",
"metadata": {},
"source": [
"### Reading the subset chart\n",
"\n",
"- Left histogram: tall bars mean many simulations share a similar freestream speed. This is input coverage, not a performance metric.\n",
"- Middle histogram: bars count cases by angle of attack. Angles far from zero are more aggressive flow conditions and may show stronger lift, separation, or numerical difficulty.\n",
"- Right scatter: each point is one case. A rectangular cloud would mean broad coverage of speed/angle combinations; diagonal bands or clusters would mean the subset only samples certain combinations.\n",
"\n",
"These plots do not say whether a simulation is accurate or converged. They only describe the local subset's operating-condition coverage.\n"
]
},
{
"cell_type": "markdown",
"id": "30b94d50",
"metadata": {},
"source": [
"## Pick one simulation to understand\n",
"\n",
"Start with index `0`, then change `SIM_INDEX` and rerun cells below. The selected case is a complete OpenFOAM run directory.\n"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "5ca9a3fd",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"airFoil2D_SST_32.137_12.122_4.854_5.202_9.247\n",
"/home/aaron/data/airfrans/data/raw/OF_dataset/airFoil2D_SST_32.137_12.122_4.854_5.202_9.247\n"
]
},
{
"data": {
"text/plain": [
"{'name': 'airFoil2D_SST_32.137_12.122_4.854_5.202_9.247',\n",
" 'turbulence': 'SST',\n",
" 'u_inf': 32.137,\n",
" 'alpha': 12.122,\n",
" 'naca_a': 4.854,\n",
" 'naca_b': 5.202,\n",
" 'naca_c': 9.247}"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"SIM_INDEX = 0\n",
"SIM_NAME = sim_names[SIM_INDEX]\n",
"SIM_DIR = DATA_DIR / SIM_NAME\n",
"print(SIM_NAME)\n",
"print(SIM_DIR)\n",
"parse_sim_name(SIM_NAME)\n"
]
},
{
"cell_type": "markdown",
"id": "5af214c1",
"metadata": {},
"source": [
"## OpenFOAM parsing helpers\n",
"\n",
"These helpers read ASCII OpenFOAM files, including `.gz` files. They intentionally parse only the pieces used for exploration: dictionary assignments, list fields, force coefficient tables, and selected boundary faces.\n",
"\n",
"OpenFOAM list files often look like:\n",
"\n",
"```text\n",
"<number of entries>\n",
"(\n",
"(value0 value1 value2)\n",
"...\n",
")\n",
"```\n",
"\n",
"The parsers below first find the declared row count, then collect numeric rows between parentheses. That count check is important: if the file format changes or we start reading the wrong section, the notebook should fail loudly instead of plotting misaligned data.\n"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "44dd2160",
"metadata": {},
"outputs": [],
"source": [
"FLOAT_RE = re.compile(r\"[-+]?(?:\\d+(?:\\.\\d*)?|\\.\\d+)(?:[eE][-+]?\\d+)?\")\n",
"INT_RE = re.compile(r\"\\d+\")\n",
"\n",
"\n",
"def open_text(path: Path):\n",
" \"\"\"Open plain text or gzip-compressed OpenFOAM text files.\"\"\"\n",
" if path.suffix == \".gz\":\n",
" return gzip.open(path, \"rt\", errors=\"replace\")\n",
" return path.open(\"rt\", errors=\"replace\")\n",
"\n",
"\n",
"def read_text(path: Path, max_bytes: int = 200_000) -> str:\n",
" \"\"\"Read a bounded preview so huge field files do not overwhelm the notebook.\"\"\"\n",
" if path.suffix == \".gz\":\n",
" with gzip.open(path, \"rb\") as stream:\n",
" data = stream.read(max_bytes)\n",
" else:\n",
" data = path.read_bytes()[:max_bytes]\n",
" return data.decode(\"utf-8\", errors=\"replace\")\n",
"\n",
"\n",
"def assignment(text: str, key: str) -> str | None:\n",
" \"\"\"Return the raw value from an OpenFOAM dictionary line like `key value;`.\"\"\"\n",
" match = re.search(rf\"^\\s*{re.escape(key)}\\s+([^;]+);\", text, flags=re.MULTILINE)\n",
" return match.group(1).strip() if match else None\n",
"\n",
"\n",
"def vector_assignment(text: str, key: str) -> np.ndarray | None:\n",
" \"\"\"Return an OpenFOAM vector dictionary value as a numeric numpy array.\"\"\"\n",
" value = assignment(text, key)\n",
" if value is None:\n",
" return None\n",
" numbers = [float(item) for item in FLOAT_RE.findall(value)]\n",
" return np.array(numbers, dtype=float)\n",
"\n",
"\n",
"def parse_foam_list(path: Path, columns: int) -> np.ndarray:\n",
" \"\"\"Parse a numeric OpenFOAM list into an array with the requested columns.\"\"\"\n",
" expected_count: int | None = None\n",
" in_values = False\n",
" rows: list[list[float]] = []\n",
"\n",
" with open_text(path) as stream:\n",
" for line in stream:\n",
" stripped = line.strip()\n",
" if not in_values:\n",
" # First standalone integer is the OpenFOAM-declared list length.\n",
" if expected_count is None and stripped.isdigit():\n",
" expected_count = int(stripped)\n",
" continue\n",
" # Values begin on the line containing only `(` after the count.\n",
" if expected_count is not None and stripped == \"(\":\n",
" in_values = True\n",
" continue\n",
" continue\n",
"\n",
" # A line containing only `)` ends the list.\n",
" if stripped == \")\":\n",
" break\n",
"\n",
" # Field rows may be scalar (`1.23`) or vector-like (`(1 2 3)`).\n",
" # The regex strips OpenFOAM punctuation and keeps the numeric payload.\n",
" numbers = [float(item) for item in FLOAT_RE.findall(stripped)]\n",
" if len(numbers) >= columns:\n",
" rows.append(numbers[:columns])\n",
"\n",
" array = np.array(rows, dtype=float)\n",
" if expected_count is not None and len(array) != expected_count:\n",
" raise ValueError(f\"{path}: parsed {len(array)} rows, expected {expected_count}\")\n",
" if columns == 1:\n",
" return array.reshape(-1)\n",
" return array\n",
"\n",
"\n",
"def parse_boundary(path: Path) -> dict[str, dict[str, int | str]]:\n",
" \"\"\"Parse named boundary patches and their face ranges from polyMesh/boundary.\"\"\"\n",
" text = read_text(path, max_bytes=500_000)\n",
" patches: dict[str, dict[str, int | str]] = {}\n",
" for name, body in re.findall(r\"\\n\\s*([A-Za-z][A-Za-z0-9_]*)\\s*\\n\\s*\\{(.*?)\\n\\s*\\}\", text, flags=re.DOTALL):\n",
" patch_type = assignment(body, \"type\") or \"\"\n",
" n_faces = assignment(body, \"nFaces\")\n",
" start_face = assignment(body, \"startFace\")\n",
" if n_faces is not None and start_face is not None:\n",
" patches[name] = {\n",
" \"type\": patch_type,\n",
" \"nFaces\": int(n_faces),\n",
" \"startFace\": int(start_face),\n",
" }\n",
" return patches\n",
"\n",
"\n",
"def parse_faces(path: Path, start_face: int, n_faces: int) -> list[list[int]]:\n",
" \"\"\"Read only the face definitions belonging to one named boundary patch.\"\"\"\n",
" expected_count: int | None = None\n",
" in_values = False\n",
" face_index = -1\n",
" selected: list[list[int]] = []\n",
" stop_face = start_face + n_faces\n",
"\n",
" with open_text(path) as stream:\n",
" for line in stream:\n",
" stripped = line.strip()\n",
" if not in_values:\n",
" if expected_count is None and stripped.isdigit():\n",
" expected_count = int(stripped)\n",
" continue\n",
" if expected_count is not None and stripped == \"(\":\n",
" in_values = True\n",
" continue\n",
" continue\n",
"\n",
" if stripped == \")\":\n",
" break\n",
" face_index += 1\n",
" if face_index < start_face:\n",
" continue\n",
" if face_index >= stop_face:\n",
" break\n",
"\n",
" # OpenFOAM face row format is `N(v0 v1 ... vN)`. The first integer is\n",
" # the number of vertices; the remaining integers index into `points`.\n",
" values = [int(item) for item in INT_RE.findall(stripped)]\n",
" if not values:\n",
" continue\n",
" selected.append(values[1:])\n",
"\n",
" if len(selected) != n_faces:\n",
" raise ValueError(f\"{path}: parsed {len(selected)} selected faces, expected {n_faces}\")\n",
" return selected\n",
"\n",
"\n",
"def load_force_coefficients(path: Path) -> tuple[list[str], np.ndarray]:\n",
" \"\"\"Load coefficient.dat and keep the header names aligned with data columns.\"\"\"\n",
" columns: list[str] = []\n",
" with path.open(\"rt\", errors=\"replace\") as stream:\n",
" for line in stream:\n",
" if line.startswith(\"# Time\"):\n",
" columns = line[1:].split()\n",
" break\n",
" data = np.loadtxt(path, comments=\"#\")\n",
" return columns, data\n",
"\n",
"\n",
"def summarize_array(name: str, values: np.ndarray) -> dict[str, float | int | str]:\n",
" \"\"\"Return robust distribution landmarks for large field arrays.\"\"\"\n",
" finite = values[np.isfinite(values)]\n",
" return {\n",
" \"name\": name,\n",
" \"count\": int(values.size),\n",
" \"finite\": int(finite.size),\n",
" \"min\": float(np.min(finite)),\n",
" \"p01\": float(np.percentile(finite, 1)),\n",
" \"mean\": float(np.mean(finite)),\n",
" \"p99\": float(np.percentile(finite, 99)),\n",
" \"max\": float(np.max(finite)),\n",
" }\n"
]
},
{
"cell_type": "markdown",
"id": "cc795a71",
"metadata": {},
"source": [
"## What this simulation says about the run\n",
"\n",
"Read solver setup and physical/run metadata from OpenFOAM dictionaries instead of guessing from filenames only.\n",
"\n",
"Important values printed below:\n",
"\n",
"- `solver`: OpenFOAM application used. `simpleFoam` is a steady incompressible RANS solver, so the time axis in coefficient plots is an iteration count, not physical seconds.\n",
"- `turbulence model`: closure used for Reynolds-averaged turbulence terms.\n",
"- `Uinf`: freestream speed used by force-coefficient normalization.\n",
"- `nu`: kinematic viscosity in $m^2/s$.\n",
"- `Re for lRef=1`: Reynolds number computed as `Re = U_inf * L / nu` with `L=1`. Larger Reynolds numbers generally mean inertia dominates viscosity more strongly.\n",
"- `dragDir` / `liftDir`: unit directions used to project total surface force into drag and lift coefficients.\n",
"- `angle from dragDir`: check that the solver dictionary agrees with the angle encoded in the simulation name.\n"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "f52501cd",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"solver: simpleFoam\n",
"turbulence model: kOmegaSST\n",
"Uinf: 32.137 m/s\n",
"nu: 1.560e-05 m^2/s\n",
"Re for lRef=1: 2.060e+06\n",
"dragDir: [0.97770268 0.20999399 0. ]\n",
"liftDir: [-0.20999399 0.97770268 0. ]\n",
"angle from dragDir: 12.122 deg\n"
]
}
],
"source": [
"control_text = read_text(SIM_DIR / \"system\" / \"controlDict\")\n",
"transport_text = read_text(SIM_DIR / \"constant\" / \"transportProperties\")\n",
"turbulence_text = read_text(SIM_DIR / \"constant\" / \"turbulenceProperties\")\n",
"\n",
"# These variables are read from solver dictionaries, not the filename. That makes\n",
"# this cell the authoritative check for the selected case's physical setup.\n",
"u_inf_config = float(assignment(control_text, \"Uinf\"))\n",
"nu = float(assignment(transport_text, \"nu\"))\n",
"application = assignment(control_text, \"application\")\n",
"turbulence_model = assignment(turbulence_text, \"RASModel\")\n",
"drag_dir = vector_assignment(control_text, \"dragDir\")\n",
"lift_dir = vector_assignment(control_text, \"liftDir\")\n",
"alpha_from_drag = math.degrees(math.atan2(drag_dir[1], drag_dir[0])) if drag_dir is not None else float(\"nan\")\n",
"re_lref1 = u_inf_config / nu\n",
"\n",
"print(f\"solver: {application}\")\n",
"print(f\"turbulence model: {turbulence_model}\")\n",
"print(f\"Uinf: {u_inf_config:.3f} m/s\")\n",
"print(f\"nu: {nu:.3e} m^2/s\")\n",
"print(f\"Re for lRef=1: {re_lref1:.3e}\")\n",
"print(f\"dragDir: {drag_dir}\")\n",
"print(f\"liftDir: {lift_dir}\")\n",
"print(f\"angle from dragDir: {alpha_from_drag:.3f} deg\")\n"
]
},
{
"cell_type": "markdown",
"id": "19b2fbea",
"metadata": {},
"source": [
"## Directory organization\n",
"\n",
"This shows the case as OpenFOAM organizes it: mesh under `constant/polyMesh`, solver dictionaries under `system`, initial fields under `0`, final fields under `40000`, and time histories under `postProcessing` / `logs`.\n",
"\n",
"The file counts and sizes help separate small configuration files from large arrays. Large files in `40000/` and `constant/polyMesh/` are usually numeric fields or mesh topology; small files in `system/` are mostly human-readable dictionaries.\n"
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "f776bc17",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[('0', 'dir', '7.45 MB', 13),\n",
" ('0.orig', 'dir', '9.61 KB', 8),\n",
" ('40000', 'dir', '42.52 MB', 29),\n",
" ('SST_32.137_12.122_(4.854, 5.202, 9.247).foam', 'file', '0 B', 1),\n",
" ('coef_convergence.png', 'file', '169.98 KB', 1),\n",
" ('constant', 'dir', '16.03 MB', 8),\n",
" ('log.blockMesh', 'file', '2.96 KB', 1),\n",
" ('log.checkMesh', 'file', '3.36 KB', 1),\n",
" ('log.decomposePar', 'file', '7.67 KB', 1),\n",
" ('log.foamLog', 'file', '607 B', 1),\n",
" ('log.foamToVTK', 'file', '1.75 KB', 1),\n",
" ('log.reconstructPar', 'file', '2.05 KB', 1),\n",
" ('log.simpleFoam', 'file', '84.13 MB', 1),\n",
" ('logs', 'dir', '19.75 MB', 35),\n",
" ('naca_(4.854, 5.202, 9.247).png', 'file', '26.79 KB', 1),\n",
" ('postProcessing', 'dir', '18.44 MB', 5),\n",
" ('residuals.png', 'file', '227.71 KB', 1),\n",
" ('system', 'dir', '395.57 KB', 7)]"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"def format_bytes(size: int) -> str:\n",
" if size >= 1_000_000_000:\n",
" return f\"{size / 1_000_000_000:.2f} GB\"\n",
" if size >= 1_000_000:\n",
" return f\"{size / 1_000_000:.2f} MB\"\n",
" if size >= 1_000:\n",
" return f\"{size / 1_000:.2f} KB\"\n",
" return f\"{size} B\"\n",
"\n",
"\n",
"def child_rows(path: Path) -> list[tuple[str, str, str, int]]:\n",
" rows = []\n",
" for child in sorted(path.iterdir(), key=lambda item: item.name):\n",
" if child.is_dir():\n",
" files = [item for item in child.rglob(\"*\") if item.is_file()]\n",
" size = sum(item.stat().st_size for item in files)\n",
" rows.append((child.name, \"dir\", format_bytes(size), len(files)))\n",
" else:\n",
" rows.append((child.name, \"file\", format_bytes(child.stat().st_size), 1))\n",
" return rows\n",
"\n",
"\n",
"child_rows(SIM_DIR)\n"
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "bead87a2",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"0 files= 13 size=7.45 MB\n",
"40000 files= 29 size=42.52 MB\n",
"constant files= 8 size=16.03 MB\n",
"system files= 7 size=395.57 KB\n",
"postProcessing files= 5 size=18.44 MB\n",
"logs files= 35 size=19.75 MB\n"
]
}
],
"source": [
"for folder in [\"0\", \"40000\", \"constant\", \"system\", \"postProcessing\", \"logs\"]:\n",
" path = SIM_DIR / folder\n",
" if path.exists():\n",
" files = [item for item in path.rglob(\"*\") if item.is_file()]\n",
" print(f\"{folder:16} files={len(files):4d} size={format_bytes(sum(item.stat().st_size for item in files))}\")\n",
" else:\n",
" print(f\"{folder:16} missing\")\n"
]
},
{
"cell_type": "markdown",
"id": "773a1932",
"metadata": {},
"source": [
"## Force coefficient history\n",
"\n",
"`postProcessing/forceCoeffs1/0/coefficient.dat` is the simulation-level history for drag, lift, and pitching-moment coefficients. These are dimensionless quantities formed by normalizing forces/moments by freestream dynamic pressure and reference geometry from `controlDict`.\n",
"\n",
"What the common columns mean:\n",
"\n",
"- `Cd`: drag coefficient. Positive drag acts along `dragDir`; smaller positive values usually mean less resistance.\n",
"- `Cl`: lift coefficient. Sign follows `liftDir`; magnitude indicates force normal to drag direction.\n",
"- `CmPitch`: pitching moment coefficient around the configured reference point. Sign indicates nose-up vs nose-down by the case convention.\n",
"- `Time`: for `simpleFoam`, this is an iteration index. It is not physical elapsed time.\n",
"\n",
"The first plot shows the whole convergence history. The second zooms into the final iterations; nearly flat traces imply the steady solve has stopped changing appreciably, while trends/oscillations would warn that final coefficients are less reliable.\n"
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "166e8dad",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"['Time', 'Cd', 'Cs', 'Cl', 'CmRoll', 'CmPitch', 'CmYaw', 'Cd(f)', 'Cd(r)', 'Cs(f)', 'Cs(r)', 'Cl(f)', 'Cl(r)']\n",
"rows: 40000 solver iterations recorded\n",
"final Cd: 0.024150 (dimensionless drag coefficient at last iteration)\n",
"final Cl: 1.675002 (dimensionless lift coefficient at last iteration)\n",
"final CmPitch: -0.520471 (dimensionless pitching-moment coefficient)\n",
"final-500 Cd span: 0.024150 .. 0.024151\n",
"final-500 Cl span: 1.675002 .. 1.675008\n"
]
}
],
"source": [
"coeff_columns, coeff_data = load_force_coefficients(SIM_DIR / \"postProcessing\" / \"forceCoeffs1\" / \"0\" / \"coefficient.dat\")\n",
"coeff_lookup = {name: idx for idx, name in enumerate(coeff_columns)}\n",
"time = coeff_data[:, coeff_lookup[\"Time\"]]\n",
"cd = coeff_data[:, coeff_lookup[\"Cd\"]]\n",
"cl = coeff_data[:, coeff_lookup[\"Cl\"]]\n",
"cm_pitch = coeff_data[:, coeff_lookup[\"CmPitch\"]]\n",
"\n",
"final_window = min(500, len(time))\n",
"print(coeff_columns)\n",
"print(f\"rows: {len(coeff_data)} solver iterations recorded\")\n",
"print(f\"final Cd: {cd[-1]:.6f} (dimensionless drag coefficient at last iteration)\")\n",
"print(f\"final Cl: {cl[-1]:.6f} (dimensionless lift coefficient at last iteration)\")\n",
"print(f\"final CmPitch: {cm_pitch[-1]:.6f} (dimensionless pitching-moment coefficient)\")\n",
"print(f\"final-{final_window} Cd span: {cd[-final_window:].min():.6f} .. {cd[-final_window:].max():.6f}\")\n",
"print(f\"final-{final_window} Cl span: {cl[-final_window:].min():.6f} .. {cl[-final_window:].max():.6f}\")\n"
]
},
{
"cell_type": "code",
"execution_count": 11,
"id": "dbd12749",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAABQgAAAGaCAYAAABDga41AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjExLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvctoD+AAAAAlwSFlzAAAPYQAAD2EBqD+naQAAzTFJREFUeJzs3Xd4U9UbB/DvTdO9W0rpglL2kr2EsgSZAgICypApDvZGcSD8GAKyVEDZCBYBBUSGgOwtQ2TKkFEotHTvNrnn90dtbJq0TdOEpu338zx9mpx77rnvuSdpT97cIQkhBIiIiIiIiIiIiKhEUhR2AERERERERERERFR4mCAkIiIiIiIiIiIqwZggJCIiIiIiIiIiKsGYICQiIiIiIiIiIirBmCAkIiIiIiIiIiIqwZggJCIiIiIiIiIiKsGYICQiIiIiIiIiIirBmCAkIiIiIiIiIiIqwZggJCIiIiIiIiIiKsGYICQqgCtXrqBly5ZwdnaGJEnYt28fAGDDhg2oVKkSlEolJEkCAHTo0AF16tTJV/tfffUVJEnC06dPTR16gcyaNQuSJCEhISHPupbaB/qPvjEqjHGzxNeKJcZkiKIaNxFRYVq1ahUkScL9+/eLxXbMpajGX1TjppKhVatWaNKkSYnZLlkmJgipyLtz5w4++OADVKlSBQ4ODnB3d0ft2rUxatQoXLt2zWzbFUKgb9++cHZ2xuPHjyGEQIcOHXD//n0MGzYM77zzDpKTkyGEMFsMprJgwQJIkoTnz58XdigalhgTEWA5r01LiYOIqCgIDQ2FJEk5/qhUqsIOUcfVq1dzjPerr77SqZ+cnIwpU6YgICAAtra2qFq1KpYuXaq37fzUNUTmF1OhoaFGt2EKlhIHEVFRxAQhFWnbtm3DSy+9hDt37uDbb79FeHg4Hj58iAULFmiO7jOXZ8+e4caNG+jWrRtcXFw05SdOnEB6ejr69u0La2trTfm+fftw+fLlfG1j5MiREEKgTJkypgr7hSsOfaCSq6i+fotq3ERE5jZmzBgIIXR+lEolhg0bBiEEAgMDCztMLWvXrtWJd+TIkTr1evfujQ0bNiAkJASxsbGYOXMmJk+ejM8++6xAdbOz1P2Ul6IaNxHRi6Is7ACIjHX16lX0798fnTp1wrZt26BQ/JfvbteuHdq2bYvPP//cbNsPDw8HANjb2xtUTkRERERkDr/99ht2796NjRs3olmzZgCAN954A2fPnsWcOXMwYsQI+Pj45LsuERGVHDyCkIqsuXPnIi0tDV9++aVWcjCTJEn49NNPtcr++OMPdO7cGW5ubrCzs0OdOnWwatUqnXXv37+PQYMGwcfHBzY2NggKCsKMGTM0p5/0798ftWvXBgAMGDBA65SPCRMmAAC8vLx0TgfJfg3C+/fvY8iQIfD394ednR1q1qyJJUuWID09HUDO1xHLKz7gv+sERkdHY/z48fDy8oKTkxO6d++OsLAwTb2RI0di0qRJOjGfOHEizzFIT0/Pte2c+vDs2TMMHz4cZcuWhb29PapUqYLJkycjOjra4JgMGcvMffD8+XOMGTMG3t7ecHZ2xvbt2yFJEnbv3q3Tp19++QWSJOGnn37Kte95bX/Lli1a16XMas+ePZAkCT///LOmLD9jmr0/Odm2bZvW68/JyQmNGzfG5s2bc+1bfuQ1lpkMfe9lZQn7MPvr15DXZkHjMGTc8oojp78d+Xnf5PW3AzB8/ImIigJ916jLvJzD06dPMWXKFJQuXRqOjo7o0qWLzmmsu3fv1vr77ejoiIYNG2L9+vVmj3379u1QKBTo2rWrVnmPHj2QlpaGnTt3GlVXn+z7aeLEiRg1ahQAICAgQNP/gwcPatYJDQ3F0KFD4evrCxsbG5QvXx4ff/wx0tLSNHUy9/WTJ08wceJElClTRnM2jiH7Nq84croG4eXLl9G1a1e4u7vDzs4OL730EpYvX65VJz+vg+fPn+O9995DuXLlYG9vj8qVK2P8+PGIjIzUqidJElq1apXrvs706NEjzf9bOzs7VK9eHQsWLNDaf6bsx86dOyFJEnbs2KETy++//w5JkvDDDz9oygo6vgCwePFiVKhQAXZ2dmjQoAGOHj2KkSNHwsnJSSeGo0ePon379nB1dYWdnR3q16+Pbdu2adXJvLbe48eP8dprr8HJyQmlS5fGxIkT9V5S4OjRo+jUqRM8PDzg7OyM4OBg7N27N9/bzYkh7QPIV7yGxGLodrP69ddf4ezsjF69eiE5Odmg/lExIYiKKE9PT1G9enWD6589e1bY2tqKTp06ib///ltERESIL774QigUCvHxxx9r6t25c0eUKlVKBAcHiwsXLoj4+Hixf/9+4ePjI/r166epd+nSJQFAbNy4UWs78+fPFwBERESEVnn79u1F7dq1Nc9v3bolPDw8RMOGDcXJkydFfHy8uH79uhg7dqw4evSoEEKIZcuWCQAiLCws3/HNnDlTABDDhw8XmzdvFjExMeLUqVOiTJkyonPnzgbFnJP8tK2vD61atRK1atUSf/75p0hJSRF37twRCxcuFEuWLDEoJkPHMjPOAQMGiO+//15ERUWJb7/9VqSnpws/Pz/RsWNHnbY7dOggvLy8RFpaWo79N2T7qampwtPTU/Tq1Utn/Z49e4rSpUtrtpHfMc3eH0PIsiyePHkiZs6cKSRJEnv27NEs0zdG+sr0MWQsDR2v7Nu0hH2obz/k9to0dRy5jVtuceiLO7/vG0Pe34aMPxGRJXj06JEAIMaMGZNjne+++04AEP/884+mLPNv7dChQ8WGDRtETEyMOHfunPD39xevvPJKjm3JsizCwsLE3LlzhUKhED/99FOu29Hnr7/+EgCEp6ensLGxEc7OzqJ58+Zi69atOnWbNGki/Pz8dMqjoqIEAPH+++8bVVcfffFn/t959OiRTv379+8Lb29v0bRpU3H+/HmRkJAgDh06JPz9/UXPnj019TL39VtvvSXWrl0roqKixMqVK3Xay23f5haHvrgvXbokHBwcRLt27cTNmzfF8+fPxaJFi4SVlZWYNGmSTmyGvA46dOggqlatKi5evCiSk5PF3bt3xZIlS8T8+fO16gEQLVu2zHE/Z7p3757w8vISdevWFcePHxdxcXHixo0bYtKkSeK3334zSz/S09OFj4+P6NKli048/fr1E+7u7iI5OVkIYZrxnTlzprCyshILFy4Uz58/F7du3RKvv/66aN++vXB0dNTa/pYtW4RCoRAffPCB+Oeff0RUVJRYtmyZUCqVYu3atZp6LVu2FLVr1xY9evQQp0+fFnFxcWLdunVCkiSxaNEirTY3bdokFAqFGDx4sLh586aIj48XJ0+e1PqsYOh29TGk/fzEa2gshm63cePGmudLliwRVlZWYsqUKUKW5Vz7RcUPE4RUJCUkJAgAokOHDgav06ZNG+Hl5SUSExO1yt9++21hbW0tnjx5IoQQ4vXXXxelSpUSUVFRWvVCQkIEAHHp0iUhRMEThF26dBFubm4iMjIyx5j1fcg3NL7MD/lffvmlVr158+YJAOLx48d5xpyT/LSdvQ9qtVpYWVmJGTNm5LqN3GIydCwz4/zf//6n08aMGTOEQqEQ9+7d05TdvXtXSJIkJkyYkGtshm5/7NixwsbGRqsPERERwsbGRkycOFFTlt8x1def/GjevLno3r275rmxCUJDx9LQ/aVvm4W9D/ObIDTnWGYft/wmCPP7vsnr/W3o+BMRWYLMBKG+n5kzZwohck8Qzp49W6u9JUuWCADi7t27eW67bdu2WnNWQxOE169fF++99564cOGCSEhIEDdu3BADBgwQAMTcuXO16laqVEnUqFFDpw1ZlgUA0bt3b6Pq6pPfBOGbb74p3NzcdP5f7dixQwAQp0+fFkL8t6+zfmmVl+z7Nr8Jwk6dOgk3NzcRGxurVXfEiBHCyspK3L9/Xys2Q14HdnZ2Ytq0aQb3IS+9evUSTk5O4tmzZznWMUc/pk6dKqysrDTzAyGEiImJEfb29mLkyJGasoKOb2xsrHBwcBADBgzQKXd1ddVKECYnJwsvLy/Rvn17nX0wbNgw4e3
"text/plain": [
"<Figure size 1300x420 with 2 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fig, axes = plt.subplots(1, 2, figsize=(13, 4.2))\n",
"axes[0].plot(time, cd, label=\"Cd: drag coefficient\")\n",
"axes[0].plot(time, cl, label=\"Cl: lift coefficient\")\n",
"axes[0].plot(time, cm_pitch, label=\"CmPitch: pitching moment\")\n",
"axes[0].set_title(\"Coefficient history over all solver iterations\")\n",
"axes[0].set_xlabel(\"simpleFoam iteration\")\n",
"axes[0].set_ylabel(\"dimensionless coefficient\")\n",
"axes[0].grid(alpha=0.25)\n",
"axes[0].legend()\n",
"\n",
"window = min(500, len(time))\n",
"axes[1].plot(time[-window:], cd[-window:], label=\"Cd\")\n",
"axes[1].plot(time[-window:], cl[-window:], label=\"Cl\")\n",
"axes[1].set_title(f\"Final {window} iterations: convergence check\")\n",
"axes[1].set_xlabel(\"simpleFoam iteration\")\n",
"axes[1].set_ylabel(\"dimensionless coefficient\")\n",
"axes[1].grid(alpha=0.25)\n",
"axes[1].legend()\n",
"fig.tight_layout()\n"
]
},
{
"cell_type": "markdown",
"id": "ce49eb43",
"metadata": {},
"source": [
"### Reading the force-coefficient plots\n",
"\n",
"- The full-history panel shows how the solver approached a steady solution from its starting fields.\n",
"- The final-window panel is the practical quality check: if `Cd` and `Cl` are almost horizontal, the final printed numbers are representative of the converged state.\n",
"- The numbers are dimensionless. They let cases with different speeds be compared more directly than raw Newton forces, because the freestream normalization removes much of the speed scaling.\n"
]
},
{
"cell_type": "markdown",
"id": "867f9511",
"metadata": {},
"source": [
"## Mesh geometry and boundary patches\n",
"\n",
"The mesh points are vertices. The volume fields below are cell-centered arrays, so vertex count and field row count differ. Boundary faces tell us where the aerofoil wall and freestream patches live.\n",
"\n",
"The mesh plots are geometry/mesh-density plots, not solution plots. Dense regions indicate where the CFD mesh has more resolution. For aerofoils, high density near the wall is expected because boundary-layer gradients are steep there.\n"
]
},
{
"cell_type": "code",
"execution_count": 12,
"id": "3b7a5b0e",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"mesh points: (520500, 3)\n"
]
},
{
"data": {
"text/plain": [
"{'aerofoil': {'type': 'wall', 'nFaces': 908, 'startFace': 516702},\n",
" 'freestream': {'type': 'patch', 'nFaces': 1624, 'startFace': 517610},\n",
" 'frontAndBack': {'type': 'empty', 'nFaces': 517968, 'startFace': 519234}}"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"points = parse_foam_list(SIM_DIR / \"constant\" / \"polyMesh\" / \"points.gz\", columns=3)\n",
"patches = parse_boundary(SIM_DIR / \"constant\" / \"polyMesh\" / \"boundary\")\n",
"print(f\"mesh points: {points.shape}\")\n",
"patches\n"
]
},
{
"cell_type": "code",
"execution_count": 13,
"id": "e7ce6dff",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAABHEAAAGaCAYAAACIU9+nAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjExLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvctoD+AAAAAlwSFlzAAAPYQAAD2EBqD+naQABAABJREFUeJzsvXmcG3d9//8a3dLoXEl732t7ba/PHE4cwpHDKTkckhByUAhX+g0U+oMCBdLSUs5QoIRSWppAIC0pBCglQJNAnDshCXYSx46P2Gt77V3vodUeum9pfn9oZzwazUgz0ujY9ef5ePixXu3MZz6fz3xmNO/XvA+KYRgGBAKBQCAQCAQCgUAgEAiEpkbT6A4QCAQCgUAgEAgEAoFAIBDKQ0QcAoFAIBAIBAKBQCAQCIRlABFxCAQCgUAgEAgEAoFAIBCWAUTEIRAIBAKBQCAQCAQCgUBYBhARh0AgEAgEAoFAIBAIBAJhGUBEHAKBQCAQCAQCgUAgEAiEZQARcQgEAoFAIBAIBAKBQCAQlgFExCEQCAQCgUAgEAgEAoFAWAYQEYdAIBAIBAKBQCAQCAQCYRlARJyzmMcffxwUReHxxx9vdFeq4uWXXwZFUXjooYca3ZUVTy3memZmBhRF4Tvf+U7ZbX0+H971rnfB4/GAoij84z/+o6JjifVfyZiqPX4p5ufnccstt8DtdnNtP//886AoCv/3f/+nuD0l41JyDggEAoGwclkp3weJRAIUReErX/lKo7tCqBOvv/46LrnkEtjtdsXPqj/84Q9BURROnjxZ8jMCoVkgIk4DOHbsGCiKAkVR+Ju/+RvRbb785S9z27z22mv17WATUo0xS1g5fOpTn8Krr76KvXv3gmEYVUWURh//s5/9LPbs2YPXXnutIWMjEAgEQmMxmUygKApXXnml6N9//etfc8+GDzzwQJ1713wEAgFQFIWvf/3rje4KoQm49dZbodPpMDExAYZhcN111zW6SwRCzdA1ugNnMzRN44EHHsDXv/51aLXagr/953/+J2iaRjQabVDvlg/nnXceGIZpdDcIdeCZZ57BFVdcgZ6enhV3/McffxyXXXZZQdsXX3wxWdsEAoFwFkHTNHbt2oWpqSl0dnYW/O3+++8nz4YyMZlM5PvzLGJ+fh4HDx7EHXfcAYfDoXj/22+/HbfffnsNekYg1AbiidNAbrjhBszMzOD3v/99wefPPvssjh8/jhtvvLFBPSMQmpO5uTmYzeYVefzZ2dmGjo1AIBAIjeeyyy6D1WrFT37yk4LPfT4fHnnkEfJsSCCIMDc3BwDkOYpw1kBEnAayZs0abN++Hffff3/B5z/+8Y+xZcsWbNmyRXS/SCSCz372sxgaGoLRaERbWxv+4i/+AvPz89w26XQaX/jCF7BmzRpYLBb09fXhPe95D44fPy7a5ve//30MDAzAZDJh+/btePnll0v2/eTJk9BoNKIhH7OzszAYDPjsZz+rqM/8HD3f/e53MTg4CK1Wi7vuugtvfvObAQA7d+7kXIm/9a1vAZDO/RGJRHDnnXdizZo1MJlMGBwcxCc/+UkEAoGazOWWLVvgdDpLzpvc9k6fPs2Nk6IoGI1GrFu3Dl/96leRyWRE5+zb3/42ent7YbPZcNNNNyEYDAIA/vmf/xl9fX0wm814+9vfjsnJyYL+8Nv41re+hZ6eHpjNZlx88cV46aWXyo5H7jwC+VDCnTt3wmq1wuv14hOf+ARSqVTZ9r/3ve+BoigkEgn827/9GzcvMzMzuP/++0FRFN54442Cfdg5/N73vidrDJUev5JzxV/ff/u3fwuKohCPxwva5v8ThhHKnW8xKj0HBAKBQKg9ZrMZN910E/7zP/+z4PMHHngAFosFN9xwg+h+2WwW3/72t7Fx40aYTCa4XC68853vxOjoaMF299xzD7Zs2QKbzYb29nZcc801+NOf/iTa5sMPP8y1t379evzmN78p2fdUKgWv1ysqNOVyOfT29uLqq69W1Gd+jp5f//rX2LRpEwwGA77+9a/D5XIBAO68807u+/L9738/AOmcONlsFt/5znewZcsWWCwWdHZ24rbbbsOpU6dqMpfvec97QFEUTp8+XXLu+OOUM+9y++h0Orm50ev16O/vx8c//nGEQqGyc/zggw9K9ldOu3L7We74r732Gq699lq4XC6YTCZs2rQJ3//+97n9b7/9dqxduxYA8Bd/8RegKArt7e3c38vtD5D8N4RlCEOoO6OjowwA5stf/jJz7733MkajkZmfn2cYhmHC4TBD0zTzne98h7n77rsZAMzevXu5faPRKHPOOecw/f39zKOPPsqEQiFm3759zHnnncds2LCBicViDMMwzOc+9znG6XQyjz32GBONRpmpqSnmpz/9KfOxj32Ma2vXrl0MAObWW29lvvKVrzA+n485fvw4c8EFFzCdnZ1MIpEoOY7LL7+c6e3tZbLZbMHn3/zmNxkAzBtvvKGoz2x/brjhBuYf//EfmZmZGebxxx9njh49yjz33HMMAOZ3v/tdUT/27NnDAGB+/etfc59FIhFmy5YtTFdXF/PrX/+aWVhYYE6ePMncfffdzL/927+pPpcMwzCbN29mHA5HyTlT0h6f+fl55qc//SljtVqZL3zhC9zn7JzdeOONzDe+8Q1mbm6OefXVV5nu7m7m5ptvZr75zW8yd911F+P3+5l9+/Yxvb29zNvf/vaCttk2rr/+euZLX/oSMzs7yxw7doy55pprGLPZzOzbt6/kXMudx9nZWaa9vZ0555xzmL179zKLi4vMj370I+aWW25hADB333132bkzGo3MRz/60YLPfvzjHzMAmMOHDxd8PjExwQBg/vVf/7Vk/8U+U3J8IeXOldj6lmpbbN3LnW+xcalxDggEAoFQG4xGI3PzzTczL7zwAgOAeemll7i/bdiwgbn99tu575Kf/OQnBfu+613vYpxOJ/Nf//VfzMLCAnPixAnm+uuvZzweDzMxMcEwDMM88MADjEajYe6//34mGAwy8/PzzKOPPsrceOONXDvT09Pcc8X/+3//jzl58iTj8/mYW265hdHr9cypU6dKjuETn/gEo9frmdnZ2YLPH330UQYA86tf/UpRn9n+vOMd72A+8IEPMGNjY8zBgweZp59+mllcXGQAMHfddVdRP+LxOPeszZLL5Zjrr7+esdlszD333MNMT08zPp+PeeCBB5jPfvazqs8lwzDMn//5nzMAuP2kUDrvcvooJBwOM4899hjT19fHvPOd7yw6ttgcy0GqXbn9LHX8vXv3MhaLhdmxYwfzxhtvMHNzc8zdd9/NaLVa5m/+5m+44xw+fJgBwPzgBz8oOL7c/X/wgx8wAJixsbGSnxEIzQIRcRoAX8QJBoOM2Wxmvve97zEMwzA/+tGPGL1ez/j9flER55/+6Z8YAMyePXsK2jx+/Dij1WqZ73//+wzDMMyFF17IXHnllSX7wT4I3HTTTQWfP/PMM0VftGI8+OCDDADmD3/4Q8Hn69atYy6++GLFfWb7s3PnzqJjKRVxvvKVrzAAmD/+8Y+S/VdzLpVQTXuf/vSnmfb2du53ds5uueWWgu2+8Y1vMBqNhnnve99b8Pk///M/MwCY06dPF7Xxjne8o2DbSCTCuN1u5tprr+U+E5trufN45513Mlqtljl27FjBdp/73OdWlIjDInWuxNa3VNti617ufIuNS41zQCAQCITawIo4DMMww8PDzIc//GGGYRhm9+7d3DONmIjz8MMPMwCYH//4xwXtRaNRpq2tjftuuf3225menp6SfWAN6pGRESaXy3Gfz87OMlqtlvniF79Ycv8DBw4wAJh//ud/Lvj8xhtvZLxeL5NKpRT1me3P0NBQ0UtDpSLOQw89xABg7r33Xsn+qzmXSlAy73L7KMX999/PAGDm5uYKji02x0oQtqvGOb7qqqsYp9PJBIPBgs/vuOMORqvVMidPnmQYRlrEkbs/EXEIyw0STtVg7HY7rr/+ei6k6sc//jF27twJj8cjuv3vfvc79Pf347zzziv4fHBwEAMDA3jmmWcAAJs3b8bjjz+OL37xizh06FDJ5G5
"text/plain": [
"<Figure size 1250x420 with 2 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"rng = np.random.default_rng(20260719)\n",
"sample_count = min(60_000, len(points))\n",
"sample_idx = rng.choice(len(points), size=sample_count, replace=False)\n",
"point_sample = points[sample_idx]\n",
"\n",
"fig, axes = plt.subplots(1, 2, figsize=(12.5, 4.2))\n",
"axes[0].scatter(point_sample[:, 0], point_sample[:, 1], s=0.2, alpha=0.25)\n",
"axes[0].set_title(\"Mesh vertices: sampled full farfield\")\n",
"axes[0].set_xlabel(\"x coordinate\")\n",
"axes[0].set_ylabel(\"y coordinate\")\n",
"axes[0].set_aspect(\"equal\", adjustable=\"box\")\n",
"axes[0].grid(alpha=0.15)\n",
"\n",
"near = point_sample[(point_sample[:, 0] > -0.5) & (point_sample[:, 0] < 1.5) & (np.abs(point_sample[:, 1]) < 0.6)]\n",
"axes[1].scatter(near[:, 0], near[:, 1], s=0.4, alpha=0.35)\n",
"axes[1].set_title(\"Mesh vertices: near aerofoil\")\n",
"axes[1].set_xlabel(\"x coordinate near chord\")\n",
"axes[1].set_ylabel(\"y coordinate\")\n",
"axes[1].set_aspect(\"equal\", adjustable=\"box\")\n",
"axes[1].grid(alpha=0.15)\n",
"fig.tight_layout()\n"
]
},
{
"cell_type": "markdown",
"id": "eca551a0",
"metadata": {},
"source": [
"## Aerofoil surface fields\n",
"\n",
"The `aerofoil` patch is the wall boundary. We map each boundary face to its face center and color those centers by final surface quantities.\n",
"\n",
"Fields plotted here:\n",
"\n",
"- `forceCoeff`: per-face contribution to force coefficient. The plotted norm combines x/y components, so color means contribution magnitude, not drag or lift sign.\n",
"- `wallShearStress`: near-wall viscous shear stress vector. Higher magnitude often marks stronger skin-friction loading or steep near-wall velocity gradients.\n",
"- `yPlus`: dimensionless wall distance of the first cell. It is a mesh/turbulence-model diagnostic: small values mean the first cell is close to the wall in viscous units. Good/bad thresholds depend on the wall treatment, so use it here as a distribution check rather than a universal pass/fail score.\n",
"\n",
"These are final-iteration values from `40000/`, so they describe the solved state, not the initialization.\n"
]
},
{
"cell_type": "code",
"execution_count": 14,
"id": "13b1a5a2",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"aerofoil faces: 908\n",
"aerofoil centers: (908, 2)\n"
]
}
],
"source": [
"aerofoil_patch = patches[\"aerofoil\"]\n",
"aerofoil_faces = parse_faces(\n",
" SIM_DIR / \"constant\" / \"polyMesh\" / \"faces.gz\",\n",
" start_face=int(aerofoil_patch[\"startFace\"]),\n",
" n_faces=int(aerofoil_patch[\"nFaces\"]),\n",
")\n",
"aerofoil_centers = np.array([points[face, :2].mean(axis=0) for face in aerofoil_faces])\n",
"print(f\"aerofoil faces: {len(aerofoil_faces)}\")\n",
"print(f\"aerofoil centers: {aerofoil_centers.shape}\")\n"
]
},
{
"cell_type": "code",
"execution_count": 15,
"id": "7766ae33",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"surface_force: (908, 3) values, one vector per aerofoil wall face\n",
"wall_shear: (908, 3) values, one vector per aerofoil wall face\n",
"y_plus: (908,) values, range 0.0025 .. 0.3785\n",
"Surface distribution summaries:\n"
]
},
{
"data": {
"text/plain": [
"[{'name': '|forceCoeff_xy|',\n",
" 'count': 908,\n",
" 'finite': 908,\n",
" 'min': 4.3631415850623045e-06,\n",
" 'p01': 2.743305896474191e-05,\n",
" 'mean': 0.002002579090723482,\n",
" 'p99': 0.0119255645132304,\n",
" 'max': 0.012571306964867877},\n",
" {'name': '|wallShearStress_xy|',\n",
" 'count': 908,\n",
" 'finite': 908,\n",
" 'min': 0.0013604693437104711,\n",
" 'p01': 0.022830534823254253,\n",
" 'mean': 9.20671527439559,\n",
" 'p99': 31.8291138984418,\n",
" 'max': 31.8967788099112},\n",
" {'name': 'yPlus',\n",
" 'count': 908,\n",
" 'finite': 908,\n",
" 'min': 0.00251439,\n",
" 'p01': 0.010309987999999999,\n",
" 'mean': 0.15439976128854627,\n",
" 'p99': 0.37796499,\n",
" 'max': 0.37849}]"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"surface_force = parse_foam_list(SIM_DIR / \"40000\" / \"forceCoeff.gz\", columns=3)\n",
"wall_shear = parse_foam_list(SIM_DIR / \"40000\" / \"wallShearStress.gz\", columns=3)\n",
"y_plus = parse_foam_list(SIM_DIR / \"40000\" / \"yPlus.gz\", columns=1)\n",
"\n",
"# Norms collapse vector components into one positive magnitude for color maps.\n",
"# Use component plots instead if you need drag/lift sign on each face.\n",
"force_norm = np.linalg.norm(surface_force[:, :2], axis=1)\n",
"shear_norm = np.linalg.norm(wall_shear[:, :2], axis=1)\n",
"\n",
"print(f\"surface_force: {surface_force.shape} values, one vector per aerofoil wall face\")\n",
"print(f\"wall_shear: {wall_shear.shape} values, one vector per aerofoil wall face\")\n",
"print(f\"y_plus: {y_plus.shape} values, range {y_plus.min():.4f} .. {y_plus.max():.4f}\")\n",
"print(\"Surface distribution summaries:\")\n",
"[\n",
" summarize_array(\"|forceCoeff_xy|\", force_norm),\n",
" summarize_array(\"|wallShearStress_xy|\", shear_norm),\n",
" summarize_array(\"yPlus\", y_plus),\n",
"]\n"
]
},
{
"cell_type": "code",
"execution_count": 16,
"id": "070e6676",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAABZoAAAGRCAYAAADsGAYiAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjExLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvctoD+AAAAAlwSFlzAAAPYQAAD2EBqD+naQABAABJREFUeJzs3Xd8E+UfB/DPJemedFE6WWXLRih7FUQUQVBQQLYs/SFDFAQFUVEUEWU4QBAcICAgm0LZAjJFZEOBMlro3k2TPL8/2oSmSdukSWlpP+/X617Qu+fuvndJ822+99xzkhBCgIiIiIiIiIiIiIiomGSlHQARERERERERERERPdlYaCYiIiIiIiIiIiIii7DQTEREREREREREREQWYaGZiIiIiIiIiIiIiCzCQjMRERERERERERERWYSFZiIiIiIiIiIiIiKyCAvNRERERERERERERGQRFpqJiIiIiIiIiIiIyCIsNBMREVGJmDhxIgYOHFjaYeiMHDkSI0eOLO0wLPbvv/9i8ODBCA0NRfPmzQEA48ePx2uvvVbi+y6t17SsvZfKi/LyO5GXRqPBd999h+effx4tWrTA7NmzTV43KSkJzZs3x2+//VbovNJgaWyWnBciIiIiUylKOwAiIiIqG+7fv4+ffvoJR48eRWxsLDw8PFC1alX06dMHnTt3Nnt7V69exZ07d0og0uK5dOlSaYdgsYSEBLRv3x59+vTB/PnzYWtrCwC4fPkyEhMTS3z/pfWalrX3UnlRHn4n8lu6dCkmTpyI5cuXo27duvDy8jJ53ezsbJw6dQoxMTGFzisNlsZmyXkhIiIiMhULzURERISVK1di7NixaNmyJYYNG4aaNWsiKSkJ27ZtQ48ePTB8+HAsXbq0tMOs8Hbt2oXExES89957qFGjhm7+kiVLoFarSzEyorJh48aNaNmyJQYPHmz2uu7u7jhx4gSCgoJKILLSZcl5ISIiIjIVC81EREQV3LZt2zB8+HCMGzcOixYt0lv27LPP4vXXX8fPP/9cStFRXtpevW5ubnrza9WqVRrhEJU59+/fR7Vq1Yq1rkKh0A1HU95Ycl6IiIiITMVCMxERUQU3ZcoU+Pr64ssvvzS6vFGjRqhXr57evH/++QdLlizBhQsXYGNjg1atWmHChAmoXLlykfszZV3tuLELFy7EvHnzcPDgQTz99NP47LPPjG4zOjoaixcvxqlTp5CRkYHatWtj9OjRaNKkiUHbjIwMfPbZZ9i/fz+cnJwwZMgQvPzyy2bHuXz5cnz//fc4cOAA7O3tAQBbt27FrFmz0L17d3z88ce6bY0aNQpCCCxbtszseLW6du2qG+qga9euUCj0/4yrV68eVq1aZXAOv/nmmyKP96OPPsKmTZsAADKZDJUqVULr1q3x5ptvwsPDo8CYjJk+fTr++usv7NmzxyDGxYsX46effsLOnTvh4eFh0X579eqFp556Su88a49l//792LNnj978O3fuYPHixfj777+RlZWFBg0a4K233kKdOnWKPCZz3rOmnO+CmBKjOefs7t27WLp0Kf7++29kZ2ejSZMmmDRpEgICAvTamfo7kZ85sZhybIX93ms0Gvzyyy/YuHEjoqOjUblyZfTu3RuDBw+GTCbD4cOH8dZbb+HGjRu4e/eurmC8cuVKNGjQoMj1gZwxj7t06YLJkyfjlVdeMeEVe6R3796oX7++3vtxzJgxOHnyJBYtWoRWrVrpzkPv3r0xY8YM9O7dGwDQqlUrqFQqAICtrS0CAgLQt29f9O/f36wYjCnqvJiz77179+Lnn3/GtWvX4Onpie7du2PUqFF6v+fh4eFYvXo1bty4AUdHR4SFheGNN96Ag4ODxcdCRERETwBBREREFdaFCxcEADF69GiT19m8ebOwsbERL7zwgti2bZtYu3ataNCggahcubK4evWqrl3Pnj1Fo0aNirVumzZtRGhoqHj++efF119/LXbs2CE++eQTo/GkpKSIGjVqiNatW4s//vhDHDx4UPzwww+iRYsW4uzZs3rbbN26tXj55ZfFokWLREREhBg7dqwAIDZu3Gh2nPv37xcAxK5du3TrDR8+XEiSJPz8/HTz0tPThZ2dnZg2bZpZ8eZ39uxZ8cYbbwgAYs+ePeLEiRO66emnnxbNmjXTa2/O8d68eVO3rb/++kv89NNPon79+qJ+/foiIyND187Ya5rfpk2bBACxadMmvfkajUZUq1ZNdO3a1Sr7DQ4OFv379zfY/4gRI4Snp6fevGPHjgl3d3fRvn17sX79erF7924xePBg4ejoKA4fPlzo8ZjznjX1fBtjaoymnrO//vpLuLu7ixYtWohffvlF7Nu3TyxYsEDUq1fPajGbGoupx1bY733//v2Fra2t+Pjjj8XevXvFp59+Kuzs7ES/fv2ERqMRiYmJ4sSJEyI4OFi0bdtWF1dqaqpJ6wshxMOHDwUAsWDBAl1MxuYZM3ToUL3f+8zMTOHg4CAkSRLTp0/Xzf/xxx8FAHH58mXdvJMnT+rijYiIEB9++KGwt7cXH374YaFxmBJbUefFlH0LIcTUqVOFTCYTb7zxhti+fbv4888/xfjx48WkSZN0baZNmyYUCoWYMmWKCA8PF7/99psICQkRLVu21Hs/EBERUfnFQjMREVEFpi0KfvHFFya1z8zMFD4+PqJVq1a64owQQsTFxQlXV1fRo0cP3bz8xUFz1m3Tpo2Qy+Viy5YtunkqlcpoTHv27BEAxLFjx/TmazQakZmZabDN8PBwvXaNGjUSrVu3NjtOpVIpnJyc9AotgYGBYtCgQQKAOH/+vBBCiF27dgkAIiIiwqx4jZk7d64AIB4+fKg3v0uXLkYLzaYcb0GioqIEALFq1SrdPFMKzdnZ2cLX11c8//zzevO1x71mzRqr7NfUQrNKpRI1a9YUTZo0MXgPde/eXTRo0KDAWIrzni3O+bYkRiEMz1l2draoVq2aaNCggcjOztZrm56ebpWYTY3FnGMr6Pd+48aNAoBYsmSJ3vrLli0TAMTatWt182rXri169uyp187U9S0pNP/yyy8CgDh37pwQQoi9e/cKSZLEwIEDRfPmzXXtXnnlFREUFFTotoQQYs6cOcLOzk73+lkSmxDGz4up+9b+7s6dO9egrfb9pL3wNn/+fL3lN2/eFHZ2dibnGCIiInqyyUq4wzQRERGVYUqlEkDOLdOmOHbsGB48eIARI0ZAkiTdfA8PD/Tr1w/h4eFIT0+3yrqVKlVCz549dT/L5XKj2/Xx8QGQMyzD3bt3dfMlSYKdnZ1e28qVK6Nr165680JDQ3Hx4kWz47SxsUGHDh2we/duAMDly5cRFRWFyZMnw9fXVzc/PDwcjo6OaN26tdnxWsqU4wWA7OxsLFu2DP369UPr1q3RvHlz9O7dG5IkGbQtikKhwJAhQ7Bjxw5ER0fr5v/444/w8PDQDRdg7f0W5O+//8a1a9cwYsQIg/fQSy+9hPPnz+PWrVtG1zX3PWvq+bYkRlPO2fHjxxEZGYnx48cbDF+SfwiD4sZsaizmnn9jv/ebN2+GXC7H8OHD9dYfMmQI7OzsdMN3FMTS9U3RtWtXSJKk93vfqFEjDBo0CKdPn0ZcXByEENi7dy/CwsL01r116xamT5+O7t274+mnn0bz5s2xYsUKZGVl4caNGxbHVhhT9r1mzRpIkoQ333zTYH3t++nXX3+FTCbDqFGj9JYHBwejRYsW2Lp1a4keBxEREZUNHKOZiIioAtOOMRsTE2NSe+3D6IKDgw2WVa1aFSqVCtHR0ahevbrF6wYFBekV9wry1FNPYf78+Zg9ezZWr16NOnXqoHPnzhg2bJjBg70CAwMN1vf09ERCQgLUajXkcrlZcYaFhWHixIm4f/8+wsPDUblyZTRq1Ahdu3ZFeHg4Jk6ciPDwcLRv315XRDYnXkuZcrwA8MILL+DYsWN477330KRJE7i4uECSJLRq1arACweFGTFiBD777DP89NNPeOedd5CYmIg//vgDr7/+ul4x3dr7NSYyMhJAztjJq1evhsi5ow8AkJCQACDnQWnGXm9z37Omnm9LYjTlnGkLtzV
"text/plain": [
"<Figure size 1450x420 with 6 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fig, axes = plt.subplots(1, 3, figsize=(14.5, 4.2))\n",
"for ax, values, title, label in [\n",
" (axes[0], force_norm, \"Surface force-coefficient magnitude\", \"dimensionless |forceCoeff_xy|\"),\n",
" (axes[1], shear_norm, \"Wall-shear magnitude\", \"|wallShearStress_xy|\"),\n",
" (axes[2], y_plus, \"First-cell wall distance\", \"yPlus\"),\n",
"]:\n",
" sc = ax.scatter(aerofoil_centers[:, 0], aerofoil_centers[:, 1], c=values, s=12, cmap=\"viridis\")\n",
" ax.set_title(title)\n",
" ax.set_xlabel(\"x along aerofoil\")\n",
" ax.set_ylabel(\"y\")\n",
" ax.set_aspect(\"equal\", adjustable=\"box\")\n",
" cbar = fig.colorbar(sc, ax=ax, shrink=0.8)\n",
" cbar.set_label(label)\n",
"fig.suptitle(\"Color shows final value on each aerofoil wall face\", y=1.03)\n",
"fig.tight_layout()\n"
]
},
{
"cell_type": "code",
"execution_count": 17,
"id": "ef843b8d",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAABAcAAAMhCAYAAACUltjnAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjExLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvctoD+AAAAAlwSFlzAAAPYQAAD2EBqD+naQABAABJREFUeJzs3Xd4U2UbBvA7o3vvQScUOihlQ0G2yBaRIUNAVByICiJucTNUVBRFEVQUPxegIKACZW9krzK6905305Xz/VEaKV1JOk6a3L/rytX25JyTJydvTnqevO/zSgRBEEBERERERERERksqdgBEREREREREJC4mB4iIiIiIiIiMHJMDREREREREREaOyQEiIiIiIiIiI8fkABEREREREZGRY3KAiIiIiIiIyMgxOUBERERERERk5JgcICIiIiIiIjJyTA4QERERERERGTkmB4iozbt06RKGDh0KW1tbSCQSbN26VeyQarh+/TqGDRsGa2trdXyff/45JBIJkpKStN7f+vXrIZFIEBcX1+i6//zzDyQSCQ4cOKB94K2sruelzXMVU1uJU0xNfZ9qc4z1/ZwgJrHbqtiP31AcYsXWFs7TK1euhEQiQVZWltihEFELYnKASA/8+++/mDRpEnx9fWFhYYFOnTph1qxZOHz4sNihtQnTp0+HXC5HYmIiBEHAhAkTxA6phgcffBAymQzJycl6GR+JqymJorakNd+n+n5OaA4NtRtjaVNERNS8mBwgEtnOnTvRr18/VFZWYseOHcjJycHvv/+OiooKDBo0CGlpaWKHqNeys7Nx5coVjB8/HnZ2dmKHU0teXh7OnDlTK76nn34agiDAy8tLxOj039y5cyEIAvz8/MQOhZqgNd+n+n5OEBvfU/XjsSEiYycXOwAiY/f+++/DysoKv/32G0xNTQEAoaGh+Pnnn9G3b19IpczhNaS6i6OFhYXIkdQtPT0dgP7GR9QaWvN9qu/nBCIiIn3Fqw4ikeXk5MDZ2VmdGLjdwoUL4erqqv57x44dkEgk6puVlRV69+6N77//vsZ21WMDU1JSsGDBAjg7O8PJyQnPP/88VCoVSktL8eyzz8LFxQU2NjaYM2cOSkpKaj1+UlISHn30UXh6esLU1BT+/v5YsmQJysrKGnxOWVlZmDdvXo1hEosWLUJ2dnaN52Zubl5r2+rneOTIkTqfz+LFi+Hu7g4TExPMnTsXQUFBAIDHHnsMEokE7u7uWh0rAEhMTMRjjz0GHx8fmJubIyQkBCtXrqzxPHU5FnPnzkVgYGCN+O683dntV9djDgBHjhxB//79YWFhAR8fH3zyySd1rqfJ61MXTY+prvuvS11jgKvbQ1paGl566SW4urrCysoK48aNq7Mbta7HtPpxkpOT8eyzz8LZ2RnW1taYMGECYmNjNYr//PnzGD9+PBwcHGBubo6wsDB8+eWX6vsXL16MZ555BgDg7e2tPrYREREAND+WEokEQ4YMaTQebd53mj52Y8e3ofepNvFooqHHAhp/PYD6zzf10eQ4aduWmtJuGmtTzfGeOnjwIMLDw2uca5o6bl6T1wZo/Hytzbn/Tncem9zc3DrP29W38+fPq7fV9Dyj6Xm6Ltq2I01iasrxio+PR5cuXdCxY0fcuHEDQPOe/4mo9TE5QCSyAQMGICYmBuvWrYMgCA2uO27cOAiCAEEQoFKpEB0djcmTJ+ORRx7BH3/8UWv9V199Ff3790dMTAw2btyINWvW4IMPPsAzzzyDfv36ISoqCps3b8amTZvwzjvv1Ng2Pj4evXr1QmRkJP78808oFAp888032LBhA2bMmNFgnLNmzcKBAwewdetWKBQK/PPPP/Dz88N3332n/QG6zQsvvIDQ0FBERkbiiy++wPr16xEZGQkA6uNXPQxD02MVGxuLnj174syZM/jpp5+QmZmJ33//HRkZGTh48GCTjkVd8VXfVq9eXWv9phzz06dPY/jw4XB3d8eVK1dw+vRpKBQKrF27tta6ur4+mh7Tlnr97/T6668jNDQUN2/exIEDB3DhwgXMmTOnxjpNOabVXnjhBfTq1QvR0dE4dOgQYmNjMWjQIOTk5DS43fnz53HXXXdBqVTixIkTSE5OxiOPPIJnnnkGL774IoCqf/ar20L1+HhBEDB8+HAArXcs66LJY2tyfBt6nza3hh5Lk9fjdneeb+qjzWukSVtqartprE01RJP31KlTpzBixAh4e3vj6tWrOHPmDAoLC+s812hK09dGk/O1tp+TDbG3t69x3q5uT35+fnB0dISHhwcAzc8z2pynG6JJO9I0Jl2P18mTJ9G3b184OjrixIkT6NSpEwBxz1lE1AwEIhJVTk6OcM899wgABBcXF2HChAnCO++8I5w6dUrjfQwfPlwYNWqU+u8PP/xQACCsWLGixnozZswQrKyshPfee6/G8lmzZgkuLi41lk2fPl2wt7cXMjMzayzfunWrAEA4fvx4vfGYm5sLr7zySoMxL1iwQDAzM6u1fPv27QIA4fDhw7Wez5IlS2qtHxkZKQAQ1q1b1+DjVbvzWE2ePFmwtrYW0tPT692mKceivvhWr14tABASExO1fpx169YJAITY2Fj1OqNHjxbc3NyEkpKSGtuOGjVKACDs379fvUyT10cbdx5TXfdf1/Oqa1l1e1i2bFmN7T/99FMBgBAdHa1e1pTXrvpx3nzzzRrLr169KkgkEuGNN95oMM4xY8YI9vb2Ql5eXo3tn3jiCUEmkwlxcXGCINTdFqo192ulzftOk8fW9PjW9z7QJp66jnFd6nssTV+Phs43ddHkOGnTlpqj3TR0X1PfUyNHjhQ8PDwEpVJZY91x48bVOtfUpSnvFU3O1/W58zyl6fnmdkVFRULv3r0FMzMz4dChQ+rlmr4PtDlP10WbdtSUc58g1P9/RWZmpvDrr78K5ubmwuzZs4XS0tIa2zX3OYuIWhd7DhCJzMHBAbt370ZkZCTefvtteHh4YOPGjejTpw/GjBmDwsJC9bqVlZVYuXIlunfvDisrqxrdRaOiomrte/To0TX+DgoKQlFRUa3lwcHByMzMREFBgXrZ9u3bMWTIEDg7O9dY9+677wYA9bc0denatSu+/vprrFmzBgkJCZofjEaMHz9e43U1PVZ///03hg8fXmP4xp2aciy0oevjCIKA/fv3Y8SIEbW6aNdVpV3X10fTY9pSr/+dxo4dW+Pv0NBQAEBMTIx6WXO8dne2u+DgYAQGBmLfvn31biMIAvbt24fhw4fD1ta2xn2TJ09GZWWlRt2vW+tY6vrYrfXeaCpdXg9NzzfavEaNtaXmaje6auw9JQgCDh48iHvuuQdmZmY11tXm/Hw7bZ6zJudrbT8nNaVSqfDggw/i9OnT+O677zBw4ED1fZq8D7Q9TzdEk3OSpu9NbY/X8uXLMX36dLz66qv4/vvvaw2JFPOcRURNx+QAkZ4ICgrCvHnzsGbNGty4cQPvv/8+/v77b7z55pvqdV566SW8+uqreOaZZxATE4OKigr1NF3l5eW19lnd5bGajY1Ng8vz8vIAAEVFRSgsLMS2bdsgl8shk8kgk8kglUrV6zY0fnDTpk0YOXIkXn75Zfj6+qJ9+/Z4/vnnNRpzKDQwtKJdu3aNbl9Nk2NVVFSEoqKiBvfb1GOhqaY8TlFREZRKJdzc3GrdV9cyXV8fTdtfU15/bdzZjqsvLHJzcwE032tX33FtaL7v4uJiKJXKGuPdq1Uv02S+8NY6lnW97xp77JZ8bzR0HtCFLq+HpucbbV6jxtpSc7UbXWnynlIqlXVenDd0wd4QTZ+zJudrQPvPSU0tWrQIW7duxbvvvovp06erl2v6PtD2PN2QxtqRNu9NbY/XDz/8ADc3Nzz00EN1xtZa5ywiahmcrYBIT73wwgt46623cPjwYfWyH374QT0W8Hb1FUeTSCRaLa9maWkJCwsLTJo0CRs3btQy8qoiWP/73/9QUVGBixcvYufOnVi+fDnOnDmj/gbIzs4OpaWlKCsrq/HNQ3Jycr37bago2J00OVaWlpawtLRs8DGbeiw01ZTHsbKygrm5uXpmhNvVtUyT16cumrY
"text/plain": [
"<Figure size 1050x800 with 3 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"surface_index = np.arange(len(y_plus))\n",
"fig, axes = plt.subplots(3, 1, figsize=(10.5, 8), sharex=True)\n",
"axes[0].plot(surface_index, force_norm)\n",
"axes[0].set_ylabel(\"|forceCoeff|\\ncontribution magnitude\")\n",
"axes[0].grid(alpha=0.25)\n",
"axes[1].plot(surface_index, shear_norm)\n",
"axes[1].set_ylabel(\"|wallShearStress|\\nviscous loading\")\n",
"axes[1].grid(alpha=0.25)\n",
"axes[2].plot(surface_index, y_plus)\n",
"axes[2].set_ylabel(\"yPlus\\nwall-distance diagnostic\")\n",
"axes[2].set_xlabel(\"aerofoil boundary face index (mesh ordering around wall)\")\n",
"axes[2].grid(alpha=0.25)\n",
"fig.suptitle(\"Same surface fields as line plots; useful for spotting localized peaks\", y=0.995)\n",
"fig.tight_layout()\n"
]
},
{
"cell_type": "markdown",
"id": "3dd1357a",
"metadata": {},
"source": [
"### Reading the surface-field plots\n",
"\n",
"- The colored aerofoil plots show where each final wall quantity is located in physical space.\n",
"- The line plots show the same values in boundary-face order. Peaks identify localized regions that may be leading/trailing edges or high-gradient wall zones; the face index is mesh ordering, not a physical distance unit.\n",
"- `forceCoeff` and `wallShearStress` are vector fields. Because these charts use vector norms, they show intensity. They do not distinguish forward/backward or upward/downward direction.\n"
]
},
{
"cell_type": "markdown",
"id": "77f42737",
"metadata": {},
"source": [
"## Final volume fields\n",
"\n",
"`40000/U.gz`, `40000/p.gz`, and turbulence fields are cell-centered OpenFOAM volume fields. These arrays describe the solved flow field over cells; they are not directly indexed by mesh vertices.\n",
"\n",
"Fields summarized below:\n",
"\n",
"- `U`: velocity vector. `|U_xy|` is speed in the 2D plane.\n",
"- `p`: OpenFOAM incompressible pressure, commonly pressure divided by density (`p/rho`), so values are in velocity-squared units rather than Pascals.\n",
"- `nut`: turbulent kinematic viscosity from the turbulence model. Larger values mean the model is adding more eddy viscosity.\n",
"- `k`: turbulent kinetic energy per unit mass. Larger values indicate stronger modeled velocity fluctuations.\n",
"\n",
"The summary table uses min, 1st percentile, mean, 99th percentile, and max so a few extreme cells do not hide the bulk distribution.\n"
]
},
{
"cell_type": "code",
"execution_count": 18,
"id": "99b5127c",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[{'name': 'speed |U_xy|',\n",
" 'count': 258984,\n",
" 'finite': 258984,\n",
" 'min': 9.274197929486949e-05,\n",
" 'p01': 0.21615380686796398,\n",
" 'mean': 32.09927286785441,\n",
" 'p99': 97.0217835728182,\n",
" 'max': 101.66421787590755},\n",
" {'name': 'Ux',\n",
" 'count': 258984,\n",
" 'finite': 258984,\n",
" 'min': -40.1626,\n",
" 'p01': -26.663384999999998,\n",
" 'mean': 24.143595454931,\n",
" 'p99': 57.23976799999999,\n",
" 'max': 64.0673},\n",
" {'name': 'Uy',\n",
" 'count': 258984,\n",
" 'finite': 258984,\n",
" 'min': -4.48142,\n",
" 'p01': -2.5176868000000003,\n",
" 'mean': 12.237414191923309,\n",
" 'p99': 95.408068,\n",
" 'max': 100.078},\n",
" {'name': 'p',\n",
" 'count': 258984,\n",
" 'finite': 258984,\n",
" 'min': -4783.75,\n",
" 'p01': -4551.2704,\n",
" 'mean': -371.1175586042671,\n",
" 'p99': 486.94471999999973,\n",
" 'max': 516.845},\n",
" {'name': 'nut',\n",
" 'count': 258984,\n",
" 'finite': 258984,\n",
" 'min': 6.05286e-19,\n",
" 'p01': 3.6575001e-14,\n",
" 'mean': 0.0006963185103278898,\n",
" 'p99': 0.015897779999999934,\n",
" 'max': 0.0205544},\n",
" {'name': 'k',\n",
" 'count': 258984,\n",
" 'finite': 258984,\n",
" 'min': 1e-15,\n",
" 'p01': 7.4761248e-10,\n",
" 'mean': 1.8431521013304757,\n",
" 'p99': 29.00798499999999,\n",
" 'max': 70.9575}]"
]
},
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"U = parse_foam_list(SIM_DIR / \"40000\" / \"U.gz\", columns=3)\n",
"p = parse_foam_list(SIM_DIR / \"40000\" / \"p.gz\", columns=1)\n",
"nut = parse_foam_list(SIM_DIR / \"40000\" / \"turbulenceProperties:nut.gz\", columns=1)\n",
"k = parse_foam_list(SIM_DIR / \"40000\" / \"turbulenceProperties:k.gz\", columns=1)\n",
"\n",
"speed = np.linalg.norm(U[:, :2], axis=1)\n",
"summary_rows = [\n",
" summarize_array(\"speed |U_xy|\", speed),\n",
" summarize_array(\"Ux\", U[:, 0]),\n",
" summarize_array(\"Uy\", U[:, 1]),\n",
" summarize_array(\"p\", p),\n",
" summarize_array(\"nut\", nut),\n",
" summarize_array(\"k\", k),\n",
"]\n",
"summary_rows\n"
]
},
{
"cell_type": "code",
"execution_count": 19,
"id": "ef2f8457",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAABNgAAANLCAYAAACT4yAfAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjExLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvctoD+AAAAAlwSFlzAAAPYQAAD2EBqD+naQABAABJREFUeJzs3Xd4VEX7//HPpoeSACG00GtoITRBpQtSAijCo/AAj6KAJQIW7BUVy9eCqFFEEWyIKIpSBZSqovSOEJBegpQUCKnz+4NfVpZswia7yW6S9+u6uC52zuw595nZ7J69d86MxRhjBAAAAAAAACBfvNwdAAAAAAAAAFCUkWADAAAAAAAAnECCDQAAAAAAAHACCTYAAAAAAADACSTYAAAAAAAAACeQYAMAAAAAAACcQIINAAAAAAAAcAIJNgAAAAAAAMAJJNgAwINkZmZq/vz5io2NdXcoKKJSU1O1ceNGLViwQBs2bFB6errmz5+vffv2Fehx3fHatXdM/oacd+rUKc2fP1+nT592qP7p06e1YsUKzZ8/XydPntSBAwc0f/58Xbx4Mc/HPnz4sObPn6/k5GSXx+kqeYkR7rN69Wr9+eef7g6j0JXU8wYAT0CCDQAK0KJFi7Ru3Tq72+Lj47MlPlJTU9WvXz/NmjUrT8cprCQKCk9++vTYsWOKjIzUf/7zH33wwQdaunSpkpKS1K9fP33zzTcFGG3+X7uuPqY74sgrT/973bBhg/r166dt27Zdte7s2bNVs2ZNPfPMM5oyZYr+/vtvzZ8/X/369dM///yT52P/9NNP6tevn06ePOnSOF0pLzGiYP3666/6/fff7W4bM2aMnn322UKOyP1K6nkDgCcgwQYABWjAgAGaMGGC3W179+7Nlvjw9vZWVFSUGjRokKfjFFYSBYUnP336xhtv6PTp09q1a5fmz5+vxx9/XL6+voqKilL9+vULMFrkRXH5ezXGaMyYMbrzzju1evVqzZ8/X+3bt1edOnUUFRWlwMBAd4eIYu7hhx/WE0884e4wAACQJPm4OwAAwL98fX01f/58d4eBImrr1q1q0KCB/Pz8rGWlS5fmNYUCcfLkScXFxalp06Y25VFRUYqKinJTVAAAAO5Bgg0APEhmZqYWLlyo8PDwbCOODh8+rP3796t06dJq3LixSpcuLUlKSEjQTz/9JEnavXu3NZlSo0YNtWjRwvp8Y4x27NihY8eOKTQ0VBEREfL29s4WQ2pqqtavX6+MjAy1bt1apUqV0pIlS1SlShVFRERIklJSUrR06VI1bdpUderU0cGDB/XXX3+pWbNmqlatmn766SelpaVJupQ0DAsLU9OmTWWxWKzHuXIf+/bt06FDh9SsWTOFhoZa6+3bt09///23wsPDVb16dbvtllPbOOLIkSPas2ePypUrp6ZNm8rf399m+9XaLT4+XqtXr1abNm1UpUoVm+f+9NNPqlatmpo3b273nA8cOKDY2FjVqVNH9erVsz7P0T7Ncvr0af3+++86ePCgAgIC7CbUGjdubD2Go3Fcfh5X68/8uFrbS9LBgwe1d+9eBQYGqlWrVgU2KsrZWAqib/NyvCv/BiUpIyNDW7Zs0T///GPtM3tOnjyp7du3W1/fjti6dat+++03SdLOnTvtvua6d++ugIAAmzJHY3JVnH///bd27Nihzp07q2zZsjbbkpOT9fPPP6tRo0Zq0KCB4uLirHNXWSwWlSpVSuHh4apatepVj7Nz504dPHhQvXv3tinP2uf111+v8uXL22xzpC0yMjK0c+dOxcXFqUaNGqpfv768vHK/ASWv55GRkaEdO3YoLi5OdevWVd26da3b9u/fr507d6pnz57y9vbWxo0bFRcXpz59+ljrHD58WH/99ZcCAgIUGRmpMmXKZDvGuXPntHv3bqWmpio8PFyVKlXKV53LLVu2TOfOnbP5YSogIEDdu3fPVvfs2bPavHmzgoKC1LJlyxzb8Ny5c9q6davS0tLUrFkzVa5cOdcYpEtznvn7++uaa67R6dOntXXrVlWvXt1mJHpWeUhISK6vXUeOn5d2cvS8AQAuYgAABcbf399ERUXZ3bZu3TojybzyyivWsuTkZCPJvPjii9ay8+fPm759+5pSpUqZjh07muuuu85UqVLFTJw40RhjzP79+03Pnj2NJBMeHm6ioqJMVFSUmTRpknUfGzZsMOHh4aZcuXKmY8eOplKlSqZmzZpm2bJlNjGtWLHCVK1a1YSGhppOnTqZRo0amdWrV5uQkBBz1113WesdP37cSDJvvvmmGTVqlGnWrJlp0KCB+eqrr4wxxtx2223WODp37mzKly9v6tevbzZv3mx3H3fddZdp2bKladasmSlVqpT57rvvTFpamhkxYoSJjIw0ERERxsfHx3zwwQc28V6tbXJz4MABc8MNNxh/f3/TunVr07ZtW1O7dm3z3Xff5andsvrxm2++yXaM4OBgc/fdd2c757feesvcf//9JiIiwrRt29ZYLBZz3333Wes50qeX27p1q4mKijLBwcGmQoUK1vpRUVHW/Vz+OnM0jiyO9Ke9164zbX/w4EHTtWtXU6pUKXPdddeZJk2amODgYDNjxoxcj5mXOFwZiyv7Ni/Hy+lv8NtvvzXVqlUz1apVM126dDGhoaGmZcuWZu/evTbn/+yzzxofHx/TuHFj07ZtW9OzZ08ze/ZsI8ksX748x3Z75513TJcuXYwkExERYfOaa9q0qZFkDh8+bPMcR2L66KOPjCTz999/uyTOjRs3Gklm8uTJ2bbNmDHDSDKrV682xhjz559/Ws+hT58+pnXr1sbb29sMGjTIpKam5hrjuHHjjL+/f7ZjzJs3z+YYeWmLX3/91dSuXdtUr17d3HDDDaZx48amadOm2fZ1JUfPwxhjvvnmGxMWFmYqVqxoOnXqZOrVq2e6dOli7bvXX3/dSDJ//PGHadOmjbn22mtN+fLljTHGnDlzxtx0003Gz8/PtGvXzjRs2NCUKlXKvPrqqzbHmDBhggkMDDQtW7Y0Xbt2NdWqVTO33nqrOXfuXJ7qXGnYsGGmXLlyNu95w4YNs25v0aKF6dmzp5k1a5Zp1KiR6dixoyldurRp06aNOXPmjM2+UlJSzP3332/8/PxMRESEad++vQkICDD33XefSUtLy7W9s44zffp0Ex4ebq677jrj5+dnHnjgAWOMMVOnTjXh4eHm+uuvN35+fmbQoEEmMzMzX8d3pJ3yct4AANciwQYABcjf39+0bdvWzJs3L9u/SZMmOZRge+mll0ypUqVsvqxeuHDBTJkyxfr47Nmz2faV5dSpU6ZixYrm2muvNWfPnrUep3///iYwMND89ddfxhhjTp48acqVK2d69uxpLly4YIwxJiEhwQwePNiUKVPGboItPDzc+oU+LS3N7Nmzx247nD9/3vTu3ds0atTI+sUiax+NGzc2P/74o7XunXfeacqVK2cee+wxmyTH6NGjTenSpc3p06fz1Db2JCYmmjp16pimTZua/fv3W8tPnDhhvv766zy1W34SbE2bNjXffvuttfydd94xkszKlSutZbn1aU5at25tOnfubFNmbz95icMee/3paGLLkbY/f/68qV+/vmnRooU5cuSItc67775rvLy8zKpVq3I8Zl4SbK6MxVV9m9fj2fsbXLp0qfHy8jIPP/ywSU9PN8YYk5SUZG644QbTsGFDa5Jl+vTpRpJ5//33rcdZs2aNiYyMvGriyhhjdu3aZSSZjz76yKb83XffzZZgczQme8krZ+Ns1aqVadGiRbbyTp06mYYNG+b63K1bt5py5cqZF154wVrmbILN0bYIDw83ffv2tUnG/PXXX2b+/Pm5xuzoeSxatMhYLBYTHR1tk3hbsWKFNXmelWDr06ePOXjwoDHGmJ07dxpjjOnVq5cpX7682bRpk/W5b775ps1rYv369UaSmTlzpk0833zzjfX14UidnLRr1y7be16WFi1amHr16plx48ZZ23nHjh3G19fXPPHEEzZ1R4wYYUqVKmXzd7pp0yZTpkwZ8/TTT+caQ4sWLUzdunXNQw89ZO2rrOTvxIkTTXR0tMnIyDD
"text/plain": [
"<Figure size 1250x830 with 4 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fig, axes = plt.subplots(2, 2, figsize=(12.5, 8.3))\n",
"for ax, values, title, xlabel in [\n",
" (axes[0, 0], speed, \"Speed distribution\", \"|U_xy| [m/s-like velocity magnitude]\"),\n",
" (axes[0, 1], p, \"Pressure distribution\", \"p/rho [velocity^2 units]\"),\n",
" (axes[1, 0], nut, \"Turbulent viscosity distribution\", \"nut [m^2/s]\"),\n",
" (axes[1, 1], k, \"Turbulent kinetic energy distribution\", \"k [m^2/s^2]\"),\n",
"]:\n",
" ax.hist(values[np.isfinite(values)], bins=80, edgecolor=\"white\")\n",
" ax.set_title(title)\n",
" ax.set_xlabel(xlabel)\n",
" ax.set_ylabel(\"number of cells\")\n",
" ax.set_yscale(\"log\")\n",
" ax.grid(alpha=0.2)\n",
"fig.suptitle(\"Histograms count final cell-centered field values across the mesh\", y=1.01)\n",
"fig.tight_layout()\n"
]
},
{
"cell_type": "markdown",
"id": "3ae9f126",
"metadata": {},
"source": [
"### Reading the volume-field histograms\n",
"\n",
"- Each bar counts cells whose final value falls in that range. The y-axis is logarithmic so both common values and rare extremes remain visible.\n",
"- These histograms do not show where cells are located. They show distribution only. To connect values to geometry, you would need cell centers and a spatial plot.\n",
"- Wide tails can indicate boundary layers, wake regions, or localized numerical/physical extremes. Use the percentile summary above before focusing on the absolute min/max.\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Raw text previews\n",
"\n",
"Use these to connect the arrays and plots back to the files on disk. Previews are capped so huge logs or fields do not blow up the notebook.\n"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"--- system/controlDict ---\n",
"/*--------------------------------*- C++ -*----------------------------------*\\\n",
"| ========= | |\n",
"| \\\\ / F ield | OpenFOAM: The Open Source CFD Toolbox |\n",
"| \\\\ / O peration | Version: v2112 |\n",
"| \\\\ / A nd | Website: www.openfoam.com |\n",
"| \\\\/ M anipulation | |\n",
"\\*---------------------------------------------------------------------------*/\n",
"FoamFile\n",
"{\n",
" version 2.0;\n",
" format ascii;\n",
" class dictionary;\n",
" object controlDict;\n",
"}\n",
"// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //\n",
"Uinf\t32.137;\n",
"\n",
"application simpleFoam;\n",
"\n",
"startFrom startTime;\n",
"\n",
"startTime 0;\n",
"\n",
"stopAt endTime;\n",
"\n",
"endTime\t40000;\n",
"\n",
"deltaT 1;\n",
"\n",
"writeControl timeStep;\n",
"\n",
"writeInterval $endTime;\n",
"\n",
"purgeWrite 0;\n",
"\n",
"writeFormat ascii;\n",
"\n",
"writePrecision 6;\n",
"\n",
"writeCompression on;\n",
"\n",
"timeFormat general;\n",
"\n",
"timePrecision 6;\n",
"\n",
"runTimeModifiable true;\n",
"\n",
"functions\n",
"{\n",
"\tforces_object\n",
"\t{\n",
"\t type forces;\n",
"\t libs (\"libforces.so\");\n",
"\n",
"\t enabled true;\n",
"\n",
"\t writeControl timeStep;\n",
"\t writeInterval $endTime;\n",
"\n",
"\t patches (\"aerofoil\");\n",
"\n",
"\t p\t\tp;\n",
"\t U\t\tU;\n",
"\t rho\trhoInf;\n",
"\n",
"\t //// Density only for incompressible flows\n",
"\t rhoInf 1.204;\n",
"\t \n",
"\t //// Centre of rotation\n",
"\t CofR (0 0 0);\n",
"\t}\n",
"\t\n",
"\tforceCoeffs1\n",
"\t{\n",
"\t // Mandatory entries\n",
"\t type forceCoeffs;\n",
"\t libs (\"libforces.so\");\n",
"\t patches (\"aerofoil\");\n",
"\n",
"\n",
"\t // Optional entries\n",
"\n",
"\t // Field names\n",
"\t p\t\tp;\n",
"\t U\t\tU;\n",
"\t rho\trhoInf;\n",
"\t \n",
"\t ////Density only for incompressible flows\n",
"\t rhoInf 1.204;\n",
"\n",
"\t // Reference pressure [Pa]\n",
"\t pRef 0;\n",
"\n",
"\t // Include porosity effects?\n",
"\t porosity no;\n",
"\n",
"\t // Store and write volume field representations of forces and moments\n",
"\t writeFields yes;\n",
"\t writeControl timeStep;\n",
"\t writeInterval $endTime;\n",
"\n",
"\t // Centre of rotation for moment calculations\n",
"\t CofR (0 0 0);\n",
"\n",
"\t // Lift direction\n",
"\t liftDir\t (-0.20999398925280716 0.9777026769308202 0);\n",
"\n",
"\t // Drag direction\n",
"\t dragDir\t (0.9777026769308202 0.20999398925280716 0);\n",
"\n",
"\t // Pitch axis\n",
"\t pitchAxis (0 0 1);\n",
"\n",
"\t // Freestream velocity magnitude [m/s]\n",
"\t magUInf $Uinf;\n",
"\n",
"\t // Reference length [m]\n",
"\t lRef 1;\n",
"\n",
"\t // Reference area [m2]\n",
"\t Aref 1;\n",
"\n",
"\t // Spatial data binning\n",
"\t // - extents given by the bounds of the input geometry\n",
"\t /*binData\n",
"\t {\n",
"\t\tnBin 20;\n",
"\t\tdirection (1 0 0);\n",
"\t\tcumulative yes;\n",
"\t }*/\n",
"\t}\n",
"\n",
" momErr\n",
" {\n",
" type momentumError;\n",
" libs (fieldFunctionObjects);\n",
" executeControl writeTime;\n",
" writeControl writeTime;\n",
" }\n",
"\n",
" contErr\n",
" {\n",
" type div;\n",
" libs (fieldFunctionObjects);\n",
" field phi;\n",
" executeControl writeTime;\n",
" writeControl writeTime;\n",
" }\n",
"\n",
"\n",
" turbulenceFields1\n",
" {\n",
" type turbulenceFields;\n",
" libs (fieldFunctionObjects);\n",
" fields\n",
" (\n",
" R\n",
" I\n",
" L\n",
" k\n",
" epsilon\n",
" omega\n",
" nut\n",
" nuEff\n",
" devReff\n",
" );\n",
"\n",
" executeControl writeTime;\n",
" writeControl writeTime;\n",
" }\n",
"\n",
" yplus\n",
" {\n",
"\ttype\t\tyPlus;\n",
"\tlibs\t\t(fieldFunctionObjects);\n",
"\n",
"\tenabled\ttrue;\n",
"\texecuteControl\twriteTime;\n",
"\twriteControl\twriteTime;\n",
" }\n",
" \n",
" wallshearstress\n",
" {\n",
" \ttype\t\twallShearStress;\n",
" \tlibs\t\t(fieldFunctionObjects);\n",
" \t\n",
" \texecuteControl\twriteTime;\n",
" \twriteControl\twriteTime;\n",
" }\n",
" \n",
" mach\n",
" {\n",
" \ttype\t\tMachNo;\n",
" \tlibs\t\t(fieldFunctionObjects);\n",
" \t\n",
" \texecuteControl\twriteTime;\n",
" \twriteControl\twriteTime;\n",
" }\n",
"}\n",
"\n",
"\n",
"// *************************************************************************\n"
]
}
],
"source": [
"print(\"--- system/controlDict ---\")\n",
"print(read_text(SIM_DIR / \"system\" / \"controlDict\", max_bytes=4_000))\n"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"--- 40000/U.gz header and first values ---\n",
"/*--------------------------------*- C++ -*----------------------------------*\\\n",
"| ========= | |\n",
"| \\\\ / F ield | OpenFOAM: The Open Source CFD Toolbox |\n",
"| \\\\ / O peration | Version: 2112 |\n",
"| \\\\ / A nd | Website: www.openfoam.com |\n",
"| \\\\/ M anipulation | |\n",
"\\*---------------------------------------------------------------------------*/\n",
"FoamFile\n",
"{\n",
" version 2.0;\n",
" format ascii;\n",
" arch \"LSB;label=32;scalar=64\";\n",
" class volVectorField;\n",
" location \"40000\";\n",
" object U;\n",
"}\n",
"// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //\n",
"\n",
"dimensions [0 1 -1 0 0 0 0];\n",
"\n",
"internalField nonuniform List<vector> \n",
"258984\n",
"(\n",
"(30.9771 6.60589 7.13234e-22)\n",
"(30.9449 6.59872 0)\n",
"(30.9331 6.59577 -4.37143e-32)\n",
"(30.924 6.59278 2.99688e-21)\n",
"(30.9226 6.5912 -2.39008e-21)\n",
"(30.9275 6.59069 -2.40176e-32)\n",
"(30.9387 6.59126 2.71343e-22)\n",
"(30.9558 6.59283 -4.89239e-22)\n",
"(30.9785 6.59534 -1.07139e-21)\n",
"(31.0064 6.59871 -2.02935e-21)\n",
"(31.0395 6.60292 2.80758e-21)\n",
"(31.0775 6.60792 -3.88218e-21)\n",
"(31.1205 6.61368 -2.44941e-21)\n",
"(31.1682 6.62018 3.2547e-21)\n",
"(31.2205 6.62735 -1.00218e-31)\n",
"(31.2768 6.63505 6.6084e-21)\n",
"(31.3352 6.64288 1.06892e-20)\n",
"(31.3911 6.64978 2.97873e-31)\n",
"(31.4321 6.65311 -1.39049e-29)\n",
"(31.4453 6.65005 0)\n",
"(31.4446 6.64357 4.23594e-28)\n",
"(31.446 6.63704 6.78282e-28)\n",
"(31.4474 6.63002 8.20923e-19)\n",
"(31.4489 6.62246 9.06432e-19)\n",
"(31.4506 6.61435 1.24796e-22)\n",
"(31.4523 6.60564 -2.22734e-18)\n",
"(31.4542 6.5963 6.06896e-19)\n",
"(31.4563 6.58629 6.63523e-19)\n",
"(31.4585 6.57555 0)\n",
"(31.4608 6.56405 -1.32388e-22)\n",
"(31.4633 6.55171 0)\n",
"(31.4659 6.5385 5.05547e-19)\n",
"(31.4687 6.52436 0)\n",
"(31.4717 6.50922 0)\n",
"(31.4749 6.49301 -7.48502e-19)\n",
"(31.4782 6.47567 -1.71535e-27)\n",
"(31.4818 6.45714 -9.07347e-23)\n",
"(31.4855 6.43732 0)\n",
"(31.4895 6.41614 -1.20768e-22)\n",
"(31.4937 6.39353 8.10999e-19)\n",
"(31.498 6.36938 -9.57919e-19)\n",
"(31.5026 6.34362 -4.55391e-27)\n",
"(31.5074 6.31615 -1.34617e-18)\n",
"(31.5124 6.28687 -7.9968e-19)\n",
"(31.5176 6.25568 9.51028e-19)\n",
"(31.523 6.22247 4.54929e-27)\n",
"(31.5286 6.18715 -6.73109e-19)\n",
"(31.5343 6.14959 -8.00407e-19)\n",
"(31.5402 6.1097 -7.64645e-27)\n",
"(31.5463 6.06735 -1.13053e-18)\n",
"(31.5524 6.02244 -1.07685e-26)\n",
"(31.5587 5.97485 1.11139e-26)\n",
"(31.565 5.92447 9.41686e-19)\n",
"(31.5712 5.87119 -1.55796e-26)\n",
"(31.5775 5.81489 2.10475e-26)\n",
"(31.5836 5.75548 -7.75313e-19)\n",
"(31.5896 5.69284 -9.12881e-19)\n",
"(31.5953 5.62689 -1.07297e-18)\n",
"(31.6007 5.55753 0)\n",
"(31.6057 5.48468 1.47569e-18)\n",
"(31.6103 5.40826 8.63075e-19)\n",
"(31.6142 5.32822 1.00778e-18)\n",
"(31.6174 5.2445 0)\n",
"(31.6198 5.15706 -1.36747e-18)\n",
"(31.6213 5.06586 -7.94401e-19)\n",
"(31.6216 4.9709 -5.20441e-26)\n",
"(31.6208 4.87218 0)\n",
"(31.6185 4.7697 -1.23312e-18)\n",
"(31.6146 4.66349 0)\n",
"(31.6091 4.5536 0)\n",
"(31.6018 4.44011 -1.88511e-18)\n",
"(31.5928 4.3231 -1.08226e-18)\n",
"(31.5821 4.20267 -1.24049e-18)\n",
"(31.5694 4.0789 -8.31105e-27)\n",
"(31.554 3.95181 9.2397e-26)\n",
"(31.5342 3.8213 0)\n",
"(31.5075 3.68715 -1.05318e-18)\n",
"(31.4716 3.54914 -6.7894e-23)\n",
"(31.4283 3.40744 2.7224e-18)\n",
"(31.3918 3.26374 -1.76631e-25)\n",
"(31.3952 3.1223 0)\n",
"(31.4677 2.98758 0)\n",
"(31.5514 2.85369 -2.75739e-25)\n",
"(31.4104 2.69109 2.94594e-26)\n",
"(30.7887 2.4651 -1.11872e-20)\n",
"(29.8417 2.2057 5.89645e-21)\n",
"(28.7783 1.94455 0)\n",
"(27.654 1.68842 -8.75387e-28)\n",
"(26.4758 1.43935 -8.61442e-26)\n",
"(25.2452 1.19962 -7.57069e-28)\n",
"(23.9674 0.971889 1.52416e-21)\n",
"(22.6519 0.758733 7.17327e-28)\n",
"(21.3102 0.562214 1.57519e-21)\n",
"(19.9537 0.383784 0)\n",
"(18.5935 0.224254 0)\n",
"(17.2397 0.0838885 0)\n",
"(15.9016 -0.0376447 -9.34507e-22)\n",
"(14.5875 -0.141296 0)\n",
"(13.3044 -0.228505 -5.15522e-22)\n",
"(12.0583 -0.301162 5.3938e-22)\n",
"(10.8529 -0.361165 5.28847e-28)\n",
"(9.69116 -0.410397 -5.71385e-28)\n",
"(8.57561 -0.450505 6.06921e-22)\n",
"(7.50773 -0.482642 3.10611e-22)\n",
"(6.48864 -0.50726 -3.15733e-22)\n",
"(5.52035 -0.524054 6.37224e-28)\n",
"(4.60703 -0.532391 -5.47e-28)\n",
"(3.75636 -0.531905 2.98549e-22)\n",
"(2.9803 -\n"
]
}
],
"source": [
"print(\"--- 40000/U.gz header and first values ---\")\n",
"print(read_text(SIM_DIR / \"40000\" / \"U.gz\", max_bytes=4_000))\n"
]
}
],
"metadata": {
"kernelspec": {
"display_name": ".venv",
"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.13.12"
}
},
"nbformat": 4,
"nbformat_minor": 5
}
Reference in a new issue Copy permalink