CFDManager/docs/theory/kbc_lbm.ipynb

Метод KBC (энтропийный Lattice Boltzmann): от формул к валидированной GPU-модели

Karlin–Bösch–Chikatamarla (KBC) — энтропийная multi-relaxation-time формулировка решёточного метода Больцмана (LBM), дающая безусловную нелинейную устойчивость при больших числах Рейнольдса без явной подсеточной модели турбулентности (работает как implicit LES).

Этот документ — рабочий разбор метода и история эволюции 2D-макета: от формул и самопроверок через классические бенчмарки (Тейлора–Грина, сдвиговый слой, каверна), блочное измельчение сетки (AMR) и границу SDF+Bouzidi — к финальной валидированной GPU-модели обтекания цилиндра solver_2x_sdf (D2Q9, KBC-N1, L0+L1, CUDA-графы) с двумя независимыми способами считывания силы.

Как устроен документ. Ноутбук не исполняет код — только теория и готовые артефакты (картинки, гифки, числа). Весь код живёт рядом и запускается на GPU-сервере:

Где	Что
demos_gpu/	GPU-прогоны всех демо и самопроверок (см. demos_gpu/README.md)
solver_2x_sdf/	финальная модель (см. solver_2x_sdf/README.md)
figures/, anim/, */out/	артефакты, встроенные ниже

Демо гоняются теми же модулями (решётка, равновесие, столкновение, перенос, AMR-связка), что и финальная модель, — каждый бенчмарк ниже одновременно валидирует компонент финальной модели.

Сводная таблица обозначений

Везде решёточные единицы: шаг сетки \Delta x = 1 (lu), шаг по времени \Delta t = 1 (ts).

Символ	Название	Смысл / формула	Диапазон
f_i	популяция	плотность частиц со скоростью \mathbf c_i	f_i>0
\mathbf c_i	дискретная скорость	направление переноса i	компоненты \in\{-1,0,1\}
w_i	вес квадратуры	\sum_i w_i=1	\{4/9,1/9,1/36\} (D2Q9)
c_s	скорость звука решётки	c_s^2=1/3	c_s=1/\sqrt3
\rho, \mathbf u	плотность, скорость	\rho=\sum_i f_i, \rho\mathbf u=\sum_i \mathbf c_i f_i	\rho\approx1, $
\tau, \beta	релаксация	\nu=c_s^2(\tau-\tfrac12), \beta=1/(2\tau)	\tau>1/2
\gamma	энтропийный стабилизатор	релаксация высших мод	\approx2
f_i^{\mathrm{eq}}	равновесие	максимум энтропии при \rho,\mathbf u	>0
k_i,s_i,h_i	разложение	кинематич./сдвиг/высшие моды	k+s+h=f
\Delta s_i,\Delta h_i	неравновесие	P_s(f{-}f^{eq}), остальное	малы
H	H-функция	H=\sum_i f_i\ln(f_i/w_i)	не растёт
$\langle X	Y\rangle$	энтропийное скал. произв.	\sum_i X_iY_i/f_i^{eq}
T,N,\Pi_{xy}	натуральные 2-е моменты	\Pi_{xx}{+}\Pi_{yy}, \Pi_{xx}{-}\Pi_{yy}, сдвиг	—
\Pi^{neq}_{\alpha\beta}	неравновесный тензор	\sum_i c_{i\alpha}c_{i\beta}(f_i-f_i^{eq})	—
\sigma_{\alpha\beta}	тензор напряжений	-p'\delta_{\alpha\beta}-(1-\tfrac{1}{2\tau})\Pi^{neq,dev}_{\alpha\beta}	—
q	доля стенки (Bouzidi)	положение стенки на линке из SDF	(0,1)
C_d, C_l, C_m	коэффициенты	\tfrac{2F_x}{U^2 D}, \tfrac{2F_y}{U^2 D}, \tfrac{2T_z}{U^2 D^2}	—
\beta_{бл}	блокировка канала	D/N_y (стеснение стенками)	\to0 — безграничный
\mathrm{Re}, \mathrm{Ma}, \mathrm{St}	подобие	UL/\nu, $	\mathbf u

---

Часть I. Теория

1. Основы LBM и решётка D2Q9

LBM отслеживает популяции f_i(\mathbf x, t) — плотность фиктивных частиц, движущихся в узле \mathbf x со скоростью \mathbf c_i. Динамика — чередование переноса (streaming) и столкновения (collision):

f_i(\mathbf x + \mathbf c_i\,\Delta t,\; t+\Delta t) = f_i(\mathbf x,t) + \Omega_i.

Решётка D2Q9: покоящаяся частица, 4 осевые и 4 диагональные скорости; c_s^2=1/3; веса w_0=4/9, осевые 1/9, диагональные 1/36 дают изотропию 4-го порядка. Макромоменты: \rho=\sum_i f_i, \rho\mathbf u=\sum_i \mathbf c_i f_i.

Самопроверка (demos_gpu/checks_gpu.py, исполняется на компонентах финальной модели — solver_2x_sdf/lattice.py): \sum w_i=1, \sum w_i\mathbf c_i=0, \sum w_i c_{i\alpha}c_{i\beta}=c_s^2\delta_{\alpha\beta}, \mathbf c_{\bar i}=-\mathbf c_i — все тождества проходят (см. demos_gpu/out/checks.txt).

2. Равновесие: полиномиальное и энтропийное

Полиномиальное (разложение максвеллиана до O(u^2), стандарт BGK): $$ f_i^{\mathrm{eq,poly}} = w_i,\rho\left[1 + \frac{\mathbf c_i\cdot\mathbf u}{c_s^2}

• \frac{(\mathbf c_i\cdot\mathbf u)^2}{2c_s^4} - \frac{\mathbf u^2}{2c_s^2}\right]. $$

Энтропийное (точный максимум дискретной H-функции при связях \rho, \rho\mathbf u) для D2Q9/D3Q27 имеет замкнутую product-form: $$ f_i^{\mathrm{eq}} = \rho,w_i\prod_{\alpha\in{x,y}} A(u_\alpha),B(u_\alpha)^{,c_{i\alpha}},\qquad A(u) = 2-\sqrt{1+3u^2},\qquad B(u) = \frac{2u+\sqrt{1+3u^2}}{1-u}. $$

