# Общее GPU-ядро для основных AMR-решателей (D2Q9 KBC) на цилиндре.
# Математика ДОСЛОВНО та же, что в проверенных CPU-прототипах (_amr_anims*.py):
#   - энтропийное равновесие (product-form), KBC-N1 столкновение f<-f-β(2Δs+γΔh),
#   - AMR-связка τ_f=2τ_c-½, R_cf=τ_f/(2τ_c), m подшагов, временна́я интерполяция,
#   - сила по Mei (обмен импульсом), вход/выход/стенки, опционально SDF+Bouzidi.
# Оптимизация под GPU (CuPy): всё векторно на device; проектор через matmul; ГУ и силы
# через where (без gather/scatter); БЕЗ пер-шаговых синхронизаций host↔device.
import os, io as _io, numpy as np

# ---- backend: CuPy (GPU) с откатом на numpy (CPU) ----
try:
    import cupy as cp
    _ = (cp.zeros(2) + 1).sum()      # реальная операция на device — ловит "cupy есть, GPU нет"
    xp = cp; GPU = True
except Exception:
    import numpy as cp               # noqa
    xp = np; GPU = False
def to_cpu(a):
    return cp.asnumpy(a) if GPU else np.asarray(a)
DTYPE = xp.float64 if os.environ.get("AMR_FP64") else xp.float32
GREL  = 1e-8    # ОТНОСИТЕЛЬНЫЙ порог вырожденности знаменателя γ (доля от ‖Δ‖²); см. solver_2x_sdf/backend.py.
                # Прежний абсолютный порог (1e-6 в fp32) срабатывал на большинстве узлов и подменял KBC на LBGK.
BACKEND = f"{'CuPy/GPU' if GPU else 'numpy/CPU'} dtype={'float64' if DTYPE == xp.float64 else 'float32'}"

# ---- решётка D2Q9 ----
Cx = [0, 1, 0, -1, 0, 1, -1, -1, 1]          # python-инты для roll/индексации
Cy = [0, 0, 1, 0, -1, 1, 1, -1, -1]
OPP = [0, 3, 4, 1, 2, 7, 8, 5, 6]
Q = 9
_W = np.array([4/9, 1/9, 1/9, 1/9, 1/9, 1/36, 1/36, 1/36, 1/36])
_CXn = np.array(Cx, float); _CYn = np.array(Cy, float)
# моментный базис и проектор на сдвиг {cx²−cy², cxcy} (KBC-N1) — считаем в numpy, грузим на device
_M = np.zeros((Q, Q)); _M[0] = 1; _M[1] = _CXn; _M[2] = _CYn; _M[3] = 3*(_CXn**2+_CYn**2)-2
_M[4] = _CXn**2-_CYn**2; _M[5] = _CXn*_CYn; _M[6] = _CXn**2*_CYn; _M[7] = _CXn*_CYn**2; _M[8] = _CXn**2*_CYn**2
_Dm = np.zeros((Q, Q)); _Dm[4, 4] = _Dm[5, 5] = 1.0
_Psn = np.linalg.inv(_M) @ _Dm @ _M
# device-константы
Wa = xp.asarray(_W, dtype=DTYPE)
CXa = xp.asarray(_CXn, dtype=DTYPE); CYa = xp.asarray(_CYn, dtype=DTYPE)
Ps = xp.asarray(_Psn, dtype=DTYPE)
CS2 = 1.0/3.0

# ---- ядро LBM/KBC (идентичная математика) ----
def feq(rho, u):
    ux = xp.clip(u[0], -0.95, 0.95); uy = xp.clip(u[1], -0.95, 0.95)
    sx = xp.sqrt(1 + 3*ux*ux); sy = xp.sqrt(1 + 3*uy*uy); base = rho*(2-sx)*(2-sy)
    Bx = ((2*ux+sx)/(1-ux))[None]**CXa[:, None, None]
    By = ((2*uy+sy)/(1-uy))[None]**CYa[:, None, None]
    return Wa[:, None, None]*base[None]*Bx*By
