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` | **создание гифок.** Тайминг относительно физического времени, палитры, нормировка, служебная надпись, кодирование. |

Рядом, вне этого разделения: `bench/` — валидационная кампания (список сценариев и драйверы под
Linux и Windows), `Dockerfile` с `docker-compose.yml` — образ для развёртывания на сервере.

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

Нужен Rust 1.75+.

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

Сборка без GPU кладёт бинарь по тому же пути, поэтому собирать её ПОСЛЕДНЕЙ нельзя: `kbc2d`
окажется перезаписан вариантом без wgpu, и `--backend gpu` будет отказывать. Порядок —
сначала `--no-default-features`, потом обычная.

На сервер удобнее ставить образом — см. [Развёртывание](#развёртывание-docker).

Пример: цилиндр 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.

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

- **поток**: скорость (м/с), направление, число Рейнольдса, решёточная скорость (число Маха);
- **постановка**: `--case channel` — обтекание тела; `taylor-green`, `shear-layer`,
  `decaying-turbulence` — периодические эталоны из статей, без тела и граничных условий;
- **сетка**: размеры домена и размер ячейки в метрах (плотность сетки), коэффициент измельчения
  вложенного патча и его границы;
- **тело**: восемь форм на выбор — `cylinder`, `square`, `diamond`, `ellipse`, `naca`,
  `triangle`, `plate`, `polygon` — плюс характерный размер, относительная толщина, угол атаки и
  положение. Профиль задаётся четырёхзначным кодом (`--naca 4412`, с кривизной), произвольный
  контур — списком вершин (`--poly "x,y;x,y;…"`), **несколько тел сразу** — списком
  `--bodies "cylinder:d=24,x=120,y=120; naca:d=48,x=260,y=120,a=8"` (сила считается по каждому
  телу отдельно, до четырёх);
- **время**: длительность прогона — либо числом шагов (`--steps`), либо прямо в СЕКУНДАХ
  физического времени (`--time`, число шагов считается как time/δt); начальное поле (однородный
  поток либо покой с разгоном), длина разгона, амплитуда и длительность стартового возмущения;
- **схема**: оператор столкновения (`kbc`/`bgk`), состав сдвиговой части (`n1`/`n2`), модель
  стенки на теле (`hrr`/`grad`/`bouzidi`/`staircase`), режим выхода, поглощающие губки перед
  выходом и после входа, бэкенд, число потоков;
- **анимация**: файл, поле (`speed`/`vorticity`/`density`/`gamma`), палитра, масштаб, шаг кадра,
  частота, скорость воспроизведения, диапазон нормировки, усреднение k×k клеток в пиксель
  (`--gif-downsample`, без него кадр с сетки 4096×2048 неподъёмен);
- **вывод**: период живых строк, три уровня подробности, число окон в отчёте о сходимости, CSV
  рядов с прореживанием (`--series-every`), x–t диаграмма осевой линии (`--xt`), метрики
  эталонных течений (`--case-csv`) и машиночитаемая сводка всего прогона (`--summary`).

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

*Шапка* — вся постановка с производными величинами: δ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`, 29 быстрых + 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
  (у эллипса для этого пришлось считать ближайшую точку итеративно: дешёвое приближение давало
  |∇φ| = 0.57 вдали от поверхности);
- **HRR-сборка сохраняет моменты**, ради которых затевалась: ρ, ρu и Π восстановленной функции
  распределения совпадают с целевыми, а при нулевой скорости она совпадает с Градовой (все
  коэффициенты 3-го порядка рекурсивно обращаются в ноль).

### Вихрь Тейлора–Грина: второй порядок сходимости (разд. 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`:

- `hrr` (умолчание) — восстановление по целевым моментам с **рекурсивной регуляризацией**
  (Malaspinas 2015; Coreixas и др., PRE 96, 033306): ряд Эрмита продолжен до 3-го порядка, а
  коэффициенты 3-го порядка не считаются по популяциям, а выражаются через 2-й рекурсивно:
  `a₃_xxy = 2u_x·a₂_xy + u_y·a₂_xx`, `a₃_xyy = 2u_y·a₂_xy + u_x·a₂_yy`. В D2Q9 `a₃_xxx` и
  `a₃_yyy` решёткой не поддерживаются и отбрасываются; множитель 1/2c_s⁶ (а не 1/6c_s⁶) учитывает
  три перестановки индексов;
