pip install "flow360>=25.10.4"

import csv
import os
import re

import flow360 as fl
import numpy as np
import matplotlib.pyplot as plt

from flow360.examples import download_benchmark_assets

download_benchmark_assets("backward_step", "ref_data")

root_asset_files = download_benchmark_assets("backward_step", "root_assets")
# The snapshot bundles the mesh with solver metadata and logs; select the
# volume mesh file (.cgns or .ugrid, possibly .zst-compressed) to load.
volume_mesh_file = next(
    f for f in root_asset_files if f.removesuffix(".zst").endswith((".cgns", ".ugrid"))
)
project = fl.Project.from_volume_mesh(volume_mesh_file, name="Backward-Facing Step")
volume_mesh = project.volume_mesh

# Case constants (Driver & Seegmiller reference conditions).
RE = 36000
MACH = 0.128
TEMPERATURE = 537 * fl.u.R
ALPHA = 0 * fl.u.deg
BETA = 0 * fl.u.deg
STATIC_PRESSURE_RATIO_OUTFLOW_W_FREESTREAM = 1.0023
STATIC_PRESSURE_RATIO_OUTFLOW_W_INFLOW = 1.0017

# Freestream turbulence intensity from NASA TMR (TI = 0.061 %). Combined with the
# SST viscosity ratio below it pins the (k, omega) inlet pair to match TMR/CFL3D,
# instead of letting Flow360 derive k from its (much smaller) default TI.
TURB_INTENSITY = {
    "sst": 0.061 * 0.01,
}

TOTAL_PRESSURE_RATIO_INFLOW = 1.011516
TOTAL_TEMPERATURE_RATIO_INFLOW = 1.003277

TURB_VISCOSITY_RATIO = {
    "sa": 0.210438,
    "sst": 0.009,
}

operating_condition = fl.AerospaceCondition.from_mach_reynolds(
    mach=MACH,
    reynolds_mesh_unit=RE,
    project_length_unit=1 * fl.u.m,
    temperature=TEMPERATURE,
    alpha=ALPHA,
    beta=BETA,
    reference_mach=MACH,
)

# Navier-Stokes solver (CFL3D/FUN3D-matched setup, locked across all cases):
# standard compressible NS with a Krylov linear solver + line search + low-Mach
# preconditioner, evaluating the equations and rebuilding the Jacobian every
# pseudo step.
ns_solver = fl.NavierStokesSolver(
    absolute_tolerance=1e-12,
    low_mach_preconditioner=True,
    equation_evaluation_frequency=1,
    update_jacobian_frequency=1,
    linear_solver=fl.KrylovLinearSolver(
        max_iterations=50,
        max_preconditioner_iterations=50,
        relative_tolerance=1e-3,
    ),
    line_search=fl.LineSearch(),
)

# Turbulence-model variations.
turbulence_sa = fl.SpalartAllmaras(
    absolute_tolerance=1e-10,
    equation_evaluation_frequency=1,
    update_jacobian_frequency=1,
)

turbulence_sst = fl.KOmegaSST(
    absolute_tolerance=1e-10,
    equation_evaluation_frequency=1,
    update_jacobian_frequency=1,
)

t_solver_variations = {
    "sa": turbulence_sa,
    "sst": turbulence_sst,
}

# Wall / symmetry boundary conditions common to every case.
bcs_base = [
    fl.SlipWall(
        surfaces=[
            volume_mesh["Block/Symm1"],
            volume_mesh["Block/Symm2"],
            volume_mesh["Block/Slip1"],
            volume_mesh["Block/Slip2"],
        ]
    ),
    fl.Wall(
        surfaces=[
            volume_mesh["Block/LowerWall"],
            volume_mesh["Block/UpperWall"],
        ]
    ),
]


