Code Implementation

from __future__ import annotations

import numpy as np


def generate_naca_4digit(
    m: float, p: float, t: float, num_points: int = 200
) -> tuple[np.ndarray, np.ndarray]:
    """
    Generate upper and lower coordinates for a NACA 4-digit airfoil.

    Uses cosine spacing to cluster points at leading and trailing edges
    for better resolution of high-curvature regions.

    Args:
        m: Maximum camber (e.g., 0.02 for NACA 2412)
        p: Position of max camber (e.g., 0.4 for NACA 2412)
        t: Maximum thickness in percent (e.g., 12 for NACA 0012)
        num_points: Number of points along the chord
    Returns:
        (upper_coords, lower_coords) as Nx2 numpy arrays
    """
    t = t / 100.0

    beta = np.linspace(0, np.pi, num_points)
    x = (1 - np.cos(beta)) / 2

    yt = (
        5 * t * (
            0.2969 * np.sqrt(x) - 0.1260 * x - 0.3516 * x**2
            + 0.2843 * x**3 - 0.1015 * x**4
        )
    )

    yc = np.zeros_like(x)
    dyc_dx = np.zeros_like(x)

    if m > 0 and p > 0:
        forward = x < p
        backward = ~forward
        yc[forward] = (m / p**2) * (2 * p * x[forward] - x[forward] ** 2)
        dyc_dx[forward] = (2 * m / p**2) * (p - x[forward])
        yc[backward] = (m / (1 - p) ** 2) * (
            (1 - 2 * p) + 2 * p * x[backward] - x[backward] ** 2
        )
        dyc_dx[backward] = (2 * m / (1 - p) ** 2) * (p - x[backward])

    theta = np.arctan(dyc_dx)

    xu = x - yt * np.sin(theta)
    yu = yc + yt * np.cos(theta)
    xl = x + yt * np.sin(theta)
    yl = yc - yt * np.cos(theta)

    upper_coords = np.column_stack((xu, yu))
    lower_coords = np.column_stack((xl, yl))

    return upper_coords, lower_coords


def save_dat_file(upper: np.ndarray, lower: np.ndarray, filename: str) -> None:
    coords = np.vstack((upper[::-1], lower[1:]))
    np.savetxt(filename, coords, fmt="%f %f")

from __future__ import annotations

from pathlib import Path
from typing import Optional

import gmsh
import numpy as np


