{ "cells": [ { "cell_type": "markdown", "id": "232f64cf", "metadata": {}, "source": [ "# AirfRANS in OpenFOAM: algorithm-first, array-visible\n", "\n", "Purpose: understand what OpenFOAM is doing in the AirfRANS steady RANS solve, without drowning in C++ framework structure.\n", "\n", "Pattern used throughout:\n", "\n", "```text\n", "OpenFOAM phase → algorithmic meaning → equations used → arrays touched → parallelization shape → small numeric probe\n", "```\n", "\n", "This is not a generic CFD derivation. It follows the repo's real OpenFOAM v14 `incompressibleFluid` path used to run a migrated AirfRANS `kOmegaSST` case.\n" ] }, { "cell_type": "markdown", "id": "7a6f30ae", "metadata": {}, "source": [ "## 0. Mental model\n", "\n", "OpenFOAM is doing a finite-volume SIMPLE solve.\n", "\n", "```text\n", "cells store unknowns\n", "faces move flux\n", "matrices encode neighbour coupling\n", "pressure repairs continuity\n", "turbulence updates effective viscosity\n", "```\n", "\n", "Core arrays:\n", "\n", "| array | meaning |\n", "|---|---|\n", "| `V[c]` | cell volume |\n", "| `C[c,3]` | cell centre |\n", "| `Sf[f,3]` | oriented face area vector |\n", "| `owner[f]`, `neighbour[f]` | face-to-cell topology |\n", "| `phi[f]` | active solve face volume flux; empty front/back faces are bookkeeping for 2D |\n", "| `U[c,3]` | mean velocity |\n", "| `p[c]` | kinematic pressure |\n", "| `k[c]`, `omega[c]`, `nut[c]` | SST turbulence state |\n", "| `diag`, `upper`, `lower`, `source` | sparse finite-volume matrix storage |\n" ] }, { "cell_type": "code", "execution_count": 1, "id": "e6138eb2", "metadata": { "execution": { "iopub.execute_input": "2026-07-24T06:27:07.329068Z", "iopub.status.busy": "2026-07-24T06:27:07.328767Z", "iopub.status.idle": "2026-07-24T06:27:09.064496Z", "shell.execute_reply": "2026-07-24T06:27:09.063766Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "{'simulation': 'airFoil2D_SST_93.213_3.79_0.418_0.0_9.665', 'cells': 283400, 'internal_face_flux_entries': 565419, 'sparse_offdiag_entries': 565419, 'boundary_patches': ['aerofoil', 'freestream', 'frontAndBack'], 'U_shape': (283400, 3), 'p_shape': (283400,), 'phi_internal_shape': (565419,)}\n" ] } ], "source": [ "# Fresh-kernel setup. Keep framework setup quarantined here.\n", "from pathlib import Path\n", "from contextlib import contextmanager\n", "import os, sys, shutil, json, re\n", "\n", "# Automation can inject an incompatible PYTHONPATH. Drop it before compiled imports.\n", "pythonpath = os.environ.pop(\"PYTHONPATH\", \"\")\n", "for entry in pythonpath.split(os.pathsep):\n", " if entry:\n", " while entry in sys.path:\n", " sys.path.remove(entry)\n", "\n", "import numpy as np\n", "\n", "\n", "def find_repo_root(start):\n", " start = Path(start).resolve()\n", " for p in [start, *start.parents]:\n", " if (p / \"pyproject.toml\").exists() and (p / \"OpenFOAM-14\").exists():\n", " return p\n", " raise RuntimeError(\"repo root not found\")\n", "\n", "\n", "@contextmanager\n", "def quiet_native_output():\n", " # Suppress OpenFOAM banner/solver logs; this notebook prints selected facts instead.\n", " sys.stdout.flush(); sys.stderr.flush()\n", " devnull = os.open(os.devnull, os.O_WRONLY)\n", " saved_out = os.dup(1)\n", " saved_err = os.dup(2)\n", " try:\n", " os.dup2(devnull, 1)\n", " os.dup2(devnull, 2)\n", " yield\n", " finally:\n", " os.dup2(saved_out, 1)\n", " os.dup2(saved_err, 2)\n", " os.close(saved_out)\n", " os.close(saved_err)\n", " os.close(devnull)\n", "\n", "\n", "def q(call, *args, **kwargs):\n", " with quiet_native_output():\n", " return call(*args, **kwargs)\n", "\n", "\n", "ROOT = find_repo_root(Path.cwd())\n", "os.chdir(ROOT)\n", "sys.path.insert(0, str(ROOT / \"scripts\"))\n", "\n", "from openfoam_env import apply_openfoam_env\n", "from prepare_airfrans_stepper_case import DEFAULT_SOURCE, prepare_case\n", "\n", "apply_openfoam_env()\n", "import foam_stepper as foam\n", "\n", "WORK = ROOT / \"tmp\" / \"airfrans_openfoam_algorithm_first_notebook\"\n", "CASE = WORK / \"airfrans_v14\"\n", "SOURCE = DEFAULT_SOURCE\n", "\n", "# Rebuild the migrated case for every full notebook execution.\n", "if CASE.exists():\n", " shutil.rmtree(CASE)\n", "WORK.mkdir(parents=True, exist_ok=True)\n", "meta = prepare_case(SOURCE, CASE, end_time=1)\n", "\n", "with quiet_native_output():\n", " stepper = foam.Case(CASE).make_stepper()\n", " mesh = stepper.mesh()\n", " fields = stepper.fields()\n", "\n", "# Plain arrays used by the rest of the notebook.\n", "V = np.asarray(mesh.V)\n", "C = np.asarray(mesh.C)\n", "Cf = np.asarray(mesh.Cf)\n", "Sf = np.asarray(mesh.Sf)\n", "magSf = np.asarray(mesh.magSf)\n", "owner = np.asarray(mesh.owner, dtype=np.int64)\n", "neighbour = np.asarray(mesh.neighbour, dtype=np.int64)\n", "patches = list(mesh.boundary)\n", "\n", "U0 = np.asarray(fields[\"U\"].internal)\n", "p0 = np.asarray(fields[\"p\"].internal)\n", "phi0 = np.asarray(fields[\"phi\"].internal)\n", "k0 = np.asarray(fields[\"k\"].internal)\n", "omega0 = np.asarray(fields[\"omega\"].internal)\n", "nut0 = np.asarray(fields[\"nut\"].internal)\n", "phi_field0 = fields[\"phi\"]\n", "\n", "print({\n", " \"simulation\": meta.simulation,\n", " \"cells\": int(V.size),\n", " \"internal_face_flux_entries\": int(phi0.size),\n", " \"sparse_offdiag_entries\": int(neighbour.size),\n", " \"boundary_patches\": [patch.name for patch in patches],\n", " \"U_shape\": U0.shape,\n", " \"p_shape\": p0.shape,\n", " \"phi_internal_shape\": phi0.shape,\n", "})\n" ] }, { "cell_type": "markdown", "id": "763c6a18", "metadata": {}, "source": [ "## 1. Case contract: what algorithm OpenFOAM has been asked to run\n", "\n", "Before thinking about kernels, read the dictionaries as the numerical contract:\n", "\n", "- `controlDict`: steady run controls and selected solver;\n", "- `fvSchemes`: discrete operators, interpolation, gradients, divergence, laplacians;\n", "- `fvSolution`: linear solvers, SIMPLE controls, relaxation;\n", "- initial fields: boundary conditions and starting values;\n", "- `momentumTransport`: RAS model choice, here `kOmegaSST`.\n", "\n", "`frontAndBack` is an OpenFOAM `empty` patch. It matters for declaring the case 2D, but it does not contribute active face-flux work in the solve view used below.\n" ] }, { "cell_type": "code", "execution_count": 2, "id": "eadd8926", "metadata": { "execution": { "iopub.execute_input": "2026-07-24T06:27:09.065740Z", "iopub.status.busy": "2026-07-24T06:27:09.065631Z", "iopub.status.idle": "2026-07-24T06:27:09.069646Z", "shell.execute_reply": "2026-07-24T06:27:09.069229Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "case\n", "{'solver': 'incompressibleFluid', 'endTime': '1', 'deltaT': '1', 'RAS_model': 'kOmegaSST', 'freestream_speed': 93.213, 'alpha_deg': 3.79, 'Re': 5975192.3}\n", "\n", "patches\n", "{'name': 'aerofoil', 'type': 'wall', 'active_faces_in_stepper_view': 1026}\n", "{'name': 'freestream', 'type': 'patch', 'active_faces_in_stepper_view': 1736}\n", "{'name': 'frontAndBack', 'type': 'empty', 'active_faces_in_stepper_view': 0}\n", "frontAndBack is an empty 2D patch; it has no active solve-flux faces here.\n", "\n", "numerics\n", "{'has_SIMPLE_block': True, 'has_div_phi_U_scheme': True, 'has_laplacian_schemes': True, 'has_relaxationFactors': True}\n" ] } ], "source": [ "def assignment(text, key, default=None):\n", " m = re.search(rf\"^\\s*{re.escape(key)}\\s+([^;]+);\", text, flags=re.MULTILINE)\n", " return default if m is None else m.group(1).strip()\n", "\n", "\n", "def contains_line(text, needle):\n", " return any(needle in line for line in text.splitlines())\n", "\n", "control = (CASE / \"system/controlDict\").read_text(errors=\"replace\")\n", "schemes = (CASE / \"system/fvSchemes\").read_text(errors=\"replace\")\n", "solution = (CASE / \"system/fvSolution\").read_text(errors=\"replace\")\n", "transport = (CASE / \"constant/momentumTransport\").read_text(errors=\"replace\")\n", "\n", "patch_rows = []\n", "for patch in patches:\n", " patch_rows.append((patch.name, patch.type, int(patch.size), int(patch.start)))\n", "\n", "print(\"case\")\n", "print({\n", " \"solver\": assignment(control, \"solver\"),\n", " \"endTime\": assignment(control, \"endTime\"),\n", " \"deltaT\": assignment(control, \"deltaT\"),\n", " \"RAS_model\": assignment(transport, \"model\"),\n", " \"freestream_speed\": meta.u_inf,\n", " \"alpha_deg\": round(meta.alpha_deg, 6),\n", " \"Re\": round(meta.reynolds, 1),\n", "})\n", "\n", "print(\"\\npatches\")\n", "for row in patch_rows:\n", " print({\"name\": row[0], \"type\": row[1], \"active_faces_in_stepper_view\": row[2]})\n", "print(\"frontAndBack is an empty 2D patch; it has no active solve-flux faces here.\")\n", "\n", "print(\"\\nnumerics\")\n", "print({\n", " \"has_SIMPLE_block\": \"SIMPLE\" in solution,\n", " \"has_div_phi_U_scheme\": contains_line(schemes, \"div(phi,U)\"),\n", " \"has_laplacian_schemes\": \"laplacianSchemes\" in schemes,\n", " \"has_relaxationFactors\": \"relaxationFactors\" in solution,\n", "})\n" ] }, { "cell_type": "markdown", "id": "943203bb", "metadata": {}, "source": [ "## 2. The actual OpenFOAM solve order\n", "\n", "For this repo's OpenFOAM v14 path, `simpleFoam` maps to `foamRun -solver incompressibleFluid`.\n", "The relevant algorithmic order is:\n", "\n", "```text\n", "pre_solve\n", "advance_time\n", "begin SIMPLE/PIMPLE iteration\n", " fv_models_correct\n", " pre_predictor\n", " momentum_transport_predictor\n", " assemble momentum terms\n", " assemble UEqn\n", " relax UEqn\n", " constrain UEqn\n", " solve momentum predictor\n", " compute pressure inputs: rAU, HbyA, phiHbyA\n", " assemble pEqn\n", " solve pEqn\n", " correct phi, p, U\n", " momentum_transport_corrector # k-omega SST: omega, k, nut\n", "end iteration\n", "post_solve\n", "```\n", "\n", "Parallel intuition:\n", "\n", "```text\n", "assembly is mostly face/cell parallel\n", "linear solves are global iterative kernels\n", "boundary conditions are patch kernels\n", "residuals/convergence are reductions\n", "phase boundaries are synchronization points\n", "```\n" ] }, { "cell_type": "code", "execution_count": 3, "id": "d13855ca", "metadata": { "execution": { "iopub.execute_input": "2026-07-24T06:27:09.070871Z", "iopub.status.busy": "2026-07-24T06:27:09.070736Z", "iopub.status.idle": "2026-07-24T06:27:09.080021Z", "shell.execute_reply": "2026-07-24T06:27:09.079575Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "{'initial_continuity_residual': {'Linf': 0.061028931972296335, 'L1': 20.391531263432356}, 'note': 'Signed internal owner/neighbour flux plus boundary patch flux into each cell.'}\n" ] } ], "source": [ "def continuity_residual_from_phi(phi_field, owner, neighbour, patches, n_cells):\n", " \"\"\"Signed flux sum per cell, including boundary patch fluxes.\"\"\"\n", " internal_phi = np.asarray(phi_field.internal)\n", " r = np.zeros(n_cells, dtype=internal_phi.dtype)\n", "\n", " np.add.at(r, owner, internal_phi[:len(owner)])\n", " np.add.at(r, neighbour, -internal_phi[:len(neighbour)])\n", "\n", " for patch in patches:\n", " values = np.asarray(phi_field.boundary[patch.name].values).reshape(-1)\n", " if values.size == 0:\n", " continue\n", " face_cells = np.asarray(patch.face_cells, dtype=np.int64)[:values.size]\n", " np.add.at(r, face_cells, values)\n", " return r\n", "\n", "\n", "def norm_report(x):\n", " x = np.asarray(x)\n", " return {\n", " \"Linf\": float(np.max(np.abs(x))) if x.size else 0.0,\n", " \"L1\": float(np.sum(np.abs(x))) if x.size else 0.0,\n", " }\n", "\n", "mass0 = continuity_residual_from_phi(phi_field0, owner, neighbour, patches, len(V))\n", "print({\n", " \"initial_continuity_residual\": norm_report(mass0),\n", " \"note\": \"Signed internal owner/neighbour flux plus boundary patch flux into each cell.\",\n", "})\n", "\n" ] }, { "cell_type": "markdown", "id": "3f5ddd40", "metadata": {}, "source": [ "## 3. Momentum predictor\n", "\n", "OpenFOAM operation, stripped to computation:\n", "\n", "```text\n", "UEqn = ddt(U) + div(phi, U) + turbulent_stress_divergence(U, nut) == sources\n", "relax UEqn\n", "apply matrix/boundary constraints\n", "solve UEqn == -grad(p)\n", "```\n", "\n", "Sparse-row view:\n", "\n", "$$\n", "A_c U_c + \\sum_{n\\in N(c)} A_{cn}U_n = b_c - V_c(\\nabla p)_c\n", "$$\n", "\n", "Arrays touched:\n", "\n", "```text\n", "read: U, p, phi, nut, V, Sf, owner, neighbour, boundary fields, fvSchemes, fvSolution\n", "write: UEqn.diag, UEqn.upper/lower, UEqn.source, possibly U\n", "```\n", "\n", "Parallel shape:\n", "\n", "```text\n", "per-face work: convection/diffusion coupling across owner-neighbour\n", "per-cell work: diagonal/source accumulation and explicit terms\n", "solver work: sparse matrix-vector iterations plus reductions\n", "hazard: face contributions scatter into two cells unless accumulation is organized carefully\n", "```\n" ] }, { "cell_type": "code", "execution_count": 4, "id": "57ce4919", "metadata": { "execution": { "iopub.execute_input": "2026-07-24T06:27:09.081243Z", "iopub.status.busy": "2026-07-24T06:27:09.081153Z", "iopub.status.idle": "2026-07-24T06:27:09.322109Z", "shell.execute_reply": "2026-07-24T06:27:09.321878Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "{'phases_reached': ['pre_solve', 'advance_time', 'begin_pimple_iteration', 'fv_models_correct', 'pre_predictor', 'momentum_transport_predict', 'assemble_momentum_terms', 'assemble_UEqn'], 'momentum_terms': ['assemble_momentum_ddt', 'assemble_momentum_div_phi_U', 'assemble_momentum_divDevSigma', 'assemble_momentum_sources', 'assemble_momentum_MRF_DDt', 'assemble_momentum_pressure_rhs_grad_p'], 'diag': (283400,), 'upper': (565419,), 'lower': (565419,), 'source': (283400, 3), 'H': (283400, 3), 'diag_first_5': [4.5674015979225775, 4.426112149473814, 4.2788585852295125, 4.134855363049338, 3.9941324145628894]}\n" ] } ], "source": [ "phase_names = []\n", "for method in [\n", " stepper.pre_solve,\n", " stepper.advance_time,\n", " stepper.begin_pimple_iteration,\n", " stepper.fv_models_correct,\n", " stepper.pre_predictor,\n", " stepper.momentum_transport_predictor,\n", "]:\n", " result = q(method)\n", " phase_names.append(result.name)\n", "\n", "terms = q(stepper.assemble_momentum_terms)\n", "UEqn_result = q(stepper.assemble_momentum_matrix)\n", "UEqn = UEqn_result.outputs[\"UEqn\"]\n", "phase_names.extend([terms.name, UEqn_result.name])\n", "\n", "A_U = np.asarray(UEqn.diag)\n", "upper_U = None if UEqn.upper is None else np.asarray(UEqn.upper)\n", "lower_U = None if UEqn.lower is None else np.asarray(UEqn.lower)\n", "b_U = np.asarray(UEqn.source)\n", "H_U = None if UEqn.H is None else np.asarray(UEqn.H.internal)\n", "\n", "print({\n", " \"phases_reached\": phase_names,\n", " \"momentum_terms\": [t[\"name\"] for t in terms.outputs[\"terms\"]],\n", " \"diag\": A_U.shape,\n", " \"upper\": None if upper_U is None else upper_U.shape,\n", " \"lower\": None if lower_U is None else lower_U.shape,\n", " \"source\": b_U.shape,\n", " \"H\": None if H_U is None else H_U.shape,\n", " \"diag_first_5\": A_U[:5].tolist(),\n", "})\n" ] }, { "cell_type": "markdown", "id": "63757c15", "metadata": {}, "source": [ "## 4. One face as kernel anatomy, not serial control flow\n", "\n", "The one-face view is useful only as a stencil exemplar.\n", "It should read as:\n", "\n", "```text\n", "this is one work item shape; many faces run in parallel\n", "```\n", "\n", "For an internal face `f`, a face kernel typically reads both adjacent cells and writes contributions to both matrix rows:\n", "\n", "```text\n", "o = owner[f]\n", "n = neighbour[f]\n", "read: U[o], U[n], C[o], C[n], Sf[f], phi[f], nut[o], nut[n]\n", "write: row contribution for o and row contribution for n\n", "```\n", "\n", "Race risk appears if many faces scatter-add into the same cell row at once.\n" ] }, { "cell_type": "code", "execution_count": 5, "id": "39d467a5", "metadata": { "execution": { "iopub.execute_input": "2026-07-24T06:27:09.323799Z", "iopub.status.busy": "2026-07-24T06:27:09.323707Z", "iopub.status.idle": "2026-07-24T06:27:09.325983Z", "shell.execute_reply": "2026-07-24T06:27:09.325697Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "{'face': 282709, 'owner_cell': 141545, 'neighbour_cell': 141546, 'Sf': [-0.009744034299999982, -0.0027874465999999654, 0.0], 'magSf': 0.010134893338729693, 'phi_before': -0.9234587503553671, 'U_owner_before': [93.00914503999995, 6.161356013756317, 0.0], 'U_neighbour_before': [93.00914503999995, 6.161356013756317, 0.0], 'center_distance': 0.001573663286642463}\n" ] } ], "source": [ "f = int(len(neighbour) // 2)\n", "o = int(owner[f])\n", "n = int(neighbour[f])\n", "\n", "face_probe = {\n", " \"face\": f,\n", " \"owner_cell\": o,\n", " \"neighbour_cell\": n,\n", " \"Sf\": Sf[f].tolist(),\n", " \"magSf\": float(magSf[f]),\n", " \"phi_before\": float(phi0[f]),\n", " \"U_owner_before\": U0[o].tolist(),\n", " \"U_neighbour_before\": U0[n].tolist(),\n", " \"center_distance\": float(np.linalg.norm(C[n] - C[o])),\n", "}\n", "print(face_probe)\n" ] }, { "cell_type": "markdown", "id": "906e541d", "metadata": {}, "source": [ "## 5. Pressure correction\n", "\n", "OpenFOAM operation, stripped to computation:\n", "\n", "```text\n", "rAU = 1 / UEqn.A\n", "HbyA = constrained(rAU * UEqn.H)\n", "phiHbyA = flux(HbyA) + time/mesh correction\n", "apply pressure boundary consistency\n", "pEqn = laplacian(rAU, p) == div(phiHbyA)\n", "solve pEqn\n", "phi = phiHbyA - pEqn.flux()\n", "p.relax()\n", "U = HbyA - rAU * grad(p)\n", "correct U boundary conditions\n", "```\n", "\n", "Why this exists:\n", "\n", "```text\n", "momentum predicts U\n", "predicted U may violate continuity\n", "pressure solve computes a flux correction\n", "corrected phi is the flux field OpenFOAM uses for continuity-error accounting\n", "U is then corrected to be consistent with that pressure/flux repair\n", "```\n", "\n", "Parallel shape:\n", "\n", "```text\n", "per-cell: rAU, HbyA\n", "per-face: interpolate HbyA/rAU and compute phiHbyA\n", "solver: pressure Laplacian iterations are global\n", "per-face: corrected phi\n", "per-cell: corrected U\n", "reduction: continuity errors\n", "```\n" ] }, { "cell_type": "code", "execution_count": 6, "id": "958d0625", "metadata": { "execution": { "iopub.execute_input": "2026-07-24T06:27:09.327451Z", "iopub.status.busy": "2026-07-24T06:27:09.327360Z", "iopub.status.idle": "2026-07-24T06:27:09.655600Z", "shell.execute_reply": "2026-07-24T06:27:09.655213Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "{'matrix_steps': ['relax_UEqn', 'constrain_UEqn', 'solve_UEqn'], 'pressure_steps': ['compute_pressure_inputs', 'assemble_pEqn'], 'rAU': (283400,), 'HbyA': (283400, 3), 'phiHbyA': (565419,), 'p_diag': (283400,), 'p_upper': (565419,), 'p_source': (283400,), 'rAU_first_5': [0.00026755862048551105, 0.00024708843305873933, 0.00022836470465003848, 0.00021125319224392247, 0.00019560221178537135]}\n" ] } ], "source": [ "relax_result = q(stepper.relax_matrix)\n", "constrain_result = q(stepper.constrain_matrix)\n", "momentum_solve = q(stepper.solve_momentum)\n", "pressure_inputs = q(stepper.compute_pressure_inputs)\n", "pEqn_result = q(stepper.assemble_pressure_matrix)\n", "pEqn = pEqn_result.outputs[\"pEqn\"]\n", "\n", "rAU = np.asarray(pressure_inputs.outputs[\"rAU\"].internal)\n", "HbyA = np.asarray(pressure_inputs.outputs[\"HbyA\"].internal)\n", "phiHbyA = np.asarray(pressure_inputs.outputs[\"phiHbyA\"].internal)\n", "A_p = np.asarray(pEqn.diag)\n", "upper_p = None if pEqn.upper is None else np.asarray(pEqn.upper)\n", "b_p = np.asarray(pEqn.source)\n", "\n", "print({\n", " \"matrix_steps\": [relax_result.name, constrain_result.name, momentum_solve.name],\n", " \"pressure_steps\": [pressure_inputs.name, pEqn_result.name],\n", " \"rAU\": rAU.shape,\n", " \"HbyA\": HbyA.shape,\n", " \"phiHbyA\": phiHbyA.shape,\n", " \"p_diag\": A_p.shape,\n", " \"p_upper\": None if upper_p is None else upper_p.shape,\n", " \"p_source\": b_p.shape,\n", " \"rAU_first_5\": rAU[:5].tolist(),\n", "})\n" ] }, { "cell_type": "code", "execution_count": 7, "id": "ab30847e", "metadata": { "execution": { "iopub.execute_input": "2026-07-24T06:27:09.656945Z", "iopub.status.busy": "2026-07-24T06:27:09.656826Z", "iopub.status.idle": "2026-07-24T06:27:12.376771Z", "shell.execute_reply": "2026-07-24T06:27:12.376091Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "{'pressure_steps': ['solve_pEqn', 'correct_velocity_pressure_flux'], 'changed_fields': ['p', 'U', 'phi'], 'delta_U_Linf': 142.80037348998195, 'delta_p_Linf': 139306.85715654117, 'delta_phi_Linf': 0.12004663628376555, 'continuity_after_pressure': {'Linf': 0.026893021480056893, 'L1': 1.0664099985737665}}\n" ] } ], "source": [ "pressure_solve = q(stepper.solve_pressure)\n", "correct_result = q(stepper.correct_velocity_pressure_flux)\n", "\n", "fields_after_pressure = stepper.fields()\n", "U_after_pressure = np.asarray(fields_after_pressure[\"U\"].internal)\n", "p_after_pressure = np.asarray(fields_after_pressure[\"p\"].internal)\n", "phi_after_pressure_field = fields_after_pressure[\"phi\"]\n", "phi_after_pressure = np.asarray(phi_after_pressure_field.internal)\n", "\n", "mass_after_pressure = continuity_residual_from_phi(phi_after_pressure_field, owner, neighbour, patches, len(V))\n", "\n", "print({\n", " \"pressure_steps\": [pressure_solve.name, correct_result.name],\n", " \"changed_fields\": correct_result.changed_fields,\n", " \"delta_U_Linf\": float(np.max(np.abs(U_after_pressure - U0))),\n", " \"delta_p_Linf\": float(np.max(np.abs(p_after_pressure - p0))),\n", " \"delta_phi_Linf\": float(np.max(np.abs(phi_after_pressure - phi0))),\n", " \"continuity_after_pressure\": norm_report(mass_after_pressure),\n", "})\n" ] }, { "cell_type": "markdown", "id": "4c384cfc", "metadata": {}, "source": [ "## 6. Turbulence corrector: k-omega SST\n", "\n", "OpenFOAM's `kOmegaSST` corrector does this computationally:\n", "\n", "```text\n", "compute grad(U)\n", "compute strain magnitude S2\n", "compute production G from nut and grad(U)\n", "update omega wall coefficients\n", "compute SST blending functions F1, F2/F23\n", "assemble and solve omega equation\n", "bound omega\n", "assemble and solve k equation\n", "bound k and omega\n", "update nut = a1*k / max(a1*omega, b1*F2*sqrt(S2))\n", "apply nut boundary conditions/constraints\n", "```\n", "\n", "For parallelization:\n", "\n", "```text\n", "grad(U): face/cell stencil\n", "production/S2: per-cell\n", "omega matrix: face/cell assembly + sparse solve\n", "k matrix: face/cell assembly + sparse solve\n", "nut update: embarrassingly per-cell, plus patch kernels for wall functions\n", "```\n" ] }, { "cell_type": "code", "execution_count": 8, "id": "d8e88122", "metadata": { "execution": { "iopub.execute_input": "2026-07-24T06:27:12.378174Z", "iopub.status.busy": "2026-07-24T06:27:12.378067Z", "iopub.status.idle": "2026-07-24T06:27:12.634993Z", "shell.execute_reply": "2026-07-24T06:27:12.634592Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "{'final_steps': ['momentum_transport_correct', 'end_pimple_iteration', 'post_solve'], 'turbulence_changed_fields': ['viscosity', 'momentumTransport'], 'delta_k_Linf': 2.9020790511255652e-08, 'delta_omega_Linf': 1718225865.8661606, 'delta_nut_Linf': 3.1199999716434326e-09, 'final_continuity_residual': {'Linf': 0.026893021480056893, 'L1': 1.0664099985737665}, 'omega_delta_note': 'Linf is wall-dominated because omega wall coefficients update strongly.'}\n" ] } ], "source": [ "turbulence_result = q(stepper.momentum_transport_corrector)\n", "end_result = q(stepper.end_pimple_iteration)\n", "post_result = q(stepper.post_solve, write=False)\n", "\n", "fields1 = stepper.fields()\n", "U1 = np.asarray(fields1[\"U\"].internal)\n", "p1 = np.asarray(fields1[\"p\"].internal)\n", "phi_field1 = fields1[\"phi\"]\n", "phi1 = np.asarray(phi_field1.internal)\n", "k1 = np.asarray(fields1[\"k\"].internal)\n", "omega1 = np.asarray(fields1[\"omega\"].internal)\n", "nut1 = np.asarray(fields1[\"nut\"].internal)\n", "\n", "mass1 = continuity_residual_from_phi(phi_field1, owner, neighbour, patches, len(V))\n", "\n", "print({\n", " \"final_steps\": [turbulence_result.name, end_result.name, post_result.name],\n", " \"turbulence_changed_fields\": turbulence_result.changed_fields,\n", " \"delta_k_Linf\": float(np.max(np.abs(k1 - k0))),\n", " \"delta_omega_Linf\": float(np.max(np.abs(omega1 - omega0))),\n", " \"delta_nut_Linf\": float(np.max(np.abs(nut1 - nut0))),\n", " \"final_continuity_residual\": norm_report(mass1),\n", " \"omega_delta_note\": \"Linf is wall-dominated because omega wall coefficients update strongly.\",\n", "})\n" ] }, { "cell_type": "markdown", "id": "a78169fd", "metadata": {}, "source": [ "## 7. Parallelization map\n", "\n", "| OpenFOAM algorithm phase | natural work items | writes | synchronization / hazard |\n", "|---|---|---|---|\n", "| face flux/interpolation | faces | `phi[f]`, face temporaries | reads owner/neighbour cells |\n", "| continuity residual | faces or cells | residual per cell | face scatter needs atomics or two-pass accumulation |\n", "| momentum assembly | faces + cells | sparse matrix rows | owner/neighbour scatter or row-wise gather |\n", "| matrix relaxation/constraints | cells + patches | matrix rows, boundary coeffs | patch-specific logic |\n", "| momentum solve | sparse rows | `U` iterations | global reductions each Krylov iteration |\n", "| pressure-input build | cells + faces | `rAU`, `HbyA`, `phiHbyA` | interpolation crosses faces |\n", "| pressure assembly/solve | faces + sparse rows | pressure matrix, `p` iterations | global solve/reductions |\n", "| flux/velocity correction | faces + cells | `phi`, `U` | pressure field must be solved first |\n", "| SST corrector | cells + faces + patches | `omega`, `k`, `nut` | two scalar solves, wall-function patches |\n", "\n", "Main lesson for GPU work:\n", "\n", "```text\n", "OpenFOAM's loop is phase-serial.\n", "Inside each phase, most assembly/update work is data-parallel.\n", "The sparse solvers and reductions are the major global synchronization points.\n", "```\n" ] }, { "cell_type": "markdown", "id": "3737d955", "metadata": {}, "source": [ "## 8. What this notebook should let you explain\n", "\n", "After working through it, you should be able to say:\n", "\n", "1. which OpenFOAM phase builds the momentum matrix;\n", "2. why pressure correction exists;\n", "3. why `phi` is central to continuity;\n", "4. where boundary conditions alter matrices/fields;\n", "5. which parts are face-parallel, cell-parallel, patch-parallel, or globally synchronized;\n", "6. why a one-face/one-cell probe is a kernel stencil example, not the serial algorithm.\n", "\n", "If a future section cannot be described as `phase → arrays → kernel shape`, it is probably framework spandrel and should be cut.\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.12.12" } }, "nbformat": 4, "nbformat_minor": 5 }