def macros(f):
    r = f.sum(0)
    return r, xp.stack([(CXa[:, None, None]*f).sum(0)/r, (CYa[:, None, None]*f).sum(0)/r])
def collide(f, fe, beta):
    dfn = f - fe
    ds = (Ps @ dfn.reshape(Q, -1)).reshape(dfn.shape)   # проектор на сдвиг (matmul, GPU-friendly)
    dh = dfn - ds; inv = 1.0/fe
    num = (ds*dh*inv).sum(0); den = (dh*dh*inv).sum(0)
    nrm = (dfn*dfn*inv).sum(0)                                # ‖Δ‖² — масштаб неравновесия узла
    ok = den > GREL*nrm
    den_safe = xp.where(ok, den, xp.asarray(1.0, DTYPE))      # избегаем 0/0 (обе ветки where считаются)
    g = xp.where(ok, 1/beta - (2 - 1/beta)*num/den_safe, xp.asarray(2.0, DTYPE))
    return f - beta*(2*ds + g[None]*dh)
def stream(f):
    fs = xp.empty_like(f)
    for i in range(Q):
        fs[i] = xp.roll(f[i], (Cy[i], Cx[i]), (0, 1))
    return fs

# ---- AMR-связка ----
def patch(pNx, pNy, ax, bx, ay, by, r):
    Wx, Wy = bx-ax, by-ay; Nfx, Nfy = r*Wx+1, r*Wy+1
    FX, FY = np.meshgrid(ax+np.arange(Nfx)/r, ay+np.arange(Nfy)/r)
    x0 = np.floor(FX).astype(np.int64); y0 = np.floor(FY).astype(np.int64)
    x1 = np.minimum(x0+1, pNx-1); y1 = np.minimum(y0+1, pNy-1)
    return dict(Nfx=Nfx, Nfy=Nfy, ax=ax, bx=bx, ay=ay, by=by, r=r,
                x0=xp.asarray(x0), y0=xp.asarray(y0), x1=xp.asarray(x1), y1=xp.asarray(y1),
                tx=xp.asarray(FX-x0, DTYPE), ty=xp.asarray(FY-y0, DTYPE),
                slx=slice(r, r*Wx, r), sly=slice(r, r*Wy, r))
def pint(fld, P):
    return (fld[..., P["y0"], P["x0"]]*(1-P["tx"])*(1-P["ty"]) + fld[..., P["y0"], P["x1"]]*P["tx"]*(1-P["ty"])
            + fld[..., P["y1"], P["x0"]]*(1-P["tx"])*P["ty"] + fld[..., P["y1"], P["x1"]]*P["tx"]*P["ty"])
def ghost(pf, P, Rcf):
    r, u = macros(pf); neq = pf - feq(r, u)
    return feq(pint(r, P), xp.stack([pint(u[0], P), pint(u[1], P)])) + Rcf*pint(neq, P)
def fill(cf, gh):
    cf[:, 0, :] = gh[:, 0, :]; cf[:, -1, :] = gh[:, -1, :]; cf[:, :, 0] = gh[:, :, 0]; cf[:, :, -1] = gh[:, :, -1]
    return cf
def restrict(cf, pf, P, Rfc, fluid):
    r, u = macros(cf); neq = cf - feq(r, u); slx, sly = P["slx"], P["sly"]
    nv = feq(r[sly, slx], u[:, sly, slx]) + Rfc*neq[:, sly, slx]
    cur = pf[:, P["ay"]+1:P["by"], P["ax"]+1:P["bx"]]
    pf[:, P["ay"]+1:P["by"], P["ax"]+1:P["bx"]] = xp.where(fluid[None], nv, cur)
    return pf

# ---- граница: (1) узловой bounce-back (без SDF); (2) SDF+Bouzidi ----
def bb_set(post, pre, solid):           # узлы тела <- развёрнутые до-столкновения (как в CPU-версии)
    for i in range(Q):
        post[i] = xp.where(solid, pre[OPP[i]], post[i])
    return post