✓ Сверено с первоисточником: форма и коэффициенты в точности совпадают с уравнениями (7)–(11) Bösch, Chikatamarla, Karlin, Phys. Rev. E 92, 043309 (2015) (восходит к Ansumali, Karlin, Öttinger, Europhys. Lett. 63, 798, 2003). Реализация — solver_2x_sdf/equilibrium.py (показатели c_{i\alpha}\in\{-1,0,1\}, поэтому вместо pow берётся B или 1/B — быстрее и совместимо с CUDA-графами).

Самопроверка (checks_gpu.py): энтропийное равновесие воспроизводит \rho и \rho\mathbf u точно (до машинного нуля выбранной точности), а второй момент при малом \mathrm{Ma} близок к \rho(c_s^2\delta_{\alpha\beta}+u_\alpha u_\beta) — тензору давления Навье–Стокса. В отличие от полинома, product-form гарантирует положительность и H-теорему.

3. Разложение популяций f = k + s + h

Сердце KBC — разбиение 9-вектора популяций по моментным подпространствам. Натуральные вторые моменты D2Q9: след T=\Pi_{xx}{+}\Pi_{yy}, разность нормальных напряжений N=\Pi_{xx}{-}\Pi_{yy}, сдвиг \Pi_{xy}; плюс третьи (Q_{xxy},Q_{xyy}) и четвёртый (\Pi_{xxyy}) моменты:

• k_i — кинематическая часть: только сохраняющиеся моменты (\rho, j_x, j_y);
• s_i — сдвиговая часть: девиаторный тензор напряжений (N и \Pi_{xy});
• h_i — высшие («ghost») моды: T, третьи и четвёртый моменты.

Строится моментная матрица M_{ai}=\phi^{(a)}(\mathbf c_i) из 9 независимых мономов; проектор на группу моментов G: P_G = M^{-1} D_G M. По построению P_k+P_s+P_h=I, разложение точное.

✓ Сверено с первоисточником: назначение s задаёт вариант KBC; вариант с s=\{N,\Pi_{xy}\} (только девиатор) — канонический «KBC-d» / N1 (Karlin, Bösch, Chikatamarla, PRE 90, 031302(R), 2014; Bösch et al., PRE 92, 043309, 2015 — там же семейство s\in\{d,\,d{+}t,\,d{+}q,\,d{+}t{+}q\}). Именно он реализован в solver_2x_sdf/lattice.py (проектор Ps на строки \{c_x^2{-}c_y^2,\,c_xc_y\}).

Самопроверка (checks_gpu.py): Ps идемпотентен, \Delta s несёт только моменты \{N,\Pi_{xy}\} (остальные строки M\Delta s — нули).

4. Оператор столкновений KBC и стабилизатор \gamma

Так как f и f^{eq} имеют одинаковые сохраняющиеся моменты, \Delta k=0 и f-f^{eq}=\Delta s+\Delta h, где \Delta s = P_s(f-f^{eq}).

Столкновение KBC (через «зеркальное состояние» f^{mirr}): $$ f_i' = (1-\beta),f_i + \beta,f_i^{mirr} ;;\Longleftrightarrow;; \boxed{,f_i' = f_i - \beta,\bigl(2,\Delta s_i + \gamma,\Delta h_i\bigr),},\qquad \nu = c_s^2\Bigl(\frac{1}{2\beta}-\frac12\Bigr);\Rightarrow;\beta=\frac{1}{2\tau}. $$ Множитель 2 перед \Delta s фиксирует сдвиговую вязкость; высшие моды релаксируются отдельным \gamma.

Стабилизатор \gamma выбирается локально (в каждом узле на каждом шаге) из условия неувеличения дискретной энтропии; аналитическое приближение корня: $$ \gamma^* = \frac{1}{\beta} - \Bigl(2-\frac{1}{\beta}\Bigr) \frac{\langle \Delta s,|,\Delta h\rangle}{\langle \Delta h,|,\Delta h\rangle}, \qquad \langle X|Y\rangle = \sum_i \frac{X_i Y_i}{f_i^{eq}}. $$

✓ Сверено с первоисточником: обе формулы (включая зеркальную форму и связь \nu(\beta)) дословно совпадают с Karlin, Bösch, Chikatamarla, PRE 90, 031302(R) (2014) и Bösch et al., PRE 92, 043309 (2015), ур. (24)–(25). Реализация — solver_2x_sdf/collision.py (с регуляризацией \gamma\to2 при \langle\Delta h|\Delta h\rangle\to0).

При \gamma=2: 2\Delta s+2\Delta h=2(f-f^{eq}) — KBC сводится к BGK с \omega=2\beta. Самопроверка (checks_gpu.py): KBC$(\gamma{=}2)\equiv$ BGK до машинного нуля; столкновение сохраняет \rho и \rho\mathbf u.

Популярно. Разложение k+s+h — три «слоя» движения: перенос массы/импульса, вязкий сдвиг (трение слоёв — задаёт вязкость и Re) и мелкая «рябь» высших мод, где на грубой сетке копится неустойчивость. KBC гасит сдвиг фиксированным темпом 2\beta (физика), а рябь — адаптивным \gamma\beta: в спокойных зонах \gamma\approx2 (как BGK), у резких градиентов \gamma отклоняется и добавляет ровно столько диссипации, сколько разрешает энтропия. Это «локальный термостат потока» — встроенная, без подгоночных констант, модель мелких масштабов (implicit LES).

5. Перенос и связь вязкости с \beta; H-функция

Перенос — сдвиг каждой популяции на её \mathbf c_i (pull-схема, solver_2x_sdf/streaming.py). Вязкость: \nu = c_s^2(\tfrac{1}{2\beta}-\tfrac12), т.е. \tau=\nu/c_s^2+\tfrac12. Число Рейнольдса \mathrm{Re}=UL/\nu; Маха \mathrm{Ma}=U/c_s\lesssim0.1 (ошибка сжимаемости \propto \mathrm{Ma}^2).

Дискретная энтропия H=\sum_i f_i\ln(f_i/w_i) — «физический компас» устойчивости: из условия «H не растёт» выведен \gamma; в демо Тейлора–Грина ниже \langle H\rangle(t) монотонно не возрастает.

6. Алгоритм одного шага KBC

для каждого узла:
  1) макро-моменты:   rho = sum f_i,   rho*u = sum c_i f_i
  2) равновесие:      f_i^eq  (энтропийное product-form)
  3) неравновесие:    Ds = P_s (f - f^eq),   Dh = (f - f^eq) - Ds
  4) стабилизатор:    gamma = 1/beta - (2 - 1/beta) <Ds|Dh> / <Dh|Dh>
  5) столкновение:    f <- f - beta (2 Ds + gamma Dh)
  6) перенос:         f_i(x + c_i) <- f_i(x)
  7) граничные условия (bounce-back / Bouzidi / Zou-He / free-slip)

