CFDManager/docs/theory/2d_solver/README.md

# kbc2d — двумерный решатель LBM D2Q9 с энтропийным столкновением KBC

Переписанный на Rust решатель обтекания тела в канале. Физика — та же, что в
`docs/theory/solver_2x_sdf` (python/CuPy), но собранная заново: с тестами против формул
первоисточников, двумя взаимозаменяемыми бэкендами и анимацией, привязанной к физическому
времени потока.

Оператор столкновения — энтропийный KBC по работам Bösch, Chikatamarla, Karlin: *Entropic
Multi-Relaxation Models for Simulation of Fluid Turbulence* (arXiv:1507.02509, двумерная
реализация) и *Entropic multi-relaxation lattice Boltzmann scheme for turbulent flows* (2024,
трёхмерная). Состав сдвиговой части выбирается ключом `--kbc-model`: `n1` — только девиатор
{N, Π_xy} (в 2D-статье KBC D), `n2` — девиатор со следом {N, Π_xy, T} (KBC C). Статьи лежат в
`docs/origins/`; ссылки на формулы в коде и ниже даны по нумерации 2D-статьи.

## Устройство

Пять файлов по ролям — каждый отвечает ровно за одно:

| файл | роль |
|---|---|
| `src/math.rs` | **математика решателя.** Решётка D2Q9, энтропийное равновесие в product-form, проектор на сдвиг, стабилизатор γ, столкновение, Zou–He, геометрия тел через SDF, перевод единиц, спектральная диагностика. Всё узловое и чистое; здесь же тесты против формул статьи. |
| `src/cpu.rs` | **бэкенд под процессор.** Раскладка AoS, обход сетки, rayon, сборка Bouzidi-линков, связка уровней AMR. Физику берёт из `math`. |
| `src/gpu.rs` | **бэкенд под видеокарту.** wgpu + WGSL (Vulkan/DX12/Metal), раскладка SoA. Физика построчно повторяет `math.rs` на f32; топологию (маски, линки, рамка патча) не дублирует, а берёт из `cpu`. |
| `src/main.rs` | **запуск и оркестрирование.** Разбор параметров, сборка постановки, цикл по шагам, живой вывод и итоговый отчёт, выгрузка рядов в CSV. Здесь же контракт `Spec` / `StepRec` / `FieldKind`, общий для обоих бэкендов. |
| `src/gif.rs` | **создание гифок.** Тайминг относительно физического времени, палитры, нормировка, служебная надпись, кодирование. |

## Сборка и запуск

Нужен Rust 1.75+.

```sh
cargo build --release                       # с GPU-бэкендом
cargo build --release --no-default-features  # только CPU (без wgpu)
cargo test --release                         # 22 быстрых теста
cargo test --release -- --include-ignored    # плюс эталоны статьи (~18 с)
```

Пример: цилиндр Re=150, гифка завихренности в реальном времени.

```sh
./target/release/kbc2d --shape cylinder --size 24 --re 150 \
    --nx 480 --ny 240 --steps 40000 --sponge-len 32 \
    --gif wake.gif --gif-field vorticity --verbose full
```

`--help` показывает все ключи, разбитые по группам: Физика, Сетка, Тело, Время, Схема,
Анимация, Вывод.

## Синхронизация анимации с физическим временем

Требование: гифка идёт с той же скоростью, что и настоящий поток, независимо от того, с какой
скоростью считает машина. Реальная производительность в тайминг не входит вообще.

В LBM скорость самой решётки жёстко равна c = δx/δt = 1, поэтому шаг по времени однозначно
определяется тем, какую **решёточную** скорость `u_lat` мы назначаем физическому потоку:

```
δt = u_lat · δx / u_phys          [с/шаг]
шагов в секунду = u_phys / (u_lat · δx)
задержка кадра  = (шагов на кадр) · δt / playback
```