def force_on(fpost, cyl, fluid):         # Mei для узловой границы: F=Σ c_i(f_i+f_ī) по линкам жидк.->тело
    Fx = 0.0; Fy = 0.0
    for i in range(1, Q):
        nb = xp.roll(cyl, (-Cy[i], -Cx[i]), (0, 1)); link = fluid & nb
        s = xp.where(link, fpost[i] + fpost[OPP[i]], xp.asarray(0.0, DTYPE)).sum()
        Fx = Fx + Cx[i]*s; Fy = Fy + Cy[i]*s
    return Fx, Fy
def build_bc(solid_np, phi_np, cyl_np):  # SDF-границу строим на host один раз, переносим на device
    fluid = ~solid_np; bc = []
    for i in range(1, Q):
        nb_solid = np.roll(solid_np, (-Cy[i], -Cx[i]), (0, 1)); mask = fluid & nb_solid
        cl = mask & np.roll(cyl_np, (-Cy[i], -Cx[i]), (0, 1))
        phinb = np.roll(phi_np, (-Cy[i], -Cx[i]), (0, 1))
        with np.errstate(divide="ignore", invalid="ignore"):
            qq = phi_np/(phi_np - phinb)
        q = np.where(mask, np.clip(qq, 0.02, 0.98), 0.5)
        xff = mask & (~np.roll(solid_np, (Cy[i], Cx[i]), (0, 1)))
        bc.append((i, OPP[i], xp.asarray(mask), xp.asarray(q, DTYPE), xp.asarray(xff), xp.asarray(cl)))
    return bc
def apply_bc(f, fpost, bc):              # Bouzidi (линейный), полностью через where (GPU)
    for (i, ib, mask, q, xff, cl) in bc:
        fi = fpost[i]; fib = fpost[ib]; fiback = xp.roll(fpost[i], (Cy[i], Cx[i]), (0, 1))
        near = mask & (q < 0.5) & xff; bad = mask & (q < 0.5) & (~xff); far = mask & (q >= 0.5)
        new = f[ib]
        new = xp.where(near, 2*q*fi + (1-2*q)*fiback, new)
        new = xp.where(far, (1/(2*q))*fi + (1-1/(2*q))*fib, new)
        new = xp.where(bad, fi, new)
        f[ib] = new
    return f
def force_bc(fpost, f, bc):              # Mei для интерполированной границы
    Fx = 0.0; Fy = 0.0
    for (i, ib, mask, q, xff, cl) in bc:
        s = xp.where(cl, fpost[i] + f[ib], xp.asarray(0.0, DTYPE)).sum()
        Fx = Fx + Cx[i]*s; Fy = Fy + Cy[i]*s
    return Fx, Fy

# ---- геометрия (ИДЕНТИЧНА CPU-версиям) ----
Nx, Ny, D, cx, cy = 150, 72, 16, 40, 36
U, Re, STEPS, ramp, FE = 0.07, 150.0, int(os.environ.get("AMR_STEPS", 8000)), 1000, 34
D0 = D; nu = U*D/Re; tau0 = nu/CS2 + 0.5
ax1, bx1, ay1, by1 = 16, 118, 8, 64
cx2a, cx2b, cy2a, cy2b = 26, 64, 18, 54
px, py = cx + 3*D, cy
def _cmask(NX, NY, ccx, ccy, R):
    Y, X = np.meshgrid(np.arange(NY), np.arange(NX), indexing="ij"); return (X-ccx)**2 + (Y-ccy)**2 <= R*R
def _csdf(NX, NY, ccx, ccy, R):
    Y, X = np.meshgrid(np.arange(NY), np.arange(NX), indexing="ij"); return np.sqrt((X-ccx)**2.0 + (Y-ccy)**2.0) - R
