1 Commits
Author SHA1 Message Date
NotBigGhostandClaude Opus 5 dd5f519c06 CLAUDE.md вне версионирования: один файл на все ветки
Файл описывает всё дерево целиком, а ветки содержат разные его части
(порт на Rust, исследование CFD, базовый C++-редактор). Отслеживаемая
копия при каждом переключении подменялась бы версией своей ветки, тогда
как нужна одна общая. Содержимое на всех ветках было идентично, так что
расхождений открепление не теряет.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01B9Gcr11JJJyf8NnWsXjDzZ
2026-09-06 16:25:03 +03:00
22 changed files with 0 additions and 15992 deletions
-16
View File
@@ -40,22 +40,6 @@ docs/theory/solver_2x_sdf/out/
# ноутбук kbc_lbm.ipynb их встраивает, и без них он нечитаем на свежем клоне # ноутбук kbc_lbm.ipynb их встраивает, и без них он нечитаем на свежем клоне
# (перегенерация требует GPU-сервера). См. заметку о них в README ветки research. # (перегенерация требует GPU-сервера). См. заметку о них в README ветки research.
# ---------------------------------------------------------------------------
# Rust (docs/theory/2d_solver — решатель KBC на Rust)
# ---------------------------------------------------------------------------
# Каталог сборки cargo: объектники, зависимости, инкрементальный кэш.
target/
# Cargo.lock для бинарного крейта обычно коммитят (воспроизводимость прогонов) —
# он НЕ игнорируется намеренно.
# Артефакты прогонов решателя: гифки, ряды, сводки. Воспроизводятся запуском.
docs/theory/2d_solver/*.gif
docs/theory/2d_solver/*.csv
docs/theory/2d_solver/out/
# Результаты валидационной кампании — гигабайты гифок и рядов. Сам список сценариев
# (bench/scenarios.json) лежит рядом с драйверами и версионируется.
docs/theory/2d_solver/bench/out/
# --------------------------------------------------------------------------- # ---------------------------------------------------------------------------
# Редакторы и ОС # Редакторы и ОС
# --------------------------------------------------------------------------- # ---------------------------------------------------------------------------
Binary file not shown.
-14
View File
@@ -1,14 +0,0 @@
# Каталог сборки cargo — сотни мегабайт, внутри образа он собирается заново.
target/
# Результаты прогонов: гифки и ряды монтируются с хоста, в образ не кладутся.
bench/out/
*.gif
*.csv
*.json
!bench/scenarios.json
# Мусор редакторов и ОС
.vs/
.vscode/
.idea/
.DS_Store
Thumbs.db
-1308
View File
File diff suppressed because it is too large Load Diff
-35
View File
@@ -1,35 +0,0 @@
[package]
name = "kbc2d"
version = "0.1.0"
edition = "2021"
rust-version = "1.75"
description = "Двумерный решатель LBM D2Q9 с энтропийным оператором столкновения KBC (модель D / N1)"
[[bin]]
name = "kbc2d"
path = "src/main.rs"
[dependencies]
clap = { version = "4", features = ["derive"] }
rayon = "1"
gif = "0.13"
bytemuck = { version = "1", features = ["derive"] }
wgpu = { version = "22", optional = true }
pollster = { version = "0.3", optional = true }
[features]
default = ["gpu"]
# GPU-бэкенд на wgpu (Vulkan/DX12/Metal). Отключается: --no-default-features
gpu = ["dep:wgpu", "dep:pollster"]
[profile.release]
opt-level = 3
lto = "thin"
codegen-units = 1
panic = "abort"
# Отдельный профиль для прогонов с проверками переполнения индексов
[profile.release-dbg]
inherits = "release"
debug = true
panic = "unwind"
-79
View File
@@ -1,79 +0,0 @@
# Образ решателя и валидационной кампании.
#
# Поддерживает два способа добраться до видеокарты, и выбирать между ними не нужно —
# доступным окажется ровно один:
#
# 1. НАСТОЯЩИЙ LINUX-СЕРВЕР с драйвером NVIDIA. Vulkan-ICD подкладывает внутрь
# NVIDIA Container Toolkit, и только если в NVIDIA_DRIVER_CAPABILITIES есть
# `graphics` — с одним `compute` wgpu не увидит ни одного адаптера. Запуск через
# docker-compose.server.yml.
#
# 2. WSL2. Драйвер NVIDIA для Linux там отсутствует как класс: карта отдаётся через
# /dev/dxg по протоколу WDDM, а нативный libGLX_nvidia про него не знает и
# enumerate возвращает ноль устройств. Зато с /dev/dxg умеет говорить dzn (Dozen) —
# драйвер Mesa, транслирующий Vulkan в D3D12. Он в образе и лежит. NVIDIA Container
# Toolkit при этом НЕ НУЖЕН: нужен проброс устройства и монтирование /usr/lib/wsl,
# где Microsoft держит libd3d12.so. Запуск через docker-compose.wsl.yml.
#
# Сборка и быстрая проверка:
# docker build -t kbc2d docs/theory/2d_solver
# docker compose -f docker-compose.wsl.yml --profile check run --rm vulkan
# ── сборка ───────────────────────────────────────────────────────────────────
# Версия тулчейна закреплена: та же, на которой решатель собирался и проверялся.
FROM rust:1.97-bookworm AS build
WORKDIR /src
# Сначала только манифесты — тогда слой с зависимостями переиспользуется, пока они не менялись.
COPY Cargo.toml Cargo.lock ./
RUN mkdir src && echo 'fn main() {}' > src/main.rs && cargo build --release --locked 2>/dev/null || true && rm -rf src
COPY src ./src
# touch нужен, чтобы cargo не спутал новые исходники с заглушкой из слоя выше
RUN touch src/main.rs && cargo build --release --locked
# ── исполнение ───────────────────────────────────────────────────────────────
# База Arch, а не debian-slim, по одной причине: dzn собирают не все дистрибутивы.
# В mesa-vulkan-drivers Debian и Ubuntu его нет (проверено по списку файлов пакета),
# в Arch он лежит отдельным пакетом vulkan-dzn той же версии Mesa, что и на хосте WSL.
FROM archlinux:base
ARG VERSION=dev
ARG REVISION=unknown
LABEL org.opencontainers.image.title="kbc2d"
LABEL org.opencontainers.image.description="Двумерный решатель LBM D2Q9 с энтропийным столкновением KBC"
LABEL org.opencontainers.image.version="${VERSION}"
LABEL org.opencontainers.image.revision="${REVISION}"
LABEL org.opencontainers.image.source="https://gitea.arseniev.info/NotBigGhost/TurbulenceCAD"
# Справка, руководства и переводы занимают в базе Arch заметно больше, чем сам
# драйвер, а в контейнере не нужны никому. NoExtract прописывается ДО установки.
RUN { echo 'NoExtract = usr/share/man/*'; \
echo 'NoExtract = usr/share/doc/*'; \
echo 'NoExtract = usr/share/info/*'; \
echo 'NoExtract = usr/share/locale/*'; \
echo 'NoExtract = usr/share/i18n/*'; } >> /etc/pacman.conf \
&& pacman -Syu --noconfirm --needed \
vulkan-icd-loader vulkan-dzn vulkan-tools python \
&& pacman -Scc --noconfirm \
&& rm -rf /var/cache/pacman/pkg/* /var/lib/pacman/sync/* \
/usr/share/man /usr/share/doc /usr/share/info /usr/share/locale
# Проверка на этапе сборки: без dzn образ бесполезен для WSL, и узнать об этом надо
# здесь, а не через час после старта кампании.
RUN ls /usr/share/vulkan/icd.d/ && ls /usr/share/vulkan/icd.d/dzn_icd*.json > /dev/null
# Для пути 1 (настоящий сервер). На WSL эти переменные ни на что не влияют.
ENV NVIDIA_VISIBLE_DEVICES=all
ENV NVIDIA_DRIVER_CAPABILITIES=compute,utility,graphics
WORKDIR /work
COPY --from=build /src/target/release/kbc2d /work/target/release/kbc2d
COPY bench /work/bench
RUN chmod +x /work/bench/run_campaign.sh
# Кампания идёт десятки часов, поэтому по умолчанию образ НЕ запускает её: случайный
# `docker run` покажет смету и выйдет.
WORKDIR /work/bench
ENTRYPOINT ["./run_campaign.sh"]
CMD ["--dry-run"]
-609
View File
@@ -1,609 +0,0 @@
# 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)
Собранный образ опубликован: **`notbigghost/kbc2d:1.2.0`** (он же `latest`, платформа `linux/amd64`).
Исходники на сервере не нужны — достаточно перенести туда один файл
`docker-compose.server.yml`:
```sh
mkdir -p ~/kbc2d && cd ~/kbc2d # сюда же ляжет ./out с результатами
# перенести docker-compose.server.yml
docker compose -f docker-compose.server.yml --profile check run --rm vulkan # карта видна?
docker compose -f docker-compose.server.yml --profile check run --rm preflight # сценарии стартуют?
docker compose -f docker-compose.server.yml --profile check run --rm calibrate # сколько MLUPS?
docker compose -f docker-compose.server.yml up -d # кампания
docker compose -f docker-compose.server.yml logs -f
```
Порядок именно такой: узнать, что карта не видна, лучше через минуту, чем через час. Замеренные
`--calibrate` числа подставляются переменными `KBC2D_GPU_MLUPS` / `KBC2D_CPU_MLUPS` — только на
оценки в часах, на счёт они не влияют.
Собрать образ самому (`docker-compose.yml` рядом делает то же самое с `build:`):
```sh
docker build -t kbc2d docs/theory/2d_solver
docker run --rm --gpus all kbc2d --calibrate
```
`ENTRYPOINT` — драйвер кампании, `CMD` по умолчанию `--dry-run`: случайный `docker run` покажет
смету и выйдет, а не запустит сточасовую задачу. В серверном compose политика перезапуска —
`on-failure`, а не `unless-stopped`: кампания завершается штатно с кодом 0, и «перезапускать
всегда» крутило бы контейнер вхолостую по кругу, тогда как падение или перезагрузку хоста
`on-failure` подхватывает, а `--resume` продолжает с места.
**Главная тонкость — 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`
без проброшенной карты отказывает явным сообщением, а не считает молча.
### WSL2 — отдельный путь
В WSL2 всё вышеописанное не работает, и не из-за настроек. **Драйвера Vulkan для Linux у
NVIDIA там нет**: карта отдаётся через `/dev/dxg` по протоколу WDDM, нативный
`libGLX_nvidia` про него не знает и перечисляет ноль устройств. Container Toolkit
подкладывать внутрь нечего, отсюда `could not select device driver "nvidia"`.
Работает другое: **dzn** (Dozen) — драйвер Mesa, транслирующий Vulkan в D3D12, он умеет
говорить с `/dev/dxg` напрямую. Он положен в образ (ради него база сменена с
`debian:bookworm-slim` на `archlinux:base`: в пакетах Mesa у Debian и Ubuntu dzn не
собирают). NVIDIA-runtime для этого пути не нужен вовсе — нужны проброс устройства и
монтирование `/usr/lib/wsl`. Запуск через `docker-compose.wsl.yml`, подробности в
[`bench/README.md`](bench/README.md).
Три вещи, которые надо знать про этот путь.
**wgpu по умолчанию прячет несоответствующие адаптеры.** dzn сообщает о себе
`conformanceVersion = 0.0.0.0`, и wgpu молча его отбрасывает — решатель докладывает, что
GPU не найден. Согласие даётся явно, переменной `WGPU_ALLOW_UNDERLYING_NONCOMPLIANT_ADAPTER=1`;
чтобы она вообще читалась, в `gpu.rs` при создании инстанса стоит
`InstanceFlags::from_build_config().with_env()`. Поведение по умолчанию не изменилось: без
переменной несоответствующие адаптеры по-прежнему скрыты.
**У dzn мал предел на размер одной привязки — 128 МиБ** против гигабайтов у нативных
драйверов. Массив популяций занимает `nx · ny · 9 · 4` байта, поэтому при одной общей
привязке потолок выходил 3.73 млн узлов, и сетки от 2048×2048 не запускались вовсе. Но
ограничена именно привязка, а не буфер (`max_buffer_size` у dzn 2047 МиБ), так что тот же
буфер теперь показывается **девятью привязками по одному направлению**: потолок
поднимается до 33.5 млн узлов. Вариант выбирается по возможностям адаптера; где предела
нет, собирается прежний общий, без `switch` в аксессорах. Оба дают одинаковые числа
(Cd 2.39486 в обоих), раздельный стоит 5.7% пропускной способности. Принудительно
включается переменной `KBC2D_SPLIT_POPULATIONS=1` — она нужна для сверки двух вариантов
на одной карте.
**Точность трансляция не портит.** Замерено на Intel Iris Xe одним и тем же прогоном
(`bench/parity.py`, цилиндр Re=20, 8000 шагов):
| | Cd | energy_end (Тейлор–Грин) |
|---|---|---|
| CPU, f64 | 2.39490 | 5.74754e-05 |
| нативный Vulkan, f32 | 2.39486 | 5.74813e-05 |
| dzn в контейнере, f32 | 2.39486 | 5.74810e-05 |
Числа dzn и нативного драйвера сходятся до 5–6 значащих цифр, и разница между ними меньше,
чем между любым из них и f64. Считает трансляция то же самое.
**Скорость она портит, но тем меньше, чем крупнее сетка** — плата почти вся приходится на
трансляцию вызова, а не счёта:
| сетка | узлов | нативно | dzn в контейнере | плата |
|---|---|---|---|---|
| 320×192 | 61 тыс. | 188.8 MLUPS | 44.5 MLUPS | 4.2× |
| 960×480 | 461 тыс. | 125.3 MLUPS | 86.0 MLUPS | 1.46× |
| 1920×960 | 1.84 млн | 128.7 MLUPS | 98.7 MLUPS | 1.30× |
Для кампании это решает дело: 95.1% её стоимости приходится на сетки крупнее 600 тыс.
узлов, а на сетки мельче 150 тыс. — 0.0%. Ожидаемое удорожание всей кампании в контейнере —
около трети, а не в разы. Проверяется профилем `calibrate`, который меряет в том числе
1920×960.
## Дальше
- σ·n-кросс-чек силы (интеграл тензора напряжений по контуру) как независимая проверка GMEM;
- согласованные начальные условия по образцу работы 2024 года: там ρ и старшие моменты
получают, решая ∂ρ/∂t + ∇·(ρu₀) = D∇²ρ до стационара, что убрало бы и остаточный
стартовый импульс от появления тела в потоке;
- подвижные и вращающиеся тела: GMEM уже записан в галилей-инвариантной форме и принимает
скорость стенки на линке, но подача этой скорости не подключена;
- разбор нерешённых 12% по Cd и предела устойчивости Re ≈ 5·10⁴ — обе задачи вынесены в
кампанию (группы B и G), выводы делать по её данным;
- больше четырёх тел в домене: сейчас сила считается по четырём вёдрам (`MAX_BODY_BUCKETS`),
геометрия при этом собирается из любого числа тел, но силы сверх четвёртого сливаются вместе.
-248
View File
@@ -1,248 +0,0 @@
# Валидационная кампания
115 прогонов, ≈90 часов на RTX 4070 Ti. Проверяет решатель по трём независимым линиям:
эталонам из статей авторов метода, литературе по обтеканию тел и внутренним инвариантам самой
схемы. Каждый прогон кладёт логи, ряды, машиночитаемую сводку и гифку в собственную папку.
## Быстрый старт на сервере
Образ опубликован, собирать ничего не нужно: **`notbigghost/kbc2d:1.2.0`**. Исходники на
сервере тоже не нужны — переносится один файл `docker-compose.server.yml`.
```sh
mkdir -p ~/kbc2d && cd ~/kbc2d # сюда же ляжет ./out с результатами
# перенести сюда docker-compose.server.yml
C=docker-compose.server.yml
docker compose -f $C --profile check run --rm vulkan # 1. карта видна из контейнера?
docker compose -f $C --profile check run --rm preflight # 2. все 115 сценариев стартуют?
docker compose -f $C --profile check run --rm calibrate # 3. сколько MLUPS на этой машине?
docker compose -f $C --profile check run --rm plan # 4. смета в часах по замеренному
docker compose -f $C up -d # 5. кампания
docker compose -f $C logs -f
```
Порядок не случайный: узнать, что карта не видна, лучше на первом шаге, чем через час счёта.
Замеренные калибровкой числа подставляются переменными окружения — они влияют только на
оценки в часах, не на счёт:
```sh
KBC2D_GPU_MLUPS=1450 KBC2D_CPU_MLUPS=40 docker compose -f $C --profile check run --rm plan
```
**Если `vulkaninfo` показывает ноль адаптеров** — почти наверняка дело в
`NVIDIA_DRIVER_CAPABILITIES`: Container Toolkit подкладывает Vulkan-ICD только при наличии
`graphics` в списке. В образе и в compose это прописано, но может быть переопределено снаружи.
Проверить, что хост вообще умеет отдавать карту:
```sh
docker run --rm --gpus all nvidia/cuda:12.4.0-base-ubuntu22.04 nvidia-smi
```
## Запуск в WSL2
Отдельный путь, потому что в WSL2 **драйвера Vulkan для Linux у NVIDIA не существует**.
Карта отдаётся через `/dev/dxg` по протоколу WDDM; нативный `libGLX_nvidia` про него не
знает и возвращает ноль устройств — загрузчик его выбрасывает. Отсюда и `could not select
device driver "nvidia"`: NVIDIA Container Toolkit тут не поможет, потому что подкладывать
внутрь нечего.
Зато с `/dev/dxg` умеет говорить **dzn** (Dozen) — драйвер Mesa, транслирующий Vulkan в
D3D12. Он лежит в образе. NVIDIA-runtime при этом **не нужен вовсе**: достаточно проброса
устройства и монтирования `/usr/lib/wsl`, где Microsoft держит `libd3d12.so`.
```sh
C=docker-compose.wsl.yml
docker compose -f $C --profile check run --rm vulkan # 1. карта видна?
docker compose -f $C --profile check run --rm parity # 2. считает ли она правильно?
docker compose -f $C --profile check run --rm calibrate # 3. и с какой скоростью?
docker compose -f $C --profile check run --rm preflight # 4. все сценарии стартуют?
docker compose -f $C --profile check run --rm plan # 5. смета в часах
docker compose -f $C up -d # 6. кампания
```
Если образа нет ни локально, ни в реестре — `docker compose -f $C build`, доступ к Docker
Hub не обязателен.
### Шаг parity обязателен
dzn сообщает о себе `conformanceVersion = 0.0.0.0`: набор тестов соответствия Vulkan он не
проходил. wgpu по этой причине по умолчанию **прячет** такие адаптеры, и решатель сообщает,
что GPU не найден. Согласие считать на непроверенном драйвере даётся явно — переменной
`WGPU_ALLOW_UNDERLYING_NONCOMPLIANT_ADAPTER=1`, она прописана в `docker-compose.wsl.yml`.
Раз соответствие не проверено вендором, его проверяем сами. `bench/parity.py` гоняет два
коротких эталона на CPU в f64 и на GPU и сверяет числа: вихрь Тейлора–Грина (есть точное
решение) и стационарное обтекание цилиндра при Re=20. Течения выбраны намеренно не
хаотические — там расхождение f32 и f64 не нарастает, поэтому заметное различие означает
проблему драйвера, а не разрядности.
Замерено на Intel Iris Xe: **dzn совпадает с нативным драйвером той же карты до 5–6
значащих цифр** (`energy_end` 3.02277e-07 против 3.02283e-07, `cd` 2.39486 против 2.39486).
Точность трансляция не портит. Но это замер на Intel; на своей карте прогоните сами — две
минуты.
### Предел 128 МиБ на привязку и как он снят
dzn объявляет `max_storage_buffer_binding_size = 128 МиБ` (2²⁷ байт), тогда как нативные
драйверы дают гигабайты. Массив популяций занимает `nx · ny · 9 · 4` байта, так что при
одной общей привязке потолок выходил **3.73 млн узлов** (примерно 1920×1920), и 14
прогонов кампании из 115 падали на создании bind group:
```
Buffer binding 0 range 150994944 exceeds `max_*_buffer_binding_size` limit 134217728
```
На эти 14 приходилось 70.7% стоимости кампании — дорогие прогоны как раз крупносеточные.
Существенно, что ограничена только **привязка**: `max_buffer_size` у dzn 2047 МиБ, то есть
буфер держать разрешено, нельзя лишь показать шейдеру его целиком. Поэтому решатель теперь
умеет показывать тот же буфер **девятью привязками**, по одному направлению в каждой.
Потолок поднимается в девять раз — до 33.5 млн узлов, чего хватает всей кампании с запасом
(самая крупная сетка в ней 4096×4096 — 16.8 млн узлов).
Вариант выбирается сам, по `max_storage_buffer_binding_size` адаптера: где предела нет,
собирается прежний общий вариант без `switch` в аксессорах. Проверить оба на одной карте
можно переменной `KBC2D_SPLIT_POPULATIONS=1` — она включает раздельные привязки
принудительно.
Замерено на Intel Iris Xe, один драйвер, два варианта привязки:
| | Cd | energy_end | MLUPS (цилиндр) |
|---|---|---|---|
| общая привязка | 2.39486 | 5.74813e-05 | 169.7 |
| девять привязок | 2.39486 | 5.74813e-05 | 160.0 |
Числа совпадают полностью; раздельный вариант стоит **5.7%** пропускной способности, и
включается только там, где без него счёт вообще невозможен.
### Чего трансляция стоит по скорости
Плата есть, но она почти вся — накладные расходы на вызов, а не на счёт, и потому падает
с ростом сетки. Замерено на Intel Iris Xe, один и тот же решатель:
| сетка | узлов | нативный Vulkan | dzn в контейнере | плата |
|---|---|---|---|---|
| 320×192 | 61 тыс. | 188.8 MLUPS | 44.5 MLUPS | **4.2×** |
| 960×480 | 461 тыс. | 125.3 MLUPS | 86.0 MLUPS | **1.46×** |
| 1920×960 | 1.84 млн | 128.7 MLUPS | 98.7 MLUPS | **1.30×** |
| 2048×2048 | 4.19 млн | 111.8 MLUPS | 67.4 MLUPS | **1.66×** |
Последняя строка стоит особняком: там уже раздельные привязки (см. ниже), и в плату
входит их `switch` в аксессорах. Ожидаемый разброс по кампании — от 1.3× до 1.7×.
Каждый шаг решателя — несколько отправок в очередь; на мелкой сетке трансляция вызова
стоит дороже самого счёта, на крупной размазывается. Для кампании это решающее
обстоятельство, потому что дорогие прогоны в ней как раз крупные:
| узлов в сетке | прогонов | доля стоимости кампании |
|---|---|---|
| < 150 тыс. | 18 | 0.0% |
| 150–600 тыс. | 47 | 4.9% |
| > 600 тыс. | 50 | **95.1%** |
То есть 95% времени кампания проводит там, где плата 1.3–1.5×, и почти не бывает там, где
она четырёхкратна. Ожидаемое удорожание по всей кампании — около трети: 90 часов
превращаются примерно в 115–120, а не в 400.
Проверьте на своей карте: профиль `calibrate` меряет три сетки, и ориентироваться надо на
**1920×960** — она представительна для кампании, а 240×120 показывает худший случай.
## Своя сборка образа
Если нужен образ из текущего состояния репозитория, а не опубликованный:
```sh
docker build -t kbc2d docs/theory/2d_solver --build-arg VERSION=dev --build-arg REVISION=$(git rev-parse --short HEAD)
```
Рядом лежит `docker-compose.yml` — то же самое через `build:`, для локальной отладки.
## Без Docker
```sh
cargo build --release # из docs/theory/2d_solver
cd bench
python preflight.py # каждый сценарий стартует на два шага
./run_campaign.sh --calibrate # замерить MLUPS этой машины
KBC2D_GPU_MLUPS=1400 ./run_campaign.sh --dry-run
./run_campaign.sh --resume
```
Под Windows то же самое делает `run_campaign.ps1` (тот же `scenarios.json`); он нужен для
локальной отладки обвязки, целевая машина — Linux.
## Ключи драйвера
| ключ | что делает |
|---|---|
| `--dry-run` | печатает смету: что, сколько шагов, сколько часов. Ничего не считает |
| `--calibrate` | три коротких прогона, замер фактических MLUPS этой машины |
| `--smoke` | три самых дешёвых прогона: проверить обвязку, а не физику |
| `--resume` | пропускает прогоны, у которых уже есть `summary.json` |
| `--group A,B` | только выбранные группы |
| `--only cyl_re150` | по подстроке идентификатора |
| `--budget-hours 24` | остановиться, когда время выйдет |
Прогон, который упал или развалился, помечается в сводке и **не останавливает кампанию**:
группы E и G специально ищут предел устойчивости, там развал — ожидаемый результат.
## Что где лежит
```
out/summary.csv сводная таблица: id, группа, статус, секунды, стоимость
out/campaign.log журнал с отметками времени
out/<id>/cmd.txt точная команда, которой прогон был запущен
out/<id>/log.txt полный вывод
out/<id>/report.txt только итоговый отчёт
out/<id>/series.csv временные ряды по шагам
out/<id>/summary.json ключевые метрики машиночитаемо — с этого удобно начинать разбор
out/<id>/*.gif анимация
out/<id>/*_case.csv энергия, энстрофия, палинстрофия (эталонные течения)
out/<id>/xt_*.csv x–t диаграммы (группа G)
```
## Группы
| | прогонов | что проверяется |
|---|---|---|
| **A** | 12 | эталоны первоисточников: Тейлор–Грин (второй порядок сходимости; при фиксированном u₀ — полка O(Ma²)), сдвиговый слой Re=3·10⁴, затухающая турбулентность |
| **B** | 16 | цилиндр против литературы: стационар Re=20/40, дорожка Re=100…300, разрешение D=16…128, блокировка Ny/D=6…32 с экстраполяцией к бесконечной среде |
| **C** | 20 | модели стенки и субсеточность: hrr / grad / bouzidi / staircase при Re=150 и 2000, плюс сдвиг тела внутри клетки на 0 / ¼ / ½ для всех четырёх |
| **D** | 14 | профили крыла: поляра NACA 0012, Re-серия, изгиб 4412, сходимость по хорде |
| **E** | 12 | сложная и множественная геометрия: тандем, решётка, перфорация, многоэлементный профиль, сверхтонкая пластина, клин, зазубренная кромка |
| **F** | 9 | поле влияния и границы домена: отступы до входа и выхода, вложенный патч |
| **G** | 14 | старт, акустика, время жизни: x–t диаграммы, губки, предел по Re, прогон на 10⁷ шагов |
| **H** | 9 | инварианты: симметрия, зеркальность, зависимость от числа Маха, расхождение f32 против f64 |
| **I** | 9 | сверхмелкие сетки 4096×2048 по всем формам и композиции из трёх тел |
## Стоимость и время
Стоимость каждого прогона хранится в **обновлениях узлов** (`nodes_per_step × steps`) — это
единственная мера, переносимая между машинами. Часы драйвер получает, поделив её на MLUPS.
Оценки по умолчанию исходят из 1200 MLUPS на GPU и 22 на CPU; **`--calibrate` обязателен**,
потому что эти числа взяты с другой машины.
Пересобрать список с другой длительностью:
```sh
python gen_scenarios.py --scale 3 # все прогоны в полтора раза длиннее (≈135 ч)
```
## Что известно заранее
- **Точность f32.** Замерено на Тейлоре–Грине: GPU совпадает с CPU/f64, пока истинная ошибка
выше ~10⁻³, и промахивается в 40 раз, когда она ниже. Поэтому исследования сходимости из
группы A идут на CPU — это указано в самих сценариях.
- **Развалы ожидаемы** в группе G (предел по Re) и у прогона `A09_shear_n512_lbgk`: LBGK при
Re=3·10⁴ обязан развалиться там, где KBC доживает — это и есть проверяемое утверждение.
- **Гифки** пишутся на всю длительность прогона, в реальном времени, 10 кадр/с. Число кадров
этим задано жёстко (у самого длинного прогона их 16 666), поэтому единственный рычаг —
размер кадра. Замерено: 0.103 байта на пиксель после LZW; отсюда бюджет
кадры×пиксели ≤ 1.5·10⁹ на гифку, но ширина не опускается ниже 480 пикселей. Итог по
кампании: **≈4.4 ГБ**, самый тяжёлый файл 226 МБ.
- **Перед запуском** имеет смысл прогнать предполётную проверку: каждый сценарий стартует на
два шага, что ловит опечатки в ключах и несовместимые сочетания до того, как кампания уйдёт
считать на девяносто часов. Именно она поймала, что вся группа I падала на пределе GPU в
65535 рабочих групп на измерение.
@@ -1,518 +0,0 @@
#!/usr/bin/env python3
# Генератор списка прогонов кампании -> scenarios.json
#
# Список порождается кодом, а не правится руками: серии по Re, по разрешению и по углу атаки
# получаются циклами, а стоимость каждого прогона считается тут же и суммируется. Запуск:
# python gen_scenarios.py # перезаписать scenarios.json и напечатать смету
# python gen_scenarios.py --scale 2 # растянуть все длительности вдвое
#
# Стоимость меряется в ОБНОВЛЕНИЯХ УЗЛОВ (nodes_per_step * steps) — единственная переносимая
# между машинами мера. Часы драйвер получает из неё, поделив на фактические MLUPS, которые
# замеряет на месте (--calibrate).
import argparse
import json
import math
import os
RUNS = []
U_LAT = 0.05 # решёточная скорость по умолчанию
DX = 0.1 # размер клетки, м (умолчание решателя)
U_PHYS = 30.0 # скорость потока, м/с (умолчание решателя)
GIF_MAX_W = 1280 # предельная ширина кадра после прореживания
GIF_MIN_W = 480 # ниже не опускаемся: нечитаемая гифка бесполезнее большой
GIF_PX_BUDGET = 1.5e9 # кадры×пиксели на гифку; при 0.103 байта на пиксель это ~150 МБ
def conv_steps(d_cells, n_conv, u_lat=U_LAT):
"""Сколько шагов нужно на n_conv конвективных времён D/U."""
return int(round(n_conv * d_cells / u_lat))
def add(rid, group, title, expect, args, nodes, steps, backend="gpu",
gif=None, nx=None, extra_out=()):
"""Записать прогон. `nodes` — обновлений узлов за шаг, `steps` — сколько шагов."""
a = list(args) + ["--steps", str(steps), "--backend", backend]
# ряды прореживаем так, чтобы их осталось порядка 200 тысяч записей
if steps > 400_000:
a += ["--series-every", str(max(1, steps // 200_000))]
if "case_csv" in extra_out:
a += ["--case-csv", f"{rid}_case.csv"]
if gif:
# Гифка идёт на всю длительность прогона в реальном времени, 10 кадр/с — число кадров
# задано физикой прогона и не обсуждается. Единственный рычаг — пиксели в кадре, и он
# нужен: замерено 0.103 байта на пиксель после LZW, так что прогон на 16 тысяч кадров
# при 960×576 весит 900 МБ, а такой файл уже нечем открыть. Поэтому кроме ограничения
# по ширине действует бюджет ПРОИЗВЕДЕНИЯ кадры×пиксели, но ширина не опускается ниже
# GIF_MIN_W: нечитаемая гифка бесполезнее большой.
ny = max(1, int(nodes) // max(1, nx or 1))
u_lat = float(a[a.index("--u-lat") + 1]) if "--u-lat" in a else U_LAT
frames = max(1, int(steps * (u_lat * DX / U_PHYS) * 10))
down = max(1, math.ceil((nx or 0) / GIF_MAX_W))
while nx and frames * (nx // down) * (ny // down) > GIF_PX_BUDGET \
and nx // (down + 1) >= GIF_MIN_W:
down += 1
a += ["--gif", gif, "--gif-field", "vorticity", "--gif-every", "auto",
"--gif-fps", "10", "--gif-speed", "1",
"--gif-downsample", str(down), "--gif-scale", "1"]
RUNS.append({
"id": rid, "group": group, "title": title, "expect": expect,
"backend": backend, "args": a,
"nodes_per_step": int(nodes), "steps": int(steps),
"cost": int(nodes) * int(steps),
"outputs": list(extra_out),
})
# ══════════════════════════════════════════════════════════════════════════════
# A. Эталоны первоисточников: периодические течения, тела и ГУ нет вовсе
# ══════════════════════════════════════════════════════════════════════════════
def group_a(scale):
NU = 0.0192
for n in (64, 128, 256, 512):
# диффузионное измельчение: nu фиксирована, u0 ~ 1/N, поэтому Re сохраняется, а
# число Маха падает вместе с сеткой — только так виден второй порядок целиком
u0 = 0.03 * 64 / n
tc = int(round(math.log(2.0) / (NU * (2 * math.pi / n) ** 2 * 17)))
add(f"A{len(RUNS)+1:02d}_tg_diff_n{n}", "A",
f"Тейлор–Грин, диффузионное измельчение, N={n} (CPU/f64)",
"порядок сходимости 2 по ряду N=64..512; ошибка на N=512 ниже 2e-4",
["--case", "taylor-green", "--nx", str(n), "--ny", str(n), "--refine", "1",
"--u-lat", f"{u0:.6f}", "--re", "100", "--case-every", str(max(1, tc // 20))],
n * n, tc, backend="cpu", extra_out=("case_csv",))
for n in (128, 256):
# при ФИКСИРОВАННОМ u0 ошибка упирается в полку O(Ma^2) и от сетки не зависит
tc = int(round(math.log(2.0) / (NU * (2 * math.pi / n) ** 2 * 17)))
add(f"A{len(RUNS)+1:02d}_tg_fixed_n{n}", "A",
f"Тейлор–Грин при фиксированном u₀, N={n} — демонстрация полки O(Ma²)",
"ошибка перестаёт падать с ростом N: это свойство слабо-сжимаемого метода",
["--case", "taylor-green", "--nx", str(n), "--ny", str(n), "--refine", "1",
"--u-lat", "0.03", "--re", f"{0.03 * n / NU:.1f}",
"--case-every", str(max(1, tc // 20))],
n * n, tc, backend="cpu", extra_out=("case_csv",))
for n in (512,):
for op, km, tag in (("kbc", "n1", "kbc_n1"), ("kbc", "n2", "kbc_n2"), ("bgk", "n1", "lbgk")):
steps = int(4 * n / 0.04 * scale)
add(f"A{len(RUNS)+1:02d}_shear_n{n}_{tag}", "A",
f"Сдвиговый слой Re=3·10⁴, N={n}, {tag}",
"KBC доживает до конца, LBGK обязан развалиться (разд. VII статьи)",
["--case", "shear-layer", "--nx", str(n), "--ny", str(n), "--refine", "1",
"--u-lat", "0.04", "--re", "30000", "--collision", op, "--kbc-model", km,
"--case-every", str(max(1, steps // 200))],
n * n, steps, gif=f"shear_{tag}.gif", nx=n, extra_out=("case_csv",))
for n, tag, ops in ((2048, "n2048", ("kbc", "bgk")), (4096, "n4096", ("kbc",))):
for op in ops:
steps = int(200_000 * scale * (2048 / n))
add(f"A{len(RUNS)+1:02d}_turb_{tag}_{op}", "A",
f"Затухающая турбулентность N={n}, {op}",
"энстрофия и палинстрофия падают монотонно; KBC держится там, где LBGK нет",
["--case", "decaying-turbulence", "--nx", str(n), "--ny", str(n), "--refine", "1",
"--u-lat", "0.02", "--re", f"{0.02 * n / 8.16e-5:.0f}", "--collision", op,
"--case-every", str(max(1, steps // 300))],
n * n, steps, gif=f"turb_{tag}_{op}.gif", nx=n, extra_out=("case_csv",))
# ══════════════════════════════════════════════════════════════════════════════
# B. Цилиндр против литературы
# ══════════════════════════════════════════════════════════════════════════════
def flatten(d):
"""Словарь ключей в плоский argv. Именно словарь, а не список: clap отвергает
повторённый ключ, а сценарии сплошь и рядом переопределяют умолчания сборщика."""
a = []
for k, v in d.items():
a += [k, str(v)]
return a
def cyl_args(d, re, nx, ny, **kw):
a = {"--shape": "cylinder", "--size": d, "--nx": nx, "--ny": ny,
"--body-x": nx // 5, "--re": re, "--refine": 1}
a.update({"--" + k.replace("_", "-"): v for k, v in kw.items()})
return flatten(a)
def group_b(scale):
# стационарные режимы: дорожки нет, Cd — одно число
for re in (20, 40):
d, nx, ny = 48, 960, 576
st = int(conv_steps(d, 60) * scale)
add(f"B{len(RUNS)-len(RUNS)+len([r for r in RUNS if r['group']=='B'])+1:02d}_cyl_re{re}",
"B", f"Цилиндр Re={re}, стационар, D=48",
"Cd выходит на постоянную; сравнение с литературой по стационарному обтеканию",
cyl_args(d, re, nx, ny, pert_amp=0), nx * ny, st,
gif=f"cyl_re{re}.gif", nx=nx)
# дорожка Кармана: длинное осреднение, это главный заход на литературу
for re in (100, 150, 200, 300):
d, nx, ny = 64, 1280, 768
st = int(conv_steps(d, 1500) * scale)
add(f"B{len([r for r in RUNS if r['group']=='B'])+1:02d}_cyl_re{re}", "B",
f"Цилиндр Re={re}, дорожка Кармана, D=64, 1500 конв. времён",
"St≈0.183 и ⟨Cd⟩≈1.33 при Re=150 после поправки на блокировку",
cyl_args(d, re, nx, ny), nx * ny, st, gif=f"cyl_re{re}.gif", nx=nx)
# разрешение тела: сходимость Cd
for d in (16, 32, 64, 128):
nx, ny = 20 * d, 12 * d
st = int(conv_steps(d, 800) * scale)
add(f"B{len([r for r in RUNS if r['group']=='B'])+1:02d}_cyl_d{d}", "B",
f"Цилиндр Re=150, разрешение D={d}",
"Cd сходится по D; экстраполяция даёт сеточно-независимое значение",
cyl_args(d, 150, nx, ny), nx * ny, st, gif=f"cyl_d{d}.gif", nx=nx)
# блокировка канала: экстраполяция к бесконечной среде — прямой заход на нерешённые 12%
for k in (6, 8, 12, 16, 24, 32):
d, nx, ny = 48, 960, k * 48
st = int(conv_steps(d, 1000) * scale)
add(f"B{len([r for r in RUNS if r['group']=='B'])+1:02d}_cyl_block{k}", "B",
f"Цилиндр Re=150, блокировка Ny/D={k}",
"экстраполяция Cd и St к нулевой блокировке; ожидание — сходимость к 1.33 и 0.183",
cyl_args(d, 150, nx, ny), nx * ny, st, gif=None, nx=nx)
# ══════════════════════════════════════════════════════════════════════════════
# C. Модели стенки и субсеточность
# ══════════════════════════════════════════════════════════════════════════════
def group_c(scale):
for wall in ("hrr", "grad", "bouzidi", "staircase"):
for re in (150, 2000):
d, nx, ny = 32, 640, 384
st = int(conv_steps(d, 400 if re > 20 else 60) * scale)
add(f"C{len([r for r in RUNS if r['group']=='C'])+1:02d}_wall_{wall}_re{re}", "C",
f"Модель стенки {wall}, Re={re}",
"все субсеточные модели обязаны сойтись к одному пределу; staircase — база",
cyl_args(d, re, nx, ny, wall=wall), nx * ny, st,
gif=f"wall_{wall}_re{re}.gif" if re == 150 else None, nx=nx)
# Чувствительность к положению тела ВНУТРИ клетки — прямая проверка субсеточности, и
# единственное измерение, которое вообще разделяет модели стенки. Двух положений для этого
# мало: по двум точкам не отличить систематический разброс от совпадения, поэтому берём
# четверти, и все четыре модели. Прогоны короткие (стационар при Re=20), вся серия — минуты.
for wall in ("hrr", "grad", "bouzidi", "staircase"):
for off in (0.0, 0.25, 0.5):
d, nx, ny = 16, 320, 192
st = int(conv_steps(d, 60) * scale)
add(f"C{len([r for r in RUNS if r['group']=='C'])+1:02d}_sub_{wall}_{int(off*100):02d}",
"C", f"Субклеточный сдвиг тела {off} клетки, стенка {wall}",
"разброс Cd по сдвигам: чем субсеточнее модель, тем он меньше",
cyl_args(d, 20, nx, ny, body_y=96 + off, wall=wall, pert_amp=0),
nx * ny, st)
# ══════════════════════════════════════════════════════════════════════════════
# D. Профили крыла
# ══════════════════════════════════════════════════════════════════════════════
def foil_args(chord, alpha, re, nx, ny, naca="0012", **kw):
a = {"--shape": "naca", "--naca": naca, "--size": chord,
"--body-angle": alpha, "--nx": nx, "--ny": ny,
"--body-x": nx // 4, "--re": re, "--refine": 1}
a.update({"--" + k.replace("_", "-"): v for k, v in kw.items()})
return flatten(a)
def group_d(scale):
# поляра: наклон Cl(alpha) и положение сваливания
for al in (0, 4, 8, 12, 16):
c, nx, ny = 128, 2048, 1024
st = int(conv_steps(c, 400) * scale)
add(f"D{len([r for r in RUNS if r['group']=='D'])+1:02d}_naca0012_a{al:02d}", "D",
f"NACA 0012, α={al}°, Re=1000, хорда 128",
"линейный участок Cl(α) и срыв; сравнение с низкорейнольдсовой литературой",
foil_args(c, al, 1000, nx, ny), nx * ny, st,
gif=f"naca0012_a{al:02d}.gif", nx=nx)
for re in (500, 2000, 10000):
c, nx, ny = 128, 2048, 1024
st = int(conv_steps(c, 400) * scale)
add(f"D{len([r for r in RUNS if r['group']=='D'])+1:02d}_naca0012_re{re}", "D",
f"NACA 0012, α=8°, Re={re}",
"перестройка следа с ростом Re; проверка устойчивости на профиле",
foil_args(c, 8, re, nx, ny), nx * ny, st,
gif=f"naca0012_re{re}.gif", nx=nx)
for al in (0, 6, 12):
c, nx, ny = 128, 2048, 1024
st = int(conv_steps(c, 400) * scale)
add(f"D{len([r for r in RUNS if r['group']=='D'])+1:02d}_naca4412_a{al:02d}", "D",
f"NACA 4412 (с изгибом), α={al}°, Re=1000",
"ненулевой Cl при α=0 — прямая проверка средней линии профиля",
foil_args(c, al, 1000, nx, ny, naca="4412"), nx * ny, st,
gif=f"naca4412_a{al:02d}.gif", nx=nx)
for c in (48, 96, 192):
nx, ny = 16 * c, 8 * c
st = int(conv_steps(c, 300) * scale)
add(f"D{len([r for r in RUNS if r['group']=='D'])+1:02d}_naca_chord{c}", "D",
f"NACA 0012, α=8°, сходимость по хорде {c}",
"Cl и Cd сходятся по разрешению хорды",
foil_args(c, 8, 1000, nx, ny), nx * ny, st, nx=nx)
# ══════════════════════════════════════════════════════════════════════════════
# E. Сложная и множественная геометрия
# ══════════════════════════════════════════════════════════════════════════════
def group_e(scale):
# тандем: литература даёт переключение режимов около L/D 3.5–4
for ld in (1.5, 3.0, 5.0):
d, nx, ny = 48, 1440, 720
gap = int(ld * d)
bodies = f"cylinder:d={d},x=360,y=360; cylinder:d={d},x={360 + gap},y=360"
st = int(conv_steps(d, 800) * scale)
add(f"E{len([r for r in RUNS if r['group']=='E'])+1:02d}_tandem_ld{int(ld*10)}", "E",
f"Тандем цилиндров, зазор L/D={ld}",
"при малом зазоре у заднего тела ОТРИЦАТЕЛЬНОЕ сопротивление; переход около 3.5–4",
["--bodies", bodies, "--nx", str(nx), "--ny", str(ny), "--re", "150", "--refine", "1"],
nx * ny, st, gif=f"tandem_ld{int(ld*10)}.gif", nx=nx)
# решётка как пористая среда
for step_d, tag in ((3, "sparse"), (2, "dense")):
d, nx, ny = 24, 1440, 720
bodies = "; ".join(
f"cylinder:d={d},x={360 + i * step_d * d},y={360 + (j - 1.5) * step_d * d}"
for i in range(4) for j in range(4))
st = int(conv_steps(d, 600) * scale)
add(f"E{len([r for r in RUNS if r['group']=='E'])+1:02d}_array_{tag}", "E",
f"Решётка 4×4 цилиндров, шаг {step_d}D — пористая среда",
"суммарное сопротивление и структура течения сквозь набор тел",
["--bodies", bodies, "--nx", str(nx), "--ny", str(ny), "--re", "150", "--refine", "1"],
nx * ny, st, gif=f"array_{tag}.gif", nx=nx)
# перфорированная пластина: набор коротких пластин со щелями
for open_frac, tag in ((0.2, "open20"), (0.5, "open50")):
nx, ny = 1200, 600
seg, gap = 40, int(40 * open_frac / (1 - open_frac))
bodies = "; ".join(
f"plate:d={seg},t=0.12,a=90,x=300,y={60 + k * (seg + gap)}"
for k in range(max(1, (ny - 120) // (seg + gap))))
st = int(conv_steps(seg, 500) * scale)
add(f"E{len([r for r in RUNS if r['group']=='E'])+1:02d}_perf_{tag}", "E",
f"Перфорированная пластина, скважность {open_frac}",
"струи в щелях и общее сопротивление; проверка множества тонких тел",
["--bodies", bodies, "--nx", str(nx), "--ny", str(ny), "--re", "500", "--refine", "1"],
nx * ny, st, gif=f"perf_{tag}.gif", nx=nx)
# многоэлементный профиль: основной + предкрылок
nx, ny = 2048, 1024
bodies = ("naca:naca=4412,d=200,a=6,x=560,y=512; "
"naca:naca=2412,d=70,a=18,x=430,y=470")
st = int(conv_steps(200, 300) * scale)
add(f"E{len([r for r in RUNS if r['group']=='E'])+1:02d}_multielem", "E",
"Многоэлементный профиль: основной 4412 плюс предкрылок",
"щель между элементами и раздельные Cl/Cd по телам",
["--bodies", bodies, "--nx", str(nx), "--ny", str(ny), "--re", "2000", "--refine", "1"],
nx * ny, st, gif="multielem.gif", nx=nx)
# сверхтонкая пластина: предел субсеточности, тело тоньше клетки
for t, tag in ((0.02, "t05"), (0.01, "t025")):
nx, ny = 1200, 600
st = int(conv_steps(50, 500) * scale)
add(f"E{len([r for r in RUNS if r['group']=='E'])+1:02d}_thin_{tag}", "E",
f"Сверхтонкая пластина, толщина {t*50:.2f} клетки",
"предел субсеточной границы: тело тоньше клетки; ждём откаты на простой отскок",
["--shape", "plate", "--size", "50", "--thickness", str(t), "--body-angle", "20",
"--nx", str(nx), "--ny", str(ny), "--re", "500", "--refine", "1"],
nx * ny, st, gif=f"thin_{tag}.gif", nx=nx)
# клин и зазубренная кромка — субклеточная деталь на многоугольнике
nx, ny = 1200, 600
st = int(conv_steps(60, 500) * scale)
add(f"E{len([r for r in RUNS if r['group']=='E'])+1:02d}_wedge", "E",
"Клин остриём против потока", "острая кромка: особая точка геометрии",
["--shape", "polygon", "--poly=-30,-18;36,0;-30,18", "--size", "60",
"--nx", str(nx), "--ny", str(ny), "--re", "500", "--refine", "1"],
nx * ny, st, gif="wedge.gif", nx=nx)
saw = ";".join(f"{-30 + 12 * k},{(-1) ** k * 6}" for k in range(6)) + ";30,-18;-30,-18"
add(f"E{len([r for r in RUNS if r['group']=='E'])+1:02d}_serrated", "E",
"Зазубренная задняя кромка", "субклеточные зубцы: как их видит граница",
["--shape", "polygon", f"--poly={saw}", "--size", "60",
"--nx", str(nx), "--ny", str(ny), "--re", "500", "--refine", "1"],
nx * ny, st, gif="serrated.gif", nx=nx)
# ══════════════════════════════════════════════════════════════════════════════
# F. Поле влияния и границы домена
# ══════════════════════════════════════════════════════════════════════════════
def group_f(scale):
d = 48
for up in (2, 6, 16):
nx, ny = int((up + 20) * d), 12 * d
st = int(conv_steps(d, 600) * scale)
add(f"F{len([r for r in RUNS if r['group']=='F'])+1:02d}_inlet{up}d", "F",
f"Отступ до входа {up}D",
"с какого отступа результат перестаёт зависеть от положения входа",
cyl_args(d, 150, nx, ny, body_x=up * d), nx * ny, st)
for down in (5, 15, 40):
nx, ny = int((6 + down) * d), 12 * d
st = int(conv_steps(d, 600) * scale)
add(f"F{len([r for r in RUNS if r['group']=='F'])+1:02d}_outlet{down}d", "F",
f"Отступ до выхода {down}D",
"с какой длины следа выход перестаёт влиять на силы",
cyl_args(d, 150, nx, ny, body_x=6 * d), nx * ny, st)
for r in (1, 2, 3):
d, nx, ny = 32, 640, 384
st = int(conv_steps(d, 600) * scale)
nodes = nx * ny
if r > 1:
ax, bx = int(nx // 5 - 1.5 * d), int(nx // 5 + 6.25 * d)
ay, by = int(ny / 2 - 2 * d), int(ny / 2 + 2 * d)
nodes += r * (r * (bx - ax) + 1) * (r * (by - ay) + 1)
add(f"F{len([r_ for r_ in RUNS if r_['group']=='F'])+1:02d}_refine{r}", "F",
f"Вложенный патч ×{r}",
"вносит ли связка уровней систематику: Cd и St обязаны совпасть",
cyl_args(d, 150, nx, ny, refine=r), nodes, st)
# ══════════════════════════════════════════════════════════════════════════════
# G. Старт, акустика, время жизни
# ══════════════════════════════════════════════════════════════════════════════
def group_g(scale):
d, nx, ny = 48, 960, 576
for re in (150, 2000, 20000):
for init in ("uniform", "rest"):
st = int(conv_steps(d, 300) * scale)
add(f"G{len([r for r in RUNS if r['group']=='G'])+1:02d}_xt_re{re}_{init}", "G",
f"x–t диаграмма, Re={re}, старт {init}",
"откуда идёт стартовое возмущение и раскачивается ли продольная мода",
cyl_args(d, re, nx, ny, init=init, sponge_len=0,
xt=f"xt_re{re}_{init}.csv", xt_every=200),
nx * ny, st, gif=f"xt_re{re}_{init}.gif", nx=nx, extra_out=("xt",))
for sp in (0, 40):
st = int(conv_steps(d, 400) * scale)
add(f"G{len([r for r in RUNS if r['group']=='G'])+1:02d}_sponge{sp}", "G",
f"Губка выхода {sp} столбцов при Re=5000",
"губка — единственное, что реально ест продольную моду",
cyl_args(d, 5000, nx, ny, sponge_len=sp), nx * ny, st)
for sp in (0, 40):
st = int(conv_steps(d, 400) * scale)
add(f"G{len([r for r in RUNS if r['group']=='G'])+1:02d}_spongein{sp}", "G",
f"Губка входа {sp} столбцов при Re=5000",
"помогает ли гасить волну до отражения от входа",
cyl_args(d, 5000, nx, ny, sponge_in=sp), nx * ny, st)
for re in (10000, 50000, 200000):
st = int(conv_steps(d, 200) * scale)
add(f"G{len([r for r in RUNS if r['group']=='G'])+1:02d}_limit_re{re}", "G",
f"Предел устойчивости, Re={re}",
"где именно гибнет счёт; часть этих прогонов обязана развалиться",
cyl_args(d, re, nx, ny), nx * ny, st)
# тест на время жизни: один очень длинный прогон
st = int(5_000_000 * scale)
add(f"G{len([r for r in RUNS if r['group']=='G'])+1:02d}_lifetime", "G",
"Тест на время жизни: 5·10⁶ шагов",
"не уплывает ли ⟨ρ⟩ и не деградирует ли решение на очень длинной дистанции",
cyl_args(d, 150, nx, ny), nx * ny, st, gif="lifetime.gif", nx=nx)
# ══════════════════════════════════════════════════════════════════════════════
# H. Инварианты и точность
# ══════════════════════════════════════════════════════════════════════════════
def group_h(scale):
d, nx, ny = 32, 640, 384
st = int(conv_steps(d, 400) * scale)
add(f"H{len([r for r in RUNS if r['group']=='H'])+1:02d}_sym_cyl", "H",
"Симметрия: цилиндр при α=0 без возмущения",
"Cl и Cm обязаны остаться ≈0 — это внутреннее свойство схемы, не литература",
cyl_args(d, 40, nx, ny, pert_amp=0), nx * ny, st)
add(f"H{len([r for r in RUNS if r['group']=='H'])+1:02d}_sym_foil", "H",
"Симметрия: NACA 0012 при α=0 без возмущения",
"Cl ≈ 0 у симметричного профиля под нулевым углом",
foil_args(64, 0, 500, 1024, 512, pert_amp=0), 1024 * 512,
int(conv_steps(64, 300) * scale))
for al, tag in ((8, "plus"), (-8, "minus")):
add(f"H{len([r for r in RUNS if r['group']=='H'])+1:02d}_mirror_{tag}", "H",
f"Зеркальность: NACA 0012 при α={al}°",
"ряды обязаны зеркалиться: Cl меняет знак, Cd совпадает",
foil_args(64, al, 500, 1024, 512), 1024 * 512,
int(conv_steps(64, 300) * scale))
for ul in (0.02, 0.05, 0.1):
st2 = int(conv_steps(d, 400, ul) * scale)
add(f"H{len([r for r in RUNS if r['group']=='H'])+1:02d}_mach{int(ul*100):02d}", "H",
f"Зависимость от числа Маха, u_lat={ul}",
"результат обязан сходиться при Ma→0; расхождение и есть сжимаемостная ошибка",
cyl_args(d, 150, nx, ny, u_lat=ul), nx * ny, st2)
# прямой замер расхождения f32 и f64 от длины прогона
for bk in ("cpu", "gpu"):
add(f"H{len([r for r in RUNS if r['group']=='H'])+1:02d}_prec_{bk}", "H",
f"Точность: одинаковая постановка на {bk}",
"расхождение f32 против f64 в зависимости от длины прогона",
cyl_args(d, 150, nx, ny), nx * ny, int(conv_steps(d, 800) * scale), backend=bk)
# ══════════════════════════════════════════════════════════════════════════════
# I. Сверхмелкие сетки по всем формам
# ══════════════════════════════════════════════════════════════════════════════
def group_i(scale):
shapes = [
("cylinder", ["--shape", "cylinder"]),
("square", ["--shape", "square"]),
("diamond", ["--shape", "diamond"]),
("ellipse", ["--shape", "ellipse", "--thickness", "0.35"]),
("naca", ["--shape", "naca", "--naca", "0012", "--body-angle", "8"]),
("triangle", ["--shape", "triangle"]),
("plate", ["--shape", "plate", "--thickness", "0.08", "--body-angle", "25"]),
("wedge", ["--shape", "polygon", "--poly=-128,-72;154,0;-128,72"]),
]
d, nx, ny = 256, 4096, 2048
st = int(conv_steps(d, 300) * scale)
for tag, extra in shapes:
add(f"I{len([r for r in RUNS if r['group']=='I'])+1:02d}_fine_{tag}", "I",
f"Сверхмелкая сетка 4096×2048, {tag}, D=256",
"держится ли субсеточная граница на предельном разрешении для этой формы",
extra + ["--size", str(d), "--nx", str(nx), "--ny", str(ny),
"--body-x", "820", "--re", "1000", "--refine", "1"],
nx * ny, st, gif=f"fine_{tag}.gif", nx=nx)
bodies = ("cylinder:d=200,x=800,y=1024; cylinder:d=200,x=1400,y=1024; "
"naca:naca=4412,d=300,a=10,x=2200,y=1024")
add(f"I{len([r for r in RUNS if r['group']=='I'])+1:02d}_fine_multi", "I",
"Сверхмелкая сетка, композиция из трёх тел",
"взаимодействие следов на предельном разрешении",
["--bodies", bodies, "--nx", str(nx), "--ny", str(ny), "--re", "1000", "--refine", "1"],
nx * ny, st, gif="fine_multi.gif", nx=nx)
def main():
ap = argparse.ArgumentParser()
ap.add_argument("--scale", type=float, default=2.0,
help="общий множитель длительности всех прогонов")
ap.add_argument("--out", default=os.path.join(os.path.dirname(__file__), "scenarios.json"))
args = ap.parse_args()
for g in (group_a, group_b, group_c, group_d, group_e, group_f, group_g, group_h, group_i):
g(args.scale)
doc = {
"meta": {
"note": "Стоимость в обновлениях узлов — единственная переносимая между машинами "
"мера. Часы драйвер получает, поделив её на фактические MLUPS, которые "
"замеряет на месте (--calibrate).",
"gif_policy": "вся длительность прогона, реальная скорость, 10 кадр/с; кадр "
"прореживается до ширины не больше 1280 пикселей",
"scale": args.scale,
"assumed_gpu_mlups": 1200,
"assumed_cpu_mlups": 22,
},
"runs": RUNS,
}
with open(args.out, "w", encoding="utf-8") as f:
json.dump(doc, f, ensure_ascii=False, indent=1)
print(f"прогонов: {len(RUNS)} файл: {args.out}\n")
print(f"{'гр':>3}{'прогонов':>10}{'обновлений узлов':>20}{'часов GPU':>12}{'часов CPU':>12}")
tot_gpu = tot_cpu = 0.0
for g in sorted({r["group"] for r in RUNS}):
rs = [r for r in RUNS if r["group"] == g]
cost = sum(r["cost"] for r in rs)
hg = sum(r["cost"] / 1.2e9 / 3600 for r in rs if r["backend"] == "gpu")
hc = sum(r["cost"] / 2.2e7 / 3600 for r in rs if r["backend"] == "cpu")
tot_gpu += hg
tot_cpu += hc
print(f"{g:>3}{len(rs):>10}{cost:>20.3e}{hg:>12.1f}{hc:>12.1f}")
print(f"{'—':>3}{len(RUNS):>10}{sum(r['cost'] for r in RUNS):>20.3e}"
f"{tot_gpu:>12.1f}{tot_cpu:>12.1f}")
print(f"\nвсего ≈ {tot_gpu + tot_cpu:.1f} ч при 1200 MLUPS на GPU и 22 на CPU")
if __name__ == "__main__":
main()
-188
View File
@@ -1,188 +0,0 @@
#!/usr/bin/env python3
"""Сверка GPU-пути с эталонным CPU-путём.
Зачем. На настоящем сервере GPU-путь идёт через драйвер NVIDIA. В WSL такого драйвера
для Linux нет, и единственный доступный Vulkan — dzn (Dozen), транслирующий вызовы в
D3D12. Он честно сообщает о себе `conformanceVersion = 0.0.0.0`, то есть набор тестов
соответствия не проходил. Арифметика f32 в обоих API описана одним и тем же IEEE 754,
но точность трансцендентных функций (exp, log, sqrt — а они в энтропийном равновесии
на каждом узле) стандартами задаётся с допуском в несколько ULP, и допуски эти разные.
Плюс энтропийный стабилизатор сравнивает знаменатель с относительным порогом GREL=1e-8:
достаточно мелкого расхождения, чтобы решение «вырожден или нет» поменялось.
Поэтому пригодность dzn нельзя принимать на веру — её надо мерить. Скрипт гоняет два
коротких эталона на CPU (f64) и на GPU и сверяет числа из summary.json.
Прогоны подобраны намеренно НЕ хаотические: вихрь Тейлора–Грина затухает гладко и имеет
точное решение, обтекание цилиндра при Re=20 стационарно. В таких течениях расхождение
f32 и f64 не нарастает, поэтому любое заметное различие означает проблему драйвера, а не
разрядности. На развитой дорожке (Re >= 100) сверять бессмысленно: там разойдётся любая
пара запусков с разной арифметикой.
Порог взят с запасом от факта: на нативном Vulkan разница Cd между CPU-f64 и GPU-f32
измерена и составила 0.007%.
python3 parity.py # сверка
python3 parity.py --tol 2e-3 # свой порог
python3 parity.py --keep out/ # оставить сводки для разбора
Код возврата 1, если хоть одно поле вышло за порог или прогон развалился.
"""
import argparse
import json
import os
import shutil
import subprocess
import sys
import tempfile
HERE = os.path.dirname(os.path.abspath(__file__))
BIN = os.environ.get("KBC2D_BIN") or os.path.join(
HERE, "..", "target", "release", "kbc2d" + (".exe" if os.name == "nt" else ""))
# Поля сводки: "fields" сверяются с порогом, "report" только печатаются.
#
# * energy_end, enstrophy_end — интегралы поля, ловят расхождение самого счёта;
# * cd — интегральная сила, накапливает ошибку по всей границе тела;
# * gamma_mean и degenerate_frac — поведение энтропийного стабилизатора, самый чуткий
# индикатор расхождений в exp/log;
# * rho_mean — сохранение массы, обязано совпадать почти побитово.
#
# Пороги откалиброваны по замеру на НАТИВНОМ драйвере Vulkan (Intel Iris Xe, Windows):
# там же, где f64 и f32 сравниваются в отсутствие какой-либо трансляции, отклонения
# составили cd 1.7e-5, rho_mean 1.5e-6, gamma_mean 6.5e-4, energy_end 3.4e-3 (вихрь
# доведён до затухания). Пороги ниже взяты с запасом к этим числам, но заметно ниже
# того, что дал бы по-настоящему сломанный драйвер.
#
# Прогоны намеренно НЕ доводятся до глубокого затухания: когда энергия падает до 1e-7,
# f32 упирается в собственный шум, и сверять становится нечего — это свойство
# разрядности, а не драйвера.
CASES = [
{
"name": "Тейлор-Грин 128x128, 500 шагов",
"args": ["--case", "taylor-green", "--nx", "128", "--ny", "128",
"--refine", "1", "--steps", "500", "--case-every", "100"],
# taylor_green_error намеренно НЕ в списке сверяемых, только в отчётных: это
# разность почти равных величин, и её относительное отклонение f32 от f64
# достигает десятков процентов на любом, в том числе нативном, драйвере.
# Замерено на нативном Vulkan: 0.0116 (f64) против 0.0218 (f32). Сверять по
# ней бессмысленно, а видеть её полезно.
"fields": {"energy_end": 1.0, "enstrophy_end": 1.0, "gamma_mean": 1.0},
"report": ["taylor_green_error"],
},
{
"name": "цилиндр Re=20, стационар, 8000 шагов",
"args": ["--shape", "cylinder", "--size", "16", "--nx", "320", "--ny", "192",
"--re", "20", "--refine", "1", "--steps", "8000"],
"fields": {"cd": 1.0, "rho_mean": 0.05, "gamma_mean": 2.0,
"degenerate_frac": 5.0},
"report": ["strouhal"],
},
]
def run(backend, args, summary_path):
"""Один прогон. Возвращает разобранную сводку либо строку с ошибкой."""
cmd = [BIN] + args + ["--backend", backend, "--verbose", "quiet",
"--report-every", "0", "--summary", summary_path]
proc = subprocess.run(cmd, capture_output=True, text=True, errors="replace")
if proc.returncode != 0:
tail = (proc.stderr or proc.stdout or "").strip().splitlines()
return None, "\n".join(tail[-6:]) or f"код возврата {proc.returncode}"
try:
with open(summary_path, encoding="utf-8") as f:
return json.load(f), None
except Exception as e:
return None, f"сводка не читается: {e}"
def deviation(ref, got):
"""Относительное отклонение; для околонулевых величин — абсолютное."""
if abs(ref) > 1e-12:
return abs(got - ref) / abs(ref)
return abs(got - ref)
def main():
ap = argparse.ArgumentParser(description="Сверка GPU-пути с CPU-эталоном")
ap.add_argument("--tol", type=float, default=1e-3,
help="базовый относительный порог (по умолчанию 1e-3)")
ap.add_argument("--keep", metavar="DIR", help="куда положить сводки прогонов")
args = ap.parse_args()
if not os.path.exists(BIN):
print(f"не найден бинарь решателя: {BIN}")
return 1
workdir = args.keep or tempfile.mkdtemp(prefix="kbc2d-parity-")
os.makedirs(workdir, exist_ok=True)
failures = []
speeds = []
for case in CASES:
print(f"\n=== {case['name']} ===")
summaries = {}
for backend in ("cpu", "gpu"):
path = os.path.join(workdir, f"{backend}.json")
data, err = run(backend, case["args"], path)
if err:
print(f" {backend}: ПРОГОН НЕ УДАЛСЯ\n {err}")
failures.append(f"{case['name']}: {backend} не отработал")
break
if data.get("blew_up"):
print(f" {backend}: счёт РАЗВАЛИЛСЯ")
failures.append(f"{case['name']}: {backend} развалился")
break
summaries[backend] = data
if len(summaries) != 2:
continue
cpu, gpu = summaries["cpu"], summaries["gpu"]
m_cpu, m_gpu = cpu.get("mlups", 0.0), gpu.get("mlups", 0.0)
speeds.append((case["name"], m_cpu, m_gpu))
ratio = f", ускорение x{m_gpu / m_cpu:.1f}" if m_cpu > 0 else ""
print(f" скорость: CPU {m_cpu:.1f} MLUPS, GPU {m_gpu:.1f} MLUPS{ratio}")
print(f" {'поле':<22}{'CPU (f64)':>16}{'GPU':>16}{'откл.':>12} порог")
for field in case.get("report", []):
if field in cpu and field in gpu:
ref, got = float(cpu[field]), float(gpu[field])
print(f" {field:<22}{ref:>16.6g}{got:>16.6g}{deviation(ref, got):>11.2e}"
f" (только к сведению)")
for field, mult in case["fields"].items():
if field not in cpu or field not in gpu:
print(f" {field:<22}{'нет в сводке':>16}")
continue
ref, got = float(cpu[field]), float(gpu[field])
dev = deviation(ref, got)
tol = args.tol * mult
mark = "" if dev <= tol else " <-- ВЫШЕ ПОРОГА"
print(f" {field:<22}{ref:>16.6g}{got:>16.6g}{dev:>11.2e} {tol:.1e}{mark}")
if dev > tol:
failures.append(f"{case['name']}: {field} отклонилось на {dev:.2e} "
f"при пороге {tol:.1e}")
print()
if speeds:
print("Скорость по прогонам:")
for name, m_cpu, m_gpu in speeds:
print(f" {name}: CPU {m_cpu:.1f} MLUPS, GPU {m_gpu:.1f} MLUPS")
print()
if failures:
print("СВЕРКА НЕ ПРОЙДЕНА:")
for f in failures:
print(f" * {f}")
print("\nGPU-путь на этом драйвере доверия не заслуживает — кампанию на нём "
"запускать нельзя.")
else:
print("Сверка пройдена: GPU считает то же, что CPU в двойной точности.")
if not args.keep:
shutil.rmtree(workdir, ignore_errors=True)
return 1 if failures else 0
if __name__ == "__main__":
sys.exit(main())
-113
View File
@@ -1,113 +0,0 @@
#!/usr/bin/env python3
# Предполётная проверка кампании: каждый сценарий запускается на два шага.
#
# python preflight.py # проверить все
# python preflight.py --group I # только группу
#
# Ловит опечатки в ключах, несовместимые сочетания и пределы железа ДО того, как кампания
# уйдёт считать на девяносто часов. Не бесплатная привычка, а окупившаяся: именно так
# выяснилось, что вся группа сверхмелких сеток падала на пределе GPU в 65535 рабочих групп
# на измерение, а девять прогонов группы F передавали --body-x дважды.
#
# Физику проверка не трогает вовсе: два шага — это ровно «запустилось и не упало».
import argparse
import json
import os
import subprocess
import sys
HERE = os.path.dirname(os.path.abspath(__file__))
BIN = os.environ.get("KBC2D_BIN") or os.path.join(
HERE, "..", "target", "release", "kbc2d" + (".exe" if os.name == "nt" else ""))
# Всё, что пишет файлы, выкидываем: проверяем запуск, а не вывод. Ключ идёт со значением,
# поэтому следующий за ним аргумент тоже пропускается.
DROP = {"--gif", "--xt", "--case-csv", "--gif-field", "--gif-downsample", "--gif-scale",
"--gif-every", "--gif-fps", "--series-every"}
# Строки, которые ничего не объясняют и только вытесняют полезные.
NOISE = (
"note: run with `RUST_BACKTRACE",
"note: Some details are omitted",
"WARNING: dzn is not a conformant",
"stack backtrace",
"Caused by:",
)
def explain(output):
"""Достать из вывода упавшего прогона то, ради чего его вообще читают.
Раньше бралась последняя строка вывода — а у паники Rust последняя строка это
«note: run with RUST_BACKTRACE=1», то есть ровно ноль сведений о причине.
Полезное лежит либо в собственном сообщении решателя, либо на строке ПОСЛЕ
«panicked at»: там текст самой паники.
"""
out = [l.rstrip() for l in output.splitlines() if l.strip()]
out = [l for l in out if not any(n in l for n in NOISE)]
if not out:
return "(без вывода)"
for l in out:
if l.lstrip().startswith("ОШИБКА"):
return l.strip()
for i, l in enumerate(out):
if "panicked at" in l:
rest = [x.strip() for x in out[i + 1:]]
return " | ".join(rest[:4]) if rest else l.strip()
return out[-1].strip()
def main():
ap = argparse.ArgumentParser()
ap.add_argument("--group", help="только эта группа, например I")
ap.add_argument("--scenarios", default=os.path.join(HERE, "scenarios.json"))
args = ap.parse_args()
runs = json.load(open(args.scenarios, encoding="utf-8"))["runs"]
if args.group:
runs = [r for r in runs if r["group"].upper() == args.group.upper()]
if not runs:
print("под выборку не попал ни один прогон")
return 0
if not os.path.exists(BIN):
print(f"не найден бинарь решателя: {BIN}")
return 1
bad = []
for i, r in enumerate(runs, 1):
argv, skip = [], False
for a in r["args"]:
if skip:
skip = False
continue
if a in DROP:
skip = True
continue
argv.append(a)
if "--steps" in argv:
argv[argv.index("--steps") + 1] = "2"
else:
argv += ["--steps", "2"]
p = subprocess.run([BIN] + argv + ["--verbose", "quiet", "--report-every", "0"],
capture_output=True, cwd=HERE)
ok = p.returncode == 0
print("." if ok else "X", end="", flush=True)
if not ok:
out = (p.stdout + p.stderr).decode("utf-8", "replace")
bad.append((r["id"], explain(out)))
if i % 40 == 0:
print(f" {i}/{len(runs)}", flush=True)
print(f"\n\nпроверено сценариев: {len(runs)}, не запустились: {len(bad)}")
for id_, msg in bad:
print(f" {id_}")
print(f" {msg}")
return 1 if bad else 0
if __name__ == "__main__":
sys.exit(main())
@@ -1,110 +0,0 @@
# Драйвер валидационной кампании для Windows. Делает ровно то же, что run_campaign.sh, и
# читает тот же scenarios.json — целевая машина Linux, а этот вариант нужен для локальной
# отладки обвязки.
#
# .\run_campaign.ps1 -DryRun
# .\run_campaign.ps1 -Calibrate
# .\run_campaign.ps1 -Smoke
# .\run_campaign.ps1 -Resume -Group A,B
# .\run_campaign.ps1 -Only cyl_re150 -BudgetHours 6
[CmdletBinding()]
param(
[switch]$DryRun,
[switch]$Resume,
[switch]$Calibrate,
[switch]$Smoke,
[string[]]$Group,
[string]$Only,
[double]$BudgetHours = 0,
[string]$Bin = "..\target\release\kbc2d.exe",
[string]$Scenarios = "scenarios.json",
[string]$Out = "out",
[double]$GpuMlups = 1200,
[double]$CpuMlups = 22
)
Set-Location $PSScriptRoot
if (-not (Test-Path $Bin)) { throw "не найден бинарь решателя: $Bin (собери cargo build --release)" }
if (-not (Test-Path $Scenarios)) { throw "не найден список сценариев: $Scenarios" }
if ($Calibrate) {
"Калибровка на этой машине (три коротких прогона)…"
# 1920x960 добавлена не для красоты: на 95% стоимости кампании сетки крупнее
# 600 тыс. узлов, и оценивать по мелким - значит занижать пропускную способность.
foreach ($c in @(@(240,120,'gpu'), @(960,480,'gpu'), @(1920,960,'gpu'), @(480,240,'cpu'))) {
$o = & $Bin --nx $c[0] --ny $c[1] --size 16 --refine 1 --steps 3000 `
--report-every 3000 --verbose full --backend $c[2] 2>$null
$m = ($o | Select-String 'MLUPS' | Select-Object -Last 1)
" $($c[0])x$($c[1]) на $($c[2]): $(if ($m) { ($m -replace '.*?([0-9.]+) MLUPS.*','$1') } else { 'не измерено' })"
}
""
"Подставь замеренное: .\run_campaign.ps1 -DryRun -GpuMlups <число> -CpuMlups <число>"
return
}
$runs = (Get-Content $Scenarios -Raw -Encoding UTF8 | ConvertFrom-Json).runs
if ($Group) { $g = $Group | ForEach-Object { $_.ToUpper() }; $runs = $runs | Where-Object { $g -contains $_.group.ToUpper() } }
if ($Only) { $runs = $runs | Where-Object { $_.id -like "*$Only*" } }
if ($Smoke) { $runs = $runs | Sort-Object cost | Select-Object -First 3 }
if (-not $runs) { "под выборку не попал ни один прогон"; return }
$totalCost = ($runs | Measure-Object -Property cost -Sum).Sum
$estH = $totalCost / ($GpuMlups * 1e6) / 3600
"прогонов: $($runs.Count) обновлений узлов: {0:e3} ≈ {1:F1} ч при $GpuMlups MLUPS" -f $totalCost, $estH
"результаты: $Out\<id>\"
""
if ($DryRun) {
"{0,-28} {1,-3} {2,-4} {3,12} {4,8} {5}" -f 'id','гр','бэк','шагов','часов','название'
foreach ($r in $runs) {
$m = if ($r.backend -eq 'cpu') { $CpuMlups } else { $GpuMlups }
"{0,-28} {1,-3} {2,-4} {3,12} {4,8:F2} {5}" -f $r.id, $r.group, $r.backend, $r.steps, ($r.cost/($m*1e6)/3600), $r.title
}
return
}
New-Item -ItemType Directory -Force $Out | Out-Null
$summary = Join-Path $Out 'summary.csv'
$log = Join-Path $Out 'campaign.log'
if (-not (Test-Path $summary)) { 'id,group,backend,status,seconds,steps,cost,title' | Out-File $summary -Encoding utf8 }
$started = Get-Date
foreach ($r in $runs) {
$dir = Join-Path $Out $r.id
if ($Resume -and (Test-Path (Join-Path $dir 'summary.json'))) { "· $($r.id) — уже посчитан, пропускаю"; continue }
if ($BudgetHours -gt 0 -and ((Get-Date) - $started).TotalHours -ge $BudgetHours) {
"бюджет $BudgetHours ч исчерпан, останавливаюсь"; break
}
New-Item -ItemType Directory -Force $dir | Out-Null
# относительные имена файлов из сценария кладём внутрь папки прогона
$argv = @()
for ($i = 0; $i -lt $r.args.Count; $i++) {
$v = $r.args[$i]
if ($i -gt 0 -and @('--gif','--xt','--case-csv') -contains $r.args[$i-1]) { $v = Join-Path $dir $v }
$argv += $v
}
$argv += @('--summary', (Join-Path $dir 'summary.json'),
'--csv', (Join-Path $dir 'series.csv'),
'--verbose','full','--report-every','5000')
"$Bin $($argv -join ' ')" | Out-File (Join-Path $dir 'cmd.txt') -Encoding utf8
"▶ $($r.id) ($($r.group), $($r.backend)) $($r.title)"
$t0 = Get-Date
& $Bin @argv 2>&1 | Out-File (Join-Path $dir 'log.txt') -Encoding utf8
$st = if ($LASTEXITCODE -eq 0) { 'ok' } else { 'fail' }
$dt = [int]((Get-Date) - $t0).TotalSeconds
if (Select-String -Path (Join-Path $dir 'log.txt') -Pattern 'РАЗВАЛИЛСЯ' -Quiet) { $st = 'blewup' }
$rep = Get-Content (Join-Path $dir 'log.txt') -Encoding utf8
$k = ($rep | Select-String 'ИТОГОВЫЙ ОТЧЁТ' | Select-Object -First 1).LineNumber
if ($k) { $rep[($k-1)..($rep.Count-1)] | Out-File (Join-Path $dir 'report.txt') -Encoding utf8 }
'{0},{1},{2},{3},{4},{5},{6},"{7}"' -f $r.id,$r.group,$r.backend,$st,$dt,$r.steps,$r.cost,$r.title |
Out-File $summary -Append -Encoding utf8
"$(Get-Date -Format o) $($r.id) $st ${dt}s" | Out-File $log -Append -Encoding utf8
" → $st за $dt с"
}
""
"готово. сводная таблица: $summary"
-152
View File
@@ -1,152 +0,0 @@
#!/usr/bin/env bash
# Драйвер валидационной кампании. Читает scenarios.json, гоняет прогоны по одному, кладёт
# логи, ряды, сводки и гифки каждого в собственную папку и ведёт общую таблицу.
#
# ./run_campaign.sh --dry-run смета: что и сколько будет считаться
# ./run_campaign.sh --calibrate замерить фактические MLUPS этой машины
# ./run_campaign.sh --smoke три коротких прогона: проверить обвязку
# ./run_campaign.sh --resume считать, пропуская уже готовое
# ./run_campaign.sh --group A,B только выбранные группы
# ./run_campaign.sh --only cyl_re150 по подстроке идентификатора
# ./run_campaign.sh --budget-hours 24 остановиться, когда время выйдет
#
# Прогон, который упал или развалился, помечается в сводке и НЕ останавливает кампанию:
# группы E и G специально ищут предел устойчивости.
set -u
cd "$(dirname "$0")"
BIN="${KBC2D_BIN:-../target/release/kbc2d}"
SCEN="${KBC2D_SCENARIOS:-scenarios.json}"
OUT="${KBC2D_OUT:-out}"
# ВНИМАНИЕ: имя GROUPS занято самим bash (список групп пользователя), присваивание в неё
# молча игнорируется, а "$GROUPS" под root разворачивается в 0 — и выборка съедает всю
# кампанию, не сказав ни слова. Отсюда GRP_SEL.
DRY=0; RESUME=0; CALIB=0; SMOKE=0; GRP_SEL=""; ONLY=""; BUDGET=""
GPU_MLUPS="${KBC2D_GPU_MLUPS:-1200}"
CPU_MLUPS="${KBC2D_CPU_MLUPS:-22}"
while [ $# -gt 0 ]; do
case "$1" in
--dry-run) DRY=1 ;;
--resume) RESUME=1 ;;
--calibrate) CALIB=1 ;;
--smoke) SMOKE=1 ;;
--group) GRP_SEL="$2"; shift ;;
--only) ONLY="$2"; shift ;;
--budget-hours) BUDGET="$2"; shift ;;
-h|--help) sed -n '2,20p' "$0"; exit 0 ;;
*) echo "неизвестный ключ: $1" >&2; exit 2 ;;
esac
shift
done
command -v python3 >/dev/null 2>&1 && PY=python3 || PY=python
[ -x "$BIN" ] || { echo "не найден бинарь решателя: $BIN (собери cargo build --release)" >&2; exit 1; }
[ -f "$SCEN" ] || { echo "не найден список сценариев: $SCEN" >&2; exit 1; }
# ── калибровка: три коротких прогона на месте, чтобы оценки в часах были не гаданием ──
if [ "$CALIB" = 1 ]; then
echo "Калибровка на этой машине (три коротких прогона)…"
# 1920x960 добавлена не для красоты: на 95% стоимости кампании сетки крупнее
# 600 тыс. узлов, и оценивать по мелким — значит занижать пропускную способность.
for spec in "240 120 gpu" "960 480 gpu" "1920 960 gpu" "480 240 cpu"; do
set -- $spec
m=$("$BIN" --nx "$1" --ny "$2" --size 16 --refine 1 --steps 3000 --report-every 3000 \
--verbose full --backend "$3" 2>/dev/null | grep -o '[0-9.]* MLUPS' | tail -1)
echo " $1x$2 на $3: ${m:-не измерено}"
done
echo
echo "Подставь замеренное в переменные окружения и запусти смету:"
echo " KBC2D_GPU_MLUPS=<число> KBC2D_CPU_MLUPS=<число> ./run_campaign.sh --dry-run"
exit 0
fi
mkdir -p "$OUT"
SUMMARY="$OUT/summary.csv"
LOG="$OUT/campaign.log"
[ -f "$SUMMARY" ] || echo "id,group,backend,status,seconds,steps,cost,title" > "$SUMMARY"
# Выборку и порядок считает python: разбирать JSON башем — верный способ ошибиться.
PLAN=$("$PY" - "$SCEN" "$GRP_SEL" "$ONLY" "$SMOKE" <<'PYEOF'
import json, sys
scen, groups, only, smoke = sys.argv[1], sys.argv[2], sys.argv[3], sys.argv[4] == "1"
runs = json.load(open(scen, encoding="utf-8"))["runs"]
if groups:
keep = {g.strip().upper() for g in groups.split(",")}
runs = [r for r in runs if r["group"].upper() in keep]
if only:
runs = [r for r in runs if only in r["id"]]
if smoke:
# три самых дешёвых прогона: обвязку проверяем, а не физику
runs = sorted(runs, key=lambda r: r["cost"])[:3]
for r in runs:
print("\t".join([r["id"], r["group"], r["backend"], str(r["cost"]),
str(r["steps"]), r["title"], "\x1f".join(r["args"])]))
PYEOF
)
[ -n "$PLAN" ] || { echo "под выборку не попал ни один прогон"; exit 0; }
total_cost=0; n=0
while IFS=$'\t' read -r id grp bk cost steps title args; do
n=$((n + 1)); total_cost=$((total_cost + cost))
done <<< "$PLAN"
est_h=$("$PY" -c "
import sys
print('%.1f' % (float(sys.argv[1]) / (float(sys.argv[2]) * 1e6) / 3600))" "$total_cost" "$GPU_MLUPS")
echo "прогонов: $n обновлений узлов: $total_cost ≈ $est_h ч при $GPU_MLUPS MLUPS"
echo "результаты: $OUT/<id>/"
echo
if [ "$DRY" = 1 ]; then
printf '%-28s %-3s %-4s %12s %10s %s\n' id гр бэк шагов часов название
while IFS=$'\t' read -r id grp bk cost steps title args; do
m=$GPU_MLUPS; [ "$bk" = cpu ] && m=$CPU_MLUPS
h=$("$PY" -c "import sys;print('%.2f'%(float(sys.argv[1])/(float(sys.argv[2])*1e6)/3600))" "$cost" "$m")
printf '%-28s %-3s %-4s %12s %10s %s\n' "$id" "$grp" "$bk" "$steps" "$h" "$title"
done <<< "$PLAN"
exit 0
fi
started=$(date +%s)
while IFS=$'\t' read -r id grp bk cost steps title args; do
dir="$OUT/$id"
if [ "$RESUME" = 1 ] && [ -f "$dir/summary.json" ]; then
echo "· $id — уже посчитан, пропускаю"
continue
fi
if [ -n "$BUDGET" ]; then
spent=$(( ($(date +%s) - started) ))
limit=$("$PY" -c "import sys;print(int(float(sys.argv[1])*3600))" "$BUDGET")
[ "$spent" -ge "$limit" ] && { echo "бюджет $BUDGET ч исчерпан, останавливаюсь"; break; }
fi
mkdir -p "$dir"
# аргументы разделены символом 0x1f: в них есть пробелы (список тел) и запятые
IFS=$'\x1f' read -r -a ARGV <<< "$args"
ARGV+=(--summary "$dir/summary.json" --csv "$dir/series.csv" --verbose full --report-every 5000)
# относительные имена файлов из сценария кладём внутрь папки прогона
for i in "${!ARGV[@]}"; do
case "${ARGV[$((i-1))]:-}" in
--gif|--xt|--case-csv) ARGV[$i]="$dir/${ARGV[$i]}" ;;
esac
done
printf '%s\n' "$BIN ${ARGV[*]}" > "$dir/cmd.txt"
echo "▶ $id ($grp, $bk) $title"
t0=$(date +%s)
if "$BIN" "${ARGV[@]}" > "$dir/log.txt" 2>&1; then st=ok; else st=fail; fi
t1=$(date +%s); dt=$((t1 - t0))
grep -q "РАЗВАЛИЛСЯ" "$dir/log.txt" 2>/dev/null && st=blewup
sed -n '/ИТОГОВЫЙ ОТЧЁТ/,$p' "$dir/log.txt" > "$dir/report.txt" 2>/dev/null
printf '%s,%s,%s,%s,%s,%s,%s,"%s"\n' "$id" "$grp" "$bk" "$st" "$dt" "$steps" "$cost" "$title" >> "$SUMMARY"
echo "$(date -Is) $id $st ${dt}s" >> "$LOG"
echo " → $st за ${dt} с"
done <<< "$PLAN"
echo
echo "готово. сводная таблица: $SUMMARY"
File diff suppressed because it is too large Load Diff
@@ -1,97 +0,0 @@
# Развёртывание кампании на сервере с NVIDIA. Этот файл САМОДОСТАТОЧЕН: он тянет готовый образ
# из реестра и исходников репозитория не требует. Скопировать на сервер достаточно его одного.
#
# mkdir -p ~/kbc2d && cd ~/kbc2d
# curl -O <ссылка на этот файл> # либо просто перенести файл руками
#
# docker compose -f docker-compose.server.yml --profile check run --rm vulkan
# docker compose -f docker-compose.server.yml --profile check run --rm preflight
# docker compose -f docker-compose.server.yml --profile check run --rm calibrate
# docker compose -f docker-compose.server.yml --profile check run --rm plan
# docker compose -f docker-compose.server.yml up -d
# docker compose -f docker-compose.server.yml logs -f
#
# Порядок именно такой. Сначала vulkan: если карта не видна, кампания молча уйдёт считать
# ничего — точнее, откажется стартовать на первом же прогоне, но узнать об этом через час
# обиднее, чем через минуту. Потом preflight (каждый сценарий стартует на два шага), потом
# calibrate (замер MLUPS этой машины), и только потом сама кампания.
#
# Результаты складываются в ./out на хосте — около 4.4 ГБ гифок и рядов за полный проход.
name: kbc2d
# ── общая часть всех сервисов ────────────────────────────────────────────────
x-kbc2d: &kbc2d
image: notbigghost/kbc2d:1.2.0
pull_policy: missing
volumes:
- ./out:/work/bench/out
environment:
# ГЛАВНОЕ МЕСТО ВСЕГО ФАЙЛА. NVIDIA Container Toolkit подкладывает внутрь Vulkan-ICD
# (nvidia_icd.json) только если в этом списке есть `graphics`. С одним `compute` wgpu не
# увидит НИ ОДНОГО адаптера, и решатель откажется стартовать с --backend gpu. В образе
# значение уже прописано, здесь оно продублировано явно — чтобы его было видно тому, кто
# читает compose, а не Dockerfile.
NVIDIA_DRIVER_CAPABILITIES: compute,utility,graphics
NVIDIA_VISIBLE_DEVICES: all
# Оценки в часах драйвер считает из этих чисел. Умолчания взяты с другой машины —
# подставьте сюда то, что напечатает профиль calibrate.
KBC2D_GPU_MLUPS: ${KBC2D_GPU_MLUPS:-1200}
KBC2D_CPU_MLUPS: ${KBC2D_CPU_MLUPS:-22}
deploy:
resources:
reservations:
devices:
- driver: nvidia
count: all
capabilities: [gpu]
services:
# ── сама кампания: единственный сервис, который поднимается по `up -d` ──────
campaign:
<<: *kbc2d
container_name: kbc2d-campaign
# --resume пропускает всё, у чего уже есть summary.json: перезапуск продолжает с места,
# а не начинает заново. Именно поэтому перезапуск здесь безопасен и дёшев.
command: ["--resume"]
# on-failure, а НЕ unless-stopped: кампания завершается штатно с кодом 0, и политика
# «перезапускать всегда» после её окончания крутила бы контейнер вхолостую по кругу.
# Падение (OOM, перезагрузка хоста) даёт ненулевой код и будет подхвачено.
restart: on-failure:5
stop_grace_period: 30s
# Девяносто часов вывода — это сотни мегабайт журнала; без ограничения он съест диск.
# Полные логи каждого прогона всё равно лежат в out/<id>/log.txt.
logging:
driver: json-file
options:
max-size: "50m"
max-file: "5"
# ── проверки перед запуском (профиль check, сами не поднимаются) ────────────
vulkan:
<<: *kbc2d
profiles: ["check"]
entrypoint: ["vulkaninfo"]
command: ["--summary"]
preflight:
<<: *kbc2d
profiles: ["check"]
entrypoint: ["python3"]
command: ["preflight.py"]
calibrate:
<<: *kbc2d
profiles: ["check"]
command: ["--calibrate"]
plan:
<<: *kbc2d
profiles: ["check"]
command: ["--dry-run"]
# ── если результаты нужны не под root ────────────────────────────────────────
# Контейнер работает от root, поэтому файлы в ./out окажутся root:root. Чтобы они
# принадлежали вам, добавьте в x-kbc2d строку
# user: "${UID}:${GID}"
# и запускайте как UID=$(id -u) GID=$(id -g) docker compose -f … up -d
# (каталог ./out при этом должен быть создан заранее и принадлежать вам).
@@ -1,114 +0,0 @@
# Запуск решателя и кампании в WSL2 — там, где драйвера NVIDIA для Linux не существует.
#
# ЧЕМ ЭТОТ ФАЙЛ ОТЛИЧАЕТСЯ ОТ docker-compose.server.yml. Тот рассчитан на настоящий
# Linux-сервер: просит у Docker устройство через драйвер `nvidia`, а Vulkan-ICD внутрь
# подкладывает NVIDIA Container Toolkit. В WSL2 этого сделать нельзя в принципе —
# у NVIDIA там нет Linux-библиотеки Vulkan, карта отдаётся по протоколу WDDM через
# /dev/dxg. Отсюда и ошибка `could not select device driver "nvidia"`.
#
# Здесь NVIDIA-runtime не запрашивается ВООБЩЕ. Нужно ровно две вещи:
# * проброс /dev/dxg — сама видеокарта;
# * монтирование /usr/lib/wsl — там Microsoft держит libd3d12.so и libdxcore.so.
# Дальше работает dzn (Dozen) из образа: драйвер Mesa, транслирующий Vulkan в D3D12.
# Ни nvidia-container-toolkit, ни nvidia-ctk, ни правки в daemon.json не требуются.
#
# ВАЖНО, ПРОЧТИТЕ ДО ЗАПУСКА КАМПАНИИ. dzn сообщает о себе conformanceVersion 0.0.0.0,
# то есть набор тестов соответствия Vulkan не проходил. Поэтому порядок такой:
#
# docker compose -f docker-compose.wsl.yml --profile check run --rm vulkan # карта видна?
# docker compose -f docker-compose.wsl.yml --profile check run --rm parity # числа те же?
# docker compose -f docker-compose.wsl.yml --profile check run --rm calibrate # скорость какая?
# docker compose -f docker-compose.wsl.yml --profile check run --rm plan
# docker compose -f docker-compose.wsl.yml up -d
# docker compose -f docker-compose.wsl.yml logs -f
#
# Шаг parity — не формальность. Он гоняет вихрь Тейлора-Грина (есть точное решение) и
# стационарное обтекание цилиндра на CPU в f64 и на GPU, и сверяет числа. Если он не
# прошёл, кампанию запускать нельзя: результаты будет нечем защищать.
#
# Результаты складываются в ./out на хосте.
name: kbc2d
# ── общая часть всех сервисов ────────────────────────────────────────────────
x-kbc2d: &kbc2d
image: notbigghost/kbc2d:1.2.0
# Контекст сборки — каталог с этим файлом. Если образа нет ни локально, ни в реестре,
# достаточно `docker compose -f docker-compose.wsl.yml build`: доступ к Docker Hub
# для запуска не обязателен.
build:
context: .
args:
VERSION: "1.2.0"
pull_policy: missing
devices:
- /dev/dxg:/dev/dxg
volumes:
- ./out:/work/bench/out
# Только на чтение: контейнеру нужны отсюда libd3d12.so, libd3d12core.so и
# libdxcore.so. Каталог наполняет сам WSL из C:\Windows\System32\lxss\lib.
- /usr/lib/wsl:/usr/lib/wsl:ro
environment:
# Без этого dzn не найдёт D3D12 и молча не даст ни одного адаптера.
LD_LIBRARY_PATH: /usr/lib/wsl/lib
# Осознанное согласие считать на драйвере, не прошедшем тесты соответствия Vulkan.
# Без этой переменной wgpu прячет dzn и решатель сообщает, что адаптер не найден.
# Прежде чем запускать кампанию, обязательно пройдите профиль parity.
WGPU_ALLOW_UNDERLYING_NONCOMPLIANT_ADAPTER: "1"
# Оценки в часах драйвер считает из этих чисел. Умолчания взяты с другой машины —
# подставьте сюда то, что напечатает профиль calibrate.
KBC2D_GPU_MLUPS: ${KBC2D_GPU_MLUPS:-1200}
KBC2D_CPU_MLUPS: ${KBC2D_CPU_MLUPS:-22}
services:
# ── сама кампания: единственный сервис, который поднимается по `up -d` ──────
campaign:
<<: *kbc2d
container_name: kbc2d-campaign
# --resume пропускает всё, у чего уже есть summary.json: перезапуск продолжает с
# места, а не начинает заново.
command: ["--resume"]
# on-failure, а НЕ unless-stopped: кампания завершается штатно с кодом 0, и политика
# «перезапускать всегда» после её окончания крутила бы контейнер вхолостую.
restart: on-failure:5
stop_grace_period: 30s
logging:
driver: json-file
options:
max-size: "50m"
max-file: "5"
# ── проверки перед запуском (профиль check, сами не поднимаются) ────────────
vulkan:
<<: *kbc2d
profiles: ["check"]
entrypoint: ["vulkaninfo"]
command: ["--summary"]
parity:
<<: *kbc2d
profiles: ["check"]
entrypoint: ["python3"]
command: ["parity.py"]
preflight:
<<: *kbc2d
profiles: ["check"]
entrypoint: ["python3"]
command: ["preflight.py"]
calibrate:
<<: *kbc2d
profiles: ["check"]
command: ["--calibrate"]
plan:
<<: *kbc2d
profiles: ["check"]
command: ["--dry-run"]
# ── если результаты нужны не под root ────────────────────────────────────────
# Контейнер работает от root, поэтому файлы в ./out окажутся root:root. Чтобы они
# принадлежали вам, добавьте в x-kbc2d строку
# user: "${UID}:${GID}"
# и запускайте как UID=$(id -u) GID=$(id -g) docker compose -f ... up -d
# (каталог ./out при этом должен быть создан заранее и принадлежать вам).
-27
View File
@@ -1,27 +0,0 @@
# Запуск кампании на сервере с NVIDIA одной командой:
# docker compose run --rm kbc2d --calibrate # сначала убедиться, что GPU виден
# docker compose up -d # кампания в фоне
# docker compose logs -f # смотреть ход
#
# Результаты складываются в ./out на хосте.
services:
kbc2d:
build: .
image: kbc2d
container_name: kbc2d
volumes:
- ./out:/work/bench/out
environment:
# Без graphics NVIDIA Container Toolkit не подложит Vulkan-ICD и wgpu не увидит карту.
NVIDIA_DRIVER_CAPABILITIES: compute,utility,graphics
deploy:
resources:
reservations:
devices:
- driver: nvidia
count: all
capabilities: [gpu]
# Кампания идёт десятки часов и переживает перезапуск: --resume пропускает готовое.
command: ["--resume"]
restart: "no"
File diff suppressed because it is too large Load Diff
-632
View File
@@ -1,632 +0,0 @@
//! БЛОК СОЗДАНИЯ ГИФОК: перевод поля в индексированный кадр, палитры, служебная надпись и —
//! главное — СИНХРОНИЗАЦИЯ АНИМАЦИИ С ФИЗИЧЕСКИМ ВРЕМЕНЕМ ПОТОКА.
//!
//! Требование: гифка идёт с той же скоростью, что и настоящий поток, независимо от того, с
//! какой скоростью считает машина. Реальная производительность (шагов в секунду wall-clock)
//! в тайминг не входит вообще — берётся только физическая длительность шага δt из `Units`:
//!
//! δt = u_lat·δx/u_phys с/шаг
//! кадр каждые S шагов ⇒ между кадрами S·δt секунд физического времени
//! задержка кадра = S·δt / playback с (playback = 1 ⇒ реальное время)
//!
//! Ограничение формата: задержка в GIF хранится в СОТЫХ долях секунды (u16), т.е. сетка 10 мс,
//! и почти все декодеры считают задержку 0–1 сс как «как можно быстрее». Поэтому реальное время
//! достижимо не при любом S: при δt = 1.7e-4 с (30 м/с, ячейка 0.1 м, u_lat = 0.05) кадр каждые
//! 10 шагов — это 600 кадр/с, чего GIF не умеет. Отсюда два режима:
//!
//! * `--gif-every N` — шаг задан жёстко; если реальное время не представимо, задержка
//! зажимается и в отчёт печатается фактический коэффициент замедления;
//! * `--gif-every auto` — шаг подбирается так, чтобы реальное время получилось точно при
//! целевой частоте `--gif-fps` (по умолчанию 30). Это режим по
//! умолчанию: синхронность важнее, чем конкретное число шагов на кадр.
use std::fs::File;
use std::io::BufWriter;
use crate::math::{Units, R};
// ─────────────────────────────────────────────────────────────────────────────
// Тайминг
// ─────────────────────────────────────────────────────────────────────────────
/// Как кадры анимации ложатся на физическое время.
#[derive(Clone, Copy, Debug)]
pub struct GifPlan {
/// Кадр каждые столько шагов симуляции.
pub stride: u64,
/// ТОЧНАЯ задержка кадра в сотых долях секунды — вообще говоря дробная.
pub exact_delay_cs: R,
/// Физическое время между соседними кадрами, с.
pub frame_dt: R,
/// Запрошенная скорость воспроизведения (1 = реальное время, 0.1 = 10× замедление).
pub requested_playback: R,
/// Фактическая, в среднем по анимации.
pub actual_playback: R,
/// Средняя частота кадров получившейся гифки.
pub fps: R,
/// Шаг подбирался автоматически.
pub auto: bool,
/// Точная задержка нецелая ⇒ задержки кадров чередуются (см. `DelayDither`).
pub dithered: bool,
/// Точная задержка меньше одной сотой ⇒ формат физически не тянет, время сжато.
pub clamped: bool,
}
impl GifPlan {
/// `stride`: Some(N) — жёстко N шагов на кадр; None — подобрать под `target_fps`.
pub fn new(units: &Units, stride: Option<u64>, target_fps: R, playback: R) -> GifPlan {
let dt = units.dt;
let playback = if playback > 0.0 { playback } else { 1.0 };
let auto = stride.is_none();
let stride = match stride {
Some(s) => s.max(1),
None => {
// хотим кадр раз в 1/target_fps секунды экранного времени, что при скорости
// playback отвечает playback/target_fps секундам физического времени
((playback / target_fps / dt).round() as i64).max(1) as u64
}
};
let frame_dt = stride as R * dt;
let exact_delay_cs = 100.0 * frame_dt / playback;
// 0 сотых большинство декодеров трактует как «как можно быстрее», поэтому пол — единица
let clamped = exact_delay_cs < 1.0;
let mean_delay_cs = if clamped { 1.0 } else { exact_delay_cs };
GifPlan {
stride,
exact_delay_cs,
frame_dt,
requested_playback: playback,
actual_playback: frame_dt / (mean_delay_cs / 100.0),
fps: 100.0 / mean_delay_cs,
auto,
dithered: !clamped && (exact_delay_cs - exact_delay_cs.round()).abs() > 1e-9,
clamped,
}
}
/// Отношение фактической скорости воспроизведения к запрошенной (1.0 = точно).
pub fn sync_error(&self) -> R {
self.actual_playback / self.requested_playback
}
/// Синхронность выдержана: средняя задержка совпадает с точной.
pub fn is_synced(&self) -> bool {
(self.sync_error() - 1.0).abs() < 0.02
}
/// Совет, как починить рассинхрон, если формат его не тянет.
pub fn advice(&self, units: &Units) -> Option<String> {
// разойтись с запрошенной скоростью можно только одним способом: точная задержка
// не влезла в минимальную единицу формата и была зажата
if !self.clamped {
return None;
}
let good = ((0.01 * self.requested_playback / units.dt).ceil() as i64).max(1);
Some(format!(
"кадр требует задержки {:.4} сотых — меньше минимальной единицы формата, поэтому \
анимация идёт в {:.2}× быстрее запрошенного. Починить: --gif-every auto (подберёт \
шаг сам), либо --gif-every {} (минимальный шаг с представимой задержкой), либо \
--gif-speed {:.4} (замедление вместо реального времени)",
self.exact_delay_cs,
self.sync_error(),
good,
self.requested_playback / self.sync_error()
))
}
}
/// Раскладка дробной задержки по целым сотым без накопления ошибки.
///
/// Формат GIF хранит задержку кадра целым числом сотых долей секунды. Точная физическая
/// задержка почти никогда не целая: 30 кадр/с — это 3⅓ сотых. Округление каждого кадра
/// по отдельности даёт систематический уход (3 вместо 3.333 — анимация на 11% быстрее
/// реальности, и расхождение копится линейно). Поэтому задержки выдаются так, чтобы
/// НАКОПЛЕННОЕ время кадров всегда совпадало с накопленным физическим с точностью до одной
/// сотой: для 3.333 получается 3, 3, 4, 3, 3, 4, … Ошибка ограничена ±5 мс и не растёт.
#[derive(Clone, Copy, Debug, Default)]
pub struct DelayDither {
exact_cum: R,
emitted_cum: i64,
}
impl DelayDither {
/// Задержка очередного кадра в сотых долях секунды.
pub fn next(&mut self, exact_delay_cs: R) -> u16 {
self.exact_cum += exact_delay_cs;
let want = self.exact_cum.round() as i64;
let d = (want - self.emitted_cum).max(1);
self.emitted_cum += d;
d.clamp(1, u16::MAX as i64) as u16
}
}
// ─────────────────────────────────────────────────────────────────────────────
// Палитры
// ─────────────────────────────────────────────────────────────────────────────
/// Индексы служебных цветов в палитре (после градиента).
const NCOLORS: usize = 250;
const IDX_BODY: u8 = 250;
const IDX_PATCH: u8 = 251;
const IDX_TEXT: u8 = 252;
const IDX_TEXT_BG: u8 = 253;
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub enum ColorMap {
/// Последовательная, воспринимаемо-равномерная — для |u| и ρ.
Viridis,
/// Яркая последовательная с широким динамическим диапазоном.
Turbo,
/// Расходящаяся сине-бело-красная — для знакопеременных полей (завихренность, γ−2).
CoolWarm,
}
impl ColorMap {
pub fn from_str(s: &str) -> Option<Self> {
Some(match s.to_ascii_lowercase().as_str() {
"viridis" => ColorMap::Viridis,
"turbo" => ColorMap::Turbo,
"coolwarm" | "rdbu" => ColorMap::CoolWarm,
_ => return None,
})
}
pub const ALL: [&'static str; 3] = ["viridis", "turbo", "coolwarm"];
/// Опорные точки градиента; между ними линейная интерполяция в sRGB.
fn anchors(&self) -> &'static [[u8; 3]] {
match self {
ColorMap::Viridis => &[
[68, 1, 84],
[72, 40, 120],
[62, 74, 137],
[49, 104, 142],
[38, 130, 142],
[31, 158, 137],
[53, 183, 121],
[109, 205, 89],
[180, 222, 44],
[253, 231, 37],
],
ColorMap::Turbo => &[
[48, 18, 59],
[70, 107, 227],
[40, 187, 236],
[49, 242, 153],
[138, 252, 61],
[210, 226, 27],
[254, 165, 45],
[239, 89, 17],
[180, 27, 1],
[122, 4, 3],
],
// Расходящаяся палитра Морланда. Опорных точек НЕЧЁТНОЕ число: только тогда
// середина диапазона (для знакопеременного поля — ноль) попадает ровно на
// нейтральный серый, а не между двумя соседними точками.
ColorMap::CoolWarm => &[
[59, 76, 192],
[98, 130, 234],
[141, 176, 254],
[184, 208, 249],
[221, 221, 221],
[245, 196, 173],
[244, 154, 123],
[222, 96, 77],
[180, 4, 38],
],
}
}
/// Палитра GIF: NCOLORS ступеней градиента плюс служебные цвета.
fn palette(&self) -> Vec<u8> {
let a = self.anchors();
let mut p = vec![0u8; 256 * 3];
for k in 0..NCOLORS {
let t = k as R / (NCOLORS - 1) as R * (a.len() - 1) as R;
let i = (t.floor() as usize).min(a.len() - 2);
let fr = t - i as R;
for c in 0..3 {
p[k * 3 + c] = (a[i][c] as R * (1.0 - fr) + a[i + 1][c] as R * fr).round() as u8;
}
}
let set = |p: &mut Vec<u8>, idx: u8, rgb: [u8; 3]| {
let o = idx as usize * 3;
p[o] = rgb[0];
p[o + 1] = rgb[1];
p[o + 2] = rgb[2];
};
set(&mut p, IDX_BODY, [25, 25, 28]);
set(&mut p, IDX_PATCH, [255, 255, 255]);
set(&mut p, IDX_TEXT, [255, 255, 255]);
set(&mut p, IDX_TEXT_BG, [0, 0, 0]);
p
}
}
// ─────────────────────────────────────────────────────────────────────────────
// Микрошрифт 5×7 для служебной надписи
// ─────────────────────────────────────────────────────────────────────────────
/// Битовая маска глифа: 7 строк по 5 бит (старший бит — левый пиксель).
fn glyph(c: char) -> [u8; 7] {
match c.to_ascii_uppercase() {
'0' => [0x0E, 0x11, 0x13, 0x15, 0x19, 0x11, 0x0E],
'1' => [0x04, 0x0C, 0x04, 0x04, 0x04, 0x04, 0x0E],
'2' => [0x0E, 0x11, 0x01, 0x02, 0x04, 0x08, 0x1F],
'3' => [0x1F, 0x02, 0x04, 0x02, 0x01, 0x11, 0x0E],
'4' => [0x02, 0x06, 0x0A, 0x12, 0x1F, 0x02, 0x02],
'5' => [0x1F, 0x10, 0x1E, 0x01, 0x01, 0x11, 0x0E],
'6' => [0x06, 0x08, 0x10, 0x1E, 0x11, 0x11, 0x0E],
'7' => [0x1F, 0x01, 0x02, 0x04, 0x08, 0x08, 0x08],
'8' => [0x0E, 0x11, 0x11, 0x0E, 0x11, 0x11, 0x0E],
'9' => [0x0E, 0x11, 0x11, 0x0F, 0x01, 0x02, 0x0C],
'A' => [0x0E, 0x11, 0x11, 0x1F, 0x11, 0x11, 0x11],
'B' => [0x1E, 0x11, 0x11, 0x1E, 0x11, 0x11, 0x1E],
'C' => [0x0E, 0x11, 0x10, 0x10, 0x10, 0x11, 0x0E],
'D' => [0x1C, 0x12, 0x11, 0x11, 0x11, 0x12, 0x1C],
'E' => [0x1F, 0x10, 0x10, 0x1E, 0x10, 0x10, 0x1F],
'F' => [0x1F, 0x10, 0x10, 0x1E, 0x10, 0x10, 0x10],
'G' => [0x0E, 0x11, 0x10, 0x17, 0x11, 0x11, 0x0F],
'H' => [0x11, 0x11, 0x11, 0x1F, 0x11, 0x11, 0x11],
'I' => [0x0E, 0x04, 0x04, 0x04, 0x04, 0x04, 0x0E],
'K' => [0x11, 0x12, 0x14, 0x18, 0x14, 0x12, 0x11],
'L' => [0x10, 0x10, 0x10, 0x10, 0x10, 0x10, 0x1F],
'M' => [0x11, 0x1B, 0x15, 0x15, 0x11, 0x11, 0x11],
'N' => [0x11, 0x19, 0x15, 0x13, 0x11, 0x11, 0x11],
'O' => [0x0E, 0x11, 0x11, 0x11, 0x11, 0x11, 0x0E],
'P' => [0x1E, 0x11, 0x11, 0x1E, 0x10, 0x10, 0x10],
'R' => [0x1E, 0x11, 0x11, 0x1E, 0x14, 0x12, 0x11],
'S' => [0x0F, 0x10, 0x10, 0x0E, 0x01, 0x01, 0x1E],
'T' => [0x1F, 0x04, 0x04, 0x04, 0x04, 0x04, 0x04],
'U' => [0x11, 0x11, 0x11, 0x11, 0x11, 0x11, 0x0E],
'V' => [0x11, 0x11, 0x11, 0x11, 0x11, 0x0A, 0x04],
'W' => [0x11, 0x11, 0x11, 0x15, 0x15, 0x1B, 0x11],
'X' => [0x11, 0x11, 0x0A, 0x04, 0x0A, 0x11, 0x11],
'Y' => [0x11, 0x11, 0x0A, 0x04, 0x04, 0x04, 0x04],
'Z' => [0x1F, 0x01, 0x02, 0x04, 0x08, 0x10, 0x1F],
'.' => [0x00, 0x00, 0x00, 0x00, 0x00, 0x0C, 0x0C],
',' => [0x00, 0x00, 0x00, 0x00, 0x0C, 0x04, 0x08],
':' => [0x00, 0x0C, 0x0C, 0x00, 0x0C, 0x0C, 0x00],
'-' => [0x00, 0x00, 0x00, 0x1F, 0x00, 0x00, 0x00],
'+' => [0x00, 0x04, 0x04, 0x1F, 0x04, 0x04, 0x00],
'=' => [0x00, 0x00, 0x1F, 0x00, 0x1F, 0x00, 0x00],
'/' => [0x01, 0x02, 0x02, 0x04, 0x08, 0x08, 0x10],
'%' => [0x19, 0x1A, 0x02, 0x04, 0x08, 0x0B, 0x13],
'(' => [0x02, 0x04, 0x08, 0x08, 0x08, 0x04, 0x02],
')' => [0x08, 0x04, 0x02, 0x02, 0x02, 0x04, 0x08],
'Ч' => [0x11, 0x11, 0x11, 0x0F, 0x01, 0x01, 0x01],
_ => [0; 7], // пробел и всё незнакомое
}
}
/// Нарисовать строку в индексированный буфер (координаты — левый верхний угол).
fn draw_text(buf: &mut [u8], w: usize, h: usize, x0: usize, y0: usize, s: &str, scale: usize) {
let gw = 6 * scale; // 5 пикселей глифа + 1 пробел
// подложка, чтобы текст читался на любом фоне
let tw = s.chars().count() * gw;
for yy in y0.saturating_sub(scale)..(y0 + 7 * scale + scale).min(h) {
for xx in x0.saturating_sub(scale)..(x0 + tw + scale).min(w) {
buf[yy * w + xx] = IDX_TEXT_BG;
}
}
for (k, ch) in s.chars().enumerate() {
let g = glyph(ch);
for (row, bits) in g.iter().enumerate() {
for col in 0..5 {
if bits & (1 << (4 - col)) == 0 {
continue;
}
for sy in 0..scale {
for sx in 0..scale {
let px = x0 + k * gw + col * scale + sx;
let py = y0 + row * scale + sy;
if px < w && py < h {
buf[py * w + px] = IDX_TEXT;
}
}
}
}
}
}
}
// ─────────────────────────────────────────────────────────────────────────────
// Диапазон нормировки
// ─────────────────────────────────────────────────────────────────────────────
/// Диапазон значений, отображаемый на палитру. Фиксируется на весь прогон: плавающая
/// автонормировка делает анимацию нечитаемой — цвет перестаёт что-либо значить.
#[derive(Clone, Copy, Debug)]
pub struct Range {
pub lo: R,
pub hi: R,
}
impl Range {
#[inline]
fn index(&self, v: R) -> u8 {
if !v.is_finite() {
return 0;
}
let t = ((v - self.lo) / (self.hi - self.lo)).clamp(0.0, 1.0);
(t * (NCOLORS - 1) as R).round() as u8
}
}
// ─────────────────────────────────────────────────────────────────────────────
// Запись гифки
// ─────────────────────────────────────────────────────────────────────────────
/// Что подписывать на кадре.
pub struct Hud {
pub u_phys: R,
pub field: &'static str,
pub show: bool,
}
pub struct GifWriter {
enc: gif::Encoder<BufWriter<File>>,
pub plan: GifPlan,
/// Размеры ИСХОДНОЙ сетки.
src_nx: usize,
src_ny: usize,
/// Во сколько раз клетки усредняются в пиксель перед отрисовкой.
down: usize,
/// Размеры картинки после прореживания.
nx: usize,
ny: usize,
scale: usize,
w: usize,
h: usize,
range: Range,
hud: Hud,
patch: Option<(usize, usize, usize, usize)>,
dither: DelayDither,
pub frames: u32,
}
impl GifWriter {
#[allow(clippy::too_many_arguments)]
pub fn create(
path: &str,
nx: usize,
ny: usize,
scale: usize,
down: usize,
cmap: ColorMap,
range: Range,
plan: GifPlan,
hud: Hud,
patch: Option<(usize, usize, usize, usize)>,
) -> std::io::Result<GifWriter> {
let scale = scale.max(1);
let down = down.max(1);
let (src_nx, src_ny) = (nx, ny);
// прореживание с округлением вверх: последний блок может быть неполным
let nx = nx.div_ceil(down);
let ny = ny.div_ceil(down);
let w = nx * scale;
let h = ny * scale;
let file = BufWriter::new(File::create(path)?);
let mut enc = gif::Encoder::new(file, w as u16, h as u16, &cmap.palette())
.map_err(std::io::Error::other)?;
enc.set_repeat(gif::Repeat::Infinite).map_err(std::io::Error::other)?;
Ok(GifWriter {
enc,
plan,
src_nx,
src_ny,
down,
nx,
ny,
scale,
w,
h,
range,
hud,
patch,
dither: DelayDither::default(),
frames: 0,
})
}
/// Добавить кадр. `field` — значения по узлам L0, `solid` — маска тела,
/// `t_phys` — физическое время этого кадра (с), `cd` — текущий коэффициент сопротивления.
pub fn push(
&mut self,
field: &[R],
solid: &[bool],
t_phys: R,
cd: Option<R>,
) -> std::io::Result<()> {
// Прореживание: блок down×down усредняется в один пиксель, а телом пиксель считается,
// если тело занимает хотя бы половину блока. Без этого кадр с сетки 4096×2048 весит
// столько, что гифка становится непригодной.
let (field, solid) = if self.down == 1 {
(field.to_vec(), solid.to_vec())
} else {
let d = self.down;
let mut fv = vec![0.0; self.nx * self.ny];
let mut sv = vec![false; self.nx * self.ny];
for gy in 0..self.ny {
for gx in 0..self.nx {
let (mut acc, mut cnt, mut sol) = (0.0, 0usize, 0usize);
for yy in gy * d..((gy + 1) * d).min(self.src_ny) {
for xx in gx * d..((gx + 1) * d).min(self.src_nx) {
let k = yy * self.src_nx + xx;
acc += field[k];
sol += solid[k] as usize;
cnt += 1;
}
}
let g = gy * self.nx + gx;
fv[g] = if cnt > 0 { acc / cnt as R } else { 0.0 };
sv[g] = cnt > 0 && 2 * sol >= cnt;
}
}
(fv, sv)
};
let (field, solid) = (&field[..], &solid[..]);
let mut buf = vec![0u8; self.w * self.h];
// строка 0 изображения — это ВЕРХ, а y = 0 решётки — низ канала: переворачиваем
for y in 0..self.ny {
let src = (self.ny - 1 - y) * self.nx;
for x in 0..self.nx {
let idx = if solid[src + x] { IDX_BODY } else { self.range.index(field[src + x]) };
for sy in 0..self.scale {
let row = (y * self.scale + sy) * self.w + x * self.scale;
for sx in 0..self.scale {
buf[row + sx] = idx;
}
}
}
}
if let Some((ax, bx, ay, by)) = self.patch {
let d = self.down;
self.draw_patch_outline(&mut buf, ax / d, bx / d, ay / d, by / d);
}
if self.hud.show {
let line = match cd {
Some(c) => format!(
"T={:.4}S U={:.1}M/S CD={:.2}",
t_phys, self.hud.u_phys, c
),
None => format!("T={:.4}S U={:.1}M/S", t_phys, self.hud.u_phys),
};
let sc = (self.scale.min(3)).max(1);
draw_text(&mut buf, self.w, self.h, 2 * sc, 2 * sc, &line, sc);
draw_text(&mut buf, self.w, self.h, 2 * sc, self.h - 9 * sc, self.hud.field, sc);
}
let mut frame = gif::Frame::from_indexed_pixels(self.w as u16, self.h as u16, buf, None);
// задержка берётся из накопителя: сумма задержек кадров отслеживает точное
// физическое время, а не округляется покадрово (см. DelayDither)
frame.delay = self.dither.next(self.plan.exact_delay_cs.max(1.0));
self.enc.write_frame(&frame).map_err(std::io::Error::other)?;
self.frames += 1;
Ok(())
}
/// Контур области измельчения — видно, где считает тонкая сетка.
fn draw_patch_outline(&self, buf: &mut [u8], ax: usize, bx: usize, ay: usize, by: usize) {
let flip = |y: usize| self.ny - 1 - y.min(self.ny - 1);
let mut put = |x: usize, y: usize| {
if x < self.nx && y < self.ny {
let py = flip(y) * self.scale;
let px = x * self.scale;
for sy in 0..self.scale {
for sx in 0..self.scale {
let o = (py + sy) * self.w + px + sx;
if o < buf.len() {
buf[o] = IDX_PATCH;
}
}
}
}
};
for x in ax..=bx.min(self.nx - 1) {
put(x, ay);
put(x, by);
}
for y in ay..=by.min(self.ny - 1) {
put(ax, y);
put(bx, y);
}
}
pub fn finish(self) -> std::io::Result<u32> {
Ok(self.frames)
}
}
#[cfg(test)]
mod tests {
use super::*;
/// Пример из постановки: 30 м/с, ячейка 0.1 м, u_lat = 1 ⇒ 300 шагов/с; кадр каждые
/// 10 шагов ⇒ 30 кадр/с. Точная задержка 3⅓ сотых — нецелая, поэтому включается дизеринг,
/// и в среднем гифка идёт ровно в реальном времени.
#[test]
fn spec_example_is_realtime() {
let u = Units::new(0.1, 30.0, 1.0, 150.0, 16.0);
assert!((u.steps_per_second() - 300.0).abs() < 1e-9);
let p = GifPlan::new(&u, Some(10), 30.0, 1.0);
assert!((p.exact_delay_cs - 10.0 / 3.0).abs() < 1e-9, "{}", p.exact_delay_cs);
assert!(p.dithered && !p.clamped);
assert!((p.fps - 30.0).abs() < 1e-9);
assert!(p.is_synced(), "коэффициент {}", p.sync_error());
}
/// Дизеринг обязан вести НАКОПЛЕННОЕ время кадров вплотную к физическому, без ухода.
/// Покадровое округление 3.333 → 3 дало бы к тысячному кадру уход в 3.3 секунды.
#[test]
fn dither_tracks_exact_time_without_drift() {
let exact = 10.0 / 3.0;
let mut d = DelayDither::default();
let mut sum = 0i64;
for k in 1..=1000 {
sum += d.next(exact) as i64;
let want = exact * k as R;
assert!(
(sum as R - want).abs() <= 0.5 + 1e-9,
"кадр {k}: накоплено {sum} сотых против точных {want:.3}"
);
}
// а наивное округление к этому моменту отстало бы на треть сотой на кадр
assert!((sum as R - exact * 1000.0).abs() < 1.0);
assert!((3000 - sum).abs() > 300, "дизеринг обязан отличаться от покадрового округления");
}
/// Тайминг обязан зависеть ТОЛЬКО от физики, не от того, как быстро считает машина:
/// один и тот же δt даёт одну и ту же задержку, вдвое больший шаг — вдвое большую.
#[test]
fn timing_ignores_wall_clock() {
let u = Units::new(0.1, 30.0, 1.0, 150.0, 16.0);
let a = GifPlan::new(&u, Some(10), 30.0, 1.0);
let b = GifPlan::new(&u, Some(10), 30.0, 1.0);
assert!((a.exact_delay_cs - b.exact_delay_cs).abs() < 1e-12);
let c = GifPlan::new(&u, Some(20), 30.0, 1.0);
assert!((c.exact_delay_cs - 2.0 * a.exact_delay_cs).abs() < 1e-12);
assert!(c.is_synced() && a.is_synced());
}
/// Авто-режим обязан вытянуть реальное время при физически корректном u_lat, где жёсткий
/// шаг «каждые 10» потребовал бы 600 кадр/с и формат бы этого не вытянул.
#[test]
fn auto_stride_restores_realtime() {
let u = Units::new(0.1, 30.0, 0.05, 150.0, 16.0);
assert!((u.steps_per_second() - 6000.0).abs() < 1e-6);
let fixed = GifPlan::new(&u, Some(10), 30.0, 1.0);
assert!(fixed.clamped, "0.167 сотых обязаны упереться в пол формата");
assert!(!fixed.is_synced(), "жёсткий шаг не должен был попасть в реальное время");
assert!(fixed.advice(&u).is_some(), "в таком режиме обязан быть совет");
let auto = GifPlan::new(&u, None, 30.0, 1.0);
assert_eq!(auto.stride, 200); // 1/30 с при 6000 шаг/с
assert!(auto.is_synced(), "коэффициент {}", auto.sync_error());
assert!((auto.fps - 30.0).abs() < 1e-9);
}
/// Замедление: playback = 0.1 растягивает то же физическое время ровно в 10 раз.
#[test]
fn slow_motion_scales_delay() {
let u = Units::new(0.1, 30.0, 1.0, 150.0, 16.0);
let fast = GifPlan::new(&u, Some(10), 30.0, 1.0);
let slow = GifPlan::new(&u, Some(10), 30.0, 0.1);
assert!((slow.exact_delay_cs - 10.0 * fast.exact_delay_cs).abs() < 1e-9);
assert!((slow.actual_playback - 0.1).abs() < 1e-9);
assert!(slow.is_synced());
}
/// Палитра обязана быть ровно 256 цветов по 3 байта и не затирать служебные индексы.
#[test]
fn palette_layout() {
for cm in [ColorMap::Viridis, ColorMap::Turbo, ColorMap::CoolWarm] {
let p = cm.palette();
assert_eq!(p.len(), 768);
let o = IDX_BODY as usize * 3;
assert_eq!(&p[o..o + 3], &[25, 25, 28]);
}
}
}
File diff suppressed because it is too large Load Diff
File diff suppressed because it is too large Load Diff
File diff suppressed because it is too large Load Diff