**Про пример из постановки.** «30 м/с, ячейка 0.1 м ⇒ 300 шагов/с» — это арифметика
δt = δx/u_phys, то есть `u_lat = 1`: поток проходит ровно ячейку за шаг. Формула
воспроизводится буквально ключом `--u-lat 1.0` (тест `units_time_scaling` это проверяет), но
физически такой режим негоден: Ma = u_lat/c_s = √3 ≈ 1.73, сверхзвук, разложение
Чепмена–Энскога не работает. Поэтому по умолчанию `u_lat = 0.05` (Ma ≈ 0.087), и те же 30 м/с
при ячейке 0.1 м дают 6000 шагов/с. Синхронность гифки выдерживается в обоих случаях — меняется
только, сколько шагов приходится на кадр.

**Две тонкости формата GIF**, обе разобраны:

1. *Задержка хранится в сотых долях секунды.* Точная физическая задержка почти никогда не целая:
   30 кадр/с — это 3⅓ сотых. Покадровое округление до 3 дало бы анимацию на 11% быстрее
   реальности, и уход копился бы линейно (к тысячному кадру — 3.3 секунды). Поэтому задержки
   выдаёт накопитель `DelayDither`: суммарное время кадров отслеживает точное физическое с
   точностью до одной сотой (3, 3, 4, 3, 3, 4, …), ошибка ограничена ±5 мс и не растёт.
2. *Меньше одной сотой не бывает.* При физически корректном `u_lat` кадр каждые 10 шагов — это
   600 кадр/с, чего формат не умеет. Поэтому режим по умолчанию `--gif-every auto` подбирает
   шаг сам под `--gif-fps` (30) так, чтобы получилось ровно реальное время. Если шаг задан
   жёстко и задержка не представима, программа печатает фактический коэффициент расхождения и
   конкретный совет, как починить.

`--gif-speed 0.1` даёт замедление в 10 раз (тоже точно, через тот же накопитель).

## Шаг — не фиксированная порция времени

Физическая длительность шага привязана к размеру клетки: δt = u_lat·δx/u_phys. Уменьшили клетку
вдвое — вдвое уменьшился и δt, и то же число шагов покроет вдвое меньше физического времени.
Отсюда `--time`: задаёте длительность в секундах, число шагов считается само.

Измельчение стоит дважды: клеток становится (1/δx)², а шагов на ту же секунду — 1/δx, итого
работы ~(1/δx)³ в двумерии. И помните, что `--size` задаётся В КЛЕТКАХ: уменьшив клетку и не
тронув `--size`, вы уменьшите тело физически.

Судить о длительности удобнее всего по **конвективным временам D/U** — их печатает шапка. Это
единственная мера, не зависящая ни от сетки, ни от выбора u_lat.

## Что параметризовано

- **поток**: скорость (м/с), направление, число Рейнольдса, решёточная скорость (число Маха);
- **сетка**: размеры домена и размер ячейки в метрах (плотность сетки), коэффициент измельчения
  вложенного патча и его границы;
- **тело**: семь форм на выбор — `cylinder`, `square`, `diamond`, `ellipse`, `naca`, `triangle`,
  `plate` — плюс характерный размер, относительная толщина, угол атаки и положение;
- **время**: длительность прогона — либо числом шагов (`--steps`), либо прямо в СЕКУНДАХ
  физического времени (`--time`, число шагов считается как time/δt); начальное поле (однородный
  поток либо покой с разгоном), длина разгона, амплитуда и длительность стартового возмущения;
- **схема**: оператор столкновения (`kbc`/`bgk`), состав сдвиговой части (`n1`/`n2`), модель
  стенки на теле (`bouzidi`/`grad`/`staircase`), режим выхода, поглощающая губка, бэкенд,
  число потоков;
- **анимация**: файл, поле (`speed`/`vorticity`/`density`/`gamma`), палитра, масштаб, шаг кадра,
  частота, скорость воспроизведения, диапазон нормировки;
- **вывод**: период живых строк, три уровня подробности, число окон в отчёте о сходимости, CSV.

## Что печатает отчёт

*Шапка* — вся постановка с производными величинами: δt, шагов на секунду, Ma, τ, физическая
вязкость, блокировка канала, геометрия патча, полный план тайминга анимации.

*Живой вывод* — шаг, физическое время, ⟨ρ⟩, max|u|, Cd, Cl, скорость счёта и ETA; на уровне
`full` дополнительно ⟨γ⟩ с размахом, доля вырожденных узлов, доля узлов с ξ < 0, MLUPS и
отношение скорости счёта к реальному времени.