# host-маски/SDF
_walls = np.zeros((Ny, Nx), bool); _walls[0, :] = True; _walls[-1, :] = True
cyl0_np = _cmask(Nx, Ny, cx, cy, D/2); solid0_np = cyl0_np | _walls
_Yg = np.arange(Ny)[:, None]*np.ones((1, Nx))
phi0_np = np.minimum(_csdf(Nx, Ny, cx, cy, D/2), np.minimum(_Yg-0.5, (Ny-1.5)-_Yg))
P1 = patch(Nx, Ny, ax1, bx1, ay1, by1, 2); P2 = patch(P1["Nfx"], P1["Nfy"], (cx2a-ax1)*2, (cx2b-ax1)*2, (cy2a-ay1)*2, (cy2b-ay1)*2, 2)
solid1_np = _cmask(P1["Nfx"], P1["Nfy"], (cx-ax1)*2, (cy-ay1)*2, D); phi1_np = _csdf(P1["Nfx"], P1["Nfy"], (cx-ax1)*2, (cy-ay1)*2, D)
solid2_np = _cmask(P2["Nfx"], P2["Nfy"], (cx-cx2a)*4, (cy-cy2a)*4, 2*D); phi2_np = _csdf(P2["Nfx"], P2["Nfy"], (cx-cx2a)*4, (cy-cy2a)*4, 2*D)
# device-маски (без SDF)
g_solid0 = xp.asarray(solid0_np); g_cyl0 = xp.asarray(cyl0_np); g_fl0 = ~g_solid0
g_solid1 = xp.asarray(solid1_np); g_fl1 = ~g_solid1
g_solid2 = xp.asarray(solid2_np); g_fl2 = ~g_solid2
# SDF-ГУ
BC0 = build_bc(solid0_np, phi0_np, cyl0_np); BC1 = build_bc(solid1_np, phi1_np, solid1_np); BC2 = build_bc(solid2_np, phi2_np, solid2_np)
# рестрикция: жидкие узлы
fl1 = xp.asarray(~solid0_np[ay1+1:by1, ax1+1:bx1])
fl2 = xp.ones((P2["by"]-P2["ay"]-1, P2["bx"]-P2["ax"]-1), bool); fl2 = xp.asarray(fl2)
# τ / масштабы
t1 = 2*tau0 - 0.5; t2 = 2*t1 - 0.5
b0 = 1/(2*tau0); b1 = 1/(2*t1); b2 = 1/(2*t2)
R01 = t1/(2*tau0); Rf01 = 1/R01; R12 = t2/(2*t1); Rf12 = 1/R12

def smoothstep(x):
    x = min(max(x, 0.0), 1.0); return x*x*(3-2*x)

