zacheryasc/openFOAM-RANS-to-GPU
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions
openFOAM-RANS-to-GPU/notebooks/airfrans_openfoam_algorithm_first.ipynb

815 lines
30 KiB
Text
Raw Permalink Normal View History

init
2026-07-25 07:09:19 +00:00
{
"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
}
Reference in a new issue Copy permalink