*Итог* — установившийся режим (St, ⟨Cd⟩, rms Cl, ⟨Cm⟩) сырой и с поправкой на блокировку, рядом
литературные значения для цилиндра; таблица сходимости по окнам с вердиктом о дрейфе массы и
насыщении; разбор стабилизатора γ; производительность.

## Состояние проверки

Опора — статьи авторов метода из `docs/origins/`, а не сторонние реализации. Ссылки на формулы
даны по нумерации 2D-статьи (arXiv:1507.02509); там, где полезнее формулировка из работы 2024
года по трёхмерной реализации, это отмечено отдельно.

### Оператор столкновения сверен с листингом статьи

Работа 2024 года приводит оператор явным пошаговым листингом (разд. IV). Реализация повторяет
его дословно: ρ, u → f^eq → s и s^eq → Δs = s − s^eq → **Δh = h − h^eq = f − f^eq − Δs** →
γ по замкнутой оценке → **f′ = f − β(2Δs + γΔh)**.

Проверено тестами (`cargo test`, 22 быстрых + 4 длинных):

- **проектор Δs** совпадает с матричным `M⁻¹·D·M` в базисе натуральных моментов (6)–(7) до 1e-13
  — для обоих составов сдвиговой части; идемпотентен, не несёт ни массы, ни импульса;
- **γ из замкнутой оценки** (ур. 17 / ур. 25 работы 2024) — корень условия максимума энтропии
  (ур. 15 / 23): невязка при γ\* более чем в 20 раз меньше, чем при γ\*±1;
- **при γ = 2 схема совпадает с LBGK** поточечно — как и заявлено под ур. (14);
- **сдвиговые моменты релаксируют ровно с 2β при любой γ** (β = 0.3, 0.6, 0.95, обе модели) —
  именно это гарантирует, что стабилизатор не трогает вязкость;
- **вязкость по ур. (5)** воспроизводится затуханием сдвиговой волны точнее 1% (τ = 0.6 и 1.0);
- **объёмная вязкость по ур. (57)**: ξ = ν при следе в сдвиговой части и ξ = c_s²(1/(γβ) − ½)
  без него;
- равновесие в product-form сохраняет ρ и ρu до 1e-13; Zou–He ставит ровно заданные скорость
  на входе и плотность на выходе; SDF всех семи форм даёт верный знак и |∇φ| = 1 ± 0.05.

### Вихрь Тейлора–Грина: второй порядок сходимости (разд. VI)

Единственное из трёх эталонных течений статьи с ТОЧНЫМ аналитическим решением, поэтому
проверяется не «похоже на чужой прогон», а прямое совпадение с формулой:

    u = ∇×[(u₀/k₂)cos(k₁x)cos(k₂y)·exp(−ν(k₁²+k₂²)t)],  k₁ = 1, k₂ = 4,
    область 0 < x,y < 2π на N×N,  Re = u₀N/ν,  полураспад t_c = ln2/[ν(k₁²+k₂²)].

Старт — приближением Града (ур. 58), как в статье. Метрика — как на рис. 1:
Σ|u_x − u_x^точн| / Σ|u_x^точн| в момент t_c.

| N | u₀ | Re | полная | амплитуда | форма |
|---|---|---|---|---|---|
| 64 | 0.03 | 100 | 6.68e-3 | 6.49e-3 | 8.65e-4 |
| 128 | 0.015 | 100 | 1.61e-3 | 1.57e-3 | 2.07e-4 |
| 256 | 0.0075 | 100 | 4.00e-4 | 3.89e-4 | 4.36e-5 |

**Порядок 2.05 и 2.01** — второй порядок статьи воспроизведён, причём отдельно по амплитуде
(скорость затухания) и по форме.

Два места, где пришлось разобраться, и оба поучительны:

1. **Давление в начальных условиях.** Течение несёт собственное поле давления порядка ρu₀²,
   находимое из ∇²p = 2ρ(ψ_xx·ψ_yy − ψ_xy²):
   `p = −(ρu₀²/4)[cos(2k₁x) + (k₁²/k₂²)cos(2k₂y)]`. Старт с ρ ≡ 1 сбрасывает эту разницу в
   акустику, которая в периодическом ящике почти не затухает и садится полкой на ошибку. Работа
   2024 года делает то же самое явно: там начальные ρ и старшие моменты получают, решая
   ∂ρ/∂t + ∇·(ρu₀) = D∇²ρ до стационара.