def make_bc_variations(turb_visc_ratio, turb_intensity=None):
    if turb_intensity is not None:
        turb_q = fl.TurbulenceQuantities(
            viscosity_ratio=turb_visc_ratio,
            turbulent_intensity=turb_intensity,
        )
    else:
        turb_q = fl.TurbulenceQuantities(viscosity_ratio=turb_visc_ratio)

    bcs_freestream_outflow = [
        fl.Freestream(
            surfaces=[volume_mesh["Block/Inlet"]],
            turbulence_quantities=turb_q,
        ),
        fl.Outflow(
            surfaces=[volume_mesh["Block/Outlet"]],
            spec=fl.Pressure(
                value=operating_condition.thermal_state.pressure
                * STATIC_PRESSURE_RATIO_OUTFLOW_W_FREESTREAM
            ),
        ),
    ]

    bcs_inflow_outflow = [
        fl.Inflow(
            surfaces=[volume_mesh["Block/Inlet"]],
            spec=fl.TotalPressure(
                value=operating_condition.thermal_state.pressure * TOTAL_PRESSURE_RATIO_INFLOW
            ),
            total_temperature=TOTAL_TEMPERATURE_RATIO_INFLOW
            * operating_condition.thermal_state.temperature,
            turbulence_quantities=turb_q,
        ),
        fl.Outflow(
            surfaces=[volume_mesh["Block/Outlet"]],
            spec=fl.Pressure(
                value=operating_condition.thermal_state.pressure
                * STATIC_PRESSURE_RATIO_OUTFLOW_W_INFLOW
            ),
        ),
    ]

    return {
        "freestream_outflow": bcs_freestream_outflow,
        "inflow_outflow": bcs_inflow_outflow,
    }

# Probe arrays: the bottom wall (before and after the step) and vertical profiles.
bottom_before_step = fl.PointArray(
    name="bottom_before_step",
    start=(-4, 0.5, 1) * fl.u.m,
    end=(0, 0.5, 1) * fl.u.m,
    number_of_points=40,
)

bottom_after_step = fl.PointArray(
    name="bottom_after_step",
    start=(1e-3, 0.5, 0) * fl.u.m,
    end=(30, 0.5, 0) * fl.u.m,
    number_of_points=300,
)

profile_0 = fl.PointArray(
    name="profile_x_-4",
    start=(-4, 0.5, 1) * fl.u.m,
    end=(-4, 0.5, 6) * fl.u.m,
    number_of_points=200,
)

profile_1 = fl.PointArray(
    name="profile_x_1",
    start=(1, 0.5, 0) * fl.u.m,
    end=(1, 0.5, 5) * fl.u.m,
    number_of_points=200,
)

profile_2 = fl.PointArray(
    name="profile_x_4",
    start=(4, 0.5, 0) * fl.u.m,
    end=(4, 0.5, 5) * fl.u.m,
    number_of_points=200,
)

profile_x6 = fl.PointArray(
    name="profile_x_6",
    start=(6, 0.5, 0) * fl.u.m,
    end=(6, 0.5, 5) * fl.u.m,
    number_of_points=200,
)

profile_3 = fl.PointArray(
    name="profile_x_10",
    start=(10, 0.5, 0) * fl.u.m,
    end=(10, 0.5, 5) * fl.u.m,
    number_of_points=200,
)

# Normalized user variables (u, v, |vel|, Reynolds shear stress, streamwise Cf).
u_norm = fl.UserVariable(
    name="u_norm",
    value=fl.solution.velocity[0] / operating_condition.reference_velocity_magnitude,
)

v_norm = fl.UserVariable(
    name="v_norm",
    value=fl.solution.velocity[2] / operating_condition.reference_velocity_magnitude,
)

mag_norm = fl.UserVariable(
    name="vel_mag_norm",
    value=fl.math.sqrt(u_norm**2 + v_norm**2),
)

nu_t = fl.solution.mut / fl.solution.density
uv_norm = fl.UserVariable(
    name="uv_norm",
    value=-nu_t
    * (fl.solution.grad_u[2] + fl.solution.grad_w[0])
    * 1000
    / operating_condition.reference_velocity_magnitude**2,
)

cf_x = fl.UserVariable(
    name="Cf_X",
    value=fl.solution.CfVec[0],
)