def run(mode, use_sdf, progress=None):
    """mode: none|amr2x|nested. Возвращает {Fx,Fy,uy (на host), frames (host), DL, mode}.
    progress(t, STEPS) — необязательный колбэк прогресса (вызывается каждый шаг)."""
    use1 = mode in ("amr2x", "nested"); use2 = (mode == "nested")
    DL = {"none": D0, "amr2x": 2*D0, "nested": 4*D0}[mode]
    f0 = feq(xp.ones((Ny, Nx), DTYPE), xp.zeros((2, Ny, Nx), DTYPE))
    f1 = feq(xp.ones((P1["Nfy"], P1["Nfx"]), DTYPE), xp.zeros((2, P1["Nfy"], P1["Nfx"]), DTYPE)) if use1 else None
    f2 = feq(xp.ones((P2["Nfy"], P2["Nfx"]), DTYPE), xp.zeros((2, P2["Nfy"], P2["Nfx"]), DTYPE)) if use2 else None
    Fx = xp.zeros(STEPS+1, DTYPE); Fy = xp.zeros(STEPS+1, DTYPE); uy = xp.zeros(STEPS+1, DTYPE)
    onecol = xp.ones((Ny, 1), DTYPE); frames = []
    for t in range(STEPS+1):
        rho, u0 = macros(f0)
        if t % 500 == 0 and not bool(xp.isfinite(rho).all()):
            print("BLEW UP", mode, t); Fx = Fx[:t]; Fy = Fy[:t]; uy = uy[:t]; break
        pre = f0.copy(); post0 = collide(f0, feq(rho, u0), b0)
        if not use_sdf:
            post0b = bb_set(post0.copy(), pre, g_solid0); f0 = stream(post0b)
        else:
            f0 = stream(post0); f0 = apply_bc(f0, post0, BC0)
        if mode == "none":
            Fx[t], Fy[t] = (force_bc(post0, f0, BC0) if use_sdf else force_on(post0, g_cyl0, g_fl0))
        uin = xp.zeros((2, Ny, 1), DTYPE); uin[0, :, 0] = U*smoothstep(t/ramp)
        f0[:, 1:-1, 0] = feq(onecol, uin)[:, 1:-1, 0]; f0[:, 1:-1, -1] = f0[:, 1:-1, -2]
        if use1:
            g1o = ghost(pre, P1, R01); g1n = ghost(f0, P1, R01)
            for s1 in range(2):
                pre1 = f1.copy(); r1, u1 = macros(f1); post1 = collide(f1, feq(r1, u1), b1)
                if not use_sdf:
                    f1 = stream(bb_set(post1.copy(), pre1, g_solid1))
                else:
                    f1 = stream(post1); f1 = apply_bc(f1, post1, BC1)
                if mode == "amr2x" and s1 == 1:
                    Fx[t], Fy[t] = (force_bc(post1, f1, BC1) if use_sdf else force_on(post1, g_solid1, g_fl1))
                f1 = fill(f1, (1-(s1+1)/2)*g1o + ((s1+1)/2)*g1n)
                if use2:
                    g2o = ghost(pre1, P2, R12); g2n = ghost(f1, P2, R12)
                    for s2 in range(2):
                        pre2 = f2.copy(); r2, u2 = macros(f2); post2 = collide(f2, feq(r2, u2), b2)
                        if not use_sdf:
                            f2 = stream(bb_set(post2.copy(), pre2, g_solid2))
                        else:
                            f2 = stream(post2); f2 = apply_bc(f2, post2, BC2)
                        if s1 == 1 and s2 == 1:
                            Fx[t], Fy[t] = (force_bc(post2, f2, BC2) if use_sdf else force_on(post2, g_solid2, g_fl2))
                        f2 = fill(f2, (1-(s2+1)/2)*g2o + ((s2+1)/2)*g2n)
                    f1 = restrict(f2, f1, P2, Rf12, fl2)
            f0 = restrict(f1, f0, P1, Rf01, fl1)
        col = f0[:, py, px]; uy[t] = (CYa*col).sum()/col.sum()           # зонд следа (без full-macros)
        if progress is not None: progress(t, STEPS)
        if use1 and t % FE == 0:
            frames.append((t, to_cpu(macros(f0)[1]), to_cpu(macros(f1)[1]),
                           (to_cpu(macros(f2)[1]) if use2 else None)))
    return dict(mode=mode, DL=DL, Fx=to_cpu(Fx), Fy=to_cpu(Fy), uy=to_cpu(uy), frames=frames)