2. **Способ измельчения.** При фиксированном u₀ ошибка упирается в полку O(Ma²), от сетки не
   зависящую (измерено: относительная ошибка формы ∝ u₀¹·⁰⁷, то есть абсолютная ∝ Ma²).
   Второй порядок виден целиком только при диффузионном измельчении — ν фиксирована, u₀ ∝ 1/N,
   тогда Re сохраняется, а Маха падает вместе с сеткой. Это свойство слабо-сжимаемого метода,
   а не реализации: **LBGK на том же тесте даёт ту же полку** (1.39/1.37/1.36e-3 против
   0.86/1.05/1.23e-3 у KBC), что согласуется с утверждением статьи «все модели работают
   практически одинаково».

### Дважды периодический сдвиговый слой (разд. VII)

Второй эталон статьи: N = 128, Re = 30000, u₀ = 0.04, κ = 80, δ = 0.05, одно конвективное
время. Отношение энстрофии к начальной сходится с fp64-значением 0.6035. Тест чувствителен
именно к тому, что важно: на испорченном (абсолютном) пороге вырожденности γ тот же прогон
даёт 0.6599, то есть +9.3%, а чистый LBGK при этих параметрах разваливается.

### Порог вырожденности γ

`GREL = 1e-8` — **относительный** порог, доля от ⟨Δ|Δ⟩, а не абсолютный. Знаменатель ⟨Δh|Δh⟩
квадратичен по неравновесию и физически мал (~1e-7…1e-9 в развитом следе), поэтому абсолютный
порог срабатывает на подавляющем большинстве узлов и молча подменяет γ на 2 — то есть гонит
чистый LBGK вместо KBC. Работа 2024 года прямо об этом: γ «далеко не постоянна», её эволюция
тесно связана с состоянием потока, и «любой MRT с γ = const не достигнет той же устойчивости».
Доля вырожденных узлов печатается в отчёте; на исправном пороге она обязана быть ~0.

Единственное исключение — самый первый шаг: поле в точности равно равновесию, Δ ≡ 0, и порог
честно срабатывает везде. На результат это не влияет: γ умножается на Δh = 0.

### Выбор состава сдвиговой части (табл. I)

Ключ `--kbc-model`:

- `n1` (умолчание) — s = {N, Π_xy}, только девиатор. 2D-статья: KBC D; 3D: KBC-N1.
- `n2` — s = {N, Π_xy, T}, девиатор со следом. 2D-статья: KBC C; 3D: KBC-N2.

По точности они неразличимы, как и заявляет статья: на одной постановке St 0.1828 у обоих,
⟨Cd⟩ 1.4479 против 1.4468, rms Cl 0.383 против 0.388.

Разница — в объёмной вязкости (ур. 57). У `n2` она фиксирована: ξ = ν. У `n1` она равна
c_s²(1/(γβ) − ½) и, поскольку измеренная ⟨γ⟩ ≈ 1.73 < 2, в среднем оказывается примерно
вчетверо БОЛЬШЕ ν. Отрицательной она бывает лишь в долях процента узлов. Практический вывод
против ожидания: `n1` демпфирует продольную акустику сильнее, и в специально испорченной
постановке (старт из покоя, губка выключена) `n1` доживает до конца с пульсацией 75% от U,
а `n2` разваливается. Поэтому умолчание — `n1`.

### Модель стенки на теле: насколько она субсеточная

Ключ `--wall`:

- `bouzidi` (умолчание) — интерполированный отскок: доля пересечения q входит в КАЖДУЮ
  восстанавливаемую популяцию, полинково;
- `grad` — условие Града (Dorschner, Bösch, Chikatamarla, Boulouchos, Karlin, JFM 801 (2016),
  разд. 2.1 и прил. B): задаются не популяции, а целевые моменты — ρ, u и тензор давлений, —
  после чего недостающие популяции собираются приближением Града (2.13);