with fl.SI_unit_system:
    outputs = [
        fl.SurfaceProbeOutput(
            probe_points=[bottom_before_step, bottom_after_step],
            output_fields=["Cp", cf_x],
            name="surface_probes_bottom",
            target_surfaces=[volume_mesh["Block/LowerWall"]],
        ),
        fl.SurfaceProbeOutput(
            probe_points=[fl.Point(location=(40, 0.5, 0) * fl.u.m, name="Cp_ref")],
            output_fields=["Cp"],
            name="Cp_shift_ref",
            target_surfaces=[volume_mesh["Block/LowerWall"]],
        ),
        fl.ProbeOutput(
            probe_points=[
                profile_0,
                profile_1,
                profile_2,
                profile_x6,
                profile_3,
            ],
            output_fields=[u_norm, v_norm, uv_norm],
        ),
        fl.ProbeOutput(
            probe_points=[fl.Point(location=(-4, 0.5, 4) * fl.u.m, name="U_ref")],
            output_fields=[u_norm, "Mach", "Cp"],
            name="ref_location",
        ),
        fl.SliceOutput(
            slices=[
                fl.Slice(
                    normal=(0, 1, 0),
                    origin=(0, 0.5, 0) * fl.u.m,
                    name="Yslice",
                )
            ],
            output_fields=[u_norm, v_norm, mag_norm, "Cp", uv_norm, "velocity"],
        ),
    ]

with fl.SI_unit_system:
    params = fl.SimulationParams(
        reference_geometry=fl.ReferenceGeometry(
            area=1 * fl.u.m**2,
            moment_center=(0, 0, 0) * fl.u.m,
            moment_length=(1, 1, 1) * fl.u.m,
        ),
        operating_condition=operating_condition,
        time_stepping=fl.Steady(
            max_steps=15000,
            CFL=fl.AdaptiveCFL(),
        ),
        outputs=outputs,
    )

SOLVER_VERSION = os.environ.get("SOLVER_VERSION_OVERRIDE", "release-25.9")

cases = []
for t_var_name, t_var in t_solver_variations.items():
    bc_variations = make_bc_variations(
        TURB_VISCOSITY_RATIO[t_var_name],
        turb_intensity=TURB_INTENSITY.get(t_var_name),
    )
    for bc_var_name, bc_var in bc_variations.items():
        case_name = f"BackwardStep_{bc_var_name}_{t_var_name}_{SOLVER_VERSION}"
        params.models = [
            fl.Fluid(
                navier_stokes_solver=ns_solver,
                turbulence_model_solver=t_var,
            ),
            *bc_var,
            *bcs_base,
        ]
        case = project.run_case(
            params=params,
            name=case_name,
            tags=[bc_var_name, t_var_name, SOLVER_VERSION],
            solver_version=SOLVER_VERSION,
        )
        print(f"Case {case_name} submitted: {case.id}")
        cases.append({"name": case_name, "case": case})

for record in cases:
    record["case"].wait()

RESULTS_DIR = "results"
REF_DATA_DIR = "ref_data"
os.makedirs(RESULTS_DIR, exist_ok=True)

STEP_HEIGHT = 1.0  # H in meters

# Shared benchmark colors: Flow360 = brand green, references = black / grey.
FLOW360_COLOR = "#00643c"
REFERENCE_COLOR = "black"

# Flow360 curves get distinct saturated colors and are drawn on top; the CFL3D
# references are thinner dashed lines underneath.
FLOW360_MODEL_COLORS = {"sa": "#1f77b4", "sst": FLOW360_COLOR}  # blue (SA), brand green (SST)
BC_LINESTYLES = {"freestream_outflow": "-", "inflow_outflow": (0, (1, 1))}
BC_LABELS = {"freestream_outflow": "freestream", "inflow_outflow": "inflow"}
MODEL_LABELS = {"sa": "SA", "sst": "SST"}

CFL3D_STYLES = {
    "sa": {"color": "#9467bd", "linestyle": "--", "linewidth": 1.2, "alpha": 0.85, "zorder": 2},
    "sst": {"color": "#d62728", "linestyle": "--", "linewidth": 1.2, "alpha": 0.85, "zorder": 2},
}
CFL3D_MARKER = "x"
CFL3D_MARKERSIZE = 6

# x/H window shown on the Cf and Cp distribution plots.
CF_CP_XLIM = (-4, 30)

# Driver & Seegmiller published experimental reattachment length (NASA TMR).
EXP_LABEL = "Experiment (Driver & Seegmiller)"
EXP_REATTACH_XH = 6.26