class MeshGenerator:
    """Gmsh-based mesh generation for 2D airfoil C-grid meshes."""

    FARFIELD_RADIUS: float = 15.0
    DOWNSTREAM_LENGTH: float = 30.0

    def __init__(self, mesh_density: float = 1.0):
        self.mesh_density = mesh_density

    @staticmethod
    def _clean_coords(dat_file: str) -> np.ndarray:
        coords = np.loadtxt(dat_file)
        if np.allclose(coords[0], coords[-1], atol=1e-5):
            return coords[:-1]
        return coords

    def generate(
        self,
        dat_file: str,
        output_su2: str,
        bl_first_layer: float = 2e-5,
        bl_ratio: float = 1.15,
        bl_thickness: float = 0.05,
        farfield_size: float = 1.2,
        airfoil_size: float = 0.003,
        quiet: bool = True,
    ) -> Path:
        gmsh.initialize()
        if quiet:
            gmsh.option.setNumber("General.Terminal", 0)

        gmsh.model.add("airfoil_cgrid")

        coords = self._clean_coords(dat_file)
        airfoil_points = []
        for x, y in coords:
            pid = gmsh.model.geo.addPoint(x, y, 0.0)
            airfoil_points.append(pid)

        airfoil_tag = gmsh.model.geo.addSpline(airfoil_points + [airfoil_points[0]])

        R = self.FARFIELD_RADIUS
        L = self.DOWNSTREAM_LENGTH

        p_bot = gmsh.model.geo.addPoint(0.0, -R, 0.0)
        p_ctr = gmsh.model.geo.addPoint(0.0, 0.0, 0.0)
        p_lef = gmsh.model.geo.addPoint(-R, 0.0, 0.0)
        p_top = gmsh.model.geo.addPoint(0.0, R, 0.0)
        p_out_top = gmsh.model.geo.addPoint(L, R, 0.0)
        p_out_bot = gmsh.model.geo.addPoint(L, -R, 0.0)

        arc_bot = gmsh.model.geo.addCircleArc(p_bot, p_ctr, p_lef)
        arc_top = gmsh.model.geo.addCircleArc(p_lef, p_ctr, p_top)
        line_top = gmsh.model.geo.addLine(p_top, p_out_top)
        line_out = gmsh.model.geo.addLine(p_out_top, p_out_bot)
        line_bot = gmsh.model.geo.addLine(p_out_bot, p_bot)

        farfield_curves = [arc_bot, arc_top, line_top, line_out, line_bot]

        farfield_loop = gmsh.model.geo.addCurveLoop(farfield_curves)
        airfoil_loop = gmsh.model.geo.addCurveLoop([airfoil_tag])
        fluid_surf = gmsh.model.geo.addPlaneSurface([farfield_loop, airfoil_loop])

        gmsh.model.geo.synchronize()

        gmsh.model.addPhysicalGroup(1, farfield_curves, name="farfield")
        gmsh.model.addPhysicalGroup(1, [airfoil_tag], name="airfoil")
        gmsh.model.addPhysicalGroup(2, [fluid_surf], name="fluid")

        dist = gmsh.model.mesh.field.add("Distance")
        gmsh.model.mesh.field.setNumbers(dist, "CurvesList", [airfoil_tag])

        thresh = gmsh.model.mesh.field.add("Threshold")
        gmsh.model.mesh.field.setNumber(thresh, "InField", dist)
        gmsh.model.mesh.field.setNumber(thresh, "SizeMin", airfoil_size * self.mesh_density)
        gmsh.model.mesh.field.setNumber(thresh, "SizeMax", farfield_size * self.mesh_density)
        gmsh.model.mesh.field.setNumber(thresh, "DistMin", 0.05)
        gmsh.model.mesh.field.setNumber(thresh, "DistMax", 2.5)

        bl = gmsh.model.mesh.field.add("BoundaryLayer")
        gmsh.model.mesh.field.setNumbers(bl, "CurvesList", [airfoil_tag])
        gmsh.model.mesh.field.setNumber(bl, "Size", bl_first_layer)
        gmsh.model.mesh.field.setNumber(bl, "Ratio", bl_ratio)
        gmsh.model.mesh.field.setNumber(bl, "Thickness", bl_thickness)
        gmsh.model.mesh.field.setAsBoundaryLayer(bl)

        gmsh.model.mesh.field.setAsBackgroundMesh(thresh)

        gmsh.option.setNumber("Mesh.Algorithm", 6)
        gmsh.option.setNumber("Mesh.RecombinationAlgorithm", 0)
        gmsh.option.setNumber("Mesh.SubdivisionAlgorithm", 0)

        gmsh.model.mesh.generate(2)

        gmsh.model.mesh.createTopology()
        gmsh.write(output_su2)

        gmsh.finalize()
        return Path(output_su2)


def generate_su2_mesh(dat_file, output_su2, mesh_density=1.0):
    gen = MeshGenerator(mesh_density=mesh_density)
    return gen.generate(dat_file, output_su2)

from __future__ import annotations

import dataclasses
import shutil
import subprocess
import tempfile
from pathlib import Path
from typing import Optional


