CFDManager/docs/theory/solver_2x_sdf/config.py

# Конфигурация модели 2×+SDF. Всё в одном месте, чтобы математику/режимы было легко менять.
# Единицы — решёточные (lu — длина воксель, ts — шаг по времени).
import os
from dataclasses import dataclass
from typing import Optional

@dataclass
class Config:
    # --- сетка / геометрия L0 (умеренно расширены: меньше блокировка + больше следа) ---
    Nx: int = 174            # домен по x (цилиндр на 40 → ~8.4D вниз по потоку до выхода)
    Ny: int = 90             # домен по y (стенки free-slip; блокировка β=D/Ny — варьируется серией)
    D: int = 16              # диаметр цилиндра на грубом уровне L0
    cx: int = 40
    cy: Optional[int] = None # центр по вертикали; по умолчанию Ny//2 (см. __post_init__)

    # --- физика ---
    U: float = 0.07          # скорость набегающего потока
    Re: float = 150.0        # число Рейнольдса по D и U

    # --- время ---
    steps: int = 100000
    ramp: int = 1000         # smoothstep-разгон входной скорости (гасит импульсный старт)
    frame_every: int = 34    # период сохранения кадров для гифки

    # --- патч уровня L1 (в координатах L0); refine=2 → разрешение Δx/2 ---
    ax1: Optional[int] = None  # по умолчанию cx − 1.5D (запас перед телом; см. __post_init__)
    bx1: Optional[int] = None  # по умолчанию cx + 6.25D (расширен в след)
    ay1: Optional[int] = None  # по умолчанию cy − 2D (центрирован по cy)
    by1: Optional[int] = None  # по умолчанию cy + 2D
    refine: int = 2          # коэффициент измельчения L1 (модель «2×» рассчитана на 2)

    # --- зонд скорости в следе: px = cx + probe_dx*D, py = cy (берётся с тонкой сетки L1) ---
    probe_dx: int = 3

    # --- стартовое возмущение (триггер вихревой дорожки): поперечный sin-импульс входа ---
    pert_amp: float = 0.1    # амплитуда u_y входа, доля от U (0 — выключить)
    pert_dur: int = 300      # длительность импульса в шагах; окно [ramp, ramp+pert_dur]

    # --- кросс-чек силы: интеграл тензора напряжений σ·n по окружности R1+δ на L1 (forces_stress.py) ---
    stress_every: int = 50   # период выборки в шагах L0 (0 — выключить кросс-чек)
    stress_ntheta: int = 256 # число точек контура по углу
    stress_delta: float = 2.0  # отступ контура от поверхности (в тонких ячейках L1);
                               # δ=2 → билинейный стенсиль не трогает Bouzidi-зону (q_max≈0.98)

    # --- ГУ выхода: режим поперечной скорости на Zou-He давление-выходе ---
    #   "zero"        — классика: uy=0 на выходе (жёстко; отражает вихри дорожки)
    #   "extrapolate" — uy берётся из соседнего столбца (нуль-градиент); ρ-якорь сохраняется
    outlet_uy: str = "zero"

    # --- губка (absorbing layer) перед выходом: плавный рост вязкости в последних sponge_len
    # столбцах L0 (smoothstep, ν → ν·sponge_nu_mult). Гасит вихри И акустику до прихода на
    # давление-выход. Зачем: канал «вход-скорость + выход-давление» — недодемпфированный
    # акустический резонатор (затухание моды ∝ ν(π/Nx)²); на длинных доменах мода раскачивается
    # стартовым транзиентом → смещённый режим ⟨ρ⟩≈1.66 или взрыв (факторное исследование,
    # эксперимент №1: F1/F2/F3/F5). 0 — выключена.
    sponge_len: int = 0
    sponge_nu_mult: float = 30.0

    # --- исполнение ---
    use_cuda_graph: bool = True

    def __post_init__(self):
        # производная геометрия по вертикали: центр и патч следуют за Ny (для серии по блокировке)
        if self.cy is None:
            self.cy = self.Ny // 2
        if self.ay1 is None:
            self.ay1 = self.cy - 2 * self.D
        if self.by1 is None:
            self.by1 = self.cy + 2 * self.D
        # производная геометрия по горизонтали: патч следует за cx (для факторов вход/выход)
        if self.ax1 is None:
            self.ax1 = self.cx - (3 * self.D) // 2
        if self.bx1 is None:
            self.bx1 = self.cx + (25 * self.D) // 4
        assert 1 <= self.ay1 < self.by1 <= self.Ny - 2, (
            f"патч L1 [{self.ay1},{self.by1}] должен лежать строго внутри (0,{self.Ny-1}) — "
            f"specular-ряды стенок не накрывать; минимально допустимое Ny ≈ {4*self.D + 6}")
        assert 2 <= self.ax1 < self.cx < self.bx1 <= self.Nx - 2, (
            f"патч L1 [{self.ax1},{self.bx1}] должен лежать строго внутри (1,{self.Nx-1}) "
            f"и накрывать тело (cx={self.cx})")
        assert self.outlet_uy in ("zero", "extrapolate"), self.outlet_uy
        if self.sponge_len:
            assert self.Nx - 1 - self.sponge_len > self.bx1, (
                f"губка (последние {self.sponge_len} столбцов) не должна накрывать патч L1 "
                f"(bx1={self.bx1}, Nx={self.Nx})")

    # ---------- производные величины ----------
    @property
    def nu(self):  return self.U * self.D / self.Re          # кинематическая вязкость (lu²/ts)
    @property
    def tau0(self): return self.nu / (1.0/3.0) + 0.5         # время релаксации на L0 (cs²=1/3)
    @property
    def tau1(self):  # связка τ при измельчении в r раз: ν одна и та же ⇒ τ_f = r(τ_c − ½) + ½
        return self.refine * (self.tau0 - 0.5) + 0.5
    @property
    def beta0(self): return 1.0 / (2.0 * self.tau0)          # β = 1/(2τ) для KBC/BGK на L0
    @property
    def beta1(self): return 1.0 / (2.0 * self.tau1)          # β на L1
    @property
    def R01(self):   # масштаб неравновесной части грубый→тонкий: f^neq ∝ τ·δt ⇒ (τ_f/τ_c)/r
        return self.tau1 / (self.refine * self.tau0)
    @property
    def Rf01(self): return 1.0 / self.R01                    # тонкий→грубый
    @property
    def DL(self):   return self.refine * self.D              # эффективный диаметр у тела (на L1)
    @property
    def px(self):   return self.cx + self.probe_dx * self.D
    @property
    def py(self):   return self.cy

