{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "ea656dfb",
   "metadata": {},
   "source": [
    "# Метод KBC (энтропийный Lattice Boltzmann): от формул к валидированной GPU-модели\n",
    "\n",
    "**Karlin–Bösch–Chikatamarla (KBC)** — энтропийная multi-relaxation-time формулировка\n",
    "решёточного метода Больцмана (LBM), дающая **безусловную нелинейную устойчивость** при больших\n",
    "числах Рейнольдса без явной подсеточной модели турбулентности (работает как *implicit LES*).\n",
    "\n",
    "Этот документ — рабочий разбор метода и **история эволюции 2D-макета**: от формул и\n",
    "самопроверок через классические бенчмарки (Тейлора–Грина, сдвиговый слой, каверна),\n",
    "блочное измельчение сетки (AMR) и границу SDF+Bouzidi — к **финальной валидированной\n",
    "GPU-модели** обтекания цилиндра `solver_2x_sdf` (D2Q9, KBC-N1, L0+L1, CUDA-графы) с\n",
    "двумя независимыми способами считывания силы.\n",
    "\n",
    "> **Как устроен документ.** Ноутбук **не исполняет код** — только теория и готовые\n",
    "> артефакты (картинки, гифки, числа). Весь код живёт рядом и запускается на GPU-сервере:\n",
    ">\n",
    "> | Где | Что |\n",
    "> |---|---|\n",
    "> | `demos_gpu/` | GPU-прогоны всех демо и самопроверок (см. `demos_gpu/README.md`) |\n",
    "> | `solver_2x_sdf/` | финальная модель (см. `solver_2x_sdf/README.md`) |\n",
    "> | `figures/`, `anim/`, `*/out/` | артефакты, встроенные ниже |\n",
    ">\n",
    "> Демо гоняются **теми же модулями** (решётка, равновесие, столкновение, перенос,\n",
    "> AMR-связка), что и финальная модель, — каждый бенчмарк ниже одновременно валидирует\n",
    "> компонент финальной модели."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "23513c20",
   "metadata": {},
   "source": [
    "## Сводная таблица обозначений\n",
    "\n",
    "Везде **решёточные единицы**: шаг сетки $\\Delta x = 1$ (lu), шаг по времени $\\Delta t = 1$ (ts).\n",
    "\n",
    "| Символ | Название | Смысл / формула | Диапазон |\n",
    "|---|---|---|---|\n",
    "| $f_i$ | популяция | плотность частиц со скоростью $\\mathbf c_i$ | $f_i>0$ |\n",
    "| $\\mathbf c_i$ | дискретная скорость | направление переноса $i$ | компоненты $\\in\\{-1,0,1\\}$ |\n",
    "| $w_i$ | вес квадратуры | $\\sum_i w_i=1$ | $\\{4/9,1/9,1/36\\}$ (D2Q9) |\n",
    "| $c_s$ | скорость звука решётки | $c_s^2=1/3$ | $c_s=1/\\sqrt3$ |\n",
    "| $\\rho$, $\\mathbf u$ | плотность, скорость | $\\rho=\\sum_i f_i$, $\\rho\\mathbf u=\\sum_i \\mathbf c_i f_i$ | $\\rho\\approx1$, $|\\mathbf u|\\ll c_s$ |\n",
    "| $\\tau$, $\\beta$ | релаксация | $\\nu=c_s^2(\\tau-\\tfrac12)$, $\\beta=1/(2\\tau)$ | $\\tau>1/2$ |\n",
    "| $\\gamma$ | **энтропийный стабилизатор** | релаксация высших мод | $\\approx2$ |\n",
    "| $f_i^{\\mathrm{eq}}$ | равновесие | максимум энтропии при $\\rho,\\mathbf u$ | $>0$ |\n",
    "| $k_i,s_i,h_i$ | разложение | кинематич./сдвиг/высшие моды | $k+s+h=f$ |\n",
    "| $\\Delta s_i,\\Delta h_i$ | неравновесие | $P_s(f{-}f^{eq})$, остальное | малы |\n",
    "| $H$ | H-функция | $H=\\sum_i f_i\\ln(f_i/w_i)$ | не растёт |\n",
    "| $\\langle X|Y\\rangle$ | энтропийное скал. произв. | $\\sum_i X_iY_i/f_i^{eq}$ | — |\n",
    "| $T,N,\\Pi_{xy}$ | натуральные 2-е моменты | $\\Pi_{xx}{+}\\Pi_{yy}$, $\\Pi_{xx}{-}\\Pi_{yy}$, сдвиг | — |\n",
    "| $\\Pi^{neq}_{\\alpha\\beta}$ | неравновесный тензор | $\\sum_i c_{i\\alpha}c_{i\\beta}(f_i-f_i^{eq})$ | — |\n",
    "| $\\sigma_{\\alpha\\beta}$ | тензор напряжений | $-p'\\delta_{\\alpha\\beta}-(1-\\tfrac{1}{2\\tau})\\Pi^{neq,dev}_{\\alpha\\beta}$ | — |\n",
    "| $q$ | доля стенки (Bouzidi) | положение стенки на линке из SDF | $(0,1)$ |\n",
    "| $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}$ | — |\n",
    "| $\\beta_{бл}$ | блокировка канала | $D/N_y$ (стеснение стенками) | $\\to0$ — безграничный |\n",
    "| $\\mathrm{Re}$, $\\mathrm{Ma}$, $\\mathrm{St}$ | подобие | $UL/\\nu$, $|\\mathbf u|/c_s$, $fD/U$ | — |"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "02c48319",
   "metadata": {},
   "source": [
    "---\n",
    "# Часть I. Теория\n",
    "\n",
    "## 1. Основы LBM и решётка D2Q9\n",
    "\n",
    "LBM отслеживает **популяции** $f_i(\\mathbf x, t)$ — плотность фиктивных частиц, движущихся в\n",
    "узле $\\mathbf x$ со скоростью $\\mathbf c_i$. Динамика — чередование **переноса** (streaming) и\n",
    "**столкновения** (collision):\n",
    "$$ f_i(\\mathbf x + \\mathbf c_i\\,\\Delta t,\\; t+\\Delta t) = f_i(\\mathbf x,t) + \\Omega_i. $$\n",
    "\n",
    "**Решётка D2Q9**: покоящаяся частица, 4 осевые и 4 диагональные скорости; $c_s^2=1/3$;\n",
    "веса $w_0=4/9$, осевые $1/9$, диагональные $1/36$ дают изотропию 4-го порядка.\n",
    "Макромоменты: $\\rho=\\sum_i f_i$, $\\rho\\mathbf u=\\sum_i \\mathbf c_i f_i$.\n",
    "\n",
    "![Стенсиль D2Q9](figures/01_stencil_d2q9.png)\n",
    "\n",
    "**Самопроверка** (`demos_gpu/checks_gpu.py`, исполняется на компонентах финальной модели —\n",
    "`solver_2x_sdf/lattice.py`): $\\sum w_i=1$, $\\sum w_i\\mathbf c_i=0$,\n",
    "$\\sum w_i c_{i\\alpha}c_{i\\beta}=c_s^2\\delta_{\\alpha\\beta}$, $\\mathbf c_{\\bar i}=-\\mathbf c_i$ —\n",
    "все тождества проходят (см. `demos_gpu/out/checks.txt`)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8a74a375",
   "metadata": {},
   "source": [
    "## 2. Равновесие: полиномиальное и энтропийное\n",
    "\n",
    "**Полиномиальное** (разложение максвеллиана до $O(u^2)$, стандарт BGK):\n",
    "$$ f_i^{\\mathrm{eq,poly}} = w_i\\,\\rho\\left[1 + \\frac{\\mathbf c_i\\cdot\\mathbf u}{c_s^2}\n",
    "   + \\frac{(\\mathbf c_i\\cdot\\mathbf u)^2}{2c_s^4} - \\frac{\\mathbf u^2}{2c_s^2}\\right]. $$\n",
    "\n",
    "**Энтропийное** (точный максимум дискретной H-функции при связях $\\rho$, $\\rho\\mathbf u$) для\n",
    "D2Q9/D3Q27 имеет замкнутую **product-form**:\n",
    "$$ f_i^{\\mathrm{eq}} = \\rho\\,w_i\\prod_{\\alpha\\in\\{x,y\\}}\n",
    "   A(u_\\alpha)\\,B(u_\\alpha)^{\\,c_{i\\alpha}},\\qquad\n",
    "   A(u) = 2-\\sqrt{1+3u^2},\\qquad B(u) = \\frac{2u+\\sqrt{1+3u^2}}{1-u}. $$\n",
    "\n",
    "✓ **Сверено с первоисточником**: форма и коэффициенты в точности совпадают с уравнениями\n",
    "(7)–(11) Bösch, Chikatamarla, Karlin, *Phys. Rev. E* **92**, 043309 (2015)\n",
    "(восходит к Ansumali, Karlin, Öttinger, *Europhys. Lett.* **63**, 798, 2003).\n",
    "Реализация — `solver_2x_sdf/equilibrium.py` (показатели $c_{i\\alpha}\\in\\{-1,0,1\\}$, поэтому\n",
    "вместо `pow` берётся $B$ или $1/B$ — быстрее и совместимо с CUDA-графами).\n",
    "\n",
    "**Самопроверка** (`checks_gpu.py`): энтропийное равновесие воспроизводит $\\rho$ и\n",
    "$\\rho\\mathbf u$ **точно** (до машинного нуля выбранной точности), а второй момент при малом\n",
    "$\\mathrm{Ma}$ близок к $\\rho(c_s^2\\delta_{\\alpha\\beta}+u_\\alpha u_\\beta)$ — тензору давления\n",
    "Навье–Стокса. В отличие от полинома, product-form гарантирует положительность и H-теорему."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c6abf6ca",
   "metadata": {},
   "source": [
    "## 3. Разложение популяций $f = k + s + h$\n",
    "\n",
    "Сердце KBC — разбиение 9-вектора популяций по **моментным подпространствам**. Натуральные\n",
    "вторые моменты D2Q9: след $T=\\Pi_{xx}{+}\\Pi_{yy}$, разность нормальных напряжений\n",
    "$N=\\Pi_{xx}{-}\\Pi_{yy}$, сдвиг $\\Pi_{xy}$; плюс третьи ($Q_{xxy},Q_{xyy}$) и четвёртый\n",
    "($\\Pi_{xxyy}$) моменты:\n",
    "\n",
    "- $k_i$ — **кинематическая** часть: только сохраняющиеся моменты ($\\rho, j_x, j_y$);\n",
    "- $s_i$ — **сдвиговая** часть: девиаторный тензор напряжений ($N$ и $\\Pi_{xy}$);\n",
    "- $h_i$ — **высшие** («ghost») моды: $T$, третьи и четвёртый моменты.\n",
    "\n",
    "Строится моментная матрица $M_{ai}=\\phi^{(a)}(\\mathbf c_i)$ из 9 независимых мономов; проектор\n",
    "на группу моментов $G$: $P_G = M^{-1} D_G M$. По построению $P_k+P_s+P_h=I$, разложение точное.\n",
    "\n",
    "✓ **Сверено с первоисточником**: назначение $s$ задаёт *вариант* KBC; вариант с\n",
    "$s=\\{N,\\Pi_{xy}\\}$ (только девиатор) — канонический «KBC-d» / **N1**\n",
    "(Karlin, Bösch, Chikatamarla, *PRE* **90**, 031302(R), 2014; Bösch et al., *PRE* **92**,\n",
    "043309, 2015 — там же семейство $s\\in\\{d,\\,d{+}t,\\,d{+}q,\\,d{+}t{+}q\\}$). Именно он\n",
    "реализован в `solver_2x_sdf/lattice.py` (проектор `Ps` на строки $\\{c_x^2{-}c_y^2,\\,c_xc_y\\}$).\n",
    "\n",
    "**Самопроверка** (`checks_gpu.py`): `Ps` идемпотентен, $\\Delta s$ несёт **только** моменты\n",
    "$\\{N,\\Pi_{xy}\\}$ (остальные строки $M\\Delta s$ — нули)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3fe4992c",
   "metadata": {},
   "source": [
    "## 4. Оператор столкновений KBC и стабилизатор $\\gamma$\n",
    "\n",
    "Так как $f$ и $f^{eq}$ имеют одинаковые сохраняющиеся моменты, $\\Delta k=0$ и\n",
    "$f-f^{eq}=\\Delta s+\\Delta h$, где $\\Delta s = P_s(f-f^{eq})$.\n",
    "\n",
    "**Столкновение KBC** (через «зеркальное состояние» $f^{mirr}$):\n",
    "$$ f_i' = (1-\\beta)\\,f_i + \\beta\\,f_i^{mirr}\n",
    "   \\;\\;\\Longleftrightarrow\\;\\;\n",
    "   \\boxed{\\,f_i' = f_i - \\beta\\,\\bigl(2\\,\\Delta s_i + \\gamma\\,\\Delta h_i\\bigr)\\,},\\qquad\n",
    "   \\nu = c_s^2\\Bigl(\\frac{1}{2\\beta}-\\frac12\\Bigr)\\;\\Rightarrow\\;\\beta=\\frac{1}{2\\tau}. $$\n",
    "Множитель **2** перед $\\Delta s$ фиксирует сдвиговую вязкость; высшие моды релаксируются\n",
    "отдельным $\\gamma$.\n",
    "\n",
    "**Стабилизатор** $\\gamma$ выбирается локально (в каждом узле на каждом шаге) из условия\n",
    "неувеличения дискретной энтропии; аналитическое приближение корня:\n",
    "$$ \\gamma^\\* = \\frac{1}{\\beta} - \\Bigl(2-\\frac{1}{\\beta}\\Bigr)\n",
    "   \\frac{\\langle \\Delta s\\,|\\,\\Delta h\\rangle}{\\langle \\Delta h\\,|\\,\\Delta h\\rangle},\n",
    "   \\qquad \\langle X|Y\\rangle = \\sum_i \\frac{X_i Y_i}{f_i^{eq}}. $$\n",
    "\n",
    "✓ **Сверено с первоисточником**: обе формулы (включая зеркальную форму и связь\n",
    "$\\nu(\\beta)$) дословно совпадают с Karlin, Bösch, Chikatamarla, *PRE* **90**, 031302(R)\n",
    "(2014) и Bösch et al., *PRE* **92**, 043309 (2015), ур. (24)–(25). Реализация —\n",
    "`solver_2x_sdf/collision.py` (с регуляризацией $\\gamma\\to2$ при\n",
    "$\\langle\\Delta h|\\Delta h\\rangle\\to0$).\n",
    "\n",
    "При $\\gamma=2$: $2\\Delta s+2\\Delta h=2(f-f^{eq})$ — KBC **сводится к BGK** с $\\omega=2\\beta$.\n",
    "**Самопроверка** (`checks_gpu.py`): KBC$(\\gamma{=}2)\\equiv$ BGK до машинного нуля; столкновение\n",
    "сохраняет $\\rho$ и $\\rho\\mathbf u$.\n",
    "\n",
    "**Популярно.** Разложение $k+s+h$ — три «слоя» движения: перенос массы/импульса, вязкий сдвиг\n",
    "(трение слоёв — задаёт вязкость и Re) и мелкая «рябь» высших мод, где на грубой сетке копится\n",
    "неустойчивость. KBC гасит сдвиг фиксированным темпом $2\\beta$ (физика), а рябь — адаптивным\n",
    "$\\gamma\\beta$: в спокойных зонах $\\gamma\\approx2$ (как BGK), у резких градиентов $\\gamma$\n",
    "отклоняется и добавляет ровно столько диссипации, сколько разрешает энтропия. Это «локальный\n",
    "термостат потока» — встроенная, без подгоночных констант, модель мелких масштабов\n",
    "(*implicit LES*)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9eba837e",
   "metadata": {},
   "source": [
    "## 5. Перенос и связь вязкости с $\\beta$; H-функция\n",
    "\n",
    "Перенос — сдвиг каждой популяции на её $\\mathbf c_i$ (pull-схема, `solver_2x_sdf/streaming.py`).\n",
    "Вязкость: $\\nu = c_s^2(\\tfrac{1}{2\\beta}-\\tfrac12)$, т.е. $\\tau=\\nu/c_s^2+\\tfrac12$.\n",
    "Число Рейнольдса $\\mathrm{Re}=UL/\\nu$; Маха $\\mathrm{Ma}=U/c_s\\lesssim0.1$ (ошибка\n",
    "сжимаемости $\\propto \\mathrm{Ma}^2$).\n",
    "\n",
    "Дискретная энтропия $H=\\sum_i f_i\\ln(f_i/w_i)$ — «физический компас» устойчивости: из условия\n",
    "«$H$ не растёт» выведен $\\gamma$; в демо Тейлора–Грина ниже $\\langle H\\rangle(t)$ монотонно\n",
    "не возрастает.\n",
    "\n",
    "## 6. Алгоритм одного шага KBC\n",
    "\n",
    "```\n",
    "для каждого узла:\n",
    "  1) макро-моменты:   rho = sum f_i,   rho*u = sum c_i f_i\n",
    "  2) равновесие:      f_i^eq  (энтропийное product-form)\n",
    "  3) неравновесие:    Ds = P_s (f - f^eq),   Dh = (f - f^eq) - Ds\n",
    "  4) стабилизатор:    gamma = 1/beta - (2 - 1/beta) <Ds|Dh> / <Dh|Dh>\n",
    "  5) столкновение:    f <- f - beta (2 Ds + gamma Dh)\n",
    "  6) перенос:         f_i(x + c_i) <- f_i(x)\n",
    "  7) граничные условия (bounce-back / Bouzidi / Zou-He / free-slip)\n",
    "```\n",
    "\n",
    "![Блок-схема шага KBC](figures/06_algorithm_flow.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a27e6548",
   "metadata": {},
   "source": [
    "---\n",
    "# Часть II. Валидационные демо\n",
    "\n",
    "Четыре классических теста; каждый гоняется скриптом из `demos_gpu/` **на тех же модулях, что\n",
    "финальная модель** (решётка/равновесие/KBC/перенос из `solver_2x_sdf`). Числа ниже цитируются\n",
    "из `demos_gpu/out/*.txt` — перегенерируются при каждом прогоне на GPU-сервере.\n",
    "\n",
    "## 7. Демо-1: вихрь Тейлора–Грина — валидация вязкости и H-теоремы\n",
    "\n",
    "Точное решение несжимаемого Навье–Стокса в периодической области:\n",
    "$u_x=-U_0\\cos(kx)\\sin(ky)$, $u_y=U_0\\sin(kx)\\cos(ky)$, $k=2\\pi/N$; кинетическая энергия\n",
    "затухает как $E(t)=E_0\\,e^{-4\\nu k^2 t}$. Прогоняем KBC ($N=96$, $\\nu=0.0125$, 10 000 шагов)\n",
    "и **измеряем** $\\nu$ по наклону $\\ln E(t)$ — критерий приёмки: отклонение $<6\\%$\n",
    "(фактически $\\lesssim1\\%$, см. `demos_gpu/out/tgv.txt`). $\\langle H\\rangle(t)$ монотонно\n",
    "не возрастает — H-теорема выполняется.\n",
    "\n",
    "![TGV: завихренность, затухание энергии, H-функция](figures/07_taylor_green.png)\n",
    "\n",
    "![Распад вихря: 10 срезов на единой шкале цвета](figures/07b_tgv_decay_slices.png)\n",
    "\n",
    "![Анимация затухания](anim/tgv_decay.gif)\n",
    "\n",
    "*Скрипт: `demos_gpu/tgv_gpu.py`.*"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0f68444b",
   "metadata": {},
   "source": [
    "## 8. Демо-2: устойчивость — двойной сдвиговый слой (BGK против KBC)\n",
    "\n",
    "Классический тест Minion–Brown: тонкие сдвиговые слои на **недоразрешённой** сетке\n",
    "($N=128$, $\\nu=1.7\\cdot10^{-4}$, $\\mathrm{Re}\\sim3\\cdot10^4$). Стандартный BGK\n",
    "(+полиномиальное равновесие) при малой вязкости разрушается (отрицательные популяции → NaN),\n",
    "а KBC за счёт адаптивного $\\gamma$ остаётся устойчивым и гладко сворачивает вихри — это и есть\n",
    "*implicit LES*. На средней панели — поле $\\gamma-2$: стабилизатор отклоняется от BGK-предела\n",
    "именно там, где градиенты резкие.\n",
    "\n",
    "![Сдвиговый слой: KBC устойчив, BGK разрушается](figures/08_shear_stability.png)\n",
    "\n",
    "![BGK (слева, замирает в момент разрушения) против KBC (справа)](anim/shear_bgk_vs_kbc.gif)\n",
    "\n",
    "*Скрипт: `demos_gpu/shear_gpu.py`; момент разрушения BGK — в `demos_gpu/out/shear.txt`.*"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d8ebf6c3",
   "metadata": {},
   "source": [
    "## 9. Демо-3: каверна с движущейся крышкой против эталона Ghia (Re=1000)\n",
    "\n",
    "Квадратная каверна, верхняя стенка движется со скоростью $U$. На неподвижных стенках —\n",
    "**bounce-back**; на крышке — bounce-back с поправкой импульса подвижной стенки\n",
    "$-2w_i\\rho\\,(\\mathbf c_i\\cdot\\mathbf u_w)/c_s^2$. Профили скорости по центральным линиям\n",
    "сравниваются с эталоном **Ghia, Ghia & Shin (1982)**. Сетка намеренно грубая ($N=110$) и\n",
    "bounce-back первого порядка, поэтому совпадение **качественное** (RMS-отклонение —\n",
    "в `demos_gpu/out/cavity.txt`); оно улучшается измельчением и ГУ второго порядка.\n",
    "\n",
    "![Каверна: линии тока и сравнение с Ghia](figures/09_cavity_ghia.png)\n",
    "\n",
    "![Становление течения в каверне](anim/cavity_transient.gif)\n",
    "\n",
    "*Скрипт: `demos_gpu/cavity_gpu.py`; профили сохраняются в `demos_gpu/out/cavity.npz`.*"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "761a536f",
   "metadata": {},
   "source": [
    "## 10. Демо-4: обтекание тел разной формы — дорожки Кармана\n",
    "\n",
    "Канал с втоком слева (Дирихле + плавный разгон), нуль-градиентным оттоком и стенками.\n",
    "Тела (цилиндр, квадрат, треугольник, профиль под углом атаки) — **узловой bounce-back** по\n",
    "маске. При $\\mathrm{Re}\\approx150$ за телами формируется **вихревая дорожка Кармана**;\n",
    "KBC держит грубую сетку устойчиво. Колебания $u_y$ в зонде за телом — мерило частоты срыва.\n",
    "\n",
    "Именно «лесенка» ступенчатой маски — мотивация перехода к **SDF + Bouzidi** в части IV.\n",
    "\n",
    "![Четыре тела: поле завихренности в конце прогона](figures/10_obstacle_shapes.png)\n",
    "\n",
    "![Зонд в следе: колебания = срыв вихрей](figures/10b_obstacle_probe.png)\n",
    "\n",
    "![Анимация всех четырёх тел](anim/obstacle_shapes.gif)\n",
    "\n",
    "*Скрипт: `demos_gpu/obstacle_gpu.py`.*"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2056c378",
   "metadata": {},
   "source": [
    "---\n",
    "# Часть III. Измельчение сетки (AMR)\n",
    "\n",
    "## 11. Блочное измельчение: зачем и как\n",
    "\n",
    "Чистый LBM живёт на **равномерной** решётке: перенос требует соседа ровно на $\\mathbf c_i$.\n",
    "Неоднородное разрешение делают **вложенными равномерными блоками** разного шага: крупный\n",
    "вдали, мелкий у тела и в следе. На стыке блоков ($\\Delta x\\to\\Delta x/2$,\n",
    "$\\Delta t\\to\\Delta t/2$ — акустическое масштабирование, $c_s$ сохраняется):\n",
    "\n",
    "- **время релаксации**: $\\beta_f = \\dfrac{1}{1+r\\,(1/\\beta_c-1)}$, что при $r=2$ эквивалентно\n",
    "  $\\boxed{\\tau_f = 2\\tau_c - \\tfrac12}$ (вязкость $\\nu$ непрерывна через стык);\n",
    "- **неравновесная часть** масштабируется: $f^{neq}_f = \\dfrac{\\beta_c}{r\\,\\beta_f}\\,f^{neq}_c\n",
    "  = \\dfrac{\\tau_f}{2\\tau_c}\\,f^{neq}_c$ (обратно — множитель $\\tfrac{2\\tau_c}{\\tau_f}$);\n",
    "- **ghost-рамка** тонкого блока заполняется равновесием по интерполированным $\\rho,\\mathbf u$\n",
    "  плюс масштабированным $f^{neq}$; между подшагами тонкого уровня граница **линейно\n",
    "  интерполируется по времени** между состояниями родителя в $t$ и $t+\\Delta t$;\n",
    "- **рестрикция**: внутренние узлы родителя обновляются с тонкого блока (каждый $r$-й узел)\n",
    "  с обратным масштабом $f^{neq}$.\n",
    "\n",
    "✓ **Сверено с первоисточниками**: формулы $\\beta_f$ и множитель $\\beta_c/(r\\beta_f)$ —\n",
    "Dorschner, Frapolli, Chikatamarla, Karlin, *PRE* **94**, 053311 (2016), ур. (34), (40)\n",
    "(KBC полностью совместим с измельчением — $\\gamma$ поузловой на каждом уровне); временна́я\n",
    "интерполяция границы между подшагами — Lagrava, Malaspinas, Latt, Chopard, *J. Comput. Phys.*\n",
    "**231**, 4808 (2012). Реализация — `solver_2x_sdf/amr.py` (`patch/ghost/fill/restrict`),\n",
    "`R01 = τ_f/(2τ_c)` в `config.py`.\n",
    "\n",
    "![Схема блочного измельчения поверх реального поля |ω|](figures/11c_refinement_schematic.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6645f68a",
   "metadata": {},
   "source": [
    "## 12. Минимальный 2-уровневый AMR на цилиндре\n",
    "\n",
    "Рабочая связка «крупная сетка + мелкий блок 2× вокруг тела» (то же тело bounce-back на обоих\n",
    "уровнях): на 1 крупный шаг — 2 мелких подшага, ghost из крупного уровня, рестрикция только по\n",
    "жидким узлам. В этой *минимальной* версии граница блока держится постоянной в течение\n",
    "подшагов (без временно́й интерполяции) — её эффект виден в бенчмарке ниже.\n",
    "\n",
    "![Минимальный AMR: крупное поле + мелкий блок 2× (жёлтая рамка)](anim/amr_cylinder.gif)\n",
    "\n",
    "*Скрипт: `demos_gpu/amr_demo_gpu.py` (часть 1) — на `amr.patch/ghost/fill/restrict`\n",
    "финальной модели.*\n",
    "\n",
    "## 13. Полный вложенный 2×+4× и численное сравнение\n",
    "\n",
    "Три уровня (коарс → 2× → 4×), рекурсивный sub-cycling 1:2:4, временна́я интерполяция границ.\n",
    "Контролируемый тест — локализованный резкий вихрь; эталон — равномерная мелкая 4× сетка.\n",
    "Работа меряется в **cell-updates** (аппаратно-независимо).\n",
    "\n",
    "![AMR: точность и фронт «точность–работа»](figures/11d_amr_comparison.png)\n",
    "\n",
    "Выводы (точные числа — `demos_gpu/out/amr_bench.txt`):\n",
    "- **вложенность эффективна**: полный 2×+4× точнее одиночного 4× при заметно меньшей работе;\n",
    "- **временна́я интерполяция** снижает ошибку при той же работе;\n",
    "- полный AMR даёт точность уровня мелкой сетки за **в разы меньшую работу**;\n",
    "- крупная сетка дёшева, но её ошибка принципиально не лечится без измельчения.\n",
    "\n",
    "![Вложенный AMR: видны три размера ячеек](anim/amr_nested.gif)\n",
    "\n",
    "## 14. AMR на большой области (исторические прогоны)\n",
    "\n",
    "Та же связка на большой области (150×72, длинный прогон) — видны грубая сетка и учащённые\n",
    "зоны: **2× (тело+след)** и **4× (у тела)**:\n",
    "\n",
    "![Одиночный 2× AMR на цилиндре](anim/amr2x_cyl.gif)\n",
    "\n",
    "![Вложенный 2×+4× AMR на цилиндре](anim/amr_nested_cyl.gif)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ed665c26",
   "metadata": {},
   "source": [
    "---\n",
    "# Часть IV. Граница SDF + Bouzidi\n",
    "\n",
    "## 15. Каталог граничных условий\n",
    "\n",
    "KBC задаёт только оператор столкновения; полному решателю нужны ГУ:\n",
    "\n",
    "- **Bounce-back** — no-slip стенка «на полпути между узлами»: $f_{\\bar i}=f_i$. Первый\n",
    "  порядок; для движущейся стенки — поправка $-2w_i\\rho(\\mathbf c_i\\cdot\\mathbf u_w)/c_s^2$.\n",
    "- **Интерполированный bounce-back (Bouzidi)** — для **кривых** стенок: доля стенки на линке\n",
    "  $q\\in(0,1)$ из SDF; две ветви:\n",
    "  $$ 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\n",
    "     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). $$\n",
    "  ✓ Bouzidi, Firdaouss, Lallemand, *Phys. Fluids* **13**, 3452 (2001); реализация —\n",
    "  `solver_2x_sdf/boundary.py` (`build_bc`/`apply_bc`).\n",
    "- **Zou–He** — заданные скорость/давление на входе/выходе (восстановление недостающих\n",
    "  популяций из $\\rho,\\mathbf u$ + bounce-back неравновесной части; второй порядок).\n",
    "  ✓ Zou & He, *Phys. Fluids* **9**, 1591 (1997); реализация — `channel_bc`. На давление-выходе\n",
    "  поперечная скорость задаётся режимом `outlet_uy`: классическое жёсткое $u_y=0$ либо\n",
    "  нуль-градиент (мягкий выход — не отражает вихри дорожки; см. факторное исследование,\n",
    "  раздел 24).\n",
    "- **Free-slip (specular)** — скользящая стенка: касательный импульс сохраняется, нормальный\n",
    "  отражается; массу не дрейфует. Реализация — `free_slip_walls`.\n",
    "- **Sponge / absorbing-слой** — плавный рост вязкости перед выходом; гасит вихри и акустику до\n",
    "  отражающей границы. В финальной модели **добавлен** по итогам факторного исследования\n",
    "  (раздел 25): без него длинные домены с парой «вход-скорость / выход-давление» работают как\n",
    "  недодемпфированный акустический резонатор. Использовался и в прежней 3D-реализации.\n",
    "- **Grad / регуляризованные** — реконструкция $f^{neq}$ через $\\Pi^{neq}$; в 2D-модели\n",
    "  не понадобились.\n",
    "\n",
    "## 16. Зачем SDF: лесенка против гладкой границы\n",
    "\n",
    "Ступенчатая маска (часть II, демо-4) приближает окружность «лесенкой»: положение стенки\n",
    "скачет на $O(\\Delta x)$, что шумит в силе и завышает $C_d$. **SDF** (signed distance field,\n",
    "$\\varphi<0$ в теле) даёт точную долю пересечения каждого линка $q=\\varphi_f/(\\varphi_f-\\varphi_s)$,\n",
    "и Bouzidi восстанавливает отражение **в правильной точке** — граница эффективно гладкая\n",
    "при том же разрешении.\n",
    "\n",
    "Сравнение реализаций без SDF (bounce-back) и с SDF+Bouzidi на одинаковой сетке:\n",
    "\n",
    "![Зонды: AMR с SDF](figures/11f_amr_sdf_probe_comparison.png)\n",
    "\n",
    "![Сводное сравнение SDF vs ступенчатая граница](figures/11g_sdf_vs_nosdf.png)\n",
    "\n",
    "![Цилиндр 2× AMR + SDF](anim/amr2x_cyl_sdf.gif)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9d235a12",
   "metadata": {},
   "source": [
    "---\n",
    "# Часть V. Эволюция GPU-макета\n",
    "\n",
    "## 17. Поколения\n",
    "\n",
    "| Поколение | Файлы | Что нового | Известные проблемы (чинились в следующем) |\n",
    "|---|---|---|---|\n",
    "| 0. CPU-прототипы | `_amr_anims.py`, `_amr_anims_sdf.py` | первая полная связка: KBC + AMR (none/2×/2×+4×) на цилиндре; ветка с SDF+Bouzidi | numpy-CPU медленный; сила в BB-ветке считалась из post-collision |\n",
    "| 1. Базовый GPU | `amr_gpu_core.py` + `amr_cylinder_gpu.py`, `amr_cylinder_sdf_gpu.py` | перенос на CuPy (~100× быстрее), оба типа границы в одном ядре | **сила `force_on`: оба члена из post-collision → Cd отрицательный** |\n",
    "| 2. Оптимизированный GPU | `amr_opt/*` | **фикс силы** (второй член из post-streaming); `feq` без `pow`; `cp.fuse`; **CUDA-графы** с рантайм-самопроверкой; GPU-only | монолитные скрипты, параметры захардкожены |\n",
    "| 3. Финальная модель | `solver_2x_sdf/*` | модульная архитектура; **Zou-He вход + давление-выход**; **free-slip стенки**; GMEM-сила с торком; **кросс-чек ∮σ·n ds**; серия по блокировке | — |\n",
    "\n",
    "## 18. Хронология багов и фиксов (чему научил макет)\n",
    "\n",
    "1. **Сила: второй член — из поля ПОСЛЕ стриминга.** В обмене импульсом\n",
    "   $F=\\sum c_i(f_i^{post}+f_{\\bar i})$ популяция $f_{\\bar i}$ — та, что *вернулась* в жидкий\n",
    "   узел после отражения. Если взять её из post-collision поля (как в поколении 1), ведущий\n",
    "   симметричный член $\\sim2\\rho w_i$ сокращается: $C_d$ выходил **отрицательным** и\n",
    "   схлопывался к нулю при измельчении. Урок: считывание силы валидируется так же строго,\n",
    "   как сам решатель.\n",
    "2. **Дрейф массы из-за ГУ выхода.** Грубый zero-gradient выход не выпускал поток: масса и\n",
    "   противодавление копились ($\\langle\\rho\\rangle\\uparrow$), поток глох, $C_d$ и rms $C_l$\n",
    "   уползали к нулю за 100k шагов. Вскрыто **диагностикой по окнам времени** (тренды\n",
    "   $\\langle\\rho\\rangle$, $C_d/\\rho$). Фикс: **Zou-He скоростной вход + давление-выход** —\n",
    "   масса заякорена ($\\langle\\rho\\rangle\\approx1.000$).\n",
    "3. **Free-slip стенки канала** (specular) вместо no-slip: нет паразитного пограничного слоя\n",
    "   на искусственных стенках, блокировка уменьшается; цилиндр остаётся no-slip (Bouzidi).\n",
    "4. **Стартовое возмущение** — короткий поперечный sin-импульс входа после разгона: вихревая\n",
    "   дорожка запускается за единицы периодов вместо тысяч шагов симметричного «застоя».\n",
    "5. **Зонд скорости — на тонком уровне L1** (меньше численного размытия вихрей) и **оценка St\n",
    "   с параболической интерполяцией пика FFT** (иначе St «застревает» на шаге сетки частот).\n",
    "6. **CUDA-графы**: захват шага требовал убрать все host→device на горячем пути —\n",
    "   `cupy.stack` заменён преаллокацией (`backend.pack`), проектор KBC переписан без cuBLAS\n",
    "   (gemm тащит alpha/beta с хоста), границы clip запечены литералами в ядро. На шаге $t=1$ —\n",
    "   **рантайм-самопроверка** граф vs eager с автооткатом.\n",
    "7. **Осреднение силы по подшагам L1 и момент $C_m$** — финальная ревизия считывания\n",
    "   (часть VI): сила меряется на каждом подшаге тонкого уровня, торк — с точкой приложения\n",
    "   на пересечении линка со стенкой."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e7b57ba6",
   "metadata": {},
   "source": [
    "---\n",
    "# Часть VI. Финальная модель `solver_2x_sdf`\n",
    "\n",
    "## 19. Архитектура\n",
    "\n",
    "Обтекание цилиндра: грубый уровень **L0** (174×90, $D=16$, $\\mathrm{Re}=150$, $U=0.07$) +\n",
    "вложенный патч **L1** (×2 — тело и ближний след). Граница цилиндра — SDF+Bouzidi (no-slip);\n",
    "вход/выход — Zou-He (скорость/давление); верх/низ — free-slip. Только GPU/CuPy, шаг целиком\n",
    "захватывается в CUDA-граф.\n",
    "\n",
    "| Модуль | Ответственность |\n",
    "|---|---|\n",
    "| `config.py` | параметры + производные ($\\tau$, $\\beta$, патч из $N_y$); env-переопределения |\n",
    "| `lattice.py` | D2Q9 + KBC-проектор $P_s$ |\n",
    "| `equilibrium.py` | product-form $f^{eq}$, макромоменты |\n",
    "| `collision.py` | KBC-N1 (единственный оператор) |\n",
    "| `streaming.py` | pull-перенос |\n",
    "| `geometry.py` | маски/SDF для L0 и L1 |\n",
    "| `boundary.py` | Bouzidi + Zou-He + free-slip |\n",
    "| `forces.py` | GMEM-сила и момент (считывание №1) |\n",
    "| `forces_stress.py` | ∮σ·n ds по контуру (считывание №2, кросс-чек) |\n",
    "| `amr.py` | связка L0↔L1 (часть III) |\n",
    "| `solver.py` | оркестратор: persistent-буферы, CUDA-граф с самопроверкой |\n",
    "| `diagnostics.py` | St (параболический пик), таблицы, кросс-чек, фиты серии по $\\beta_{бл}$ |\n",
    "| `run.py` / `run_blockage.py` / `run_factors.py` | одиночный прогон / серия по блокировке / факторное исследование (реестр `out/factors.csv`) |\n",
    "\n",
    "## 20. Считывание силы №1: галилей-инвариантный обмен импульсом (GMEM)\n",
    "\n",
    "Для каждого линка жидкость→тело (Bouzidi-маска):\n",
    "$$ \\Delta\\mathbf F = (\\mathbf c_i - \\mathbf u_w)\\,f_i^{post}(\\mathbf x_f)\n",
    "   - (\\mathbf c_{\\bar i} - \\mathbf u_w)\\,f_{\\bar i}^{stream}(\\mathbf x_f), $$\n",
    "где $f_i^{post}$ — популяция, уходящая в стенку (после столкновения), $f_{\\bar i}^{stream}$ —\n",
    "вернувшаяся (после стриминга/Bouzidi), $\\mathbf u_w$ — скорость стенки (для неподвижного\n",
    "цилиндра $0$, но члены с $\\mathbf u_w$ делают оценку **галилей-инвариантной** и готовой к\n",
    "подвижным телам — цель SimV4). ✓ Wen, Zhang, Tu, Wang, Fang, *J. Comput. Phys.* **266**,\n",
    "161 (2014).\n",
    "\n",
    "**Момент (торк)**: $T_z=\\sum_{links} (\\mathbf r_w\\times\\Delta\\mathbf F)_z$ с точкой приложения\n",
    "на **пересечении линка со стенкой** $\\mathbf r_w=\\mathbf r_f+q\\,\\mathbf c_i$ ($q$ уже лежит в\n",
    "структуре Bouzidi). $C_m = 2T_z/(U^2D_{L1}^2)$; для кругового цилиндра $\\langle C_m\\rangle\\approx0$ —\n",
    "встроенный **контроль симметрии** считывания.\n",
    "\n",
    "Сила меряется на **каждом** подшаге L1 и усредняется (анти-алиасинг быстрых мод).\n",
    "\n",
    "## 21. Считывание силы №2: баланс импульса контрольного объёма (кросс-чек)\n",
    "\n",
    "Независимая оценка по замкнутому контуру (окружность $R+\\delta$, $\\delta=2$ тонкие ячейки —\n",
    "вне зоны Bouzidi-реконструкции). Полный баланс импульса контрольного объёма:\n",
    "$$ \\mathbf F_{тела} = \\oint \\bigl[\\sigma\\cdot\\mathbf n\n",
    "   - \\rho\\,\\mathbf u\\,(\\mathbf u\\cdot\\mathbf n)\\bigr]\\, ds\n",
    "   \\;-\\;\\underbrace{\\frac{d}{dt}\\!\\int_{кольцо}\\rho\\mathbf u\\,dV}_{\\to\\,0\\ в\\ среднем},\\qquad\n",
    "   \\sigma_{\\alpha\\beta} = -p'\\,\\delta_{\\alpha\\beta}\n",
    "   - \\Bigl(1-\\frac{1}{2\\tau_1}\\Bigr)\\,\\Pi^{neq,dev}_{\\alpha\\beta}. $$\n",
    "Тонкости:\n",
    "- давление без константы ($\\oint p_0\\mathbf n\\,ds\\equiv0$);\n",
    "- строго **девиаторная** проекция $\\Pi^{neq}=\\sum_i \\mathbf c_i\\mathbf c_i(f_i-f_i^{eq})$ —\n",
    "  в KBC след релаксирует со стабилизатором $\\gamma\\beta$, а не $1/\\tau$, тогда как\n",
    "  сдвиг-моменты — с точным $2\\beta=1/\\tau$;\n",
    "- **конвективный член** $-\\rho\\mathbf u(\\mathbf u\\cdot\\mathbf n)$ зануляется только при\n",
    "  $\\delta\\to0$ (no-slip); в первой версии кросс-чека он был опущен, что давало\n",
    "  **систематический сдвиг $+4\\%$ по $\\langle C_d\\rangle$ на всех сетках** (см. таблицы\n",
    "  разделов 22–23) — диагностическая ценность: смещение одинаково при любом $N_y$, т.е. это\n",
    "  свойство эстиматора, а не течения. Добавлен в `forces_stress.py`;\n",
    "- сравнение — по средним ($\\langle C_d\\rangle$, rms $C_l$): мгновенные ряды сдвинуты по фазе\n",
    "  инерцией кольца $R..R{+}\\delta$ (член $d/dt$).\n",
    "\n",
    "**Согласие двух независимых методов = считывание силы корректно** — расхождение с\n",
    "литературой тогда лежит в физике постановки, что и проверяют эксперименты разделов 23–24."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ccac6a82",
   "metadata": {},
   "source": [
    "## 22. Результаты основного прогона (174×90, 100k шагов, 2026-06-09)\n",
    "\n",
    "Прогон на GPU-сервере: CuPy float32, CUDA-графы активны (self-check $\\Delta=0$), 271 с на\n",
    "100k шагов L0 (вместе с L1-подшагами и гифкой).\n",
    "\n",
    "![Финальная модель: |ω| с патчем L1](solver_2x_sdf/out/cyl_2x_sdf.gif)\n",
    "\n",
    "![Ряды Cd/Cl/Cm, кросс-чек, спектр, ⟨ρ⟩(t)](solver_2x_sdf/out/series_2x_sdf.png)\n",
    "\n",
    "**Валидация** ($\\beta_{бл}=D/N_y=0.178$, осреднение по второй половине ряда):\n",
    "\n",
    "| | raw | ×(1−β)$^k$ | лит. безгранич. (Re≈150) |\n",
    "|---|---|---|---|\n",
    "| St | 0.221 | **0.181** | **0.183** ✓ (1%) |\n",
    "| $C_d$ | 1.799 | 1.216 | 1.330 |\n",
    "| rms $C_l$ | 0.7051 | 0.4767 | ~0.3 |\n",
    "| $\\langle C_m\\rangle$ | **−0.00003** | | 0 ✓ |\n",
    "\n",
    "**Кросс-чек считывания** (версия ещё без конвективного члена — см. раздел 21):\n",
    "\n",
    "| метод | $\\langle C_d\\rangle$ | rms $C_l$ | $\\langle C_m\\rangle$ |\n",
    "|---|---|---|---|\n",
    "| GMEM (обмен импульсом) | 1.799 | 0.7051 | −0.00003 |\n",
    "| ∮σ·n ds (контур R+δ) | 1.872 | 0.7252 | −0.00003 |\n",
    "| расхождение | **4.1%** | 2.8% | — |\n",
    "\n",
    "Диагностика сходимости: $\\langle\\rho\\rangle = 1.001$ стабильна (масса заякорена); оконные\n",
    "средние $C_d$ в окнах 6–10: 1.797, 1.813, 1.795, 1.797, 1.791 — плато; rms $C_l$ насыщен.\n",
    "\n",
    "**Выводы по считыванию** (всё подтверждается):\n",
    "1. Два независимых метода согласуются в пределах 4% (систематика +4% объяснена опущенным\n",
    "   конвективным членом и закрыта — раздел 21).\n",
    "2. $\\langle C_m\\rangle\\approx-3\\cdot10^{-5}$ при масштабе $C_l\\sim0.7$ — симметрия в норме.\n",
    "3. $C_d=1.799$ воспроизводит предыдущую итерацию бит-в-бит: осреднение силы по подшагам L1\n",
    "   не сместило среднее (регрессии нет, ряды стали глаже).\n",
    "4. Следовательно $C_d\\approx1.8$ — **истинная сила в смоделированной постановке**; вопрос\n",
    "   к расхождению с литературой — вопрос к *постановке*, а не к считыванию.\n",
    "\n",
    "## 23. Эксперимент №0: серия по боковой блокировке — гипотеза ОПРОВЕРГНУТА\n",
    "\n",
    "**Гипотеза**: завышение $C_d$/rms $C_l$ объясняется боковым стеснением канала; тогда при\n",
    "$D=16$ фикс. и $N_y\\uparrow$ величины должны монотонно стремиться к безграничной литературе,\n",
    "а экстраполяция $\\beta\\to0$ — попасть в $C_d\\approx1.33$, rms $C_l\\approx0.3$.\n",
    "\n",
    "**Постановка**: $N_y\\in\\{90,128,180\\}$ → $\\beta\\in\\{0.178,0.125,0.089\\}$; центр и патч\n",
    "выводятся из $N_y$; всё остальное фиксировано. 3×100k шагов (по ~282 с).\n",
    "\n",
    "**Результат** (`python run_blockage.py`, 2026-06-09):\n",
    "\n",
    "| $N_y$ | β | St | St(1−β) | $\\langle C_d\\rangle$ | rms $C_l$ | $\\langle C_m\\rangle$ | $C_d$ σ·n | dCd |\n",
    "|---|---|---|---|---|---|---|---|---|\n",
    "| 90 | 0.178 | 0.221 | 0.181 | 1.799 | 0.7051 | −0.00003 | 1.872 | 4.1% |\n",
    "| 128 | 0.125 | 0.196 | 0.171 | 1.771 | 0.6734 | −0.00007 | 1.838 | 3.8% |\n",
    "| 180 | 0.089 | 0.193 | 0.176 | **1.870** | **0.8354** | +0.00043 | 1.937 | 3.6% |\n",
    "\n",
    "Экстраполяция $\\beta\\to0$ (интерсепт линейного / квадратичного фита, спред = неопределённость):\n",
    "\n",
    "| | лин. | квадр. | спред | лит. |\n",
    "|---|---|---|---|---|\n",
    "| $C_d(0)$ | 1.905 | 1.855 | 0.050 | 1.33 |\n",
    "| rms $C_l(0)$ | 0.910 | 0.818 | 0.091 | ~0.3 |\n",
    "| St(0) | 0.161 | **0.181** | 0.020 | **0.183** |\n",
    "\n",
    "![Серия по блокировке: Cd(β), rms Cl(β), St(β) с экстраполяцией](solver_2x_sdf/out/blockage_extrapolation.png)\n",
    "\n",
    "**Анализ:**\n",
    "1. $C_d(\\beta)$ **немонотонен** (1.80 → 1.77 → 1.87), rms $C_l$ при самом широком канале даже\n",
    "   вырос — экстраполяция «$C_d(0)\\approx1.9$» физически абсурдна и означает одно: **β — не\n",
    "   управляющий параметр**. Поправка $(1-\\beta)^2$ из предыдущей итерации была совпадением на\n",
    "   одной точке ($N_y=90$).\n",
    "2. Это согласуется с классической теорией продувок: поправки Аллена–Винченти\n",
    "   (solid + wake blockage) при $\\beta=0.178$ дают завышение лишь ~15–18%, при $\\beta=0.089$ —\n",
    "   ~10%. Боковое стеснение объясняло максимум половину превышения даже на исходной сетке.\n",
    "3. **St — единственная величина с правильным поведением**: raw 0.221→0.196→0.193 монотонно\n",
    "   стремится к 0.183 сверху; квадратичная экстраполяция даёт **0.181 ≈ лит. 0.183**.\n",
    "   Кинематика вихреобразования (частота, эффективное Re) верна; завышены именно силовые\n",
    "   (давленческие) амплитуды.\n",
    "4. **Улика — поведение $N_y=180$**: низкочастотная модуляция с размахом оконных $C_d$\n",
    "   1.68–1.98 и rms $C_l$ 0.58–0.96 (период ~20–30k шагов; на 50k установившегося режима —\n",
    "   всего 2–3 цикла, статистика среднего не набрана; $\\langle C_m\\rangle$ на порядок больше,\n",
    "   чем на узких каналах). Чем шире канал, тем свободнее меандрирует дорожка — и тем сильнее\n",
    "   она бьётся о **жёсткое условие $u_y=0$ на выходе**, порождая отражения и обратную связь.\n",
    "5. **Диагноз**: доминируют **продольные ограничения**, которые серия держала константой:\n",
    "   вход 2.5D до тела (жёсткий равномерный профиль), выход 8.3D (Zou-He давление + $u_y=0$ —\n",
    "   отражает вихри), стык L1→L0 на 6.25D. Их изолирует эксперимент №1.\n",
    "\n",
    "## 24. Эксперимент №1: факторное исследование продольных границ\n",
    "\n",
    "**Метод**: каждый подозреваемый фактор варьируется **отдельно** при фиксированной боковой\n",
    "блокировке ($N_y=90$, $\\beta=0.178$) и 100k шагах — изоляция вкладов входа, выхода и\n",
    "выходного ГУ. Реализация: `config.py` выводит x-границы патча из `cx` (патч едет за телом),\n",
    "`boundary.channel_bc` получил режим `outlet_uy`:\n",
    "- `\"zero\"` — классический Zou-He: $u_y=0$ жёстко (как во всех прогонах выше);\n",
    "- `\"extrapolate\"` — $u_y$ нуль-градиент (берётся из столбца $x{=}-2$), якорь $\\rho=1$\n",
    "  сохраняется; алгебраически: $f_6 \\mathrel{+}= \\tfrac12\\rho u_y$, $f_7 \\mathrel{-}= \\tfrac12\\rho u_y$\n",
    "  при тех же $f_3$ и $u_x$ (балансы массы и импульса выполняются точно).\n",
    "\n",
    "| Фактор | Nx | cx | вход | выход | $u_y$-выход | что изолирует |\n",
    "|---|---|---|---|---|---|---|\n",
    "| B0_baseline | 174 | 40 | 2.5D | 8.3D | zero | контроль (+ кросс-чек с конвективным членом) |\n",
    "| F1_outlet_far | 300 | 40 | 2.5D | 16.2D | zero | близость выхода |\n",
    "| F2_inlet_far | 230 | 96 | 6.0D | 8.3D | zero | близость входа |\n",
    "| F3_inout_far | 390 | 96 | 6.0D | 18.3D | zero | оба продольных расстояния |\n",
    "| F4_outlet_uy | 174 | 40 | 2.5D | 8.3D | extrapolate | кинематическое отражение вихрей ГУ |\n",
    "| F5_combo | 390 | 96 | 6.0D | 18.3D | extrapolate | «чистая» постановка |\n",
    "\n",
    "Метрики каждого прогона: St, $\\langle C_d\\rangle$, rms $C_l$, $\\langle C_m\\rangle$, кросс-чек\n",
    "dCd (уже с конвективным членом), **std оконных средних** (мера НЧ-модуляции), $\\langle\\rho\\rangle$.\n",
    "Каждый прогон фиксируется строкой в накопительном реестре `solver_2x_sdf/out/factors.csv`\n",
    "(полная конфигурация + результаты + время) — документация всех изменений постановки.\n",
    "\n",
    "**Результаты** (`python run_factors.py`, 2026-06-10, 6×100k шагов по ~280 с):\n",
    "\n",
    "| Фактор | вход | выход | $u_y$-выход | St(1−β) | $\\langle C_d\\rangle$ | rms $C_l$ | dCd кросс | $\\langle\\rho\\rangle_{фин}$ | модуляция | статус |\n",
    "|---|---|---|---|---|---|---|---|---|---|---|\n",
    "| B0_baseline | 2.5D | 8.3D | zero | 0.181 | 1.799 | 0.705 | **1.0%** | 1.000 | 0.024 | ok |\n",
    "| F1_outlet_far | 2.5D | 16.2D | zero | 0.163 | 2.838 | 6.794 | 1.6% | **1.663** | 0.225 | режим смещён |\n",
    "| F2_inlet_far | 6.0D | 8.3D | zero | 0.138 | 2.229 | 6.120 | 9.7% | **1.655** | 0.306 | режим смещён (срыв на ~75k) |\n",
    "| F3_inout_far | 6.0D | 18.3D | zero | 0.160 | 2.577 | 5.348 | 2.0% | **1.660** | 0.269 | режим смещён (срыв на ~20k) |\n",
    "| **F4_outlet_uy** | 2.5D | 8.3D | extrapolate | **0.183** | **1.754** | **0.557** | **0.6%** | **1.0000** | **0.023** | **ok ✓** |\n",
    "| F5_combo | 6.0D | 18.3D | extrapolate | — | — | — | — | NaN | — | **взрыв на 4.5k** |\n",
    "\n",
    "![Факторное исследование: Cd, rms Cl и НЧ-модуляция по факторам](solver_2x_sdf/out/factors_summary.png)\n",
    "\n",
    "**Анализ:**\n",
    "1. **Конвективный член подтверждён**: кросс-чек B0 сжался с 4.1% до **1.0%** (F4 — 0.6%) при\n",
    "   неизменной физике. Считывание силы теперь согласовано двумя методами на уровне ~1%.\n",
    "2. **F4 (мягкий $u_y$-выход) — выигрыш по всем метрикам**: rms $C_l$ 0.705→0.557 (−21%),\n",
    "   $\\langle C_d\\rangle$ 1.799→1.754, $\\langle\\rho\\rangle=1.0000$ (точнее базы), нулевая\n",
    "   НЧ-модуляция (окна 6–10: $C_d$ 1.753–1.755, идеальное плато) и **St·(1−β)=0.1828 ≈ лит.\n",
    "   0.183 точно**. Жёсткое $u_y=0$ действительно кинематически отражало вихри.\n",
    "3. **Главное открытие — длинные домены без демпфера невалидны.** Все три (F1/F2/F3, $u_y=0$)\n",
    "   уходят на смещённую ветвь $\\langle\\rho\\rangle\\approx1.66$ (F1 и F3 — за ~20k шагов, F2 —\n",
    "   на ~75k), где $C_d$, $C_l$ бессмысленны для сравнения (у F2 разваливается и кросс-чек:\n",
    "   9.7% — поле у контура нестационарно «звенит»). Механизм: канал «вход-скорость +\n",
    "   выход-давление» — **недодемпфированный акустический резонатор**; затухание продольной моды\n",
    "   $\\propto\\nu(\\pi/N_x)^2$ падает в ~5 раз при $N_x$ 174→390. Стартовый транзиент раскачивает\n",
    "   моду, вход с фиксированной скоростью качает массу $\\propto\\rho_{in}$ (положительная обратная\n",
    "   связь) — система садится на новую ветвь равновесия потоков массы. **F5** — чистая картина\n",
    "   роста: $\\langle\\rho\\rangle$ по окнам осциллирует 1.04→0.93→1.05 с периодом\n",
    "   $\\approx 2N_x/c_s\\approx1350$ шагов (точно акустический период туда-обратно) → NaN на 4.5k.\n",
    "4. При этом rms $u_y$ зонда нормальный (0.033–0.046) во **всех** прогонах — гидродинамика\n",
    "   дорожки жива; ломаются акустика и массовый баланс. Это ретроспективно объясняет и «загадку»\n",
    "   эксперимента №0: $N_y=180$ был зачатком той же моды (НЧ-модуляция), но поперечное расширение\n",
    "   менее опасно продольного.\n",
    "\n",
    "**Урок**: продольное удлинение требует **поглощающего слоя** — стандартное лекарство из\n",
    "каталога ГУ (раздел 15), упоминавшееся ещё в прежней 3D-реализации.\n",
    "\n",
    "## 25. Эксперимент №2: губка (absorbing layer) перед выходом\n",
    "\n",
    "Реализация: в последних `sponge_len` столбцах L0 вязкость плавно (smoothstep) растёт до\n",
    "$\\nu\\cdot$`sponge_nu_mult` (по умолчанию ×30 — локальное Re падает до ~5, вихри и акустика\n",
    "диссипируют до прихода на выход). Технически $\\beta$ на L0 становится **полем** $\\beta(x)$\n",
    "(KBC это допускает: $\\gamma$ и так поузловой; CUDA-граф не страдает — поле предвычислено);\n",
    "патч L1 губку не видит (гарантировано assert'ом). Ручки: `AMR_SPONGE_LEN`, `AMR_SPONGE_NU`.\n",
    "\n",
    "| Фактор | геометрия | $u_y$-выход | губка | что проверяет |\n",
    "|---|---|---|---|---|\n",
    "| F6_sponge | 174×90 (база) | extrapolate | 32 | валидация губки против F4 (не портит ли короткий домен) |\n",
    "| F7_long_sponge | 390×90, вход 6D / выход 18.3D | extrapolate | 32 | **целевая «чистая» постановка** |\n",
    "| F8_long_sp_uy0 | 390×90 | zero | 32 | изоляция: достаточно ли губки без мягкого $u_y$ |\n",
    "\n",
    "> ⏳ *Числа появятся после `python run_factors.py` с факторами F6–F8 (3×100k ≈ 15 мин).*\n",
    "\n",
    "**Критерии**: у всех $\\langle\\rho\\rangle\\approx1.000$ и нет НЧ-модуляции; F7 — ожидание\n",
    "$\\langle C_d\\rangle\\approx1.5$–1.6 (остаётся боковая блокировка β=0.178: по Аллену–Винченти\n",
    "+15–18% к безграничному 1.33) и rms $C_l<0.45$. Если критерии выполнены — **эксперимент №3**:\n",
    "повторная серия по β ($N_y=90/128/180$) на конфигурации F7 → финальная экстраполяция к\n",
    "безграничному пределу. Остаточное расхождение после неё — фактор разрешения (D-рефайнмент)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0afa973f",
   "metadata": {},
   "source": [
    "---\n",
    "# Часть VII. Перенос на D3Q27 и связь с целевой реализацией\n",
    "\n",
    "Вся математика KBC переносится в 3D **без изменений по сути**: решётка D3Q27 (веса\n",
    "$8/27, 2/27, 1/54, 1/216$ — изотропия 4-го порядка, предпочтительна для KBC/турбулентности),\n",
    "моментный базис 27×27, сдвиговая часть $s$ — девиаторный тензор (5 независимых компонент).\n",
    "Формула $\\gamma$ и оператор $f\\leftarrow f-\\beta(2\\Delta s+\\gamma\\Delta h)$ идентичны.\n",
    "\n",
    "![Решётка D3Q27](figures/11_d3q27.png)\n",
    "\n",
    "Чем настоящий KBC отличается от прежнего подхода в 3D-движке (poly eq + TRT + α-лимитер +\n",
    "Smagorinsky):\n",
    "\n",
    "| Аспект | Настоящий KBC | Прежний подход |\n",
    "|---|---|---|\n",
    "| Равновесие | энтропийное product-form | полиномиальное $O(u^2)$ |\n",
    "| Столкновение | энтропийный MRT, $k/s/h$ | TRT (симм./антисимм.) |\n",
    "| Стабилизация | $\\gamma$ из энтропии | α-лимитер положительности |\n",
    "| H-теорема | встроена в вывод $\\gamma$ | контролируется отдельно |\n",
    "| Турбулентность | implicit LES (параметр-free) | явный Smagorinsky |\n",
    "| Решётка | D3Q27 | D3Q27 (совпадает) |\n",
    "\n",
    "**Дорожная карта переноса** (всё отработано на 2D-макете): (1) product-form $f^{eq}$;\n",
    "(2) проекторы $P_s$ для D3Q27 (предвычислить); (3) per-node $\\gamma$; (4) столкновение;\n",
    "(5) SDF+Bouzidi и GMEM-сила с $\\mathbf u_w$ на линке (подвижные тела) + кросс-чек ∮σ·n dA по\n",
    "SDF-изоповерхности; (6) AMR-связка $\\tau_f=2\\tau_c-\\tfrac12$ с временно́й интерполяцией."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d0691787",
   "metadata": {},
   "source": [
    "---\n",
    "# Часть VIII. Сверка формул с первоисточниками\n",
    "\n",
    "Все формулы, ранее помеченные «⚠ сверить с оригиналом», проверены по первоисточникам\n",
    "(плюс численная самопроверка каждой в `demos_gpu/checks_gpu.py`):\n",
    "\n",
    "| Формула | Источник | Статус |\n",
    "|---|---|---|\n",
    "| 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) | ✓ дословно |\n",
    "| $f' = f-\\beta(2\\Delta s+\\gamma\\Delta h)$, $\\beta=1/(2\\tau)$, mirror-форма | Karlin/Bösch/Chikatamarla, PRE **90**, 031302(R) (2014) | ✓ дословно |\n",
    "| $\\gamma^\\*=\\frac1\\beta-(2-\\frac1\\beta)\\frac{\\langle\\Delta s|\\Delta h\\rangle}{\\langle\\Delta h|\\Delta h\\rangle}$, $\\langle X|Y\\rangle=\\sum X_iY_i/f_i^{eq}$ | там же; PRE **92**, ур. (24)–(25) | ✓ дословно |\n",
    "| Состав $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) | ✓ |\n",
    "| 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) | ✓ дословно |\n",
    "| Временна́я интерполяция ghost-границы | Lagrava/Malaspinas/Latt/Chopard, JCP **231**, 4808 (2012) | ✓ |\n",
    "| Bouzidi: ветви $q<\\frac12$ / $q\\ge\\frac12$ | Bouzidi/Firdaouss/Lallemand, Phys. Fluids **13**, 3452 (2001) | ✓ |\n",
    "| Zou-He: скоростной вход / давление-выход | Zou & He, Phys. Fluids **9**, 1591 (1997) | ✓ |\n",
    "| 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) | ✓ |\n",
    "| $\\sigma^v=-(1-\\frac{1}{2\\tau})\\Pi^{neq}$ (Чепмен–Энског) | стандарт (Krüger et al., *The Lattice Boltzmann Method*, 2017) | ✓ |\n",
    "\n",
    "**Литература:**\n",
    "- I. V. Karlin, F. Bösch, S. S. Chikatamarla, *Gibbs' principle for the lattice-kinetic theory\n",
    "  of fluid dynamics*, Phys. Rev. E **90**, 031302(R) (2014).\n",
    "- F. Bösch, S. S. Chikatamarla, I. V. Karlin, *Entropic multirelaxation lattice Boltzmann\n",
    "  models for turbulent flows*, Phys. Rev. E **92**, 043309 (2015) (arXiv:1507.02518).\n",
    "- S. Ansumali, I. V. Karlin, H. C. Öttinger, *Minimal entropic kinetic models for\n",
    "  hydrodynamics*, Europhys. Lett. **63**, 798 (2003).\n",
    "- B. Dorschner, F. Bösch, S. S. Chikatamarla, K. Boulouchos, I. V. Karlin, *Entropic\n",
    "  multi-relaxation time lattice Boltzmann model for complex flows*, J. Fluid Mech. **801**,\n",
    "  623–651 (2016).\n",
    "- B. Dorschner, N. Frapolli, S. S. Chikatamarla, I. V. Karlin, *Grid refinement for entropic\n",
    "  lattice Boltzmann models*, Phys. Rev. E **94**, 053311 (2016) (arXiv:1608.06915).\n",
    "- D. Lagrava, O. Malaspinas, J. Latt, B. Chopard, *Advances in multi-domain lattice Boltzmann\n",
    "  grid refinement*, J. Comput. Phys. **231**, 4808 (2012).\n",
    "- M. Bouzidi, M. Firdaouss, P. Lallemand, *Momentum transfer of a Boltzmann-lattice fluid with\n",
    "  boundaries*, Phys. Fluids **13**, 3452 (2001).\n",
    "- Q. Zou, X. He, *On pressure and velocity boundary conditions for the lattice Boltzmann BGK\n",
    "  model*, Phys. Fluids **9**, 1591 (1997).\n",
    "- B. Wen, C. Zhang, Y. Tu, C. Wang, H. Fang, *Galilean invariant fluid–solid interfacial\n",
    "  dynamics in lattice Boltzmann simulations*, J. Comput. Phys. **266**, 161 (2014).\n",
    "- Z. Guo, C. Zheng, B. Shi, *Discrete lattice effects on the forcing term in the lattice\n",
    "  Boltzmann method*, Phys. Rev. E **65**, 046308 (2002).\n",
    "- U. Ghia, K. N. Ghia, C. T. Shin, *High-Re solutions for incompressible flow using the\n",
    "  Navier-Stokes equations and a multigrid method*, J. Comput. Phys. **48**, 387 (1982)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a8cc8c0d",
   "metadata": {},
   "source": [
    "---\n",
    "# Часть IX. Итоги\n",
    "\n",
    "- **KBC** = энтропийный MRT-LBM: сохраняющиеся моды не трогаем, сдвиг релаксируем фиксированным\n",
    "  $2\\beta$ (вязкость), высшие моды — адаптивным $\\gamma$ из условия неувеличения энтропии.\n",
    "  При $\\gamma=2$ — в точности BGK. Все формулы сверены с первоисточниками и самопроверены\n",
    "  численно.\n",
    "- **Бенчмарки** (на компонентах финальной модели): вязкость по TGV (<6%, факт ≲1%), H-теорема,\n",
    "  устойчивость на сдвиговом слое (BGK разрушается — KBC держит), каверна ≈ Ghia, дорожки\n",
    "  Кармана на четырёх телах.\n",
    "- **Измельчение сетки**: связка $\\tau_f=2\\tau_c-\\tfrac12$ + масштабирование $f^{neq}$ +\n",
    "  временна́я интерполяция ghost; вложенный 2×+4× даёт точность мелкой сетки за в разы меньшую\n",
    "  работу.\n",
    "- **Финальная модель** `solver_2x_sdf`: KBC-N1 + SDF/Bouzidi + Zou-He + free-slip + CUDA-графы.\n",
    "  Главные уроки эволюции: второй член силы — из post-streaming; дрейф массы лечится\n",
    "  давлением-выходом; зонд и сила — на тонком уровне; графы требуют чистоты горячего пути.\n",
    "- **Считывание сил ДОКАЗАНО корректным** (разделы 22, 24): два независимых метода — GMEM с\n",
    "  торком и баланс импульса ∮[σ·n−ρu(u·n)]ds — после добавления конвективного члена согласуются\n",
    "  на уровне **0.6–1.0%** по $\\langle C_d\\rangle$; $\\langle C_m\\rangle\\approx0$; масса стабильна.\n",
    "- **Исследование факторов постановки** (разделы 23–25):\n",
    "  - эксперимент №0 (серия по β) **опроверг** гипотезу бокового стеснения: $C_d(\\beta)$\n",
    "    немонотонен; при этом **St(β→0)=0.181 ≈ лит. 0.183** — кинематика верна;\n",
    "  - эксперимент №1 (6 факторов): **мягкий $u_y$-выход (F4) улучшил всё** (rms $C_l$ −21%,\n",
    "    St·(1−β)=0.183 точно, идеальное плато) — жёсткое $u_y=0$ отражало вихри. **Длинные домены\n",
    "    без демпфера невалидны**: пара «вход-скорость / выход-давление» — недодемпфированный\n",
    "    акустический резонатор (затухание ∝ν(π/Nx)²); масса накачивается до смещённой ветви\n",
    "    ⟨ρ⟩≈1.66 либо взрыв (F5, период осцилляций ⟨ρ⟩ ≈ 2Nx/c_s — прямая улика);\n",
    "  - эксперимент №2 (в работе): **губка перед выходом** (β(x)-поле на L0) — F6/F7/F8; затем\n",
    "    эксперимент №3 — серия по β на чистой постановке F7 → финальная экстраполяция.\n",
    "- **Реестр всех прогонов** — `solver_2x_sdf/out/factors.csv`: полная конфигурация + метрики\n",
    "  каждого изменения постановки (готовый материал для отдельного исследования факторов).\n",
    "- **Дальше**: эксперименты №2–3 → финальные числа; при остаточном превышении — фактор\n",
    "  разрешения (D-рефайнмент). Затем перенос отработанной схемы на D3Q27/GPU в основной движок\n",
    "  (часть VII)."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "version": "3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}