@dataclasses.dataclass
class SU2Config:
    angle_of_attack: float = 0.0
    mach_number: float = 0.15
    reynolds_number: float = 1_000_000
    reynolds_length: float = 1.0
    iterations: int = 2000
    cfl_number: float = 0.5
    cfl_adapt: bool = True
    turbulence_model: str = "SA"
    muscl_flow: str = "NO"
    slope_limiter_flow: str = "VENKATAKRISHNAN"
    conv_num_method_flow: str = "ROE"
    screen_output: str = "ITER, RMS_DENSITY, LIFT, DRAG"

    def write(self, path: Path) -> None:
        tag = str(int(self.angle_of_attack))
        lines = [
            f"% ------ CONFIG FILE (auto-generated by SU2Config) ------",
            f"",
            f"% --- SOLVER DEFINITIONS ---",
            f"SOLVER= RANS",
            f"KIND_TURB_MODEL= {self.turbulence_model}",
            f"MATH_PROBLEM= DIRECT",
            f"RESTART_SOL= NO",
            f"",
            f"% --- FREESTREAM DEFINITIONS ---",
            f"MACH_NUMBER= {self.mach_number}",
            f"AOA= {self.angle_of_attack:.2f}",
            f"REYNOLDS_NUMBER= {self.reynolds_number}",
            f"REYNOLDS_LENGTH= {self.reynolds_length}",
            f"",
            f"% --- THERMODYNAMIC PROPERTIES ---",
            f"FLUID_MODEL= STANDARD_AIR",
            f"GAMMA_VALUE= 1.4",
            f"GAS_CONSTANT= 287.05",
            f"FREESTREAM_TEMPERATURE= 288.15",
            f"FREESTREAM_PRESSURE= 101325.0",
            f"",
            f"% --- BOUNDARY CONDITIONS ---",
            f"MARKER_FAR= ( farfield )",
            f"MARKER_HEATFLUX= ( airfoil, 0.0 )",
            f"MARKER_MONITORING= ( airfoil )",
            f"",
            f"% --- NUMERICAL METHODS ---",
            f"NUM_METHOD_GRAD= GREEN_GAUSS",
            f"CONV_NUM_METHOD_FLOW= {self.conv_num_method_flow}",
            f"CONV_NUM_METHOD_TURB= SCALAR_UPWIND",
            f"MUSCL_FLOW= {self.muscl_flow}",
            f"SLOPE_LIMITER_FLOW= {self.slope_limiter_flow}",
            f"MUSCL_TURB= NO",
            f"TIME_DISCRE_FLOW= EULER_IMPLICIT",
            f"",
            f"% --- ITERATIVE CONTROLS ---",
            f"ITER= {self.iterations}",
            f"CFL_NUMBER= {self.cfl_number}",
            f"CFL_ADAPT= {'YES' if self.cfl_adapt else 'NO'}",
            f"CFL_ADAPT_PARAM= ( 0.01, 5.0, 1.1, 25.0 )",
            f"",
            f"% --- CONVERGENCE CRITERIA ---",
            f"CONV_FIELD= RMS_DENSITY",
            f"CONV_RESIDUAL_MINVAL= -8",
            f"CONV_STARTITER= 10",
            f"",
            f"% --- LINEAR SOLVER ---",
            f"LINEAR_SOLVER= FGMRES",
            f"LINEAR_SOLVER_PREC= ILU",
            f"LINEAR_SOLVER_ERROR= 1E-6",
            f"LINEAR_SOLVER_ITER= 10",
            f"",
            f"% --- MULTIGRID ---",
            f"MGLEVEL= 0",
            f"",
            f"% --- OUTPUT ---",
            f"MESH_FORMAT= SU2",
            f"MESH_FILENAME= mesh.su2",
            f"TABULAR_FORMAT= CSV",
            f"CONV_FILENAME= history",
            f"VOLUME_FILENAME= flow_results_{tag}",
            f"SURFACE_FILENAME= surface_flow_{tag}",
            f"OUTPUT_WRT_FREQ= {self.iterations}",
            f"HISTORY_WRT_FREQ_INNER= 1",
            f"SCREEN_OUTPUT= ( {self.screen_output} )",
            f"HISTORY_OUTPUT= ( INNER_ITER, RMS_RES, AERO_COEFF )",
            f"OUTPUT_FILES= (RESTART, PARAVIEW)",
            f"",
        ]
        path.write_text("\n".join(lines))


@dataclasses.dataclass
class SU2Results:
    cd: float = 0.0
    cl: float = 0.0
    cmz: float = 0.0
    converged: bool = False
    iterations: int = 0
    history: list[dict] = dataclasses.field(default_factory=list)

    def to_dict(self) -> dict:
        return dataclasses.asdict(self)