def from_env(**overrides):
    """Создать Config с учётом переменных окружения и явных overrides.
    Геометрия (AMR_NY/AMR_NX/AMR_CX) попадает в overrides ДО конструктора —
    производные cy/ay1/by1/ax1/bx1 строятся в __post_init__."""
    if "Ny" not in overrides and os.environ.get("AMR_NY"):
        overrides["Ny"] = int(os.environ["AMR_NY"])
    if "Nx" not in overrides and os.environ.get("AMR_NX"):
        overrides["Nx"] = int(os.environ["AMR_NX"])
    if "cx" not in overrides and os.environ.get("AMR_CX"):
        overrides["cx"] = int(os.environ["AMR_CX"])
    if "outlet_uy" not in overrides and os.environ.get("AMR_OUTLET_UY"):
        overrides["outlet_uy"] = os.environ["AMR_OUTLET_UY"]
    if "sponge_len" not in overrides and os.environ.get("AMR_SPONGE_LEN"):
        overrides["sponge_len"] = int(os.environ["AMR_SPONGE_LEN"])
    if "sponge_nu_mult" not in overrides and os.environ.get("AMR_SPONGE_NU"):
        overrides["sponge_nu_mult"] = float(os.environ["AMR_SPONGE_NU"])
    cfg = Config(**overrides)
    if os.environ.get("AMR_STEPS"):
        cfg.steps = int(os.environ["AMR_STEPS"])
    if "AMR_CUDAGRAPH" in os.environ:
        cfg.use_cuda_graph = (os.environ["AMR_CUDAGRAPH"] != "0")
    if os.environ.get("AMR_STRESS_EVERY"):
        cfg.stress_every = int(os.environ["AMR_STRESS_EVERY"])
    if os.environ.get("AMR_STRESS_NTHETA"):
        cfg.stress_ntheta = int(os.environ["AMR_STRESS_NTHETA"])
    if os.environ.get("AMR_STRESS_DELTA"):
        cfg.stress_delta = float(os.environ["AMR_STRESS_DELTA"])
    return cfg