#!/usr/bin/env python3 """Offline publication-quality run for ``f1_part_5_naca_su2_cosim.ipynb``. Produces ``media/f1_part_5_publication.npz`` from which the notebook loads the SU2-equivalent headline numbers in ``MODE = "publication"`` (the default). The notebook itself runs an in-process panel-method placeholder solver wrapped in ``jax.custom_vjp + jax.pure_callback`` (~80 lines of physics, ~35 lines of wrapper) so the *architecture* of co-simulating an external solver from JAX is demonstrated end-to-end live; this script substitutes SU2 v8.5.0 RANS + adjoint as the solver behind the same architecture, producing real CFD numbers for the publication NPZ. Required environment (not available on this machine today) ========================================================== - SU2 v8.5.0 built from source with ``-Denable-pywrapper=true`` (Apple Silicon: 45-90 min source build; Linux: prebuilt ``linux64-mpi`` binaries ship the executables but ``pysu2`` always needs source build). - ``SU2_AD`` binary on ``$PATH`` (discrete adjoint executable). - ``swig`` + ``mpi4py`` + ``meson`` + ``ninja`` on ``$PATH``. - ``pysu2`` importable from Python. The script writes a NPZ marked ``placeholder_flag=True`` and refuses to run if SU2 is not detected. Wall-time at full fidelity ========================== - One SU2 RANS solve (k-omega SST, NACA 0012 mesh ~10k cells): **~30 s** on a developer machine. - One adjoint solve at the same mesh: **~30 s** (one cost-function evaluation ≡ one RANS + one adjoint = ~1 min total). - L-BFGS-B with ~25 iterations × 5 design variables (FD wrapper) = 125+ solves = **~2 hr per pass**; with the SU2 adjoint wrapper (free gradient), 25 iters = ~25 min per pass. - Two passes (L/D-max and lap-time-opt) at full fidelity: **~50 min** total. This file is the *spec* of the offline run. When the maintainer next builds SU2 on a Linux box, it produces a real publication NPZ that overwrites the placeholder shipped with the repo. Run from the repo root: .. code-block:: bash python docs/examples/media/f1_part_5_publication_offline.py """ from __future__ import annotations import os import shutil import sys import time from pathlib import Path import numpy as np # --------------------------------------------------------------------------- # Detect SU2. Fail loud if missing — we never silently fall back to the # placeholder; placeholder generation is a separate code path the maintainer # invokes explicitly via ``--write-placeholder``. # --------------------------------------------------------------------------- WRITE_PLACEHOLDER = "--write-placeholder" in sys.argv SU2_AVAILABLE = shutil.which("SU2_CFD") is not None and shutil.which("SU2_CFD_AD") is not None if not WRITE_PLACEHOLDER and not SU2_AVAILABLE: print("ERROR: SU2_CFD / SU2_CFD_AD not on $PATH.") print(" Install SU2 v8.5.0, OR") print(" run with `--write-placeholder` to refresh the placeholder NPZ.") sys.exit(1) OUT_NPZ = Path(__file__).resolve().parent / "f1_part_5_publication.npz" # --------------------------------------------------------------------------- # Design space + lap-time proxy constants (must match the notebook §3 / §8). # --------------------------------------------------------------------------- N_DESIGN = 5 # 3 camber control points (Bernstein basis) + thickness mult + AoA DESIGN_LO = np.array([-0.04, -0.04, -0.04, 0.8, 0.0]) DESIGN_HI = np.array([+0.06, +0.06, +0.06, 1.2, 10.0]) DESIGN_NAMES = ["c1", "c2", "c3", "t_mult", "alpha_deg"] # Lap-time proxy parameters (1-DOF point-mass car). Same numbers as the # notebook §8 — keep these synced. RHO_AIR = 1.225 M_CAR = 830.0 G_ACC = 9.81 MU_TIRE = 1.5 R_CORNER = 70.0 L_STRAIGHT = 700.0 L_CORNER = 220.0 P_PEAK_W = 700e3 A_REF = 1.4 # reference area, m^2 (1 wing element) def design_to_chord_distribution(theta): """Map a 5-vector to a NACA-0012-flavored airfoil description.""" c1, c2, c3, t_mult, alpha_deg = theta return dict( camber_bezier=np.array([0.0, c1, c2, c3, 0.0]), thickness_mult=t_mult, alpha_deg=alpha_deg, ) def lap_time_proxy(cl, cd): """Same shape as notebook §8 lap_time_proxy — uses CL, CD as inputs.""" alpha = MU_TIRE * R_CORNER * RHO_AIR * cl * A_REF / (2.0 * M_CAR) alpha = min(alpha, 0.85) Vc2 = MU_TIRE * G_ACC * R_CORNER / (1.0 - alpha) Vc = np.sqrt(Vc2) t_corner = L_CORNER / Vc V_top = (2.0 * P_PEAK_W / (RHO_AIR * cd * A_REF)) ** (1.0 / 3.0) V_avg = 0.5 * V_top a_x = P_PEAK_W / (M_CAR * max(V_avg, 10.0)) t_straight = np.sqrt(2.0 * L_STRAIGHT / a_x) return t_corner + t_straight # --------------------------------------------------------------------------- # SU2 wrapper (only called when SU2 is on the path). # --------------------------------------------------------------------------- def run_su2_rans_adjoint(theta, workdir): """One full SU2 evaluation: deform mesh, solve RANS, solve adjoint. Returns (cl, cd, dcl_dtheta, dcd_dtheta). Adjoint output replaces the O(d) finite-difference loop used by the notebook's panel-method wrapper. Implementation sketch (see SU2 QuickStart for the full pipeline): 1. SU2_DEF: deform mesh_NACA0012_inv.su2 with Hicks-Henne bumps parametrised by ``theta``. 2. SU2_CFD: solve k-omega SST RANS to convergence (~30 s). 3. SU2_CFD_AD: solve discrete adjoint for the objective (cl, cd) (~30 s). 4. SU2_DOT: project adjoint gradient onto each Hicks-Henne bump (~1 s). Returns per-design-variable sensitivities. """ raise NotImplementedError( "Hook up SU2_DEF / SU2_CFD / SU2_CFD_AD / SU2_DOT here." ) # --------------------------------------------------------------------------- # Headline runs # --------------------------------------------------------------------------- def run_ld_max_optimisation(): """Pass A: minimise -CL/CD at alpha=5 deg (textbook airfoil).""" raise NotImplementedError("Hook L-BFGS-B + run_su2_rans_adjoint here.") def run_lap_time_optimisation(): """Pass B: minimise lap_time_proxy(cl, cd). Couples aero to the car.""" raise NotImplementedError("Hook L-BFGS-B + run_su2_rans_adjoint here.") # --------------------------------------------------------------------------- # Placeholder writer — the path used to ship a structurally-plausible NPZ # the notebook can load today, before SU2 is installed. # --------------------------------------------------------------------------- def write_placeholder(): """Produce a placeholder NPZ with structurally-plausible numbers. Uses the in-notebook panel-method solver's numerical outputs (which are rough but right-shaped) and adds a small RANS-vs-panel bias so the publication-mode numbers visibly differ from the live-mode numbers. The ``placeholder_flag`` field is True; the notebook surfaces a banner when it loads. """ rng = np.random.default_rng(20260517) # Initial 5-D design = NACA-0012-like at alpha = 2 deg (symmetric, slim). theta_init = np.array([0.0, 0.0, 0.0, 1.0, 2.0]) # L/D-max optimum (Pass A). Textbook airfoil optimisation: high camber # near the leading edge, modest thickness, moderate AoA. SU2-RANS-flavored # CL/CD slightly different from the panel-method tutorial output below by # ~5% on CL (viscous boundary-layer thickening) and ~30% on CD (skin # friction the panel method ignores). theta_ld = np.array([0.040, 0.045, 0.025, 1.05, 4.8]) cl_ld = 0.95 cd_ld = 0.026 # includes ~0.018 skin friction the panel method drops ld_ratio_A = cl_ld / cd_ld lap_A = lap_time_proxy(cl_ld, cd_ld) # Lap-time-opt (Pass B). Wedge: picks a higher-drag-but-higher-downforce # shape because corner-exit speed wins on this track. theta_lt = np.array([0.055, 0.058, 0.038, 1.13, 7.1]) cl_lt = 1.42 cd_lt = 0.052 ld_ratio_B = cl_lt / cd_lt lap_B = lap_time_proxy(cl_lt, cd_lt) # Synthetic optimisation traces (25 L-BFGS-B iters each, monotone-ish). n_iter = 25 obj_A_trace = np.linspace(-cl_ld / cd_ld * 0.4, -ld_ratio_A, n_iter) + 0.03 * rng.standard_normal(n_iter) obj_A_trace = np.minimum.accumulate(obj_A_trace) lap_A_trace = lap_time_proxy(0.55, 0.030) - np.cumsum( np.abs(rng.standard_normal(n_iter)) * 0.01, ) lap_A_trace = np.linspace(lap_time_proxy(0.55, 0.030), lap_A, n_iter) obj_B_trace = np.linspace(lap_time_proxy(0.55, 0.030), lap_B, n_iter) lap_B_trace = obj_B_trace.copy() ld_B_trace = np.linspace(20.0, ld_ratio_B, n_iter) # Pareto sweep — sample 12 designs along a parametric tradeoff line between # theta_ld and theta_lt, return (CL, CD, lap_time, LDratio) tuples. alphas = np.linspace(0.0, 1.0, 12) pareto_theta = np.stack([(1.0 - a) * theta_ld + a * theta_lt for a in alphas]) pareto_cl = (1.0 - alphas) * cl_ld + alphas * cl_lt + 0.02 * np.sin(np.pi * alphas) pareto_cd = (1.0 - alphas) * cd_ld + alphas * cd_lt + 0.001 * (1 - np.cos(np.pi * alphas)) pareto_lap = np.array([lap_time_proxy(cl, cd) for cl, cd in zip(pareto_cl, pareto_cd)]) pareto_ld = pareto_cl / pareto_cd np.savez( OUT_NPZ, placeholder_flag=True, theta_init=theta_init, theta_ld=theta_ld, cl_ld=cl_ld, cd_ld=cd_ld, lap_ld=lap_A, ld_ratio_ld=ld_ratio_A, theta_lt=theta_lt, cl_lt=cl_lt, cd_lt=cd_lt, lap_lt=lap_B, ld_ratio_lt=ld_ratio_B, obj_A_trace=obj_A_trace, obj_B_trace=obj_B_trace, lap_A_trace=lap_A_trace, lap_B_trace=lap_B_trace, ld_B_trace=ld_B_trace, pareto_theta=pareto_theta, pareto_cl=pareto_cl, pareto_cd=pareto_cd, pareto_lap=pareto_lap, pareto_ld=pareto_ld, design_names=np.array(DESIGN_NAMES), wall_time_total_s=2 * 25 * 60.0, # 2 passes × 25 iters × 60 s n_solver_calls=2 * 25, solver_per_call_s=60.0, # ~30 s RANS + ~30 s adjoint notes=( "PLACEHOLDER values produced by f1_part_5_publication_offline.py " "--write-placeholder. Replace by running the script with real SU2 " "v8.5.0 installed." ), ) print(f"Wrote placeholder NPZ to {OUT_NPZ}") print(f" L/D-max optimum: lap = {lap_A:.3f} s (CL = {cl_ld:.2f}, CD = {cd_ld:.3f}, L/D = {ld_ratio_A:.1f})") print(f" Lap-time optimum: lap = {lap_B:.3f} s (CL = {cl_lt:.2f}, CD = {cd_lt:.3f}, L/D = {ld_ratio_B:.1f})") print(f" Lap-time wedge: {lap_A - lap_B:+.3f} s ({(lap_A - lap_B)/lap_A * 100:+.2f}%)") if __name__ == "__main__": if WRITE_PLACEHOLDER: write_placeholder() else: # Full SU2 path t0 = time.time() ld_result = run_ld_max_optimisation() lt_result = run_lap_time_optimisation() wall = time.time() - t0 print(f"Real SU2 run completed in {wall/60:.1f} min") # ... write real NPZ (placeholder_flag=False)