class SU2Solver:
    def __init__(self, su2_cfd: str = "SU2_CFD", workdir: Optional[Path] = None):
        self.su2_cfd = su2_cfd
        self.workdir = workdir or Path(tempfile.mkdtemp())

    def run(self, config: SU2Config, mesh_file: Path, output_dir: Path,
            timeout: Optional[int] = None) -> SU2Results:
        output_dir = Path(output_dir).resolve()
        output_dir.mkdir(parents=True, exist_ok=True)

        cfg_path = output_dir / "config.cfg"
        config.write(cfg_path)

        mesh_local = output_dir / "mesh.su2"
        if mesh_file.resolve() != mesh_local.resolve():
            shutil.copy2(mesh_file, mesh_local)

        cmd = [self.su2_cfd, str(cfg_path)]
        print(f"  Executing SU2_CFD (up to {config.iterations} iterations)...")

        try:
            proc = subprocess.run(cmd, cwd=str(output_dir),
                                  capture_output=True, text=True, timeout=timeout)
        except subprocess.TimeoutExpired:
            print("  SU2_CFD timed out")
            return SU2Results()

        return self._parse_results(proc, output_dir)

    def _normalize_header(self, raw_header: str) -> list[str]:
        return [h.strip().strip('"').strip() for h in raw_header.split(",")]

    def _parse_screen_cl_cd(self, stdout: str) -> tuple[Optional[float], Optional[float]]:
        lines = stdout.strip().split("\n")
        for line in reversed(lines):
            parts = [p.strip() for p in line.split("|")]
            if len(parts) >= 4:
                try:
                    return float(parts[-2]), float(parts[-1])
                except (ValueError, IndexError):
                    continue
        return None, None

    def _parse_results(self, proc: subprocess.CompletedProcess, workdir: Path) -> SU2Results:
        results = SU2Results()
        stdout, stderr = proc.stdout, proc.stderr

        if proc.returncode != 0:
            print(f"  [SU2 exited with code {proc.returncode}]")
            for line in (stderr or "").strip().split("\n")[-5:]:
                if line.strip():
                    print(f"  STDERR: {line.strip()}")
            return results

        results.converged = "Convergence reached" in stdout or "convergence" in stdout.lower()

        col_map = {"ITER": "ITER", "Inner_Iter": "ITER", "INNER_ITER": "ITER",
                    "DRAG": "DRAG", "CD": "DRAG", "cd": "DRAG",
                    "LIFT": "LIFT", "CL": "LIFT", "cl": "LIFT",
                    "MOMENT": "MOMENT", "CMZ": "MOMENT", "Cmz": "MOMENT"}

        hist_file = workdir / "history.csv"
        if hist_file.exists():
            lines = hist_file.read_text().strip().split("\n")
            if len(lines) > 1:
                header = self._normalize_header(lines[0])
                n_cols = len(header)
                last_entry = None
                for line in lines[1:]:
                    vals = [v.strip() for v in line.split(",")]
                    if len(vals) == n_cols:
                        entry = {col_map.get(k, k): v for k, v in zip(header, vals)}
                        results.history.append(entry)
                        last_entry = entry
                if last_entry:
                    try:
                        results.iterations = int(float(last_entry.get("ITER", 0)))
                        cd = last_entry.get("DRAG")
                        cl = last_entry.get("LIFT")
                        if cd is not None: results.cd = float(cd)
                        if cl is not None: results.cl = float(cl)
                    except (ValueError, TypeError):
                        pass

        if abs(results.cd) < 1e-12 or abs(results.cl) > 1e3:
            cl, cd = self._parse_screen_cl_cd(stdout)
            if cl is not None:
                results.cl = cl
                results.cd = cd

        return results

from __future__ import annotations

# Post-processing module is intentionally empty.
# Velocity, pressure, mesh, and Cp visualizations are created manually in ParaView.
# Convergence and aggregate plots are generated by physics/analysis.py.

from __future__ import annotations

from pathlib import Path
from typing import Optional

import numpy as np
from physics.solver import SU2Results

BG_COLOR = "#282a36"
CARD_BG = "#44475a"
FG_COLOR = "#f8f8f2"
PINK = "#ff79c6"
PURPLE = "#bd93f9"
CYAN = "#8be9fd"
GREEN = "#50fa7b"
YELLOW = "#f1fa8c"
COMMENT = "#6272a4"


def _setup_matplotlib():
    import matplotlib
    matplotlib.use("Agg")
    import matplotlib.pyplot as plt
    plt.rcParams.update({
        "figure.facecolor": BG_COLOR, "axes.facecolor": BG_COLOR,
        "axes.edgecolor": CARD_BG, "axes.labelcolor": FG_COLOR,
        "text.color": FG_COLOR, "xtick.color": COMMENT, "ytick.color": COMMENT,
        "grid.alpha": 0.15, "grid.color": FG_COLOR,
        "legend.facecolor": CARD_BG, "legend.labelcolor": FG_COLOR,
    })
    return plt


