# Конфигурация модели 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