---

Часть II. Валидационные демо

Четыре классических теста; каждый гоняется скриптом из demos_gpu/ на тех же модулях, что финальная модель (решётка/равновесие/KBC/перенос из solver_2x_sdf). Числа ниже цитируются из demos_gpu/out/*.txt — перегенерируются при каждом прогоне на GPU-сервере.

7. Демо-1: вихрь Тейлора–Грина — валидация вязкости и H-теоремы

Точное решение несжимаемого Навье–Стокса в периодической области: u_x=-U_0\cos(kx)\sin(ky), u_y=U_0\sin(kx)\cos(ky), k=2\pi/N; кинетическая энергия затухает как E(t)=E_0\,e^{-4\nu k^2 t}. Прогоняем KBC (N=96, \nu=0.0125, 10 000 шагов) и измеряем \nu по наклону \ln E(t) — критерий приёмки: отклонение <6\% (фактически \lesssim1\%, см. demos_gpu/out/tgv.txt). \langle H\rangle(t) монотонно не возрастает — H-теорема выполняется.

Скрипт: demos_gpu/tgv_gpu.py.

8. Демо-2: устойчивость — двойной сдвиговый слой (BGK против KBC)

Классический тест Minion–Brown: тонкие сдвиговые слои на недоразрешённой сетке (N=128, \nu=1.7\cdot10^{-4}, \mathrm{Re}\sim3\cdot10^4). Стандартный BGK (+полиномиальное равновесие) при малой вязкости разрушается (отрицательные популяции → NaN), а KBC за счёт адаптивного \gamma остаётся устойчивым и гладко сворачивает вихри — это и есть implicit LES. На средней панели — поле \gamma-2: стабилизатор отклоняется от BGK-предела именно там, где градиенты резкие.

Скрипт: demos_gpu/shear_gpu.py; момент разрушения BGK — в demos_gpu/out/shear.txt.

9. Демо-3: каверна с движущейся крышкой против эталона Ghia (Re=1000)

Квадратная каверна, верхняя стенка движется со скоростью U. На неподвижных стенках — bounce-back; на крышке — bounce-back с поправкой импульса подвижной стенки -2w_i\rho\,(\mathbf c_i\cdot\mathbf u_w)/c_s^2. Профили скорости по центральным линиям сравниваются с эталоном Ghia, Ghia & Shin (1982). Сетка намеренно грубая (N=110) и bounce-back первого порядка, поэтому совпадение качественное (RMS-отклонение — в demos_gpu/out/cavity.txt); оно улучшается измельчением и ГУ второго порядка.

Скрипт: demos_gpu/cavity_gpu.py; профили сохраняются в demos_gpu/out/cavity.npz.

10. Демо-4: обтекание тел разной формы — дорожки Кармана

Канал с втоком слева (Дирихле + плавный разгон), нуль-градиентным оттоком и стенками. Тела (цилиндр, квадрат, треугольник, профиль под углом атаки) — узловой bounce-back по маске. При \mathrm{Re}\approx150 за телами формируется вихревая дорожка Кармана; KBC держит грубую сетку устойчиво. Колебания u_y в зонде за телом — мерило частоты срыва.

Именно «лесенка» ступенчатой маски — мотивация перехода к SDF + Bouzidi в части IV.

Скрипт: demos_gpu/obstacle_gpu.py.

---

Часть III. Измельчение сетки (AMR)

11. Блочное измельчение: зачем и как

Чистый LBM живёт на равномерной решётке: перенос требует соседа ровно на \mathbf c_i. Неоднородное разрешение делают вложенными равномерными блоками разного шага: крупный вдали, мелкий у тела и в следе. На стыке блоков (\Delta x\to\Delta x/2, \Delta t\to\Delta t/2 — акустическое масштабирование, c_s сохраняется):

• время релаксации: \beta_f = \dfrac{1}{1+r\,(1/\beta_c-1)}, что при r=2 эквивалентно \boxed{\tau_f = 2\tau_c - \tfrac12} (вязкость \nu непрерывна через стык);
• неравновесная часть масштабируется: $f^{neq}_f = \dfrac{\beta_c}{r,\beta_f},f^{neq}_c = \dfrac{\tau_f}{2\tau_c},f^{neq}_c$ (обратно — множитель \tfrac{2\tau_c}{\tau_f});
• ghost-рамка тонкого блока заполняется равновесием по интерполированным \rho,\mathbf u плюс масштабированным f^{neq}; между подшагами тонкого уровня граница линейно интерполируется по времени между состояниями родителя в t и t+\Delta t;
• рестрикция: внутренние узлы родителя обновляются с тонкого блока (каждый $r$-й узел) с обратным масштабом f^{neq}.

✓ Сверено с первоисточниками: формулы \beta_f и множитель \beta_c/(r\beta_f) — Dorschner, Frapolli, Chikatamarla, Karlin, PRE 94, 053311 (2016), ур. (34), (40) (KBC полностью совместим с измельчением — \gamma поузловой на каждом уровне); временна́я интерполяция границы между подшагами — Lagrava, Malaspinas, Latt, Chopard, J. Comput. Phys. 231, 4808 (2012). Реализация — solver_2x_sdf/amr.py (patch/ghost/fill/restrict), R01 = τ_f/(2τ_c) в config.py.

12. Минимальный 2-уровневый AMR на цилиндре

Рабочая связка «крупная сетка + мелкий блок 2× вокруг тела» (то же тело bounce-back на обоих уровнях): на 1 крупный шаг — 2 мелких подшага, ghost из крупного уровня, рестрикция только по жидким узлам. В этой минимальной версии граница блока держится постоянной в течение подшагов (без временно́й интерполяции) — её эффект виден в бенчмарке ниже.

Скрипт: demos_gpu/amr_demo_gpu.py (часть 1) — на amr.patch/ghost/fill/restrict финальной модели.

13. Полный вложенный 2×+4× и численное сравнение

Три уровня (коарс → 2× → 4×), рекурсивный sub-cycling 1:2:4, временна́я интерполяция границ. Контролируемый тест — локализованный резкий вихрь; эталон — равномерная мелкая 4× сетка. Работа меряется в cell-updates (аппаратно-независимо).

Выводы (точные числа — demos_gpu/out/amr_bench.txt):

• вложенность эффективна: полный 2×+4× точнее одиночного 4× при заметно меньшей работе;
• временна́я интерполяция снижает ошибку при той же работе;
• полный AMR даёт точность уровня мелкой сетки за в разы меньшую работу;
• крупная сетка дёшева, но её ошибка принципиально не лечится без измельчения.

14. AMR на большой области (исторические прогоны)

Та же связка на большой области (150×72, длинный прогон) — видны грубая сетка и учащённые зоны: 2× (тело+след) и 4× (у тела):

---

Часть IV. Граница SDF + Bouzidi

15. Каталог граничных условий

KBC задаёт только оператор столкновения; полному решателю нужны ГУ:

• Bounce-back — no-slip стенка «на полпути между узлами»: f_{\bar i}=f_i. Первый порядок; для движущейся стенки — поправка -2w_i\rho(\mathbf c_i\cdot\mathbf u_w)/c_s^2.
• Интерполированный bounce-back (Bouzidi) — для кривых стенок: доля стенки на линке q\in(0,1) из SDF; две ветви: $$ q<\tfrac12:;; f_{\bar i}(\mathbf x_f) = 2q,f_i^{post}(\mathbf x_f) + (1-2q),f_i^{post}(\mathbf x_f-\mathbf c_i);\qquad q\ge\tfrac12:;; f_{\bar i}(\mathbf x_f) = \tfrac{1}{2q},f_i^{post}(\mathbf x_f) + \bigl(1-\tfrac{1}{2q}\bigr)f_{\bar i}^{post}(\mathbf x_f). $$ ✓ Bouzidi, Firdaouss, Lallemand, Phys. Fluids 13, 3452 (2001); реализация — solver_2x_sdf/boundary.py (build_bc/apply_bc).
• Zou–He — заданные скорость/давление на входе/выходе (восстановление недостающих популяций из \rho,\mathbf u + bounce-back неравновесной части; второй порядок). ✓ Zou & He, Phys. Fluids 9, 1591 (1997); реализация — channel_bc. На давление-выходе поперечная скорость задаётся режимом outlet_uy: классическое жёсткое u_y=0 либо нуль-градиент (мягкий выход — не отражает вихри дорожки; см. факторное исследование, раздел 24).
• Free-slip (specular) — скользящая стенка: касательный импульс сохраняется, нормальный отражается; массу не дрейфует. Реализация — free_slip_walls.
• Sponge / absorbing-слой — плавный рост вязкости перед выходом; гасит вихри и акустику до отражающей границы. В финальной модели добавлен по итогам факторного исследования (раздел 25): без него длинные домены с парой «вход-скорость / выход-давление» работают как недодемпфированный акустический резонатор. Использовался и в прежней 3D-реализации.
• Grad / регуляризованные — реконструкция f^{neq} через \Pi^{neq}; в 2D-модели не понадобились.

16. Зачем SDF: лесенка против гладкой границы

Ступенчатая маска (часть II, демо-4) приближает окружность «лесенкой»: положение стенки скачет на O(\Delta x), что шумит в силе и завышает C_d. SDF (signed distance field, \varphi<0 в теле) даёт точную долю пересечения каждого линка q=\varphi_f/(\varphi_f-\varphi_s), и Bouzidi восстанавливает отражение в правильной точке — граница эффективно гладкая при том же разрешении.

Сравнение реализаций без SDF (bounce-back) и с SDF+Bouzidi на одинаковой сетке:

---

Часть V. Эволюция GPU-макета

17. Поколения

Поколение	Файлы	Что нового	Известные проблемы (чинились в следующем)
0. CPU-прототипы	_amr_anims.py, _amr_anims_sdf.py	первая полная связка: KBC + AMR (none/2×/2×+4×) на цилиндре; ветка с SDF+Bouzidi	numpy-CPU медленный; сила в BB-ветке считалась из post-collision
1. Базовый GPU	amr_gpu_core.py + amr_cylinder_gpu.py, amr_cylinder_sdf_gpu.py	перенос на CuPy (~100× быстрее), оба типа границы в одном ядре	сила force_on: оба члена из post-collision → Cd отрицательный
2. Оптимизированный GPU	amr_opt/*	фикс силы (второй член из post-streaming); feq без pow; cp.fuse; CUDA-графы с рантайм-самопроверкой; GPU-only	монолитные скрипты, параметры захардкожены
3. Финальная модель	solver_2x_sdf/*	модульная архитектура; Zou-He вход + давление-выход; free-slip стенки; GMEM-сила с торком; кросс-чек ∮σ·n ds; серия по блокировке	—

18. Хронология багов и фиксов (чему научил макет)

1. Сила: второй член — из поля ПОСЛЕ стриминга. В обмене импульсом F=\sum c_i(f_i^{post}+f_{\bar i}) популяция f_{\bar i} — та, что вернулась в жидкий узел после отражения. Если взять её из post-collision поля (как в поколении 1), ведущий симметричный член \sim2\rho w_i сокращается: C_d выходил отрицательным и схлопывался к нулю при измельчении. Урок: считывание силы валидируется так же строго, как сам решатель.
2. Дрейф массы из-за ГУ выхода. Грубый zero-gradient выход не выпускал поток: масса и противодавление копились (\langle\rho\rangle\uparrow), поток глох, C_d и rms C_l уползали к нулю за 100k шагов. Вскрыто диагностикой по окнам времени (тренды \langle\rho\rangle, C_d/\rho). Фикс: Zou-He скоростной вход + давление-выход — масса заякорена (\langle\rho\rangle\approx1.000).
3. Free-slip стенки канала (specular) вместо no-slip: нет паразитного пограничного слоя на искусственных стенках, блокировка уменьшается; цилиндр остаётся no-slip (Bouzidi).
4. Стартовое возмущение — короткий поперечный sin-импульс входа после разгона: вихревая дорожка запускается за единицы периодов вместо тысяч шагов симметричного «застоя».
5. Зонд скорости — на тонком уровне L1 (меньше численного размытия вихрей) и оценка St с параболической интерполяцией пика FFT (иначе St «застревает» на шаге сетки частот).
6. CUDA-графы: захват шага требовал убрать все host→device на горячем пути — cupy.stack заменён преаллокацией (backend.pack), проектор KBC переписан без cuBLAS (gemm тащит alpha/beta с хоста), границы clip запечены литералами в ядро. На шаге t=1 — рантайм-самопроверка граф vs eager с автооткатом.
7. Осреднение силы по подшагам L1 и момент $C_m$ — финальная ревизия считывания (часть VI): сила меряется на каждом подшаге тонкого уровня, торк — с точкой приложения на пересечении линка со стенкой.

---

Часть VI. Финальная модель solver_2x_sdf

19. Архитектура

Обтекание цилиндра: грубый уровень L0 (174×90, D=16, \mathrm{Re}=150, U=0.07) + вложенный патч L1 (×2 — тело и ближний след). Граница цилиндра — SDF+Bouzidi (no-slip); вход/выход — Zou-He (скорость/давление); верх/низ — free-slip. Только GPU/CuPy, шаг целиком захватывается в CUDA-граф.

Модуль	Ответственность
config.py	параметры + производные (\tau, \beta, патч из N_y); env-переопределения
lattice.py	D2Q9 + KBC-проектор P_s
equilibrium.py	product-form f^{eq}, макромоменты
collision.py	KBC-N1 (единственный оператор)
streaming.py	pull-перенос
geometry.py	маски/SDF для L0 и L1
boundary.py	Bouzidi + Zou-He + free-slip
forces.py	GMEM-сила и момент (считывание №1)
forces_stress.py	∮σ·n ds по контуру (считывание №2, кросс-чек)
amr.py	связка L0↔L1 (часть III)
solver.py	оркестратор: persistent-буферы, CUDA-граф с самопроверкой
diagnostics.py	St (параболический пик), таблицы, кросс-чек, фиты серии по \beta_{бл}
run.py / run_blockage.py / run_factors.py	одиночный прогон / серия по блокировке / факторное исследование (реестр out/factors.csv)

20. Считывание силы №1: галилей-инвариантный обмен импульсом (GMEM)

Для каждого линка жидкость→тело (Bouzidi-маска): $$ \Delta\mathbf F = (\mathbf c_i - \mathbf u_w),f_i^{post}(\mathbf x_f)

• (\mathbf c_{\bar i} - \mathbf u_w),f_{\bar i}^{stream}(\mathbf x_f), $$ где f_i^{post} — популяция, уходящая в стенку (после столкновения), f_{\bar i}^{stream} — вернувшаяся (после стриминга/Bouzidi), \mathbf u_w — скорость стенки (для неподвижного цилиндра 0, но члены с \mathbf u_w делают оценку галилей-инвариантной и готовой к подвижным телам — цель SimV4). ✓ Wen, Zhang, Tu, Wang, Fang, J. Comput. Phys. 266, 161 (2014).

Момент (торк): T_z=\sum_{links} (\mathbf r_w\times\Delta\mathbf F)_z с точкой приложения на пересечении линка со стенкой \mathbf r_w=\mathbf r_f+q\,\mathbf c_i (q уже лежит в структуре Bouzidi). C_m = 2T_z/(U^2D_{L1}^2); для кругового цилиндра \langle C_m\rangle\approx0 — встроенный контроль симметрии считывания.

Сила меряется на каждом подшаге L1 и усредняется (анти-алиасинг быстрых мод).

21. Считывание силы №2: баланс импульса контрольного объёма (кросс-чек)

Независимая оценка по замкнутому контуру (окружность R+\delta, \delta=2 тонкие ячейки — вне зоны Bouzidi-реконструкции). Полный баланс импульса контрольного объёма: $$ \mathbf F_{тела} = \oint \bigl[\sigma\cdot\mathbf n

• \rho,\mathbf u,(\mathbf u\cdot\mathbf n)\bigr], ds ;-;\underbrace{\frac{d}{dt}!\int_{кольцо}\rho\mathbf u,dV}{\to,0\ в\ среднем},\qquad \sigma{\alpha\beta} = -p',\delta_{\alpha\beta}
• \Bigl(1-\frac{1}{2\tau_1}\Bigr),\Pi^{neq,dev}_{\alpha\beta}. $$ Тонкости:
• давление без константы (\oint p_0\mathbf n\,ds\equiv0);
• строго девиаторная проекция \Pi^{neq}=\sum_i \mathbf c_i\mathbf c_i(f_i-f_i^{eq}) — в KBC след релаксирует со стабилизатором \gamma\beta, а не 1/\tau, тогда как сдвиг-моменты — с точным 2\beta=1/\tau;
• конвективный член -\rho\mathbf u(\mathbf u\cdot\mathbf n) зануляется только при \delta\to0 (no-slip); в первой версии кросс-чека он был опущен, что давало систематический сдвиг +4\% по \langle C_d\rangle на всех сетках (см. таблицы разделов 22–23) — диагностическая ценность: смещение одинаково при любом N_y, т.е. это свойство эстиматора, а не течения. Добавлен в forces_stress.py;
• сравнение — по средним (\langle C_d\rangle, rms C_l): мгновенные ряды сдвинуты по фазе инерцией кольца R..R{+}\delta (член d/dt).

Согласие двух независимых методов = считывание силы корректно — расхождение с литературой тогда лежит в физике постановки, что и проверяют эксперименты разделов 23–24.

22. Результаты основного прогона (174×90, 100k шагов, 2026-06-09)

Прогон на GPU-сервере: CuPy float32, CUDA-графы активны (self-check \Delta=0), 271 с на 100k шагов L0 (вместе с L1-подшагами и гифкой).

Валидация (\beta_{бл}=D/N_y=0.178, осреднение по второй половине ряда):

	raw	×(1−β)^k	лит. безгранич. (Re≈150)
St	0.221	0.181	0.183 ✓ (1%)
C_d	1.799	1.216	1.330
rms C_l	0.7051	0.4767	~0.3
\langle C_m\rangle	−0.00003		0 ✓

Кросс-чек считывания (версия ещё без конвективного члена — см. раздел 21):

метод	\langle C_d\rangle	rms C_l	\langle C_m\rangle
GMEM (обмен импульсом)	1.799	0.7051	−0.00003
∮σ·n ds (контур R+δ)	1.872	0.7252	−0.00003
расхождение	4.1%	2.8%	—

Диагностика сходимости: \langle\rho\rangle = 1.001 стабильна (масса заякорена); оконные средние C_d в окнах 6–10: 1.797, 1.813, 1.795, 1.797, 1.791 — плато; rms C_l насыщен.

Выводы по считыванию (всё подтверждается):

1. Два независимых метода согласуются в пределах 4% (систематика +4% объяснена опущенным конвективным членом и закрыта — раздел 21).
2. \langle C_m\rangle\approx-3\cdot10^{-5} при масштабе C_l\sim0.7 — симметрия в норме.
3. C_d=1.799 воспроизводит предыдущую итерацию бит-в-бит: осреднение силы по подшагам L1 не сместило среднее (регрессии нет, ряды стали глаже).
4. Следовательно C_d\approx1.8 — истинная сила в смоделированной постановке; вопрос к расхождению с литературой — вопрос к постановке, а не к считыванию.

23. Эксперимент №0: серия по боковой блокировке — гипотеза ОПРОВЕРГНУТА

Гипотеза: завышение $C_d$/rms C_l объясняется боковым стеснением канала; тогда при D=16 фикс. и N_y\uparrow величины должны монотонно стремиться к безграничной литературе, а экстраполяция \beta\to0 — попасть в C_d\approx1.33, rms C_l\approx0.3.

Постановка: N_y\in\{90,128,180\} → \beta\in\{0.178,0.125,0.089\}; центр и патч выводятся из N_y; всё остальное фиксировано. 3×100k шагов (по ~282 с).

Результат (python run_blockage.py, 2026-06-09):

N_y	β	St	St(1−β)	\langle C_d\rangle	rms C_l	\langle C_m\rangle	C_d σ·n	dCd
90	0.178	0.221	0.181	1.799	0.7051	−0.00003	1.872	4.1%
128	0.125	0.196	0.171	1.771	0.6734	−0.00007	1.838	3.8%
180	0.089	0.193	0.176	1.870	0.8354	+0.00043	1.937	3.6%

Экстраполяция \beta\to0 (интерсепт линейного / квадратичного фита, спред = неопределённость):

	лин.	квадр.	спред	лит.
C_d(0)	1.905	1.855	0.050	1.33
rms C_l(0)	0.910	0.818	0.091	~0.3
St(0)	0.161	0.181	0.020	0.183

Анализ:

1. C_d(\beta) немонотонен (1.80 → 1.77 → 1.87), rms C_l при самом широком канале даже вырос — экстраполяция «$C_d(0)\approx1.9$» физически абсурдна и означает одно: β — не управляющий параметр. Поправка (1-\beta)^2 из предыдущей итерации была совпадением на одной точке (N_y=90).
2. Это согласуется с классической теорией продувок: поправки Аллена–Винченти (solid + wake blockage) при \beta=0.178 дают завышение лишь ~15–18%, при \beta=0.089 — ~10%. Боковое стеснение объясняло максимум половину превышения даже на исходной сетке.
3. St — единственная величина с правильным поведением: raw 0.221→0.196→0.193 монотонно стремится к 0.183 сверху; квадратичная экстраполяция даёт 0.181 ≈ лит. 0.183. Кинематика вихреобразования (частота, эффективное Re) верна; завышены именно силовые (давленческие) амплитуды.
4. Улика — поведение $N_y=180$: низкочастотная модуляция с размахом оконных C_d 1.68–1.98 и rms C_l 0.58–0.96 (период ~20–30k шагов; на 50k установившегося режима — всего 2–3 цикла, статистика среднего не набрана; \langle C_m\rangle на порядок больше, чем на узких каналах). Чем шире канал, тем свободнее меандрирует дорожка — и тем сильнее она бьётся о жёсткое условие u_y=0 на выходе, порождая отражения и обратную связь.
5. Диагноз: доминируют продольные ограничения, которые серия держала константой: вход 2.5D до тела (жёсткий равномерный профиль), выход 8.3D (Zou-He давление + u_y=0 — отражает вихри), стык L1→L0 на 6.25D. Их изолирует эксперимент №1.

24. Эксперимент №1: факторное исследование продольных границ

Метод: каждый подозреваемый фактор варьируется отдельно при фиксированной боковой блокировке (N_y=90, \beta=0.178) и 100k шагах — изоляция вкладов входа, выхода и выходного ГУ. Реализация: config.py выводит x-границы патча из cx (патч едет за телом), boundary.channel_bc получил режим outlet_uy:

• "zero" — классический Zou-He: u_y=0 жёстко (как во всех прогонах выше);
• "extrapolate" — u_y нуль-градиент (берётся из столбца x{=}-2), якорь \rho=1 сохраняется; алгебраически: f_6 \mathrel{+}= \tfrac12\rho u_y, f_7 \mathrel{-}= \tfrac12\rho u_y при тех же f_3 и u_x (балансы массы и импульса выполняются точно).

Фактор	Nx	cx	вход	выход	$u_y$-выход	что изолирует
B0_baseline	174	40	2.5D	8.3D	zero	контроль (+ кросс-чек с конвективным членом)
F1_outlet_far	300	40	2.5D	16.2D	zero	близость выхода
F2_inlet_far	230	96	6.0D	8.3D	zero	близость входа
F3_inout_far	390	96	6.0D	18.3D	zero	оба продольных расстояния
F4_outlet_uy	174	40	2.5D	8.3D	extrapolate	кинематическое отражение вихрей ГУ
F5_combo	390	96	6.0D	18.3D	extrapolate	«чистая» постановка

Метрики каждого прогона: St, \langle C_d\rangle, rms C_l, \langle C_m\rangle, кросс-чек dCd (уже с конвективным членом), std оконных средних (мера НЧ-модуляции), \langle\rho\rangle. Каждый прогон фиксируется строкой в накопительном реестре solver_2x_sdf/out/factors.csv (полная конфигурация + результаты + время) — документация всех изменений постановки.

Результаты (python run_factors.py, 2026-06-10, 6×100k шагов по ~280 с):

Фактор	вход	выход	$u_y$-выход	St(1−β)	\langle C_d\rangle	rms C_l	dCd кросс	\langle\rho\rangle_{фин}	модуляция	статус
B0_baseline	2.5D	8.3D	zero	0.181	1.799	0.705	1.0%	1.000	0.024	ok
F1_outlet_far	2.5D	16.2D	zero	0.163	2.838	6.794	1.6%	1.663	0.225	режим смещён
F2_inlet_far	6.0D	8.3D	zero	0.138	2.229	6.120	9.7%	1.655	0.306	режим смещён (срыв на ~75k)
F3_inout_far	6.0D	18.3D	zero	0.160	2.577	5.348	2.0%	1.660	0.269	режим смещён (срыв на ~20k)
F4_outlet_uy	2.5D	8.3D	extrapolate	0.183	1.754	0.557	0.6%	1.0000	0.023	ok ✓
F5_combo	6.0D	18.3D	extrapolate	—	—	—	—	NaN	—	взрыв на 4.5k

Анализ:

1. Конвективный член подтверждён: кросс-чек B0 сжался с 4.1% до 1.0% (F4 — 0.6%) при неизменной физике. Считывание силы теперь согласовано двумя методами на уровне ~1%.
2. F4 (мягкий $u_y$-выход) — выигрыш по всем метрикам: rms C_l 0.705→0.557 (−21%), \langle C_d\rangle 1.799→1.754, \langle\rho\rangle=1.0000 (точнее базы), нулевая НЧ-модуляция (окна 6–10: C_d 1.753–1.755, идеальное плато) и St·(1−β)=0.1828 ≈ лит. 0.183 точно. Жёсткое u_y=0 действительно кинематически отражало вихри.
3. Главное открытие — длинные домены без демпфера невалидны. Все три (F1/F2/F3, u_y=0) уходят на смещённую ветвь \langle\rho\rangle\approx1.66 (F1 и F3 — за ~20k шагов, F2 — на ~75k), где C_d, C_l бессмысленны для сравнения (у F2 разваливается и кросс-чек: 9.7% — поле у контура нестационарно «звенит»). Механизм: канал «вход-скорость + выход-давление» — недодемпфированный акустический резонатор; затухание продольной моды \propto\nu(\pi/N_x)^2 падает в ~5 раз при N_x 174→390. Стартовый транзиент раскачивает моду, вход с фиксированной скоростью качает массу \propto\rho_{in} (положительная обратная связь) — система садится на новую ветвь равновесия потоков массы. F5 — чистая картина роста: \langle\rho\rangle по окнам осциллирует 1.04→0.93→1.05 с периодом \approx 2N_x/c_s\approx1350 шагов (точно акустический период туда-обратно) → NaN на 4.5k.
4. При этом rms u_y зонда нормальный (0.033–0.046) во всех прогонах — гидродинамика дорожки жива; ломаются акустика и массовый баланс. Это ретроспективно объясняет и «загадку» эксперимента №0: N_y=180 был зачатком той же моды (НЧ-модуляция), но поперечное расширение менее опасно продольного.

Урок: продольное удлинение требует поглощающего слоя — стандартное лекарство из каталога ГУ (раздел 15), упоминавшееся ещё в прежней 3D-реализации.

25. Эксперимент №2: губка (absorbing layer) перед выходом

Реализация: в последних sponge_len столбцах L0 вязкость плавно (smoothstep) растёт до $\nu\cdot$sponge_nu_mult (по умолчанию ×30 — локальное Re падает до ~5, вихри и акустика диссипируют до прихода на выход). Технически \beta на L0 становится полем \beta(x) (KBC это допускает: \gamma и так поузловой; CUDA-граф не страдает — поле предвычислено); патч L1 губку не видит (гарантировано assert'ом). Ручки: AMR_SPONGE_LEN, AMR_SPONGE_NU.

Фактор	геометрия	$u_y$-выход	губка	что проверяет
F6_sponge	174×90 (база)	extrapolate	32	валидация губки против F4 (не портит ли короткий домен)
F7_long_sponge	390×90, вход 6D / выход 18.3D	extrapolate	32	целевая «чистая» постановка
F8_long_sp_uy0	390×90	zero	32	изоляция: достаточно ли губки без мягкого u_y

⏳ Числа появятся после python run_factors.py с факторами F6–F8 (3×100k ≈ 15 мин).

Критерии: у всех \langle\rho\rangle\approx1.000 и нет НЧ-модуляции; F7 — ожидание $\langle C_d\rangle\approx1.5$–1.6 (остаётся боковая блокировка β=0.178: по Аллену–Винченти +15–18% к безграничному 1.33) и rms C_l<0.45. Если критерии выполнены — эксперимент №3: повторная серия по β (N_y=90/128/180) на конфигурации F7 → финальная экстраполяция к безграничному пределу. Остаточное расхождение после неё — фактор разрешения (D-рефайнмент).

---

Часть VII. Перенос на D3Q27 и связь с целевой реализацией

Вся математика KBC переносится в 3D без изменений по сути: решётка D3Q27 (веса 8/27, 2/27, 1/54, 1/216 — изотропия 4-го порядка, предпочтительна для KBC/турбулентности), моментный базис 27×27, сдвиговая часть s — девиаторный тензор (5 независимых компонент). Формула \gamma и оператор f\leftarrow f-\beta(2\Delta s+\gamma\Delta h) идентичны.

Чем настоящий KBC отличается от прежнего подхода в 3D-движке (poly eq + TRT + α-лимитер + Smagorinsky):

Аспект	Настоящий KBC	Прежний подход
Равновесие	энтропийное product-form	полиномиальное O(u^2)
Столкновение	энтропийный MRT, k/s/h	TRT (симм./антисимм.)
Стабилизация	\gamma из энтропии	α-лимитер положительности
H-теорема	встроена в вывод \gamma	контролируется отдельно
Турбулентность	implicit LES (параметр-free)	явный Smagorinsky
Решётка	D3Q27	D3Q27 (совпадает)

Дорожная карта переноса (всё отработано на 2D-макете): (1) product-form f^{eq}; (2) проекторы P_s для D3Q27 (предвычислить); (3) per-node \gamma; (4) столкновение; (5) SDF+Bouzidi и GMEM-сила с \mathbf u_w на линке (подвижные тела) + кросс-чек ∮σ·n dA по SDF-изоповерхности; (6) AMR-связка \tau_f=2\tau_c-\tfrac12 с временно́й интерполяцией.

---

Часть VIII. Сверка формул с первоисточниками

Все формулы, ранее помеченные «⚠ сверить с оригиналом», проверены по первоисточникам (плюс численная самопроверка каждой в demos_gpu/checks_gpu.py):

Формула	Источник	Статус
Product-form f^{eq}: A(u)=2-\sqrt{1+3u^2}, B(u)=\frac{2u+\sqrt{1+3u^2}}{1-u}	Bösch/Chikatamarla/Karlin, PRE 92, 043309 (2015), ур. (7)–(11); Ansumali/Karlin/Öttinger, EPL 63, 798 (2003)	✓ дословно
f' = f-\beta(2\Delta s+\gamma\Delta h), \beta=1/(2\tau), mirror-форма	Karlin/Bösch/Chikatamarla, PRE 90, 031302(R) (2014)	✓ дословно
$\gamma^*=\frac1\beta-(2-\frac1\beta)\frac{\langle\Delta s	\Delta h\rangle}{\langle\Delta h	\Delta h\rangle}$, $\langle X
Состав s=\{N,\Pi_{xy}\} (девиатор; «KBC-d»/N1)	PRE 92 (семейство s\in\{d,d{+}t,d{+}q,d{+}t{+}q\}); Dorschner et al., JFM 801, 623 (2016)	✓
AMR: \beta_f=\frac{1}{1+r(1/\beta_c-1)} ⇔ \tau_f=2\tau_c-\frac12; f^{neq}\cdot\frac{\beta_c}{r\beta_f}=\frac{\tau_f}{2\tau_c}f^{neq}	Dorschner/Frapolli/Chikatamarla/Karlin, PRE 94, 053311 (2016), ур. (34), (40)	✓ дословно
Временна́я интерполяция ghost-границы	Lagrava/Malaspinas/Latt/Chopard, JCP 231, 4808 (2012)	✓
Bouzidi: ветви q<\frac12 / q\ge\frac12	Bouzidi/Firdaouss/Lallemand, Phys. Fluids 13, 3452 (2001)	✓
Zou-He: скоростной вход / давление-выход	Zou & He, Phys. Fluids 9, 1591 (1997)	✓
GMEM: (\mathbf c_i-\mathbf u_w)f_i^{post}-(\mathbf c_{\bar i}-\mathbf u_w)f_{\bar i}^{stream}	Wen/Zhang/Tu/Wang/Fang, JCP 266, 161 (2014)	✓
\sigma^v=-(1-\frac{1}{2\tau})\Pi^{neq} (Чепмен–Энског)	стандарт (Krüger et al., The Lattice Boltzmann Method, 2017)	✓

Литература:

• I. V. Karlin, F. Bösch, S. S. Chikatamarla, Gibbs' principle for the lattice-kinetic theory of fluid dynamics, Phys. Rev. E 90, 031302(R) (2014).
• F. Bösch, S. S. Chikatamarla, I. V. Karlin, Entropic multirelaxation lattice Boltzmann models for turbulent flows, Phys. Rev. E 92, 043309 (2015) (arXiv:1507.02518).
• S. Ansumali, I. V. Karlin, H. C. Öttinger, Minimal entropic kinetic models for hydrodynamics, Europhys. Lett. 63, 798 (2003).
• B. Dorschner, F. Bösch, S. S. Chikatamarla, K. Boulouchos, I. V. Karlin, Entropic multi-relaxation time lattice Boltzmann model for complex flows, J. Fluid Mech. 801, 623–651 (2016).
• B. Dorschner, N. Frapolli, S. S. Chikatamarla, I. V. Karlin, Grid refinement for entropic lattice Boltzmann models, Phys. Rev. E 94, 053311 (2016) (arXiv:1608.06915).
• D. Lagrava, O. Malaspinas, J. Latt, B. Chopard, Advances in multi-domain lattice Boltzmann grid refinement, J. Comput. Phys. 231, 4808 (2012).
• M. Bouzidi, M. Firdaouss, P. Lallemand, Momentum transfer of a Boltzmann-lattice fluid with boundaries, Phys. Fluids 13, 3452 (2001).
• Q. Zou, X. He, On pressure and velocity boundary conditions for the lattice Boltzmann BGK model, Phys. Fluids 9, 1591 (1997).
• B. Wen, C. Zhang, Y. Tu, C. Wang, H. Fang, Galilean invariant fluid–solid interfacial dynamics in lattice Boltzmann simulations, J. Comput. Phys. 266, 161 (2014).
• Z. Guo, C. Zheng, B. Shi, Discrete lattice effects on the forcing term in the lattice Boltzmann method, Phys. Rev. E 65, 046308 (2002).
• U. Ghia, K. N. Ghia, C. T. Shin, High-Re solutions for incompressible flow using the Navier-Stokes equations and a multigrid method, J. Comput. Phys. 48, 387 (1982).

---

Часть IX. Итоги

• KBC = энтропийный MRT-LBM: сохраняющиеся моды не трогаем, сдвиг релаксируем фиксированным 2\beta (вязкость), высшие моды — адаптивным \gamma из условия неувеличения энтропии. При \gamma=2 — в точности BGK. Все формулы сверены с первоисточниками и самопроверены численно.
• Бенчмарки (на компонентах финальной модели): вязкость по TGV (<6%, факт ≲1%), H-теорема, устойчивость на сдвиговом слое (BGK разрушается — KBC держит), каверна ≈ Ghia, дорожки Кармана на четырёх телах.
• Измельчение сетки: связка \tau_f=2\tau_c-\tfrac12 + масштабирование f^{neq} + временна́я интерполяция ghost; вложенный 2×+4× даёт точность мелкой сетки за в разы меньшую работу.
• Финальная модель solver_2x_sdf: KBC-N1 + SDF/Bouzidi + Zou-He + free-slip + CUDA-графы. Главные уроки эволюции: второй член силы — из post-streaming; дрейф массы лечится давлением-выходом; зонд и сила — на тонком уровне; графы требуют чистоты горячего пути.
• Считывание сил ДОКАЗАНО корректным (разделы 22, 24): два независимых метода — GMEM с торком и баланс импульса ∮[σ·n−ρu(u·n)]ds — после добавления конвективного члена согласуются на уровне 0.6–1.0% по \langle C_d\rangle; \langle C_m\rangle\approx0; масса стабильна.
• Исследование факторов постановки (разделы 23–25):
  • эксперимент №0 (серия по β) опроверг гипотезу бокового стеснения: C_d(\beta) немонотонен; при этом St(β→0)=0.181 ≈ лит. 0.183 — кинематика верна;
  • эксперимент №1 (6 факторов): мягкий $u_y$-выход (F4) улучшил всё (rms C_l −21%, St·(1−β)=0.183 точно, идеальное плато) — жёсткое u_y=0 отражало вихри. Длинные домены без демпфера невалидны: пара «вход-скорость / выход-давление» — недодемпфированный акустический резонатор (затухание ∝ν(π/Nx)²); масса накачивается до смещённой ветви ⟨ρ⟩≈1.66 либо взрыв (F5, период осцилляций ⟨ρ⟩ ≈ 2Nx/c_s — прямая улика);
  • эксперимент №2 (в работе): губка перед выходом (β(x)-поле на L0) — F6/F7/F8; затем эксперимент №3 — серия по β на чистой постановке F7 → финальная экстраполяция.
• Реестр всех прогонов — solver_2x_sdf/out/factors.csv: полная конфигурация + метрики каждого изменения постановки (готовый материал для отдельного исследования факторов).
• Дальше: эксперименты №2–3 → финальные числа; при остаточном превышении — фактор разрешения (D-рефайнмент). Затем перенос отработанной схемы на D3Q27/GPU в основной движок (часть VII).