Начальный коммит: Vulkan-редактор SimVulcan + исследование KBC-LBM
Состояние на момент заведения репозитория. C++ приложение (src/, shaders/, tests/) — минимальный редактор 3D-моделей на Vulkan 1.3: орбитальная камера, три опорные сетки через начало координат, загрузка .obj с режимами отображения. Весь Vulkan изолирован в src/vk/. Исследование (docs/) — оригинальные статьи по KBC (docs/origins) и Python-решатель D2Q9 KBC-N1 с AMR 2x и SDF+Bouzidi (docs/theory). В решателе перед коммитом исправлены дефекты, найденные сверкой с первоисточниками: относительный порог знаменателя энтропийного стабилизатора (абсолютный вырождал KBC в LBGK на 77-99% узлов), заворот вход/выход в углах домена, диагностика средней плотности по фиктивным узлам тела, зашитый refine=2. Подробности — docs/theory/solver_2x_sdf/README.md. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
This commit is contained in:
Claude Opus 5
@@ -0,0 +1,160 @@
|
||||
# solver_2x_sdf — модель цилиндра 2×+SDF, разложенная по компонентам
|
||||
|
||||
Самодостаточный LBM-решатель (D2Q9, KBC) для обтекания цилиндра: грубый уровень **L0**
|
||||
плюс один вложенный патч **L1** (измельчение ×2). Граница цилиндра — **SDF + Bouzidi (no-slip)**;
|
||||
вход/выход — **Zou-He** (скорость/давление); верх/низ канала — **free-slip (specular)**.
|
||||
Только GPU/CuPy — **CPU-фолбэка нет**. Запуск — из этой папки: `python run.py`
|
||||
(серия по блокировке: `python run_blockage.py`).
|
||||
|
||||
## Карта файлов (каждый компонент изолирован и легко меняется)
|
||||
|
||||
| Файл | Ответственность |
|
||||
|---|---|
|
||||
| `backend.py` | жёсткое требование CuPy/GPU, `DTYPE`, device-скаляры, `to_cpu` |
|
||||
| `config.py` | все параметры модели + производные (τ, β, R, DL, cy/патч из Ny); `from_env()` |
|
||||
| `lattice.py` | D2Q9 (c_i, w_i, OPP) + KBC-проектор `Ps` на сдвиг-моменты |
|
||||
| `equilibrium.py` | `feq` (product-form, без `pow`), `macros` (ρ, u) |
|
||||
| `collision.py` | **столкновение**: `kbc_collide` (энтропийный KBC-N1) — единственный оператор |
|
||||
| `streaming.py` | **перенос**: `stream` (пул-схема через roll) |
|
||||
| `geometry.py` | маски тела/стенок и SDF для L0 и L1 |
|
||||
| `boundary.py` | `build_bc`/`apply_bc` (SDF+Bouzidi) + `channel_bc` (Zou-He) + `free_slip_walls` |
|
||||
| `forces.py` | **силы**: галилей-инвариантный обмен импульсом `force_gmem` (Fx, Fy, **торк Tz**) + `coefficients` (Cd, Cl, **Cm**) |
|
||||
| `forces_stress.py` | **кросс-чек силы**: ∮σ·n ds по контуру R+δ вокруг тела (независимый эстиматор) |
|
||||
| `amr.py` | связка L0↔L1: `patch`/`ghost`/`fill`/`restrict` |
|
||||
| `solver.py` | оркестратор: persistent-буферы, `step()`, `run()` (+ CUDA-граф) |
|
||||
| `diagnostics.py` | `analyze`, St (параболическая интерполяция пика), кросс-чек, фиты серии по β |
|
||||
| `visualization.py` | гифка |ω|, png-панель рядов (`series_2x_sdf.png`), прогресс-бар |
|
||||
| `run.py` | точка входа (один прогон) |
|
||||
| `run_blockage.py` | серия по блокировке Ny↑ → β↓, экстраполяция β→0, сводный npz+png |
|
||||
| `run_factors.py` | факторное исследование постановки (вход/выход/uy-выход); реестр `out/factors.csv` |
|
||||
|
||||
## Где какая математика (для ревизии)
|
||||
|
||||
- **Столкновение** — `collision.py`. KBC-N1: `f ← f − β(2Δs + γΔh)`, где `Δs = Ps·(f−feq)` —
|
||||
проекция на сдвиг-моменты, `γ` — энтропийный лимитер. Чтобы поэкспериментировать:
|
||||
переключить `cfg.collision="bgk"`, либо добавить TRT новой функцией и зарегистрировать в `get`.
|
||||
- **Перенос** — `streaming.py`. Чистая пул-схема; на физику влияет только корректность сдвигов.
|
||||
- **Силы** — `forces.py`. Обмен импульсом по линкам тела:
|
||||
`F = Σ_links c_i (f_i^{после столкн.} + f_ī^{после стриминга})`. Второй член — из поля ПОСЛЕ
|
||||
стриминга (вернувшаяся Bouzidi-популяция); иначе ведущий симметричный член сокращается и Cd врёт.
|
||||
Сила/момент меряются на **каждом** подшаге L1 и усредняются (анти-алиасинг; средний Cd не меняется).
|
||||
**Торк** `Tz` — вокруг центра тела, точка приложения — пересечение линка со стенкой `r_f + q·c_i`
|
||||
(Bouzidi-доля q уже лежит в BC); `Cm = 2Tz/(U²·DL²)`, для цилиндра ⟨Cm⟩≈0 — контроль симметрии.
|
||||
- **Кросс-чек силы** — `forces_stress.py`. Независимый эстиматор — баланс импульса
|
||||
контрольного объёма: `F = ∮[σ·n − ρu(u·n)] ds` по окружности R+δ (δ=2 тонкие ячейки);
|
||||
`σ = −p′·I − (1−1/(2τ₁))·Π^neq_dev` (строго девиаторная проекция: след Π^neq в KBC
|
||||
релаксирует с γβ, а не 1/τ). Конвективный член −ρu(u·n) обязателен при δ>0 — без него
|
||||
систематика +4% по ⟨Cd⟩. Считается раз в `stress_every` шагов **вне CUDA-графа**; сравнение
|
||||
по средним (⟨Cd⟩, rms Cl) — мгновенные ряды сдвинуты по фазе кольцом R..R+δ.
|
||||
Согласие двух методов = считывание силы корректно.
|
||||
- **Серия по блокировке** — `run_blockage.py`. D фиксирован, Ny ∈ {90,128,180} → β ∈ {0.178,
|
||||
0.125, 0.089}; `cy` и y-границы патча выводятся из Ny автоматически (`config.__post_init__`).
|
||||
Экстраполяция к β→0 двумя фитами (a+b·β и a+b·β²; спред интерсептов — неопределённость) —
|
||||
прямое сравнение с безграничной литературой без ручной поправки (1−β)².
|
||||
|
||||
## ⚠ Сверка с первоисточниками (2026-08-14): найден дефект, все числа ниже требуют перепрогона
|
||||
|
||||
Сверка реализации со статьями в `docs/origins` (Karlin/Bösch/Chikatamarla, PRE **90**, 031302(R)
|
||||
2014 — принцип; Bösch/Chikatamarla/Karlin, PRE **92**, 043309 2015 — алгоритм; arXiv:1507.02509 —
|
||||
2D-реализация). **Ядро KBC-N1 подтверждено точным**: проектор `Ps` совпадает с аналитической
|
||||
s-частью из ур.(10) 2D-статьи (ранг 2, идемпотентен, несёт ровно {N, Π_xy}); равновесие —
|
||||
точный энтропийный максимизатор (сходится с численной максимизацией S при фикс. ρ, ρu);
|
||||
γ по замкнутой формуле ур.(17)/(39) совпадает с корнем точного ур.(15); сдвиг-моменты
|
||||
релаксируют ровно с 2β при любом γ; вязкость по затуханию сдвиговой волны воспроизводится
|
||||
с ошибкой <0.1%.
|
||||
|
||||
**Дефект — порог `GEPS` в `kbc_collide`.** Он был АБСОЛЮТНЫМ (1e-6 в float32, режим по умолчанию),
|
||||
тогда как `den = ⟨Δh|Δh⟩` квадратична по неравновесию и при разрешённой физике равна ~1e-7…1e-9.
|
||||
Замерено на боевой конфигурации: откат `γ←2` срабатывал на **77–99% узлов**, то есть решатель
|
||||
считал **чистый LBGK**, а не KBC. На эталоне 2D-статьи (двойной периодический сдвиговый слой,
|
||||
Re=30000, N=128, разд. VII) это давало **+9.3%** по энстрофии на t=t_c (0.6599 против эталонных
|
||||
0.6035); после исправления fp32 даёт 0.6038. Порог сделан относительным (`B.GREL`, доля от ‖Δ‖²);
|
||||
на тех же полях он теперь не срабатывает нигде (min den/‖Δ‖² ≈ 1.5e-3 против порога 1e-8).
|
||||
|
||||
**Следствие: все количественные результаты ниже (Cd, rms Cl, St, серия по блокировке,
|
||||
факторное исследование №0/№1) получены на дефектной ветке и подлежат перепрогону.**
|
||||
Качественные выводы про ГУ (мода акустического резонатора, роль конвективного члена
|
||||
в кросс-чеке) дефектом не затрагиваются.
|
||||
|
||||
Исправлено там же:
|
||||
- **углы домена**: `stream()` на `roll` периодичен по обеим осям, а Zou-He стоял на срезе `1:-1` —
|
||||
в 4 угловых узлах вход был напрямую связан с выходом (`f[1,0,0] == post[1,0,-1]`, величина ~0.12).
|
||||
Порядок изменён на `free_slip_walls` → `channel_bc` по всему столбцу; заворот закрыт полностью.
|
||||
- **диагностика ⟨ρ⟩** считалась по всему домену, включая фиктивные узлы внутри цилиндра
|
||||
(смещение −5.4e-4 — ровно тот разряд, в котором ⟨ρ⟩ и цитируется). Теперь только по жидкости.
|
||||
- **`refine` был зашит как 2** в `tau1`, `R01` и радиусе тела на L1/контуре σ·n. Обобщено
|
||||
(при r=2 значения бит-в-бит прежние; при r=3 вязкость уровней теперь согласована).
|
||||
|
||||
**Отдельно — не дефект, а следствие выбора модели.** В KBC-N1 след тензора напряжений T отнесён
|
||||
к h-части, поэтому объёмная вязкость плавает: ξ = c_s²(1/(γβ) − ½) (ур.36 статьи 2015, ур.57
|
||||
2D-статьи). Замер на боевой конфигурации: γ ∈ [−4.1, 4.6], медиана ξ/ν ≈ 17, и на **2.1% узлов
|
||||
ξ < 0** — акустическое АНТИдемпфирование. Это прямо релевантно диагнозу «недодемпфированный
|
||||
акустический резонатор» ниже: KBC-N2 (T в s-части) даёт ξ = ν фиксированно и является
|
||||
каноническим средством именно от этой проблемы. Модель не менялась — решение за автором.
|
||||
|
||||
## Статус: считывание силы валидировано; идёт факторное исследование постановки
|
||||
|
||||
Оператор столкновения зафиксирован — **KBC-N1 с энтропийным стабилизатором** (канонический
|
||||
Karlin/Bösch/Chikatamarla 2014), без замен (никаких BGK/TRT).
|
||||
**Полная документация исследования — в ноутбуке `../kbc_lbm.ipynb` (части VI, разделы 22–24).**
|
||||
|
||||
**Считывание силы доказано корректным** (прогон 100k, 174×90, 2026-06-09):
|
||||
- кросс-чек GMEM vs ∮σ·n: dCd=4.1% (систематика объяснена конвективным членом −ρu(u·n),
|
||||
добавлен в `forces_stress.py` — ожидаемое расхождение ≲1%);
|
||||
- ⟨Cm⟩=−0.00003 (≈0 — симметрия), ⟨ρ⟩=1.001 стабильна, Cd=1.799 воспроизводит прежнюю
|
||||
итерацию бит-в-бит (осреднение по подшагам не сместило среднее).
|
||||
|
||||
**Эксперимент №0 (серия по блокировке) — гипотеза бокового стеснения ОПРОВЕРГНУТА:**
|
||||
Cd(β) немонотонен: 1.799 / 1.771 / **1.870** при Ny=90/128/180; rms Cl на широком канале
|
||||
вырос (0.835), НЧ-модуляция (окна Cd 1.68–1.98). При этом **St(β→0)=0.181 ≈ лит. 0.183** —
|
||||
кинематика верна, завышены давленческие амплитуды. Диагноз: продольные границы
|
||||
(вход 2.5D с жёстким профилем; выход 8.3D Zou-He с жёстким uy=0 — отражает вихри).
|
||||
|
||||
**Эксперимент №1 (прогнан 2026-06-10) — факторное исследование** `run_factors.py`, 6 прогонов
|
||||
при Ny=90 фикс. Итоги (полный разбор — ноутбук, раздел 24; реестр — `out/factors.csv`):
|
||||
- конвективный член подтверждён: кросс-чек B0 4.1%→**1.0%**, F4 — **0.6%**;
|
||||
- **F4 (мягкий uy-выход) — лучший по всем метрикам**: Cd 1.799→1.754, rms Cl 0.705→**0.557**,
|
||||
St·(1−β)=**0.183** (точно лит.), ⟨ρ⟩=1.0000, нулевая НЧ-модуляция;
|
||||
- **длинные домены без демпфера НЕвалидны** (F1/F2/F3 → смещённый режим ⟨ρ⟩≈1.66; F5 → NaN):
|
||||
пара «вход-скорость / выход-давление» — недодемпфированный акустический резонатор
|
||||
(затухание моды ∝ν(π/Nx)²; период осцилляций ⟨ρ⟩ у F5 ≈ 2Nx/c_s — прямая улика).
|
||||
|
||||
**Эксперимент №2 (в работе) — губка перед выходом**: плавный рост ν (×30, smoothstep) в
|
||||
последних `sponge_len` столбцах L0; β(x) — поле (graph-safe), патч L1 губки не видит.
|
||||
Факторы: F6_sponge (база+uy-extrap+губка 32 — валидация против F4), F7_long_sponge
|
||||
(390×90, вход 6D/выход 18.3D — целевая постановка), F8_long_sp_uy0 (изоляция вклада губки).
|
||||
Критерии: ⟨ρ⟩≈1.000 у всех; F7: Cd≈1.5–1.6, rms Cl<0.45. Затем эксперимент №3 — серия по β
|
||||
на конфигурации F7.
|
||||
|
||||
| прогон | команда | артефакты в `out/` |
|
||||
|---|---|---|
|
||||
| смоук | `AMR_STEPS=2000 python run.py` | graph self-check OK, Cm≈0, кросс-чек ≲1% (с конв. членом) |
|
||||
| основной | `python run.py` | `cyl_2x_sdf.gif`, `series_2x_sdf.png`, `probe_2x_sdf.npz` |
|
||||
| серия по β | `python run_blockage.py` | `probe_blockage_Ny*.npz`, `blockage_summary.npz`, `blockage_extrapolation.png` |
|
||||
| **факторы №2** | `AMR_FACTORS=F6_sponge,F7_long_sponge,F8_long_sp_uy0 python run_factors.py` | `factor_*.npz`, **`factors.csv`** (реестр; при смене схемы старый → `factors_vN.csv`), `factors_summary.png` |
|
||||
|
||||
Env-ручки: `AMR_NX`, `AMR_CX` (продольная геометрия; x-границы патча L1 выводятся из cx),
|
||||
`AMR_OUTLET_UY=zero|extrapolate`, `AMR_SPONGE_LEN`/`AMR_SPONGE_NU` (губка),
|
||||
`AMR_FACTORS=...` (подмножество факторов).
|
||||
|
||||
**Что привело сюда (хронология фиксов):**
|
||||
1. Зонд скорости → L1; оценка St с параболической интерполяцией пика.
|
||||
2. Сила → галилей-инвариантный `force_gmem` (Wen 2014), готов к подвижным телам.
|
||||
3. Стартовое возмущение (поперечный sin-импульс входа) — запуск дорожки за единицы циклов.
|
||||
4. Диагностика по окнам (`convergence_report`, трекинг ⟨ρ⟩) — вскрыла дрейф массы.
|
||||
5. **Ключевой фикс:** грубое ГУ выхода (zero-gradient) копило массу → поток глох (Cd/rms→0 за 100k).
|
||||
Заменено на **Zou-He скоростной вход + давление-выход** (`boundary.channel_bc`) — масса заякорена.
|
||||
6. **Free-slip стенки** (`boundary.free_slip_walls`, specular) + умеренно расширенный домен (меньше
|
||||
блокировка). Цилиндр остаётся no-slip.
|
||||
7. Cm + осреднение по подшагам + кросс-чек ∮σ·n ds + серия по блокировке (текущая итерация).
|
||||
|
||||
**Прямой raw-матч на одной сетке непрактичен** (нужен `β<0.03 → Ny≈500+`); серия из трёх умеренных
|
||||
Ny с экстраполяцией даёт безпоправочное сравнение за ~3 стоимости одного прогона (патч L1 —
|
||||
основная работа — от Ny не зависит).
|
||||
|
||||
## CUDA Graphs
|
||||
|
||||
По умолчанию ВКЛ (`cfg.use_cuda_graph`, env `AMR_CUDAGRAPH=0` — выкл). **Захватываются успешно** (после
|
||||
замены `cupy.stack`→`backend.pack`, отказа от cuBLAS в проекторе KBC и clip-литералов в ядре —
|
||||
всё это убрало host→device при capture). На шаге t=1 — рантайм-самопроверка (граф vs eager); при
|
||||
расхождении/сбое — откат на eager с печатью трейсбэка.
|
||||
@@ -0,0 +1,54 @@
|
||||
# AMR-связка L0 (грубый) ↔ L1 (тонкий, измельчение r).
|
||||
#
|
||||
# Идея: на тонком патче считаем r подшагов по времени на каждый шаг L0. Границы патча («ghost»)
|
||||
# заполняем интерполяцией с грубого уровня (по пространству + по времени между «до» и «после» шага L0),
|
||||
# а внутренность патча после подшагов проецируем обратно на грубый уровень (рестрикция).
|
||||
# Неравновесная часть масштабируется множителями Rcf (грубый→тонкий) и Rfc (тонкий→грубый),
|
||||
# т.к. при смене τ масштаб неравновесия меняется.
|
||||
import backend as B
|
||||
import lattice as L
|
||||
from equilibrium import feq, macros
|
||||
|
||||
cp = B.cp; xp = B.xp; DTYPE = B.DTYPE
|
||||
|
||||
def patch(pNx, pNy, ax, bx, ay, by, r):
|
||||
"""Описание патча: размеры (Nfy,Nfx), индексы/веса билинейной интерполяции с грубого уровня,
|
||||
срезы для рестрикции каждого r-го узла."""
|
||||
Wx, Wy = bx - ax, by - ay
|
||||
Nfx, Nfy = r*Wx + 1, r*Wy + 1
|
||||
FX, FY = cp.meshgrid(ax + cp.arange(Nfx)/r, ay + cp.arange(Nfy)/r) # координаты тонких узлов в L0
|
||||
x0 = cp.floor(FX).astype(cp.int64); y0 = cp.floor(FY).astype(cp.int64)
|
||||
x1 = cp.minimum(x0 + 1, pNx - 1); y1 = cp.minimum(y0 + 1, pNy - 1)
|
||||
return dict(Nfx=int(Nfx), Nfy=int(Nfy), ax=ax, bx=bx, ay=ay, by=by, r=r,
|
||||
x0=x0, y0=y0, x1=x1, y1=y1,
|
||||
tx=(FX - x0).astype(DTYPE), ty=(FY - y0).astype(DTYPE),
|
||||
slx=slice(r, r*Wx, r), sly=slice(r, r*Wy, r))
|
||||
|
||||
def pint(fld, P):
|
||||
"""Билинейная интерполяция поля fld (грубого) в узлы тонкого патча."""
|
||||
return (fld[..., P["y0"], P["x0"]]*(1 - P["tx"])*(1 - P["ty"])
|
||||
+ fld[..., P["y0"], P["x1"]]*P["tx"]*(1 - P["ty"])
|
||||
+ fld[..., P["y1"], P["x0"]]*(1 - P["tx"])*P["ty"]
|
||||
+ fld[..., P["y1"], P["x1"]]*P["tx"]*P["ty"])
|
||||
|
||||
def ghost(pf, P, Rcf):
|
||||
"""Граничные значения тонкого уровня из грубого поля pf: равновесие по интерполированным ρ,u
|
||||
плюс масштабированная неравновесная часть Rcf·neq.
|
||||
u в feq передаём кортежем (без cupy.stack) — H2D ломает CUDA-graph capture."""
|
||||
r, u = macros(pf); neq = pf - feq(r, u)
|
||||
return feq(pint(r, P), (pint(u[0], P), pint(u[1], P))) + Rcf*pint(neq, P)
|
||||
|
||||
def fill(cf, gh):
|
||||
"""Записать ghost-значения gh в рамку (1 ряд) тонкого поля cf."""
|
||||
cf[:, 0, :] = gh[:, 0, :]; cf[:, -1, :] = gh[:, -1, :]
|
||||
cf[:, :, 0] = gh[:, :, 0]; cf[:, :, -1] = gh[:, :, -1]
|
||||
return cf
|
||||
|
||||
def restrict(cf, pf, P, Rfc, fluid):
|
||||
"""Спроецировать тонкое поле cf обратно на грубое pf (каждый r-й узел), масштабируя neq на Rfc.
|
||||
Меняется только внутренняя жидкая часть перекрытия (маска fluid)."""
|
||||
r, u = macros(cf); neq = cf - feq(r, u); slx, sly = P["slx"], P["sly"]
|
||||
nv = feq(r[sly, slx], u[:, sly, slx]) + Rfc*neq[:, sly, slx]
|
||||
cur = pf[:, P["ay"]+1:P["by"], P["ax"]+1:P["bx"]]
|
||||
pf[:, P["ay"]+1:P["by"], P["ax"]+1:P["bx"]] = xp.where(fluid[None], nv, cur)
|
||||
return pf
|
||||
@@ -0,0 +1,40 @@
|
||||
# Бэкенд: ТОЛЬКО GPU/CuPy. CPU-фолбэка нет (по требованию). Если CuPy/GPU недоступны — жёсткая ошибка.
|
||||
# Здесь же — общие device-константы (DTYPE, скаляры 0/1/2) и помощник переноса на host для вывода.
|
||||
import os
|
||||
import cupy as cp
|
||||
|
||||
try:
|
||||
_probe = (cp.zeros(2) + 1).sum(); cp.cuda.Device().synchronize() # реальная операция на device
|
||||
except Exception as e:
|
||||
raise RuntimeError("solver_2x_sdf требует рабочую GPU + CuPy (CPU-фолбэк удалён намеренно).") from e
|
||||
|
||||
xp = cp # единый алиас вычислительного бэкенда
|
||||
DTYPE = cp.float64 if os.environ.get("AMR_FP64") else cp.float32
|
||||
# Порог вырожденности знаменателя γ. ОТНОСИТЕЛЬНЫЙ (доля от ‖Δ‖² = ⟨Δ|Δ⟩), а НЕ абсолютный.
|
||||
# Почему: den = ⟨Δh|Δh⟩ — квадратичная по неравновесию величина, у развитого следа она
|
||||
# ~1e-7…1e-9 при вполне разрешённой физике. Прежний абсолютный порог 1e-6 (fp32) срабатывал
|
||||
# на 77–99% узлов и молча подменял γ на 2, т.е. гнал ЧИСТЫЙ LBGK вместо KBC (замерено:
|
||||
# энстрофия сдвигового слоя Re=3e4 уходила на +9%). den — сумма НЕотрицательных слагаемых,
|
||||
# её относительная точность ~eps типа, поэтому дробный порог одинаков для fp32/fp64.
|
||||
# 1e-8 на шесть порядков ниже наблюдаемого минимума den/‖Δ‖² (~5e-4) — срабатывает только на
|
||||
# истинном вырождении (Δh ≡ 0), где γ всё равно ни на что не влияет, т.к. умножается на Δh.
|
||||
GREL = 1e-8
|
||||
BACKEND = f"CuPy/GPU dtype={'float64' if DTYPE == cp.float64 else 'float32'}"
|
||||
|
||||
# device-скаляры — не создаём host→device на горячем пути (важно для CUDA-graph capture)
|
||||
ZERO = cp.zeros((), DTYPE)
|
||||
ONE = cp.asarray(1.0, DTYPE)
|
||||
TWO = cp.asarray(2.0, DTYPE)
|
||||
|
||||
def to_cpu(a):
|
||||
"""Перенос на host (numpy) — только для диагностики/визуализации, не в солвере."""
|
||||
return cp.asnumpy(a)
|
||||
|
||||
def pack(comps):
|
||||
"""Собрать кортеж одинаковых по форме массивов в один (len,*shape) БЕЗ cupy.stack.
|
||||
cupy.stack/concatenate грузит на device host-массив указателей (H2D) → запрещено при CUDA-graph
|
||||
capture. Здесь: преаллокация + присваивание срезов — только device→device."""
|
||||
out = xp.empty((len(comps),) + comps[0].shape, dtype=comps[0].dtype)
|
||||
for i, c in enumerate(comps):
|
||||
out[i] = c
|
||||
return out
|
||||
@@ -0,0 +1,104 @@
|
||||
# Граничные условия: (1) тело/стенки — SDF + интерполированный отскок Bouzidi; (2) вход/выход канала.
|
||||
#
|
||||
# Bouzidi (линейный): для линка из жидкого узла x_f в твёрдый по направлению i вводится
|
||||
# доля пересечения q = |x_f→стенка| / |x_f→x_solid| ∈ (0,1):
|
||||
# q < 1/2 (есть «дальний» жидкий сосед x_f − c_i): f_ī = 2q·f_i^post + (1−2q)·f_i^post(дальний)
|
||||
# q ≥ 1/2: f_ī = (1/2q)·f_i^post + (1−1/2q)·f_ī^post
|
||||
# Здесь f_i^post — после столкновения, f_ī — то, что вернётся в x_f после переноса.
|
||||
import backend as B
|
||||
import lattice as L
|
||||
|
||||
cp = B.cp; xp = B.xp; DTYPE = B.DTYPE
|
||||
Q = L.Q; Cx = L.Cx; Cy = L.Cy; OPP = L.OPP
|
||||
|
||||
def build_bc(solid, phi, body):
|
||||
"""Сборка структуры ГУ (один раз). Для каждого направления i (1..8):
|
||||
mask — жидкий узел, чей сосед по +c_i твёрдый;
|
||||
q — доля пересечения по SDF (клип [0.02,0.98]);
|
||||
xff — есть ли «дальний» жидкий сосед (для линейной формулы при q<1/2);
|
||||
cl — подмаска линков именно ТЕЛА body (для расчёта силы)."""
|
||||
fluid = ~solid; bc = []
|
||||
for i in range(1, Q):
|
||||
nb_solid = cp.roll(solid, (-Cy[i], -Cx[i]), (0, 1)); mask = fluid & nb_solid
|
||||
cl = mask & cp.roll(body, (-Cy[i], -Cx[i]), (0, 1))
|
||||
phinb = cp.roll(phi, (-Cy[i], -Cx[i]), (0, 1))
|
||||
qq = phi / (phi - phinb) # деление; nan/inf маскируются ниже
|
||||
q = xp.where(mask, cp.clip(qq, 0.02, 0.98), cp.asarray(0.5, DTYPE))
|
||||
xff = mask & (~cp.roll(solid, (Cy[i], Cx[i]), (0, 1)))
|
||||
bc.append((i, OPP[i], mask, q.astype(DTYPE), xff, cl))
|
||||
return bc
|
||||
|
||||
def apply_bc(f, fpost, bc):
|
||||
"""Применить Bouzidi ко всем линкам (полностью через where, GPU-friendly).
|
||||
f — поле ПОСЛЕ стриминга (его и правим); fpost — поле ПОСЛЕ столкновения (до стриминга)."""
|
||||
for (i, ib, mask, q, xff, cl) in bc:
|
||||
fi = fpost[i]; fib = fpost[ib]
|
||||
fiback = cp.roll(fpost[i], (Cy[i], Cx[i]), (0, 1)) # популяция «дальнего» соседа
|
||||
near = mask & (q < 0.5) & xff
|
||||
bad = mask & (q < 0.5) & (~xff) # нет дальнего соседа → простой отскок
|
||||
far = mask & (q >= 0.5)
|
||||
new = f[ib]
|
||||
new = xp.where(near, 2*q*fi + (1 - 2*q)*fiback, new)
|
||||
new = xp.where(far, (1/(2*q))*fi + (1 - 1/(2*q))*fib, new)
|
||||
new = xp.where(bad, fi, new)
|
||||
f[ib] = new
|
||||
return f
|
||||
|
||||
def channel_bc(f, uin, rho_out=1.0, outlet_uy="zero"):
|
||||
"""Корректные ГУ канала (Zou & He 1997), нумерация D2Q9: 1=E,2=N,3=W,4=S,5=NE,6=NW,7=SW,8=SE.
|
||||
ВХОД (запад, x=0): скоростной Zou-He — задаём u=(ux,uy) из uin, плотность ρ ПЛАВАЕТ; достраиваем
|
||||
приходящие извне популяции 1,5,8 из известных 0,2,3,4,6,7.
|
||||
ВЫХОД (восток, x=-1): давление Zou-He — задаём ρ=rho_out, скорость ux плавает; достраиваем 3,6,7.
|
||||
Поперечная скорость на выходе — по outlet_uy:
|
||||
"zero" — классический Zou-He: uy=0 (жёстко). Вихри дорожки приходят на выход с
|
||||
uy~±0.3U → принудительное зануление отражает их назад к телу (завышает
|
||||
rms Cl и модулирует Cd — см. факторное исследование в ноутбуке);
|
||||
"extrapolate" — uy берётся из соседнего столбца x=-2 (нуль-градиент поперечной скорости);
|
||||
ρ-якорь (масса) сохраняется, кинематическое отражение снимается.
|
||||
Зачем давление-выход вообще: грубый zero-gradient по f не выпускал поток → масса/противодавление
|
||||
копились, поток глох (⟨ρ⟩↑, Cd/rmsCl→0). Давление-выход якорит массу, скоростной вход её не
|
||||
пересоздаёт → дрейфа нет. Всё на срезах (graph-safe, без host→device).
|
||||
|
||||
ПОРЯДОК ВЫЗОВА: строго ПОСЛЕ free_slip_walls и покрывая ВЕСЬ столбец (включая ряды 0 и Ny−1).
|
||||
Причина: stream() сделан через roll и потому периодичен по ОБЕИМ осям. Заворот по y снимают
|
||||
стенки, заворот по x — этот ГУ. Если оставить углы (ряды 0 и Ny−1) без Zou-He, в них популяции
|
||||
с cx=+1 на входе приходят прямо из столбца выхода (замерено: f[1,0,0] == post[1,0,-1], величина
|
||||
~0.12) — вход и выход оказываются физически связаны в 4 узлах. Известные популяции, по которым
|
||||
здесь считаются ρ (вход) и ux (выход), x-заворота не содержат ни в одном ряду, а те, что
|
||||
заворот содержат, — это ровно неизвестные, которые ГУ и перезаписывает; поэтому такой порядок
|
||||
закрывает углы полностью."""
|
||||
# --- вход (запад): скоростной Zou-He ---
|
||||
ux = uin[0, :, 0]; uy = uin[1, :, 0]
|
||||
f0 = f[0, :, 0]; f2 = f[2, :, 0]; f3 = f[3, :, 0]
|
||||
f4 = f[4, :, 0]; f6 = f[6, :, 0]; f7 = f[7, :, 0]
|
||||
rho = (f0 + f2 + f4 + 2.0*(f3 + f6 + f7)) / (1.0 - ux)
|
||||
f[1, :, 0] = f3 + (2.0/3.0)*rho*ux
|
||||
f[5, :, 0] = f7 - 0.5*(f2 - f4) + (1.0/6.0)*rho*ux + 0.5*rho*uy
|
||||
f[8, :, 0] = f6 + 0.5*(f2 - f4) + (1.0/6.0)*rho*ux - 0.5*rho*uy
|
||||
# --- выход (восток): давление Zou-He ---
|
||||
g0 = f[0, :, -1]; g1 = f[1, :, -1]; g2 = f[2, :, -1]
|
||||
g4 = f[4, :, -1]; g5 = f[5, :, -1]; g8 = f[8, :, -1]
|
||||
uxo = (g0 + g2 + g4 + 2.0*(g1 + g5 + g8)) / rho_out - 1.0
|
||||
if outlet_uy == "extrapolate":
|
||||
# нуль-градиент поперечной скорости: uy с предвыходного столбца (он полностью известен)
|
||||
col = f[:, :, -2]
|
||||
uyo = ((col[2] + col[5] + col[6]) - (col[4] + col[7] + col[8])) / col.sum(0)
|
||||
else:
|
||||
uyo = 0.0
|
||||
f[3, :, -1] = g1 - (2.0/3.0)*rho_out*uxo
|
||||
f[6, :, -1] = g8 - 0.5*(g2 - g4) - (1.0/6.0)*rho_out*uxo + 0.5*rho_out*uyo
|
||||
f[7, :, -1] = g5 + 0.5*(g2 - g4) - (1.0/6.0)*rho_out*uxo - 0.5*rho_out*uyo
|
||||
return f
|
||||
|
||||
def free_slip_walls(f):
|
||||
"""FREE-SLIP (скользящие) стенки канала: верх/низ — зеркальное (specular) отражение.
|
||||
Сохраняет касательный (x) импульс (нет трения, нет пограничного слоя) и нормальный заворачивает;
|
||||
масса сохраняется (копирование популяций) → дрейфа не вносит. Применять ПОСЛЕ стриминга.
|
||||
Нумерация: 2=(0,+1) 4=(0,-1) 5=(+1,+1) 6=(-1,+1) 7=(-1,-1) 8=(+1,-1).
|
||||
Нижняя стенка (row 0): достраиваем вверх-идущие из вниз-идущих (зеркало по y, x сохранён).
|
||||
Верхняя стенка (row -1): достраиваем вниз-идущие из вверх-идущих.
|
||||
ВНИМАНИЕ: это только искусственные стенки канала; цилиндр остаётся no-slip (Bouzidi).
|
||||
ПОРЯДОК: вызывать ПЕРЕД channel_bc (см. его docstring — так закрываются углы домена)."""
|
||||
f[2, 0, :] = f[4, 0, :]; f[5, 0, :] = f[8, 0, :]; f[6, 0, :] = f[7, 0, :]
|
||||
f[4, -1, :] = f[2, -1, :]; f[7, -1, :] = f[6, -1, :]; f[8, -1, :] = f[5, -1, :]
|
||||
return f
|
||||
@@ -0,0 +1,42 @@
|
||||
# Оператор СТОЛКНОВЕНИЯ — KBC-N1 (энтропийный), единственный оператор модели.
|
||||
#
|
||||
# KBC-N1 (Karlin/Bösch/Chikatamarla 2014):
|
||||
# Δ = f − feq (неравновесная часть)
|
||||
# Δs = Ps·Δ (проекция на сдвиг-моменты {cx²−cy², cxcy})
|
||||
# Δh = Δ − Δs (остальная часть)
|
||||
# γ = 1/β − (2 − 1/β)·⟨Δs,Δh⟩ / ⟨Δh,Δh⟩, ⟨a,b⟩ = Σ_i a_i b_i / feq_i (H-теорема)
|
||||
# f_post = f − β(2Δs + γΔh), β = 1/(2τ)
|
||||
# γ — энтропийный стабилизатор: подстраивает диссипацию неhydro-мод под энтропийное ограничение
|
||||
# (поэтому KBC устойчив при τ→1/2). Сдвиг-моменты (вязкость) релаксируют с правильным 2β = 1/τ,
|
||||
# то есть физическая вязкость задаётся только τ — стабилизатор её не меняет.
|
||||
#
|
||||
# Порядок шагов дословно совпадает с алгоритмом статьи (Bösch/Chikatamarla/Karlin, PRE 92,
|
||||
# 043309 (2015), разд. III, листинг после ур.39): ρ,u → feq → Δs → Δh = f − feq − Δs → γ → релаксация.
|
||||
# Вырожденный случай ⟨Δh|Δh⟩→0 отсекается ОТНОСИТЕЛЬНЫМ порогом B.GREL (см. backend.py):
|
||||
# абсолютный порог здесь недопустим — den квадратична по неравновесию и физически мала.
|
||||
import backend as B
|
||||
import lattice as L
|
||||
|
||||
cp = B.cp; xp = B.xp
|
||||
Q = L.Q; Ps = L.Ps; GREL = B.GREL; ONE = B.ONE; TWO = B.TWO
|
||||
|
||||
def _project_shift(dfn):
|
||||
"""Δs = Ps·Δ без cuBLAS (gemm при CUDA-graph capture тащит alpha/beta host→device).
|
||||
Реализовано как Σ_k Ps[i,k]·Δ[k] через broadcast+редукцию — только обычные ядра."""
|
||||
dfn2d = dfn.reshape(Q, -1) # (Q, N)
|
||||
return (Ps[:, :, None] * dfn2d[None]).sum(1).reshape(dfn.shape)
|
||||
|
||||
def kbc_collide(f, fe, beta):
|
||||
dfn = f - fe
|
||||
ds = _project_shift(dfn)
|
||||
dh = dfn - ds; inv = 1.0 / fe
|
||||
num = (ds*dh*inv).sum(0); den = (dh*dh*inv).sum(0)
|
||||
nrm = (dfn*dfn*inv).sum(0) # ‖Δ‖² — масштаб неравновесия узла
|
||||
ok = den > GREL*nrm # относительный порог, не абсолютный
|
||||
den_safe = xp.where(ok, den, ONE) # избегаем 0/0 (обе ветки where считаются)
|
||||
binv = 1.0 / beta
|
||||
g = xp.where(ok, binv - (2.0 - binv)*num/den_safe, TWO)
|
||||
return f - beta*(2.0*ds + g[None]*dh)
|
||||
|
||||
# единственный оператор столкновения модели
|
||||
collide = kbc_collide
|
||||
@@ -0,0 +1,138 @@
|
||||
# Конфигурация модели 2×+SDF. Всё в одном месте, чтобы математику/режимы было легко менять.
|
||||
# Единицы — решёточные (lu — длина воксель, ts — шаг по времени).
|
||||
import os
|
||||
from dataclasses import dataclass
|
||||
from typing import Optional
|
||||
|
||||
@dataclass
|
||||
class Config:
|
||||
# --- сетка / геометрия L0 (умеренно расширены: меньше блокировка + больше следа) ---
|
||||
Nx: int = 174 # домен по x (цилиндр на 40 → ~8.4D вниз по потоку до выхода)
|
||||
Ny: int = 90 # домен по y (стенки free-slip; блокировка β=D/Ny — варьируется серией)
|
||||
D: int = 16 # диаметр цилиндра на грубом уровне L0
|
||||
cx: int = 40
|
||||
cy: Optional[int] = None # центр по вертикали; по умолчанию Ny//2 (см. __post_init__)
|
||||
|
||||
# --- физика ---
|
||||
U: float = 0.07 # скорость набегающего потока
|
||||
Re: float = 150.0 # число Рейнольдса по D и U
|
||||
|
||||
# --- время ---
|
||||
steps: int = 100000
|
||||
ramp: int = 1000 # smoothstep-разгон входной скорости (гасит импульсный старт)
|
||||
frame_every: int = 34 # период сохранения кадров для гифки
|
||||
|
||||
# --- патч уровня L1 (в координатах L0); refine=2 → разрешение Δx/2 ---
|
||||
ax1: Optional[int] = None # по умолчанию cx − 1.5D (запас перед телом; см. __post_init__)
|
||||
bx1: Optional[int] = None # по умолчанию cx + 6.25D (расширен в след)
|
||||
ay1: Optional[int] = None # по умолчанию cy − 2D (центрирован по cy)
|
||||
by1: Optional[int] = None # по умолчанию cy + 2D
|
||||
refine: int = 2 # коэффициент измельчения L1 (модель «2×» рассчитана на 2)
|
||||
|
||||
# --- зонд скорости в следе: px = cx + probe_dx*D, py = cy (берётся с тонкой сетки L1) ---
|
||||
probe_dx: int = 3
|
||||
|
||||
# --- стартовое возмущение (триггер вихревой дорожки): поперечный sin-импульс входа ---
|
||||
pert_amp: float = 0.1 # амплитуда u_y входа, доля от U (0 — выключить)
|
||||
pert_dur: int = 300 # длительность импульса в шагах; окно [ramp, ramp+pert_dur]
|
||||
|
||||
# --- кросс-чек силы: интеграл тензора напряжений σ·n по окружности R1+δ на L1 (forces_stress.py) ---
|
||||
stress_every: int = 50 # период выборки в шагах L0 (0 — выключить кросс-чек)
|
||||
stress_ntheta: int = 256 # число точек контура по углу
|
||||
stress_delta: float = 2.0 # отступ контура от поверхности (в тонких ячейках L1);
|
||||
# δ=2 → билинейный стенсиль не трогает Bouzidi-зону (q_max≈0.98)
|
||||
|
||||
# --- ГУ выхода: режим поперечной скорости на Zou-He давление-выходе ---
|
||||
# "zero" — классика: uy=0 на выходе (жёстко; отражает вихри дорожки)
|
||||
# "extrapolate" — uy берётся из соседнего столбца (нуль-градиент); ρ-якорь сохраняется
|
||||
outlet_uy: str = "zero"
|
||||
|
||||
# --- губка (absorbing layer) перед выходом: плавный рост вязкости в последних sponge_len
|
||||
# столбцах L0 (smoothstep, ν → ν·sponge_nu_mult). Гасит вихри И акустику до прихода на
|
||||
# давление-выход. Зачем: канал «вход-скорость + выход-давление» — недодемпфированный
|
||||
# акустический резонатор (затухание моды ∝ ν(π/Nx)²); на длинных доменах мода раскачивается
|
||||
# стартовым транзиентом → смещённый режим ⟨ρ⟩≈1.66 или взрыв (факторное исследование,
|
||||
# эксперимент №1: F1/F2/F3/F5). 0 — выключена.
|
||||
sponge_len: int = 0
|
||||
sponge_nu_mult: float = 30.0
|
||||
|
||||
# --- исполнение ---
|
||||
use_cuda_graph: bool = True
|
||||
|
||||
def __post_init__(self):
|
||||
# производная геометрия по вертикали: центр и патч следуют за Ny (для серии по блокировке)
|
||||
if self.cy is None:
|
||||
self.cy = self.Ny // 2
|
||||
if self.ay1 is None:
|
||||
self.ay1 = self.cy - 2 * self.D
|
||||
if self.by1 is None:
|
||||
self.by1 = self.cy + 2 * self.D
|
||||
# производная геометрия по горизонтали: патч следует за cx (для факторов вход/выход)
|
||||
if self.ax1 is None:
|
||||
self.ax1 = self.cx - (3 * self.D) // 2
|
||||
if self.bx1 is None:
|
||||
self.bx1 = self.cx + (25 * self.D) // 4
|
||||
assert 1 <= self.ay1 < self.by1 <= self.Ny - 2, (
|
||||
f"патч L1 [{self.ay1},{self.by1}] должен лежать строго внутри (0,{self.Ny-1}) — "
|
||||
f"specular-ряды стенок не накрывать; минимально допустимое Ny ≈ {4*self.D + 6}")
|
||||
assert 2 <= self.ax1 < self.cx < self.bx1 <= self.Nx - 2, (
|
||||
f"патч L1 [{self.ax1},{self.bx1}] должен лежать строго внутри (1,{self.Nx-1}) "
|
||||
f"и накрывать тело (cx={self.cx})")
|
||||
assert self.outlet_uy in ("zero", "extrapolate"), self.outlet_uy
|
||||
if self.sponge_len:
|
||||
assert self.Nx - 1 - self.sponge_len > self.bx1, (
|
||||
f"губка (последние {self.sponge_len} столбцов) не должна накрывать патч L1 "
|
||||
f"(bx1={self.bx1}, Nx={self.Nx})")
|
||||
|
||||
# ---------- производные величины ----------
|
||||
@property
|
||||
def nu(self): return self.U * self.D / self.Re # кинематическая вязкость (lu²/ts)
|
||||
@property
|
||||
def tau0(self): return self.nu / (1.0/3.0) + 0.5 # время релаксации на L0 (cs²=1/3)
|
||||
@property
|
||||
def tau1(self): # связка τ при измельчении в r раз: ν одна и та же ⇒ τ_f = r(τ_c − ½) + ½
|
||||
return self.refine * (self.tau0 - 0.5) + 0.5
|
||||
@property
|
||||
def beta0(self): return 1.0 / (2.0 * self.tau0) # β = 1/(2τ) для KBC/BGK на L0
|
||||
@property
|
||||
def beta1(self): return 1.0 / (2.0 * self.tau1) # β на L1
|
||||
@property
|
||||
def R01(self): # масштаб неравновесной части грубый→тонкий: f^neq ∝ τ·δt ⇒ (τ_f/τ_c)/r
|
||||
return self.tau1 / (self.refine * self.tau0)
|
||||
@property
|
||||
def Rf01(self): return 1.0 / self.R01 # тонкий→грубый
|
||||
@property
|
||||
def DL(self): return self.refine * self.D # эффективный диаметр у тела (на L1)
|
||||
@property
|
||||
def px(self): return self.cx + self.probe_dx * self.D
|
||||
@property
|
||||
def py(self): return self.cy
|
||||
|
||||
def from_env(**overrides):
|
||||
"""Создать Config с учётом переменных окружения и явных overrides.
|
||||
Геометрия (AMR_NY/AMR_NX/AMR_CX) попадает в overrides ДО конструктора —
|
||||
производные cy/ay1/by1/ax1/bx1 строятся в __post_init__."""
|
||||
if "Ny" not in overrides and os.environ.get("AMR_NY"):
|
||||
overrides["Ny"] = int(os.environ["AMR_NY"])
|
||||
if "Nx" not in overrides and os.environ.get("AMR_NX"):
|
||||
overrides["Nx"] = int(os.environ["AMR_NX"])
|
||||
if "cx" not in overrides and os.environ.get("AMR_CX"):
|
||||
overrides["cx"] = int(os.environ["AMR_CX"])
|
||||
if "outlet_uy" not in overrides and os.environ.get("AMR_OUTLET_UY"):
|
||||
overrides["outlet_uy"] = os.environ["AMR_OUTLET_UY"]
|
||||
if "sponge_len" not in overrides and os.environ.get("AMR_SPONGE_LEN"):
|
||||
overrides["sponge_len"] = int(os.environ["AMR_SPONGE_LEN"])
|
||||
if "sponge_nu_mult" not in overrides and os.environ.get("AMR_SPONGE_NU"):
|
||||
overrides["sponge_nu_mult"] = float(os.environ["AMR_SPONGE_NU"])
|
||||
cfg = Config(**overrides)
|
||||
if os.environ.get("AMR_STEPS"):
|
||||
cfg.steps = int(os.environ["AMR_STEPS"])
|
||||
if "AMR_CUDAGRAPH" in os.environ:
|
||||
cfg.use_cuda_graph = (os.environ["AMR_CUDAGRAPH"] != "0")
|
||||
if os.environ.get("AMR_STRESS_EVERY"):
|
||||
cfg.stress_every = int(os.environ["AMR_STRESS_EVERY"])
|
||||
if os.environ.get("AMR_STRESS_NTHETA"):
|
||||
cfg.stress_ntheta = int(os.environ["AMR_STRESS_NTHETA"])
|
||||
if os.environ.get("AMR_STRESS_DELTA"):
|
||||
cfg.stress_delta = float(os.environ["AMR_STRESS_DELTA"])
|
||||
return cfg
|
||||
@@ -0,0 +1,140 @@
|
||||
# Диагностика по собранным рядам сил/скорости. numpy здесь — только host-постобработка малых рядов
|
||||
# (FFT/статистика), это НЕ вычислительный бэкенд солвера.
|
||||
import numpy as np
|
||||
import backend as B
|
||||
|
||||
def strouhal(uy_host, D, U, pad=8192):
|
||||
"""Число Струхаля из ряда поперечной скорости в следе.
|
||||
УЛУЧШЕНО: параболическая (3-точечная) интерполяция пика спектра → суб-биновая точность.
|
||||
(Раньше брался ближайший FFT-бин: шаг по St ≈ D/(U·pad) ≈ 0.028, из-за чего St 'застревал'.)"""
|
||||
sig = uy_host[len(uy_host)//2:] # установившийся режим — вторая половина
|
||||
sig = sig - sig.mean()
|
||||
amp = np.abs(np.fft.rfft(sig, n=pad))
|
||||
if len(amp) <= 3:
|
||||
return 0.0, np.fft.rfftfreq(pad, 1.0), amp
|
||||
k = 1 + int(np.argmax(amp[1:])) # индекс пика (без DC)
|
||||
if 1 <= k < len(amp) - 1: # параболическая интерполяция вершины
|
||||
a0, a1, a2 = amp[k-1], amp[k], amp[k+1]
|
||||
denom = (a0 - 2*a1 + a2)
|
||||
delta = 0.5*(a0 - a2)/denom if denom != 0 else 0.0
|
||||
else:
|
||||
delta = 0.0
|
||||
f_peak = (k + delta) / pad
|
||||
return f_peak * D / U, np.fft.rfftfreq(pad, 1.0), amp
|
||||
|
||||
def analyze(res, cfg):
|
||||
Fx = B.to_cpu(res["Fx"]); Fy = B.to_cpu(res["Fy"]); uy = B.to_cpu(res["uy"])
|
||||
DL, U, D = res["DL"], cfg.U, cfg.D
|
||||
n = len(uy); h = slice(n//2, n)
|
||||
Cd = 2.0*Fx/(U**2*DL); Cl = 2.0*Fy/(U**2*DL)
|
||||
St, freqs, amp = strouhal(uy, D, U)
|
||||
a = dict(Cd=float(Cd[h].mean()), Clrms=float(Cl[h].std()), St=float(St),
|
||||
uyrms=float(uy[h].std()), Cd_s=Cd, Cl_s=Cl, uy=uy, freqs=freqs, amp=amp, h0=n//2)
|
||||
if "Tz" in res: # момент: Cm = 2Tz/(U²·DL²); для цилиндра ⟨Cm⟩≈0
|
||||
Cm = 2.0*B.to_cpu(res["Tz"])/(U**2*DL*DL)
|
||||
a["Cm"] = float(Cm[h].mean()); a["Cmrms"] = float(Cm[h].std()); a["Cm_s"] = Cm
|
||||
K = int(res.get("stress_every", 0) or 0) # кросс-чек σ·n (редкая сетка по времени)
|
||||
if K and len(res.get("Fx_st", ())) > 3:
|
||||
Fxs = B.to_cpu(res["Fx_st"]); Fys = B.to_cpu(res["Fy_st"]); Tzs = B.to_cpu(res["Tz_st"])
|
||||
m = len(Fxs); hs = slice(m//2, m)
|
||||
Cds = 2.0*Fxs/(U**2*DL); Cls = 2.0*Fys/(U**2*DL); Cms = 2.0*Tzs/(U**2*DL*DL)
|
||||
a.update(Cd_st=float(Cds[hs].mean()), Clrms_st=float(Cls[hs].std()),
|
||||
Cm_st=float(Cms[hs].mean()), Cd_st_s=Cds, Cl_st_s=Cls, stress_every=K)
|
||||
return a
|
||||
|
||||
def print_table(a, cfg):
|
||||
# поправка на блокировку канала (приведение к скорости в зазоре U/(1−β)):
|
||||
# St — кинематическая → ×(1−β); Cd, rms Cl ~ скорость² → ×(1−β)²
|
||||
beta = cfg.D / cfg.Ny
|
||||
St_c = a['St'] * (1 - beta)
|
||||
Cd_c = a['Cd'] * (1 - beta)**2
|
||||
Cl_c = a['Clrms'] * (1 - beta)**2
|
||||
cm = a.get('Cm')
|
||||
cms = f"{cm:>10.5f}" if cm is not None else f"{'—':>10}"
|
||||
print("\n=== Модель 2×+SDF — установившийся режим ===")
|
||||
print(f"{'версия':18}{'D у тела':>9}{'St':>8}{'<Cd>':>9}{'rms Cl':>9}{'<Cm>':>10}{'rms u_y':>9}")
|
||||
print(f"{'2×+SDF (raw)':18}{cfg.DL:>9}{a['St']:>8.3f}{a['Cd']:>9.3f}{a['Clrms']:>9.4f}{cms}{a['uyrms']:>9.4f}")
|
||||
print(f"{'2×+SDF (×блок.)':18}{cfg.DL:>9}{St_c:>8.3f}{Cd_c:>9.3f}{Cl_c:>9.4f}{'—':>10}{'—':>9}")
|
||||
print(f"{'лит. Re≈150':18}{'—':>9}{0.183:>8.3f}{1.330:>9.3f}{'~0.3':>9}{0.0:>10.1f}{'':>9}")
|
||||
print(f"(блокировка β=D/Ny={beta:.3f}; поправка: St×(1−β), Cd/rmsCl×(1−β)². "
|
||||
"Лит. — безграничный цилиндр; ⟨Cm⟩≈0 — контроль симметрии; см. README)")
|
||||
|
||||
def print_crosscheck(a):
|
||||
"""Сравнение двух НЕЗАВИСИМЫХ считываний силы: GMEM (обмен импульсом по линкам) vs ∮σ·n ds
|
||||
(интеграл тензора напряжений по контуру R+δ). Сходство — прямое свидетельство корректности
|
||||
считывания; сравниваем средние (мгновенные ряды σ·n сдвинуты по фазе кольцом R..R+δ)."""
|
||||
if "Cd_st" not in a:
|
||||
print("\n[кросс-чек σ·n: выключен (stress_every=0) или слишком мало точек]")
|
||||
return
|
||||
dcd = abs(a["Cd_st"] - a["Cd"]) / max(abs(a["Cd"]), 1e-12)
|
||||
dcl = abs(a["Clrms_st"] - a["Clrms"]) / max(abs(a["Clrms"]), 1e-12)
|
||||
cm = a.get("Cm", float("nan"))
|
||||
print("\n=== Кросс-чек считывания силы: GMEM vs ∮σ·n ds (установившийся режим) ===")
|
||||
print(f"{'метод':26}{'<Cd>':>9}{'rms Cl':>9}{'<Cm>':>10}")
|
||||
print(f"{'GMEM (обмен импульсом)':26}{a['Cd']:>9.3f}{a['Clrms']:>9.4f}{cm:>10.5f}")
|
||||
print(f"{'∮σ·n ds (контур R+δ)':26}{a['Cd_st']:>9.3f}{a['Clrms_st']:>9.4f}{a['Cm_st']:>10.5f}")
|
||||
verdict = "СОГЛАСОВАНО" if (dcd < 0.05 and dcl < 0.10) else "ПРОВЕРИТЬ"
|
||||
print(f"расхождение: dCd={dcd*100:.1f}% (порог 5%), d(rmsCl)={dcl*100:.1f}% (порог 10%) → {verdict}")
|
||||
|
||||
def fit_blockage(betas, vals):
|
||||
"""Экстраполяция величины серии к β→0 двумя моделями: val ≈ a + b·β и val ≈ a + b·β².
|
||||
(Теория solid-blockage даёт ведущий член O(β²), практика поправки (1−β)² содержит и линейный.)
|
||||
Спред интерсептов |a_lin − a_quad| — оценка модельной неопределённости экстраполяции."""
|
||||
b = np.asarray(betas, float); v = np.asarray(vals, float)
|
||||
sl, al = np.polyfit(b, v, 1)
|
||||
sq, aq = np.polyfit(b**2, v, 1)
|
||||
return dict(lin=(float(al), float(sl)), quad=(float(aq), float(sq)),
|
||||
spread=float(abs(al - aq)))
|
||||
|
||||
def print_blockage_table(rows, fits):
|
||||
"""Сводка серии по блокировке. rows — список dict(Ny, beta, St, Cd, Clrms, Cm, Cd_st);
|
||||
fits — dict('Cd'/'Clrms'/'St' → результат fit_blockage)."""
|
||||
print("\n=== Серия по блокировке: D фикс., Ny варьируется ===")
|
||||
print(f"{'Ny':>5}{'β=D/Ny':>9}{'St':>8}{'St(1−β)':>9}{'<Cd>':>9}{'rms Cl':>9}{'<Cm>':>10}{'Cd σ·n':>9}")
|
||||
for r in rows:
|
||||
cdst = f"{r['Cd_st']:>9.3f}" if r.get("Cd_st") is not None else f"{'—':>9}"
|
||||
print(f"{r['Ny']:>5}{r['beta']:>9.3f}{r['St']:>8.3f}{r['St']*(1-r['beta']):>9.3f}"
|
||||
f"{r['Cd']:>9.3f}{r['Clrms']:>9.4f}{r['Cm']:>10.5f}{cdst}")
|
||||
print("экстраполяция β→0 (интерсепт лин. фита / квадр. фита; спред — неопределённость):")
|
||||
print(f" Cd(0) = {fits['Cd']['lin'][0]:.3f} / {fits['Cd']['quad'][0]:.3f} "
|
||||
f"(спред {fits['Cd']['spread']:.3f}; лит. безгранич. 1.33)")
|
||||
print(f" rms Cl(0) = {fits['Clrms']['lin'][0]:.3f} / {fits['Clrms']['quad'][0]:.3f} "
|
||||
f"(спред {fits['Clrms']['spread']:.3f}; лит. ~0.3)")
|
||||
print(f" St(0) = {fits['St']['lin'][0]:.3f} / {fits['St']['quad'][0]:.3f} "
|
||||
f"(спред {fits['St']['spread']:.3f}; лит. 0.183)")
|
||||
|
||||
def convergence_report(res, cfg, nwin=10):
|
||||
"""Эволюция по окнам времени: видно, НАСЫЩЕНО ли решение или ДРЕЙФУЕТ.
|
||||
Ключевое: ⟨ρ⟩ (средняя плотность) и Cd, нормированный на реальную ⟨ρ⟩ (дрейф-устойчивый Cd).
|
||||
Если ⟨ρ⟩ уплывает от 1 — это дрейф массы из-за ГУ входа/выхода, а Cd_raw искажён нормировкой на ρ=1."""
|
||||
Fx = B.to_cpu(res["Fx"]); Fy = B.to_cpu(res["Fy"]); uy = B.to_cpu(res["uy"]); rho = B.to_cpu(res["rho"])
|
||||
n = len(Fx); U = cfg.U; DL = res["DL"]
|
||||
if n < nwin * 2:
|
||||
return
|
||||
Cd = 2.0*Fx/(U**2*DL); Cl = 2.0*Fy/(U**2*DL)
|
||||
w = n // nwin
|
||||
print("\n=== Сходимость по времени / дрейф (по окнам) ===")
|
||||
print(f"{'окно':>4}{'шаги':>16}{'<ρ>':>8}{'<Cd>raw':>9}{'<Cd>/ρ':>9}{'rms Cl':>9}{'rms u_y':>9}")
|
||||
rows = []
|
||||
for k in range(nwin):
|
||||
s = slice(k*w, (k+1)*w if k < nwin-1 else n)
|
||||
rm = float(rho[s].mean()); cdr = float(Cd[s].mean())
|
||||
cdn = cdr / rm if rm != 0 else float("nan")
|
||||
clr = float(Cl[s].std()); uyr = float(uy[s].std())
|
||||
rows.append((rm, cdr, cdn, clr, uyr))
|
||||
print(f"{k+1:>4}{f'{k*w}-{(k+1)*w}':>16}{rm:>8.3f}{cdr:>9.3f}{cdn:>9.3f}{clr:>9.4f}{uyr:>9.4f}")
|
||||
# вердикт по последним двум окнам
|
||||
(rm1, cdr1, cdn1, clr1, _), (rm0, _, _, clr0, _) = rows[-1], rows[-2]
|
||||
drho = rm1 - rm0; dcl = clr1 - clr0
|
||||
verdict = []
|
||||
if abs(rm1 - 1.0) > 0.02:
|
||||
# «стабильна, но не там»: система села на смещённую ветвь (акустическая накачка массы;
|
||||
# ⟨Cd⟩/Cl на ней НЕвалидны для сравнения с литературой) — см. факторное исследование
|
||||
verdict.append(f"плотность: ⟨ρ⟩={rm1:.3f} → РЕЖИМ СМЕЩЁН (⟨ρ⟩ далеко от 1; прогон невалиден)")
|
||||
else:
|
||||
verdict.append(f"плотность: ⟨ρ⟩={rm1:.3f}, Δ за окно={drho:+.4f} → "
|
||||
+ ("ДРЕЙФ массы (правь ГУ выхода: давление/Zou-He)" if abs(drho) > 1e-3 else "стабильна"))
|
||||
verdict.append(f"rms Cl: Δ за окно={dcl:+.4f} → "
|
||||
+ ("ещё растёт (не насыщено)" if dcl > 1e-3 else "насыщено/стабильно"))
|
||||
verdict.append(f"дрейф-устойчивый Cd (⟨Cd⟩/ρ) в последнем окне = {cdn1:.3f}")
|
||||
print("вердикт: " + "; ".join(verdict))
|
||||
@@ -0,0 +1,45 @@
|
||||
# Равновесие и макропеременные.
|
||||
#
|
||||
# Используется ЭНТРОПИЙНОЕ равновесие в product-form (как в KBC), а не усечённый полином:
|
||||
# feq_i = w_i · ρ · Π_d (2 − sqrt(1+3 u_d²)) · ((2 u_d + sqrt(1+3 u_d²))/(1 − u_d))^{c_{i,d}}
|
||||
# Показатели c_{i,d} ∈ {−1,0,+1}, поэтому степень НЕ нужна — берём q или 1/q (быстрее и совместимо
|
||||
# с CUDA-graph: нет pow и нет host→device для границ clip).
|
||||
import backend as B
|
||||
import lattice as L
|
||||
|
||||
cp = B.cp; xp = B.xp
|
||||
CXa = L.CXa; CYa = L.CYa
|
||||
|
||||
def _feq9_body(rho, ux, uy):
|
||||
ux = cp.minimum(cp.maximum(ux, -0.95), 0.95) # clip скоростей литералами (запекается в кернел)
|
||||
uy = cp.minimum(cp.maximum(uy, -0.95), 0.95)
|
||||
sx = cp.sqrt(1.0 + 3.0*ux*ux); sy = cp.sqrt(1.0 + 3.0*uy*uy)
|
||||
base = rho * (2.0 - sx) * (2.0 - sy)
|
||||
qx = (2.0*ux + sx) / (1.0 - ux); qy = (2.0*uy + sy) / (1.0 - uy)
|
||||
ix = 1.0/qx; iy = 1.0/qy
|
||||
return (base*(4.0/9.0), # ( 0, 0)
|
||||
base*(1.0/9.0)*qx, base*(1.0/9.0)*qy, # (+1,0) (0,+1)
|
||||
base*(1.0/9.0)*ix, base*(1.0/9.0)*iy, # (-1,0) (0,-1)
|
||||
base*(1.0/36.0)*qx*qy, base*(1.0/36.0)*ix*qy, # (+1,+1) (-1,+1)
|
||||
base*(1.0/36.0)*ix*iy, base*(1.0/36.0)*qx*iy) # (-1,-1) (+1,-1)
|
||||
|
||||
_feq9 = cp.fuse()(_feq9_body) # одно слитное ядро (если cp.fuse доступен)
|
||||
_FUSE_OK = [True]
|
||||
|
||||
def feq(rho, u):
|
||||
"""Равновесные популяции (Q,Ny,Nx) по плотности rho и скорости u (индексируется [0],[1]).
|
||||
Сборка через B.pack (без cupy.stack) — иначе H2D ломает CUDA-graph capture."""
|
||||
if _FUSE_OK[0]:
|
||||
try:
|
||||
return B.pack(_feq9(rho, u[0], u[1]))
|
||||
except Exception as e: # откат на не-fused (тоже без pow), один раз
|
||||
print(f"[equilibrium] cp.fuse отключён ({type(e).__name__}: {e}); обычный путь")
|
||||
_FUSE_OK[0] = False
|
||||
return B.pack(_feq9_body(rho, u[0], u[1]))
|
||||
|
||||
def macros(f):
|
||||
"""Макропеременные: ρ = Σ_i f_i, u = (Σ_i c_i f_i)/ρ. Возвращает (rho, u=(2,Ny,Nx))."""
|
||||
r = f.sum(0)
|
||||
ux = (CXa[:, None, None] * f).sum(0) / r
|
||||
uy = (CYa[:, None, None] * f).sum(0) / r
|
||||
return r, B.pack((ux, uy))
|
||||
@@ -0,0 +1,52 @@
|
||||
# СИЛА и МОМЕНТ на теле — галилей-инвариантный обмен импульсом (GMEM, Wen et al. 2014).
|
||||
#
|
||||
# Для каждого линка жидкость→тело по направлению i (маска cl):
|
||||
# ΔF = (c_i − u_w)·f_i^{после столкновения}(x_f) − (c_ī − u_w)·f_ī^{после стриминга}(x_f)
|
||||
# где u_w — скорость стенки на линке (для неподвижного цилиндра = 0). Т.к. c_ī = −c_i, по компонентам:
|
||||
# F_x += (Cx_i − u_wx)·f_i + (Cx_i + u_wx)·f_ī
|
||||
# F_y += (Cy_i − u_wy)·f_i + (Cy_i + u_wy)·f_ī
|
||||
#
|
||||
# МОМЕНТ (торк) вокруг центра тела: T_z = Σ_links (r_w × ΔF)_z, где точка приложения —
|
||||
# ПЕРЕСЕЧЕНИЕ линка со стенкой r_w = r_f + q·c_i (q — Bouzidi-доля, уже лежит в bc).
|
||||
# Узел жидкости дал бы ошибку плеча до 1 ячейки; q бесплатен. Для кругового цилиндра
|
||||
# ⟨Cm⟩ ≈ 0 — это контроль симметрии считывания; для подвижных/вращающихся тел (цель SimV4) —
|
||||
# рабочая величина.
|
||||
#
|
||||
# Свойства:
|
||||
# * f_i — популяция, уходящая в стенку (после столкновения); f_ī — вернувшаяся (после стриминга/Bouzidi).
|
||||
# f_ī берётся из поля ПОСЛЕ стриминга — иначе ведущий симметричный член ~2ρw_i сокращается и Cd врёт.
|
||||
# * При u_w=0 сводится к классическому Σ c_i (f_i + f_ī); члены с u_w делают оценку галилей-инвариантной
|
||||
# и ГОТОВОЙ к подвижным/вращающимся телам (цель проекта SimV4) — туда передаётся u_w на линке.
|
||||
#
|
||||
# Независимый кросс-чек этой оценки — интегрирование тензора напряжений по контуру вокруг тела:
|
||||
# см. forces_stress.py (σ = −p′·I + σ^v из неравновесных моментов; вызывается редко, вне CUDA-графа).
|
||||
import backend as B
|
||||
import lattice as L
|
||||
|
||||
xp = B.xp; ZERO = B.ZERO
|
||||
Cx = L.Cx; Cy = L.Cy
|
||||
|
||||
def force_gmem(fpost, fnew, bc, rxy=None, u_wall=(0.0, 0.0)):
|
||||
"""GMEM по линкам тела (маска cl). Возвращает (Fx, Fy, Tz) — device-скаляры.
|
||||
rxy=(rx,ry) — координаты узлов решётки относительно центра тела (для торка); None → Tz=0.
|
||||
u_wall — скорость стенки. Graph-safe: только where/арифметика с литералами."""
|
||||
uwx, uwy = u_wall
|
||||
Fx = ZERO; Fy = ZERO; Tz = ZERO
|
||||
if rxy is not None:
|
||||
rx, ry = rxy
|
||||
for (i, ib, mask, q, xff, cl) in bc:
|
||||
fi = fpost[i]; fib = fnew[ib]
|
||||
dFx = xp.where(cl, (Cx[i] - uwx)*fi + (Cx[i] + uwx)*fib, ZERO)
|
||||
dFy = xp.where(cl, (Cy[i] - uwy)*fi + (Cy[i] + uwy)*fib, ZERO)
|
||||
Fx = Fx + dFx.sum(); Fy = Fy + dFy.sum()
|
||||
if rxy is not None: # плечо до точки пересечения линка со стенкой
|
||||
Tz = Tz + ((rx + q*Cx[i])*dFy - (ry + q*Cy[i])*dFx).sum()
|
||||
return Fx, Fy, Tz
|
||||
|
||||
def coefficients(Fx, Fy, U, DL, Tz=None):
|
||||
"""Безразмерные коэффициенты: Cd = 2Fx/(U²·DL), Cl = 2Fy/(U²·DL), Cm = 2Tz/(U²·DL²).
|
||||
DL — диаметр у тела (на L1); сила/плечо тоже в единицах L1, так что нормировка согласована."""
|
||||
Cd = 2.0*Fx/(U**2*DL); Cl = 2.0*Fy/(U**2*DL)
|
||||
if Tz is None:
|
||||
return Cd, Cl
|
||||
return Cd, Cl, 2.0*Tz/(U**2*DL*DL)
|
||||
@@ -0,0 +1,80 @@
|
||||
# НЕЗАВИСИМЫЙ КРОСС-ЧЕК силы/момента: баланс импульса контрольного объёма, ограниченного
|
||||
# замкнутым контуром (окружностью) вокруг тела на тонком уровне L1.
|
||||
#
|
||||
# F_тела = ∮ [σ·n − ρ·u·(u·n)] ds − d/dt ∫_кольцо ρu dV, σ = −p′·I + σ^v
|
||||
#
|
||||
# * давление: p = ρ·c_s²; интегрируем p′ = (ρ − ⟨ρ⟩_контур)·c_s² — константа давления по замкнутому
|
||||
# контуру даёт ∮ p0·n ds ≡ 0, а вычитание убирает катастрофическое сокращение в float32;
|
||||
# * вязкая часть из неравновесных моментов (Чепмен–Энског): σ^v_αβ = −(1 − 1/(2τ))·Π^neq_αβ,
|
||||
# Π^neq_αβ = Σ_i c_iα c_iβ (f_i − f_i^eq).
|
||||
# ВАЖНО (KBC): след Π^neq релаксирует не с 1/τ, а с γβ (поле-зависимый стабилизатор), поэтому
|
||||
# берём строго ДЕВИАТОРНУЮ проекцию (след исключён) — сдвиг-моменты в KBC релаксируют с точным
|
||||
# 2β = 1/τ, для них префактор корректен. Bulk-часть в силу трения и не входит.
|
||||
# * КОНВЕКТИВНЫЙ член −ρu(u·n): поток импульса через контур. Зануляется только при δ→0 (no-slip
|
||||
# на поверхности); при δ=2 скорость на контуре уже заметна, и без этого члена кросс-чек давал
|
||||
# систематический сдвиг +4% по ⟨Cd⟩ на всех сетках (факторное исследование, эксперимент №0).
|
||||
# * член d/dt ∫ρu dV по кольцу R..R1+δ НЕ считается: в среднем по времени он зануляется,
|
||||
# на мгновенных рядах даёт малый сдвиг фазы → кросс-чек ведётся по средним (⟨Cd⟩, rms Cl).
|
||||
# * контур — окружность радиуса R1+δ (δ = cfg.stress_delta тонких ячеек, по умолчанию 2):
|
||||
# билинейный стенсиль выборки не задевает первое жидкое кольцо у стенки, где популяции
|
||||
# реконструированы Bouzidi (q_max ≈ 0.98 < 1).
|
||||
#
|
||||
# Вызывается РЕДКО (каждые cfg.stress_every шагов) в Solver.run() ВНЕ step()/CUDA-графа.
|
||||
import math
|
||||
import backend as B
|
||||
import lattice as L
|
||||
import equilibrium as E
|
||||
|
||||
cp = B.cp; xp = B.xp; DTYPE = B.DTYPE
|
||||
CS2 = L.CS2
|
||||
|
||||
def build_probe(cfg):
|
||||
"""Прекомпьют контура (один раз, вне захвата): N=stress_ntheta точек на окружности R1+δ
|
||||
вокруг центра тела на L1; индексы/веса билинейной интерполяции, нормали, плечи, элемент дуги."""
|
||||
r = cfg.refine
|
||||
R1 = r * float(cfg.D) / 2.0 # радиус цилиндра на L1 (тонкие ячейки)
|
||||
Rp = R1 + float(cfg.stress_delta)
|
||||
ccx = (cfg.cx - cfg.ax1) * r; ccy = (cfg.cy - cfg.ay1) * r
|
||||
N = int(cfg.stress_ntheta)
|
||||
th = (2.0 * math.pi / N) * cp.arange(N)
|
||||
nx = cp.cos(th).astype(DTYPE); ny = cp.sin(th).astype(DTYPE)
|
||||
px = ccx + Rp * nx; py = ccy + Rp * ny
|
||||
x0 = cp.floor(px).astype(cp.int64); y0 = cp.floor(py).astype(cp.int64)
|
||||
# контур обязан лежать внутри патча L1 (он центрирован вокруг тела — выполняется с запасом)
|
||||
Nfx = (cfg.bx1 - cfg.ax1) * r + 1; Nfy = (cfg.by1 - cfg.ay1) * r + 1
|
||||
assert int(x0.min()) >= 0 and int(x0.max()) + 1 <= Nfx - 1, "контур σ·n вышел за патч L1 по x"
|
||||
assert int(y0.min()) >= 0 and int(y0.max()) + 1 <= Nfy - 1, "контур σ·n вышел за патч L1 по y"
|
||||
return dict(x0=x0, y0=y0, x1=x0 + 1, y1=y0 + 1,
|
||||
tx=(px - x0).astype(DTYPE), ty=(py - y0).astype(DTYPE),
|
||||
nx=nx, ny=ny, rx=(Rp * nx), ry=(Rp * ny), ds=2.0 * math.pi * Rp / N)
|
||||
|
||||
def _sample(fld, pr):
|
||||
"""Билинейная выборка 2D-поля в точках контура (аналог amr.pint, но по 1D-набору точек)."""
|
||||
return (fld[pr["y0"], pr["x0"]] * (1 - pr["tx"]) * (1 - pr["ty"])
|
||||
+ fld[pr["y0"], pr["x1"]] * pr["tx"] * (1 - pr["ty"])
|
||||
+ fld[pr["y1"], pr["x0"]] * (1 - pr["tx"]) * pr["ty"]
|
||||
+ fld[pr["y1"], pr["x1"]] * pr["tx"] * pr["ty"])
|
||||
|
||||
def force_stress(F1, pr, tau1):
|
||||
"""Сила/момент на тело из баланса импульса ∮[σ·n − ρu(u·n)]ds по контуру pr на состоянии F1
|
||||
(L1, после стриминга). Возвращает (Fx, Fy, Tz) — device-скаляры. tau1 — релаксация на L1."""
|
||||
rho, u = E.macros(F1)
|
||||
fneq = F1 - E.feq(rho, u)
|
||||
CX = L.CXa[:, None, None]; CY = L.CYa[:, None, None]
|
||||
Pxx = (CX * CX * fneq).sum(0)
|
||||
Pyy = (CY * CY * fneq).sum(0)
|
||||
Pxy = (CX * CY * fneq).sum(0)
|
||||
rho_s = _sample(rho, pr)
|
||||
ux_s = _sample(u[0], pr); uy_s = _sample(u[1], pr)
|
||||
Pxx_s = _sample(Pxx, pr); Pyy_s = _sample(Pyy, pr); Pxy_s = _sample(Pxy, pr)
|
||||
p = (rho_s - rho_s.mean()) * CS2 # p′: давление без константы (см. шапку)
|
||||
pref = 1.0 - 1.0 / (2.0 * tau1)
|
||||
tr = 0.5 * (Pxx_s + Pyy_s) # девиаторная проекция (след исключён)
|
||||
sxx = -pref * (Pxx_s - tr); syy = -pref * (Pyy_s - tr); sxy = -pref * Pxy_s
|
||||
un = ux_s * pr["nx"] + uy_s * pr["ny"] # нормальная скорость на контуре
|
||||
flux = rho_s * un # ρ(u·n): конвективный вынос импульса
|
||||
tx = (-p + sxx) * pr["nx"] + sxy * pr["ny"] - flux * ux_s # t = σ·n − ρu(u·n)
|
||||
ty = sxy * pr["nx"] + (-p + syy) * pr["ny"] - flux * uy_s
|
||||
Fx = tx.sum() * pr["ds"]; Fy = ty.sum() * pr["ds"]
|
||||
Tz = (pr["rx"] * ty - pr["ry"] * tx).sum() * pr["ds"] # давление в торк не входит (n ∥ r)
|
||||
return Fx, Fy, Tz
|
||||
@@ -0,0 +1,47 @@
|
||||
# Геометрия: маски тела/стенок и знаковое поле расстояний (SDF) для L0 и патча L1.
|
||||
# Чистый CuPy. Соглашение по SDF: phi>0 в жидкости, phi<0 в теле; |phi| ≈ расстоянию до границы.
|
||||
import backend as B
|
||||
import lattice as L
|
||||
import amr
|
||||
|
||||
cp = B.cp; xp = B.xp; DTYPE = B.DTYPE
|
||||
|
||||
def _cmask(NX, NY, ccx, ccy, R):
|
||||
"""Булева маска круга радиуса R (центр ccx,ccy) на сетке (NY,NX)."""
|
||||
Y, X = cp.meshgrid(cp.arange(NY), cp.arange(NX), indexing="ij")
|
||||
return (X - ccx)**2 + (Y - ccy)**2 <= R*R
|
||||
|
||||
def _csdf(NX, NY, ccx, ccy, R):
|
||||
"""SDF круга: sqrt((x−ccx)²+(y−ccy)²) − R (отрицательно внутри)."""
|
||||
Y, X = cp.meshgrid(cp.arange(NY), cp.arange(NX), indexing="ij")
|
||||
return cp.sqrt((X - ccx).astype(DTYPE)**2 + (Y - ccy).astype(DTYPE)**2) - R
|
||||
|
||||
class Geometry:
|
||||
"""Контейнер масок/SDF/патча. Поля *_np — host-копии для визуализации."""
|
||||
pass
|
||||
|
||||
def build(cfg):
|
||||
Nx, Ny, D, cx, cy = cfg.Nx, cfg.Ny, cfg.D, cfg.cx, cfg.cy
|
||||
r = cfg.refine
|
||||
ax1, bx1, ay1, by1 = cfg.ax1, cfg.bx1, cfg.ay1, cfg.by1
|
||||
g = Geometry()
|
||||
|
||||
# --- L0: только цилиндр как твёрдое тело (стенки канала верх/низ — free-slip, не solid) ---
|
||||
g.cyl0 = _cmask(Nx, Ny, cx, cy, D/2.0)
|
||||
g.solid0 = g.cyl0 # стенки обрабатываются specular-отражением, см. boundary.free_slip_walls
|
||||
g.phi0 = _csdf(Nx, Ny, cx, cy, D/2.0) # SDF только цилиндра (для Bouzidi тела)
|
||||
|
||||
# --- патч L1 (измельчение r) над областью [ax1,bx1]×[ay1,by1] ---
|
||||
g.P1 = amr.patch(Nx, Ny, ax1, bx1, ay1, by1, r)
|
||||
Nfx, Nfy = g.P1["Nfx"], g.P1["Nfy"]
|
||||
# цилиндр на L1: центр и радиус масштабируются на r (диаметр r·D → радиус r·D/2)
|
||||
R1 = r * D / 2.0
|
||||
g.solid1 = _cmask(Nfx, Nfy, (cx - ax1)*r, (cy - ay1)*r, R1)
|
||||
g.phi1 = _csdf(Nfx, Nfy, (cx - ax1)*r, (cy - ay1)*r, R1)
|
||||
|
||||
# маска жидких узлов для рестрикции L1→L0 (внутренняя часть патча на L0)
|
||||
g.fl1 = ~g.solid0[ay1+1:by1, ax1+1:bx1]
|
||||
|
||||
# host-копии для визуализации
|
||||
g.solid0_np = B.to_cpu(g.solid0); g.solid1_np = B.to_cpu(g.solid1)
|
||||
return g
|
||||
@@ -0,0 +1,34 @@
|
||||
# Решётка D2Q9 и KBC-проектор. Чистый CuPy; всё считается один раз при импорте.
|
||||
#
|
||||
# Скорости e_i и веса w_i (стандартная нумерация D2Q9):
|
||||
# 0:( 0, 0) 1:(+1,0) 2:(0,+1) 3:(-1,0) 4:(0,-1)
|
||||
# 5:(+1,+1) 6:(-1,+1) 7:(-1,-1) 8:(+1,-1)
|
||||
# OPP[i] — индекс противоположного направления.
|
||||
import backend as B
|
||||
|
||||
cp = B.cp; xp = B.xp; DTYPE = B.DTYPE
|
||||
|
||||
Q = 9
|
||||
Cx = [0, 1, 0, -1, 0, 1, -1, -1, 1] # python-инты (для roll/индексации)
|
||||
Cy = [0, 0, 1, 0, -1, 1, 1, -1, -1]
|
||||
OPP = [0, 3, 4, 1, 2, 7, 8, 5, 6]
|
||||
_Wlist = [4/9, 1/9, 1/9, 1/9, 1/9, 1/36, 1/36, 1/36, 1/36]
|
||||
CS2 = 1.0 / 3.0 # квадрат скорости звука решётки
|
||||
|
||||
# device-векторы
|
||||
Wa = cp.asarray(_Wlist, DTYPE)
|
||||
CXa = cp.asarray(Cx, DTYPE)
|
||||
CYa = cp.asarray(Cy, DTYPE)
|
||||
|
||||
def _build_projector():
|
||||
"""Проектор Ps на сдвиг-моменты {cx²−cy², cx·cy} (ядро KBC-N1).
|
||||
Ps = M⁻¹ · diag(оставить только эти 2 момента) · M, где M — моментный базис."""
|
||||
cxv = cp.asarray(Cx, cp.float64); cyv = cp.asarray(Cy, cp.float64)
|
||||
M = cp.zeros((Q, Q), cp.float64)
|
||||
M[0] = 1.0; M[1] = cxv; M[2] = cyv; M[3] = 3.0*(cxv**2 + cyv**2) - 2.0
|
||||
M[4] = cxv**2 - cyv**2; M[5] = cxv*cyv
|
||||
M[6] = cxv**2*cyv; M[7] = cxv*cyv**2; M[8] = cxv**2*cyv**2
|
||||
Dm = cp.zeros((Q, Q), cp.float64); Dm[4, 4] = Dm[5, 5] = 1.0
|
||||
return (cp.linalg.inv(M) @ Dm @ M)
|
||||
|
||||
Ps = _build_projector().astype(DTYPE) # (Q,Q) проектор на сдвиг
|
||||
@@ -0,0 +1,53 @@
|
||||
# Точка входа модели 2×+SDF. Запуск из этой папки:
|
||||
# python run.py
|
||||
# AMR_STEPS=2000 python run.py # короткий прогон
|
||||
# AMR_CUDAGRAPH=0 python run.py # без CUDA-графов
|
||||
# AMR_FP64=1 python run.py # двойная точность
|
||||
import os, time
|
||||
import config as C
|
||||
import backend as B
|
||||
import solver as S
|
||||
import diagnostics as Diag
|
||||
import visualization as Viz
|
||||
|
||||
HERE = os.path.dirname(os.path.abspath(__file__))
|
||||
OUT = os.path.join(HERE, "out"); os.makedirs(OUT, exist_ok=True)
|
||||
|
||||
def main():
|
||||
cfg = C.from_env()
|
||||
print(f"[2×+SDF] backend={B.BACKEND}; домен {cfg.Nx}x{cfg.Ny}; D={cfg.D}; steps={cfg.steps}; "
|
||||
f"collision=kbc; force=gmem; "
|
||||
f"graphs={'on' if cfg.use_cuda_graph else 'off'}; "
|
||||
f"L1 {cfg.refine}× = {(cfg.bx1-cfg.ax1)*cfg.refine+1}x{(cfg.by1-cfg.ay1)*cfg.refine+1}; "
|
||||
f"вход {cfg.cx/cfg.D:.1f}D / выход {(cfg.Nx-1-cfg.cx)/cfg.D:.1f}D; outlet_uy={cfg.outlet_uy}")
|
||||
|
||||
t0 = time.perf_counter()
|
||||
sim = S.Solver(cfg)
|
||||
res = sim.run(Viz.make_progress("2x+sdf"))
|
||||
print(" время:", round(time.perf_counter() - t0, 1), "s")
|
||||
|
||||
a = Diag.analyze(res, cfg)
|
||||
Diag.print_table(a, cfg)
|
||||
Diag.print_crosscheck(a)
|
||||
Diag.convergence_report(res, cfg)
|
||||
|
||||
Viz.save_gif(res, sim.geo, cfg, os.path.join(OUT, "cyl_2x_sdf.gif"))
|
||||
Viz.save_series_panel(res, a, cfg, os.path.join(OUT, "series_2x_sdf.png"))
|
||||
import numpy as np
|
||||
np.savez(os.path.join(OUT, "probe_2x_sdf.npz"),
|
||||
# сводные числа (самоописываемость для ноутбука: он только читает артефакты)
|
||||
St=a["St"], Cd=a["Cd"], Clrms=a["Clrms"], uyrms=a["uyrms"],
|
||||
Cm=a.get("Cm", float("nan")),
|
||||
Cd_st=a.get("Cd_st", float("nan")), Clrms_st=a.get("Clrms_st", float("nan")),
|
||||
Cm_st=a.get("Cm_st", float("nan")),
|
||||
Ny=cfg.Ny, beta=cfg.D/cfg.Ny, D=cfg.D, U=cfg.U, Re=cfg.Re, steps=cfg.steps,
|
||||
Nx=cfg.Nx, cx=cfg.cx, outlet_uy=cfg.outlet_uy,
|
||||
stress_every=res["stress_every"],
|
||||
# ряды
|
||||
Cl=a["Cl_s"], uy=a["uy"], Fx=B.to_cpu(res["Fx"]), Fy=B.to_cpu(res["Fy"]),
|
||||
Tz=B.to_cpu(res["Tz"]), rho=B.to_cpu(res["rho"]),
|
||||
Fx_st=B.to_cpu(res["Fx_st"]), Fy_st=B.to_cpu(res["Fy_st"]), Tz_st=B.to_cpu(res["Tz_st"]))
|
||||
print("DONE (2×+SDF)")
|
||||
|
||||
if __name__ == "__main__":
|
||||
main()
|
||||
@@ -0,0 +1,95 @@
|
||||
# Серия по блокировке: D фиксирован, Ny растёт → β=D/Ny падает; экстраполяция β→0 даёт
|
||||
# прямое сравнение с БЕЗГРАНИЧНОЙ литературой (Re=150: Cd≈1.33, St≈0.183, rms Cl≈0.3)
|
||||
# без ручной поправки (1−β)². Запуск из этой папки:
|
||||
# python run_blockage.py
|
||||
# AMR_NY_LIST=90,128,180 python run_blockage.py # свой список Ny
|
||||
# AMR_STEPS=20000 AMR_NY_LIST=128 python run_blockage.py # смоук
|
||||
# Стоимость прогона ≈ как у run.py: патч L1 (основная работа) от Ny не зависит (высота 4D фикс.).
|
||||
import os, time
|
||||
import numpy as np
|
||||
import config as C
|
||||
import backend as B
|
||||
import solver as S
|
||||
import diagnostics as Diag
|
||||
import visualization as Viz
|
||||
|
||||
HERE = os.path.dirname(os.path.abspath(__file__))
|
||||
OUT = os.path.join(HERE, "out"); os.makedirs(OUT, exist_ok=True)
|
||||
|
||||
LIT = dict(Cd=1.33, Clrms=0.3, St=0.183) # безграничный цилиндр, Re≈150
|
||||
|
||||
def run_one(Ny):
|
||||
cfg = C.from_env(Ny=Ny, frame_every=0) # серия без кадров/гифок (только ряды)
|
||||
print(f"\n[блокировка] Ny={cfg.Ny} (β={cfg.D/cfg.Ny:.3f}); домен {cfg.Nx}x{cfg.Ny}; "
|
||||
f"патч y=[{cfg.ay1},{cfg.by1}]; steps={cfg.steps}")
|
||||
t0 = time.perf_counter()
|
||||
sim = S.Solver(cfg)
|
||||
res = sim.run(Viz.make_progress(f"Ny={cfg.Ny}"))
|
||||
print(" время:", round(time.perf_counter() - t0, 1), "s")
|
||||
a = Diag.analyze(res, cfg)
|
||||
Diag.print_table(a, cfg)
|
||||
Diag.print_crosscheck(a)
|
||||
Diag.convergence_report(res, cfg)
|
||||
np.savez(os.path.join(OUT, f"probe_blockage_Ny{cfg.Ny}.npz"),
|
||||
St=a["St"], Cd=a["Cd"], Clrms=a["Clrms"], uyrms=a["uyrms"],
|
||||
Cm=a.get("Cm", float("nan")),
|
||||
Cd_st=a.get("Cd_st", float("nan")), Clrms_st=a.get("Clrms_st", float("nan")),
|
||||
Ny=cfg.Ny, beta=cfg.D/cfg.Ny, D=cfg.D, U=cfg.U, Re=cfg.Re, steps=cfg.steps,
|
||||
Cl=a["Cl_s"], uy=a["uy"], rho=B.to_cpu(res["rho"]),
|
||||
Fx=B.to_cpu(res["Fx"]), Fy=B.to_cpu(res["Fy"]), Tz=B.to_cpu(res["Tz"]))
|
||||
row = dict(Ny=cfg.Ny, beta=cfg.D/cfg.Ny, St=a["St"], Cd=a["Cd"], Clrms=a["Clrms"],
|
||||
Cm=a.get("Cm", float("nan")), Cd_st=a.get("Cd_st"))
|
||||
# освободить VRAM перед следующим прогоном (буферы, BC, захваченный CUDA-граф)
|
||||
del sim, res, a
|
||||
B.cp.get_default_memory_pool().free_all_blocks()
|
||||
return row
|
||||
|
||||
def save_extrapolation_figure(rows, fits, path):
|
||||
import matplotlib; matplotlib.use("Agg")
|
||||
import matplotlib.pyplot as plt
|
||||
betas = np.array([r["beta"] for r in rows]); bb = np.linspace(0, betas.max()*1.05, 100)
|
||||
panels = [("Cd", "⟨Cd⟩", LIT["Cd"], "лит. 1.33"),
|
||||
("Clrms", "rms Cl", LIT["Clrms"], "лит. ~0.3"),
|
||||
("St", "St", LIT["St"], "лит. 0.183")]
|
||||
fig, axs = plt.subplots(1, 3, figsize=(13.5, 4.2), dpi=110)
|
||||
for ax, (key, lab, lit, litlab) in zip(axs, panels):
|
||||
vals = np.array([r[key] for r in rows]); ft = fits[key]
|
||||
ax.plot(betas, vals, "o", ms=7, color="#d32f2f", label="серия (raw)")
|
||||
al, sl = ft["lin"]; aq, sq = ft["quad"]
|
||||
ax.plot(bb, al + sl*bb, "-", lw=1.2, color="#1976d2", label=f"a+b·β → {al:.3f}")
|
||||
ax.plot(bb, aq + sq*bb**2, "--", lw=1.2, color="#388e3c", label=f"a+b·β² → {aq:.3f}")
|
||||
ax.axhline(lit, color="gray", lw=1.0, ls=":", label=litlab)
|
||||
if key == "St": # для St полезен и кинематически поправленный вид
|
||||
ax.plot(betas, vals*(1-betas), "s", ms=6, mfc="none", color="#7b1fa2", label="St·(1−β)")
|
||||
ax.set_xlabel("β = D/Ny"); ax.set_ylabel(lab); ax.set_xlim(0, bb[-1])
|
||||
ax.legend(fontsize=8); ax.set_title(f"{lab}(β): экстраполяция β→0", fontsize=10)
|
||||
fig.suptitle("Серия по блокировке: стеснение каналом — свойство ПОСТАНОВКИ, не считывания силы",
|
||||
fontsize=11)
|
||||
fig.tight_layout(rect=(0, 0, 1, 0.93))
|
||||
fig.savefig(path, bbox_inches="tight"); plt.close(fig)
|
||||
print("saved", os.path.basename(path))
|
||||
|
||||
def main():
|
||||
ny_list = [int(s) for s in os.environ.get("AMR_NY_LIST", "90,128,180").split(",")]
|
||||
print(f"[серия по блокировке] backend={B.BACKEND}; Ny={ny_list}")
|
||||
rows = [run_one(Ny) for Ny in ny_list]
|
||||
if len(rows) < 2:
|
||||
print("\n(одна точка — фит/экстраполяция пропущены)")
|
||||
return
|
||||
betas = [r["beta"] for r in rows]
|
||||
fits = {key: Diag.fit_blockage(betas, [r[key] for r in rows]) for key in ("Cd", "Clrms", "St")}
|
||||
Diag.print_blockage_table(rows, fits)
|
||||
np.savez(os.path.join(OUT, "blockage_summary.npz"),
|
||||
Ny=np.array([r["Ny"] for r in rows]), beta=np.array(betas),
|
||||
St=np.array([r["St"] for r in rows]), Cd=np.array([r["Cd"] for r in rows]),
|
||||
Clrms=np.array([r["Clrms"] for r in rows]), Cm=np.array([r["Cm"] for r in rows]),
|
||||
Cd_st=np.array([r["Cd_st"] if r["Cd_st"] is not None else np.nan for r in rows]),
|
||||
Cd0_lin=fits["Cd"]["lin"][0], Cd0_quad=fits["Cd"]["quad"][0],
|
||||
Cl0_lin=fits["Clrms"]["lin"][0], Cl0_quad=fits["Clrms"]["quad"][0],
|
||||
St0_lin=fits["St"]["lin"][0], St0_quad=fits["St"]["quad"][0],
|
||||
lit_Cd=LIT["Cd"], lit_Clrms=LIT["Clrms"], lit_St=LIT["St"])
|
||||
save_extrapolation_figure(rows, fits, os.path.join(OUT, "blockage_extrapolation.png"))
|
||||
print("DONE (серия по блокировке)")
|
||||
|
||||
if __name__ == "__main__":
|
||||
main()
|
||||
@@ -0,0 +1,184 @@
|
||||
# ФАКТОРНОЕ ИССЛЕДОВАНИЕ постановки: какие границы расчётной области отвечают за завышение
|
||||
# Cd / rms Cl относительно безграничной литературы (Re=150: Cd≈1.33, rms Cl≈0.3, St≈0.183).
|
||||
#
|
||||
# Контекст (эксперимент №0 — серия по боковой блокировке, см. ноутбук, раздел 23): Cd(β)
|
||||
# немонотонен (1.80 → 1.77 → 1.87 при β = 0.178 → 0.125 → 0.089) → боковое стеснение НЕ
|
||||
# доминирующий фактор. Подозреваемые — продольные границы, которые серия держала константой:
|
||||
# вход 2.5D (жёсткий равномерный профиль), выход 8.3D (Zou-He давление с жёстким uy=0,
|
||||
# отражает вихри дорожки), и их связка.
|
||||
#
|
||||
# Здесь каждый фактор варьируется ОТДЕЛЬНО при фиксированной боковой блокировке (Ny=90):
|
||||
# изоляция вклада входа, выхода и выходного ГУ. Все результаты накапливаются в out/factors.csv
|
||||
# (одна строка на прогон, с полной конфигурацией) — реестр для документирования исследования.
|
||||
#
|
||||
# Запуск: python run_factors.py
|
||||
# AMR_FACTORS=B0_baseline,F4_outlet_uy python run_factors.py # подмножество факторов
|
||||
# AMR_STEPS=20000 python run_factors.py # смоук
|
||||
import os, csv, time
|
||||
import numpy as np
|
||||
import config as C
|
||||
import backend as B
|
||||
import solver as S
|
||||
import diagnostics as Diag
|
||||
import visualization as Viz
|
||||
|
||||
HERE = os.path.dirname(os.path.abspath(__file__))
|
||||
OUT = os.path.join(HERE, "out"); os.makedirs(OUT, exist_ok=True)
|
||||
LIT = dict(Cd=1.33, Clrms=0.3, St=0.183)
|
||||
|
||||
# (имя, overrides конфига, что изолирует)
|
||||
# Эксперимент №1 (B0–F5, прогнан 2026-06-10): итоги — мягкий uy-выход (F4) выигрывает по всем
|
||||
# метрикам; ВСЕ длинные домены без демпфера уходят в смещённый режим ⟨ρ⟩≈1.66 (акустический
|
||||
# резонатор «вход-скорость/выход-давление» + накачка массы входом) или взрываются (F5).
|
||||
# Эксперимент №2 (F6–F8): губка перед выходом — лечит акустику, открывает длинные домены.
|
||||
FACTORS = [
|
||||
("B0_baseline", dict(),
|
||||
"база 174x90: вход 2.5D, выход 8.3D, uy=0 (контроль; кросс-чек уже с конвективным членом)"),
|
||||
("F1_outlet_far", dict(Nx=300),
|
||||
"выход 8.3D -> 16.2D (вход 2.5D как в базе) — изолирует близость выхода"),
|
||||
("F2_inlet_far", dict(Nx=230, cx=96),
|
||||
"вход 2.5D -> 6.0D (выход 8.3D как в базе) — изолирует близость входа"),
|
||||
("F3_inout_far", dict(Nx=390, cx=96),
|
||||
"вход 6.0D + выход 18.3D — оба продольных фактора вместе"),
|
||||
("F4_outlet_uy", dict(outlet_uy="extrapolate"),
|
||||
"мягкий выход: uy нуль-градиент вместо uy=0 (геометрия базовая) — изолирует отражение вихрей"),
|
||||
("F5_combo", dict(Nx=390, cx=96, outlet_uy="extrapolate"),
|
||||
"вход 6.0D + выход 18.3D + мягкий uy-выход — всё вместе"),
|
||||
# --- эксперимент №2: губка (absorbing layer) ---
|
||||
("F6_sponge", dict(outlet_uy="extrapolate", sponge_len=32),
|
||||
"база 174x90 + мягкий uy + ГУБКА 32 столбца (nu x30) — валидация губки против F4"),
|
||||
("F7_long_sponge", dict(Nx=390, cx=96, outlet_uy="extrapolate", sponge_len=32),
|
||||
"вход 6D + выход 18.3D + мягкий uy + губка — целевая «чистая» постановка"),
|
||||
("F8_long_sp_uy0", dict(Nx=390, cx=96, sponge_len=32),
|
||||
"длинный домен + губка, но ЖЁСТКОЕ uy=0 — изолирует вклад губки от uy-режима"),
|
||||
]
|
||||
|
||||
CSV_FIELDS = ["when", "name", "status", "steps", "Nx", "Ny", "cx", "in_D", "out_D", "outlet_uy",
|
||||
"sponge_len", "beta", "St", "St_corr", "Cd", "Clrms", "Cm", "Cd_st", "Clrms_st",
|
||||
"dCd_pct", "dCl_pct", "rho_last", "cd_win_std", "cl_win_std", "wall_s", "desc"]
|
||||
|
||||
def window_stats(res, cfg, nwin=10):
|
||||
"""Срединная изменчивость: std оконных средних Cd и оконных rms Cl (окна 2..nwin,
|
||||
первое окно с разгоном пропускается). Метрика низкочастотной модуляции дорожки."""
|
||||
Fx = B.to_cpu(res["Fx"]); Fy = B.to_cpu(res["Fy"])
|
||||
Cd = 2.0*Fx/(cfg.U**2*res["DL"]); Cl = 2.0*Fy/(cfg.U**2*res["DL"])
|
||||
n = len(Cd); w = n // nwin
|
||||
cds = [float(Cd[k*w:(k+1)*w].mean()) for k in range(1, nwin)]
|
||||
cls = [float(Cl[k*w:(k+1)*w].std()) for k in range(1, nwin)]
|
||||
return float(np.std(cds)), float(np.std(cls))
|
||||
|
||||
def run_factor(name, overrides, desc):
|
||||
cfg = C.from_env(frame_every=0, **overrides)
|
||||
in_D = cfg.cx / cfg.D; out_D = (cfg.Nx - 1 - cfg.cx) / cfg.D
|
||||
print(f"\n=== Фактор {name}: {desc}")
|
||||
print(f" домен {cfg.Nx}x{cfg.Ny}; вход {in_D:.1f}D / выход {out_D:.1f}D; "
|
||||
f"outlet_uy={cfg.outlet_uy}; губка={cfg.sponge_len or '—'}; "
|
||||
f"патч x=[{cfg.ax1},{cfg.bx1}] y=[{cfg.ay1},{cfg.by1}]; steps={cfg.steps}")
|
||||
t0 = time.perf_counter()
|
||||
sim = S.Solver(cfg)
|
||||
res = sim.run(Viz.make_progress(name[:8]))
|
||||
wall = time.perf_counter() - t0
|
||||
print(" время:", round(wall, 1), "s")
|
||||
status = "ok" if len(res["Fx"]) == cfg.steps + 1 else f"BLEW_UP@{len(res['Fx'])}"
|
||||
a = Diag.analyze(res, cfg)
|
||||
Diag.print_table(a, cfg)
|
||||
Diag.print_crosscheck(a)
|
||||
Diag.convergence_report(res, cfg)
|
||||
cd_ws, cl_ws = window_stats(res, cfg)
|
||||
rho_last = float(B.to_cpu(res["rho"])[-1])
|
||||
np.savez(os.path.join(OUT, f"factor_{name}.npz"),
|
||||
St=a["St"], Cd=a["Cd"], Clrms=a["Clrms"], uyrms=a["uyrms"],
|
||||
Cm=a.get("Cm", float("nan")),
|
||||
Cd_st=a.get("Cd_st", float("nan")), Clrms_st=a.get("Clrms_st", float("nan")),
|
||||
Nx=cfg.Nx, Ny=cfg.Ny, cx=cfg.cx, D=cfg.D, U=cfg.U, Re=cfg.Re, steps=cfg.steps,
|
||||
outlet_uy=cfg.outlet_uy, beta=cfg.D/cfg.Ny,
|
||||
Cl=a["Cl_s"], uy=a["uy"], rho=B.to_cpu(res["rho"]),
|
||||
Fx=B.to_cpu(res["Fx"]), Fy=B.to_cpu(res["Fy"]), Tz=B.to_cpu(res["Tz"]))
|
||||
beta = cfg.D / cfg.Ny
|
||||
dcd = abs(a.get("Cd_st", float("nan")) - a["Cd"]) / a["Cd"] * 100.0
|
||||
dcl = abs(a.get("Clrms_st", float("nan")) - a["Clrms"]) / max(a["Clrms"], 1e-12) * 100.0
|
||||
row = dict(when=time.strftime("%Y-%m-%d %H:%M:%S"), name=name, status=status, steps=cfg.steps,
|
||||
Nx=cfg.Nx, Ny=cfg.Ny, cx=cfg.cx, in_D=round(in_D, 2), out_D=round(out_D, 2),
|
||||
outlet_uy=cfg.outlet_uy, sponge_len=cfg.sponge_len, beta=round(beta, 4),
|
||||
St=round(a["St"], 4), St_corr=round(a["St"]*(1-beta), 4),
|
||||
Cd=round(a["Cd"], 4), Clrms=round(a["Clrms"], 4),
|
||||
Cm=round(a.get("Cm", float("nan")), 6),
|
||||
Cd_st=round(a.get("Cd_st", float("nan")), 4),
|
||||
Clrms_st=round(a.get("Clrms_st", float("nan")), 4),
|
||||
dCd_pct=round(dcd, 2), dCl_pct=round(dcl, 2),
|
||||
rho_last=round(rho_last, 5), cd_win_std=round(cd_ws, 4),
|
||||
cl_win_std=round(cl_ws, 4), wall_s=round(wall, 1), desc=desc)
|
||||
del sim, res, a
|
||||
B.cp.get_default_memory_pool().free_all_blocks()
|
||||
return row
|
||||
|
||||
def append_csv(rows):
|
||||
path = os.path.join(OUT, "factors.csv")
|
||||
if os.path.exists(path):
|
||||
with open(path, encoding="utf-8") as fh:
|
||||
header = fh.readline().strip()
|
||||
if header != ",".join(CSV_FIELDS):
|
||||
# схема реестра изменилась — старый файл сохраняем (история не теряется)
|
||||
k = 1
|
||||
while os.path.exists(os.path.join(OUT, f"factors_v{k}.csv")):
|
||||
k += 1
|
||||
os.rename(path, os.path.join(OUT, f"factors_v{k}.csv"))
|
||||
print(f"[реестр] схема CSV изменилась — прежний файл сохранён как factors_v{k}.csv")
|
||||
new = not os.path.exists(path)
|
||||
with open(path, "a", newline="", encoding="utf-8") as fh:
|
||||
wr = csv.DictWriter(fh, fieldnames=CSV_FIELDS)
|
||||
if new:
|
||||
wr.writeheader()
|
||||
for r in rows:
|
||||
wr.writerow(r)
|
||||
print("CSV:", path, f"(+{len(rows)} строк)")
|
||||
|
||||
def save_summary_figure(rows, path):
|
||||
import matplotlib; matplotlib.use("Agg")
|
||||
import matplotlib.pyplot as plt
|
||||
names = [r["name"] for r in rows]
|
||||
xs = np.arange(len(rows))
|
||||
fig, axs = plt.subplots(1, 3, figsize=(15, 4.6))
|
||||
ax = axs[0]
|
||||
ax.bar(xs, [r["Cd"] for r in rows], color="#1f6feb", label="⟨Cd⟩ GMEM")
|
||||
ax.plot(xs, [r["Cd_st"] for r in rows], "o", color="#ff8f00", label="⟨Cd⟩ ∮[σ·n−ρu(u·n)]ds")
|
||||
ax.axhline(LIT["Cd"], color="#d1495b", ls="--", label=f"лит. безгранич. {LIT['Cd']}")
|
||||
ax.set_ylabel("⟨Cd⟩"); ax.set_title("Сопротивление по факторам")
|
||||
ax = axs[1]
|
||||
ax.bar(xs, [r["Clrms"] for r in rows], color="#2e7d32")
|
||||
ax.axhline(LIT["Clrms"], color="#d1495b", ls="--", label=f"лит. ~{LIT['Clrms']}")
|
||||
ax.set_ylabel("rms Cl"); ax.set_title("Амплитуда подъёмной силы")
|
||||
ax = axs[2]
|
||||
ax.bar(xs, [r["cd_win_std"] for r in rows], color="#7b1fa2")
|
||||
ax.set_ylabel("std оконных ⟨Cd⟩"); ax.set_title("НЧ-модуляция (меньше — стабильнее)")
|
||||
for ax in axs:
|
||||
ax.set_xticks(xs)
|
||||
ax.set_xticklabels([f"{r['name']}\nвх {r['in_D']}D / вых {r['out_D']}D\n"
|
||||
f"uy:{r['outlet_uy']} губ:{r['sponge_len']}" for r in rows], fontsize=7)
|
||||
ax.legend(fontsize=8)
|
||||
fig.suptitle("Факторное исследование постановки (Ny=90, β=0.178 фикс., 100k шагов)", fontsize=12)
|
||||
fig.tight_layout(rect=(0, 0, 1, 0.94))
|
||||
fig.savefig(path, bbox_inches="tight"); plt.close(fig)
|
||||
print("saved", os.path.basename(path))
|
||||
|
||||
def main():
|
||||
sel = os.environ.get("AMR_FACTORS")
|
||||
todo = FACTORS if not sel else [f for f in FACTORS if f[0] in sel.split(",")]
|
||||
print(f"[факторы] backend={B.BACKEND}; прогонов: {len(todo)}")
|
||||
rows = [run_factor(*f) for f in todo]
|
||||
print("\n=== Сводка факторного исследования ===")
|
||||
hdr = (f"{'фактор':15}{'вх,D':>6}{'вых,D':>7}{'uy':>12}{'губка':>6}{'St':>7}{'<Cd>':>8}"
|
||||
f"{'rms Cl':>8}{'Cd σ·n':>8}{'dCd%':>6}{'<ρ>фин':>8}{'модул.':>8} статус")
|
||||
print(hdr)
|
||||
for r in rows:
|
||||
print(f"{r['name']:15}{r['in_D']:>6}{r['out_D']:>7}{r['outlet_uy']:>12}{r['sponge_len']:>6}"
|
||||
f"{r['St']:>7.3f}{r['Cd']:>8.3f}{r['Clrms']:>8.4f}{r['Cd_st']:>8.3f}{r['dCd_pct']:>6.1f}"
|
||||
f"{r['rho_last']:>8.4f}{r['cd_win_std']:>8.4f} {r['status']}")
|
||||
print(f"{'лит. безгранич.':15}{'':>6}{'':>7}{'':>12}{'':>6}{LIT['St']:>7.3f}{LIT['Cd']:>8.3f}{LIT['Clrms']:>8.4f}")
|
||||
append_csv(rows)
|
||||
if len(rows) > 1:
|
||||
save_summary_figure(rows, os.path.join(OUT, "factors_summary.png"))
|
||||
print("DONE (факторы)")
|
||||
|
||||
if __name__ == "__main__":
|
||||
main()
|
||||
@@ -0,0 +1,191 @@
|
||||
# Оркестратор модели 2×+SDF: держит persistent-буферы, собирает один шаг из компонентов и гоняет цикл.
|
||||
# Один шаг = замкнутая step() над фикс. буферами без host→device (готова к захвату в CUDA-граф).
|
||||
import math
|
||||
import backend as B
|
||||
import lattice as L
|
||||
import geometry as G
|
||||
import equilibrium as E
|
||||
import collision as Coll
|
||||
import streaming as Strm
|
||||
import boundary as Bnd
|
||||
import forces as Fz
|
||||
import forces_stress as FzS
|
||||
import amr as Amr
|
||||
|
||||
cp = B.cp; xp = B.xp; DTYPE = B.DTYPE; to_cpu = B.to_cpu
|
||||
CYa = L.CYa
|
||||
|
||||
def _capture(step_fn):
|
||||
"""Захват одного шага в CUDA-граф. При сбое гарантированно выводит поток из режима записи."""
|
||||
cp.cuda.Device().synchronize()
|
||||
s = cp.cuda.Stream(non_blocking=True)
|
||||
with s:
|
||||
s.begin_capture()
|
||||
try:
|
||||
step_fn()
|
||||
except Exception:
|
||||
try: s.end_capture()
|
||||
except Exception: pass
|
||||
raise
|
||||
return s.end_capture()
|
||||
|
||||
class Solver:
|
||||
def __init__(self, cfg):
|
||||
self.cfg = cfg
|
||||
self.collide = Coll.collide # KBC-N1 (единственный оператор столкновения)
|
||||
self.force = Fz.force_gmem # галилей-инвариантный обмен импульсом
|
||||
self.geo = G.build(cfg)
|
||||
self.P1 = self.geo.P1
|
||||
self.beta1 = cfg.beta1
|
||||
if cfg.sponge_len:
|
||||
# губка: τ(x) = τ0 + (mult−1)·ν·s(x)/cs², s — smoothstep 0→1 в последних sponge_len
|
||||
# столбцах. β становится ПОЛЕМ (Nx,) → broadcast в kbc_collide; патч L1 губки не видит
|
||||
# (гарантировано assert'ом в config), там остаётся скалярный beta1.
|
||||
xs = cp.arange(cfg.Nx, dtype=DTYPE)
|
||||
s = cp.clip((xs - (cfg.Nx - 1 - cfg.sponge_len)) / float(cfg.sponge_len), 0.0, 1.0)
|
||||
s = s*s*(3.0 - 2.0*s)
|
||||
tau_x = cfg.tau0 + (cfg.sponge_nu_mult - 1.0) * cfg.nu / (1.0/3.0) * s
|
||||
self.beta0 = (1.0 / (2.0 * tau_x)).astype(DTYPE)[None, :] # (1,Nx) → broadcast по y
|
||||
else:
|
||||
self.beta0 = cfg.beta0
|
||||
self.R01 = cfg.R01; self.Rf01 = cfg.Rf01
|
||||
self.refine = cfg.refine
|
||||
# ЗОНД скорости в следе: физическая точка (px,py) в координатах L0.
|
||||
# Берём её с ТОНКОЙ сетки L1 (меньше численного размытия вихрей → чище St / rms u_y),
|
||||
# если точка внутри патча; иначе fallback на L0. Индексы фиксированы → graph-safe.
|
||||
self.probe_l1 = (cfg.ax1 <= cfg.px <= cfg.bx1) and (cfg.ay1 <= cfg.py <= cfg.by1)
|
||||
if self.probe_l1:
|
||||
self.pfx = (cfg.px - cfg.ax1) * cfg.refine
|
||||
self.pfy = (cfg.py - cfg.ay1) * cfg.refine
|
||||
else:
|
||||
self.pfx = cfg.px; self.pfy = cfg.py
|
||||
# ГУ (SDF/Bouzidi) на L0 (цилиндр+стенки) и L1 (цилиндр)
|
||||
self.BC0 = Bnd.build_bc(self.geo.solid0, self.geo.phi0, self.geo.cyl0)
|
||||
self.BC1 = Bnd.build_bc(self.geo.solid1, self.geo.phi1, self.geo.solid1)
|
||||
self.fl1 = self.geo.fl1
|
||||
# persistent-состояние
|
||||
Ny, Nx = cfg.Ny, cfg.Nx; Nfy, Nfx = self.P1["Nfy"], self.P1["Nfx"]
|
||||
self.F0 = E.feq(xp.ones((Ny, Nx), DTYPE), xp.zeros((2, Ny, Nx), DTYPE)).copy()
|
||||
self.F1 = E.feq(xp.ones((Nfy, Nfx), DTYPE), xp.zeros((2, Nfy, Nfx), DTYPE)).copy()
|
||||
self.onecol = xp.ones((Ny, 1), DTYPE)
|
||||
self.uin = xp.zeros((2, Ny, 1), DTYPE) # вход; обновляется СНАРУЖИ захвата
|
||||
self.Fx_s = xp.zeros((), DTYPE); self.Fy_s = xp.zeros((), DTYPE); self.uy_s = xp.zeros((), DTYPE)
|
||||
self.Tz_s = xp.zeros((), DTYPE)
|
||||
# координаты узлов L1 относительно центра тела — плечи торка (прекомпьют вне захвата)
|
||||
ccx1 = (cfg.cx - cfg.ax1) * cfg.refine; ccy1 = (cfg.cy - cfg.ay1) * cfg.refine
|
||||
Yf, Xf = cp.meshgrid(cp.arange(Nfy), cp.arange(Nfx), indexing="ij")
|
||||
self.rxy1 = ((Xf - ccx1).astype(DTYPE), (Yf - ccy1).astype(DTYPE))
|
||||
# контур кросс-чека силы σ·n (forces_stress; считается редко, вне CUDA-графа)
|
||||
self.sprobe = FzS.build_probe(cfg) if cfg.stress_every else None
|
||||
# разгон входа на device: U*smoothstep(t/ramp)
|
||||
ts = cp.minimum(cp.arange(cfg.steps + 1) / cfg.ramp, 1.0)
|
||||
self.ramp_u = (cfg.U * (ts*ts*(3 - 2*ts))).astype(DTYPE)
|
||||
# стартовое возмущение: поперечный sin-импульс u_y входа в окне [ramp, ramp+pert_dur] (триггер схода)
|
||||
uyp = cp.zeros(cfg.steps + 1, DTYPE)
|
||||
if cfg.pert_dur > 0 and cfg.pert_amp != 0.0:
|
||||
k = cp.arange(cfg.steps + 1)
|
||||
inw = (k >= cfg.ramp) & (k < cfg.ramp + cfg.pert_dur)
|
||||
ph = (k - cfg.ramp).astype(DTYPE) / float(cfg.pert_dur)
|
||||
uyp = cp.where(inw, (cfg.pert_amp * cfg.U) * cp.sin(math.pi * ph), 0.0).astype(DTYPE)
|
||||
self.ramp_uy = uyp
|
||||
# ⟨ρ⟩ — диагностика дрейфа массы. Считать ТОЛЬКО по жидкости: узлы внутри цилиндра —
|
||||
# фиктивные (их популяции не имеют физического смысла, ГУ их не трогает), а их плотность
|
||||
# уезжает независимо. При D=12 на 120×60 это 1.6% узлов и смещение ⟨ρ⟩ на −5.4e-4 —
|
||||
# ровно тот разряд, в котором ⟨ρ⟩ и цитируется в отчётах (1.0000 / 1.001).
|
||||
self.fluid0 = (~self.geo.solid0).astype(DTYPE)
|
||||
self._ncells = float(self.fluid0.sum())
|
||||
|
||||
def step(self):
|
||||
F0 = self.F0; F1 = self.F1; P1 = self.P1
|
||||
feq = E.feq; macros = E.macros; stream = Strm.stream
|
||||
# ----- уровень 0 -----
|
||||
rho, u0 = macros(F0)
|
||||
pre = F0.copy()
|
||||
post0 = self.collide(F0, feq(rho, u0), self.beta0)
|
||||
f0 = stream(post0); f0 = Bnd.apply_bc(f0, post0, self.BC0)
|
||||
# порядок важен: стенки снимают y-заворот roll, затем Zou-He закрывает x-заворот НА ВЕСЬ
|
||||
# столбец, включая угловые узлы (иначе вход и выход связаны напрямую — см. boundary.py)
|
||||
f0 = Bnd.free_slip_walls(f0) # верх/низ — скользящие стенки (specular)
|
||||
f0 = Bnd.channel_bc(f0, self.uin, outlet_uy=self.cfg.outlet_uy) # Zou-He вход/выход
|
||||
# ----- уровень 1: refine подшагов по времени с временно́й интерполяцией ghost -----
|
||||
g1o = Amr.ghost(pre, P1, self.R01); g1n = Amr.ghost(f0, P1, self.R01)
|
||||
f1 = F1
|
||||
# сила/момент: меряем на КАЖДОМ подшаге L1 и усредняем (анти-алиасинг; раньше — мгновенная
|
||||
# на последнем подшаге). Средний Cd не меняется, ряды становятся глаже.
|
||||
fxa = B.ZERO; fya = B.ZERO; tza = B.ZERO
|
||||
for s1 in range(self.refine):
|
||||
r1, u1 = macros(f1); post1 = self.collide(f1, feq(r1, u1), self.beta1)
|
||||
f1 = stream(post1); f1 = Bnd.apply_bc(f1, post1, self.BC1)
|
||||
fx, fy, tz = self.force(post1, f1, self.BC1, self.rxy1)
|
||||
fxa = fxa + fx; fya = fya + fy; tza = tza + tz
|
||||
w = (s1 + 1) / self.refine # временна́я интерполяция ghost: o→n
|
||||
f1 = Amr.fill(f1, (1 - w)*g1o + w*g1n)
|
||||
self.Fx_s[...] = fxa / self.refine; self.Fy_s[...] = fya / self.refine
|
||||
self.Tz_s[...] = tza / self.refine
|
||||
f0 = Amr.restrict(f1, f0, P1, self.Rf01, self.fl1)
|
||||
F1[...] = f1
|
||||
# зонд скорости в следе (с тонкой сетки L1, если точка внутри патча) + запись состояния
|
||||
psrc = f1 if self.probe_l1 else f0
|
||||
col = psrc[:, self.pfy, self.pfx]; self.uy_s[...] = (CYa*col).sum() / col.sum()
|
||||
F0[...] = f0
|
||||
|
||||
def run(self, progress=None):
|
||||
cfg = self.cfg
|
||||
Fx = xp.zeros(cfg.steps + 1, DTYPE); Fy = xp.zeros(cfg.steps + 1, DTYPE); uy = xp.zeros(cfg.steps + 1, DTYPE)
|
||||
Tz = xp.zeros(cfg.steps + 1, DTYPE) # момент (торк) вокруг центра тела
|
||||
rho = xp.zeros(cfg.steps + 1, DTYPE) # ⟨ρ⟩(t) — диагностика дрейфа массы
|
||||
# кросс-чек σ·n: редкая выборка (каждые stress_every шагов), вне CUDA-графа
|
||||
nst = (cfg.steps // cfg.stress_every + 1) if cfg.stress_every else 0
|
||||
Fx_st = xp.zeros(nst, DTYPE); Fy_st = xp.zeros(nst, DTYPE); Tz_st = xp.zeros(nst, DTYPE)
|
||||
kst = 0
|
||||
frames = []
|
||||
graph = None; graph_failed = False
|
||||
def snap(): return (self.F0.copy(), self.F1.copy())
|
||||
def restore(s): self.F0[...] = s[0]; self.F1[...] = s[1]
|
||||
def maxdiff(s): return max(float(xp.abs(s[0] - self.F0).max()), float(xp.abs(s[1] - self.F1).max()))
|
||||
for t in range(cfg.steps + 1):
|
||||
self.uin[0, :, 0] = self.ramp_u[t]
|
||||
self.uin[1, :, 0] = self.ramp_uy[t] # поперечный импульс-триггер (0 вне окна возмущения)
|
||||
did = False
|
||||
if t == 1 and cfg.use_cuda_graph and not graph_failed:
|
||||
# захват + РАНТАЙМ-САМОПРОВЕРКА (граф vs eager на одном шаге); при расхождении — откат
|
||||
try:
|
||||
graph = _capture(self.step)
|
||||
before = snap(); graph.launch(); gres = snap()
|
||||
restore(before); self.step()
|
||||
d = maxdiff(gres)
|
||||
if d > 1e-3:
|
||||
print(f"\n[graph] self-check разошёлся (Δ={d:.2e}) — отключаю графы, работаю eager")
|
||||
graph = None; graph_failed = True
|
||||
else:
|
||||
print(f"\n[graph] self-check OK (Δ={d:.2e}) — CUDA graphs активны")
|
||||
did = True
|
||||
except Exception as e:
|
||||
print(f"\n[graph] недоступны ({type(e).__name__}: {e}) — работаю eager")
|
||||
import traceback as _tb; _tb.print_exc()
|
||||
graph = None; graph_failed = True
|
||||
try: cp.cuda.Device().synchronize()
|
||||
except Exception: pass
|
||||
if not did:
|
||||
if graph is not None and t >= 1:
|
||||
graph.launch()
|
||||
else:
|
||||
self.step()
|
||||
Fx[t] = self.Fx_s; Fy[t] = self.Fy_s; Tz[t] = self.Tz_s; uy[t] = self.uy_s
|
||||
rho[t] = (self.F0.sum(0) * self.fluid0).sum() / self._ncells # ⟨ρ⟩ по жидкости (дрейф массы)
|
||||
if cfg.stress_every and t % cfg.stress_every == 0:
|
||||
# eager-вычисление на текущем F1 после launch/step — упорядочено на стриме,
|
||||
# сам захваченный граф не трогаем (паттерн как у rho[t] выше)
|
||||
sfx, sfy, stz = FzS.force_stress(self.F1, self.sprobe, cfg.tau1)
|
||||
Fx_st[kst] = sfx; Fy_st[kst] = sfy; Tz_st[kst] = stz; kst += 1
|
||||
if progress is not None:
|
||||
progress(t, cfg.steps)
|
||||
if t % 500 == 0 and not bool(xp.isfinite(self.F0).all()):
|
||||
print("BLEW UP", t)
|
||||
Fx = Fx[:t]; Fy = Fy[:t]; Tz = Tz[:t]; uy = uy[:t]; rho = rho[:t]
|
||||
Fx_st = Fx_st[:kst]; Fy_st = Fy_st[:kst]; Tz_st = Tz_st[:kst]
|
||||
break
|
||||
if cfg.frame_every and t % cfg.frame_every == 0:
|
||||
frames.append((t, to_cpu(E.macros(self.F0)[1]), to_cpu(E.macros(self.F1)[1])))
|
||||
return dict(Fx=Fx, Fy=Fy, Tz=Tz, uy=uy, rho=rho, frames=frames, DL=cfg.DL,
|
||||
Fx_st=Fx_st, Fy_st=Fy_st, Tz_st=Tz_st, stress_every=cfg.stress_every)
|
||||
@@ -0,0 +1,17 @@
|
||||
# ПЕРЕНОС (streaming): f_i(x, t+1) = f_i(x − c_i, t).
|
||||
# Реализован пул-схемой через cp.roll по осям (y, x). 9 сдвигов в новый буфер.
|
||||
#
|
||||
# Заметка для ревизии: roll создаёт по ядру на направление (9 ядер). Это НЕ влияет на физику,
|
||||
# только на число запусков ядер. Главный способ убрать накладные расходы — CUDA-graph (см. solver.py),
|
||||
# а не переписывание переноса. Альтернатива (один gather-кернел по предвыч. индексам) — точка расширения.
|
||||
import backend as B
|
||||
import lattice as L
|
||||
|
||||
xp = B.xp
|
||||
Q = L.Q; Cx = L.Cx; Cy = L.Cy
|
||||
|
||||
def stream(f):
|
||||
fs = xp.empty_like(f)
|
||||
for i in range(Q):
|
||||
fs[i] = xp.roll(f[i], (Cy[i], Cx[i]), (0, 1))
|
||||
return fs
|
||||
@@ -0,0 +1,92 @@
|
||||
# Визуализация (host): гифка |ω| с наложением тонкого патча L1 и прогресс-бар.
|
||||
# matplotlib/PIL работают только с host-массивами — кадры уже перенесены на host в solver.run().
|
||||
import os
|
||||
import time as _time
|
||||
import numpy as np
|
||||
|
||||
def make_progress(label, width=32, every=0.1):
|
||||
"""ASCII-шкала прогресса по выполненным шагам из заданных (обновление на месте через \\r)."""
|
||||
st = {"t0": _time.perf_counter(), "last": 0.0}
|
||||
def cb(t, total):
|
||||
now = _time.perf_counter(); done = (t >= total)
|
||||
if not done and (now - st["last"] < every):
|
||||
return
|
||||
st["last"] = now; frac = (t + 1) / (total + 1)
|
||||
filled = int(round(width * frac)); bar = "#"*filled + "-"*(width - filled)
|
||||
el = now - st["t0"]; eta = (el/frac - el) if frac > 0 else 0.0
|
||||
print(f"\r {label:8}[{bar}] {frac*100:5.1f}% ({t}/{total}) {el:5.1f}s ETA {eta:5.1f}s",
|
||||
end=("\n" if done else ""), flush=True)
|
||||
return cb
|
||||
|
||||
def _vort(u):
|
||||
return (np.roll(u[1], -1, 1) - np.roll(u[1], 1, 1))*0.5 - (np.roll(u[0], -1, 0) - np.roll(u[0], 1, 0))*0.5
|
||||
|
||||
def save_series_panel(res, a, cfg, path):
|
||||
"""Готовая png-панель для ноутбука (он код не исполняет — только вставляет картинки):
|
||||
(1) ряды Cd(t)/Cl(t) + редкие точки кросс-чека ∮σ·n ds; (2) Cm(t) — контроль симметрии;
|
||||
(3) спектр зонда u_y с пиком St; (4) ⟨ρ⟩(t) — контроль дрейфа массы."""
|
||||
import backend as B
|
||||
import matplotlib; matplotlib.use("Agg")
|
||||
import matplotlib.pyplot as plt
|
||||
Cd = a["Cd_s"]; Cl = a["Cl_s"]; n = len(Cd); h0 = a["h0"]; t = np.arange(n)
|
||||
fig, axs = plt.subplots(2, 2, figsize=(12.5, 7.2), dpi=110)
|
||||
ax = axs[0, 0]
|
||||
ax.plot(t, Cd, lw=0.7, color="#d32f2f", label=f"Cd (GMEM), ⟨Cd⟩={a['Cd']:.3f}")
|
||||
ax.plot(t, Cl, lw=0.7, color="#1976d2", label=f"Cl (GMEM), rms={a['Clrms']:.3f}")
|
||||
if "Cd_st_s" in a:
|
||||
K = a["stress_every"]; ts = np.arange(len(a["Cd_st_s"])) * K
|
||||
ax.plot(ts, a["Cd_st_s"], ".", ms=3, color="#ff8f00", label=f"Cd ∮σ·n, ⟨⟩={a['Cd_st']:.3f}")
|
||||
ax.plot(ts, a["Cl_st_s"], ".", ms=3, color="#4fc3f7", label="Cl ∮σ·n")
|
||||
ax.axvline(h0, color="gray", lw=0.8, ls="--")
|
||||
ax.set_xlabel("шаг L0"); ax.set_ylabel("Cd, Cl"); ax.legend(fontsize=8, loc="upper right")
|
||||
ax.set_title("Сила: GMEM (линии) vs ∮σ·n ds (точки)", fontsize=10)
|
||||
ax = axs[0, 1]
|
||||
if "Cm_s" in a:
|
||||
ax.plot(t, a["Cm_s"], lw=0.7, color="#388e3c")
|
||||
ax.axhline(0.0, color="gray", lw=0.8)
|
||||
ax.set_title(f"Cm(t): ⟨Cm⟩={a['Cm']:.5f} (контроль симметрии, ожид. ≈0)", fontsize=10)
|
||||
ax.set_xlabel("шаг L0"); ax.set_ylabel("Cm")
|
||||
ax = axs[1, 0]
|
||||
ax.plot(a["freqs"]*cfg.D/cfg.U, a["amp"], lw=0.9, color="#7b1fa2")
|
||||
ax.axvline(a["St"], color="#d32f2f", lw=1.0, ls="--", label=f"St={a['St']:.3f}")
|
||||
ax.set_xlim(0, 0.6); ax.set_xlabel("St = f·D/U"); ax.set_ylabel("|FFT(u_y)|")
|
||||
ax.legend(fontsize=9); ax.set_title("Спектр зонда в следе (вторая половина ряда)", fontsize=10)
|
||||
ax = axs[1, 1]
|
||||
rho = B.to_cpu(res["rho"]) if not isinstance(res["rho"], np.ndarray) else res["rho"]
|
||||
ax.plot(np.arange(len(rho)), rho, lw=0.9, color="#5d4037")
|
||||
ax.axhline(1.0, color="gray", lw=0.8)
|
||||
ax.set_xlabel("шаг L0"); ax.set_ylabel("⟨ρ⟩")
|
||||
ax.set_title(f"⟨ρ⟩(t): дрейф массы (финально {float(rho[-1]):.4f})", fontsize=10)
|
||||
fig.suptitle(f"2×+SDF · Ny={cfg.Ny} (β={cfg.D/cfg.Ny:.3f}) · Re={cfg.Re:g} · {len(rho)-1} шагов",
|
||||
fontsize=11)
|
||||
fig.tight_layout(rect=(0, 0, 1, 0.96))
|
||||
os.makedirs(os.path.dirname(path), exist_ok=True)
|
||||
fig.savefig(path, bbox_inches="tight"); plt.close(fig)
|
||||
print("saved", os.path.basename(path))
|
||||
|
||||
def save_gif(res, geo, cfg, path, fps=24):
|
||||
import matplotlib; matplotlib.use("Agg")
|
||||
import matplotlib.pyplot as plt, matplotlib.patches as mpatches
|
||||
from PIL import Image
|
||||
s0 = geo.solid0_np; s1 = geo.solid1_np
|
||||
Nx, Ny = cfg.Nx, cfg.Ny; ax1, bx1, ay1, by1 = cfg.ax1, cfg.bx1, cfg.ay1, cfg.by1
|
||||
frames = res["frames"]
|
||||
vm = np.nanpercentile(np.abs(np.where(s0, np.nan, _vort(frames[len(frames)//2][1]))), 99.0)
|
||||
cm = plt.get_cmap("inferno").copy(); cm.set_bad("white"); pil = []
|
||||
for (t, u0, u1) in frames:
|
||||
fig = plt.figure(figsize=(11.2, 5.7), dpi=96); ax = fig.add_subplot(111)
|
||||
ax.imshow(np.abs(np.where(s0, np.nan, _vort(u0))), origin="lower", cmap=cm, vmin=0, vmax=vm,
|
||||
extent=(0, Nx, 0, Ny), interpolation="nearest")
|
||||
ax.imshow(np.abs(np.where(s1, np.nan, _vort(u1)))[1:-1, 1:-1], origin="lower", cmap=cm, vmin=0, vmax=vm,
|
||||
extent=(ax1+0.5, bx1-0.5, ay1+0.5, by1-0.5), interpolation="nearest")
|
||||
ax.add_patch(mpatches.Rectangle((0.2, 0.2), Nx-0.4, Ny-0.4, fill=False, edgecolor="#4fc3f7", lw=1.6))
|
||||
ax.text(1.5, Ny-4, "L0: Δx", color="#4fc3f7", fontsize=9, weight="bold")
|
||||
ax.add_patch(mpatches.Rectangle((ax1, ay1), bx1-ax1, by1-ay1, fill=False, edgecolor="#aed581", lw=2.2))
|
||||
ax.text(ax1+1, by1-3.5, "L1: Δx/2 (тело+след)", color="#aed581", fontsize=9, weight="bold")
|
||||
ax.set_title(f"2×+SDF · |\\omega| · t={t}", fontsize=11)
|
||||
ax.set_xlabel("x [lu]"); ax.set_ylabel("y [lu]"); ax.set_xlim(0, Nx); ax.set_ylim(0, Ny)
|
||||
buf = __import__("io").BytesIO(); fig.savefig(buf, format="png", bbox_inches="tight"); buf.seek(0); plt.close(fig)
|
||||
pil.append(Image.open(buf).convert("RGB").convert("P", palette=Image.ADAPTIVE))
|
||||
os.makedirs(os.path.dirname(path), exist_ok=True)
|
||||
pil[0].save(path, save_all=True, append_images=pil[1:], duration=int(1000/fps), loop=0, optimize=True)
|
||||
print("saved", os.path.basename(path), len(frames), "кадров")
|
||||
Reference in New Issue
Block a user