- `grad` — то же самое с обрывом ряда на тензоре давлений: условие Града (Dorschner, Bösch,
  Chikatamarla, Boulouchos, Karlin, JFM 801 (2016), разд. 2.1 и прил. B). Задаются не популяции,
  а целевые моменты — ρ, u и Π, — после чего недостающие популяции собираются приближением
  Града (2.13);
- `bouzidi` — интерполированный отскок: доля пересечения q входит в КАЖДУЮ восстанавливаемую
  популяцию, полинково;
- `staircase` — простой отскок, q игнорируется. Не для счёта: это база сравнения, показывающая,
  сколько именно даёт субсеточность.

Целевые моменты у `hrr` и `grad` одни и те же — (B 1) и (B 3) прил. B JFM 801; отличается только
то, до какого порядка восстанавливается функция распределения по этим моментам.

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

| модель | разброс Cd | Cd |
|---|---|---|
| `staircase` | 0.98% | 2.482–2.506 |
| `grad` | 0.62% | 2.449–2.464 |
| `hrr` | 0.61% | 2.450–2.464 |
| `bouzidi` | **0.14%** | 2.459–2.463 |

**HRR не улучшает разрешение геометрии и не должен** — 0.61% против 0.62% у Града. Это следует
из устройства обеих схем: третий порядок Эрмита уточняет ВОССТАНОВЛЕНИЕ популяций по моментам, а
положение стенки входит в моментные схемы совсем другим местом. Обе моментные схемы оказываются
ровно между ступенькой и 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 при τ → ½.

**Умолчание — `hrr`, и это решение временное.** Оно принято по устройству схемы (третий порядок
Эрмита фильтрует высокочастотный мусор у стенки, чего обрыв на Π не делает), а не по здешним
измерениям: на стационарном цилиндре при Re = 20 отличить `hrr` от `grad` нельзя вовсе. Вопрос
ставит ребром группа C кампании — там обе моментные схемы и Bouzidi гоняются на Re = 20, 150 и
2000. Если данные не подтвердят преимущества HRR на турбулентном режиме, умолчанием станет
`bouzidi`, у которого измеренное разрешение геометрии вчетверо лучше.

Все четыре модели работают на обоих бэкендах и совпадают между ними до 0.007% по Cd. Это
специально проверяется: раньше GPU при `--wall staircase` молча считал по Bouzidi, и обнаружилось
это только потому, что две модели дали побитово одинаковый результат там, где обязаны были
разойтись.

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

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

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

**Предел на размер сетки снят.** У GPU есть жёсткий предел `maxComputeWorkgroupsPerDimension`
= 65535, а диспетчеризация была одномерной: при 64 узлах на рабочую группу это упирало сетку в
4.2 миллиона узлов, то есть примерно 2048×2048. Всё, что крупнее, падало ошибкой валидации —
не считало медленно, а не запускалось вовсе. Теперь диспетчеризация двумерная, а линейный
индекс собирается в шейдере (`lin()`/`wlin()`); отображение «рабочая группа → узлы» при этом
остаётся ровно линейным, поэтому редукции ничего не заметили. Проверено до 4096×4096, паритет
с CPU не сдвинулся ни в одной цифре.

**Где именно кончается f32.** Прямой замер на Тейлоре–Грине, где ошибка известна точно:

| N | CPU, f64 | GPU, f32 |
|---|---|---|
| 64 | 9.83e-3 | 9.80e-3 |
| 128 | 2.40e-3 | 3.29e-3 |
| 256 | 5.97e-4 | 2.37e-2 |

