CFDManager/docs/theory/2d_solver

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+.

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, потом обычная.

На сервер удобнее ставить образом — см. Развёртывание.

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

./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. Каждый прогон кладёт логи, ряды, машиночитаемую сводку и гифку на всю свою длительность в собственную папку.

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:

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:):

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.

Три вещи, которые надо знать про этот путь.

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), геометрия при этом собирается из любого числа тел, но силы сверх четвёртого сливаются вместе.