- `staircase` — простой отскок, q игнорируется. Не для счёта: это база сравнения, показывающая,
  сколько именно даёт субсеточность.

**Субсеточность — измеренная.** Прямой тест: сдвигаем тело внутри клетки и смотрим, насколько
поедет Cd. У по-настоящему субсеточной границы ответ не должен зависеть от того, где тело стоит
относительно узлов (Re = 20, D = 16, стационар, сдвиги 0…½ клетки):

| модель | разброс Cd | Cd |
|---|---|---|
| `staircase` | 1.11% | 2.68–2.71 |
| `grad` | 0.64% | 2.64–2.66 |
| `bouzidi` | **0.19%** | 2.653–2.658 |

Град оказывается ровно между ступенькой и Bouzidi, и это следует из его устройства: положение
стенки входит туда ТОЛЬКО через целевую скорость (B 1) — одну усреднённую по узлу величину.
Целевая плотность (B 3) — обычная сумма отскочивших и известных популяций, без q вовсе; тензор
давлений — конечные разности по решётке, тоже без q. Плюс все недостающие популяции узла
собираются из ОДНОГО набора моментов, так что полинковая направленность теряется. Bouzidi же
подставляет свою q в каждую популяцию отдельно. Ступенчатой поверхность у Града не становится,
но геометрия у него разрешена заметно грубее.

**Сходимость по разрешению тела.** Физическая постановка фиксирована (домен 15D × 10D,
блокировка 0.1, Re = 20), меняется только число клеток на диаметр:

| D | `bouzidi` | `grad` |
|---|---|---|
| 8 | 2.581 | 2.618 |
| 16 | 2.529 | 2.534 |
| 32 | **2.521** | **2.521** |

Обе модели состоятельны и сходятся к одному пределу с наблюдаемым порядком ≈2.7; к D = 32 они
неразличимы. Но на грубой сетке Град заметно хуже: ошибка при D = 8 равна 0.097 против 0.060.
Для сравнения, `staircase` при D = 16 даёт 2.69 — то есть +6.7% к пределу, тогда как обе
субсеточные модели держатся в пределах +0.4%.

**Зачем тогда Град.** Его преимущество в статье — не геометрическая точность, а устойчивость
на турбулентных режимах (авторы пишут, что интерполяционные схемы «ограничены низкими числами
Рейнольдса, поскольку на границе возникают паразитные скачки») и естественная форма для
подвижных стенок: скорость стенки входит в целевые значения, а не отдельной поправкой.
В здешней канальной постановке преимущества по устойчивости воспроизвести не удалось: при росте
Re обе модели теряют счёт на одном и том же значении (Re ≈ 5·10⁴ при теле в 16 клеток), то есть
ограничивает не стенка, а что-то другое — вероятнее всего Zou–He при τ → ½. Поэтому умолчание —
`bouzidi`, а `grad` стоит держать в виду для будущих подвижных тел.

Ограничение реализации: `--wall grad` пока только на процессорном бэкенде; GPU при таком выборе
отказывается запускаться явно, а не считает молча по Bouzidi.

### Паритет бэкендов и согласованность уровней

CPU (f64) и GPU (f32) на одной постановке совпадают до 4–5 значащих цифр шаг в шаг: ⟨ρ⟩
1.04933 против 1.04934, Cd 2.339 против 2.338, ⟨γ⟩ 1.2645 против 1.2646. На Intel Iris Xe GPU
даёт ≈105 MLUPS против ≈18 MLUPS у процессора.

Один и тот же случай с патчем ×2 и вовсе без измельчения (`--refine 1`) даёт St 0.1951 против
0.1970 и ⟨Cd⟩ 1.956 против 1.943 — связка уровней систематики не вносит.

### Обтекание цилиндра против литературы

Постановка 480×240, D = 24, Re = 150, блокировка β = D/Ny = 0.1, умолчания решателя:

| величина | сырое | с поправкой на блокировку | литература (безгранич. цилиндр) |
|---|---|---|---|
| St | 0.1828 | 0.1645 | 0.183 |
| ⟨Cd⟩ | 1.448 | 1.173 | 1.33 |
| rms Cl | 0.383 | 0.310 | ~0.30 |
| ⟨Cm⟩ | 0.00002 | — | 0 (симметрия) |