Пока истинная ошибка выше ~10⁻³, f32 идёт с f64 вровень; ниже — промахивается на порядок и
больше. Практический вывод, заложенный в кампанию: исследования сходимости считаются на CPU,
всё остальное — на GPU. Двойной точности на GPU здесь быть не может в принципе: **в WGSL типа
`f64` не существует**, поэтому wgpu не выразит её ни на каком железе; локальная Iris Xe вдобавок
сообщает `shaderFloat64 = false`, а на потребительских NVIDIA f64 идёт в 1/64 от f32 — то есть
медленнее, чем CPU.

Что удалось выжать вместо точности — производительность. Два изменения:

1. **Батчинг чтения.** Раньше после каждого шага делался `map_async` + `poll(Wait)` ради 48 байт
   статистики: на 240×120 счёт упирался в 863 шаг/с при том, что сам счёт занимал 0.27 мс из
   1.16. Теперь итоги копятся в кольце на 128 слотов, синхронизация — раз в батч: **6715 шаг/с**,
   в 7.8 раза быстрее, при неизменном пошаговом интерфейсе снаружи.
2. **Компенсированное суммирование** (Кэхена–Ноймайера) в редукциях и в сумме сил. Наивная сумма
   по 10⁵–10⁷ узлам съедает ~log₂N бит мантиссы — именно там f32 терял основную точность.

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