def plot_convergence(history: list[dict], save_path: str, aoa: int = 0) -> Optional[str]:
    if not history:
        return None
    import matplotlib.pyplot as plt
    _setup_matplotlib()

    n = len(history)
    iters = np.arange(1, n + 1)
    rho_res, rhoU_res, rhoV_res, nu_res = [], [], [], []
    cl_hist, cd_hist = [], []
    for entry in history:
        rho_res.append(float(entry.get("rms[Rho]", float("nan"))))
        rhoU_res.append(float(entry.get("rms[RhoU]", float("nan"))))
        rhoV_res.append(float(entry.get("rms[RhoV]", float("nan"))))
        nu_res.append(float(entry.get("rms[nu]", float("nan"))))
        cl_hist.append(float(entry.get("LIFT", float("nan"))))
        cd_hist.append(float(entry.get("DRAG", float("nan"))))

    fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 8))
    ax1.plot(iters, rho_res, label="RMS Density", color=CYAN, linewidth=1)
    ax1.plot(iters, rhoU_res, label="RMS RhoU", color=PINK, linewidth=1)
    ax1.plot(iters, rhoV_res, label="RMS RhoV", color=GREEN, linewidth=1)
    ax1.plot(iters, nu_res, label="RMS nu", color=YELLOW, linewidth=1)
    ax1.set_ylabel("Log10(Residual)")
    ax1.set_title("Convergence History - RMS Residuals")
    ax1.legend()
    ax1.grid(True, alpha=0.15)

    ax2.plot(iters, cl_hist, label="$C_l$", color=PURPLE, linewidth=1.5)
    ax2.plot(iters, cd_hist, label="$C_d$", color=PINK, linewidth=1.5)
    ax2.set_xlabel("Iteration")
    ax2.set_ylabel("Coefficient")
    ax2.set_title("Force Coefficient Convergence")
    ax2.legend()
    ax2.grid(True, alpha=0.15)

    fig.tight_layout()
    out = f"{save_path}/convergence_{aoa}.png"
    fig.savefig(out, dpi=150, bbox_inches="tight")
    plt.close(fig)
    return out


def plot_cl_alpha(results: list[tuple[float, SU2Results]], save_path: str,
                  experimental: Optional[list[tuple[float, float]]] = None) -> Optional[str]:
    if not results:
        return None
    import matplotlib.pyplot as plt
    _setup_matplotlib()

    angles = [r[0] for r in results]
    cl_vals = [r[1].cl for r in results]

    fig, ax = plt.subplots(figsize=(8, 6))
    ax.plot(angles, cl_vals, "o-", color=PURPLE, linewidth=2, markersize=8, label="SU2 RANS (SA)")
    if experimental:
        exp_angles, exp_cl = zip(*experimental)
        ax.plot(exp_angles, exp_cl, "s--", color=COMMENT, linewidth=1.5, markersize=5, label="Experimental")
    ax.set_xlabel("Angle of Attack (deg)")
    ax.set_ylabel("Lift Coefficient $C_l$")
    ax.set_title("Lift Curve - NACA 0012, Re = 1x10^6, M = 0.15")
    ax.legend()
    ax.grid(True, alpha=0.15)
    fig.tight_layout()
    out = f"{save_path}/cl_vs_alpha.png"
    fig.savefig(out, dpi=150, bbox_inches="tight")
    plt.close(fig)
    return out


def plot_cd_alpha(results: list[tuple[float, SU2Results]], save_path: str,
                  experimental: Optional[list[tuple[float, float]]] = None) -> Optional[str]:
    if not results:
        return None
    import matplotlib.pyplot as plt
    _setup_matplotlib()

    angles = [r[0] for r in results]
    cd_vals = [r[1].cd for r in results]

    fig, ax = plt.subplots(figsize=(8, 6))
    ax.plot(angles, cd_vals, "o-", color=PINK, linewidth=2, markersize=8, label="SU2 RANS (SA)")
    if experimental:
        exp_angles, exp_cd = zip(*experimental)
        ax.plot(exp_angles, exp_cd, "s--", color=COMMENT, linewidth=1.5, markersize=5, label="Experimental")
    ax.set_xlabel("Angle of Attack (deg)")
    ax.set_ylabel("Drag Coefficient $C_d$")
    ax.set_title("Drag Polar - NACA 0012, Re = 1x10^6, M = 0.15")
    ax.legend()
    ax.grid(True, alpha=0.15)
    fig.tight_layout()
    out = f"{save_path}/cd_vs_alpha.png"
    fig.savefig(out, dpi=150, bbox_inches="tight")
    plt.close(fig)
    return out


