CFDManager/docs/theory/solver_2x_sdf/run_factors.py

# ФАКТОРНОЕ ИССЛЕДОВАНИЕ постановки: какие границы расчётной области отвечают за завышение
# Cd / rms Cl относительно безграничной литературы (Re=150: Cd≈1.33, rms Cl≈0.3, St≈0.183).
#
# Контекст (эксперимент №0 — серия по боковой блокировке, см. ноутбук, раздел 23): Cd(β)
# немонотонен (1.80 → 1.77 → 1.87 при β = 0.178 → 0.125 → 0.089) → боковое стеснение НЕ
# доминирующий фактор. Подозреваемые — продольные границы, которые серия держала константой:
#   вход 2.5D (жёсткий равномерный профиль), выход 8.3D (Zou-He давление с жёстким uy=0,
#   отражает вихри дорожки), и их связка.
#
# Здесь каждый фактор варьируется ОТДЕЛЬНО при фиксированной боковой блокировке (Ny=90):
# изоляция вклада входа, выхода и выходного ГУ. Все результаты накапливаются в out/factors.csv
# (одна строка на прогон, с полной конфигурацией) — реестр для документирования исследования.
#
# Запуск: python run_factors.py
#   AMR_FACTORS=B0_baseline,F4_outlet_uy python run_factors.py   # подмножество факторов
#   AMR_STEPS=20000 python run_factors.py                        # смоук
import os, csv, time
import numpy as np
import config as C
import backend as B
import solver as S
import diagnostics as Diag
import visualization as Viz

HERE = os.path.dirname(os.path.abspath(__file__))
OUT = os.path.join(HERE, "out"); os.makedirs(OUT, exist_ok=True)
LIT = dict(Cd=1.33, Clrms=0.3, St=0.183)

# (имя, overrides конфига, что изолирует)
# Эксперимент №1 (B0–F5, прогнан 2026-06-10): итоги — мягкий uy-выход (F4) выигрывает по всем
# метрикам; ВСЕ длинные домены без демпфера уходят в смещённый режим ⟨ρ⟩≈1.66 (акустический
# резонатор «вход-скорость/выход-давление» + накачка массы входом) или взрываются (F5).
# Эксперимент №2 (F6–F8): губка перед выходом — лечит акустику, открывает длинные домены.
FACTORS = [
    ("B0_baseline",   dict(),
     "база 174x90: вход 2.5D, выход 8.3D, uy=0 (контроль; кросс-чек уже с конвективным членом)"),
    ("F1_outlet_far", dict(Nx=300),
     "выход 8.3D -> 16.2D (вход 2.5D как в базе) — изолирует близость выхода"),
    ("F2_inlet_far",  dict(Nx=230, cx=96),
     "вход 2.5D -> 6.0D (выход 8.3D как в базе) — изолирует близость входа"),
    ("F3_inout_far",  dict(Nx=390, cx=96),
     "вход 6.0D + выход 18.3D — оба продольных фактора вместе"),
    ("F4_outlet_uy",  dict(outlet_uy="extrapolate"),
     "мягкий выход: uy нуль-градиент вместо uy=0 (геометрия базовая) — изолирует отражение вихрей"),
    ("F5_combo",      dict(Nx=390, cx=96, outlet_uy="extrapolate"),
     "вход 6.0D + выход 18.3D + мягкий uy-выход — всё вместе"),
    # --- эксперимент №2: губка (absorbing layer) ---
    ("F6_sponge",     dict(outlet_uy="extrapolate", sponge_len=32),
     "база 174x90 + мягкий uy + ГУБКА 32 столбца (nu x30) — валидация губки против F4"),
    ("F7_long_sponge", dict(Nx=390, cx=96, outlet_uy="extrapolate", sponge_len=32),
     "вход 6D + выход 18.3D + мягкий uy + губка — целевая «чистая» постановка"),
    ("F8_long_sp_uy0", dict(Nx=390, cx=96, sponge_len=32),
     "длинный домен + губка, но ЖЁСТКОЕ uy=0 — изолирует вклад губки от uy-режима"),
]

CSV_FIELDS = ["when", "name", "status", "steps", "Nx", "Ny", "cx", "in_D", "out_D", "outlet_uy",
              "sponge_len", "beta", "St", "St_corr", "Cd", "Clrms", "Cm", "Cd_st", "Clrms_st",
              "dCd_pct", "dCl_pct", "rho_last", "cd_win_std", "cl_win_std", "wall_s", "desc"]