# Downstream profile stations: color encodes the station, line style the source/BC.
PROFILE_STATIONS = {
    1: {"monitor": "profile_x_1", "color": "#d62728"},   # red
    4: {"monitor": "profile_x_4", "color": "#1f77b4"},   # blue
    6: {"monitor": "profile_x_6", "color": "#2ca02c"},   # green
    10: {"monitor": "profile_x_10", "color": "#000000"}, # black
}
PROFILE_SERIES_STYLES = {
    "exp": {"linestyle": "none", "marker": "s", "markersize": 5,
            "markerfacecolor": "none", "zorder": 6},
    "cfl3d": {"linestyle": "--", "linewidth": 1.3, "zorder": 3},
    "freestream": {"linestyle": "-", "linewidth": 1.6, "zorder": 5},
    "inflow": {"linestyle": "-.", "linewidth": 1.6, "zorder": 4},
}
MODEL_TITLES = {"sa": "SA", "sst": "k-$\\omega$ SSTm"}


def setup_plot_style():
    plt.rcParams.update({
        "figure.figsize": (10, 7),
        "font.size": 11,
        "lines.linewidth": 1.5,
        "axes.grid": True,
        "grid.alpha": 0.3,
    })


setup_plot_style()

def load_cf_exp():
    rows = []
    with open(os.path.join(REF_DATA_DIR, "cf.exp.dat")) as f:
        for line in f:
            line = line.strip()
            if not line or line.startswith("#") or line.startswith("variables"):
                continue
            rows.append([float(p) for p in line.split()])
    arr = np.array(rows)
    return {"x_H": arr[:, 0], "cf": arr[:, 1], "error": arr[:, 2]}


def load_cp_exp():
    rows = []
    in_bottom = False
    with open(os.path.join(REF_DATA_DIR, "cp.expnew.dat")) as f:
        for line in f:
            line = line.strip()
            if not line or line.startswith("#") or line.startswith("variables"):
                continue
            if line.startswith("zone"):
                in_bottom = "bottom wall" in line
                continue
            if in_bottom:
                rows.append([float(p) for p in line.split()])
    arr = np.array(rows)
    return {"x_H": arr[:, 0], "cp": arr[:, 1], "cp_orig": arr[:, 2]}


def load_profiles_exp():
    """Load experimental profile data, keyed by x/H station."""
    profiles = {}
    current_xH = None
    rows = []
    col_names = ["ID", "y", "u", "v", "uu", "vv", "uv", "uuu", "uvv", "vuu", "vvv"]
    with open(os.path.join(REF_DATA_DIR, "profiles.exp.dat")) as f:
        for line in f:
            line = line.strip()
            if not line or line.startswith("#") or line.startswith("variables"):
                continue
            if line.startswith("zone"):
                if current_xH is not None and rows:
                    arr = np.array(rows)
                    profiles[current_xH] = {col_names[i]: arr[:, i] for i in range(arr.shape[1])}
                match = re.search(r"x/H=([-\d.]+)", line)
                current_xH = float(match.group(1)) if match else None
                rows = []
                continue
            rows.append([float(p) for p in line.split()])
    if current_xH is not None and rows:
        arr = np.array(rows)
        profiles[current_xH] = {col_names[i]: arr[:, i] for i in range(arr.shape[1])}
    return profiles


def load_cfl3d_xy(filename):
    """Load a simple two-column CFL3D file (cf or cp)."""
    rows = []
    with open(os.path.join(REF_DATA_DIR, filename)) as f:
        for line in f:
            line = line.strip()
            if not line or line.startswith("#") or line.lower().startswith("variables"):
                continue
            rows.append([float(p) for p in line.split()])
    arr = np.array(rows)
    return arr[:, 0], arr[:, 1]  # x/H, value


def load_cfl3d_zones(filename):
    """Load a zone-based CFL3D file (vel or uv). Returns dict keyed by x/H float."""
    profiles = {}
    current_xH = None
    rows = []
    with open(os.path.join(REF_DATA_DIR, filename)) as f:
        for line in f:
            line = line.strip()
            if not line or line.startswith("#") or line.upper().startswith("VARIABLES"):
                continue
            if line.upper().startswith("ZONE"):
                if current_xH is not None and rows:
                    arr = np.array(rows)
                    profiles[current_xH] = {"y_H": arr[:, 1], "value": arr[:, 2]}
                match = re.search(r"x/H=(?:approx\s*)?([-\d.]+)", line)
                current_xH = float(match.group(1)) if match else None
                rows = []
                continue
            rows.append([float(p) for p in line.split()])
    if current_xH is not None and rows:
        arr = np.array(rows)
        profiles[current_xH] = {"y_H": arr[:, 1], "value": arr[:, 2]}
    return profiles