Постановка 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`** — нуль-градиент поперечной скорости; жёсткий ноль отражает вихри
  дорожки обратно к телу.

**Что именно помогает — разделено измерением** (480×240, D=24, Re=150, 30 000 шагов):

| старт | губка, столбцов | пульсация u′/U | St | ⟨Cd⟩ | rms Cl |
|---|---|---|---|---|---|
| uniform | 0 | 1.12% | 0.1839 | 1.483 | 0.404 |
| uniform | 40 | 1.03% | 0.1840 | 1.483 | 0.404 |
| uniform | 120 | 0.83% | 0.1839 | 1.482 | 0.403 |
| rest | 0 | **74.91%** | 0.1534 | 2.009 | 0.807 |
| rest | 40 | **74.90%** | 0.1539 | 2.010 | 0.808 |
| rest | 120 | **74.99%** | 0.1542 | 2.008 | 0.808 |

Читается однозначно: **весь эффект даёт однородный старт**, 74.9% → 1.12%, причём при
полностью выключенной губке. Губка снимает только остаток — 1.12% → 1.03% → 0.83%, — и на
силы с частотой схода не влияет вовсе (St 0.1839 во всех трёх строках).

**При старте из покоя губка не помогает совсем.** Это не осечка реализации, а свойство моды:
губка поднимает вязкость на последних столбцах, а там у стоячей четвертьволновой моды **узел
давления и пучность скорости** — то самое место, где повышенная вязкость её почти не трогает.
Бегущую волну такая губка съедает, стоячую — нет. Убрать моду можно только не возбуждая её.

Заодно видно, ЧЕМ платит неверный старт: St 0.153 вместо 0.184, ⟨Cd⟩ 2.01 вместо 1.48,
rms Cl 0.81 вместо 0.40 — то есть 75-процентная продольная пульсация ломает не косметику, а
все три величины, ради которых постановка и считается.

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

**Инструмент разбора — `--xt`.** Каждые `--xt-every` шагов пишется срез ⟨ρ⟩ и u_x вдоль осевой
линии, по строке на срез. Наклон полос на такой диаграмме прямо даёт скорость распространения:
звук (±c_s) или конвекция (U). Стоячие узлы видны как вертикальные линии постоянной фазы, и это
сразу отличает резонанс от неустойчивости самого граничного условия, привязанной к столбцу x = 0
и никуда не бегущей.

Есть и **губка после входа** — `--sponge-in <столбцов>`, по умолчанию выключена. Гасит продольные
волны до того, как они отразятся от входа-поршня. Осмысленна на высоких Re, где стартовая волна
перестаёт быть безобидной полоской; цена — искажение профиля прямо на входе, поэтому включать её
надо осознанно, а не «на всякий случай».

Что **не помогает** и оставлено только для повторной проверки — `--init-taper`. Замерено:
сглаживание стартовой скорости у тела давит возмущение плотности на первом шаге восьмикратно
(2.54·10⁻² → 3.13·10⁻³), но пик ЗА ПРОГОН при этом даже подрастает (до 3.22·10⁻²). Возмущение
просто переносится во времени: поток всё равно обязан разогнаться вокруг тела, и энергия этого
переходного процесса задана физикой, а не гладкостью начального поля. Умолчание — 0.

## Валидационная кампания

`bench/` — 115 прогонов на ≈90 часов GPU, разложенных по девяти группам: эталоны первоисточников,
цилиндр против литературы, модели стенки, профили крыла, сложная и множественная геометрия,
границы домена, старт и время жизни, внутренние инварианты, сверхмелкие сетки до 4096×2048.
Каждый прогон кладёт логи, ряды, машиночитаемую сводку и гифку на всю свою длительность в
собственную папку.

```sh
cd bench
python preflight.py                # каждый сценарий стартует на два шага: ловит опечатки
./run_campaign.sh --calibrate      # замерить MLUPS этой машины: оценки в часах иначе гадание
./run_campaign.sh --dry-run        # смета: что, сколько шагов, сколько часов
./run_campaign.sh --resume         # считать, пропуская уже готовое
```

Предполётную проверку стоит гонять всерьёз: она поймала, что вся группа сверхмелких сеток
падала на пределе GPU (см. ниже), а девять прогонов передавали `--body-x` дважды. Оба отказа
проявились бы только на сервере, часов через двадцать после старта кампании.

Подробности — в `bench/README.md`: раскладка выходных файлов, таблица групп, модель стоимости и
список того, что известно заранее (какие прогоны обязаны развалиться и почему).

## Развёртывание (Docker)

```sh
docker build -t kbc2d docs/theory/2d_solver
docker run --rm --gpus all kbc2d --calibrate                        # проверить, что GPU виден
docker run -d --gpus all -v "$PWD/out:/work/bench/out" kbc2d --resume
```

`ENTRYPOINT` — драйвер кампании, `CMD` по умолчанию `--dry-run`: случайный `docker run` покажет
смету и выйдет, а не запустит сточасовую задачу.

**Главная тонкость — Vulkan внутри контейнера.** NVIDIA Container Toolkit подкладывает
Vulkan-ICD (`nvidia_icd.json`) только если в `NVIDIA_DRIVER_CAPABILITIES` есть `graphics`;
с одним `compute` wgpu не увидит ни одного адаптера. В образе это прописано, но может быть
переопределено снаружи, поэтому `vulkan-tools` лежит внутрь: первым делом на сервере стоит
выполнить `docker run --rm --gpus all --entrypoint vulkaninfo kbc2d --summary`.

Проверено локально: образ собирается, кампания внутри него проходит смоук с монтированием
результатов на хост, физика совпадает с хостовой до последней цифры (ошибка Тейлора–Грина
9.829e-3 при N=64 и 2.401e-3 при N=128 — те же значения, что вне контейнера), а `--backend gpu`
без проброшенной карты отказывает явным сообщением, а не считает молча. GPU-путь в контейнере
проверяется только на машине с картой.

## Дальше

- σ·n-кросс-чек силы (интеграл тензора напряжений по контуру) как независимая проверка GMEM;
- согласованные начальные условия по образцу работы 2024 года: там ρ и старшие моменты
  получают, решая ∂ρ/∂t + ∇·(ρu₀) = D∇²ρ до стационара, что убрало бы и остаточный
  стартовый импульс от появления тела в потоке;
- подвижные и вращающиеся тела: GMEM уже записан в галилей-инвариантной форме и принимает
  скорость стенки на линке, но подача этой скорости не подключена;
- разбор нерешённых 12% по Cd и предела устойчивости Re ≈ 5·10⁴ — обе задачи вынесены в
  кампанию (группы B и G), выводы делать по её данным;
- больше четырёх тел в домене: сейчас сила считается по четырём вёдрам (`MAX_BODY_BUCKETS`),
  геометрия при этом собирается из любого числа тел, но силы сверх четвёртого сливаются вместе.