# ---- анализ (на host: FFT малых рядов) ----
def analyze(res):
    Fx, Fy, uy, DL = res["Fx"], res["Fy"], res["uy"], res["DL"]
    n = len(uy); h = slice(n//2, n); Cd = 2*Fx/(U**2*DL); Cl = 2*Fy/(U**2*DL)
    sig = uy[h] - uy[h].mean(); npad = 8192
    freqs = np.fft.rfftfreq(npad, 1.0); amp = np.abs(np.fft.rfft(sig, n=npad))
    St = (freqs[1+int(np.argmax(amp[1:]))] if len(amp) > 2 else 0.0)*D0/U
    return dict(Cd=float(Cd[h].mean()), Clrms=float(Cl[h].std()), St=float(St), uyrms=float(uy[h].std()),
                Cl_s=Cl, uy=uy, freqs=freqs, amp=amp, h0=n//2)

# ---- визуализация (matplotlib на CPU; кадры уже host) ----
def _vort_np(u):
    return (np.roll(u[1], -1, 1)-np.roll(u[1], 1, 1))*0.5 - (np.roll(u[0], -1, 0)-np.roll(u[0], 1, 0))*0.5
def save_gif(res, name, anim_dir, sdf_tag, fps=24):
    import matplotlib; matplotlib.use("Agg")
    import matplotlib.pyplot as plt, matplotlib.patches as mpatches
    from PIL import Image
    nested = (res["mode"] == "nested"); frames = res["frames"]
    vm = np.nanpercentile(np.abs(np.where(solid0_np, np.nan, _vort_np(frames[len(frames)//2][1]))), 99.0)
    cm = plt.get_cmap("inferno").copy(); cm.set_bad("white"); pil = []
    for (t, u0, u1, u2) in frames:
        fig = plt.figure(figsize=(11.2, 5.7), dpi=96); ax = fig.add_subplot(111)
        ax.imshow(np.abs(np.where(solid0_np, np.nan, _vort_np(u0))), origin="lower", cmap=cm, vmin=0, vmax=vm,
                  extent=(0, Nx, 0, Ny), interpolation="nearest")
        ax.imshow(np.abs(np.where(solid1_np, np.nan, _vort_np(u1)))[1:-1, 1:-1], origin="lower", cmap=cm, vmin=0, vmax=vm,
                  extent=(ax1+0.5, bx1-0.5, ay1+0.5, by1-0.5), interpolation="nearest")
        if nested and u2 is not None:
            ax.imshow(np.abs(np.where(solid2_np, np.nan, _vort_np(u2)))[1:-1, 1:-1], origin="lower", cmap=cm, vmin=0, vmax=vm,
                      extent=(cx2a+0.25, cx2b-0.25, cy2a+0.25, cy2b-0.25), interpolation="nearest")
        ax.add_patch(mpatches.Rectangle((0.2, 0.2), Nx-0.4, Ny-0.4, fill=False, edgecolor="#4fc3f7", lw=1.6))
        ax.text(1.5, Ny-4, "уровень 0: Δx (крупный)", color="#4fc3f7", fontsize=9, weight="bold")
        ax.add_patch(mpatches.Rectangle((ax1, ay1), bx1-ax1, by1-ay1, fill=False, edgecolor="#aed581", lw=2.2))
        ax.text(ax1+1, by1-3.5, "уровень 1: Δx/2 (тело+след)", color="#aed581", fontsize=9, weight="bold")
        if nested:
            ax.add_patch(mpatches.Rectangle((cx2a, cy2a), cx2b-cx2a, cy2b-cy2a, fill=False, edgecolor="#ff8a65", lw=2.2))
            ax.text(cx2a+1, cy2b-3.5, "уровень 2: Δx/4 (у тела)", color="#ff8a65", fontsize=9, weight="bold")
        base = "Вложенный AMR" if nested else "AMR (одиночный блок 2×)"
        ax.set_title(f"{base}{sdf_tag} · поверх $|\\omega|$ · t={t}", fontsize=11)
        ax.set_xlabel("x [lu]"); ax.set_ylabel("y [lu]"); ax.set_xlim(0, Nx); ax.set_ylim(0, Ny)
        buf = _io.BytesIO(); fig.savefig(buf, format="png", bbox_inches="tight"); buf.seek(0); plt.close(fig)
        pil.append(Image.open(buf).convert("RGB").convert("P", palette=Image.ADAPTIVE))
    path = os.path.join(anim_dir, name)
    pil[0].save(path, save_all=True, append_images=pil[1:], duration=int(1000/fps), loop=0, optimize=True)
    print("saved", name, len(frames), "кадров")

LAB = ["без AMR (D=16)", "AMR 2× (D=32)", "AMR 2×+4× (D=64)"]
def print_table(AS, header):
    print(f"\n=== {header} ===")
    print(f"{'версия':22}{'D у тела':>9}{'St':>8}{'<Cd>':>9}{'rms Cl':>9}{'rms u_y':>9}")
    for lab, DL, a in zip(LAB, [16, 32, 64], AS):
        print(f"{lab:22}{DL:>9}{a['St']:>8.3f}{a['Cd']:>9.3f}{a['Clrms']:>9.4f}{a['uyrms']:>9.4f}")
    print(f"{'литература Re≈150':22}{'—':>9}{0.183:>8.3f}{1.330:>9.3f}{'~0.3':>9}{'':>9}")
    print("(лит.: St≈0.18, Cd≈1.3 для цилиндра при Re≈150; Williamson 1996, Henderson 1995)")
def probe_figure(AS, R2, path, suptitle):
    import matplotlib; matplotlib.use("Agg")
    import matplotlib.pyplot as plt, matplotlib.patches as mpatches
    cols = ["#9e9e9e", "#1f6feb", "#d1495b"]
    fig = plt.figure(figsize=(13.5, 8.2)); gs = fig.add_gridspec(2, 2, hspace=0.32, wspace=0.22)
    axA = fig.add_subplot(gs[0, 0]); mid = R2["frames"][len(R2["frames"])//2]
    vmf = np.nanpercentile(np.abs(np.where(solid0_np, np.nan, _vort_np(mid[1]))), 99)
    cmf = plt.get_cmap("inferno").copy(); cmf.set_bad("white")
    axA.imshow(np.abs(np.where(solid0_np, np.nan, _vort_np(mid[1]))), origin="lower", cmap=cmf, vmin=0, vmax=vmf,
               extent=(0, Nx, 0, Ny), interpolation="nearest")
    axA.add_patch(mpatches.Rectangle((ax1, ay1), bx1-ax1, by1-ay1, fill=False, edgecolor="#aed581", lw=1.8))
    axA.add_patch(mpatches.Rectangle((cx2a, cy2a), cx2b-cx2a, cy2b-cy2a, fill=False, edgecolor="#ff8a65", lw=1.8))
    axA.set_title("$|\\omega|$ + рамки зон 2×/4×"); axA.set_xlabel("x [lu]"); axA.set_ylabel("y [lu]")
    axB = fig.add_subplot(gs[0, 1])
    for a, lab, c in zip(AS, LAB, cols): axB.plot(a["Cl_s"][a["h0"]:], lw=0.8, color=c, label=lab)
    axB.set_title("Подъёмная сила $C_l(t)$ (зонд силы)"); axB.set_xlabel("шаг"); axB.set_ylabel("$C_l$")
    axB.legend(fontsize=8); axB.grid(alpha=0.3)
    axC = fig.add_subplot(gs[1, 0])
    for a, lab, c in zip(AS, LAB, cols):
        sp = a["amp"]/a["amp"][1:].max(); axC.plot(a["freqs"]*D0/U, sp, lw=1.0, color=c, label=lab)
    axC.axvline(0.183, color="k", ls="--", alpha=0.6, label="лит. St≈0.18"); axC.set_xlim(0, 0.6)
    axC.set_title("Спектр $u_y$ в следе → St"); axC.set_xlabel("St = f·D/U"); axC.set_ylabel("норм. ампл.")
    axC.legend(fontsize=8); axC.grid(alpha=0.3)
    axD = fig.add_subplot(gs[1, 1]); x = np.arange(3); w = 0.35
    axD.bar(x-w/2, [a["St"] for a in AS], w, color="#1f6feb", label="St")
    axD.bar(x+w/2, [a["Cd"] for a in AS], w, color="#fb8c00", label="<Cd>")
    axD.axhline(0.183, color="#1f6feb", ls="--", alpha=0.6); axD.axhline(1.33, color="#fb8c00", ls="--", alpha=0.6)
    axD.set_xticks(x); axD.set_xticklabels(["без AMR", "2×", "2×+4×"]); axD.set_title("St и <Cd> (пунктир — лит.)")
    axD.legend(fontsize=8); axD.grid(alpha=0.3, axis="y")
    fig.suptitle(suptitle, fontweight="bold"); fig.savefig(path, dpi=110, bbox_inches="tight"); plt.close(fig)
    print("saved", os.path.basename(path))