def load_all_cfl3d():
    """Load all CFL3D reference data for SA and SST."""
    cfl3d = {}
    for turb in ("sa", "sst"):
        uv_zones = load_cfl3d_zones(f"backstep_cfl3d_uv_{turb}.dat")
        # Auto-detect uv scaling: if peak |value| < 0.5, data is in u'v'/Uref^2 and needs *1000.
        peak = max(np.max(np.abs(d["value"])) for d in uv_zones.values())
        if peak < 0.5:
            for d in uv_zones.values():
                d["value"] = d["value"] * 1000
        cfl3d[turb] = {
            "cf_x": load_cfl3d_xy(f"backstep_cfl3d_cf_{turb}.dat"),
            "cp_x": load_cfl3d_xy(f"backstep_cfl3d_cp_{turb}.dat"),
            "vel": load_cfl3d_zones(f"backstep_cfl3d_vel_{turb}.dat"),
            "uv": uv_zones,
        }
    return cfl3d

def parse_case(case_name):
    """Return (bc_key, model_key) parsed from BackwardStep_<bc>_<model>_<solver>."""
    bc = "inflow_outflow" if "inflow_outflow" in case_name else "freestream_outflow"
    model = "sst" if re.search(r"_sst(_|$)", case_name) else "sa"
    return bc, model


def short_label(case_name):
    bc, model = parse_case(case_name)
    return f"Flow360 {MODEL_LABELS[model]} ({BC_LABELS[bc]})"


def extract_case_data(case):
    """Extract monitor DataFrames from an (already-finished) in-kernel case object."""
    monitors = case.results.monitors
    return {
        "surface_probes": monitors.get_monitor_by_name("surface_probes_bottom").as_dataframe(),
        "cp_shift_ref": monitors.get_monitor_by_name("Cp_shift_ref").as_dataframe(),
        "ref_location": monitors.get_monitor_by_name("ref_location").as_dataframe(),
        "profiles": monitors.get_monitor_by_name("Probe output").as_dataframe(),
    }


def average_last_fraction(df, fraction=0.1):
    """Average the last `fraction` of rows, excluding step columns."""
    step_cols = {"pseudo_step", "physical_step", "time", "time_units"}
    value_cols = [c for c in df.columns if c not in step_cols]
    n = max(1, int(len(df) * fraction))
    return df[value_cols].tail(n).mean()


def extract_probe_points(avg_series, prefix):
    """Extract per-probe-point data from averaged monitor columns ({prefix}_{idx}_{field})."""
    pattern = re.compile(re.escape(prefix) + r"_(\d+)_(.+)")
    points = {}
    for col in avg_series.index:
        m = pattern.match(col)
        if m:
            points.setdefault(int(m.group(1)), {})[m.group(2)] = avg_series[col]
    return points


def probe_points_to_arrays(points, x_field, value_field):
    indices = sorted(points.keys())
    x_vals, v_vals = [], []
    for idx in indices:
        p = points[idx]
        if x_field in p and value_field in p:
            x_vals.append(p[x_field])
            v_vals.append(p[value_field])
    return np.array(x_vals), np.array(v_vals)


def get_surface_data(surf_avg, field, array_prefix):
    prefix = f"surface_probes_bottom_{array_prefix}"
    points = extract_probe_points(surf_avg, prefix)
    return probe_points_to_arrays(points, "monitorLocation_x", field)


def get_profile_data(prof_avg, station_name, field):
    prefix = f"Probe output_{station_name}"
    points = extract_probe_points(prof_avg, prefix)
    return probe_points_to_arrays(points, "monitorLocation_z", field)


def get_bottom_cf(cdata):
    """Combined (x, Cf_X) along the bottom wall, before and after the step."""
    x_all, cf_all = [], []
    for array_name in ("bottom_before_step", "bottom_after_step"):
        x, cf = get_surface_data(cdata["surf_avg"], "Cf_X", array_name)
        if len(x) > 0:
            x_all.append(x)
            cf_all.append(cf)
    if not x_all:
        return np.array([]), np.array([])
    return np.concatenate(x_all), np.concatenate(cf_all)