Поправка: St×(1−β), Cd и rms Cl ×(1−β)². Сырое St и скорректированный rms Cl ложатся на
литературу; Cd после поправки ниже на 12%. ⟨Cm⟩ ≈ 0 — контроль симметрии считывания силы.

**Формы тел ведут себя физично.** Прогон на каждую форму (260×130, размер 20, угол атаки 12°)
даёт ожидаемый порядок сопротивления: профиль 0.44, эллипс 0.45, пластина 0.60, цилиндр 1.33,
ромб 1.65, квадрат 2.38, треугольник 2.72. ⟨Cm⟩ ≈ 0 **только** у круга (−0.0002), которому
угол атаки безразличен, а у несимметричных под углом тел он ненулевой (профиль +0.31,
эллипс +0.15, пластина +0.13).

### Мелкие отличия от питоновского прототипа

Решатель писался заново, не как порт, но пара мест разошлась с `solver_2x_sdf` намеренно:
внутри тела здесь не считается столкновение (эти популяции фиктивны — Bouzidi перекрывает всё,
что могло бы прийти из тела в жидкость; побочно статистика γ собирается строго по жидкости), а
рестрикция дополнительно пропускает узлы, у которых тонкий узел-источник лежит внутри тела.
Процессорный бэкенд работает в f64, GPU-бэкенд — в f32.

## Продольная акустика канала: почему поток может «дышать»

Самая заметная ловушка этой постановки, и её стоит понимать до первого запуска.

**Граничное условие с заданной скоростью на входе акустически есть жёсткий поршень**: оно
отражает продольные волны с коэффициентом +1. Выход по давлению — наоборот, открытый конец.
Вместе они делают из канала четвертьволновый резонатор с пучностью давления на входе и узлом
на выходе:

    период основной моды  =  4·Nx/c_s  шагов
    затухание вязкостью   ~  ν(π/2Nx)²  — на длинном домене практически ноль

Измерено на 480×240: период пульсации ⟨ρ⟩ **3332 шага** против расчётных 4·Nx/c_s = 3325
(0.2%), первый ноль автокорреляции на 827 шагах = ровно Nx/c_s (четверть периода). За 30000
шагов амплитуда упала на 3.8% — то есть мода не гаснет вообще. Ничто её не подкачивает; это
звон от старта, запертый в почти без потерь резонаторе.

Отсюда умолчания:

- **`--init uniform`** — домен сразу заполнен набегающим потоком, вход включён на полную.
  Старт из покоя (`--init rest`) разгоняет весь столб жидкости и закачивает моду; при разгоне
  за 1000 шагов против акустического пробега 831 шаг это для звука удар.
- **губка перед выходом включена и подобрана по домену** (`nx/12`, не меньше 16 столбцов,
  с запасом до патча). Отключается `--sponge-len 0`.
- **`--outlet extrapolate`** — нуль-градиент поперечной скорости; жёсткий ноль отражает вихри
  дорожки обратно к телу.

Эффект умолчаний на той же постановке:

| | размах ⟨ρ⟩ | пульсация u′/U |
|---|---|---|
| старт из покоя, без губки | 0.1307 | **75%** |
| умолчания | 0.0017 | **1.0%** |

Заодно выправились и сами числа: St 0.183 вместо 0.194, rms Cl 0.38 вместо 0.84, ⟨ρ⟩ 1.0000.

Отчёт печатает период моды, время её вязкого затухания и измеренную пульсацию в конце прогона,
с явным предупреждением, если она превысила 5% от U.

## Дальше

- σ·n-кросс-чек силы (интеграл тензора напряжений по контуру) как независимая проверка GMEM;
- согласованные начальные условия по образцу работы 2024 года: там ρ и старшие моменты
  получают, решая ∂ρ/∂t + ∇·(ρu₀) = D∇²ρ до стационара, что убрало бы и остаточный
  стартовый импульс от появления тела в потоке;
- подвижные и вращающиеся тела: GMEM уже записан в галилей-инвариантной форме и принимает
  скорость стенки на линке, но подача этой скорости не подключена;
- несколько тел одновременно (сейчас — одно тело на выбор из семи форм).