diff --git a/docs/theory/2d_solver/README.md b/docs/theory/2d_solver/README.md index f02454c..948530e 100644 --- a/docs/theory/2d_solver/README.md +++ b/docs/theory/2d_solver/README.md @@ -5,11 +5,12 @@ первоисточников, двумя взаимозаменяемыми бэкендами и анимацией, привязанной к физическому времени потока. -Модель столкновения — **KBC D** по таблице I работы Bösch, Chikatamarla, Karlin, *Entropic -Multi-Relaxation Models for Simulation of Fluid Turbulence* (arXiv:1507.02509); в трёхмерных -работах тех же авторов она называется **KBC-N1**. Сдвиговая часть `s` несёт натуральные -моменты {N, Π_xy}, всё остальное (T, Q_xyy, Q_yxx, A) уходит в `h`. Статьи лежат в -`docs/origins/`; ссылки на формулы в коде даны по нумерации 2D-статьи. +Оператор столкновения — энтропийный 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-статьи. ## Устройство @@ -30,8 +31,8 @@ Multi-Relaxation Models for Simulation of Fluid Turbulence* (arXiv:1507.02509); ```sh cargo build --release # с GPU-бэкендом cargo build --release --no-default-features # только CPU (без wgpu) -cargo test --release # 21 быстрый тест -cargo test --release -- --ignored # плюс эталонный сдвиговый слой (~3 с) +cargo test --release # 22 быстрых теста +cargo test --release -- --include-ignored # плюс эталоны статьи (~18 с) ``` Пример: цилиндр Re=150, гифка завихренности в реальном времени. @@ -89,9 +90,10 @@ cargo test --release -- --ignored # плюс эталонный сд вложенного патча и его границы; - **тело**: семь форм на выбор — `cylinder`, `square`, `diamond`, `ellipse`, `naca`, `triangle`, `plate` — плюс характерный размер, относительная толщина, угол атаки и положение; -- **время**: число шагов, длина разгона, амплитуда и длительность стартового возмущения; -- **схема**: оператор столкновения (`kbc`/`bgk`), режим выхода, поглощающая губка, бэкенд, - число потоков; +- **время**: число шагов, начальное поле (однородный поток либо покой с разгоном), длина + разгона, амплитуда и длительность стартового возмущения; +- **схема**: оператор столкновения (`kbc`/`bgk`), состав сдвиговой части (`n1`/`n2`), режим + выхода, поглощающая губка, бэкенд, число потоков; - **анимация**: файл, поле (`speed`/`vorticity`/`density`/`gamma`), палитра, масштаб, шаг кадра, частота, скорость воспроизведения, диапазон нормировки; - **вывод**: период живых строк, три уровня подробности, число окон в отчёте о сходимости, CSV. @@ -111,116 +113,185 @@ cargo test --release -- --ignored # плюс эталонный сд ## Состояние проверки -**Ядро схемы проверено против формул статей** (`cargo test`, 21 тест + 1 длинный): +Опора — статьи авторов метода из `docs/origins/`, а не сторонние реализации. Ссылки на формулы +даны по нумерации 2D-статьи (arXiv:1507.02509); там, где полезнее формулировка из работы 2024 +года по трёхмерной реализации, это отмечено отдельно. -- проектор Δs, выписанный аналитически из представления популяций через натуральные моменты - (ур. 10), совпадает с матричным `M⁻¹·diag(…,1,1,…)·M` до 1e-13; он идемпотентен и не несёт - ни массы, ни импульса; -- равновесие в product-form сохраняет ρ и ρu до 1e-13; -- γ из замкнутой оценки (ур. 17) — корень условия критической точки энтропии (ур. 15): невязка - при γ\* более чем в 20 раз меньше, чем при γ\*±1; -- при γ = 2 схема совпадает с LBGK поточечно; -- **сдвиговые моменты релаксируют ровно с 2β при любой γ** (проверено при β = 0.3, 0.6, 0.95) — - это и есть гарантия того, что стабилизатор не трогает вязкость; -- затухание сдвиговой волны даёт ν из ур. (5) с погрешностью < 1% при τ = 0.6 и τ = 1.0; -- Zou–He ставит ровно заданные скорость на входе и плотность на выходе; -- SDF всех семи форм: знак верен внутри и снаружи, |∇φ| = 1 ± 0.05 на контрольном кольце. +### Оператор столкновения сверен с листингом статьи -**Сквозная сверка с эталоном.** Дважды периодический сдвиговый слой — один из трёх бенчмарков -2D-статьи (N=128, Re=30000, u₀=0.04, κ=80, δ=0.05, одно конвективное время). Отношение -энстрофии к начальной сходится с fp64-эталоном питоновского решателя **0.6035**. Тест -чувствителен именно к тому, что важно: на испорченном (абсолютном) пороге вырожденности γ тот -же прогон давал 0.6599, то есть +9.3%, а чистый LBGK при этих параметрах разваливается. +Работа 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)**. -**Паритет бэкендов.** CPU (f64) и GPU (f32) на одной постановке совпадают до 4–5 значащих -цифр шаг в шаг: ⟨ρ⟩ 1.04933 против 1.04934, Cd 2.339 против 2.338, ⟨γ⟩ 1.2645 против 1.2646. -На Intel Iris Xe GPU даёт ≈82 MLUPS против ≈18 MLUPS у процессора. +Проверено тестами (`cargo test`, 22 быстрых + 4 длинных): -**Согласованность уровней AMR.** Один и тот же случай, посчитанный с патчем ×2 и вовсе без -измельчения (`--refine 1`), даёт St 0.1951 против 0.1970 и ⟨Cd⟩ 1.956 против 1.943 — расхождение -в пределах 1–3%. То есть связка уровней (рамка, подшаги, рестрикция) не вносит систематики. +- **проектор Δ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. -**Формы тел ведут себя физично.** Один короткий прогон на каждую форму (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) — -то есть и плечо, и знак момента считаются осмысленно. +### Вихрь Тейлора–Грина: второй порядок сходимости (разд. VI) -### Сверка канального случая с питоновским решателем +Единственное из трёх эталонных течений статьи с ТОЧНЫМ аналитическим решением, поэтому +проверяется не «похоже на чужой прогон», а прямое совпадение с формулой: -Постановка ровно та же, что в `solver_2x_sdf` по умолчанию: 174×90, D=16, cx=40, Re=150, -U=0.07, патч ×2 на [16,140]×[13,77], жёсткий ноль поперечной скорости на выходе, 100 000 шагов. -Слева — что даёт этот решатель, справа — что записано в README питоновского. + 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₂²)]. -| величина | kbc2d | solver_2x_sdf | | -|---|---|---|---| -| ⟨ρ⟩ | 1.001 | 1.001 | совпало | -| ⟨Cm⟩ | +0.00001 | −0.00003 | оба ≈ 0 — симметрия считывания в порядке | -| rms Cl (с попр.) | 0.571 | 0.557 | расхождение 2.5% | -| ⟨Cd⟩ (сырой) | 1.956 | 1.799 | **на 9% выше** | -| St·(1−β) | 0.161 | 0.183 | **на 12% ниже** | +Старт — приближением Града (ур. 58), как в статье. Метрика — как на рис. 1: +Σ|u_x − u_x^точн| / Σ|u_x^точн| в момент t_c. -Расхождение по Cd и St не объяснено. Что про него известно: +| 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 | -- оно **не** от связки уровней: без измельчения картина та же (см. выше); -- оно **не** от ядра схемы: сдвиговый слой из статьи воспроизводится с эталоном, вязкость точна - до 1%, сдвиговые моменты релаксируют ровно с 2β; -- после поправки на блокировку **этот** Cd = 1.32 ближе к литературным 1.33 (0.6%), чем - питоновский 1.22 (8.6% ниже), а вот питоновский St ближе к литературным 0.183; -- rms Cl у обоих решателей вдвое выше литературной ~0.3 — это известная особенность постановки, - разобранная в `solver_2x_sdf/README.md` (продольные границы завышают давленческие амплитуды). +**Порядок 2.05 и 2.01** — второй порядок статьи воспроизведён, причём отдельно по амплитуде +(скорость затухания) и по форме. -Развести это можно только прогоном обоих кодов бок о бок на одной машине; питоновская версия -требует GPU и CuPy, которых на машине сборки нет, поэтому сравнение сделано с числами, -записанными в её README. +Два места, где пришлось разобраться, и оба поучительны: -### Известные отличия от питоновского решателя +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), что согласуется с утверждением статьи «все модели работают + практически одинаково». -- **Внутри тела столкновение не считается.** Эти популяции фиктивны: Bouzidi перекрывает всё, - что могло бы прийти из тела в жидкость, так что на физику они не влияют. Питоновская версия - считает их наравне со всеми. Побочный эффект — статистика γ здесь собирается строго по - жидкости. -- **Рестрикция дополнительно пропускает узлы, у которых тонкий узел-источник лежит внутри - тела.** В питоновской версии маска строится только по грубому уровню. Случай краевой (тело на - обоих уровнях одно и то же), но здесь он закрыт явно: из фиктивного узла в жидкий переносить - нечего. -- **Точность.** Процессорный бэкенд работает в f64 (питоновский по умолчанию в f32, - переключается переменной `AMR_FP64`); GPU-бэкенд — в f32, потому что в WGSL нет двойной - точности. +### Дважды периодический сдвиговый слой (разд. 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. Доля вырожденных узлов печатается в отчёте — на исправном -пороге она обязана быть ~0. +`GREL = 1e-8` — **относительный** порог, доля от ⟨Δ|Δ⟩, а не абсолютный. Знаменатель ⟨Δh|Δh⟩ +квадратичен по неравновесию и физически мал (~1e-7…1e-9 в развитом следе), поэтому абсолютный +порог срабатывает на подавляющем большинстве узлов и молча подменяет γ на 2 — то есть гонит +чистый LBGK вместо KBC. Работа 2024 года прямо об этом: γ «далеко не постоянна», её эволюция +тесно связана с состоянием потока, и «любой MRT с γ = const не достигнет той же устойчивости». +Доля вырожденных узлов печатается в отчёте; на исправном пороге она обязана быть ~0. Единственное исключение — самый первый шаг: поле в точности равно равновесию, Δ ≡ 0, и порог -честно срабатывает на всех узлах. На γ это не влияет, потому что она умножается на Δh = 0. +честно срабатывает везде. На результат это не влияет: γ умножается на Δh = 0. -## Акустика канала +### Выбор состава сдвиговой части (табл. I) -Пара «вход по скорости / выход по давлению» — недодемпфированный акустический резонатор: -затухание продольной моды идёт как ν(π/Nx)², то есть на длинном домене её почти ничто не гасит. -Это разобрано в факторном исследовании питоновского решателя (`solver_2x_sdf/README.md`, -эксперимент №1): длинные домены без демпфера дают смещённый режим ⟨ρ⟩ ≈ 1.66 или развал счёта. +Ключ `--kbc-model`: -Отсюда два решения в умолчаниях: +- `n1` (умолчание) — s = {N, Π_xy}, только девиатор. 2D-статья: KBC D; 3D: KBC-N1. +- `n2` — s = {N, Π_xy, T}, девиатор со следом. 2D-статья: KBC C; 3D: KBC-N2. -- `--outlet extrapolate` стоит по умолчанию (в том исследовании — лучший вариант по всем - метрикам: жёсткий ноль поперечной скорости отражает вихри дорожки обратно к телу); -- при `nx ≥ 250` и выключенной губке программа печатает предупреждение с готовым рецептом - (`--sponge-len 32`). +По точности они неразличимы, как и заявляет статья: на одной постановке St 0.1828 у обоих, +⟨Cd⟩ 1.4479 против 1.4468, rms Cl 0.383 против 0.388. -Мгновенное ⟨ρ⟩ при этом всё равно осциллирует вокруг единицы — это сама акустическая мода. -Смотреть надо на **оконные средние** в таблице сходимости: именно они должны стоять на 1.000. +Разница — в объёмной вязкости (ур. 57). У `n2` она фиксирована: ξ = ν. У `n1` она равна +c_s²(1/(γβ) − ½) и, поскольку измеренная ⟨γ⟩ ≈ 1.73 < 2, в среднем оказывается примерно +вчетверо БОЛЬШЕ ν. Отрицательной она бывает лишь в долях процента узлов. Практический вывод +против ожидания: `n1` демпфирует продольную акустику сильнее, и в специально испорченной +постановке (старт из покоя, губка выключена) `n1` доживает до конца с пульсацией 75% от U, +а `n2` разваливается. Поэтому умолчание — `n1`. + +### Паритет бэкендов и согласованность уровней + +CPU (f64) и GPU (f32) на одной постановке совпадают до 4–5 значащих цифр шаг в шаг: ⟨ρ⟩ +1.04933 против 1.04934, Cd 2.339 против 2.338, ⟨γ⟩ 1.2645 против 1.2646. На Intel Iris Xe GPU +даёт ≈105 MLUPS против ≈18 MLUPS у процессора. + +Один и тот же случай с патчем ×2 и вовсе без измельчения (`--refine 1`) даёт St 0.1951 против +0.1970 и ⟨Cd⟩ 1.956 против 1.943 — связка уровней систематики не вносит. + +### Обтекание цилиндра против литературы + +Постановка 480×240, D = 24, Re = 150, блокировка β = D/Ny = 0.1, умолчания решателя: + +| величина | сырое | с поправкой на блокировку | литература (безгранич. цилиндр) | +|---|---|---|---| +| St | 0.1828 | 0.1645 | 0.183 | +| ⟨Cd⟩ | 1.448 | 1.173 | 1.33 | +| rms Cl | 0.383 | 0.310 | ~0.30 | +| ⟨Cm⟩ | 0.00002 | — | 0 (симметрия) | + +Поправка: St×(1−β), Cd и rms Cl ×(1−β)². Сырое St и скорректированный rms Cl ложатся на +литературу; Cd после поправки ниже на 12%. ⟨Cm⟩ ≈ 0 — контроль симметрии считывания силы. + +**Формы тел ведут себя физично.** Прогон на каждую форму (260×130, размер 20, угол атаки 12°) +даёт ожидаемый порядок сопротивления: профиль 0.44, эллипс 0.45, пластина 0.60, цилиндр 1.33, +ромб 1.65, квадрат 2.38, треугольник 2.72. ⟨Cm⟩ ≈ 0 **только** у круга (−0.0002), которому +угол атаки безразличен, а у несимметричных под углом тел он ненулевой (профиль +0.31, +эллипс +0.15, пластина +0.13). + +### Мелкие отличия от питоновского прототипа + +Решатель писался заново, не как порт, но пара мест разошлась с `solver_2x_sdf` намеренно: +внутри тела здесь не считается столкновение (эти популяции фиктивны — Bouzidi перекрывает всё, +что могло бы прийти из тела в жидкость; побочно статистика γ собирается строго по жидкости), а +рестрикция дополнительно пропускает узлы, у которых тонкий узел-источник лежит внутри тела. +Процессорный бэкенд работает в f64, GPU-бэкенд — в f32. + +## Продольная акустика канала: почему поток может «дышать» + +Самая заметная ловушка этой постановки, и её стоит понимать до первого запуска. + +**Граничное условие с заданной скоростью на входе акустически есть жёсткий поршень**: оно +отражает продольные волны с коэффициентом +1. Выход по давлению — наоборот, открытый конец. +Вместе они делают из канала четвертьволновый резонатор с пучностью давления на входе и узлом +на выходе: + + период основной моды = 4·Nx/c_s шагов + затухание вязкостью ~ ν(π/2Nx)² — на длинном домене практически ноль + +Измерено на 480×240: период пульсации ⟨ρ⟩ **3332 шага** против расчётных 4·Nx/c_s = 3325 +(0.2%), первый ноль автокорреляции на 827 шагах = ровно Nx/c_s (четверть периода). За 30000 +шагов амплитуда упала на 3.8% — то есть мода не гаснет вообще. Ничто её не подкачивает; это +звон от старта, запертый в почти без потерь резонаторе. + +Отсюда умолчания: + +- **`--init uniform`** — домен сразу заполнен набегающим потоком, вход включён на полную. + Старт из покоя (`--init rest`) разгоняет весь столб жидкости и закачивает моду; при разгоне + за 1000 шагов против акустического пробега 831 шаг это для звука удар. +- **губка перед выходом включена и подобрана по домену** (`nx/12`, не меньше 16 столбцов, + с запасом до патча). Отключается `--sponge-len 0`. +- **`--outlet extrapolate`** — нуль-градиент поперечной скорости; жёсткий ноль отражает вихри + дорожки обратно к телу. + +Эффект умолчаний на той же постановке: + +| | размах ⟨ρ⟩ | пульсация u′/U | +|---|---|---| +| старт из покоя, без губки | 0.1307 | **75%** | +| умолчания | 0.0017 | **1.0%** | + +Заодно выправились и сами числа: St 0.183 вместо 0.194, rms Cl 0.38 вместо 0.84, ⟨ρ⟩ 1.0000. + +Отчёт печатает период моды, время её вязкого затухания и измеренную пульсацию в конце прогона, +с явным предупреждением, если она превысила 5% от U. ## Дальше -- σ·n-кросс-чек силы (интеграл тензора напряжений по контуру) как независимая проверка GMEM — - в питоновской версии есть, здесь пока нет; +- σ·n-кросс-чек силы (интеграл тензора напряжений по контуру) как независимая проверка GMEM; +- согласованные начальные условия по образцу работы 2024 года: там ρ и старшие моменты + получают, решая ∂ρ/∂t + ∇·(ρu₀) = D∇²ρ до стационара, что убрало бы и остаточный + стартовый импульс от появления тела в потоке; - подвижные и вращающиеся тела: GMEM уже записан в галилей-инвариантной форме и принимает скорость стенки на линке, но подача этой скорости не подключена; - несколько тел одновременно (сейчас — одно тело на выбор из семи форм). diff --git a/docs/theory/2d_solver/src/cpu.rs b/docs/theory/2d_solver/src/cpu.rs index c213d8a..fc399e5 100644 --- a/docs/theory/2d_solver/src/cpu.rs +++ b/docs/theory/2d_solver/src/cpu.rs @@ -113,8 +113,11 @@ pub struct Level { } impl Level { - fn new(nx: usize, ny: usize, beta: Vec, geom: Geom) -> Self { - let f0 = math::feq(1.0, 0.0, 0.0); + /// `u0` — скорость, которой заполняется поле на старте. Заполнять сразу набегающим + /// потоком принципиально: старт из покоя разгоняет весь столб жидкости и закачивает в + /// канал продольную акустическую моду, которую вязкость потом почти не гасит. + fn new(nx: usize, ny: usize, beta: Vec, geom: Geom, u0: (R, R)) -> Self { + let f0 = math::feq(1.0, u0.0, u0.1); let n = nx * ny; Level { nx, ny, f: vec![f0; n], post: vec![f0; n], gamma: vec![2.0; n], beta, geom } } @@ -122,7 +125,7 @@ impl Level { /// Столкновение по всем жидким узлам. Внутри тела не считаем: эти популяции фиктивны /// (ГУ Bouzidi перезаписывает всё, что могло бы прийти из тела в жидкость), а счёт там /// только жжёт такты и способен родить NaN при экстремальных режимах. - fn collide(&mut self, op: Collision) -> Stats { + fn collide(&mut self, op: Collision, model: math::KbcModel) -> Stats { let nx = self.nx; let beta = &self.beta; let solid = &self.geom.solid; @@ -141,11 +144,11 @@ impl Level { } let b = beta[node % nx]; let k: Kbc = match op { - Collision::Kbc => math::collide_node(p, b), + Collision::Kbc => math::collide_node(p, b, model), Collision::Bgk => math::collide_node_bgk(p, b), }; *g = k.gamma; - Stats::of(k, b) + Stats::of(k, b, model) }) .reduce(|| Stats::EMPTY, Stats::merge) } @@ -288,14 +291,14 @@ impl Stats { xi_negative: 0, }; #[inline] - fn of(k: Kbc, beta: R) -> Stats { + fn of(k: Kbc, beta: R, model: math::KbcModel) -> Stats { Stats { n: 1, gsum: k.gamma, gmin: k.gamma, gmax: k.gamma, degenerate: k.degenerate as u64, - xi_negative: (math::xi_of_gamma(k.gamma, beta) < 0.0) as u64, + xi_negative: (math::xi_of_gamma(k.gamma, beta, model) < 0.0) as u64, } } fn merge(a: Stats, b: Stats) -> Stats { @@ -545,7 +548,14 @@ impl Sim { }) .collect(); - let l0 = Level::new(nx, ny, beta0, geom0); + // стартовое поле: либо сразу набегающий поток, либо покой + let u0 = if spec.init_uniform { + let (sn, cs) = spec.flow_angle.sin_cos(); + (spec.units.u_lat * cs, spec.units.u_lat * sn) + } else { + (0.0, 0.0) + }; + let l0 = Level::new(nx, ny, beta0, geom0, u0); let fluid_count = l0.geom.solid.iter().filter(|s| !**s).count() as R; let (l1, patch) = if spec.refine > 1 { @@ -561,7 +571,7 @@ impl Sim { let beta1 = 1.0 / (2.0 * tau1); let r01 = tau1 / (r as R * tau0); let patch = Patch::new(&spec, &l0.geom, &geom1.solid, r01); - (Some(Level::new(nfx, nfy, vec![beta1; nfx], geom1)), Some(patch)) + (Some(Level::new(nfx, nfy, vec![beta1; nfx], geom1, u0)), Some(patch)) } else { (None, None) }; @@ -594,12 +604,19 @@ impl Sim { /// Скорость на входе в момент t: разгон smoothstep плюс окно поперечного возмущения. /// - /// Разгон гасит импульсный старт (мгновенное включение входа шлёт по домену ударную волну). - /// Возмущение — короткий поперечный импульс, сбивающий симметрию: без него дорожка Кармана - /// заводится только на численном шуме и стартует на порядок позже. + /// При старте из однородного потока разгон не нужен и пропускается: поле и вход и так + /// согласованы, а плавный разгон поверх согласованного поля сам стал бы рассогласованием. + /// При старте из покоя разгон обязателен — мгновенное включение входа шлёт по домену + /// ударную волну. Возмущение — короткий поперечный импульс, сбивающий симметрию: без него + /// дорожка Кармана заводится только на численном шуме и стартует на порядок позже. fn inlet(&self, t: u64) -> (R, R) { let sp = &self.spec; - let u = sp.units.u_lat * math::smoothstep(t as R / sp.ramp.max(1) as R); + let ramp = if sp.init_uniform { + 1.0 + } else { + math::smoothstep(t as R / sp.ramp.max(1) as R) + }; + let u = sp.units.u_lat * ramp; let (s, c) = sp.flow_angle.sin_cos(); let mut uy = u * s; if sp.pert_dur > 0 && t >= sp.ramp && t < sp.ramp + sp.pert_dur { @@ -617,7 +634,8 @@ impl Sim { // ── уровень 0 ── self.pre.copy_from_slice(&self.l0.f); - let stats0 = self.l0.collide(sp_collision); + let model = self.spec.kbc_model; + let stats0 = self.l0.collide(sp_collision, model); self.l0.stream(); self.l0.apply_bouzidi(); self.l0.free_slip_walls(); @@ -629,7 +647,7 @@ impl Sim { p.ghost_values(&self.pre, &mut self.gh_old); p.ghost_values(&self.l0.f, &mut self.gh_new); for s in 0..p.r { - l1.collide(sp_collision); + l1.collide(sp_collision, model); l1.stream(); l1.apply_bouzidi(); // силу снимаем на КАЖДОМ подшаге и усредняем — мгновенное значение на @@ -767,15 +785,30 @@ mod tests { body_cx: 0.0, body_cy: 0.0, }, + (0.0, 0.0), ) } /// Голый периодический шаг: столкновение + перенос, без единого ГУ. /// Нужен, чтобы отделить ядро схемы от граничных условий. fn periodic_step(f: &mut [[R; Q]], tmp: &mut [[R; Q]], nx: usize, ny: usize, beta: R) { + periodic_step_op(f, tmp, nx, ny, beta, Collision::Kbc) + } + + fn periodic_step_op( + f: &mut [[R; Q]], + tmp: &mut [[R; Q]], + nx: usize, + ny: usize, + beta: R, + op: Collision, + ) { for (n, c) in f.iter().enumerate() { tmp[n] = *c; - math::collide_node(&mut tmp[n], beta); + match op { + Collision::Kbc => math::collide_node(&mut tmp[n], beta, math::KbcModel::N1), + Collision::Bgk => math::collide_node_bgk(&mut tmp[n], beta), + }; } for y in 0..ny { for x in 0..nx { @@ -865,6 +898,183 @@ mod tests { } } + /// Диагностика: из чего складывается ошибка Тейлора–Грина. Сетка и вязкость ЗАФИКСИРОВАНЫ, + /// меняется только амплитуда u₀. Дискретизационная ошибка от u₀ почти не зависит (задача + /// в этом пределе линейна), а сжимаемостная идёт как Ma² ∝ u₀². Наклон и покажет, что + /// доминирует. + #[test] + #[ignore = "диагностика; запуск: cargo test --release -- --ignored taylor_green_error_scaling --nocapture"] + fn taylor_green_error_scaling() { + let n = 128usize; + let nu = 0.0384; + let beta = math::beta_of_nu(nu); + println!(" N={n}, ν={nu} (τ={:.4}) — меняем только u₀", 1.0 / (2.0 * beta)); + let mut prev: Option<(R, R)> = None; + for &u0 in &[0.04, 0.02, 0.01, 0.005] { + let e = taylor_green_error(n, u0, nu, beta); + let slope = prev.map(|(p_u, p_e): (R, R)| (p_e / e.shape).ln() / (p_u / u0).ln()); + match slope { + Some(s) => println!( + " u₀={u0:<7} Ma={:.4} форма {:.3e} ампл {:.3e} наклон формы по u₀: {s:.2}", + u0 / math::CS2.sqrt(), e.shape, e.amplitude), + None => println!(" u₀={u0:<7} Ma={:.4} форма {:.3e} ампл {:.3e}", + u0 / math::CS2.sqrt(), e.shape, e.amplitude), + } + prev = Some((u0, e.shape)); + } + println!(" наклон ≈2 ⇒ правит сжимаемость (Ma²); ≈0 ⇒ правит дискретизация"); + } + + /// Один прогон вихря Тейлора–Грина до полураспада; возвращает Σ|u_x−точн|/Σ|точн|. + /// + /// Начальное состояние ставится ПОЛНОСТЬЮ согласованным: скорость, давление и неравновесная + /// часть. Давление здесь не константа — течение несёт собственное поле порядка ρu₀², + /// которое находится из ∇²p = 2ρ(ψ_xx·ψ_yy − ψ_xy²): + /// + /// p = −(ρu₀²/4)[cos(2k₁x) + (k₁²/k₂²)·cos(2k₂y)], δρ = p/c_s². + /// + /// Если стартовать с ρ ≡ 1, эта разница уходит в акустику, которая в периодическом ящике + /// почти не затухает и садится полкой на ошибку скорости, ломая порядок сходимости. + fn taylor_green_error(n: usize, u0: R, nu: R, beta: R) -> TgError { + taylor_green_error_op(n, u0, nu, beta, Collision::Kbc) + } + + fn taylor_green_error_op(n: usize, u0: R, nu: R, beta: R, op: Collision) -> TgError { + let (k1, k2) = (1.0, 4.0); + let tau = 1.0 / (2.0 * beta); + let kk1 = 2.0 * std::f64::consts::PI * k1 / n as R; + let kk2 = 2.0 * std::f64::consts::PI * k2 / n as R; + let decay = nu * (kk1 * kk1 + kk2 * kk2); + let tc = ((2.0_f64).ln() / decay).round() as usize; + let exact = |x: usize, y: usize, e: R| -> (R, R) { + let (a, b) = (kk1 * x as R, kk2 * y as R); + (-u0 * a.cos() * b.sin() * e, (k1 / k2) * u0 * a.sin() * b.cos() * e) + }; + let mut f = vec![[0.0; Q]; n * n]; + for y in 0..n { + for x in 0..n { + let (ux, uy) = exact(x, y, 1.0); + let (a, b) = (kk1 * x as R, kk2 * y as R); + let p = -(u0 * u0 / 4.0) + * ((2.0 * a).cos() + (k1 * k1) / (k2 * k2) * (2.0 * b).cos()); + let rho = 1.0 + p / math::CS2; + let dxux = u0 * kk1 * a.sin() * b.sin(); + let dyux = -u0 * kk2 * a.cos() * b.cos(); + let dxuy = (k1 / k2) * u0 * kk1 * a.cos() * b.cos(); + let dyuy = -(k1 / k2) * u0 * kk2 * a.sin() * b.sin(); + // Π = ρc_s²δ + ρuu + Π⁽¹⁾, Π⁽¹⁾ = −τρc_s²(∂_αu_β + ∂_βu_α) — ур. (53) + let pref = -tau * rho * math::CS2; + let pxx = rho * math::CS2 + rho * ux * ux + pref * 2.0 * dxux; + let pyy = rho * math::CS2 + rho * uy * uy + pref * 2.0 * dyuy; + let pxy = rho * ux * uy + pref * (dyux + dxuy); + f[y * n + x] = math::grad_init(rho, ux, uy, pxx, pxy, pyy); + } + } + let mut tmp = f.clone(); + for _ in 0..tc { + periodic_step_op(&mut f, &mut tmp, n, n, beta, op); + } + let e = (-decay * tc as R).exp(); + let (mut num, mut den) = (0.0, 0.0); + // заодно раскладываем ошибку: наилучшая подгонка амплитуды к точной форме + let (mut dot, mut nrm2) = (0.0, 0.0); + for y in 0..n { + for x in 0..n { + let got = math::macros(&f[y * n + x]).1; + let want = exact(x, y, e).0; + num += (got - want).abs(); + den += want.abs(); + dot += got * want; + nrm2 += want * want; + } + } + let amp = dot / nrm2; // 1.0 = амплитуда совпала + let (mut snum, mut sden) = (0.0, 0.0); + for y in 0..n { + for x in 0..n { + let got = math::macros(&f[y * n + x]).1; + let want = amp * exact(x, y, e).0; + snum += (got - want).abs(); + sden += want.abs(); + } + } + TgError { total: num / den, amplitude: (amp - 1.0).abs(), shape: snum / sden } + } + + /// Разложение ошибки эталона: полная, вклад амплитуды (скорость затухания) и вклад формы. + struct TgError { + total: R, + amplitude: R, + shape: R, + } + + /// ВИХРЬ ТЕЙЛОРА–ГРИНА — первый эталон 2D-статьи (разд. VI). Единственное из трёх течений + /// статьи, у которого есть ТОЧНОЕ аналитическое решение, поэтому проверяется не «похоже на + /// чужой прогон», а прямое совпадение с формулой и заявленный статьёй ВТОРОЙ ПОРЯДОК + /// сходимости (рис. 1a: Re = 100, u₀ = 0.03, N ∈ {64, 128, 256}). + /// + /// Постановка ровно по статье: + /// 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₂²)]. + /// Волновые числа переводятся в решёточные: K = 2πk/N (узел — единица длины). + /// Старт — приближением Града (ур. 58), как в статье. + /// Метрика — та же, что на рис. 1: Σ|u_x − u_x^точн| / Σ|u_x^точн| в момент t_c. + #[test] + #[ignore = "долгий (до 256², t_c растёт как N²); запуск: cargo test --release -- --ignored"] + fn taylor_green_converges_at_second_order() { + // ДИФФУЗИОННОЕ ИЗМЕЛЬЧЕНИЕ: ν фиксирована, u₀ ∝ 1/N. Тогда Re = u₀N/ν сохраняется + // (все три сетки считают ОДНО И ТО ЖЕ течение), а число Маха падает как 1/N — вместе + // с ним падает и сжимаемостная ошибка метода. Только при таком измельчении второй + // порядок виден целиком; при фиксированном u₀ ошибка упирается в полку O(Ma²), которая + // от сетки не зависит вовсе (см. taylor_green_error_scaling). + let nu = 0.0192; + let beta = math::beta_of_nu(nu); + let mut errs = Vec::new(); + for &n in &[64usize, 128, 256] { + let u0 = 0.03 * 64.0 / n as R; + let e = taylor_green_error(n, u0, nu, beta); + println!( + " N={n:>4} u₀={u0:.5} Re={:.0} полная {:.3e} амплитуда {:.3e} форма {:.3e}", + u0 * n as R / nu, e.total, e.amplitude, e.shape + ); + errs.push((n, e)); + } + for w in errs.windows(2) { + let (n0, e0) = (&w[0].0, &w[0].1); + let (n1, e1) = (&w[1].0, &w[1].1); + let ord = |a: R, b: R| (a / b).ln() / (*n1 as R / *n0 as R).ln(); + let (pt, pa, ps) = ( + ord(e0.total, e1.total), + ord(e0.amplitude, e1.amplitude), + ord(e0.shape, e1.shape), + ); + println!(" порядок {n0}→{n1}: полная {pt:.2} амплитуда {pa:.2} форма {ps:.2}"); + assert!( + (1.7..2.6).contains(&pt), + "порядок {n0}→{n1} = {pt:.2}, статья (разд. VI) заявляет второй" + ); + } + } + + /// Диагностика: тот же вихрь Тейлора–Грина оператором LBGK. Статья (рис. 1) утверждает, + /// что на этом течении LBGK и все варианты KBC идут практически одинаково; если полка по + /// форме есть у обоих — она от схемы D2Q9, а не от энтропийного стабилизатора. + #[test] + #[ignore = "диагностика; запуск: cargo test --release -- --ignored taylor_green_kbc_vs_bgk --nocapture"] + fn taylor_green_kbc_vs_bgk() { + let (u0, re) = (0.03, 100.0); + for &n in &[64usize, 128, 256] { + let nu = u0 * n as R / re; + let beta = math::beta_of_nu(nu); + let k = taylor_green_error_op(n, u0, nu, beta, Collision::Kbc); + let b = taylor_green_error_op(n, u0, nu, beta, Collision::Bgk); + println!( + " N={n:>4} KBC: форма {:.3e} ампл {:.3e} LBGK: форма {:.3e} ампл {:.3e}", + k.shape, k.amplitude, b.shape, b.amplitude + ); + } + } + /// СКВОЗНАЯ СВЕРКА С ЭТАЛОНОМ. Дважды периодический сдвиговый слой — один из трёх /// бенчмарков 2D-статьи (разд. VIII). Постановка: N=128, Re=30000, u0=0.04, κ=80, δ=0.05, /// одно конвективное время t_c = N/u0 шагов. diff --git a/docs/theory/2d_solver/src/gpu.rs b/docs/theory/2d_solver/src/gpu.rs index e0d037f..2d6fc6b 100644 --- a/docs/theory/2d_solver/src/gpu.rs +++ b/docs/theory/2d_solver/src/gpu.rs @@ -49,7 +49,8 @@ struct Dyn { ux_in: f32, uy_in: f32, rho_out: f32, - _pad0: f32, + /// 0 — сдвиговая часть только девиатор (N1), 1 — девиатор со следом (N2) + kbc_model: u32, outlet_extrap: u32, nparts: u32, refine: u32, @@ -126,7 +127,7 @@ struct Results { const SHADER: &str = r#" // ───── структуры (обязаны совпадать с gpu.rs) ───── struct LevelParams { n:u32, nx:u32, ny:u32, nlinks:u32, bcx:f32, bcy:f32, probe_node:u32, flags:u32 }; -struct Dyn { ux_in:f32, uy_in:f32, rho_out:f32, p0:f32, outlet_extrap:u32, nparts:u32, refine:u32, collision:u32 }; +struct Dyn { ux_in:f32, uy_in:f32, rho_out:f32, kbc_model:u32, outlet_extrap:u32, nparts:u32, refine:u32, collision:u32 }; struct Substep { idx:u32, p0:u32, p1:u32, p2:u32 }; struct AmrParams { nghost:u32, nrestrict:u32, ccount:u32, fcount:u32, r01:f32, rfc:f32, w:f32, pad:f32 }; struct GLink { node:u32, far:u32, i:u32, ib:u32, kind:u32, q:f32, p0:u32, p1:u32 }; @@ -186,11 +187,18 @@ fn macros9(fv: array) -> vec3 { return vec3(rho, mx/rho, my/rho); } -// Проекция неравновесия на сдвиговую часть модели KBC D: моменты N = M20-M02 и Pi_xy = M11 -fn shift9(d: array) -> array { +// Проекция неравновесия на сдвиговую часть. model=0: {N, Pi_xy} (N1); model=1: плюс след T (N2) +fn shift9(d: array, model: u32) -> array { let a = 0.25*(d[1] + d[3] - d[2] - d[4]); let b = 0.25*(d[5] + d[7] - d[6] - d[8]); - return array(0.0, a, -a, a, -a, b, -b, b, -b); + var s = array(0.0, a, -a, a, -a, b, -b, b, -b); + if (model == 1u) { + let dt = d[1] + d[2] + d[3] + d[4] + 2.0*(d[5] + d[6] + d[7] + d[8]); + let q = 0.25*dt; + s[0] = s[0] - dt; + s[1] = s[1] + q; s[2] = s[2] + q; s[3] = s[3] + q; s[4] = s[4] + q; + } + return s; } fn load9(base: ptr>, off: u32, n: u32) { @@ -200,7 +208,7 @@ fn load9(base: ptr>, off: u32, n: u32) { // Возвращает (пост-столкновительные популяции, gamma, признак вырождения) struct CollOut { fv: array, gamma: f32, degen: f32 }; -fn collide9(fin: array, b: f32, op: u32) -> CollOut { +fn collide9(fin: array, b: f32, op: u32, model: u32) -> CollOut { var out: CollOut; // WGSL разрешает переменный индекс только по памяти (var), а не по значению (let/параметр), // поэтому всё, что индексируется в цикле, кладётся в var @@ -215,7 +223,7 @@ fn collide9(fin: array, b: f32, op: u32) -> CollOut { } var d: array; for (var i = 0u; i < 9u; i = i + 1u) { d[i] = fv[i] - fe[i]; } - var ds = shift9(d); + var ds = shift9(d, model); var num = 0.0; var den = 0.0; var nrm = 0.0; for (var i = 0u; i < 9u; i = i + 1u) { let inv = 1.0 / fe[i]; @@ -273,7 +281,7 @@ fn k_collide(@builtin(global_invocation_id) gid: vec3) { gam[nd] = 2.0; return; } - var r = collide9(fv, beta[nd % P.nx], D.collision); + var r = collide9(fv, beta[nd % P.nx], D.collision, D.kbc_model); for (var i = 0u; i < 9u; i = i + 1u) { post[i*P.n + nd] = r.fv[i]; } gam[nd] = r.gamma; } @@ -422,7 +430,9 @@ fn k_stats1(@builtin(global_invocation_id) gid: vec3, p.degen = select(0.0, 1.0, g == 2.0); // объёмная вязкость модели D: xi = cs^2 (1/(gamma*beta) - 1/2) let b = beta[nd % P.nx]; - let xi = CS2*(1.0/(g*b) - 0.5); + // при следе в сдвиговой части (N2) объёмная вязкость равна сдвиговой и всегда > 0 + var xi = CS2*(1.0/(b + b) - 0.5); + if (D.kbc_model == 0u) { xi = CS2*(1.0/(g*b) - 0.5); } p.xineg = select(0.0, 1.0, xi < 0.0); } // зонд следа снимается с того уровня, на котором он лежит @@ -596,8 +606,8 @@ impl GpuLevel { } } -fn soa_equilibrium(n: usize) -> Vec { - let fe = math::feq(1.0, 0.0, 0.0); +fn soa_equilibrium(n: usize, u0: (R, R)) -> Vec { + let fe = math::feq(1.0, u0.0, u0.1); let mut v = vec![0.0f32; 9 * n]; for i in 0..9 { for k in 0..n { @@ -617,9 +627,10 @@ fn make_level( beta: &[f32], probe_node: u32, flags: u32, + u0: (R, R), ) -> GpuLevel { let n = nx * ny; - let init = soa_equilibrium(n); + let init = soa_equilibrium(n, u0); let mkf = |label: &str| { device.create_buffer_init(&wgpu::util::BufferInitDescriptor { label: Some(label), @@ -861,6 +872,14 @@ impl Sim { }) .collect(); + // стартовое поле: либо сразу набегающий поток, либо покой (см. cpu.rs) + let u0 = if spec.init_uniform { + let (sn, cs) = spec.flow_angle.sin_cos(); + (spec.units.u_lat * cs, spec.units.u_lat * sn) + } else { + (0.0, 0.0) + }; + // ── зонд ── let (px, py) = spec.probe; let probe_on_fine = matches!(spec.patch, Some((ax, bx, ay, by)) @@ -879,6 +898,7 @@ impl Sim { &beta0, probe0 as u32, flags0, + u0, ); let pre = device.create_buffer(&wgpu::BufferDescriptor { @@ -907,7 +927,7 @@ impl Sim { let probe1 = if probe_on_fine { ((py - ay) * r) * nfx + (px - ax) * r } else { 0 }; let flags1 = if probe_on_fine { 2 } else { 0 }; let lvl1 = - make_level(&device, &bgl_level, nfx, nfy, &geom1, &beta1, probe1 as u32, flags1); + make_level(&device, &bgl_level, nfx, nfy, &geom1, &beta1, probe1 as u32, flags1, u0); let ghosts: Vec = patch .ghosts() @@ -1053,7 +1073,12 @@ impl Sim { fn inlet(&self, t: u64) -> (R, R) { let sp = &self.spec; - let u = sp.units.u_lat * math::smoothstep(t as R / sp.ramp.max(1) as R); + let ramp = if sp.init_uniform { + 1.0 + } else { + math::smoothstep(t as R / sp.ramp.max(1) as R) + }; + let u = sp.units.u_lat * ramp; let (s, c) = sp.flow_angle.sin_cos(); let mut uy = u * s; if sp.pert_dur > 0 && t >= sp.ramp && t < sp.ramp + sp.pert_dur { @@ -1074,7 +1099,7 @@ impl Sim { ux_in: ux_in as f32, uy_in: uy_in as f32, rho_out: 1.0, - _pad0: 0.0, + kbc_model: (self.spec.kbc_model == math::KbcModel::N2) as u32, outlet_extrap: self.spec.outlet_extrapolate as u32, nparts: self.l0.parts_count, refine: refine as u32, diff --git a/docs/theory/2d_solver/src/main.rs b/docs/theory/2d_solver/src/main.rs index 964d883..0dccff5 100644 --- a/docs/theory/2d_solver/src/main.rs +++ b/docs/theory/2d_solver/src/main.rs @@ -88,9 +88,13 @@ pub struct Spec { pub pert_amp: R, pub pert_dur: u64, pub outlet_extrapolate: bool, + /// Начальное поле уже несёт набегающий поток (иначе — покой с разгоном входа). + pub init_uniform: bool, pub sponge_len: usize, pub sponge_mult: R, pub collision: Collision, + /// Что входит в сдвиговую часть s (см. `math::KbcModel`). + pub kbc_model: math::KbcModel, /// Узел зонда следа в координатах L0. pub probe: (usize, usize), } @@ -175,20 +179,32 @@ struct Cli { /// Длительность стартового импульса, шагов #[arg(long, default_value_t = 300, help_heading = "Время")] pert_dur: u64, + /// Начальное поле: uniform — домен сразу заполнен набегающим потоком и вход включён на + /// полную (умолчание); rest — покой, вход разгоняется за --ramp шагов. Старт из покоя + /// закачивает в канал продольную акустическую моду, которая почти не затухает. + #[arg(long, default_value = "uniform", value_parser = ["uniform", "rest"], + help_heading = "Время")] + init: String, // ── численная схема ── /// Оператор столкновения #[arg(long, default_value = "kbc", value_parser = ["kbc", "bgk"], help_heading = "Схема")] collision: String, + /// Состав сдвиговой части KBC (табл. I 2D-статьи): n1 — только девиатор {N, Π_xy} + /// (KBC D), объёмная вязкость гуляет вместе с γ и может стать отрицательной; n2 — + /// девиатор со следом {N, Π_xy, T} (KBC C), объёмная вязкость фиксирована ξ = ν. + #[arg(long, default_value = "n1", value_parser = math::KbcModel::ALL, help_heading = "Схема")] + kbc_model: String, /// Поперечная скорость на выходе: extrapolate — нуль-градиент (умолчание: в факторном /// исследовании эталонного решателя это лучший вариант по всем метрикам), zero — жёсткий /// ноль (классический Zou–He; отражает вихри дорожки назад к телу и завышает rms Cl) #[arg(long, default_value = "extrapolate", value_parser = ["zero", "extrapolate"], help_heading = "Схема")] outlet: String, - /// Длина поглощающей губки перед выходом, столбцов (0 — выключена) - #[arg(long, default_value_t = 0, help_heading = "Схема")] - sponge_len: usize, + /// Длина поглощающей губки перед выходом, столбцов. Не задана — подбирается по домену; + /// 0 — выключить. Губка гасит продольную моду, которую вязкость сама не гасит. + #[arg(long, help_heading = "Схема")] + sponge_len: Option, /// Во сколько раз губка поднимает вязкость #[arg(long, default_value_t = 30.0, help_heading = "Схема")] sponge_mult: R, @@ -322,17 +338,38 @@ fn build_spec(cli: &Cli) -> Result { if !(2 <= ax && ax < bx && bx <= cli.nx - 2) { return Err(format!("патч по x [{ax},{bx}] обязан лежать внутри (1,{})", cli.nx - 1)); } - if cli.sponge_len > 0 && cli.nx - 1 - cli.sponge_len <= bx { - return Err(format!( - "губка (последние {} столбцов) накрывает патч (bx={bx})", - cli.sponge_len - )); - } Some((ax, bx, ay, by)) } else { None }; + // Губка перед выходом. Задача канала «вход по скорости + выход по давлению» акустически + // есть четвертьволновая труба: вход отражает продольные волны как жёсткий поршень, выход — + // как открытый конец. Собственная мода затухает как ν(π/2Nx)², то есть практически не + // затухает, и любой стартовый удар остаётся в домене навсегда. Губка — единственное, что + // её реально ест, поэтому по умолчанию она включена и подобрана под длину домена. + let sponge_len = match cli.sponge_len { + Some(v) => { + if v > 0 { + if let Some((_, bx, _, _)) = patch { + if cli.nx - 1 - v <= bx { + return Err(format!("губка (последние {v} столбцов) накрывает патч (bx={bx})")); + } + } + } + v + } + None => { + let want = (cli.nx / 12).max(16); + // губка не должна доставать до патча (и вообще до тела) + let limit = match patch { + Some((_, bx, _, _)) => cli.nx.saturating_sub(bx + 2), + None => cli.nx / 4, + }; + want.min(limit) + } + }; + let probe = ( ((cx + 3.0 * cli.size).round() as usize).min(cli.nx - 2), (cy.round() as usize).min(cli.ny - 2), @@ -352,9 +389,11 @@ fn build_spec(cli: &Cli) -> Result { pert_amp: cli.pert_amp, pert_dur: cli.pert_dur, outlet_extrapolate: cli.outlet == "extrapolate", - sponge_len: cli.sponge_len, + init_uniform: cli.init == "uniform", + sponge_len, sponge_mult: cli.sponge_mult, collision: if cli.collision == "bgk" { Collision::Bgk } else { Collision::Kbc }, + kbc_model: math::KbcModel::from_str(&cli.kbc_model).ok_or("неизвестная модель KBC")?, probe, }) } @@ -618,21 +657,52 @@ fn print_header(spec: &Spec, cli: &Cli, plan: &gif::GifPlan, field: FieldKind) { } None => println!(" L1 {:>12} (измельчение выключено)", "нет"), } - println!(" оператор {:>12}", if spec.collision == Collision::Kbc { "KBC-D" } else { "LBGK" }); + println!(" оператор {:>12}", + if spec.collision == Collision::Kbc { + match spec.kbc_model { + math::KbcModel::N1 => "KBC-N1 (D)", + math::KbcModel::N2 => "KBC-N2 (C)", + } + } else { + "LBGK" + }); println!(" выход по u_y {:>12} губка {} столбцов ×{:.0}", if spec.outlet_extrapolate { "extrapolate" } else { "zero" }, spec.sponge_len, spec.sponge_mult); + // Продольная акустика — не мелочь оформления, а то, что видно глазом: при заметной + // амплитуде поток «дышит» рывками. Печатаем период моды и то, чем она гасится. + { + let cs = math::CS2.sqrt(); + let period = 4.0 * spec.nx as R / cs; // четвертьволновая труба: lambda = 4L + let nu = math::nu_of_beta(spec.beta0); + let k = std::f64::consts::PI / (2.0 * spec.nx as R); + let efold = 1.0 / (nu * k * k); + println!(" +── Продольная акустика канала ──────────────────────────────────────────────"); + println!(" вход по скорости отражает звук как жёсткий поршень, выход по давлению — как"); + println!(" открытый конец: канал работает четвертьволновым резонатором."); + println!(" период основной моды {:>10.0} шагов = {:.4} с", period, period * spec.units.dt); + println!(" затухание вязкостью {:>10.0} шагов на е-кратное — то есть сама не гаснет", efold); + print!(" чем гасится: "); + if spec.init_uniform { + print!("старт из однородного потока (моду не возбуждаем)"); + } else { + print!("РАЗГОН ИЗ ПОКОЯ за {} шагов", spec.ramp); + if (spec.ramp as R) < period { + print!(" — короче периода моды, будет звон!"); + } + } + if spec.sponge_len > 0 { + println!(" + губка {} столбцов ×{:.0}", spec.sponge_len, spec.sponge_mult); + } else { + println!("; губка ВЫКЛЮЧЕНА"); + } + } + // Пара «вход по скорости / выход по давлению» — недодемпфированный акустический резонатор: // затухание продольной моды ~ ν(π/Nx)², то есть на длинном домене её почти ничто не гасит, // и стартовый транзиент способен раскачать её до смещённого режима или до развала счёта. - if spec.sponge_len == 0 && spec.nx >= 250 { - println!(" - ⚠ домен длинный ({} столбцов), губка выключена. Продольная акустическая", spec.nx); - println!(" мода затухает как ν(π/Nx)² и на такой длине почти не гасится: возможен уход"); - println!(" ⟨ρ⟩ от единицы (режим смещён, числа несопоставимы) вплоть до развала счёта."); - println!(" Смотрите ⟨ρ⟩ в отчёте о сходимости; лечится ключом --sponge-len 32."); - } if cli.gif.is_some() { println!("\n── Анимация (синхронизация с физическим временем) ───────────────────────────"); @@ -791,6 +861,31 @@ fn final_report( println!("rms Cl насыщен ({:+.4}); дрейф-устойчивый Cd = {cdn1:.3}", clr1 - clr0); } + // ── измеренная продольная пульсация ── + // Колебание средней плотности — это и есть та самая «поршневая» мода. Переводим его в + // амплитуду скорости: для бегущей акустической волны u' = (δρ/ρ)·c_s. + { + let cs = math::CS2.sqrt(); + let last = &rho[n - n / 10..]; + let lo = last.iter().cloned().fold(R::INFINITY, R::min); + let hi = last.iter().cloned().fold(R::NEG_INFINITY, R::max); + let mean = last.iter().sum::() / last.len() as R; + let u_ac = 0.5 * (hi - lo) / mean * cs; + let ratio = u_ac / u.u_lat; + println!(" +── Продольная пульсация (последняя десятая часть прогона) ───────────────────"); + println!(" размах ⟨ρ⟩ = {:.5} вокруг {mean:.5} ⇒ амплитуда скорости u' = {:.2}% от U", + hi - lo, ratio * 100.0); + if ratio > 0.05 { + println!(" ⚠ ПОТОК ЗАМЕТНО ПУЛЬСИРУЕТ: это продольная мода канала, а не физика следа."); + println!(" Лечится стартом из однородного потока (--init uniform) и губкой"); + println!(" (--sponge-len {}). Числа Cd/St на таком прогоне недостоверны.", + (spec.nx / 12).max(16)); + } else { + println!(" пульсация в пределах нормы — продольная мода не раскачана"); + } + } + // ── диагностика KBC ── if spec.collision == Collision::Kbc { let (g_m, g_s) = math::mean_std(&recs[h..].iter().map(|r| r.gamma_mean).collect::>()); diff --git a/docs/theory/2d_solver/src/math.rs b/docs/theory/2d_solver/src/math.rs index c6ee540..d95364f 100644 --- a/docs/theory/2d_solver/src/math.rs +++ b/docs/theory/2d_solver/src/math.rs @@ -86,6 +86,33 @@ pub fn feq(rho: R, ux: R, uy: R) -> [R; Q] { ] } +/// Приближение Града — ур. (58) 2D-статьи. Восстанавливает популяции по ρ, u и ПОЛНОМУ +/// тензору давлений Π (то есть вместе с неравновесной частью Π⁽¹⁾): +/// +/// f*_i = W_i [ ρ + ρu_α v_iα/c_s² + (1/2c_s⁴)(Π_αβ − ρc_s²δ_αβ)(v_iα v_iβ − c_s²δ_αβ) ] +/// +/// Нужно для корректного СТАРТА: инициализация одним равновесием обнуляет Π⁽¹⁾, схема тратит +/// несколько шагов на её восстановление, и этот стартовый дефект портит наблюдаемый порядок +/// сходимости. Статья инициализирует эталонные течения именно так. +/// +/// По построению сохраняет ρ и ρu точно (третий момент весов обнуляется по симметрии). +/// +/// В самом решателе пока не используется: канальная постановка стартует с однородного потока, +/// где Π⁽¹⁾ = 0 и приближение Града вырождается в равновесие. Нужна для эталонных течений, +/// у которых стартовое поле имеет ненулевые градиенты. +#[cfg_attr(not(test), allow(dead_code))] +#[inline] +pub fn grad_init(rho: R, ux: R, uy: R, pxx: R, pxy: R, pyy: R) -> [R; Q] { + let axx = pxx - rho * CS2; + let ayy = pyy - rho * CS2; + std::array::from_fn(|i| { + let cx = CX[i] as R; + let cy = CY[i] as R; + let quad = axx * (cx * cx - CS2) + 2.0 * pxy * cx * cy + ayy * (cy * cy - CS2); + W[i] * (rho + rho * (ux * cx + uy * cy) / CS2 + quad / (2.0 * CS2 * CS2)) + }) +} + /// Макропеременные: ρ = Σ f_i, u = (Σ c_i f_i)/ρ. #[inline] pub fn macros(f: &[R; Q]) -> (R, R, R) { @@ -99,22 +126,58 @@ pub fn macros(f: &[R; Q]) -> (R, R, R) { // KBC: проектор на сдвиг, энтропийный стабилизатор, столкновение // ───────────────────────────────────────────────────────────────────────────── -/// Δs = Ps·Δ — проекция неравновесия на сдвиговую часть модели KBC D: натуральные моменты -/// N = M20 − M02 и Π_xy = M11. +/// Состав сдвиговой части s — то, что релаксирует жёстко с 2β. Всё остальное попадает в h и +/// отдаётся энтропийному стабилизатору γ. Выбор берётся из табл. I 2D-статьи; в 3D-работах +/// те же варианты называются N1 и N2. +#[derive(Clone, Copy, Debug, PartialEq, Eq)] +pub enum KbcModel { + /// s = {N, Π_xy} — только девиатор напряжений. 2D-статья: **KBC D**; 3D: **KBC-N1**. + /// След T уходит в h, поэтому объёмная вязкость ξ = c_s²(1/(γβ) − ½) гуляет вместе с γ + /// и способна уходить в минус — локальное антизатухание звука (ур. 57). + N1, + /// s = {N, Π_xy, T} — девиатор плюс след. 2D-статья: **KBC C**; 3D: **KBC-N2**. + /// Объёмная вязкость фиксирована и положительна: ξ = ν (ур. 57). Полезно там, где в + /// задаче есть слабозатухающая акустика. + N2, +} + +impl KbcModel { + pub fn from_str(s: &str) -> Option { + Some(match s.to_ascii_lowercase().as_str() { + "n1" | "d" => KbcModel::N1, + "n2" | "c" => KbcModel::N2, + _ => return None, + }) + } + pub const ALL: [&'static str; 2] = ["n1", "n2"]; +} + +/// Δs = Ps·Δ — проекция неравновесия на сдвиговую часть выбранной модели. /// -/// Аналитическая форма, прямо из представления популяций через натуральные моменты (10): -/// вклад N сидит только в f(σ,0) (+ρN/4) и f(0,λ) (−ρN/4), вклад Π_xy — только в f(σ,λ) -/// (+σλ·ρΠ_xy/4). Поскольку n_i и p_i линейны по своим моментам, а ρ у f и feq одинакова, -/// разности берутся покомпонентно. Эквивалентно матричному Ps = M⁻¹·diag(0,0,0,0,1,1,0,0,0)·M -/// (проверяется тестом `projector_matches_moment_matrix`), но без матрицы 9×9. +/// Аналитическая форма, прямо из представления популяций через натуральные моменты (ур. 10). +/// Вклад N = M20 − M02 сидит только в f(σ,0) (+ρN/4) и f(0,λ) (−ρN/4); вклад Π_xy = M11 — +/// только в f(σ,λ) (+σλ·ρΠ_xy/4); вклад следа T = M20 + M02 — в f(0,0) (−ρT) и в f(σ,0), +/// f(0,λ) (по +ρT/4). Все три вклада линейны по своим моментам, а ρ у f и feq одинакова, +/// поэтому разности берутся покомпонентно. Эквивалентно матричному +/// Ps = M⁻¹·diag(оставить нужные моменты)·M, но без матрицы 9×9. #[inline] -pub fn project_shift(d: &[R; Q]) -> [R; Q] { +pub fn project_shift(d: &[R; Q], model: KbcModel) -> [R; Q] { // ρΔN = Σ (c_x² − c_y²)Δ_i и ρΔΠ_xy = Σ c_x c_y Δ_i let dn = d[1] + d[3] - d[2] - d[4]; let dp = d[5] + d[7] - d[6] - d[8]; let a = 0.25 * dn; let b = 0.25 * dp; - [0.0, a, -a, a, -a, b, -b, b, -b] + let mut s = [0.0, a, -a, a, -a, b, -b, b, -b]; + if let KbcModel::N2 = model { + // ρΔT = Σ (c_x² + c_y²)Δ_i + let dt = d[1] + d[2] + d[3] + d[4] + 2.0 * (d[5] + d[6] + d[7] + d[8]); + s[0] -= dt; + let q = 0.25 * dt; + for v in s.iter_mut().take(5).skip(1) { + *v += q; + } + } + s } /// Что вернуло столкновение помимо новых популяций — нужно для диагностики KBC. @@ -137,7 +200,7 @@ pub struct Kbc { /// При γ = 2 совпадает с LBGK. Сдвиговые моменты релаксируют с точным 2β = 1/τ при любом γ, /// поэтому кинематическая вязкость (5) от стабилизатора не зависит. #[inline] -pub fn collide_node(f: &mut [R; Q], beta: R) -> Kbc { +pub fn collide_node(f: &mut [R; Q], beta: R, model: KbcModel) -> Kbc { let (rho, ux, uy) = macros(f); let fe = feq(rho, ux, uy); @@ -145,7 +208,7 @@ pub fn collide_node(f: &mut [R; Q], beta: R) -> Kbc { for i in 0..Q { d[i] = f[i] - fe[i]; } - let ds = project_shift(&d); + let ds = project_shift(&d, model); let (mut num, mut den, mut nrm) = (0.0, 0.0, 0.0); for i in 0..Q { @@ -197,8 +260,13 @@ pub fn beta_of_nu(nu: R) -> R { /// Отрицательная ξ означает локальное антизатухание акустической моды — поэтому в отчёте /// печатается доля узлов с ξ < 0 (см. report в main.rs). #[inline] -pub fn xi_of_gamma(gamma: R, beta: R) -> R { - CS2 * (1.0 / (gamma * beta) - 0.5) +pub fn xi_of_gamma(gamma: R, beta: R, model: KbcModel) -> R { + match model { + // след в h — ξ зависит от стабилизатора и может стать отрицательной + KbcModel::N1 => CS2 * (1.0 / (gamma * beta) - 0.5), + // след в s — ξ равна сдвиговой вязкости и всегда положительна + KbcModel::N2 => nu_of_beta(beta), + } } // ───────────────────────────────────────────────────────────────────────────── @@ -737,17 +805,18 @@ mod tests { use super::*; /// Проектор Δs, построенный аналитически, обязан совпадать с матричным - /// Ps = M⁻¹·diag(…,1,1,…)·M, где M — базис натуральных моментов из (6)–(7). + /// Ps = M⁻¹·D·M, где M — базис натуральных моментов (6)–(7), а D оставляет ровно те + /// моменты, которые входят в сдвиговую часть выбранной модели (табл. I). #[test] fn projector_matches_moment_matrix() { - // строки M: 1, cx, cy, 3(cx²+cy²)−2, cx²−cy², cx·cy, cx²cy, cx cy², cx²cy² + // строки M: 1, cx, cy, T=cx²+cy², N=cx²−cy², Π_xy=cx·cy, cx²cy, cx cy², cx²cy² let mut m = [[0.0f64; Q]; Q]; for i in 0..Q { let (x, y) = (CX[i] as f64, CY[i] as f64); m[0][i] = 1.0; m[1][i] = x; m[2][i] = y; - m[3][i] = 3.0 * (x * x + y * y) - 2.0; + m[3][i] = x * x + y * y; m[4][i] = x * x - y * y; m[5][i] = x * y; m[6][i] = x * x * y; @@ -755,29 +824,62 @@ mod tests { m[8][i] = x * x * y * y; } let minv = invert9(&m); - // Ps = M⁻¹ · D · M, D оставляет строки 4 (N) и 5 (Π_xy) - let mut ps = [[0.0f64; Q]; Q]; - for a in 0..Q { - for b in 0..Q { - ps[a][b] = minv[a][4] * m[4][b] + minv[a][5] * m[5][b]; - } - } - // сверяем на случайных (детерминированных) векторах - let mut seed = 12345u64; - let mut rnd = || { - seed = seed.wrapping_mul(6364136223846793005).wrapping_add(1442695040888963407); - ((seed >> 33) as f64 / (1u64 << 31) as f64) - 1.0 - }; - for _ in 0..200 { - let d: [f64; Q] = std::array::from_fn(|_| rnd()); - let got = project_shift(&d); + for (model, keep) in + [(KbcModel::N1, vec![4usize, 5]), (KbcModel::N2, vec![3usize, 4, 5])] + { + let mut ps = [[0.0f64; Q]; Q]; for a in 0..Q { - let want: f64 = (0..Q).map(|b| ps[a][b] * d[b]).sum(); - assert!((got[a] - want).abs() < 1e-13, "i={a}: {} vs {}", got[a], want); + for b in 0..Q { + ps[a][b] = keep.iter().map(|&k| minv[a][k] * m[k][b]).sum(); + } + } + let mut seed = 12345u64; + let mut rnd = || { + seed = seed.wrapping_mul(6364136223846793005).wrapping_add(1442695040888963407); + ((seed >> 33) as f64 / (1u64 << 31) as f64) - 1.0 + }; + for _ in 0..200 { + // проектор сравниваем на ФИЗИЧЕСКИ допустимых Δ: у них моменты массы и + // импульса нулевые (и f, и feq несут одни и те же ρ, ρu) + let mut d: [f64; Q] = std::array::from_fn(|_| rnd()); + let mass: f64 = d.iter().sum(); + let mx: f64 = (0..Q).map(|i| CX[i] as f64 * d[i]).sum(); + let my: f64 = (0..Q).map(|i| CY[i] as f64 * d[i]).sum(); + for i in 0..Q { + d[i] -= W[i] * (mass + 3.0 * (CX[i] as f64 * mx + CY[i] as f64 * my)); + } + let got = project_shift(&d, model); + for a in 0..Q { + let want: f64 = (0..Q).map(|b| ps[a][b] * d[b]).sum(); + assert!( + (got[a] - want).abs() < 1e-13, + "{model:?} i={a}: {} vs {want}", + got[a] + ); + } } } } + /// Объёмная вязкость по ур. (57): у модели со следом в сдвиговой части она равна + /// сдвиговой и не зависит от стабилизатора, у модели без следа — зависит. + #[test] + fn bulk_viscosity_follows_model() { + let beta = 1.0 / (2.0 * 0.7); + let nu = nu_of_beta(beta); + for &g in &[0.5, 1.0, 1.7, 3.0] { + assert!( + (xi_of_gamma(g, beta, KbcModel::N2) - nu).abs() < 1e-15, + "N2: ξ обязана равняться ν при любой γ" + ); + let want = CS2 * (1.0 / (g * beta) - 0.5); + assert!((xi_of_gamma(g, beta, KbcModel::N1) - want).abs() < 1e-15); + } + // у N1 ξ уходит в минус при достаточно большой γ — то самое антизатухание звука + let g_neg = 2.0 / (beta * CS2) * CS2; + assert!(xi_of_gamma(g_neg * 2.0, beta, KbcModel::N1) < 0.0); + } + /// Ps обязан быть идемпотентным (это проектор) и не трогать сохраняющиеся моменты. #[test] fn projector_is_idempotent_and_conservative() { @@ -788,8 +890,9 @@ mod tests { }; for _ in 0..100 { let d: [f64; Q] = std::array::from_fn(|_| rnd()); - let s1 = project_shift(&d); - let s2 = project_shift(&s1); + for model in [KbcModel::N1, KbcModel::N2] { + let s1 = project_shift(&d, model); + let s2 = project_shift(&s1, model); for i in 0..Q { assert!((s1[i] - s2[i]).abs() < 1e-14); } @@ -798,6 +901,7 @@ mod tests { let mx: f64 = (0..Q).map(|i| CX[i] as f64 * s1[i]).sum(); let my: f64 = (0..Q).map(|i| CY[i] as f64 * s1[i]).sum(); assert!(mx.abs() < 1e-13 && my.abs() < 1e-13); + } } } @@ -823,7 +927,7 @@ mod tests { let (rho, ux, uy) = macros(&f); let fe = feq(rho, ux, uy); let d: [f64; Q] = std::array::from_fn(|i| f[i] - fe[i]); - let ds = project_shift(&d); + let ds = project_shift(&d, KbcModel::N1); let dh: [f64; Q] = std::array::from_fn(|i| d[i] - ds[i]); // невязка условия (15) @@ -835,7 +939,7 @@ mod tests { .sum() }; let mut fc = f; - let info = collide_node(&mut fc, beta); + let info = collide_node(&mut fc, beta, KbcModel::N1); let r0 = resid(info.gamma).abs(); let r1 = resid(info.gamma + 1.0).abs().min(resid(info.gamma - 1.0).abs()); assert!(r0 < 0.05 * r1, "невязка (15) при γ*: {r0:.3e}, на ±1 от него: {r1:.3e}"); @@ -850,10 +954,10 @@ mod tests { let (rho, ux, uy) = macros(&f); let fe = feq(rho, ux, uy); let d: [f64; Q] = std::array::from_fn(|i| f[i] - fe[i]); - let ds = project_shift(&d); + let ds = project_shift(&d, KbcModel::N1); let mut pure_shear: [f64; Q] = std::array::from_fn(|i| fe[i] + ds[i]); let mut bgk = pure_shear; - collide_node(&mut pure_shear, beta); + collide_node(&mut pure_shear, beta, KbcModel::N1); collide_node_bgk(&mut bgk, beta); for i in 0..Q { assert!((pure_shear[i] - bgk[i]).abs() < 1e-12, "i={i}"); @@ -864,20 +968,26 @@ mod tests { /// что вязкость (5) задаётся только β, а стабилизатор её не трогает. #[test] fn shear_relaxes_at_two_beta_regardless_of_gamma() { - for &beta in &[0.3, 0.6, 0.95] { - let mut f = perturbed_state(0.8); - let (rho, ux, uy) = macros(&f); - let fe = feq(rho, ux, uy); - let d0: [f64; Q] = std::array::from_fn(|i| f[i] - fe[i]); - let s0 = project_shift(&d0); - collide_node(&mut f, beta); - let (r1, x1, y1) = macros(&f); - let fe1 = feq(r1, x1, y1); - let d1: [f64; Q] = std::array::from_fn(|i| f[i] - fe1[i]); - let s1 = project_shift(&d1); - for i in 0..Q { - let want = (1.0 - 2.0 * beta) * s0[i]; - assert!((s1[i] - want).abs() < 1e-10, "β={beta} i={i}: {} vs {want}", s1[i]); + for model in [KbcModel::N1, KbcModel::N2] { + for &beta in &[0.3, 0.6, 0.95] { + let mut f = perturbed_state(0.8); + let (rho, ux, uy) = macros(&f); + let fe = feq(rho, ux, uy); + let d0: [f64; Q] = std::array::from_fn(|i| f[i] - fe[i]); + let s0 = project_shift(&d0, model); + collide_node(&mut f, beta, model); + let (r1, x1, y1) = macros(&f); + let fe1 = feq(r1, x1, y1); + let d1: [f64; Q] = std::array::from_fn(|i| f[i] - fe1[i]); + let s1 = project_shift(&d1, model); + for i in 0..Q { + let want = (1.0 - 2.0 * beta) * s0[i]; + assert!( + (s1[i] - want).abs() < 1e-10, + "{model:?} β={beta} i={i}: {} vs {want}", + s1[i] + ); + } } } } @@ -888,7 +998,7 @@ mod tests { for &beta in &[0.25, 0.5, 0.99] { let mut f = perturbed_state(0.7); let before = macros(&f); - collide_node(&mut f, beta); + collide_node(&mut f, beta, KbcModel::N1); let after = macros(&f); assert!((before.0 - after.0).abs() < 1e-13); assert!((before.1 - after.1).abs() < 1e-12);