CFDManager/docs/theory/solver_2x_sdf/diagnostics.py

# Диагностика по собранным рядам сил/скорости. numpy здесь — только host-постобработка малых рядов
# (FFT/статистика), это НЕ вычислительный бэкенд солвера.
import numpy as np
import backend as B

def strouhal(uy_host, D, U, pad=8192):
    """Число Струхаля из ряда поперечной скорости в следе.
    УЛУЧШЕНО: параболическая (3-точечная) интерполяция пика спектра → суб-биновая точность.
    (Раньше брался ближайший FFT-бин: шаг по St ≈ D/(U·pad) ≈ 0.028, из-за чего St 'застревал'.)"""
    sig = uy_host[len(uy_host)//2:]                      # установившийся режим — вторая половина
    sig = sig - sig.mean()
    amp = np.abs(np.fft.rfft(sig, n=pad))
    if len(amp) <= 3:
        return 0.0, np.fft.rfftfreq(pad, 1.0), amp
    k = 1 + int(np.argmax(amp[1:]))                      # индекс пика (без DC)
    if 1 <= k < len(amp) - 1:                            # параболическая интерполяция вершины
        a0, a1, a2 = amp[k-1], amp[k], amp[k+1]
        denom = (a0 - 2*a1 + a2)
        delta = 0.5*(a0 - a2)/denom if denom != 0 else 0.0
    else:
        delta = 0.0
    f_peak = (k + delta) / pad
    return f_peak * D / U, np.fft.rfftfreq(pad, 1.0), amp

def analyze(res, cfg):
    Fx = B.to_cpu(res["Fx"]); Fy = B.to_cpu(res["Fy"]); uy = B.to_cpu(res["uy"])
    DL, U, D = res["DL"], cfg.U, cfg.D
    n = len(uy); h = slice(n//2, n)
    Cd = 2.0*Fx/(U**2*DL); Cl = 2.0*Fy/(U**2*DL)
    St, freqs, amp = strouhal(uy, D, U)
    a = dict(Cd=float(Cd[h].mean()), Clrms=float(Cl[h].std()), St=float(St),
             uyrms=float(uy[h].std()), Cd_s=Cd, Cl_s=Cl, uy=uy, freqs=freqs, amp=amp, h0=n//2)
    if "Tz" in res:                                       # момент: Cm = 2Tz/(U²·DL²); для цилиндра ⟨Cm⟩≈0
        Cm = 2.0*B.to_cpu(res["Tz"])/(U**2*DL*DL)
        a["Cm"] = float(Cm[h].mean()); a["Cmrms"] = float(Cm[h].std()); a["Cm_s"] = Cm
    K = int(res.get("stress_every", 0) or 0)              # кросс-чек σ·n (редкая сетка по времени)
    if K and len(res.get("Fx_st", ())) > 3:
        Fxs = B.to_cpu(res["Fx_st"]); Fys = B.to_cpu(res["Fy_st"]); Tzs = B.to_cpu(res["Tz_st"])
        m = len(Fxs); hs = slice(m//2, m)
        Cds = 2.0*Fxs/(U**2*DL); Cls = 2.0*Fys/(U**2*DL); Cms = 2.0*Tzs/(U**2*DL*DL)
        a.update(Cd_st=float(Cds[hs].mean()), Clrms_st=float(Cls[hs].std()),
                 Cm_st=float(Cms[hs].mean()), Cd_st_s=Cds, Cl_st_s=Cls, stress_every=K)
    return a

def print_table(a, cfg):
    # поправка на блокировку канала (приведение к скорости в зазоре U/(1−β)):
    #   St — кинематическая → ×(1−β);  Cd, rms Cl ~ скорость² → ×(1−β)²
    beta = cfg.D / cfg.Ny
    St_c = a['St'] * (1 - beta)
    Cd_c = a['Cd'] * (1 - beta)**2
    Cl_c = a['Clrms'] * (1 - beta)**2
    cm = a.get('Cm')
    cms = f"{cm:>10.5f}" if cm is not None else f"{'—':>10}"
    print("\n=== Модель 2×+SDF — установившийся режим ===")
    print(f"{'версия':18}{'D у тела':>9}{'St':>8}{'<Cd>':>9}{'rms Cl':>9}{'<Cm>':>10}{'rms u_y':>9}")
    print(f"{'2×+SDF (raw)':18}{cfg.DL:>9}{a['St']:>8.3f}{a['Cd']:>9.3f}{a['Clrms']:>9.4f}{cms}{a['uyrms']:>9.4f}")
    print(f"{'2×+SDF (×блок.)':18}{cfg.DL:>9}{St_c:>8.3f}{Cd_c:>9.3f}{Cl_c:>9.4f}{'—':>10}{'—':>9}")
    print(f"{'лит. Re≈150':18}{'—':>9}{0.183:>8.3f}{1.330:>9.3f}{'~0.3':>9}{0.0:>10.1f}{'':>9}")
    print(f"(блокировка β=D/Ny={beta:.3f}; поправка: St×(1−β), Cd/rmsCl×(1−β)². "
          "Лит. — безграничный цилиндр; ⟨Cm⟩≈0 — контроль симметрии; см. README)")

def print_crosscheck(a):
    """Сравнение двух НЕЗАВИСИМЫХ считываний силы: GMEM (обмен импульсом по линкам) vs ∮σ·n ds
    (интеграл тензора напряжений по контуру R+δ). Сходство — прямое свидетельство корректности
    считывания; сравниваем средние (мгновенные ряды σ·n сдвинуты по фазе кольцом R..R+δ)."""
    if "Cd_st" not in a:
        print("\n[кросс-чек σ·n: выключен (stress_every=0) или слишком мало точек]")
        return
    dcd = abs(a["Cd_st"] - a["Cd"]) / max(abs(a["Cd"]), 1e-12)
    dcl = abs(a["Clrms_st"] - a["Clrms"]) / max(abs(a["Clrms"]), 1e-12)
    cm = a.get("Cm", float("nan"))
    print("\n=== Кросс-чек считывания силы: GMEM vs ∮σ·n ds (установившийся режим) ===")
    print(f"{'метод':26}{'<Cd>':>9}{'rms Cl':>9}{'<Cm>':>10}")
    print(f"{'GMEM (обмен импульсом)':26}{a['Cd']:>9.3f}{a['Clrms']:>9.4f}{cm:>10.5f}")
    print(f"{'∮σ·n ds (контур R+δ)':26}{a['Cd_st']:>9.3f}{a['Clrms_st']:>9.4f}{a['Cm_st']:>10.5f}")
    verdict = "СОГЛАСОВАНО" if (dcd < 0.05 and dcl < 0.10) else "ПРОВЕРИТЬ"
    print(f"расхождение: dCd={dcd*100:.1f}% (порог 5%), d(rmsCl)={dcl*100:.1f}% (порог 10%) → {verdict}")

def fit_blockage(betas, vals):
    """Экстраполяция величины серии к β→0 двумя моделями: val ≈ a + b·β и val ≈ a + b·β².
    (Теория solid-blockage даёт ведущий член O(β²), практика поправки (1−β)² содержит и линейный.)
    Спред интерсептов |a_lin − a_quad| — оценка модельной неопределённости экстраполяции."""
    b = np.asarray(betas, float); v = np.asarray(vals, float)
    sl, al = np.polyfit(b, v, 1)
    sq, aq = np.polyfit(b**2, v, 1)
    return dict(lin=(float(al), float(sl)), quad=(float(aq), float(sq)),
                spread=float(abs(al - aq)))

def print_blockage_table(rows, fits):
    """Сводка серии по блокировке. rows — список dict(Ny, beta, St, Cd, Clrms, Cm, Cd_st);
    fits — dict('Cd'/'Clrms'/'St' → результат fit_blockage)."""
    print("\n=== Серия по блокировке: D фикс., Ny варьируется ===")
    print(f"{'Ny':>5}{'β=D/Ny':>9}{'St':>8}{'St(1−β)':>9}{'<Cd>':>9}{'rms Cl':>9}{'<Cm>':>10}{'Cd σ·n':>9}")
    for r in rows:
        cdst = f"{r['Cd_st']:>9.3f}" if r.get("Cd_st") is not None else f"{'—':>9}"
        print(f"{r['Ny']:>5}{r['beta']:>9.3f}{r['St']:>8.3f}{r['St']*(1-r['beta']):>9.3f}"
              f"{r['Cd']:>9.3f}{r['Clrms']:>9.4f}{r['Cm']:>10.5f}{cdst}")
    print("экстраполяция β→0 (интерсепт лин. фита / квадр. фита; спред — неопределённость):")
    print(f"  Cd(0)     = {fits['Cd']['lin'][0]:.3f} / {fits['Cd']['quad'][0]:.3f} "
          f"(спред {fits['Cd']['spread']:.3f};  лит. безгранич. 1.33)")
    print(f"  rms Cl(0) = {fits['Clrms']['lin'][0]:.3f} / {fits['Clrms']['quad'][0]:.3f} "
          f"(спред {fits['Clrms']['spread']:.3f};  лит. ~0.3)")
    print(f"  St(0)     = {fits['St']['lin'][0]:.3f} / {fits['St']['quad'][0]:.3f} "
          f"(спред {fits['St']['spread']:.3f};  лит. 0.183)")

def convergence_report(res, cfg, nwin=10):
    """Эволюция по окнам времени: видно, НАСЫЩЕНО ли решение или ДРЕЙФУЕТ.
    Ключевое: ⟨ρ⟩ (средняя плотность) и Cd, нормированный на реальную ⟨ρ⟩ (дрейф-устойчивый Cd).
    Если ⟨ρ⟩ уплывает от 1 — это дрейф массы из-за ГУ входа/выхода, а Cd_raw искажён нормировкой на ρ=1."""
    Fx = B.to_cpu(res["Fx"]); Fy = B.to_cpu(res["Fy"]); uy = B.to_cpu(res["uy"]); rho = B.to_cpu(res["rho"])
    n = len(Fx); U = cfg.U; DL = res["DL"]
    if n < nwin * 2:
        return
    Cd = 2.0*Fx/(U**2*DL); Cl = 2.0*Fy/(U**2*DL)
    w = n // nwin
    print("\n=== Сходимость по времени / дрейф (по окнам) ===")
    print(f"{'окно':>4}{'шаги':>16}{'<ρ>':>8}{'<Cd>raw':>9}{'<Cd>/ρ':>9}{'rms Cl':>9}{'rms u_y':>9}")
    rows = []
    for k in range(nwin):
        s = slice(k*w, (k+1)*w if k < nwin-1 else n)
        rm = float(rho[s].mean()); cdr = float(Cd[s].mean())
        cdn = cdr / rm if rm != 0 else float("nan")
        clr = float(Cl[s].std()); uyr = float(uy[s].std())
        rows.append((rm, cdr, cdn, clr, uyr))
        print(f"{k+1:>4}{f'{k*w}-{(k+1)*w}':>16}{rm:>8.3f}{cdr:>9.3f}{cdn:>9.3f}{clr:>9.4f}{uyr:>9.4f}")
    # вердикт по последним двум окнам
    (rm1, cdr1, cdn1, clr1, _), (rm0, _, _, clr0, _) = rows[-1], rows[-2]
    drho = rm1 - rm0; dcl = clr1 - clr0
    verdict = []
    if abs(rm1 - 1.0) > 0.02:
        # «стабильна, но не там»: система села на смещённую ветвь (акустическая накачка массы;
        # ⟨Cd⟩/Cl на ней НЕвалидны для сравнения с литературой) — см. факторное исследование
        verdict.append(f"плотность: ⟨ρ⟩={rm1:.3f} → РЕЖИМ СМЕЩЁН (⟨ρ⟩ далеко от 1; прогон невалиден)")
    else:
        verdict.append(f"плотность: ⟨ρ⟩={rm1:.3f}, Δ за окно={drho:+.4f} → "
                       + ("ДРЕЙФ массы (правь ГУ выхода: давление/Zou-He)" if abs(drho) > 1e-3 else "стабильна"))
    verdict.append(f"rms Cl: Δ за окно={dcl:+.4f} → "
                   + ("ещё растёт (не насыщено)" if dcl > 1e-3 else "насыщено/стабильно"))
    verdict.append(f"дрейф-устойчивый Cd (⟨Cd⟩/ρ) в последнем окне = {cdn1:.3f}")
    print("вердикт: " + "; ".join(verdict))