def renormalize(data):
    ref_avg = average_last_fraction(data["ref_location"])
    shift_avg = average_last_fraction(data["cp_shift_ref"])

    u_norm_ref = ref_avg[[c for c in ref_avg.index if c.endswith("_u_norm")][0]]
    cp_ref = ref_avg[[c for c in ref_avg.index if c.endswith("_Cp")][0]]
    cp_shift_raw = shift_avg[[c for c in shift_avg.index if c.endswith("_Cp")][0]]
    cp_shift = (cp_shift_raw - cp_ref) / u_norm_ref**2

    return {
        "u_norm_ref": u_norm_ref,
        "cp_ref": cp_ref,
        "cp_shift_raw": cp_shift_raw,
        "cp_shift": cp_shift,
        "surf_avg": average_last_fraction(data["surface_probes"]),
        "prof_avg": average_last_fraction(data["profiles"]),
    }


def flow360_style(case_name):
    bc, model = parse_case(case_name)
    return {
        "color": FLOW360_MODEL_COLORS[model],
        "linestyle": BC_LINESTYLES[bc],
        "linewidth": 2.2,
        "zorder": 5,
    }


def markevery_in_window(x, xlim=CF_CP_XLIM, target=20):
    """Integer ``markevery`` giving ~``target`` markers within the visible ``xlim``."""
    x = np.asarray(x)
    n_visible = int(np.count_nonzero((x >= xlim[0]) & (x <= xlim[1])))
    return max(1, n_visible // target)

exp_cf = load_cf_exp()
exp_cp = load_cp_exp()
exp_profiles = load_profiles_exp()
cfl3d = load_all_cfl3d()
print("Loaded CFL3D reference data (SA, SST)")

all_case_data = {}
for record in cases:
    case_name = record["name"]
    print(f"Processing: {case_name}")
    data = extract_case_data(record["case"])
    all_case_data[case_name] = renormalize(data)

def plot_cf(all_case_data, exp_cf, cfl3d):
    fig, ax = plt.subplots()
    for case_name, cdata in all_case_data.items():
        u_ref = cdata["u_norm_ref"]
        x_all, cf_all = get_bottom_cf(cdata)
        if len(x_all) == 0:
            continue
        cf_renorm = cf_all / u_ref**2
        sort_idx = np.argsort(x_all)
        xs, ys = x_all[sort_idx] / STEP_HEIGHT, cf_renorm[sort_idx]
        ax.plot(xs, ys, label=short_label(case_name), **flow360_style(case_name))

    ax.errorbar(exp_cf["x_H"], exp_cf["cf"], yerr=exp_cf["error"],
                fmt="o", color=REFERENCE_COLOR, markersize=4, capsize=2, label="Exp", zorder=6)
    for turb, tdata in cfl3d.items():
        x_H, cf = tdata["cf_x"]
        ax.plot(x_H, cf, label=f"CFL3D {turb.upper()}",
                marker=CFL3D_MARKER, markersize=CFL3D_MARKERSIZE,
                markevery=markevery_in_window(x_H), **CFL3D_STYLES[turb])

    ax.set_xlabel("x / H")
    ax.set_ylabel("Cf")
    ax.set_xlim(CF_CP_XLIM)
    ax.set_title("Skin Friction Coefficient")
    ax.legend(fontsize=7, loc="best", ncol=2)
    fig.tight_layout()
    fig.savefig(os.path.join(RESULTS_DIR, "cf_vs_x.png"), dpi=200)
    return fig


plot_cf(all_case_data, exp_cf, cfl3d)
plt.show()

def plot_cp(all_case_data, exp_cp, cfl3d):
    fig, ax = plt.subplots()
    for case_name, cdata in all_case_data.items():
        u_ref = cdata["u_norm_ref"]
        cp_ref = cdata["cp_ref"]
        cp_shift = cdata["cp_shift"]
        x_all, cp_all = [], []
        for array_name in ("bottom_before_step", "bottom_after_step"):
            x, cp = get_surface_data(cdata["surf_avg"], "Cp", array_name)
            if len(x) > 0:
                x_all.append(x)
                cp_all.append(cp)
        if not x_all:
            continue
        x_all = np.concatenate(x_all)
        cp_all = np.concatenate(cp_all)
        cp_renorm = (cp_all - cp_ref) / u_ref**2 - cp_shift
        sort_idx = np.argsort(x_all)
        xs, ys = x_all[sort_idx] / STEP_HEIGHT, cp_renorm[sort_idx]
        ax.plot(xs, ys, label=short_label(case_name), **flow360_style(case_name))

    ax.plot(exp_cp["x_H"], exp_cp["cp"], "o", color=REFERENCE_COLOR, markersize=4, label="Exp", zorder=6)
    for turb, tdata in cfl3d.items():
        x_H, cp = tdata["cp_x"]
        ax.plot(x_H, cp, label=f"CFL3D {turb.upper()}",
                marker=CFL3D_MARKER, markersize=CFL3D_MARKERSIZE,
                markevery=markevery_in_window(x_H), **CFL3D_STYLES[turb])

    ax.set_xlabel("x / H")
    ax.set_ylabel("Cp")
    ax.set_xlim(CF_CP_XLIM)
    ax.set_title("Pressure Coefficient")
    ax.legend(fontsize=7, loc="best", ncol=2)
    fig.tight_layout()
    fig.savefig(os.path.join(RESULTS_DIR, "cp_vs_x.png"), dpi=200)
    return fig


plot_cp(all_case_data, exp_cp, cfl3d)
plt.show()

def _profile_legend(ax):
    """Stack two legends upper-left: station colors on top, source/BC styles below."""
    from matplotlib.lines import Line2D
    station_handles = [
        Line2D([0], [0], color=s["color"], lw=2, label=f"x/H = {xH}")
        for xH, s in PROFILE_STATIONS.items()
    ]
    source_handles = [
        Line2D([0], [0], color="0.35", label=label, **PROFILE_SERIES_STYLES[key])
        for key, label in (("exp", "Experiment"),
                           ("cfl3d", "CFL3D"),
                           ("freestream", "Flow360 - Freestream"),
                           ("inflow", "Flow360 - Subsonic Inflow"))
    ]
    leg_stations = ax.legend(handles=station_handles, loc="upper left",
                             fontsize=8, framealpha=0.9, title="Station")
    ax.add_artist(leg_stations)
    ax.figure.canvas.draw()
    bbox = leg_stations.get_window_extent().transformed(ax.transAxes.inverted())
    ax.legend(handles=source_handles, loc="upper left",
              bbox_to_anchor=(bbox.x0, bbox.y0), fontsize=8, framealpha=0.9)


def _plot_profile_set(all_case_data, exp_profiles, cfl3d, *, field, exp_key, cfl3d_key,
                      renorm, xlabel, xlim, title, fname_base):
    """Produce one profile figure per turbulence model (SA, SST) for a quantity."""
    figs = []
    for turb in ("sa", "sst"):
        fig, ax = plt.subplots(figsize=(7, 8))
        cases_by_bc = {}
        for case_name, cdata in all_case_data.items():
            bc, model = parse_case(case_name)
            if model == turb:
                cases_by_bc[bc] = (case_name, cdata)

        for xH, sinfo in PROFILE_STATIONS.items():
            color = sinfo["color"]
            if xH in exp_profiles and exp_key in exp_profiles[xH]:
                ax.plot(exp_profiles[xH][exp_key], exp_profiles[xH]["y"],
                        color=color, **PROFILE_SERIES_STYLES["exp"])
            ref = cfl3d[turb][cfl3d_key]
            if xH in ref:
                d = ref[xH]
                ax.plot(d["value"], d["y_H"], color=color, **PROFILE_SERIES_STYLES["cfl3d"])
            for bc_key, style_key in (("freestream_outflow", "freestream"),
                                      ("inflow_outflow", "inflow")):
                if bc_key not in cases_by_bc:
                    continue
                case_name, cdata = cases_by_bc[bc_key]
                z, val = get_profile_data(cdata["prof_avg"], sinfo["monitor"], field)
                if len(z) == 0:
                    continue
                order = np.argsort(z)
                zs = z[order] / STEP_HEIGHT
                vs = renorm(val[order], cdata)
                ax.plot(vs, zs, color=color, **PROFILE_SERIES_STYLES[style_key])

        ax.set_xlabel(xlabel)
        ax.set_ylabel("y / H")
        ax.set_xlim(xlim)
        ax.set_ylim((0, 3))
        ax.set_title(f"{title} - {MODEL_TITLES[turb]} turbulence model")
        fig.tight_layout()
        _profile_legend(ax)
        fig.savefig(os.path.join(RESULTS_DIR, f"{fname_base}_{turb}.png"), dpi=200)
        figs.append(fig)
    return figs

_plot_profile_set(
    all_case_data, exp_profiles, cfl3d,
    field="u_norm", exp_key="u", cfl3d_key="vel",
    renorm=lambda val, cdata: val / cdata["u_norm_ref"],
    xlabel="u / Uref", xlim=(-0.4, 1.2),
    title="Axial velocity profiles", fname_base="velocity_profiles",
)
plt.show()

_plot_profile_set(
    all_case_data, exp_profiles, cfl3d,
    field="uv_norm", exp_key="uv", cfl3d_key="uv",
    renorm=lambda val, cdata: val / cdata["u_norm_ref"]**2,
    xlabel="1000 u'v' / Uref^2", xlim=(-15, 5),
    title="Turbulent shear stress profiles", fname_base="shear_stress_profiles",
)
plt.show()

def reattachment_xH(x, cf, x_min=0.05):
    """First downstream (x>x_min) sign change of Cf from negative to positive, x/H."""
    order = np.argsort(x)
    x, cf = x[order], cf[order]
    mask = x > x_min
    x, cf = x[mask], cf[mask]
    for i in range(len(cf) - 1):
        if cf[i] < 0 and cf[i + 1] >= 0:
            return (x[i] - cf[i] * (x[i + 1] - x[i]) / (cf[i + 1] - cf[i])) / STEP_HEIGHT
    return float("nan")


def compute_reattachment_table(all_case_data, cfl3d):
    cfl3d_reattach = {
        turb: round(reattachment_xH(cfl3d[turb]["cf_x"][0] * STEP_HEIGHT, cfl3d[turb]["cf_x"][1]), 3)
        for turb in ("sa", "sst")
    }
    rows = []
    for case_name, cdata in all_case_data.items():
        bc, model = parse_case(case_name)
        x, cf = get_bottom_cf(cdata)
        xr = round(reattachment_xH(x, cf), 3)
        rows.append({
            "source": "Flow360",
            "turbulence_model": MODEL_LABELS[model],
            "inlet_bc": BC_LABELS[bc],
            "reattachment_xH": xr,
            "diff_vs_cfl3d_xH": round(abs(xr - cfl3d_reattach[model]), 3),
        })
    for turb in ("sa", "sst"):
        rows.append({
            "source": "CFL3D",
            "turbulence_model": MODEL_LABELS[turb],
            "inlet_bc": "-",
            "reattachment_xH": cfl3d_reattach[turb],
            "diff_vs_cfl3d_xH": "-",
        })
    rows.append({
        "source": EXP_LABEL,
        "turbulence_model": "-",
        "inlet_bc": "-",
        "reattachment_xH": EXP_REATTACH_XH,
        "diff_vs_cfl3d_xH": "-",
    })
    return rows


reattach_rows = compute_reattachment_table(all_case_data, cfl3d)

fieldnames = ["source", "turbulence_model", "inlet_bc", "reattachment_xH", "diff_vs_cfl3d_xH"]
with open(os.path.join(RESULTS_DIR, "reattachment.csv"), "w", newline="") as f:
    writer = csv.DictWriter(f, fieldnames=fieldnames)
    writer.writeheader()
    writer.writerows(reattach_rows)

print(f"{'Source':<32} {'Model':<6} {'Inlet BC':<12} {'x_r/H':<8}")
for r in reattach_rows:
    print(f"{r['source']:<32} {r['turbulence_model']:<6} {r['inlet_bc']:<12} {r['reattachment_xH']}")

Backward-Facing Step: Separation and Reattachment¶

Requirements¶

Imports¶

Input data¶

Load project¶

Physics setup¶

Outputs setup¶

Simulation Params¶

Submit cases¶

Wait for completion¶

Postprocessing¶

Setup¶

Reference data loaders¶

Extract and renormalize the case results¶

Skin-friction coefficient along the bottom wall¶

Pressure coefficient along the bottom wall¶

Downstream profiles¶

Axial velocity profiles (SA and k-ω SSTm)¶

Reynolds shear-stress profiles (SA and k-ω SSTm)¶

Reattachment length table¶