def plot_drag_polar(results: list[tuple[float, SU2Results]], save_path: str,
                     experimental: Optional[list[tuple[float, float]]] = None) -> Optional[str]:
    if not results:
        return None
    import matplotlib.pyplot as plt
    _setup_matplotlib()

    cl_vals = [r[1].cl for r in results]
    cd_vals = [r[1].cd for r in results]
    angles = [r[0] for r in results]

    fig, ax = plt.subplots(figsize=(8, 6))
    sc = ax.scatter(cd_vals, cl_vals, c=angles, cmap="coolwarm", s=100, zorder=5)
    ax.plot(cd_vals, cl_vals, color=COMMENT, linewidth=1, alpha=0.5, zorder=3)
    cbar = fig.colorbar(sc, ax=ax, label="AoA (deg)")
    if experimental:
        exp_cd, exp_cl = zip(*experimental)
        ax.plot(exp_cd, exp_cl, "s--", color=COMMENT, linewidth=1.5, markersize=5, label="Experimental")
    for cd, cl, aoa in zip(cd_vals, cl_vals, angles):
        ax.annotate(f"{aoa}deg", (cd, cl), xytext=(5, 5),
                    textcoords="offset points", color=FG_COLOR, fontsize=9)
    ax.set_xlabel("Drag Coefficient $C_d$")
    ax.set_ylabel("Lift Coefficient $C_l$")
    ax.set_title("Drag Polar - NACA 0012, Re = 1x10^6, M = 0.15")
    ax.legend()
    ax.grid(True, alpha=0.15)
    fig.tight_layout()
    out = f"{save_path}/drag_polar.png"
    fig.savefig(out, dpi=150, bbox_inches="tight")
    plt.close(fig)
    return out

from __future__ import annotations

import shutil
from pathlib import Path

from physics.geometry import generate_naca_4digit, save_dat_file
from physics.mesher import MeshGenerator
from physics.solver import SU2Config, SU2Solver
from physics.analysis import (
    plot_convergence,
    plot_cl_alpha,
    plot_cd_alpha,
    plot_drag_polar,
)
from physics.validate import get_cl_alpha, get_cd_alpha, get_drag_polar

NACA_NAME = "NACA 0012"
NACA_PARAMS = (0, 0, 12)
ANGLES_OF_ATTACK = [0, 4, 8, 12, 16]
BASE_DIR = Path("./output/cfd/naca0012")
DOCS_IMG = Path("./docs/assets/images/naca0012")

REGIME_LABELS = {
    0: "Symmetric Baseline",
    4: "Linear Lift",
    8: "High Lift",
    12: "Onset of Stall",
    16: "Deep Stall",
}


def main():
    if BASE_DIR.exists():
        shutil.rmtree(BASE_DIR)
    BASE_DIR.mkdir(parents=True)

    solver = SU2Solver()
    mesher = MeshGenerator(mesh_density=1.0)

    all_results: list[tuple[float, SU2Results]] = []

    for aoa in ANGLES_OF_ATTACK:
        print(f"\n{'='*60}")
        print(f"  AoA = {aoa} - {REGIME_LABELS[aoa]}")
        print(f"{'='*60}")

        aoa_dir = BASE_DIR / f"aoa_{aoa}"
        aoa_dir.mkdir(parents=True, exist_ok=True)

        print("  [1/4] Generating geometry...")
        upper, lower = generate_naca_4digit(*NACA_PARAMS)
        dat_path = aoa_dir / "airfoil.dat"
        save_dat_file(upper, lower, str(dat_path))

        print("  [2/4] Generating C-grid mesh with boundary layers...")
        mesh_path = aoa_dir / "mesh.su2"
        mesher.generate(str(dat_path), str(mesh_path))

        print("  [3/4] Running SU2 CFD solver...")
        config = SU2Config(angle_of_attack=aoa)
        results = solver.run(config, mesh_path, aoa_dir, timeout=600)

        if results.history:
            print(f"    Iterations: {results.iterations}  CL={results.cl:.6f}  CD={results.cd:.6f}  Converged: {results.converged}")
        else:
            print("    WARNING: No convergence history parsed")

        all_results.append((aoa, results))

        print("  [4/4] Rendering convergence plot...")
        aoa_img_dir = DOCS_IMG / f"aoa_{aoa}"
        aoa_img_dir.mkdir(parents=True, exist_ok=True)
        plot_convergence(results.history, str(aoa_img_dir), aoa=aoa)

    print(f"\n{'='*60}")
    print("  Generating aggregate Cl/Cd curves...")
    print(f"{'='*60}")

    experimental_cl = get_cl_alpha()
    experimental_cd = get_cd_alpha()
    experimental_polar = get_drag_polar()

    DOCS_IMG.mkdir(parents=True, exist_ok=True)
    plot_cl_alpha(all_results, str(DOCS_IMG), experimental=experimental_cl)
    plot_cd_alpha(all_results, str(DOCS_IMG), experimental=experimental_cd)
    plot_drag_polar(all_results, str(DOCS_IMG), experimental=experimental_polar)

    print(f"\n  Done. Open docs/index.html in your browser.")


if __name__ == "__main__":
    main()