def window_stats(res, cfg, nwin=10):
    """Срединная изменчивость: std оконных средних Cd и оконных rms Cl (окна 2..nwin,
    первое окно с разгоном пропускается). Метрика низкочастотной модуляции дорожки."""
    Fx = B.to_cpu(res["Fx"]); Fy = B.to_cpu(res["Fy"])
    Cd = 2.0*Fx/(cfg.U**2*res["DL"]); Cl = 2.0*Fy/(cfg.U**2*res["DL"])
    n = len(Cd); w = n // nwin
    cds = [float(Cd[k*w:(k+1)*w].mean()) for k in range(1, nwin)]
    cls = [float(Cl[k*w:(k+1)*w].std()) for k in range(1, nwin)]
    return float(np.std(cds)), float(np.std(cls))

def run_factor(name, overrides, desc):
    cfg = C.from_env(frame_every=0, **overrides)
    in_D = cfg.cx / cfg.D; out_D = (cfg.Nx - 1 - cfg.cx) / cfg.D
    print(f"\n=== Фактор {name}: {desc}")
    print(f"    домен {cfg.Nx}x{cfg.Ny}; вход {in_D:.1f}D / выход {out_D:.1f}D; "
          f"outlet_uy={cfg.outlet_uy}; губка={cfg.sponge_len or '—'}; "
          f"патч x=[{cfg.ax1},{cfg.bx1}] y=[{cfg.ay1},{cfg.by1}]; steps={cfg.steps}")
    t0 = time.perf_counter()
    sim = S.Solver(cfg)
    res = sim.run(Viz.make_progress(name[:8]))
    wall = time.perf_counter() - t0
    print("  время:", round(wall, 1), "s")
    status = "ok" if len(res["Fx"]) == cfg.steps + 1 else f"BLEW_UP@{len(res['Fx'])}"
    a = Diag.analyze(res, cfg)
    Diag.print_table(a, cfg)
    Diag.print_crosscheck(a)
    Diag.convergence_report(res, cfg)
    cd_ws, cl_ws = window_stats(res, cfg)
    rho_last = float(B.to_cpu(res["rho"])[-1])
    np.savez(os.path.join(OUT, f"factor_{name}.npz"),
             St=a["St"], Cd=a["Cd"], Clrms=a["Clrms"], uyrms=a["uyrms"],
             Cm=a.get("Cm", float("nan")),
             Cd_st=a.get("Cd_st", float("nan")), Clrms_st=a.get("Clrms_st", float("nan")),
             Nx=cfg.Nx, Ny=cfg.Ny, cx=cfg.cx, D=cfg.D, U=cfg.U, Re=cfg.Re, steps=cfg.steps,
             outlet_uy=cfg.outlet_uy, beta=cfg.D/cfg.Ny,
             Cl=a["Cl_s"], uy=a["uy"], rho=B.to_cpu(res["rho"]),
             Fx=B.to_cpu(res["Fx"]), Fy=B.to_cpu(res["Fy"]), Tz=B.to_cpu(res["Tz"]))
    beta = cfg.D / cfg.Ny
    dcd = abs(a.get("Cd_st", float("nan")) - a["Cd"]) / a["Cd"] * 100.0
    dcl = abs(a.get("Clrms_st", float("nan")) - a["Clrms"]) / max(a["Clrms"], 1e-12) * 100.0
    row = dict(when=time.strftime("%Y-%m-%d %H:%M:%S"), name=name, status=status, steps=cfg.steps,
               Nx=cfg.Nx, Ny=cfg.Ny, cx=cfg.cx, in_D=round(in_D, 2), out_D=round(out_D, 2),
               outlet_uy=cfg.outlet_uy, sponge_len=cfg.sponge_len, beta=round(beta, 4),
               St=round(a["St"], 4), St_corr=round(a["St"]*(1-beta), 4),
               Cd=round(a["Cd"], 4), Clrms=round(a["Clrms"], 4),
               Cm=round(a.get("Cm", float("nan")), 6),
               Cd_st=round(a.get("Cd_st", float("nan")), 4),
               Clrms_st=round(a.get("Clrms_st", float("nan")), 4),
               dCd_pct=round(dcd, 2), dCl_pct=round(dcl, 2),
               rho_last=round(rho_last, 5), cd_win_std=round(cd_ws, 4),
               cl_win_std=round(cl_ws, 4), wall_s=round(wall, 1), desc=desc)
    del sim, res, a
    B.cp.get_default_memory_pool().free_all_blocks()
    return row

