Пристеночная функция Spalding, многоуровневый домен, группы J и K
Решатель: - --wall-function spalding поверх bouzidi/grad/hrr: моментные схемы берут модельную скорость узла и u_τ²/ν в тензор давлений, Bouzidi — скольжение стенки, согласованное по напряжению (в вязком подслое тождественно ноль); - --levels: цепочка вложенных уровней с рекурсивным шагом на CPU и GPU, одиночный патч — её частный случай, прежние постановки побитово те же; - --sponge-side, β по узлам; в сводке y⁺, cl_mean, rho_probe_rms. - Исправлено: на GPU при --refine >= 3 сила лишних подшагов терялась за краем буфера, Cd занижался в r/2 раза. Кампания: группы J (65, схемы стенки с функцией) и K (25, внешний домен), 205 прогонов на 5 долей по ≈27 ч; bench/compare.py сводит J и K против эталонов; parity.py сверяет Spalding и три уровня. Образ 1.4.0. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
This commit is contained in:
Claude Opus 5.5
+124
-11
@@ -311,8 +311,11 @@ c_s²(1/(γβ) − ½) и, поскольку измеренная ⟨γ⟩ ≈
|
||||
низкими числами Рейнольдса, поскольку на границе возникают паразитные скачки») и естественная
|
||||
форма для подвижных стенок: скорость стенки входит в целевые значения, а не отдельной поправкой.
|
||||
В здешней канальной постановке преимущества по устойчивости воспроизвести не удалось: при росте
|
||||
Re обе модели теряют счёт на одном и том же значении (Re ≈ 5·10⁴ при теле в 16 клеток), то есть
|
||||
ограничивает не стенка, а что-то другое — вероятнее всего Zou–He при τ → ½.
|
||||
Re обе модели теряют счёт на одном и том же значении, то есть ограничивает не стенка, а
|
||||
что-то другое — вероятнее всего Zou–He при τ → ½. (Сама оценка «Re ≈ 5·10⁴ при теле в 16
|
||||
клеток» была сделана по коротким прогонам; длинные дают предел τ − ½ ≈ 6·10⁻⁴, то есть при
|
||||
D = 16 уже Re ≈ 4000, — см. раздел «Дальше». Вывод «стенка не ограничивает» он не меняет:
|
||||
Bouzidi и HRR со Spalding разваливаются при одних и тех же Re.)
|
||||
|
||||
**Умолчание — `hrr`, и это решение временное.** Оно принято по устройству схемы (третий порядок
|
||||
Эрмита фильтрует высокочастотный мусор у стенки, чего обрыв на Π не делает), а не по здешним
|
||||
@@ -326,6 +329,59 @@ Re обе модели теряют счёт на одном и том же зн
|
||||
это только потому, что две модели дали побитово одинаковый результат там, где обязаны были
|
||||
разойтись.
|
||||
|
||||
### Пристеночная функция Spalding
|
||||
|
||||
Ключ `--wall-function none|spalding` (умолчание `none` — всё прежнее не меняется побитово)
|
||||
ставит поверх `bouzidi`, `grad` или `hrr` закон стенки Spalding (1961): одна формула на весь
|
||||
пристеночный слой, от вязкого подслоя (u⁺ = y⁺) до логарифмической зоны, κ = 0.41, B = 5.2.
|
||||
С `staircase` не сочетается: ступенька не знает, где стенка.
|
||||
|
||||
Как устроено:
|
||||
|
||||
- у каждого граничного узла есть нормаль (градиент SDF), расстояние до стенки y_w и точка
|
||||
отбора скорости на одну клетку дальше по нормали (Malaspinas, Sagaut 2014); скорость там —
|
||||
билинейно по полю на момент t. Геометрия строится один раз в `cpu::Geom` и загружается на
|
||||
GPU как есть;
|
||||
- по касательной скорости в точке отбора Ньютоном находится u_τ, по нему — модельная скорость
|
||||
узла u_node = u_τ·u⁺(y⁺_w) и напряжение τ_w = ρu_τ²;
|
||||
- **`grad`/`hrr`**: касательная составляющая целевой скорости (B 1) заменяется на u_node, а в
|
||||
тензоре давлений нормальная производная касательной скорости — на u_τ²/ν. Плотность (B 3) и
|
||||
сборка популяций не меняются;
|
||||
- **`bouzidi`**: стенка получает скорость скольжения u_s = u_node − u_τ²·y_w/ν и входит в
|
||||
отскок поправкой подвижной стенки. Скольжение согласовано по **напряжению**: отскок передаёт
|
||||
стенке ρν(u_node − u_s)/y_w = τ_w. В вязком подслое u⁺ = y⁺, и u_s тождественно ноль; в
|
||||
буферной и логарифмической зонах u_s < 0 — стенка «отступает», добирая трение до τ_w.
|
||||
Модуль скольжения зажат скоростью в точке отбора.
|
||||
|
||||
Первая версия согласовывала скольжение по **скорости** (прямая через стенку, узел и точку
|
||||
отбора). В решателе без турбулентной вязкости это даёт большое положительное скольжение —
|
||||
стенка становится почти свободной: цилиндр Re = 10⁴, D = 16 выдал Cd = 0.016. Вторая брала
|
||||
разрешённую скорость узла вместо модельной и при Re = 20 расходилась с прилипанием на 2 %:
|
||||
у искривлённой стенки узел не лежит на прямой от стенки к точке отбора. Нынешняя формула
|
||||
свободна от обоих изъянов.
|
||||
|
||||
**Сведение к прилипанию при малом y⁺** (цилиндр Re = 20, D = 16, 320×192, y⁺ ≈ 0.25, GPU):
|
||||
|
||||
| модель | прилипание | Spalding | разница |
|
||||
|---|---|---|---|
|
||||
| `bouzidi` | 2.38252 | 2.38252 | 0 (до 7-го знака) |
|
||||
| `grad` | 2.38774 | 2.38989 | +0.09 % |
|
||||
| `hrr` | 2.38835 | 2.39053 | +0.09 % |
|
||||
|
||||
У моментных схем остаточные 0.09 % — от того, что целевая скорость строится теперь по нормали
|
||||
к стенке, а не усреднением по линкам (B 1).
|
||||
|
||||
**Закон стенки описывает турбулентный слой, а при умеренных Re он ламинарный.** На цилиндре при
|
||||
Re = 2000 (D = 16, y⁺ ≈ 5, местами до 19, 16 тыс. шагов) функция заметно двигает только Bouzidi:
|
||||
Cd 2.10 против 1.56 без неё; у `grad` 1.59 против 1.65, у `hrr` 1.70 против 1.69. Bouzidi
|
||||
навязывает модельное трение через обмен импульсом жёстко, моментные схемы — через тензор
|
||||
давлений одного слоя узлов, мягче. Какая схема ближе к эталону — меряет группа J кампании.
|
||||
|
||||
Каждый прогон канала пишет в сводку y⁺ первого узла на самом тонком уровне (`yplus_mean`,
|
||||
`yplus_max`, долю узлов с y⁺ > 30 `yplus_frac_log`, `u_tau_mean`) — два десятка снимков во
|
||||
второй половине прогона, с функцией и без неё. Без y⁺ сравнение схем стенки не
|
||||
интерпретировать.
|
||||
|
||||
### Паритет бэкендов и согласованность уровней
|
||||
|
||||
CPU (f64) и GPU (f32) на одной постановке совпадают до 4–5 значащих цифр шаг в шаг: ⟨ρ⟩
|
||||
@@ -335,6 +391,18 @@ CPU (f64) и GPU (f32) на одной постановке совпадают
|
||||
Один и тот же случай с патчем ×2 и вовсе без измельчения (`--refine 1`) даёт St 0.1951 против
|
||||
0.1970 и ⟨Cd⟩ 1.956 против 1.943 — связка уровней систематики не вносит.
|
||||
|
||||
**Исправлено: сила на GPU при `--refine 3` и выше.** Рабочие слоты силы в буфере итогов были
|
||||
рассчитаны ровно на два подшага тонкого уровня; сила третьего и следующих писалась за край
|
||||
буфера и терялась, а делилась сумма всё равно на r. Cd на GPU выходил занижен в r/2 раза:
|
||||
два тела при `--refine 3` — 1.316 против 1.972 на CPU. Теперь слотов столько, сколько подшагов
|
||||
самого тонкого уровня на шаг L0: 1.9733 против 1.9716. Это затрагивало прогон `F09_refine3`
|
||||
кампании и любой ручной запуск с `--refine ≥ 3` на GPU; при `--refine 2` и без патча
|
||||
результат прежний побитово.
|
||||
|
||||
`bench/parity.py` сверяет CPU и GPU на четырёх случаях: Тейлор–Грин, стационарный цилиндр,
|
||||
цилиндр с `grad` + Spalding и три уровня ×2×2. Замерено на Iris Xe (Vulkan): Cd совпадает до
|
||||
1.7·10⁻⁵ (Spalding) и 7.3·10⁻⁵ (три уровня), y⁺ — до 6.4·10⁻⁵.
|
||||
|
||||
**Предел на размер сетки снят.** У GPU есть жёсткий предел `maxComputeWorkgroupsPerDimension`
|
||||
= 65535, а диспетчеризация была одномерной: при 64 узлах на рабочую группу это упирало сетку в
|
||||
4.2 миллиона узлов, то есть примерно 2048×2048. Всё, что крупнее, падало ошибкой валидации —
|
||||
@@ -470,15 +538,56 @@ rms Cl 0.81 вместо 0.40 — то есть 75-процентная прод
|
||||
просто переносится во времени: поток всё равно обязан разогнаться вокруг тела, и энергия этого
|
||||
переходного процесса задана физикой, а не гладкостью начального поля. Умолчание — 0.
|
||||
|
||||
## Многоуровневый домен
|
||||
|
||||
`--levels "r:up,down,side; r:up,down,side; …"` строит цепочку вложенных уровней вместо одного
|
||||
патча: L0 — самый грубый, на весь домен, каждый следующий в r раз тоньше родителя и лежит
|
||||
внутри него. `up`, `down`, `side` — отступы границ уровня от центра тела в его диаметрах
|
||||
(вверх по потоку, вниз, в стороны). Так строится внешний домен с грубыми буферами: большой и
|
||||
дешёвый L0, к телу сетка мельчает ступенями. `--levels` перекрывает `--refine`/`--patch`;
|
||||
одиночный патч — частный случай цепочки из одного звена, и прежние постановки дают **побитово
|
||||
тот же** результат на обоих бэкендах (сверено на пяти постановках: без патча, ×2, ×3, два тела,
|
||||
периодический случай; единственное исключение — исправленная сила при `--refine 3` на GPU).
|
||||
|
||||
```sh
|
||||
# L0 в 4 раза грубее тела; ступени ×2: 4D вверх, 15D вниз, ±4D, затем 1.5D, 10D следа, ±2D
|
||||
./target/release/kbc2d --shape cylinder --size 8 --nx 480 --ny 320 --body-x 160 --re 150 \
|
||||
--levels "2:4,15,4; 2:1.5,10,2" --wall bouzidi --steps 96000 --backend gpu
|
||||
```
|
||||
|
||||
Важно: `--size`, `--nx`, `--ny`, `--steps` и Re задаются **в клетках и шагах L0**. Тело в 8
|
||||
клеток L0 при цепочке ×2×2 — это 32 клетки на самом тонком уровне; отчёт печатает D и τ
|
||||
каждого уровня. τ растёт к тонким уровням (τ_k = r_k(τ_{k−1} − ½) + ½, ν одна на всех), так
|
||||
что ближе всего к пределу устойчивости τ → ½ именно L0.
|
||||
|
||||
Шаг рекурсивный: уровень делает шаг, затем его потомок — r подшагов с временной интерполяцией
|
||||
рамки, затем рестрикция обратно; сила снимается на самом тонком уровне на каждом его подшаге.
|
||||
Стенки канала, вход, выход и губки есть только у L0. На GPU это та же рекурсия в одном
|
||||
командном буфере; ядра связки не менялись.
|
||||
|
||||
Проверено: однородный поток без тела — неподвижная точка через три уровня с точностью 10⁻¹²;
|
||||
стационарный цилиндр Re = 20 при одном скачке ×4 и при цепочке ×2×2 с тем же разрешением тела
|
||||
даёт Cd 2.2104 и 2.2156 (0.24 %); CPU и GPU на трёх уровнях совпадают до 7.3·10⁻⁵ по Cd.
|
||||
|
||||
**Боковые губки** `--sponge-side N` — N рядов L0 у верхней и нижней стенок с той же плавной
|
||||
надбавкой вязкости, что у входа и выхода; β теперь поле по узлам, а не по столбцам. Нужны,
|
||||
чтобы отличать отражения от боковых стенок от отражений на стыках уровней.
|
||||
|
||||
В сводку добавлены `levels` (цепочка с размерами), `nodes_per_step` (обновлений узлов за шаг L0
|
||||
по всем уровням), `wall_function`, `sponge_side`, `rho_probe_rms` (пульсация плотности в зонде
|
||||
во второй половине — мера паразитной акустики) и `cl_mean`; в `series.csv` — столбец
|
||||
`rho_probe`.
|
||||
|
||||
## Валидационная кампания
|
||||
|
||||
`bench/` — 115 прогонов на ≈90 машинных часов, разложенных по девяти группам: эталоны
|
||||
`bench/` — 205 прогонов на ≈101 машинный час, разложенных по одиннадцати группам: эталоны
|
||||
первоисточников, цилиндр против литературы, модели стенки, профили крыла, сложная и
|
||||
множественная геометрия, границы домена, старт и время жизни, внутренние инварианты,
|
||||
сверхмелкие сетки до 4096×2048. Из них 108 идут на GPU (≈87 часов) и 7 — на CPU (≈3 часа,
|
||||
исследования сходимости группы A). Каждый прогон кладёт в собственную папку `cmd.txt`,
|
||||
`log.txt`, `report.txt`, `series.csv` и `summary.json`; гифку на всю свою длительность
|
||||
пишут **59 прогонов из 115**, остальным `--gif` не передаётся вовсе.
|
||||
сверхмелкие сетки до 4096×2048, схемы стенки с пристеночной функцией (J) и внешний домен с
|
||||
грубыми буферами (K). Из них 198 идут на GPU и 7 — на CPU (≈3 часа, исследования сходимости
|
||||
группы A). Каждый прогон кладёт в собственную папку `cmd.txt`, `log.txt`, `report.txt`,
|
||||
`series.csv` и `summary.json`; гифку на всю свою длительность пишут **59 прогонов из 205** —
|
||||
группы J и K гифок не пишут. Итоги J и K сводит `bench/compare.py`.
|
||||
|
||||
```sh
|
||||
cd bench
|
||||
@@ -497,7 +606,7 @@ python preflight.py # каждый сценарий старту
|
||||
|
||||
## Развёртывание (Docker)
|
||||
|
||||
Собранный образ опубликован: **`notbigghost/kbc2d:1.3.0`** (он же `latest`, платформа `linux/amd64`).
|
||||
Собранный образ опубликован: **`notbigghost/kbc2d:1.4.0`** (он же `latest`, платформа `linux/amd64`).
|
||||
Исходники на сервере не нужны — достаточно перенести туда один файл
|
||||
`docker-compose.server.yml`:
|
||||
|
||||
@@ -525,7 +634,7 @@ docker run --rm --gpus all kbc2d --calibrate
|
||||
|
||||
### Пять долей на пяти машинах
|
||||
|
||||
Кампания разложена на **пять равных по времени долей** (≈24 ч каждая, поле `shard` в
|
||||
Кампания разложена на **пять равных по времени долей** (≈27 ч каждая, поле `shard` в
|
||||
`scenarios.json`), чтобы считать её одновременно на пяти чужих ПК с видеокартой под Windows
|
||||
или WSL2. На машину переносятся два файла — `docker-compose.shards.yml` и `shard.bat` — и
|
||||
выполняется `shard.bat 3` (или `docker compose -f docker-compose.shards.yml up -d shard3`).
|
||||
@@ -618,7 +727,11 @@ GPU не найден. Согласие даётся явно, переменн
|
||||
стартовый импульс от появления тела в потоке;
|
||||
- подвижные и вращающиеся тела: GMEM уже записан в галилей-инвариантной форме и принимает
|
||||
скорость стенки на линке, но подача этой скорости не подключена;
|
||||
- разбор нерешённых 12% по Cd и предела устойчивости Re ≈ 5·10⁴ — обе задачи вынесены в
|
||||
кампанию (группы B и G), выводы делать по её данным;
|
||||
- разбор нерешённых 12% по Cd и предела устойчивости — обе задачи вынесены в кампанию
|
||||
(группы B и G), выводы делать по её данным. Оценка предела Re ≈ 5·10⁴ при D = 16 была сделана
|
||||
по коротким прогонам и не подтвердилась: в длинных (цилиндр 20D×12D, Bouzidi) Re = 10⁴
|
||||
разваливается при D = 16 и 32 на 13 и 27 тысячах шагов, а держится при τ − ½ ≥ 6·10⁻⁴
|
||||
(Re = 2000 при D = 8, 5000 при D = 32, 10⁴ при D = 64); пристеночная функция предел не
|
||||
сдвигает;
|
||||
- больше четырёх тел в домене: сейчас сила считается по четырём вёдрам (`MAX_BODY_BUCKETS`),
|
||||
геометрия при этом собирается из любого числа тел, но силы сверх четвёртого сливаются вместе.
|
||||
|
||||
@@ -1,10 +1,10 @@
|
||||
# Валидационная кампания
|
||||
|
||||
115 прогонов, ≈90 машинных часов: 108 на GPU (≈87 часов на RTX 4070 Ti) и 7 на CPU (≈3 часа —
|
||||
исследования сходимости группы A, где нужен f64). Проверяет решатель по трём независимым линиям:
|
||||
205 прогонов в одиннадцати группах, ≈101 машинный час: 198 на GPU (≈98 часов при 1200 MLUPS)
|
||||
и 7 на CPU (≈3 часа — исследования сходимости группы A, где нужен f64). Проверяет решатель по трём независимым линиям:
|
||||
эталонам из статей авторов метода, литературе по обтеканию тел и внутренним инвариантам самой
|
||||
схемы. Каждый прогон кладёт в собственную папку `cmd.txt`, `log.txt`, `report.txt`, `series.csv`
|
||||
и `summary.json`; гифку пишут **59 прогонов из 115** — остальным `--gif` не передаётся.
|
||||
и `summary.json`; гифку пишут **59 прогонов из 205** — остальным `--gif` не передаётся.
|
||||
|
||||
## Пять долей: чужие ПК под Windows / WSL2
|
||||
|
||||
@@ -44,11 +44,11 @@ Desktop с бэкендом WSL2 (так по умолчанию) либо Docke
|
||||
|
||||
| доля | прогонов | ≈ часов |
|
||||
|---|---|---|
|
||||
| 1 | 22 | 24.0 |
|
||||
| 2 | 23 | 24.0 |
|
||||
| 3 | 23 | 24.0 |
|
||||
| 4 | 23 | 24.0 |
|
||||
| 5 | 24 | 24.0 |
|
||||
| 1 | 40 | 26.7 |
|
||||
| 2 | 40 | 26.7 |
|
||||
| 3 | 41 | 26.7 |
|
||||
| 4 | 42 | 26.7 |
|
||||
| 5 | 42 | 26.7 |
|
||||
|
||||
Часы — оценка при 1200 MLUPS на крупной сетке **с поправкой на размер сетки**: на dzn мелкие
|
||||
сетки считаются медленнее на узел, а сверхмелкие из группы I идут через раздельные привязки
|
||||
@@ -66,7 +66,7 @@ Desktop с бэкендом WSL2 (так по умолчанию) либо Docke
|
||||
|
||||
## Быстрый старт на сервере
|
||||
|
||||
Образ опубликован, собирать ничего не нужно: **`notbigghost/kbc2d:1.3.0`**. Исходники на
|
||||
Образ опубликован, собирать ничего не нужно: **`notbigghost/kbc2d:1.4.0`**. Исходники на
|
||||
сервере тоже не нужны — переносится один файл `docker-compose.server.yml`.
|
||||
|
||||
```sh
|
||||
@@ -75,7 +75,7 @@ mkdir -p ~/kbc2d && cd ~/kbc2d # сюда же ляжет ./out с ре
|
||||
|
||||
C=docker-compose.server.yml
|
||||
docker compose -f $C --profile check run --rm vulkan # 1. карта видна из контейнера?
|
||||
docker compose -f $C --profile check run --rm preflight # 2. все 115 сценариев стартуют?
|
||||
docker compose -f $C --profile check run --rm preflight # 2. все 205 сценариев стартуют?
|
||||
docker compose -f $C --profile check run --rm calibrate # 3. сколько MLUPS на этой машине?
|
||||
docker compose -f $C --profile check run --rm plan # 4. смета в часах по замеренному
|
||||
docker compose -f $C up -d # 5. кампания
|
||||
@@ -277,6 +277,38 @@ out/<id>/xt_*.csv x–t диаграммы (группа G)
|
||||
| **G** | 14 | старт, акустика, время жизни: x–t диаграммы, губки, предел по Re, прогон на 10⁷ шагов |
|
||||
| **H** | 9 | инварианты: симметрия, зеркальность, зависимость от числа Маха, расхождение f32 против f64 |
|
||||
| **I** | 9 | сверхмелкие сетки 4096×2048 по всем формам и композиции из трёх тел |
|
||||
| **J** | 65 | схемы стенки с пристеночной функцией: Bouzidi без неё против Bouzidi, Grad и HRR со Spalding — цилиндр Re=150 / 2000 (D=8…64) и 5000 (D=32, 64), NACA 0012 α=4° Re=2000 (хорда 24…96) и 10⁴ (48, 96); эталон каждой серии — Bouzidi при удвоенном разрешении |
|
||||
| **K** | 25 | внешний домен: грубые буферы до тела, после и по бокам, переход разрешения скачком ×4 / ×8 против ступеней ×2, длина тонкого следа, ширина промежуточного уровня, боковые губки — против однородной тонкой сетки на домене 30D+60D, ±30D |
|
||||
|
||||
### Группы J и K: как читать
|
||||
|
||||
Обе группы отвечают на практические вопросы, и обе меряются против эталона внутри группы, а не
|
||||
против литературы — для турбулентных Re и для домена её просто нет.
|
||||
|
||||
**J** — что точнее: Bouzidi (лучшая геометрия) или моментные схемы с пристеночной функцией
|
||||
Spalding. Bouzidi + Spalding добавлен, чтобы отделить вклад функции от вклада схемы. Сетка
|
||||
однородная, разрешение задано числом клеток на тело. Re выбраны с оглядкой на устойчивость:
|
||||
предел здесь задаёт не модель стенки, а τ − ½ = 3uD/Re — при τ − ½ < 6·10⁻⁴ счёт
|
||||
разваливается у всех схем одинаково (замерено: цилиндр Re = 10⁴ при D = 16 и 32), и такой
|
||||
прогон сравнивал бы не схемы, а предел устойчивости. Отсюда Re = 5000 вместо 10⁴ у цилиндра и
|
||||
хорды от 48 у профиля при Re = 10⁴. В каждой сводке — y⁺ первого узла: при y⁺ < 5 функция
|
||||
обязана почти совпасть с прилипанием, разница возможна только там, где первая клетка дальше.
|
||||
|
||||
**K** — как строить внешний домен (ключ `--levels` решателя). Тело везде одно: цилиндр с
|
||||
D = 32 на самом тонком уровне. База — L0 в 4 раза грубее тела, домен 20D до тела, 40D после,
|
||||
±20D, две ступени ×2. Серии меняют по одному параметру. Кроме St, ⟨Cd⟩ и rms Cl сравнивается
|
||||
`rho_probe_rms` — пульсация плотности в зонде, мера паразитной акустики от отражений на границах
|
||||
и стыках уровней, — и стоимость: конфигурации K обходятся в 1–2 % стоимости эталона.
|
||||
|
||||
Сводка обеих групп:
|
||||
|
||||
```sh
|
||||
python compare.py # Markdown-таблицы против эталонов, читает out/
|
||||
python compare.py --out ../all_out --csv jk.csv # по сведённым результатам пяти долей
|
||||
```
|
||||
|
||||
Не посчитанный прогон помечается «нет», развалившийся — «развал»; под каждой серией J — какая
|
||||
схема ближе всех к эталону на каждом разрешении (по |ΔCd| у цилиндра, по |ΔCl| у профиля).
|
||||
|
||||
## Стоимость и время
|
||||
|
||||
@@ -298,6 +330,12 @@ python gen_scenarios.py --scale 3 # все прогоны в полтора
|
||||
группы A идут на CPU — это указано в самих сценариях.
|
||||
- **Развалы ожидаемы** в группе G (предел по Re) и у прогона `A09_shear_n512_lbgk`: LBGK при
|
||||
Re=3·10⁴ обязан развалиться там, где KBC доживает — это и есть проверяемое утверждение.
|
||||
Длинные прогоны канала показали предел около τ − ½ ≈ 6·10⁻⁴ (τ − ½ = 3uD/Re), а не
|
||||
Re ≈ 5·10⁴ при D = 16, как оценивалось раньше по коротким прогонам: при этом пределе
|
||||
`G11_limit_re10000` (τ − ½ = 7.2·10⁻⁴) на грани, а `G12` и `G13` обязаны развалиться.
|
||||
- **`F09_refine3` до версии 1.4.0 считался на GPU с заниженным в 1.5 раза Cd** — сила третьего
|
||||
подшага терялась за краем буфера. Исправлено; результаты этого прогона из образов до 1.4.0
|
||||
не использовать.
|
||||
- **Гифки** пишутся на всю длительность прогона, в реальном времени, 10 кадр/с. Число кадров
|
||||
этим задано жёстко (у самого длинного прогона их 16 666), поэтому единственный рычаг —
|
||||
размер кадра. Замерено: 0.103 байта на пиксель после LZW; отсюда бюджет
|
||||
|
||||
@@ -0,0 +1,194 @@
|
||||
#!/usr/bin/env python3
|
||||
"""Сводные таблицы групп J и K против их эталонов.
|
||||
|
||||
python compare.py # Markdown в консоль, читает out/
|
||||
python compare.py --out ../out_all # другой каталог результатов (например, сведённый
|
||||
# из пяти долей)
|
||||
python compare.py --csv compare.csv # те же строки в CSV
|
||||
|
||||
Группа J — схемы стенки. Каждая серия (тело, Re) имеет эталон: Bouzidi без пристеночной
|
||||
функции при удвоенном разрешении. Для каждого прогона печатается отклонение St, ⟨Cd⟩,
|
||||
rms Cl (у профиля — ещё ⟨Cl⟩) от эталона и y⁺ первого узла; в конце серии — какая схема
|
||||
ближе к эталону на каждом разрешении.
|
||||
|
||||
Группа K — внешний домен. Эталон каждого Re — однородная тонкая сетка на самом большом
|
||||
домене (K…_ref). Для каждой конфигурации — отклонения тех же величин, пульсация плотности в
|
||||
зонде (паразитная акустика) и стоимость (обновлений узлов и фактические секунды).
|
||||
|
||||
Прогон без summary.json (не посчитан) или развалившийся показывается строкой с пометкой —
|
||||
дыры в таблице видны сразу, а не теряются.
|
||||
"""
|
||||
|
||||
import argparse
|
||||
import csv
|
||||
import json
|
||||
import os
|
||||
import re
|
||||
import sys
|
||||
|
||||
HERE = os.path.dirname(os.path.abspath(__file__))
|
||||
|
||||
|
||||
def load(out_dir, scen_path):
|
||||
runs = json.load(open(scen_path, encoding="utf-8"))["runs"]
|
||||
res = {}
|
||||
for r in runs:
|
||||
p = os.path.join(out_dir, r["id"], "summary.json")
|
||||
if os.path.exists(p):
|
||||
try:
|
||||
res[r["id"]] = json.load(open(p, encoding="utf-8"))
|
||||
except json.JSONDecodeError:
|
||||
res[r["id"]] = None
|
||||
return runs, res
|
||||
|
||||
|
||||
def dev(v, ref):
|
||||
"""Относительное отклонение в процентах; None, если сравнивать нечего."""
|
||||
if v is None or ref is None or ref == 0:
|
||||
return None
|
||||
return 100.0 * (v - ref) / abs(ref)
|
||||
|
||||
|
||||
def fmt(v, spec="{:+.2f}"):
|
||||
return "—" if v is None else spec.format(v)
|
||||
|
||||
|
||||
def status(s):
|
||||
if s is None:
|
||||
return "нет"
|
||||
return "развал" if s.get("blew_up") else "ok"
|
||||
|
||||
|
||||
# ── J ────────────────────────────────────────────────────────────────────────
|
||||
|
||||
J_RE = re.compile(r"^J\d+_(cyl|naca)_re(\d+)_(?:(?:d|c)(\d+)_(\w+)|ref_(?:d|c)(\d+))$")
|
||||
|
||||
|
||||
def table_j(runs, res):
|
||||
series = {}
|
||||
for r in runs:
|
||||
m = J_RE.match(r["id"])
|
||||
if not m:
|
||||
continue
|
||||
body, re_, n, scheme, nref = m.groups()
|
||||
key = (body, int(re_))
|
||||
s = series.setdefault(key, {"ref": None, "rows": []})
|
||||
if nref:
|
||||
s["ref"] = r["id"]
|
||||
else:
|
||||
s["rows"].append((int(n), scheme, r["id"]))
|
||||
|
||||
rows, lines = [], []
|
||||
for (body, re_), s in sorted(series.items()):
|
||||
ref = res.get(s["ref"]) if s["ref"] else None
|
||||
name = "цилиндр" if body == "cyl" else "NACA 0012, α=4°"
|
||||
unit = "D" if body == "cyl" else "хорда"
|
||||
lines.append(f"\n### J: {name}, Re = {re_}\n")
|
||||
lines.append(f"Эталон `{s['ref']}`: " + (
|
||||
"не посчитан" if ref is None else
|
||||
f"St {fmt(ref.get('strouhal'), '{:.4f}')}, Cd {fmt(ref.get('cd'), '{:.4f}')}, "
|
||||
f"Cl {fmt(ref.get('cl_mean'), '{:.4f}')}, rms Cl {fmt(ref.get('cl_rms'), '{:.4f}')}"))
|
||||
lines.append("")
|
||||
lines.append(f"| {unit} | схема | статус | ΔSt % | ΔCd % | ΔCl % | Δrms Cl % | y⁺ ср | y⁺ макс |")
|
||||
lines.append("|---|---|---|---|---|---|---|---|---|")
|
||||
best = {}
|
||||
for n, scheme, rid in sorted(s["rows"]):
|
||||
v = res.get(rid)
|
||||
ok = v is not None and not v.get("blew_up")
|
||||
g = (lambda k: v.get(k)) if ok else (lambda k: None)
|
||||
rf = (lambda k: ref.get(k)) if ref else (lambda k: None)
|
||||
d_st, d_cd = dev(g("strouhal"), rf("strouhal")), dev(g("cd"), rf("cd"))
|
||||
d_cl, d_clr = dev(g("cl_mean"), rf("cl_mean")), dev(g("cl_rms"), rf("cl_rms"))
|
||||
lines.append(
|
||||
f"| {n} | {scheme} | {status(v)} | {fmt(d_st)} | {fmt(d_cd)} | "
|
||||
f"{fmt(d_cl if body == 'naca' else None)} | {fmt(d_clr)} | "
|
||||
f"{fmt(g('yplus_mean'), '{:.2f}')} | {fmt(g('yplus_max'), '{:.1f}')} |")
|
||||
rows.append({"group": "J", "body": body, "re": re_, "n": n, "scheme": scheme,
|
||||
"id": rid, "status": status(v), "d_st": d_st, "d_cd": d_cd,
|
||||
"d_cl": d_cl, "d_cl_rms": d_clr, "yplus_mean": g("yplus_mean"),
|
||||
"yplus_max": g("yplus_max")})
|
||||
# для профиля главное — подъёмная сила, для цилиндра — сопротивление
|
||||
score = d_cl if body == "naca" else d_cd
|
||||
if score is not None:
|
||||
b = best.get(n)
|
||||
if b is None or abs(score) < abs(b[1]):
|
||||
best[n] = (scheme, score)
|
||||
if best:
|
||||
what = "|ΔCl|" if body == "naca" else "|ΔCd|"
|
||||
lines.append("")
|
||||
lines.append(f"Ближе всех к эталону по {what}: " + "; ".join(
|
||||
f"{unit} {n} — **{sc}** ({sv:+.2f} %)" for n, (sc, sv) in sorted(best.items())))
|
||||
return rows, lines
|
||||
|
||||
|
||||
# ── K ────────────────────────────────────────────────────────────────────────
|
||||
|
||||
K_RE = re.compile(r"^K\d+_re(\d+)_(\w+)$")
|
||||
|
||||
|
||||
def table_k(runs, res):
|
||||
by_re = {}
|
||||
meta = {r["id"]: r for r in runs}
|
||||
for r in runs:
|
||||
m = K_RE.match(r["id"])
|
||||
if m:
|
||||
by_re.setdefault(int(m.group(1)), []).append((m.group(2), r["id"]))
|
||||
|
||||
rows, lines = [], []
|
||||
for re_, items in sorted(by_re.items()):
|
||||
ref_id = next((rid for tag, rid in items if tag == "ref"), None)
|
||||
ref = res.get(ref_id) if ref_id else None
|
||||
lines.append(f"\n### K: внешний домен, Re = {re_}\n")
|
||||
lines.append(f"Эталон `{ref_id}`: " + (
|
||||
"не посчитан" if ref is None else
|
||||
f"St {fmt(ref.get('strouhal'), '{:.4f}')}, Cd {fmt(ref.get('cd'), '{:.4f}')}, "
|
||||
f"rms Cl {fmt(ref.get('cl_rms'), '{:.4f}')}, ρ′ {fmt(ref.get('rho_probe_rms'), '{:.2e}')}"))
|
||||
lines.append("")
|
||||
lines.append("| конфигурация | статус | ΔSt % | ΔCd % | Δrms Cl % | ρ′ зонда | "
|
||||
"стоимость, отн. эталона | время, с |")
|
||||
lines.append("|---|---|---|---|---|---|---|---|")
|
||||
ref_cost = meta[ref_id]["cost"] if ref_id else None
|
||||
for tag, rid in items:
|
||||
if tag == "ref":
|
||||
continue
|
||||
v = res.get(rid)
|
||||
ok = v is not None and not v.get("blew_up")
|
||||
g = (lambda k: v.get(k)) if ok else (lambda k: None)
|
||||
rf = (lambda k: ref.get(k)) if ref else (lambda k: None)
|
||||
d_st, d_cd = dev(g("strouhal"), rf("strouhal")), dev(g("cd"), rf("cd"))
|
||||
d_clr = dev(g("cl_rms"), rf("cl_rms"))
|
||||
cost = meta[rid]["cost"] / ref_cost if ref_cost else None
|
||||
lines.append(
|
||||
f"| {tag} | {status(v)} | {fmt(d_st)} | {fmt(d_cd)} | {fmt(d_clr)} | "
|
||||
f"{fmt(g('rho_probe_rms'), '{:.2e}')} | {fmt(cost, '{:.3f}')} | "
|
||||
f"{fmt(g('wall_time_s'), '{:.0f}')} |")
|
||||
rows.append({"group": "K", "re": re_, "config": tag, "id": rid,
|
||||
"status": status(v), "d_st": d_st, "d_cd": d_cd, "d_cl_rms": d_clr,
|
||||
"rho_probe_rms": g("rho_probe_rms"), "cost_rel": cost,
|
||||
"wall_time_s": g("wall_time_s")})
|
||||
return rows, lines
|
||||
|
||||
|
||||
def main():
|
||||
ap = argparse.ArgumentParser(description="Сводка групп J и K против эталонов")
|
||||
ap.add_argument("--out", default=os.path.join(HERE, "out"), help="каталог результатов")
|
||||
ap.add_argument("--scenarios", default=os.path.join(HERE, "scenarios.json"))
|
||||
ap.add_argument("--csv", help="записать строки таблиц в CSV")
|
||||
args = ap.parse_args()
|
||||
|
||||
runs, res = load(args.out, args.scenarios)
|
||||
rj, lj = table_j(runs, res)
|
||||
rk, lk = table_k(runs, res)
|
||||
print("\n".join(["# Сравнение схем стенки (J) и конфигураций домена (K)"] + lj + lk))
|
||||
if args.csv:
|
||||
keys = sorted({k for row in rj + rk for k in row})
|
||||
with open(args.csv, "w", newline="", encoding="utf-8") as f:
|
||||
w = csv.DictWriter(f, fieldnames=keys)
|
||||
w.writeheader()
|
||||
for row in rj + rk:
|
||||
w.writerow(row)
|
||||
print(f"\nCSV: {args.csv}", file=sys.stderr)
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
main()
|
||||
@@ -485,6 +485,207 @@ def group_i(scale):
|
||||
nx * ny, st, gif="fine_multi.gif", nx=nx)
|
||||
|
||||
|
||||
# ══════════════════════════════════════════════════════════════════════════════
|
||||
# J. Схемы стенки с пристеночной функцией
|
||||
# ══════════════════════════════════════════════════════════════════════════════
|
||||
# Главный практический вопрос: что точнее — Bouzidi (лучшая геометрия, без пристеночной
|
||||
# функции) или моментные схемы с пристеночной функцией Spalding, и как это меняется с Re и
|
||||
# разрешением. Bouzidi + Spalding добавлен, чтобы отделить вклад функции от вклада схемы.
|
||||
# Сетка однородная (--refine 1): разрешение задаётся прямо числом клеток на тело, без
|
||||
# вклада связки уровней. Эталон каждой серии — Bouzidi без функции при удвоенном разрешении.
|
||||
J_SCHEMES = [
|
||||
("bouzidi", "none"),
|
||||
("bouzidi", "spalding"),
|
||||
("grad", "spalding"),
|
||||
("hrr", "spalding"),
|
||||
]
|
||||
|
||||
|
||||
def jid():
|
||||
return f"J{len([r for r in RUNS if r['group'] == 'J']) + 1:02d}"
|
||||
|
||||
|
||||
def wf_tag(wall, wf):
|
||||
return wall if wf == "none" else f"{wall}_wf"
|
||||
|
||||
|
||||
# Разрешения по Re. Предел устойчивости здесь задаёт не модель стенки, а τ − ½ = 3·u·D/Re на
|
||||
# грубейшем уровне: замерено на цилиндре 20D×12D (обычный Bouzidi, длинные прогоны), что при
|
||||
# τ − ½ = 2.4·10⁻⁴ и 4.8·10⁻⁴ (Re = 10⁴, D = 16 и 32) счёт разваливается — на 13 и 27 тысячах
|
||||
# шагов, у HRR со Spalding так же, — а при 6·10⁻⁴ и выше (Re = 2000 при D = 8, Re = 5000 при
|
||||
# D = 32, Re = 10⁴ при D = 64) держится. Поэтому большой Re цилиндра — 5000, а разрешения
|
||||
# выбраны так, чтобы τ − ½ ≥ 6·10⁻⁴. Прогоны ниже предела схемы бы не сравнивали, а лишь
|
||||
# повторяли предел устойчивости, который меряет группа G.
|
||||
J_CYL_D = {150: (8, 16, 32, 64), 2000: (8, 16, 32, 64), 5000: (32, 64)}
|
||||
J_FOIL_C = {2000: (24, 48, 96), 10000: (48, 96)}
|
||||
|
||||
|
||||
def group_j(scale):
|
||||
# цилиндр: домен 20D × 12D, как сходимость по D в группе B
|
||||
for re in (150, 2000, 5000):
|
||||
for d in J_CYL_D[re]:
|
||||
nx, ny = 20 * d, 12 * d
|
||||
st = int(conv_steps(d, 300) * scale)
|
||||
for wall, wf in J_SCHEMES:
|
||||
add(f"{jid()}_cyl_re{re}_d{d}_{wf_tag(wall, wf)}", "J",
|
||||
f"Цилиндр Re={re}, D={d}, стенка {wall}, пристеночная функция {wf}",
|
||||
"отклонение St, ⟨Cd⟩, rms Cl от эталона серии (Bouzidi при D=128); y⁺ в сводке",
|
||||
cyl_args(d, re, nx, ny, wall=wall, wall_function=wf), nx * ny, st)
|
||||
d = 128
|
||||
nx, ny = 20 * d, 12 * d
|
||||
add(f"{jid()}_cyl_re{re}_ref_d{d}", "J",
|
||||
f"Цилиндр Re={re}, ЭТАЛОН серии: Bouzidi без функции, D={d}",
|
||||
"сеточно сошедшееся значение, к которому меряются схемы при D=8…64",
|
||||
cyl_args(d, re, nx, ny, wall="bouzidi", wall_function="none"), nx * ny,
|
||||
int(conv_steps(d, 300) * scale))
|
||||
|
||||
# профиль: NACA 0012 под углом 4°, домен 12c × 6c
|
||||
for re in (2000, 10000):
|
||||
for c in J_FOIL_C[re]:
|
||||
nx, ny = 12 * c, 6 * c
|
||||
st = int(conv_steps(c, 200) * scale)
|
||||
for wall, wf in J_SCHEMES:
|
||||
add(f"{jid()}_naca_re{re}_c{c}_{wf_tag(wall, wf)}", "J",
|
||||
f"NACA 0012, α=4°, Re={re}, хорда {c}, стенка {wall}, функция {wf}",
|
||||
"отклонение Cl, Cd от эталона серии (Bouzidi при хорде 192); y⁺ в сводке",
|
||||
foil_args(c, 4, re, nx, ny, wall=wall, wall_function=wf), nx * ny, st)
|
||||
c = 192
|
||||
nx, ny = 12 * c, 6 * c
|
||||
add(f"{jid()}_naca_re{re}_ref_c{c}", "J",
|
||||
f"NACA 0012, α=4°, Re={re}, ЭТАЛОН серии: Bouzidi без функции, хорда {c}",
|
||||
"сеточно сошедшееся значение для серии по хорде",
|
||||
foil_args(c, 4, re, nx, ny, wall="bouzidi", wall_function="none"), nx * ny,
|
||||
int(conv_steps(c, 200) * scale))
|
||||
|
||||
|
||||
# ══════════════════════════════════════════════════════════════════════════════
|
||||
# K. Внешний домен: грубые буферы, переходы разрешения, губки
|
||||
# ══════════════════════════════════════════════════════════════════════════════
|
||||
# Тело одно и то же во всех прогонах — цилиндр с D = 32 клетки на САМОМ ТОНКОМ уровне.
|
||||
# Меняется только то, что вокруг: размеры домена до тела, после и по бокам; как сетка
|
||||
# грубеет к краю (одним скачком или ступенями ×2); сколько тонкой сетки остаётся в следе;
|
||||
# есть ли губки на боковых кромках. Эталон — однородная тонкая сетка на самом большом
|
||||
# домене. Метрики: St, ⟨Cd⟩, rms Cl против эталона, пульсация плотности в зонде
|
||||
# (rho_probe_rms — паразитная акустика) и стоимость.
|
||||
K_D = 32 # клеток на тело на самом тонком уровне
|
||||
K_CONV = 300 # конвективных времён D/U до умножения на scale
|
||||
|
||||
|
||||
def _rround(x):
|
||||
"""Округление как у Rust f64::round — половина от нуля, а не к чётному."""
|
||||
return math.floor(x + 0.5) if x >= 0 else -math.floor(-x + 0.5)
|
||||
|
||||
|
||||
def level_chain(levels, cx, cy, d, nx, ny):
|
||||
"""Тот же расчёт границ уровней, что `parse_levels` в main.rs. Возвращает обновлений
|
||||
узлов за шаг L0 по всем уровням и строку для --levels."""
|
||||
pnx, pny, scale, ox, oy = nx, ny, 1.0, 0.0, 0.0
|
||||
total, sub = nx * ny, 1
|
||||
for r, up, down, side in levels:
|
||||
def to(x, o):
|
||||
return _rround((x - o) * scale)
|
||||
ax = int(max(to(cx - up * d, ox), 2.0))
|
||||
bx = int(max(min(to(cx + down * d, ox), pnx - 2.0), 0.0))
|
||||
ay = int(max(to(cy - side * d, oy), 1.0))
|
||||
by = int(max(min(to(cy + side * d, oy), pny - 2.0), 0.0))
|
||||
assert ax + 1 < bx and ay + 1 < by, (levels, ax, bx, ay, by)
|
||||
ox += ax / scale
|
||||
oy += ay / scale
|
||||
scale *= r
|
||||
pnx, pny = r * (bx - ax) + 1, r * (by - ay) + 1
|
||||
sub *= r
|
||||
total += sub * pnx * pny
|
||||
spec = "; ".join(f"{r}:{up:g},{down:g},{side:g}" for r, up, down, side in levels)
|
||||
return total, spec
|
||||
|
||||
|
||||
def kid():
|
||||
return f"K{len([r for r in RUNS if r['group'] == 'K']) + 1:02d}"
|
||||
|
||||
|
||||
def k_run(tag, title, expect, re, scale, up, down, side, levels, sponge_side_d=0.0):
|
||||
"""Прогон группы K. `up`, `down`, `side` — размеры домена L0 в диаметрах; `levels` —
|
||||
цепочка [(r, up, down, side), …] от внешнего уровня к внутреннему, отступы в D.
|
||||
Пустая цепочка — однородная тонкая сетка."""
|
||||
s = 1
|
||||
for lv in levels:
|
||||
s *= lv[0]
|
||||
d0 = K_D / s # тело в клетках L0
|
||||
nx = int(round((up + down) * d0))
|
||||
ny = int(round(2 * side * d0))
|
||||
cx, cy = up * d0, ny / 2
|
||||
args = {"--shape": "cylinder", "--size": f"{d0:g}", "--nx": nx, "--ny": ny,
|
||||
"--body-x": f"{cx:g}", "--re": re, "--wall": "bouzidi"}
|
||||
if levels:
|
||||
nodes, spec = level_chain(levels, cx, cy, d0, nx, ny)
|
||||
args["--levels"] = spec
|
||||
else:
|
||||
nodes = nx * ny
|
||||
args["--refine"] = 1
|
||||
if sponge_side_d > 0:
|
||||
args["--sponge-side"] = int(round(sponge_side_d * d0))
|
||||
# шагов L0 на одно и то же физическое время: тело на L0 в s раз мельче
|
||||
st = int(conv_steps(d0, K_CONV) * scale)
|
||||
add(f"{kid()}_re{re}_{tag}", "K", title, expect, flatten(args), nodes, st)
|
||||
|
||||
|
||||
# Базовая конфигурация: L0 в 4 раза грубее тела, к телу две ступени ×2. Промежуточный
|
||||
# уровень — 4D вверх, 15D вниз, 4D в стороны; тонкий — 1.5D вверх, 10D следа, 2D в стороны.
|
||||
K_MID = (2, 4, 15, 4)
|
||||
|
||||
|
||||
def k_fine(wake=10):
|
||||
return (2, 1.5, wake, 2)
|
||||
|
||||
|
||||
def group_k(scale):
|
||||
for re in (150, 2000):
|
||||
k_run("ref", f"Re={re}, ЭТАЛОН: однородная тонкая сетка 30D+60D, ±30D",
|
||||
"значения, к которым меряются все конфигурации буферов", re, scale,
|
||||
30, 60, 30, [])
|
||||
k_run("base", f"Re={re}, база: домен 20D+40D, ±20D, ступени ×2×2",
|
||||
"отклонение от эталона при умеренных буферах", re, scale,
|
||||
20, 40, 20, [K_MID, k_fine()])
|
||||
# переход: одно и то же суммарное огрубление — скачком или ступенями
|
||||
k_run("jump4", f"Re={re}, огрубление ×4 одним скачком",
|
||||
"скачок разрешения против ступеней ×2×2 (база)", re, scale,
|
||||
20, 40, 20, [(4, 1.5, 10, 2)])
|
||||
k_run("jump8", f"Re={re}, огрубление ×8 одним скачком",
|
||||
"скачок ×8 против цепочки ×2×2×2", re, scale,
|
||||
20, 40, 20, [(8, 1.5, 10, 2)])
|
||||
k_run("chain8", f"Re={re}, огрубление ×8 цепочкой ×2×2×2",
|
||||
"плавный переход при сильном огрублении", re, scale,
|
||||
20, 40, 20, [(2, 8, 25, 8), K_MID, k_fine()])
|
||||
k_run("sponge5", f"Re={re}, боковые губки 5D",
|
||||
"гасят ли губки отражения от боковых стенок", re, scale,
|
||||
20, 40, 20, [K_MID, k_fine()], sponge_side_d=5)
|
||||
if re != 150:
|
||||
continue
|
||||
# остальные серии — только Re=150: там есть литература и чистая дорожка
|
||||
for side in (5, 10, 30):
|
||||
k_run(f"side{side}", f"Re={re}, боковой буфер ±{side}D",
|
||||
"с какого размера бока перестают влиять (блокировка + отражения)", re, scale,
|
||||
20, 40, side, [K_MID, k_fine()])
|
||||
for up in (5, 10, 30):
|
||||
k_run(f"up{up}", f"Re={re}, буфер до тела {up}D",
|
||||
"с какого отступа вход перестаёт влиять", re, scale,
|
||||
up, 40, 20, [(2, min(4, up - 1), 15, 4), k_fine()])
|
||||
for down in (15, 30, 60):
|
||||
k_run(f"down{down}", f"Re={re}, буфер после тела {down}D",
|
||||
"с какой длины выход перестаёт влиять", re, scale,
|
||||
20, down, 20, [(2, 4, min(15, down - 2), 4), k_fine(min(10, down - 4))])
|
||||
for wake in (5, 20):
|
||||
k_run(f"wake{wake}", f"Re={re}, тонкий след {wake}D за телом",
|
||||
"сколько тонкой сетки нужно в следе до первого огрубления", re, scale,
|
||||
20, 40, 20, [(2, 4, max(15, wake + 5), 4), k_fine(wake)])
|
||||
k_run("midwide", f"Re={re}, широкий промежуточный уровень (8D вверх, 25D вниз, 8D вбок)",
|
||||
"влияет ли ширина ступени перехода", re, scale,
|
||||
20, 40, 20, [(2, 8, 25, 8), k_fine()])
|
||||
k_run("sponge2", f"Re={re}, боковые губки 2D",
|
||||
"ширина губки: 2D против 5D", re, scale,
|
||||
20, 40, 20, [K_MID, k_fine()], sponge_side_d=2)
|
||||
|
||||
|
||||
def est_hours(r):
|
||||
"""Оценка времени прогона в часах для раскладки по долям (см. DZN_MLUPS)."""
|
||||
if r["backend"] == "cpu":
|
||||
@@ -526,7 +727,8 @@ def main():
|
||||
ap.add_argument("--out", default=os.path.join(os.path.dirname(__file__), "scenarios.json"))
|
||||
args = ap.parse_args()
|
||||
|
||||
for g in (group_a, group_b, group_c, group_d, group_e, group_f, group_g, group_h, group_i):
|
||||
for g in (group_a, group_b, group_c, group_d, group_e, group_f, group_g, group_h, group_i,
|
||||
group_j, group_k):
|
||||
g(args.scale)
|
||||
shard_hours = assign_shards(RUNS, args.shards)
|
||||
|
||||
|
||||
@@ -79,6 +79,26 @@ CASES = [
|
||||
"degenerate_frac": 5.0},
|
||||
"report": ["strouhal"],
|
||||
},
|
||||
# Пристеночная функция переписана в WGSL отдельно от math.rs (закон Spalding, Ньютон в
|
||||
# f32), поэтому сверяется отдельно. При Re=20 y⁺ < 1 и функция почти не меняет
|
||||
# результат — проверяется именно совпадение двух реализаций, а не физика.
|
||||
{
|
||||
"name": "цилиндр Re=20, Град + Spalding, 8000 шагов",
|
||||
"args": ["--shape", "cylinder", "--size", "16", "--nx", "320", "--ny", "192",
|
||||
"--re", "20", "--refine", "1", "--steps", "8000",
|
||||
"--wall", "grad", "--wall-function", "spalding"],
|
||||
"fields": {"cd": 1.0, "rho_mean": 0.05, "yplus_mean": 1.0},
|
||||
"report": ["strouhal"],
|
||||
},
|
||||
# Три уровня: рекурсивная связка уровней и сквозная нумерация подшагов силы на GPU.
|
||||
{
|
||||
"name": "цилиндр Re=20, три уровня ×2×2, 3000 шагов",
|
||||
"args": ["--shape", "cylinder", "--size", "8", "--nx", "240", "--ny", "160",
|
||||
"--re", "20", "--levels", "2:4,10,4; 2:1.5,4,1.5", "--wall", "bouzidi",
|
||||
"--steps", "3000"],
|
||||
"fields": {"cd": 1.0, "rho_mean": 0.05},
|
||||
"report": ["strouhal"],
|
||||
},
|
||||
]
|
||||
|
||||
|
||||
|
||||
File diff suppressed because it is too large
Load Diff
@@ -22,7 +22,7 @@ name: kbc2d
|
||||
|
||||
# ── общая часть всех сервисов ────────────────────────────────────────────────
|
||||
x-kbc2d: &kbc2d
|
||||
image: notbigghost/kbc2d:1.3.0
|
||||
image: notbigghost/kbc2d:1.4.0
|
||||
pull_policy: missing
|
||||
volumes:
|
||||
- ./out:/work/bench/out
|
||||
|
||||
@@ -15,7 +15,7 @@
|
||||
# * драйвер видеокарты под Windows — Vulkan внутри контейнера транслируется в D3D12 (dzn);
|
||||
# * Docker Desktop с бэкендом WSL2 (по умолчанию так и есть) либо Docker внутри WSL2.
|
||||
#
|
||||
# Доли посчитаны так, чтобы на ОДИНАКОВЫХ картах закончиться одновременно (≈24 ч каждая при
|
||||
# Доли посчитаны так, чтобы на ОДИНАКОВЫХ картах закончиться одновременно (≈27 ч каждая при
|
||||
# 1200 MLUPS, с поправкой на размер сетки — см. bench/gen_scenarios.py). На разных картах
|
||||
# закончатся в разное время. Результаты — в ./out рядом с этим файлом; собрать их с пяти
|
||||
# машин = скопировать все out/ в одну папку (каталоги прогонов не пересекаются, у каждой доли
|
||||
@@ -26,13 +26,13 @@
|
||||
name: kbc2d
|
||||
|
||||
x-kbc2d: &kbc2d
|
||||
image: notbigghost/kbc2d:1.3.0
|
||||
image: notbigghost/kbc2d:1.4.0
|
||||
# Если образа нет ни локально, ни в реестре — `docker compose -f docker-compose.shards.yml
|
||||
# build` соберёт его из каталога с этим файлом (нужны исходники).
|
||||
build:
|
||||
context: .
|
||||
args:
|
||||
VERSION: "1.3.0"
|
||||
VERSION: "1.4.0"
|
||||
pull_policy: missing
|
||||
devices:
|
||||
# сама видеокарта (WDDM); есть в любом WSL2, в том числе внутри Docker Desktop
|
||||
|
||||
@@ -32,14 +32,14 @@ name: kbc2d
|
||||
|
||||
# ── общая часть всех сервисов ────────────────────────────────────────────────
|
||||
x-kbc2d: &kbc2d
|
||||
image: notbigghost/kbc2d:1.3.0
|
||||
image: notbigghost/kbc2d:1.4.0
|
||||
# Контекст сборки — каталог с этим файлом. Если образа нет ни локально, ни в реестре,
|
||||
# достаточно `docker compose -f docker-compose.wsl.yml build`: доступ к Docker Hub
|
||||
# для запуска не обязателен.
|
||||
build:
|
||||
context: .
|
||||
args:
|
||||
VERSION: "1.3.0"
|
||||
VERSION: "1.4.0"
|
||||
pull_policy: missing
|
||||
devices:
|
||||
- /dev/dxg:/dev/dxg
|
||||
|
||||
+386
-167
@@ -8,8 +8,10 @@
|
||||
|
||||
use rayon::prelude::*;
|
||||
|
||||
use crate::math::{self, Kbc, Link, LinkKind, Scene, WallModel, MAX_BODY_BUCKETS, R, CX, CY, OPP, Q};
|
||||
use crate::{Case, Collision, FieldKind, Spec, StepRec};
|
||||
use crate::math::{
|
||||
self, Kbc, Link, LinkKind, Scene, WallFunction, WallModel, MAX_BODY_BUCKETS, R, CX, CY, OPP, Q, W,
|
||||
};
|
||||
use crate::{Case, Collision, FieldKind, PatchSpec, Spec, StepRec};
|
||||
|
||||
// ─────────────────────────────────────────────────────────────────────────────
|
||||
// Стартовые поля
|
||||
@@ -186,6 +188,8 @@ pub struct Geom {
|
||||
pub links: Vec<Link>,
|
||||
/// Индекс по граничным узлам: условие Града ставится сразу на весь узел, а не полинково.
|
||||
pub wall_nodes: Vec<WallNode>,
|
||||
/// Геометрия пристеночной функции, параллельно `wall_nodes` (тот же индекс).
|
||||
pub wf_nodes: Vec<WfNode>,
|
||||
/// Плечи для момента: координаты узлов относительно центра тела.
|
||||
pub body_cx: R,
|
||||
pub body_cy: R,
|
||||
@@ -214,8 +218,9 @@ impl Geom {
|
||||
}
|
||||
let links = build_links(nx, ny, &solid, &phi, &owner);
|
||||
let wall_nodes = build_wall_nodes(&links);
|
||||
let wf_nodes = build_wf_nodes(nx, ny, &solid, &wall_nodes, scene);
|
||||
let (bcx, bcy) = scene.center();
|
||||
Geom { nx, ny, solid, links, wall_nodes, body_cx: bcx, body_cy: bcy }
|
||||
Geom { nx, ny, solid, links, wall_nodes, wf_nodes, body_cx: bcx, body_cy: bcy }
|
||||
}
|
||||
}
|
||||
|
||||
@@ -320,6 +325,82 @@ pub struct WallNode {
|
||||
pub count: u8,
|
||||
}
|
||||
|
||||
/// Геометрия пристеночной функции одного граничного узла (Malaspinas, Sagaut 2014):
|
||||
/// нормаль к стенке, расстояние от неё до узла и до точки отбора скорости, отстоящей от узла
|
||||
/// на одну клетку вдоль нормали, и билинейный стенсиль этой точки. Строится один раз; GPU
|
||||
/// загружает её как есть — топология, как и линки, существует в единственном экземпляре.
|
||||
#[derive(Clone, Copy, Debug)]
|
||||
pub struct WfNode {
|
||||
/// Единичная нормаль, направленная из тела в жидкость (∇φ).
|
||||
pub nx: R,
|
||||
pub ny: R,
|
||||
/// Расстояние от стенки до граничного узла.
|
||||
pub y_w: R,
|
||||
/// Расстояние от стенки до точки отбора.
|
||||
pub y_s: R,
|
||||
/// Узлы и веса билинейного стенсиля точки отбора. Твёрдые узлы стенсиля исключены, веса
|
||||
/// перенормированы; если жидкого веса почти нет, все веса нулевые — функция в этом
|
||||
/// узле не применяется, остаётся прилипание.
|
||||
pub st: [u32; 4],
|
||||
pub w: [R; 4],
|
||||
}
|
||||
|
||||
impl WfNode {
|
||||
pub fn valid(&self) -> bool {
|
||||
self.w.iter().any(|v| *v > 0.0)
|
||||
}
|
||||
}
|
||||
|
||||
fn build_wf_nodes(
|
||||
nx: usize,
|
||||
ny: usize,
|
||||
solid: &[bool],
|
||||
wall_nodes: &[WallNode],
|
||||
scene: &Scene,
|
||||
) -> Vec<WfNode> {
|
||||
// SDF тел точный, поэтому нормаль — его градиент центральной разностью
|
||||
let h = 0.25;
|
||||
wall_nodes
|
||||
.iter()
|
||||
.map(|wn| {
|
||||
let node = wn.node as usize;
|
||||
let (x, y) = ((node % nx) as R, (node / nx) as R);
|
||||
let gx = scene.sdf(x + h, y) - scene.sdf(x - h, y);
|
||||
let gy = scene.sdf(x, y + h) - scene.sdf(x, y - h);
|
||||
let gn = (gx * gx + gy * gy).sqrt();
|
||||
let (ex, ey) = if gn > 1e-12 { (gx / gn, gy / gn) } else { (0.0, 1.0) };
|
||||
let y_w = scene.sdf(x, y).max(1e-3);
|
||||
let (sx, sy) = (x + ex, y + ey);
|
||||
let y_s = scene.sdf(sx, sy);
|
||||
let mut out = WfNode { nx: ex, ny: ey, y_w, y_s, st: [0; 4], w: [0.0; 4] };
|
||||
let (x0, y0) = (sx.floor(), sy.floor());
|
||||
if y_s <= y_w || x0 < 0.0 || y0 < 0.0 || x0 + 1.0 >= nx as R || y0 + 1.0 >= ny as R {
|
||||
return out;
|
||||
}
|
||||
let (tx, ty) = (sx - x0, sy - y0);
|
||||
let (x0, y0) = (x0 as usize, y0 as usize);
|
||||
let idx = [y0 * nx + x0, y0 * nx + x0 + 1, (y0 + 1) * nx + x0, (y0 + 1) * nx + x0 + 1];
|
||||
let wts = [(1.0 - tx) * (1.0 - ty), tx * (1.0 - ty), (1.0 - tx) * ty, tx * ty];
|
||||
let mut sum = 0.0;
|
||||
for k in 0..4 {
|
||||
out.st[k] = idx[k] as u32;
|
||||
if !solid[idx[k]] {
|
||||
out.w[k] = wts[k];
|
||||
sum += wts[k];
|
||||
}
|
||||
}
|
||||
if sum < 0.5 {
|
||||
out.w = [0.0; 4];
|
||||
} else {
|
||||
for v in out.w.iter_mut() {
|
||||
*v /= sum;
|
||||
}
|
||||
}
|
||||
out
|
||||
})
|
||||
.collect()
|
||||
}
|
||||
|
||||
fn build_wall_nodes(links: &[Link]) -> Vec<WallNode> {
|
||||
let mut out: Vec<WallNode> = Vec::new();
|
||||
for (k, l) in links.iter().enumerate() {
|
||||
@@ -345,7 +426,7 @@ pub struct Level {
|
||||
pub post: Vec<[R; Q]>,
|
||||
/// Энтропийный стабилизатор поузлово (для диагностики и картинки).
|
||||
pub gamma: Vec<R>,
|
||||
/// β по столбцам x (губка делает его полем); длина nx.
|
||||
/// β по узлам (губки делают его полем); длина nx·ny.
|
||||
pub beta: Vec<R>,
|
||||
pub geom: Geom,
|
||||
}
|
||||
@@ -372,7 +453,6 @@ impl Level {
|
||||
/// (ГУ Bouzidi перезаписывает всё, что могло бы прийти из тела в жидкость), а счёт там
|
||||
/// только жжёт такты и способен родить NaN при экстремальных режимах.
|
||||
fn collide(&mut self, op: Collision, model: math::KbcModel) -> Stats {
|
||||
let nx = self.nx;
|
||||
let beta = &self.beta;
|
||||
let solid = &self.geom.solid;
|
||||
let post = &mut self.post;
|
||||
@@ -388,7 +468,7 @@ impl Level {
|
||||
*g = 2.0;
|
||||
return Stats::EMPTY;
|
||||
}
|
||||
let b = beta[node % nx];
|
||||
let b = beta[node];
|
||||
let k: Kbc = match op {
|
||||
Collision::Kbc => math::collide_node(p, b, model),
|
||||
Collision::Bgk => math::collide_node_bgk(p, b),
|
||||
@@ -426,6 +506,67 @@ impl Level {
|
||||
}
|
||||
}
|
||||
|
||||
/// Пристеночная функция граничного узла `wi` по текущему полю: скорость в точке отбора —
|
||||
/// билинейно по пост-столкновительному полю (то есть на момент t, как и у Града).
|
||||
fn wall_fn_at(&self, wi: usize) -> Option<math::WallFnOut> {
|
||||
let g = &self.geom.wf_nodes[wi];
|
||||
if !g.valid() {
|
||||
return None;
|
||||
}
|
||||
let (mut ux, mut uy) = (0.0, 0.0);
|
||||
for k in 0..4 {
|
||||
if g.w[k] > 0.0 {
|
||||
let (a, b) = self.u_at_t(g.st[k] as usize);
|
||||
ux += g.w[k] * a;
|
||||
uy += g.w[k] * b;
|
||||
}
|
||||
}
|
||||
let node = self.geom.wall_nodes[wi].node as usize;
|
||||
let nu = math::nu_of_beta(self.beta[node]);
|
||||
math::wall_function(ux, uy, g.nx, g.ny, g.y_w, g.y_s, nu)
|
||||
}
|
||||
|
||||
/// Bouzidi с пристеночной функцией: стенка получает касательную скорость скольжения,
|
||||
/// согласованную по напряжению (см. `math::wall_slip`), и входит в отскок стандартной
|
||||
/// поправкой подвижной стенки
|
||||
/// δ = 2W_iρ(c_ī·u_s)/c_s² (Lallemand, Luo 2003); для q ≥ ½ поправка делится на 2q.
|
||||
fn apply_bouzidi_wf(&mut self) {
|
||||
for wi in 0..self.geom.wall_nodes.len() {
|
||||
let wn = self.geom.wall_nodes[wi];
|
||||
let node = wn.node as usize;
|
||||
let (usx, usy) = match self.wall_fn_at(wi) {
|
||||
Some(w) => {
|
||||
let g = &self.geom.wf_nodes[wi];
|
||||
let nu = math::nu_of_beta(self.beta[node]);
|
||||
let s = math::wall_slip(&w, g.y_w, nu);
|
||||
(s * w.ex, s * w.ey)
|
||||
}
|
||||
None => (0.0, 0.0),
|
||||
};
|
||||
let rho = math::macros(&self.post[node]).0;
|
||||
let lo = wn.first as usize;
|
||||
for l in &self.geom.links[lo..lo + wn.count as usize] {
|
||||
let i = l.i as usize;
|
||||
let ib = l.ib as usize;
|
||||
let fi = self.post[node][i];
|
||||
// c_ī = −c_i
|
||||
let delta =
|
||||
-2.0 * W[i] * rho * (CX[i] as R * usx + CY[i] as R * usy) / math::CS2;
|
||||
let v = match l.kind {
|
||||
LinkKind::Near => {
|
||||
2.0 * l.q * fi + (1.0 - 2.0 * l.q) * self.post[l.far as usize][i] + delta
|
||||
}
|
||||
LinkKind::Far => {
|
||||
let h = 1.0 / (2.0 * l.q);
|
||||
h * fi + (1.0 - h) * self.post[node][ib] + h * delta
|
||||
}
|
||||
LinkKind::Simple => fi + delta,
|
||||
};
|
||||
self.f[node][ib] = v;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/// Интерполированный отскок Bouzidi по всем линкам тела.
|
||||
/// `f` — поле ПОСЛЕ переноса (его правим), `post` — ПОСЛЕ столкновения (до переноса).
|
||||
fn apply_bouzidi(&mut self) {
|
||||
@@ -459,8 +600,7 @@ impl Level {
|
||||
///
|
||||
/// после чего недостающие популяции собираются приближением Града (2.13). Скорости
|
||||
/// соседей и градиенты берутся с прошлого шага — см. `u_prev`.
|
||||
fn apply_moment_wall(&mut self, third_order: bool) {
|
||||
let nx = self.nx;
|
||||
fn apply_moment_wall(&mut self, third_order: bool, wf: WallFunction) {
|
||||
// стенка неподвижна; для подвижного тела сюда пойдёт её скорость на линке,
|
||||
// и добавится динамическая часть плотности (B 4)
|
||||
let (uwx, uwy) = (0.0, 0.0);
|
||||
@@ -497,7 +637,21 @@ impl Level {
|
||||
rho += if missing[i] { self.post[node][OPP[i]] } else { self.f[node][i] };
|
||||
}
|
||||
|
||||
let (dudx, dudy, dvdx, dvdy) = self.grad_u_at_t(node);
|
||||
let (mut dudx, mut dudy, mut dvdx, mut dvdy) = self.grad_u_at_t(node);
|
||||
// Пристеночная функция: касательная составляющая целевой скорости и нормальная
|
||||
// производная касательной скорости в тензоре давлений берутся из закона стенки,
|
||||
// нормальная составляющая скорости и плотность (B 3) — прежние.
|
||||
if wf.is_on() {
|
||||
if let Some(w) = self.wall_fn_at(wi) {
|
||||
let g = &self.geom.wf_nodes[wi];
|
||||
let un = ux * g.nx + uy * g.ny;
|
||||
ux = un * g.nx + w.u_node * w.ex;
|
||||
uy = un * g.ny + w.u_node * w.ey;
|
||||
let nu = math::nu_of_beta(self.beta[node]);
|
||||
(dudx, dudy, dvdx, dvdy) =
|
||||
math::wall_function_gradient((dudx, dudy, dvdx, dvdy), &w, g.nx, g.ny, nu);
|
||||
}
|
||||
}
|
||||
let g = math::moment_wall(
|
||||
rho,
|
||||
ux,
|
||||
@@ -506,7 +660,7 @@ impl Level {
|
||||
dudy,
|
||||
dvdx,
|
||||
dvdy,
|
||||
self.beta[node % nx],
|
||||
self.beta[node],
|
||||
third_order,
|
||||
);
|
||||
for i in 0..Q {
|
||||
@@ -523,7 +677,7 @@ impl Level {
|
||||
/// что было в f до переноса, то есть u(x, t) — именно то, что требует прил. B. Отдельное
|
||||
/// хранилище «поля предыдущего шага» при этом не нужно.
|
||||
#[inline]
|
||||
fn u_at_t(&self, node: usize) -> (R, R) {
|
||||
pub fn u_at_t(&self, node: usize) -> (R, R) {
|
||||
let (_, a, b) = math::macros(&self.post[node]);
|
||||
(a, b)
|
||||
}
|
||||
@@ -553,11 +707,12 @@ impl Level {
|
||||
}
|
||||
|
||||
/// Замкнуть недостающие популяции выбранной моделью стенки.
|
||||
fn apply_wall(&mut self, model: WallModel) {
|
||||
fn apply_wall(&mut self, model: WallModel, wf: WallFunction) {
|
||||
match model {
|
||||
WallModel::Bouzidi if wf.is_on() => self.apply_bouzidi_wf(),
|
||||
WallModel::Bouzidi => self.apply_bouzidi(),
|
||||
WallModel::Grad => self.apply_moment_wall(false),
|
||||
WallModel::Hrr => self.apply_moment_wall(true),
|
||||
WallModel::Grad => self.apply_moment_wall(false, wf),
|
||||
WallModel::Hrr => self.apply_moment_wall(true, wf),
|
||||
WallModel::Staircase => self.apply_staircase(),
|
||||
}
|
||||
}
|
||||
@@ -729,12 +884,7 @@ pub struct Ghost {
|
||||
/// Неравновесная часть при смене уровня масштабируется: f^neq ∝ τ·δt, поэтому
|
||||
/// коэффициент грубый→тонкий равен R01 = τ_f/(r·τ_c), обратно — 1/R01.
|
||||
pub struct Patch {
|
||||
pub ax: usize,
|
||||
pub bx: usize,
|
||||
pub ay: usize,
|
||||
pub by: usize,
|
||||
pub r: usize,
|
||||
pub nfx: usize,
|
||||
pub r01: R,
|
||||
ghosts: Vec<Ghost>,
|
||||
/// Пары (узел L0, узел L1) для рестрикции — только внутренние жидкие узлы перекрытия.
|
||||
@@ -753,9 +903,9 @@ impl Patch {
|
||||
&self.restrict
|
||||
}
|
||||
|
||||
pub fn new(spec: &Spec, coarse: &Geom, fine_solid: &[bool], r01: R) -> Patch {
|
||||
let (ax, bx, ay, by) = spec.patch.expect("патч запрошен без границ");
|
||||
let r = spec.refine;
|
||||
pub fn new(ps: &PatchSpec, coarse: &Geom, fine_solid: &[bool], r01: R) -> Patch {
|
||||
let (ax, bx, ay, by) = (ps.ax, ps.bx, ps.ay, ps.by);
|
||||
let r = ps.r;
|
||||
let nfx = r * (bx - ax) + 1;
|
||||
let nfy = r * (by - ay) + 1;
|
||||
let cnx = coarse.nx;
|
||||
@@ -806,7 +956,7 @@ impl Patch {
|
||||
}
|
||||
// рамка — ровно периметр тонкого поля: два ряда по nfx плюс два столбца без углов
|
||||
assert_eq!(ghosts.len(), 2 * nfx + 2 * (nfy - 2));
|
||||
Patch { ax, bx, ay, by, r, nfx, r01, ghosts, restrict }
|
||||
Patch { r, r01, ghosts, restrict }
|
||||
}
|
||||
|
||||
/// Ghost-значения из грубого поля: равновесие по интерполированным ρ, u плюс
|
||||
@@ -885,31 +1035,72 @@ impl Patch {
|
||||
|
||||
pub struct Sim {
|
||||
pub spec: Spec,
|
||||
pub l0: Level,
|
||||
pub l1: Option<Level>,
|
||||
pub patch: Option<Patch>,
|
||||
/// Состояние L0 до столкновения — «старый» край для временной интерполяции рамки.
|
||||
pre: Vec<[R; Q]>,
|
||||
gh_old: Vec<[R; Q]>,
|
||||
gh_new: Vec<[R; Q]>,
|
||||
/// Уровни от L0 (весь домен) к самому тонкому.
|
||||
pub levels: Vec<Level>,
|
||||
/// `patches[k]` связывает уровень k с уровнем k + 1.
|
||||
pub patches: Vec<Patch>,
|
||||
/// Состояние уровня k до столкновения — «старый» край для временной интерполяции рамки
|
||||
/// уровня k + 1. Заведено только у уровней, у которых есть потомок.
|
||||
pre: Vec<Vec<[R; Q]>>,
|
||||
gh_old: Vec<Vec<[R; Q]>>,
|
||||
gh_new: Vec<Vec<[R; Q]>>,
|
||||
fluid_count: R,
|
||||
/// Зонд следа: уровень и узел на нём — самый тонкий уровень, накрывающий точку.
|
||||
probe_level: usize,
|
||||
probe_node: usize,
|
||||
probe_on_fine: bool,
|
||||
pub step_index: u64,
|
||||
}
|
||||
|
||||
/// Сводка пристеночной функции по граничным узлам самого тонкого уровня.
|
||||
#[derive(Clone, Copy, Debug, Default)]
|
||||
pub struct YPlus {
|
||||
/// Граничных узлов, где функция определена (есть точка отбора и касательный поток).
|
||||
pub nodes: usize,
|
||||
pub mean: R,
|
||||
pub max: R,
|
||||
/// Доля узлов с y⁺ > 30 — там первая клетка лежит в логарифмической зоне.
|
||||
pub frac_log: R,
|
||||
pub u_tau_mean: R,
|
||||
}
|
||||
|
||||
/// y⁺ первого узла по полю скорости `u(node)` уровня с геометрией `geom` и вязкостью `nu`.
|
||||
/// Общая для обоих бэкендов: GPU скачивает поле и считает здесь же.
|
||||
pub fn yplus_stats(geom: &Geom, nu: R, u: impl Fn(usize) -> (R, R)) -> YPlus {
|
||||
let mut y = YPlus::default();
|
||||
for g in &geom.wf_nodes {
|
||||
if !g.valid() {
|
||||
continue;
|
||||
}
|
||||
let (mut ux, mut uy) = (0.0, 0.0);
|
||||
for k in 0..4 {
|
||||
if g.w[k] > 0.0 {
|
||||
let (a, b) = u(g.st[k] as usize);
|
||||
ux += g.w[k] * a;
|
||||
uy += g.w[k] * b;
|
||||
}
|
||||
}
|
||||
if let Some(w) = math::wall_function(ux, uy, g.nx, g.ny, g.y_w, g.y_s, nu) {
|
||||
y.nodes += 1;
|
||||
y.mean += w.y_plus;
|
||||
y.max = y.max.max(w.y_plus);
|
||||
y.u_tau_mean += w.u_tau;
|
||||
if w.y_plus > 30.0 {
|
||||
y.frac_log += 1.0;
|
||||
}
|
||||
}
|
||||
}
|
||||
if y.nodes > 0 {
|
||||
let n = y.nodes as R;
|
||||
y.mean /= n;
|
||||
y.u_tau_mean /= n;
|
||||
y.frac_log /= n;
|
||||
}
|
||||
y
|
||||
}
|
||||
|
||||
impl Sim {
|
||||
pub fn new(spec: Spec) -> Sim {
|
||||
let (nx, ny) = (spec.nx, spec.ny);
|
||||
let geom0 = Geom::build(nx, ny, &spec.scene);
|
||||
|
||||
let beta0 = math::beta_profile(
|
||||
nx,
|
||||
spec.beta0,
|
||||
spec.sponge_in,
|
||||
spec.sponge_len,
|
||||
spec.sponge_mult,
|
||||
);
|
||||
let lv = spec.levels();
|
||||
|
||||
// стартовое поле: либо сразу набегающий поток, либо покой
|
||||
let u0 = if spec.init_uniform {
|
||||
@@ -918,67 +1109,78 @@ impl Sim {
|
||||
} else {
|
||||
(0.0, 0.0)
|
||||
};
|
||||
let init0 = initial_field(&Init {
|
||||
case: spec.case,
|
||||
nx,
|
||||
ny,
|
||||
u0,
|
||||
beta: spec.beta0,
|
||||
scene: Some(&spec.scene),
|
||||
taper: if spec.init_uniform { spec.init_taper } else { 0.0 },
|
||||
});
|
||||
let l0 = Level::new(nx, ny, beta0, geom0, init0);
|
||||
let fluid_count = l0.geom.solid.iter().filter(|s| !**s).count() as R;
|
||||
|
||||
let (l1, patch) = if spec.refine > 1 {
|
||||
let (ax, bx, ay, by) = spec.patch.expect("refine > 1 требует патч");
|
||||
let r = spec.refine;
|
||||
let scene1 = spec.scene.refined(r as R, ax as R, ay as R);
|
||||
let nfx = r * (bx - ax) + 1;
|
||||
let nfy = r * (by - ay) + 1;
|
||||
let geom1 = Geom::build(nfx, nfy, &scene1);
|
||||
// τ_f = r(τ_c − ½) + ½ ⇒ одинаковая ν на обоих уровнях
|
||||
let tau0 = 1.0 / (2.0 * spec.beta0);
|
||||
let tau1 = r as R * (tau0 - 0.5) + 0.5;
|
||||
let beta1 = 1.0 / (2.0 * tau1);
|
||||
let r01 = tau1 / (r as R * tau0);
|
||||
let patch = Patch::new(&spec, &l0.geom, &geom1.solid, r01);
|
||||
let init1 = initial_field(&Init {
|
||||
let mut levels: Vec<Level> = Vec::with_capacity(lv.len());
|
||||
let mut patches: Vec<Patch> = Vec::with_capacity(spec.patches.len());
|
||||
for (k, g) in lv.iter().enumerate() {
|
||||
let scene = spec.level_scene(g);
|
||||
let geom = Geom::build(g.nx, g.ny, &scene);
|
||||
// губки живут только на L0; тонкие уровни однородны, β_k = 1/(2τ_k)
|
||||
let (beta, beta_init) = if k == 0 {
|
||||
(
|
||||
math::beta_field(
|
||||
g.nx,
|
||||
g.ny,
|
||||
spec.beta0,
|
||||
spec.sponge_in,
|
||||
spec.sponge_len,
|
||||
spec.sponge_side,
|
||||
spec.sponge_mult,
|
||||
),
|
||||
spec.beta0,
|
||||
)
|
||||
} else {
|
||||
let b = 1.0 / (2.0 * g.tau);
|
||||
(vec![b; g.nx * g.ny], b)
|
||||
};
|
||||
let init = initial_field(&Init {
|
||||
case: spec.case,
|
||||
nx: nfx,
|
||||
ny: nfy,
|
||||
nx: g.nx,
|
||||
ny: g.ny,
|
||||
u0,
|
||||
beta: 1.0 / (2.0 * tau1),
|
||||
scene: Some(&scene1),
|
||||
taper: if spec.init_uniform { spec.init_taper * r as R } else { 0.0 },
|
||||
beta: beta_init,
|
||||
scene: Some(&scene),
|
||||
taper: if spec.init_uniform { spec.init_taper * g.scale } else { 0.0 },
|
||||
});
|
||||
(Some(Level::new(nfx, nfy, vec![beta1; nfx], geom1, init1)), Some(patch))
|
||||
} else {
|
||||
(None, None)
|
||||
};
|
||||
|
||||
// зонд следа: берём с тонкой сетки, если точка внутри патча (меньше численного
|
||||
// размытия вихрей → чище спектр и St), иначе с грубой
|
||||
let (px, py) = spec.probe;
|
||||
let (probe_node, probe_on_fine) = match &patch {
|
||||
Some(p) if px >= p.ax && px <= p.bx && py >= p.ay && py <= p.by => {
|
||||
(((py - p.ay) * p.r) * p.nfx + (px - p.ax) * p.r, true)
|
||||
if k > 0 {
|
||||
// f^neq ∝ τ·δt ⇒ коэффициент грубый→тонкий R01 = τ_f/(r·τ_c)
|
||||
let ps = &spec.patches[k - 1];
|
||||
let r01 = g.tau / (ps.r as R * lv[k - 1].tau);
|
||||
patches.push(Patch::new(ps, &levels[k - 1].geom, &geom.solid, r01));
|
||||
}
|
||||
_ => (py * nx + px, false),
|
||||
};
|
||||
levels.push(Level::new(g.nx, g.ny, beta, geom, init));
|
||||
}
|
||||
let fluid_count = levels[0].geom.solid.iter().filter(|s| !**s).count() as R;
|
||||
|
||||
let n = nx * ny;
|
||||
// зонд следа: с самого тонкого уровня, накрывающего точку (меньше численного
|
||||
// размытия вихрей → чище спектр и St), иначе с грубого
|
||||
let (mut px, mut py) = spec.probe;
|
||||
let mut probe_level = 0;
|
||||
for (k, p) in spec.patches.iter().enumerate() {
|
||||
if px >= p.ax && px <= p.bx && py >= p.ay && py <= p.by {
|
||||
px = (px - p.ax) * p.r;
|
||||
py = (py - p.ay) * p.r;
|
||||
probe_level = k + 1;
|
||||
} else {
|
||||
break;
|
||||
}
|
||||
}
|
||||
let probe_node = py * lv[probe_level].nx + px;
|
||||
|
||||
let nlev = levels.len();
|
||||
let pre = (0..nlev)
|
||||
.map(|k| if k + 1 < nlev { vec![[0.0; Q]; lv[k].nx * lv[k].ny] } else { Vec::new() })
|
||||
.collect();
|
||||
Sim {
|
||||
spec,
|
||||
l0,
|
||||
l1,
|
||||
patch,
|
||||
pre: vec![[0.0; Q]; n],
|
||||
gh_old: Vec::new(),
|
||||
gh_new: Vec::new(),
|
||||
levels,
|
||||
patches,
|
||||
pre,
|
||||
gh_old: vec![Vec::new(); nlev],
|
||||
gh_new: vec![Vec::new(); nlev],
|
||||
fluid_count,
|
||||
probe_level,
|
||||
probe_node,
|
||||
probe_on_fine,
|
||||
step_index: 0,
|
||||
}
|
||||
}
|
||||
@@ -1007,69 +1209,79 @@ impl Sim {
|
||||
(u * c, uy)
|
||||
}
|
||||
|
||||
/// Один шаг уровня k, а за ним — r_k подшагов его потомка (рекурсивно) с временной
|
||||
/// интерполяцией рамки и проекцией обратно. Сила снимается на самом тонком уровне на
|
||||
/// каждом его подшаге и копится в `fb`; `nf` — сколько подшагов в неё вошло.
|
||||
fn advance_level(
|
||||
&mut self,
|
||||
k: usize,
|
||||
inlet: (R, R),
|
||||
fb: &mut [[R; 3]; MAX_BODY_BUCKETS],
|
||||
nf: &mut usize,
|
||||
) -> Stats {
|
||||
let has_child = k < self.patches.len();
|
||||
let channel = self.spec.case == Case::Channel;
|
||||
let (wall, wf) = (self.spec.wall, self.spec.wall_fn);
|
||||
if has_child {
|
||||
self.pre[k].copy_from_slice(&self.levels[k].f);
|
||||
}
|
||||
let stats = self.levels[k].collide(self.spec.collision, self.spec.kbc_model);
|
||||
self.levels[k].stream();
|
||||
// Эталонные течения статей периодичны по обеим осям и ГУ не имеют вовсе: перенос
|
||||
// уже периодичен, поэтому достаточно ничего не накладывать. Стенки канала, вход и
|
||||
// выход есть только у L0 — тонкие уровни лежат строго внутри.
|
||||
if channel {
|
||||
self.levels[k].apply_wall(wall, wf);
|
||||
if k == 0 {
|
||||
self.levels[0].free_slip_walls();
|
||||
self.levels[0].channel_bc(inlet.0, inlet.1, 1.0, self.spec.outlet_extrapolate);
|
||||
}
|
||||
}
|
||||
if has_child {
|
||||
self.patches[k].ghost_values(&self.pre[k], &mut self.gh_old[k]);
|
||||
self.patches[k].ghost_values(&self.levels[k].f, &mut self.gh_new[k]);
|
||||
let r = self.patches[k].r;
|
||||
for s in 0..r {
|
||||
self.advance_level(k + 1, inlet, fb, nf);
|
||||
let w = (s + 1) as R / r as R;
|
||||
self.patches[k].fill(&mut self.levels[k + 1].f, &self.gh_old[k], &self.gh_new[k], w);
|
||||
}
|
||||
let (lo, hi) = self.levels.split_at_mut(k + 1);
|
||||
self.patches[k].restrict_to(&hi[0].f, &mut lo[k].f);
|
||||
} else {
|
||||
// силу снимаем на КАЖДОМ подшаге и усредняем — мгновенное значение на
|
||||
// последнем подшаге даёт лишний шум в рядах при том же среднем
|
||||
let g = self.levels[k].force();
|
||||
for b in 0..MAX_BODY_BUCKETS {
|
||||
for c in 0..3 {
|
||||
fb[b][c] += g[b][c];
|
||||
}
|
||||
}
|
||||
*nf += 1;
|
||||
}
|
||||
stats
|
||||
}
|
||||
|
||||
pub fn step(&mut self) -> StepRec {
|
||||
let t = self.step_index;
|
||||
let (ux_in, uy_in) = self.inlet(t);
|
||||
let sp_collision = self.spec.collision;
|
||||
let outlet_extrap = self.spec.outlet_extrapolate;
|
||||
let inlet = self.inlet(t);
|
||||
|
||||
// ── уровень 0 ──
|
||||
self.pre.copy_from_slice(&self.l0.f);
|
||||
let model = self.spec.kbc_model;
|
||||
let wall = self.spec.wall;
|
||||
let stats0 = self.l0.collide(sp_collision, model);
|
||||
self.l0.stream();
|
||||
// Эталонные течения статей периодичны по обеим осям и ГУ не имеют вовсе: перенос
|
||||
// уже периодичен, поэтому достаточно ничего не накладывать.
|
||||
if self.spec.case == Case::Channel {
|
||||
self.l0.apply_wall(wall);
|
||||
self.l0.free_slip_walls();
|
||||
self.l0.channel_bc(ux_in, uy_in, 1.0, outlet_extrap);
|
||||
}
|
||||
|
||||
// ── уровень 1: r подшагов с временной интерполяцией рамки ──
|
||||
let mut fb = [[0.0 as R; 3]; MAX_BODY_BUCKETS];
|
||||
if let (Some(l1), Some(p)) = (self.l1.as_mut(), self.patch.as_ref()) {
|
||||
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, model);
|
||||
l1.stream();
|
||||
if self.spec.case == Case::Channel {
|
||||
l1.apply_wall(wall);
|
||||
}
|
||||
// силу снимаем на КАЖДОМ подшаге и усредняем — мгновенное значение на
|
||||
// последнем подшаге даёт лишний шум в рядах при том же среднем
|
||||
let g = l1.force();
|
||||
for k in 0..MAX_BODY_BUCKETS {
|
||||
for c in 0..3 {
|
||||
fb[k][c] += g[k][c];
|
||||
}
|
||||
}
|
||||
let w = (s + 1) as R / p.r as R;
|
||||
p.fill(&mut l1.f, &self.gh_old, &self.gh_new, w);
|
||||
let mut nf = 0usize;
|
||||
let stats0 = self.advance_level(0, inlet, &mut fb, &mut nf);
|
||||
let inv = 1.0 / nf.max(1) as R;
|
||||
for b in fb.iter_mut() {
|
||||
for c in b.iter_mut() {
|
||||
*c *= inv;
|
||||
}
|
||||
let inv = 1.0 / p.r as R;
|
||||
for k in 0..MAX_BODY_BUCKETS {
|
||||
for c in 0..3 {
|
||||
fb[k][c] *= inv;
|
||||
}
|
||||
}
|
||||
p.restrict_to(&l1.f, &mut self.l0.f);
|
||||
} else {
|
||||
fb = self.l0.force();
|
||||
}
|
||||
|
||||
// ── диагностика ──
|
||||
let probe = if self.probe_on_fine {
|
||||
self.l1.as_ref().unwrap().macros_at(self.probe_node)
|
||||
} else {
|
||||
self.l0.macros_at(self.probe_node)
|
||||
};
|
||||
let probe = self.levels[self.probe_level].macros_at(self.probe_node);
|
||||
|
||||
let solid = &self.l0.geom.solid;
|
||||
let (rho_sum, max_u) = self
|
||||
.l0
|
||||
let l0 = &self.levels[0];
|
||||
let solid = &l0.geom.solid;
|
||||
let (rho_sum, max_u) = l0
|
||||
.f
|
||||
.par_iter()
|
||||
.enumerate()
|
||||
@@ -1082,10 +1294,10 @@ impl Sim {
|
||||
|
||||
self.step_index += 1;
|
||||
let (mut fx, mut fy, mut tz) = (0.0, 0.0, 0.0);
|
||||
for k in 0..MAX_BODY_BUCKETS {
|
||||
fx += fb[k][0];
|
||||
fy += fb[k][1];
|
||||
tz += fb[k][2];
|
||||
for b in &fb {
|
||||
fx += b[0];
|
||||
fy += b[1];
|
||||
tz += b[2];
|
||||
}
|
||||
StepRec {
|
||||
step: t,
|
||||
@@ -1094,6 +1306,7 @@ impl Sim {
|
||||
tz,
|
||||
body: fb,
|
||||
uy_probe: probe.2,
|
||||
rho_probe: probe.0,
|
||||
rho_mean: rho_sum / self.fluid_count,
|
||||
max_u,
|
||||
gamma_mean: stats0.gamma_mean(),
|
||||
@@ -1104,22 +1317,27 @@ impl Sim {
|
||||
}
|
||||
}
|
||||
|
||||
/// y⁺ первого узла на самом тонком уровне по текущему полю.
|
||||
pub fn yplus(&self) -> YPlus {
|
||||
let l = self.levels.last().unwrap();
|
||||
let tau = self.spec.levels().last().unwrap().tau;
|
||||
let nu = math::CS2 * (tau - 0.5);
|
||||
yplus_stats(&l.geom, nu, |n| l.u_at_t(n))
|
||||
}
|
||||
|
||||
/// Характерный размер тела в единицах того уровня, где снимается сила.
|
||||
pub fn force_ref_size(&self) -> R {
|
||||
match &self.patch {
|
||||
Some(p) => self.spec.scene.ref_size() * p.r as R,
|
||||
None => self.spec.scene.ref_size(),
|
||||
}
|
||||
self.spec.scene.ref_size() * self.spec.levels().last().unwrap().scale
|
||||
}
|
||||
|
||||
/// Поле для картинки: (значение, маска тела) на сетке L0.
|
||||
pub fn sample_field(&self, kind: FieldKind) -> (Vec<R>, &[bool]) {
|
||||
let (nx, ny) = (self.l0.nx, self.l0.ny);
|
||||
let (nx, ny) = (self.levels[0].nx, self.levels[0].ny);
|
||||
let mut out = vec![0.0; nx * ny];
|
||||
match kind {
|
||||
FieldKind::Speed => {
|
||||
for n in 0..nx * ny {
|
||||
let (_, ux, uy) = self.l0.macros_at(n);
|
||||
let (_, ux, uy) = self.levels[0].macros_at(n);
|
||||
out[n] = (ux * ux + uy * uy).sqrt();
|
||||
}
|
||||
}
|
||||
@@ -1128,7 +1346,7 @@ impl Sim {
|
||||
let mut ux = vec![0.0; nx * ny];
|
||||
let mut uy = vec![0.0; nx * ny];
|
||||
for n in 0..nx * ny {
|
||||
let (_, a, b) = self.l0.macros_at(n);
|
||||
let (_, a, b) = self.levels[0].macros_at(n);
|
||||
ux[n] = a;
|
||||
uy[n] = b;
|
||||
}
|
||||
@@ -1146,21 +1364,21 @@ impl Sim {
|
||||
}
|
||||
FieldKind::Density => {
|
||||
for n in 0..nx * ny {
|
||||
out[n] = self.l0.macros_at(n).0;
|
||||
out[n] = self.levels[0].macros_at(n).0;
|
||||
}
|
||||
}
|
||||
FieldKind::Gamma => out.copy_from_slice(&self.l0.gamma),
|
||||
FieldKind::Gamma => out.copy_from_slice(&self.levels[0].gamma),
|
||||
}
|
||||
(out, &self.l0.geom.solid)
|
||||
(out, &self.levels[0].geom.solid)
|
||||
}
|
||||
|
||||
/// Полное поле скорости уровня L0 — для метрик эталонных течений и радиуса влияния.
|
||||
pub fn sample_velocity(&self) -> (Vec<R>, Vec<R>) {
|
||||
let n = self.l0.nx * self.l0.ny;
|
||||
let n = self.levels[0].nx * self.levels[0].ny;
|
||||
let mut ux = Vec::with_capacity(n);
|
||||
let mut uy = Vec::with_capacity(n);
|
||||
for k in 0..n {
|
||||
let (_, a, b) = self.l0.macros_at(k);
|
||||
let (_, a, b) = self.levels[0].macros_at(k);
|
||||
ux.push(a);
|
||||
uy.push(b);
|
||||
}
|
||||
@@ -1171,12 +1389,12 @@ impl Sim {
|
||||
/// срезов складывается x–t диаграмма, по которой видно, бежит возмущение со скоростью
|
||||
/// звука или конвекции и есть ли стоячие узлы.
|
||||
pub fn sample_centerline(&self) -> (Vec<R>, Vec<R>) {
|
||||
let y = self.l0.ny / 2;
|
||||
let nx = self.l0.nx;
|
||||
let y = self.levels[0].ny / 2;
|
||||
let nx = self.levels[0].nx;
|
||||
let mut rho = Vec::with_capacity(nx);
|
||||
let mut ux = Vec::with_capacity(nx);
|
||||
for x in 0..nx {
|
||||
let (r, a, _) = self.l0.macros_at(y * nx + x);
|
||||
let (r, a, _) = self.levels[0].macros_at(y * nx + x);
|
||||
rho.push(r);
|
||||
ux.push(a);
|
||||
}
|
||||
@@ -1185,7 +1403,7 @@ impl Sim {
|
||||
|
||||
/// Есть ли в поле NaN/inf — признак развала счёта.
|
||||
pub fn is_finite(&self) -> bool {
|
||||
self.l0.f.par_iter().all(|c| c.iter().all(|v| v.is_finite()))
|
||||
self.levels[0].f.par_iter().all(|c| c.iter().all(|v| v.is_finite()))
|
||||
}
|
||||
}
|
||||
|
||||
@@ -1201,13 +1419,14 @@ mod tests {
|
||||
Level::new(
|
||||
nx,
|
||||
ny,
|
||||
vec![0.5; nx],
|
||||
vec![0.5; nx * ny],
|
||||
Geom {
|
||||
nx,
|
||||
ny,
|
||||
solid: vec![false; nx * ny],
|
||||
links: Vec::new(),
|
||||
wall_nodes: Vec::new(),
|
||||
wf_nodes: Vec::new(),
|
||||
body_cx: 0.0,
|
||||
body_cy: 0.0,
|
||||
},
|
||||
|
||||
@@ -382,7 +382,8 @@ pub struct GifWriter {
|
||||
h: usize,
|
||||
range: Range,
|
||||
hud: Hud,
|
||||
patch: Option<(usize, usize, usize, usize)>,
|
||||
/// Контуры патчей измельчения в клетках L0, от внешнего к внутреннему.
|
||||
patches: Vec<(usize, usize, usize, usize)>,
|
||||
dither: DelayDither,
|
||||
pub frames: u32,
|
||||
}
|
||||
@@ -399,7 +400,7 @@ impl GifWriter {
|
||||
range: Range,
|
||||
plan: GifPlan,
|
||||
hud: Hud,
|
||||
patch: Option<(usize, usize, usize, usize)>,
|
||||
patches: Vec<(usize, usize, usize, usize)>,
|
||||
) -> std::io::Result<GifWriter> {
|
||||
let scale = scale.max(1);
|
||||
let down = down.max(1);
|
||||
@@ -426,7 +427,7 @@ impl GifWriter {
|
||||
h,
|
||||
range,
|
||||
hud,
|
||||
patch,
|
||||
patches,
|
||||
dither: DelayDither::default(),
|
||||
frames: 0,
|
||||
})
|
||||
@@ -483,7 +484,7 @@ impl GifWriter {
|
||||
}
|
||||
}
|
||||
}
|
||||
if let Some((ax, bx, ay, by)) = self.patch {
|
||||
for &(ax, bx, ay, by) in &self.patches {
|
||||
let d = self.down;
|
||||
self.draw_patch_outline(&mut buf, ax / d, bx / d, ay / d, by / d);
|
||||
}
|
||||
|
||||
+495
-281
File diff suppressed because it is too large
Load Diff
@@ -97,6 +97,9 @@ pub struct StepRec {
|
||||
/// сваливаются в последнее, поэтому сумма по вёдрам всегда точна.
|
||||
pub body: [[R; 3]; MAX_BODY_BUCKETS],
|
||||
pub uy_probe: R,
|
||||
/// Плотность в том же зонде: её пульсации после установления — мера паразитной
|
||||
/// акустики (отражений от границ и от стыков уровней).
|
||||
pub rho_probe: R,
|
||||
pub rho_mean: R,
|
||||
pub max_u: R,
|
||||
pub gamma_mean: R,
|
||||
@@ -112,9 +115,9 @@ pub struct Spec {
|
||||
pub nx: usize,
|
||||
pub ny: usize,
|
||||
pub scene: Scene,
|
||||
pub refine: usize,
|
||||
/// Границы патча измельчения в координатах L0: (ax, bx, ay, by).
|
||||
pub patch: Option<(usize, usize, usize, usize)>,
|
||||
/// Цепочка вложенных патчей измельчения, от внешнего к внутреннему. Пусто — один
|
||||
/// уровень L0. Каждый патч задан в координатах СВОЕГО РОДИТЕЛЯ (для первого — L0).
|
||||
pub patches: Vec<PatchSpec>,
|
||||
pub units: Units,
|
||||
pub beta0: R,
|
||||
pub steps: u64,
|
||||
@@ -132,17 +135,116 @@ pub struct Spec {
|
||||
/// Длина губки после входа, столбцов.
|
||||
pub sponge_in: usize,
|
||||
pub sponge_mult: R,
|
||||
/// Боковые губки у верхней и нижней стенок канала, рядов L0 (0 — нет).
|
||||
pub sponge_side: usize,
|
||||
pub collision: Collision,
|
||||
/// Что входит в сдвиговую часть s (см. `math::KbcModel`).
|
||||
pub kbc_model: math::KbcModel,
|
||||
/// Чем замыкаются недостающие популяции на теле (см. `math::WallModel`).
|
||||
pub wall: math::WallModel,
|
||||
/// Пристеночная функция поверх модели стенки (см. `math::WallFunction`).
|
||||
pub wall_fn: math::WallFunction,
|
||||
/// Постановка задачи.
|
||||
pub case: Case,
|
||||
/// Узел зонда следа в координатах L0.
|
||||
pub probe: (usize, usize),
|
||||
}
|
||||
|
||||
/// Патч измельчения: прямоугольник узлов РОДИТЕЛЬСКОГО уровня [ax, bx] × [ay, by] и
|
||||
/// коэффициент r. Тонкий уровень имеет r·(bx − ax) + 1 узлов по x, его узел (0, 0) совпадает
|
||||
/// с узлом (ax, ay) родителя, и за один шаг родителя он делает r подшагов.
|
||||
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
|
||||
pub struct PatchSpec {
|
||||
pub r: usize,
|
||||
pub ax: usize,
|
||||
pub bx: usize,
|
||||
pub ay: usize,
|
||||
pub by: usize,
|
||||
}
|
||||
|
||||
/// Геометрия уровня, выведенная из цепочки патчей.
|
||||
#[derive(Clone, Copy, Debug)]
|
||||
pub struct LevelGeo {
|
||||
pub nx: usize,
|
||||
pub ny: usize,
|
||||
/// Клеток этого уровня на клетку L0.
|
||||
pub scale: R,
|
||||
/// Положение узла (0, 0) уровня в координатах L0.
|
||||
pub ox: R,
|
||||
pub oy: R,
|
||||
/// Время релаксации уровня: τ_k = r_k(τ_{k−1} − ½) + ½ — одна и та же ν на всех уровнях.
|
||||
pub tau: R,
|
||||
}
|
||||
|
||||
impl Spec {
|
||||
/// Все уровни, от L0 к самому тонкому.
|
||||
pub fn levels(&self) -> Vec<LevelGeo> {
|
||||
let tau0 = 1.0 / (2.0 * self.beta0);
|
||||
let mut out = vec![LevelGeo {
|
||||
nx: self.nx,
|
||||
ny: self.ny,
|
||||
scale: 1.0,
|
||||
ox: 0.0,
|
||||
oy: 0.0,
|
||||
tau: tau0,
|
||||
}];
|
||||
for p in &self.patches {
|
||||
let par = *out.last().unwrap();
|
||||
out.push(LevelGeo {
|
||||
nx: p.r * (p.bx - p.ax) + 1,
|
||||
ny: p.r * (p.by - p.ay) + 1,
|
||||
scale: par.scale * p.r as R,
|
||||
ox: par.ox + p.ax as R / par.scale,
|
||||
oy: par.oy + p.ay as R / par.scale,
|
||||
tau: p.r as R * (par.tau - 0.5) + 0.5,
|
||||
});
|
||||
}
|
||||
out
|
||||
}
|
||||
|
||||
/// Сцена в координатах уровня `g`.
|
||||
pub fn level_scene(&self, g: &LevelGeo) -> Scene {
|
||||
if g.scale == 1.0 && g.ox == 0.0 && g.oy == 0.0 {
|
||||
self.scene.clone()
|
||||
} else {
|
||||
self.scene.refined(g.scale, g.ox, g.oy)
|
||||
}
|
||||
}
|
||||
|
||||
/// Подшагов самого тонкого уровня на один шаг L0.
|
||||
pub fn substeps(&self) -> usize {
|
||||
self.patches.iter().map(|p| p.r).product()
|
||||
}
|
||||
|
||||
/// Обновлений узлов за один шаг L0 по всем уровням — мера стоимости.
|
||||
pub fn nodes_per_step(&self) -> f64 {
|
||||
let mut sub = 1usize;
|
||||
let mut total = 0usize;
|
||||
for (k, g) in self.levels().iter().enumerate() {
|
||||
if k > 0 {
|
||||
sub *= self.patches[k - 1].r;
|
||||
}
|
||||
total += sub * g.nx * g.ny;
|
||||
}
|
||||
total as f64
|
||||
}
|
||||
|
||||
/// Контуры патчей в клетках L0 — для рамок на гифке.
|
||||
pub fn patch_boxes_l0(&self) -> Vec<(usize, usize, usize, usize)> {
|
||||
let lv = self.levels();
|
||||
self.patches
|
||||
.iter()
|
||||
.enumerate()
|
||||
.map(|(k, p)| {
|
||||
let par = &lv[k];
|
||||
let x = |v: usize| (par.ox + v as R / par.scale).round() as usize;
|
||||
let y = |v: usize| (par.oy + v as R / par.scale).round() as usize;
|
||||
(x(p.ax), x(p.bx), y(p.ay), y(p.by))
|
||||
})
|
||||
.collect()
|
||||
}
|
||||
}
|
||||
|
||||
// ─────────────────────────────────────────────────────────────────────────────
|
||||
// Параметры запуска
|
||||
// ─────────────────────────────────────────────────────────────────────────────
|
||||
@@ -197,6 +299,13 @@ struct Cli {
|
||||
/// Границы патча в координатах L0: ax,bx,ay,by (по умолчанию строятся вокруг тела)
|
||||
#[arg(long, help_heading = "Сетка")]
|
||||
patch: Option<String>,
|
||||
/// ЦЕПОЧКА вложенных уровней вместо одного патча: "r:up,down,side; r:up,down,side; …" —
|
||||
/// от внешнего к внутреннему. r — во сколько раз уровень тоньше родителя; up, down, side —
|
||||
/// отступы его границ от центра тела в диаметрах тела (вверх по потоку, вниз, в стороны).
|
||||
/// Так строится внешний домен с грубыми буферами: L0 крупный и большой, к телу сетка
|
||||
/// мельчает ступенями. Перекрывает --refine и --patch.
|
||||
#[arg(long, help_heading = "Сетка")]
|
||||
levels: Option<String>,
|
||||
|
||||
// ── тело ──
|
||||
/// Форма обтекаемого тела
|
||||
@@ -282,6 +391,13 @@ struct Cli {
|
||||
#[arg(long, default_value = "hrr", value_parser = math::WallModel::ALL,
|
||||
help_heading = "Схема")]
|
||||
wall: String,
|
||||
/// Пристеночная функция поверх модели стенки: none — обычное прилипание; spalding — закон
|
||||
/// Spalding для неразрешённого пограничного слоя (для grad/hrr через целевую скорость и
|
||||
/// тензор давлений, для bouzidi через скорость скольжения стенки). В вязком подслое сама
|
||||
/// сводится к прилипанию. С --wall staircase не сочетается.
|
||||
#[arg(long, default_value = "none", value_parser = math::WallFunction::ALL,
|
||||
help_heading = "Схема")]
|
||||
wall_function: String,
|
||||
/// Состав сдвиговой части KBC (табл. I 2D-статьи): n1 — только девиатор {N, Π_xy}
|
||||
/// (KBC D), объёмная вязкость гуляет вместе с γ и может стать отрицательной; n2 —
|
||||
/// девиатор со следом {N, Π_xy, T} (KBC C), объёмная вязкость фиксирована ξ = ν.
|
||||
@@ -301,6 +417,10 @@ struct Cli {
|
||||
/// волны как жёсткий поршень; губка у входа гасит их до отражения.
|
||||
#[arg(long, default_value_t = 0, help_heading = "Схема")]
|
||||
sponge_in: usize,
|
||||
/// Боковые губки у верхней и нижней стенок канала, рядов L0 (0 — выключены). Гасят
|
||||
/// волны, бегущие поперёк потока, до отражения от зеркальных стенок.
|
||||
#[arg(long, default_value_t = 0, help_heading = "Схема")]
|
||||
sponge_side: usize,
|
||||
/// Во сколько раз губка поднимает вязкость
|
||||
#[arg(long, default_value_t = 30.0, help_heading = "Схема")]
|
||||
sponge_mult: R,
|
||||
@@ -455,6 +575,57 @@ fn parse_poly(v: &str) -> Result<Vec<[R; 2]>, String> {
|
||||
Ok(out)
|
||||
}
|
||||
|
||||
/// Разбор --levels: "r:up,down,side; …" от внешнего уровня к внутреннему. Отступы — от центра
|
||||
/// тела в его диаметрах; граница каждого уровня пересчитывается в координаты родителя и
|
||||
/// зажимается внутрь него (не ближе двух узлов к краю по x и одного по y — как у одиночного
|
||||
/// патча: рамка тонкого уровня не должна лечь на рамку или ряды стенок родителя).
|
||||
fn parse_levels(
|
||||
s: &str,
|
||||
cx: R,
|
||||
cy: R,
|
||||
d: R,
|
||||
nx: usize,
|
||||
ny: usize,
|
||||
) -> Result<Vec<PatchSpec>, String> {
|
||||
let mut out = Vec::new();
|
||||
// родитель: размеры, масштаб и начало в координатах L0
|
||||
let (mut pnx, mut pny, mut scale, mut ox, mut oy) = (nx, ny, 1.0 as R, 0.0 as R, 0.0 as R);
|
||||
for item in s.split(';').map(str::trim).filter(|z| !z.is_empty()) {
|
||||
let (rs, rest) = item
|
||||
.split_once(':')
|
||||
.ok_or(format!("--levels: в «{item}» ожидалось r:up,down,side"))?;
|
||||
let r: usize = rs.trim().parse().map_err(|e| format!("--levels, r в «{item}»: {e}"))?;
|
||||
if r < 2 {
|
||||
return Err(format!("--levels: r = {r} в «{item}», уровень обязан быть тоньше родителя"));
|
||||
}
|
||||
let v: Result<Vec<R>, _> = rest.split(',').map(|t| t.trim().parse::<R>()).collect();
|
||||
let v = v.map_err(|e| format!("--levels, отступы в «{item}»: {e}"))?;
|
||||
if v.len() != 3 || v.iter().any(|z| !(*z > 0.0)) {
|
||||
return Err(format!("--levels: в «{item}» нужны три положительных отступа up,down,side"));
|
||||
}
|
||||
let to = |x: R, o: R| ((x - o) * scale).round();
|
||||
let ax = to(cx - v[0] * d, ox).max(2.0) as usize;
|
||||
let bx = to(cx + v[1] * d, ox).min(pnx as R - 2.0).max(0.0) as usize;
|
||||
let ay = to(cy - v[2] * d, oy).max(1.0) as usize;
|
||||
let by = to(cy + v[2] * d, oy).min(pny as R - 2.0).max(0.0) as usize;
|
||||
if ax + 1 >= bx || ay + 1 >= by {
|
||||
return Err(format!(
|
||||
"--levels: уровень «{item}» не помещается в родителя {pnx}×{pny} (получилось x [{ax},{bx}], y [{ay},{by}])"
|
||||
));
|
||||
}
|
||||
out.push(PatchSpec { r, ax, bx, ay, by });
|
||||
ox += ax as R / scale;
|
||||
oy += ay as R / scale;
|
||||
scale *= r as R;
|
||||
pnx = r * (bx - ax) + 1;
|
||||
pny = r * (by - ay) + 1;
|
||||
}
|
||||
if out.is_empty() {
|
||||
return Err("--levels: список уровней пуст".into());
|
||||
}
|
||||
Ok(out)
|
||||
}
|
||||
|
||||
fn build_spec(cli: &Cli) -> Result<Spec, String> {
|
||||
let case = Case::from_str(&cli.case).ok_or("неизвестная постановка")?;
|
||||
if cli.nx < 16 || cli.ny < 16 {
|
||||
@@ -532,11 +703,22 @@ fn build_spec(cli: &Cli) -> Result<Spec, String> {
|
||||
));
|
||||
}
|
||||
|
||||
if case.is_periodic() && cli.refine > 1 {
|
||||
let wall = math::WallModel::from_str(&cli.wall).ok_or("неизвестная модель стенки")?;
|
||||
let wall_fn =
|
||||
math::WallFunction::from_str(&cli.wall_function).ok_or("неизвестная пристеночная функция")?;
|
||||
if wall_fn.is_on() && wall == math::WallModel::Staircase {
|
||||
return Err("--wall-function не сочетается с --wall staircase: ступенька не знает, где стенка, а пристеночной функции нужно расстояние до неё"
|
||||
.into());
|
||||
}
|
||||
if case.is_periodic() && (cli.refine > 1 || cli.levels.is_some()) {
|
||||
return Err("эталонные течения периодичны и патча измельчения не имеют: --refine 1".into());
|
||||
}
|
||||
// патч измельчения: по умолчанию охватывает тело и ближний след
|
||||
let patch = if cli.refine > 1 {
|
||||
// Цепочка патчей. --levels строит её по отступам от тела; без него — прежний одиночный
|
||||
// патч из --refine/--patch, который по умолчанию охватывает тело и ближний след.
|
||||
let patches: Vec<PatchSpec> = if let Some(lv) = &cli.levels {
|
||||
let (bcx, bcy) = scene.center();
|
||||
parse_levels(lv, bcx, bcy, scene.ref_size(), cli.nx, cli.ny)?
|
||||
} else if cli.refine > 1 {
|
||||
let d = cli.size;
|
||||
let (ax, bx, ay, by) = match &cli.patch {
|
||||
Some(s) => {
|
||||
@@ -560,8 +742,7 @@ fn build_spec(cli: &Cli) -> Result<Spec, String> {
|
||||
// на них ГУ работает по всему ряду, а рамка патча их бы перезаписала
|
||||
if !(1 <= ay && ay < by && by <= cli.ny - 2) {
|
||||
return Err(format!(
|
||||
"патч по y [{ay},{by}] обязан лежать внутри (0,{}) и не трогать ряды стенок; \
|
||||
при теле {d} ячеек минимальное ny ≈ {}",
|
||||
"патч по y [{ay},{by}] обязан лежать внутри (0,{}) и не трогать ряды стенок; при теле {d} ячеек минимальное ny ≈ {}",
|
||||
cli.ny - 1,
|
||||
(4.0 * d) as usize + 6
|
||||
));
|
||||
@@ -569,10 +750,26 @@ fn build_spec(cli: &Cli) -> Result<Spec, String> {
|
||||
if !(2 <= ax && ax < bx && bx <= cli.nx - 2) {
|
||||
return Err(format!("патч по x [{ax},{bx}] обязан лежать внутри (1,{})", cli.nx - 1));
|
||||
}
|
||||
Some((ax, bx, ay, by))
|
||||
vec![PatchSpec { r: cli.refine, ax, bx, ay, by }]
|
||||
} else {
|
||||
None
|
||||
Vec::new()
|
||||
};
|
||||
// Губки живут на L0 и не должны доставать до первого патча: внутри патча L0 перезаписывается
|
||||
// рестрикцией, и губка на нём просто бы не действовала.
|
||||
let patch = patches.first().map(|p| (p.ax, p.bx, p.ay, p.by));
|
||||
if cli.sponge_side > 0 {
|
||||
if 2 * cli.sponge_side + 4 >= cli.ny {
|
||||
return Err(format!("боковые губки ({} рядов) перекрывают весь канал", cli.sponge_side));
|
||||
}
|
||||
if let Some((_, _, ay, by)) = patch {
|
||||
if cli.sponge_side + 1 >= ay || by + cli.sponge_side + 1 >= cli.ny - 1 {
|
||||
return Err(format!(
|
||||
"боковые губки ({} рядов) достают до патча (ay={ay}, by={by})",
|
||||
cli.sponge_side
|
||||
));
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Губка перед выходом. Задача канала «вход по скорости + выход по давлению» акустически
|
||||
// есть четвертьволновая труба: вход отражает продольные волны как жёсткий поршень, выход —
|
||||
@@ -651,8 +848,7 @@ fn build_spec(cli: &Cli) -> Result<Spec, String> {
|
||||
nx: cli.nx,
|
||||
ny: cli.ny,
|
||||
scene,
|
||||
refine: 1,
|
||||
patch: None,
|
||||
patches: Vec::new(),
|
||||
units,
|
||||
beta0,
|
||||
steps,
|
||||
@@ -666,9 +862,11 @@ fn build_spec(cli: &Cli) -> Result<Spec, String> {
|
||||
sponge_len: 0,
|
||||
sponge_in: 0,
|
||||
sponge_mult: 1.0,
|
||||
sponge_side: 0,
|
||||
collision: if cli.collision == "bgk" { Collision::Bgk } else { Collision::Kbc },
|
||||
kbc_model: math::KbcModel::from_str(&cli.kbc_model).ok_or("неизвестная модель KBC")?,
|
||||
wall: math::WallModel::from_str(&cli.wall).ok_or("неизвестная модель стенки")?,
|
||||
wall,
|
||||
wall_fn,
|
||||
case,
|
||||
probe: (cli.nx / 2, cli.ny / 2),
|
||||
});
|
||||
@@ -678,8 +876,7 @@ fn build_spec(cli: &Cli) -> Result<Spec, String> {
|
||||
nx: cli.nx,
|
||||
ny: cli.ny,
|
||||
scene,
|
||||
refine: cli.refine,
|
||||
patch,
|
||||
patches,
|
||||
units,
|
||||
beta0,
|
||||
steps,
|
||||
@@ -693,10 +890,12 @@ fn build_spec(cli: &Cli) -> Result<Spec, String> {
|
||||
sponge_len,
|
||||
sponge_in: cli.sponge_in,
|
||||
sponge_mult: cli.sponge_mult,
|
||||
sponge_side: cli.sponge_side,
|
||||
case,
|
||||
collision: if cli.collision == "bgk" { Collision::Bgk } else { Collision::Kbc },
|
||||
kbc_model: math::KbcModel::from_str(&cli.kbc_model).ok_or("неизвестная модель KBC")?,
|
||||
wall: math::WallModel::from_str(&cli.wall).ok_or("неизвестная модель стенки")?,
|
||||
wall,
|
||||
wall_fn,
|
||||
probe,
|
||||
})
|
||||
}
|
||||
@@ -762,6 +961,13 @@ impl Backend {
|
||||
Backend::Gpu(s) => s.is_finite(),
|
||||
}
|
||||
}
|
||||
fn yplus(&mut self) -> cpu::YPlus {
|
||||
match self {
|
||||
Backend::Cpu(s) => s.yplus(),
|
||||
#[cfg(feature = "gpu")]
|
||||
Backend::Gpu(s) => s.yplus(),
|
||||
}
|
||||
}
|
||||
fn force_ref_size(&self) -> R {
|
||||
match self {
|
||||
Backend::Cpu(s) => s.force_ref_size(),
|
||||
@@ -864,7 +1070,7 @@ fn run(cli: Cli) -> Result<(), String> {
|
||||
range,
|
||||
plan,
|
||||
hud,
|
||||
spec.patch,
|
||||
spec.patch_boxes_l0(),
|
||||
)
|
||||
.map_err(|e| format!("не удалось создать {path}: {e}"))?,
|
||||
)
|
||||
@@ -890,13 +1096,7 @@ fn run(cli: Cli) -> Result<(), String> {
|
||||
let mut recs: Vec<StepRec> = Vec::with_capacity(spec.steps as usize + 1);
|
||||
let t_start = Instant::now();
|
||||
let mut blew_up = false;
|
||||
let nodes_per_step = (spec.nx * spec.ny
|
||||
+ spec
|
||||
.patch
|
||||
.map(|(ax, bx, ay, by)| {
|
||||
spec.refine * (spec.refine * (bx - ax) + 1) * (spec.refine * (by - ay) + 1)
|
||||
})
|
||||
.unwrap_or(0)) as f64;
|
||||
let nodes_per_step = spec.nodes_per_step();
|
||||
|
||||
let mut extras = Extras::default();
|
||||
let series_every = cli.series_every.max(1);
|
||||
@@ -930,6 +1130,15 @@ fn run(cli: Cli) -> Result<(), String> {
|
||||
if need_report {
|
||||
live_line(&spec, &rec, d_ref, t_start.elapsed().as_secs_f64(), nodes_per_step, full);
|
||||
}
|
||||
// y⁺ на теле: два десятка снимков во второй половине прогона. На GPU каждый снимок —
|
||||
// скачивание самого тонкого уровня, поэтому редко.
|
||||
if spec.case == Case::Channel && !spec.scene.bodies.is_empty() {
|
||||
let half = spec.steps / 2;
|
||||
let every = (spec.steps / 40).max(1);
|
||||
if t >= half && (t - half) % every == 0 {
|
||||
extras.yplus.push(back.yplus());
|
||||
}
|
||||
}
|
||||
if spec.case.is_periodic() && t % cli.case_every.max(1) == 0 {
|
||||
back.flush(&mut scratch);
|
||||
let (vx, vy) = back.sample_velocity();
|
||||
@@ -1057,6 +1266,8 @@ struct Extras {
|
||||
influence: Option<(R, R, R)>,
|
||||
/// Относительная ошибка поля u_x против точного решения Тейлора–Грина.
|
||||
tg_error: Option<R>,
|
||||
/// Снимки y⁺ первого узла на теле во второй половине прогона (канал).
|
||||
yplus: Vec<cpu::YPlus>,
|
||||
}
|
||||
|
||||
/// Интегральные метрики периодического течения: энергия, энстрофия, палинстрофия.
|
||||
@@ -1183,15 +1394,20 @@ fn print_header(spec: &Spec, cli: &Cli, plan: &gif::GifPlan, field: FieldKind) {
|
||||
println!("\n── Сетка ───────────────────────────────────────────────────────────────────");
|
||||
println!(" L0 {:>12} {} узлов", format!("{}×{}", spec.nx, spec.ny),
|
||||
spec.nx * spec.ny);
|
||||
match spec.patch {
|
||||
Some((ax, bx, ay, by)) => {
|
||||
let (nfx, nfy) = (spec.refine * (bx - ax) + 1, spec.refine * (by - ay) + 1);
|
||||
let tau1 = spec.refine as R * (tau - 0.5) + 0.5;
|
||||
println!(" L1 (×{}) {:>12} область x∈[{ax},{bx}] y∈[{ay},{by}], τ₁={:.4}",
|
||||
spec.refine, format!("{nfx}×{nfy}"), tau1);
|
||||
println!(" {:>12} {} подшагов на шаг L0", "", spec.refine);
|
||||
}
|
||||
None => println!(" L1 {:>12} (измельчение выключено)", "нет"),
|
||||
let _ = tau;
|
||||
let lv = spec.levels();
|
||||
if lv.len() == 1 {
|
||||
println!(" L1 {:>12} (измельчение выключено)", "нет");
|
||||
}
|
||||
for (k, p) in spec.patches.iter().enumerate() {
|
||||
let g = &lv[k + 1];
|
||||
let (ax, bx, ay, by) = (p.ax, p.bx, p.ay, p.by);
|
||||
println!(" L{} (×{}) {:>12} область L{k} x∈[{ax},{bx}] y∈[{ay},{by}], τ={:.4}, D={:.1}",
|
||||
k + 1, p.r, format!("{}×{}", g.nx, g.ny), g.tau, spec.scene.ref_size() * g.scale);
|
||||
}
|
||||
if !spec.patches.is_empty() {
|
||||
println!(" {:>12} {} подшагов тонкого уровня на шаг L0, {:.3e} обновлений узлов на шаг L0",
|
||||
"", spec.substeps(), spec.nodes_per_step());
|
||||
}
|
||||
println!(" оператор {:>12}",
|
||||
if spec.collision == Collision::Kbc {
|
||||
@@ -1574,12 +1790,12 @@ fn final_report(
|
||||
|
||||
fn write_csv(path: &str, recs: &[StepRec], spec: &Spec, d_ref: R) -> std::io::Result<()> {
|
||||
let mut f = std::io::BufWriter::new(std::fs::File::create(path)?);
|
||||
writeln!(f, "step,t_phys_s,cd,cl,cm,uy_probe,rho_mean,max_u,gamma_mean,gamma_min,gamma_max,degenerate_frac,xi_negative_frac")?;
|
||||
writeln!(f, "step,t_phys_s,cd,cl,cm,uy_probe,rho_mean,max_u,gamma_mean,gamma_min,gamma_max,degenerate_frac,xi_negative_frac,rho_probe")?;
|
||||
let s = spec.units.u_lat * spec.units.u_lat * d_ref;
|
||||
for r in recs {
|
||||
writeln!(
|
||||
f,
|
||||
"{},{:.9e},{:.6e},{:.6e},{:.6e},{:.6e},{:.9},{:.6e},{:.6},{:.6},{:.6},{:.6},{:.6}",
|
||||
"{},{:.9e},{:.6e},{:.6e},{:.6e},{:.6e},{:.9},{:.6e},{:.6},{:.6},{:.6},{:.6},{:.6},{:.9}",
|
||||
r.step,
|
||||
spec.units.time_of_step(r.step),
|
||||
2.0 * r.fx / s,
|
||||
@@ -1592,7 +1808,8 @@ fn write_csv(path: &str, recs: &[StepRec], spec: &Spec, d_ref: R) -> std::io::Re
|
||||
r.gamma_min,
|
||||
r.gamma_max,
|
||||
r.degenerate_frac,
|
||||
r.xi_negative_frac
|
||||
r.xi_negative_frac,
|
||||
r.rho_probe
|
||||
)?;
|
||||
}
|
||||
Ok(())
|
||||
@@ -1623,7 +1840,19 @@ fn write_summary(
|
||||
writeln!(f, " \"backend\": \"{}\",", cli.backend)?;
|
||||
writeln!(f, " \"wall\": \"{}\", \"kbc_model\": \"{}\", \"collision\": \"{}\",",
|
||||
cli.wall, cli.kbc_model, cli.collision)?;
|
||||
writeln!(f, " \"nx\": {}, \"ny\": {}, \"refine\": {},", spec.nx, spec.ny, spec.refine)?;
|
||||
writeln!(f, " \"nx\": {}, \"ny\": {}, \"refine\": {},", spec.nx, spec.ny, spec.substeps())?;
|
||||
let lv = spec.levels();
|
||||
let levels: Vec<String> = spec
|
||||
.patches
|
||||
.iter()
|
||||
.zip(lv.iter().skip(1))
|
||||
.map(|(p, g)| {
|
||||
format!("{{\"r\": {}, \"nx\": {}, \"ny\": {}, \"box\": [{}, {}, {}, {}]}}",
|
||||
p.r, g.nx, g.ny, p.ax, p.bx, p.ay, p.by)
|
||||
})
|
||||
.collect();
|
||||
writeln!(f, " \"levels\": [{}], \"nodes_per_step\": {}, \"wall_function\": \"{}\", \"sponge_side\": {},",
|
||||
levels.join(", "), spec.nodes_per_step(), spec.wall_fn.name(), spec.sponge_side)?;
|
||||
writeln!(f, " \"bodies\": {}, \"body_size_cells\": {},",
|
||||
spec.scene.bodies.len(), q(spec.scene.ref_size()))?;
|
||||
writeln!(f, " \"re\": {}, \"u_lat\": {}, \"mach\": {}, \"tau\": {},",
|
||||
@@ -1655,10 +1884,11 @@ fn write_summary(
|
||||
let rho: Vec<R> = recs[h..].iter().map(|r| r.rho_mean).collect();
|
||||
let (st, _, _) = math::strouhal(&uy, u.d_lat, u.u_lat, 8192);
|
||||
let (cdm, _) = math::mean_std(&cd);
|
||||
let (_, clr) = math::mean_std(&cl);
|
||||
let (clm, clr) = math::mean_std(&cl);
|
||||
let (cmm, _) = math::mean_std(&cm);
|
||||
let (rm, _) = math::mean_std(&rho);
|
||||
let bl = spec.scene.frontal_extent() / spec.ny as R;
|
||||
writeln!(f, " \"cl_mean\": {},", q(clm))?;
|
||||
writeln!(f, " \"strouhal\": {}, \"cd\": {}, \"cl_rms\": {}, \"cm\": {},",
|
||||
q(st), q(cdm), q(clr), q(cmm))?;
|
||||
writeln!(f, " \"blockage\": {}, \"strouhal_corr\": {}, \"cd_corr\": {}, \"cl_rms_corr\": {},",
|
||||
@@ -1690,6 +1920,18 @@ fn write_summary(
|
||||
math::mean_std(&recs[h..].iter().map(|r| r.degenerate_frac).collect::<Vec<_>>());
|
||||
let (xn, _) =
|
||||
math::mean_std(&recs[h..].iter().map(|r| r.xi_negative_frac).collect::<Vec<_>>());
|
||||
// Пульсация плотности в зонде во второй половине: мера паразитной акустики —
|
||||
// волн, отражённых от границ домена и от стыков уровней.
|
||||
let (_, rp) = math::mean_std(&recs[h..].iter().map(|r| r.rho_probe).collect::<Vec<_>>());
|
||||
writeln!(f, " \"rho_probe_rms\": {},", q(rp))?;
|
||||
if !extras.yplus.is_empty() {
|
||||
let k = extras.yplus.len() as R;
|
||||
let avg = |g: &dyn Fn(&cpu::YPlus) -> R| extras.yplus.iter().map(g).sum::<R>() / k;
|
||||
let ymax = extras.yplus.iter().map(|y| y.max).fold(0.0, R::max);
|
||||
writeln!(f, " \"yplus_mean\": {}, \"yplus_max\": {}, \"yplus_frac_log\": {}, \"u_tau_mean\": {}, \"yplus_nodes\": {},",
|
||||
q(avg(&|y| y.mean)), q(ymax), q(avg(&|y| y.frac_log)),
|
||||
q(avg(&|y| y.u_tau_mean)), extras.yplus[0].nodes)?;
|
||||
}
|
||||
writeln!(f, " \"gamma_mean\": {}, \"degenerate_frac\": {}, \"xi_negative_frac\": {}",
|
||||
q(gm), q(dg), q(xn))?;
|
||||
} else {
|
||||
@@ -1698,3 +1940,100 @@ fn write_summary(
|
||||
writeln!(f, "}}")?;
|
||||
Ok(())
|
||||
}
|
||||
|
||||
// ─────────────────────────────────────────────────────────────────────────────
|
||||
// Тесты постановки и многоуровневой сетки
|
||||
// ─────────────────────────────────────────────────────────────────────────────
|
||||
|
||||
#[cfg(test)]
|
||||
mod tests {
|
||||
use super::*;
|
||||
|
||||
fn spec_of(args: &[&str]) -> Spec {
|
||||
let mut v = vec!["kbc2d"];
|
||||
v.extend_from_slice(args);
|
||||
build_spec(&Cli::try_parse_from(v).expect("ключи")).expect("постановка")
|
||||
}
|
||||
|
||||
/// Цепочка уровней: каждый следующий в r раз тоньше, лежит внутри родителя, а τ растёт
|
||||
/// так, что ν одинакова на всех уровнях: ν = c_s²(τ − ½)/scale в единицах L0.
|
||||
#[test]
|
||||
fn levels_chain_is_consistent() {
|
||||
let s = spec_of(&[
|
||||
"--shape", "cylinder", "--size", "8", "--nx", "240", "--ny", "160", "--re", "20",
|
||||
"--levels", "2:4,10,4; 2:1.5,4,1.5",
|
||||
]);
|
||||
let lv = s.levels();
|
||||
assert_eq!(lv.len(), 3);
|
||||
assert_eq!(s.substeps(), 4);
|
||||
let nu0 = math::CS2 * (lv[0].tau - 0.5);
|
||||
for (k, p) in s.patches.iter().enumerate() {
|
||||
let (par, g) = (&lv[k], &lv[k + 1]);
|
||||
assert!(p.ax >= 2 && p.bx + 2 <= par.nx && p.ay >= 1 && p.by + 2 <= par.ny);
|
||||
assert_eq!(g.nx, p.r * (p.bx - p.ax) + 1);
|
||||
assert_eq!(g.scale, par.scale * p.r as R);
|
||||
let nu = math::CS2 * (g.tau - 0.5) / g.scale;
|
||||
assert!((nu - nu0).abs() < 1e-14 * nu0, "ν уровня {} = {nu}, L0 = {nu0}", k + 1);
|
||||
}
|
||||
// тело лежит внутри самого тонкого уровня
|
||||
let (cx, cy) = s.scene.center();
|
||||
let f = lv[2];
|
||||
let (x, y) = ((cx - f.ox) * f.scale, (cy - f.oy) * f.scale);
|
||||
assert!(x > 0.0 && x < f.nx as R && y > 0.0 && y < f.ny as R);
|
||||
}
|
||||
|
||||
/// Однородный поток без тела — неподвижная точка всей схемы: столкновения, переноса,
|
||||
/// стенок, входа/выхода и связки трёх уровней. Любая ошибка в рамке, рестрикции или
|
||||
/// пересчёте неравновесия между уровнями сдвинула бы его с места.
|
||||
#[test]
|
||||
fn uniform_flow_is_fixed_point_across_levels() {
|
||||
let mut s = spec_of(&[
|
||||
"--shape", "cylinder", "--size", "6", "--nx", "120", "--ny", "80", "--re", "50",
|
||||
"--levels", "2:3,8,3; 2:1.5,4,1.5", "--sponge-len", "0",
|
||||
]);
|
||||
s.scene.bodies.clear();
|
||||
let u = s.units.u_lat;
|
||||
let mut sim = cpu::Sim::new(s);
|
||||
for _ in 0..200 {
|
||||
sim.step();
|
||||
}
|
||||
for (k, l) in sim.levels.iter().enumerate() {
|
||||
let mut worst: R = 0.0;
|
||||
for n in 0..l.nx * l.ny {
|
||||
let (rho, ux, uy) = l.macros_at(n);
|
||||
worst = worst.max((rho - 1.0).abs()).max((ux - u).abs()).max(uy.abs());
|
||||
}
|
||||
assert!(worst < 1e-12, "уровень {k}: отклонение {worst:.3e}");
|
||||
}
|
||||
}
|
||||
|
||||
/// Три уровня (×2, ×2) против одного патча ×4 с тем же разрешением тела — стационарное
|
||||
/// обтекание при Re = 20 обязано дать тот же Cd. Различие — только в огрублении снаружи.
|
||||
#[test]
|
||||
#[ignore = "длинный прогон (~10 с): cargo test --release -- --ignored levels_match"]
|
||||
fn levels_match_single_patch_on_steady_cylinder() {
|
||||
let common = [
|
||||
"--shape", "cylinder", "--size", "6", "--nx", "180", "--ny", "120", "--re", "20",
|
||||
"--wall", "bouzidi", "--steps", "2500",
|
||||
];
|
||||
let cd = |extra: &[&str]| {
|
||||
let mut a: Vec<&str> = common.to_vec();
|
||||
a.extend_from_slice(extra);
|
||||
let s = spec_of(&a);
|
||||
let mut sim = cpu::Sim::new(s.clone());
|
||||
let d = sim.force_ref_size();
|
||||
let mut acc = Vec::new();
|
||||
for t in 0..s.steps {
|
||||
let r = sim.step();
|
||||
if t > s.steps * 4 / 5 {
|
||||
acc.push(math::coefficients(r.fx, r.fy, r.tz, s.units.u_lat, d).0);
|
||||
}
|
||||
}
|
||||
acc.iter().sum::<R>() / acc.len() as R
|
||||
};
|
||||
let one = cd(&["--levels", "4:2.5,6,2.5"]);
|
||||
let two = cd(&["--levels", "2:2.5,6,2.5; 2:1.5,4,1.5"]);
|
||||
println!("Cd: один патч ×4 {one:.5}, цепочка ×2×2 {two:.5}");
|
||||
assert!((one - two).abs() < 0.01 * one, "Cd {one} против {two}");
|
||||
}
|
||||
}
|
||||
|
||||
@@ -290,6 +290,43 @@ pub fn beta_profile(nx: usize, beta0: R, sponge_in: usize, sponge_out: usize, mu
|
||||
.collect()
|
||||
}
|
||||
|
||||
/// Поле β по ВСЕМ узлам уровня (узел = y·nx + x): губки у входа и выхода, как в
|
||||
/// `beta_profile`, плюс боковые губки шириной `sponge_side` рядов у верхней и нижней стенок.
|
||||
/// Губки не складываются — в каждом узле берётся сильнейшая. При `sponge_side = 0` значения
|
||||
/// совпадают с `beta_profile` побитово.
|
||||
pub fn beta_field(
|
||||
nx: usize,
|
||||
ny: usize,
|
||||
beta0: R,
|
||||
sponge_in: usize,
|
||||
sponge_out: usize,
|
||||
sponge_side: usize,
|
||||
mult: R,
|
||||
) -> Vec<R> {
|
||||
let cols = beta_profile(nx, beta0, sponge_in, sponge_out, mult);
|
||||
if sponge_side == 0 || mult == 1.0 {
|
||||
let mut out = Vec::with_capacity(nx * ny);
|
||||
for _ in 0..ny {
|
||||
out.extend_from_slice(&cols);
|
||||
}
|
||||
return out;
|
||||
}
|
||||
let nu = nu_of_beta(beta0);
|
||||
let strength = |b: R| if b == beta0 { 0.0 } else { (nu_of_beta(b) / nu - 1.0) / (mult - 1.0) };
|
||||
let w = sponge_side as R;
|
||||
let mut out = Vec::with_capacity(nx * ny);
|
||||
for y in 0..ny {
|
||||
let lo = smoothstep((w - y as R) / w);
|
||||
let hi = smoothstep((y as R - (ny - 1) as R + w) / w);
|
||||
let sy = lo.max(hi);
|
||||
for &bc in &cols {
|
||||
// где боковая губка не сильнее продольной, остаётся значение столбца как есть
|
||||
out.push(if sy <= strength(bc) { bc } else { beta_of_nu(nu * (1.0 + (mult - 1.0) * sy)) });
|
||||
}
|
||||
}
|
||||
out
|
||||
}
|
||||
|
||||
/// Кинематическая вязкость по β, формула (5): ν = c_s²(1/(2β) − 1/2).
|
||||
#[inline]
|
||||
pub fn nu_of_beta(beta: R) -> R {
|
||||
@@ -985,6 +1022,195 @@ pub fn grad_target_velocity_term(q: R, u_far: R, u_wall: R) -> R {
|
||||
(q * u_far + u_wall) / (1.0 + q)
|
||||
}
|
||||
|
||||
// ─────────────────────────────────────────────────────────────────────────────
|
||||
// Пристеночная функция
|
||||
// ─────────────────────────────────────────────────────────────────────────────
|
||||
|
||||
/// Пристеночная функция поверх модели стенки: чем задаётся профиль скорости между стенкой и
|
||||
/// первым жидким узлом, когда сетка пограничный слой не разрешает.
|
||||
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
|
||||
pub enum WallFunction {
|
||||
/// Нет: стенка — обычное прилипание выбранной модели.
|
||||
None,
|
||||
/// Закон Spalding (1961) — ОДНА формула на весь пристеночный слой, от вязкого подслоя
|
||||
/// (u⁺ = y⁺) до логарифмической зоны. Поэтому он безопасен и там, где первая клетка
|
||||
/// лежит в подслое: тогда профиль сам сводится к линейному, то есть к прилипанию.
|
||||
Spalding,
|
||||
}
|
||||
|
||||
impl WallFunction {
|
||||
pub fn from_str(s: &str) -> Option<Self> {
|
||||
Some(match s.to_ascii_lowercase().as_str() {
|
||||
"none" | "off" => WallFunction::None,
|
||||
"spalding" => WallFunction::Spalding,
|
||||
_ => return None,
|
||||
})
|
||||
}
|
||||
pub const ALL: [&'static str; 2] = ["none", "spalding"];
|
||||
pub fn is_on(&self) -> bool {
|
||||
!matches!(self, WallFunction::None)
|
||||
}
|
||||
pub fn name(&self) -> &'static str {
|
||||
match self {
|
||||
WallFunction::None => "none",
|
||||
WallFunction::Spalding => "spalding",
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/// Постоянные закона стенки (Spalding 1961, классические значения).
|
||||
pub const KAPPA: R = 0.41;
|
||||
pub const LOG_B: R = 5.2;
|
||||
|
||||
/// y⁺(u⁺) по Spalding:
|
||||
///
|
||||
/// y⁺ = u⁺ + e^{−κB}·[e^{κu⁺} − 1 − κu⁺ − (κu⁺)²/2 − (κu⁺)³/6]
|
||||
///
|
||||
/// Возвращает (y⁺, dy⁺/du⁺). Функция монотонно растёт и выпукла при u⁺ > 0.
|
||||
#[inline]
|
||||
pub fn spalding_yplus(up: R) -> (R, R) {
|
||||
let k = KAPPA * up;
|
||||
let e = (-KAPPA * LOG_B).exp();
|
||||
let ek = k.exp();
|
||||
let y = up + e * (ek - 1.0 - k - 0.5 * k * k - k * k * k / 6.0);
|
||||
let dy = 1.0 + e * KAPPA * (ek - 1.0 - k - 0.5 * k * k);
|
||||
(y, dy)
|
||||
}
|
||||
|
||||
/// u⁺(y⁺) — обращение закона Spalding методом Ньютона.
|
||||
///
|
||||
/// Старт берётся из нужной асимптоты: в подслое u⁺ ≈ y⁺, в логарифмической зоне
|
||||
/// u⁺ ≈ ln(y⁺)/κ + B. Функция y⁺(u⁺) выпукла, поэтому Ньютон от старта не меньше корня
|
||||
/// сходится монотонно; на всякий случай итерация зажата снизу нулём.
|
||||
pub fn spalding_uplus(yp: R) -> R {
|
||||
if yp <= 0.0 {
|
||||
return 0.0;
|
||||
}
|
||||
let mut u = if yp < 11.0 { yp } else { (yp.ln() / KAPPA + LOG_B).max(1e-6) };
|
||||
for _ in 0..50 {
|
||||
let (y, dy) = spalding_yplus(u);
|
||||
let du = (y - yp) / dy;
|
||||
u = (u - du).max(0.0);
|
||||
if du.abs() <= 1e-12 * (1.0 + u) {
|
||||
break;
|
||||
}
|
||||
}
|
||||
u
|
||||
}
|
||||
|
||||
/// Динамическая скорость u_τ по касательной скорости `u` в точке на расстоянии `y` от стенки
|
||||
/// и кинематической вязкости `nu` (всё в решёточных единицах).
|
||||
///
|
||||
/// Неизвестным удобнее брать x = u⁺ = u/u_τ: тогда y⁺ = y·u/(ν·x), и уравнение
|
||||
/// x·y⁺(x) = Re_y, где Re_y = y·u/ν, — с монотонной левой частью. Старт x = √Re_y — точный
|
||||
/// корень в подслое (там y⁺(x) = x).
|
||||
pub fn spalding_utau(u: R, y: R, nu: R) -> R {
|
||||
if u <= 1e-14 || y <= 0.0 {
|
||||
return 0.0;
|
||||
}
|
||||
let re_y = y * u / nu;
|
||||
let mut x = re_y.sqrt();
|
||||
for _ in 0..60 {
|
||||
let (g, dg) = spalding_yplus(x);
|
||||
let h = x * g - re_y;
|
||||
let dh = g + x * dg;
|
||||
let dx = h / dh;
|
||||
x = (x - dx).max(1e-9);
|
||||
if dx.abs() <= 1e-12 * (1.0 + x) {
|
||||
break;
|
||||
}
|
||||
}
|
||||
u / x
|
||||
}
|
||||
|
||||
/// Итог пристеночной функции для одного граничного узла.
|
||||
#[derive(Clone, Copy, Debug, Default)]
|
||||
pub struct WallFnOut {
|
||||
/// Модельная касательная скорость в самом граничном узле (вдоль `e`).
|
||||
pub u_node: R,
|
||||
/// Динамическая скорость.
|
||||
pub u_tau: R,
|
||||
/// y⁺ граничного узла: y_w·u_τ/ν.
|
||||
pub y_plus: R,
|
||||
/// Единичное касательное направление течения у стенки.
|
||||
pub ex: R,
|
||||
pub ey: R,
|
||||
/// Касательная скорость в точке отбора.
|
||||
pub u_sample: R,
|
||||
}
|
||||
|
||||
/// Пристеночная функция в одном узле: по скорости (usx, usy) в точке отбора на расстоянии
|
||||
/// `y_s` от стенки и нормали (nx, ny) — u_τ и модельная касательная скорость в самом узле на
|
||||
/// расстоянии `y_w`. `None`, если касательной скорости в точке отбора нет вовсе.
|
||||
#[allow(clippy::too_many_arguments)]
|
||||
pub fn wall_function(
|
||||
usx: R,
|
||||
usy: R,
|
||||
nx: R,
|
||||
ny: R,
|
||||
y_w: R,
|
||||
y_s: R,
|
||||
nu: R,
|
||||
) -> Option<WallFnOut> {
|
||||
let un = usx * nx + usy * ny;
|
||||
let tx = usx - un * nx;
|
||||
let ty = usy - un * ny;
|
||||
let ut = (tx * tx + ty * ty).sqrt();
|
||||
if ut <= 1e-14 {
|
||||
return None;
|
||||
}
|
||||
let u_tau = spalding_utau(ut, y_s, nu);
|
||||
let y_plus = y_w * u_tau / nu;
|
||||
let u_node = u_tau * spalding_uplus(y_plus);
|
||||
Some(WallFnOut { u_node, u_tau, y_plus, ex: tx / ut, ey: ty / ut, u_sample: ut })
|
||||
}
|
||||
|
||||
/// Скорость скольжения стенки для полинковой схемы (Bouzidi), согласованная по НАПРЯЖЕНИЮ:
|
||||
/// отскок от подвижной стенки передаёт ей импульс ≈ ρν(u_узла − u_s)/y_w, и чтобы при
|
||||
/// модельной скорости узла это было модельным τ_w = ρu_τ², стенка должна двигаться со
|
||||
/// скоростью
|
||||
///
|
||||
/// u_s = u_node − u_τ²·y_w/ν = u_τ·(u⁺(y⁺_w) − y⁺_w), зажато в [−u_sample, u_sample].
|
||||
///
|
||||
/// В вязком подслое u⁺ = y⁺ и скольжение ТОЖДЕСТВЕННО ноль — обычное прилипание, побитово.
|
||||
/// В буферной и логарифмической зонах u⁺ < y⁺, и стенка «отступает» против потока (u_s < 0),
|
||||
/// добирая трение до τ_w. Это та же пара (u_node, τ_w), что моментные схемы ставят через
|
||||
/// целевую скорость и тензор давлений, — так сравнение схем честное.
|
||||
///
|
||||
/// Разрешённую скорость узла вместо модельной брать нельзя: у искривлённой стенки она не
|
||||
/// лежит на прямой от стенки к точке отбора, и скольжение переставало обращаться в ноль
|
||||
/// даже при y⁺ ≪ 1 (цилиндр Re = 20: −2 % по Cd против прилипания).
|
||||
///
|
||||
/// Почему не согласование по СКОРОСТИ (прямая через стенку, узел и точку отбора): решатель
|
||||
/// не имеет турбулентной вязкости, и такое скольжение выходит положительным и большим —
|
||||
/// стенка становится почти свободной, исчезают трение и отрыв. Замерено: цилиндр Re = 10⁴,
|
||||
/// D = 16 давал Cd = 0.016. Зажим по модулю скорости в точке отбора держит скорость стенки
|
||||
/// в пределах скоростей самого течения — иначе при y⁺ ≫ 30 она выходила бы за число Маха.
|
||||
#[inline]
|
||||
pub fn wall_slip(w: &WallFnOut, y_w: R, nu: R) -> R {
|
||||
(w.u_node - w.u_tau * w.u_tau * y_w / nu).clamp(-w.u_sample, w.u_sample)
|
||||
}
|
||||
|
||||
/// Поправка градиента скорости граничного узла под пристеночную функцию: нормальная
|
||||
/// производная касательной скорости заменяется модельной u_τ²/ν, остальные компоненты
|
||||
/// остаются конечно-разностными. Добавка имеет вид d·(e⊗n), поэтому дивергенцию не меняет
|
||||
/// (e ⊥ n). Возвращает исправленные (∂u/∂x, ∂u/∂y, ∂v/∂x, ∂v/∂y).
|
||||
#[allow(clippy::too_many_arguments)]
|
||||
#[inline]
|
||||
pub fn wall_function_gradient(
|
||||
g: (R, R, R, R),
|
||||
w: &WallFnOut,
|
||||
nx: R,
|
||||
ny: R,
|
||||
nu: R,
|
||||
) -> (R, R, R, R) {
|
||||
let (dudx, dudy, dvdx, dvdy) = g;
|
||||
let (ex, ey) = (w.ex, w.ey);
|
||||
let g_en = ex * (dudx * nx + dudy * ny) + ey * (dvdx * nx + dvdy * ny);
|
||||
let d = w.u_tau * w.u_tau / nu - g_en;
|
||||
(dudx + d * ex * nx, dudy + d * ex * ny, dvdx + d * ey * nx, dvdy + d * ey * ny)
|
||||
}
|
||||
|
||||
// ─────────────────────────────────────────────────────────────────────────────
|
||||
// Сила на теле
|
||||
// ─────────────────────────────────────────────────────────────────────────────
|
||||
@@ -1529,4 +1755,82 @@ mod tests {
|
||||
}
|
||||
inv
|
||||
}
|
||||
|
||||
/// Закон Spalding в обеих асимптотах: в вязком подслое u⁺ = y⁺, в логарифмической зоне
|
||||
/// u⁺ = ln(y⁺)/κ + B. Формула одна, асимптоты — её свойства, а не ветки кода.
|
||||
#[test]
|
||||
fn spalding_matches_both_asymptotes() {
|
||||
for yp in [0.01, 0.1, 1.0] {
|
||||
let up = spalding_uplus(yp);
|
||||
assert!((up - yp).abs() < 1e-3 * yp.max(1e-3), "подслой: y⁺={yp}, u⁺={up}");
|
||||
}
|
||||
for yp in [300.0, 1000.0, 5000.0] {
|
||||
let up = spalding_uplus(yp);
|
||||
let log = yp.ln() / KAPPA + LOG_B;
|
||||
assert!((up - log).abs() < 0.02 * log, "log-зона: y⁺={yp}, u⁺={up}, лог {log}");
|
||||
}
|
||||
}
|
||||
|
||||
/// Обращение точное и монотонное: y⁺(u⁺(y⁺)) = y⁺ на всём диапазоне.
|
||||
#[test]
|
||||
fn spalding_inverse_roundtrip_and_monotone() {
|
||||
let mut prev = 0.0;
|
||||
let mut yp = 1e-3;
|
||||
while yp < 1e5 {
|
||||
let up = spalding_uplus(yp);
|
||||
assert!(up > prev, "u⁺ обязан расти: y⁺={yp}");
|
||||
prev = up;
|
||||
let (back, _) = spalding_yplus(up);
|
||||
assert!((back - yp).abs() < 1e-9 * yp.max(1.0), "y⁺={yp} → u⁺={up} → {back}");
|
||||
yp *= 1.7;
|
||||
}
|
||||
}
|
||||
|
||||
/// u_τ по скорости в точке: подставленная обратно, она даёт ту же скорость. В подслое
|
||||
/// u_τ = √(ν·u/y) — линейный профиль, то есть ровно прилипание.
|
||||
#[test]
|
||||
fn spalding_utau_consistent() {
|
||||
let nu = 1e-3;
|
||||
for (u, y) in [(0.05, 1.5), (0.05, 30.0), (1e-4, 1.0), (0.1, 200.0)] {
|
||||
let ut = spalding_utau(u, y, nu);
|
||||
let up = spalding_uplus(y * ut / nu);
|
||||
assert!((ut * up - u).abs() < 1e-9 * u, "u={u} y={y}: u_τ·u⁺ = {}", ut * up);
|
||||
}
|
||||
let (u, y) = (1e-5, 0.5);
|
||||
let ut = spalding_utau(u, y, nu);
|
||||
assert!((ut - (nu * u / y).sqrt()).abs() < 1e-6 * ut, "подслой: u_τ = {ut}");
|
||||
assert_eq!(spalding_utau(0.0, 1.0, nu), 0.0);
|
||||
}
|
||||
|
||||
/// В вязком подслое пристеночная функция обязана свестись к прилипанию: модельная
|
||||
/// скорость узла — линейная интерполяция от стенки, скольжение — ноль, а модельный
|
||||
/// градиент — тот же наклон прямой.
|
||||
#[test]
|
||||
fn wall_function_reduces_to_no_slip_in_sublayer() {
|
||||
let nu = 0.05; // низкий Re: y⁺ первого узла заведомо < 1
|
||||
let (y_w, y_s) = (0.4, 1.4);
|
||||
let (ux, uy) = (0.01, 0.0); // стенка снизу, нормаль +y
|
||||
let w = wall_function(ux, uy, 0.0, 1.0, y_w, y_s, nu).unwrap();
|
||||
assert!(w.y_plus < 1.0, "y⁺ = {}", w.y_plus);
|
||||
let lin = ux * y_w / y_s;
|
||||
assert!((w.u_node - lin).abs() < 1e-3 * lin, "u_node {} против {lin}", w.u_node);
|
||||
assert!(wall_slip(&w, y_w, nu).abs() < 1e-3 * lin);
|
||||
let g = wall_function_gradient((0.0, 0.0, 0.0, 0.0), &w, 0.0, 1.0, nu);
|
||||
assert!((g.1 - ux / y_s).abs() < 1e-3 * ux / y_s, "∂u/∂y = {}", g.1);
|
||||
assert_eq!((g.0, g.2, g.3), (0.0, 0.0, 0.0));
|
||||
}
|
||||
|
||||
/// В логарифмической зоне разрешённый градиент меньше модельного напряжения, и стенка
|
||||
/// отступает против потока: скольжение отрицательно, но по модулю не больше скорости
|
||||
/// в точке отбора.
|
||||
#[test]
|
||||
fn wall_function_slips_in_log_region() {
|
||||
let nu = 1e-5;
|
||||
let (y_w, y_s) = (0.5, 1.5);
|
||||
let w = wall_function(0.05, 0.0, 0.0, 1.0, y_w, y_s, nu).unwrap();
|
||||
assert!(w.y_plus > 30.0, "y⁺ = {}", w.y_plus);
|
||||
let s = wall_slip(&w, y_w, nu);
|
||||
assert!(s < 0.0 && s >= -w.u_sample, "скольжение {s}, u_sample {}", w.u_sample);
|
||||
assert!(wall_function(0.0, 0.03, 0.0, 1.0, y_w, y_s, nu).is_none(), "нормальный поток");
|
||||
}
|
||||
}
|
||||
|
||||
Reference in New Issue
Block a user