def append_csv(rows):
    path = os.path.join(OUT, "factors.csv")
    if os.path.exists(path):
        with open(path, encoding="utf-8") as fh:
            header = fh.readline().strip()
        if header != ",".join(CSV_FIELDS):
            # схема реестра изменилась — старый файл сохраняем (история не теряется)
            k = 1
            while os.path.exists(os.path.join(OUT, f"factors_v{k}.csv")):
                k += 1
            os.rename(path, os.path.join(OUT, f"factors_v{k}.csv"))
            print(f"[реестр] схема CSV изменилась — прежний файл сохранён как factors_v{k}.csv")
    new = not os.path.exists(path)
    with open(path, "a", newline="", encoding="utf-8") as fh:
        wr = csv.DictWriter(fh, fieldnames=CSV_FIELDS)
        if new:
            wr.writeheader()
        for r in rows:
            wr.writerow(r)
    print("CSV:", path, f"(+{len(rows)} строк)")

def save_summary_figure(rows, path):
    import matplotlib; matplotlib.use("Agg")
    import matplotlib.pyplot as plt
    names = [r["name"] for r in rows]
    xs = np.arange(len(rows))
    fig, axs = plt.subplots(1, 3, figsize=(15, 4.6))
    ax = axs[0]
    ax.bar(xs, [r["Cd"] for r in rows], color="#1f6feb", label="⟨Cd⟩ GMEM")
    ax.plot(xs, [r["Cd_st"] for r in rows], "o", color="#ff8f00", label="⟨Cd⟩ ∮[σ·n−ρu(u·n)]ds")
    ax.axhline(LIT["Cd"], color="#d1495b", ls="--", label=f"лит. безгранич. {LIT['Cd']}")
    ax.set_ylabel("⟨Cd⟩"); ax.set_title("Сопротивление по факторам")
    ax = axs[1]
    ax.bar(xs, [r["Clrms"] for r in rows], color="#2e7d32")
    ax.axhline(LIT["Clrms"], color="#d1495b", ls="--", label=f"лит. ~{LIT['Clrms']}")
    ax.set_ylabel("rms Cl"); ax.set_title("Амплитуда подъёмной силы")
    ax = axs[2]
    ax.bar(xs, [r["cd_win_std"] for r in rows], color="#7b1fa2")
    ax.set_ylabel("std оконных ⟨Cd⟩"); ax.set_title("НЧ-модуляция (меньше — стабильнее)")
    for ax in axs:
        ax.set_xticks(xs)
        ax.set_xticklabels([f"{r['name']}\nвх {r['in_D']}D / вых {r['out_D']}D\n"
                            f"uy:{r['outlet_uy']} губ:{r['sponge_len']}" for r in rows], fontsize=7)
        ax.legend(fontsize=8)
    fig.suptitle("Факторное исследование постановки (Ny=90, β=0.178 фикс., 100k шагов)", fontsize=12)
    fig.tight_layout(rect=(0, 0, 1, 0.94))
    fig.savefig(path, bbox_inches="tight"); plt.close(fig)
    print("saved", os.path.basename(path))

def main():
    sel = os.environ.get("AMR_FACTORS")
    todo = FACTORS if not sel else [f for f in FACTORS if f[0] in sel.split(",")]
    print(f"[факторы] backend={B.BACKEND}; прогонов: {len(todo)}")
    rows = [run_factor(*f) for f in todo]
    print("\n=== Сводка факторного исследования ===")
    hdr = (f"{'фактор':15}{'вх,D':>6}{'вых,D':>7}{'uy':>12}{'губка':>6}{'St':>7}{'<Cd>':>8}"
           f"{'rms Cl':>8}{'Cd σ·n':>8}{'dCd%':>6}{'<ρ>фин':>8}{'модул.':>8}  статус")
    print(hdr)
    for r in rows:
        print(f"{r['name']:15}{r['in_D']:>6}{r['out_D']:>7}{r['outlet_uy']:>12}{r['sponge_len']:>6}"
              f"{r['St']:>7.3f}{r['Cd']:>8.3f}{r['Clrms']:>8.4f}{r['Cd_st']:>8.3f}{r['dCd_pct']:>6.1f}"
              f"{r['rho_last']:>8.4f}{r['cd_win_std']:>8.4f}  {r['status']}")
    print(f"{'лит. безгранич.':15}{'':>6}{'':>7}{'':>12}{'':>6}{LIT['St']:>7.3f}{LIT['Cd']:>8.3f}{LIT['Clrms']:>8.4f}")
    append_csv(rows)
    if len(rows) > 1:
        save_summary_figure(rows, os.path.join(OUT, "factors_summary.png"))
    print("DONE (факторы)")

if __name__ == "__main__":
    main()