NotBigGhostandClaude Opus 5
11ff7b79b4
Состояние на момент заведения репозитория. 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>
1644 lines
98 KiB
Python
1644 lines
98 KiB
Python
# ---
|
||
# jupyter:
|
||
# jupytext:
|
||
# text_representation:
|
||
# extension: .py
|
||
# format_name: percent
|
||
# format_version: '1.3'
|
||
# kernelspec:
|
||
# display_name: Python 3
|
||
# language: python
|
||
# name: python3
|
||
# ---
|
||
|
||
# %% [markdown]
|
||
# # Метод KBC (энтропийный Lattice Boltzmann): математика, физика и рабочий разбор
|
||
#
|
||
# **Karlin–Bösch–Chikatamarla (KBC)** — энтропийная multi-relaxation-time (MRT)
|
||
# формулировка решёточного метода Больцмана (Lattice Boltzmann Method, **LBM**),
|
||
# которая обеспечивает **безусловную нелинейную устойчивость** при больших числах
|
||
# Рейнольдса без явной модели подсеточной турбулентности (она работает как
|
||
# *implicit LES*).
|
||
#
|
||
# Этот ноутбук — подготовка к реализации: мы разбираем метод «на бумаге»,
|
||
# **выводим** ключевые формулы, **численно их самопроверяем** и строим рабочий
|
||
# 2D-солвер (решётка D2Q9) с валидацией. В конце — перенос на целевую решётку
|
||
# D3Q27 и сравнение с подходом «полиномиальное равновесие + TRT + α-лимитер +
|
||
# Smagorinsky», который применялся ранее (это и есть «что менять для перехода к
|
||
# настоящему KBC»).
|
||
#
|
||
# > **О источниках.** Формулы по возможности взяты из оригинальных статей
|
||
# > (Karlin, Bösch, Chikatamarla, *Phys. Rev. E* **90**, 031302(R), 2014;
|
||
# > Bösch, Chikatamarla, Karlin, *Phys. Rev. E* **92**, 043309, 2015; Guo et al.,
|
||
# > *Phys. Rev. E* **65**, 046308, 2002). Где доступ был только к вторичным
|
||
# > источникам, формула помечена **⚠ сверить с оригиналом**, а в конце собран
|
||
# > список таких мест. Главный страж корректности здесь — **численная
|
||
# > самопроверка** каждой формулы.
|
||
|
||
# %% [markdown]
|
||
# ## Сводная таблица обозначений
|
||
#
|
||
# Везде используются **решёточные единицы** (lattice units): шаг сетки
|
||
# $\Delta x = 1$ (lu — *lattice unit*) и шаг по времени $\Delta t = 1$ (ts — *time
|
||
# step*). Переход к СИ описан в разделе «Единицы и обезразмеривание».
|
||
#
|
||
# | Символ | Название | Физический смысл | Единицы (lu/ts) | Единицы (СИ) | Диапазон |
|
||
# |---|---|---|---|---|---|
|
||
# | $f_i$ | популяция | плотность частиц со скоростью $\mathbf c_i$ | $f_i$ безразм. (доля $\rho$) | кг·с³/м⁶ | $f_i>0$ |
|
||
# | $\mathbf c_i$ | дискретная скорость | направление переноса $i$ | lu/ts | м/с | компоненты $\in\{-1,0,1\}$ |
|
||
# | $w_i$ | вес квадратуры | вес направления, $\sum_i w_i=1$ | безразм. | безразм. | $\{4/9,1/9,1/36\}$ (D2Q9) |
|
||
# | $c_s$ | скорость звука решётки | $c_s^2=1/3$ для D2Q9/D3Q27 | lu/ts | м/с | $c_s=1/\sqrt3$ |
|
||
# | $\rho$ | плотность | $\rho=\sum_i f_i$ | безразм. (опорная 1) | кг/м³ | $\approx 1$ |
|
||
# | $\mathbf u$ | скорость | $\rho\mathbf u=\sum_i \mathbf c_i f_i$ | lu/ts | м/с | $\lvert\mathbf u\rvert\ll c_s$ |
|
||
# | $\mathbf j$ | импульс | $\mathbf j=\rho\mathbf u$ | — | — | — |
|
||
# | $\tau$ | время релаксации | задаёт вязкость | ts | с | $\tau>1/2$ |
|
||
# | $\nu$ | кинематич. вязкость | $\nu=c_s^2(\tau-\tfrac12)$ | lu²/ts | м²/с | $>0$ |
|
||
# | $\beta$ | параметр сдвиговой релаксации | $\beta=1/(2\tau)$ | безразм. | — | $0<\beta<1$ |
|
||
# | $\gamma$ | **энтропийный стабилизатор** | релаксация высших мод | безразм. | — | $\approx 2$ |
|
||
# | $f_i^{\mathrm{eq}}$ | равновесие | максимум энтропии при связях $\rho,\mathbf u$ | как $f_i$ | — | $>0$ |
|
||
# | $k_i,s_i,h_i$ | части разложения | кинематич./сдвиг/высшие моды | как $f_i$ | — | — |
|
||
# | $\Delta s_i,\Delta h_i$ | неравновесные части | $s_i-s_i^{eq}$, $h_i-h_i^{eq}$ | как $f_i$ | — | малы |
|
||
# | $H$ | H-функция (энтропия) | $H=\sum_i f_i\ln(f_i/w_i)$ | безразм. | — | убывает |
|
||
# | $\langle X|Y\rangle$ | энтропийное скал. произв. | $\sum_i X_iY_i/f_i^{eq}$ | — | — | — |
|
||
# | $\mathrm{Ma}$ | число Маха | $\mathrm{Ma}=\lvert\mathbf u\rvert/c_s$ | безразм. | — | $\lesssim 0.3$ |
|
||
# | $\mathrm{Re}$ | число Рейнольдса | $\mathrm{Re}=UL/\nu$ | безразм. | — | задача |
|
||
# | $S_{\alpha\beta}$ | тензор скоростей деформации | $\tfrac12(\partial_\alpha u_\beta+\partial_\beta u_\alpha)$ | 1/ts | 1/с | — |
|
||
# | $\nu_t$ | вихревая вязкость (LES) | $\nu_t=(C_s\Delta)^2\lvert S\rvert$ | lu²/ts | м²/с | $\ge 0$ |
|
||
# | $C_s$ | константа Смагоринского | — | безразм. | — | $0.1$–$0.17$ |
|
||
# | $q$ | доля стенки | положение стенки между узлами | безразм. | — | $[0,1]$ |
|
||
|
||
# %% [markdown]
|
||
# ## 0. Стиль графиков и общие настройки
|
||
#
|
||
# Подключаем библиотеки и задаём единый «симпатичный» стиль matplotlib. Код
|
||
# написан так, чтобы файл работал **и как ноутбук** (графики встраиваются), **и
|
||
# как обычный скрипт** `python kbc_lbm.py` (графики сохраняются в `figures/` для
|
||
# валидации физики).
|
||
|
||
# %%
|
||
import numpy as np
|
||
import matplotlib
|
||
|
||
# Определяем, выполняемся ли мы внутри Jupyter. В скрипте используем неинтерактивный
|
||
# бэкенд Agg и сохраняем PNG; в ноутбуке полагаемся на inline-вывод nbconvert.
|
||
try:
|
||
get_ipython() # есть только в IPython/Jupyter
|
||
IN_NOTEBOOK = True
|
||
except NameError:
|
||
IN_NOTEBOOK = False
|
||
if not IN_NOTEBOOK:
|
||
matplotlib.use("Agg")
|
||
# Консоль Windows может быть в cp1251 — выводим в UTF-8, чтобы кириллица и
|
||
# греческие символы в print() не роняли скрипт (в ноутбуке это не нужно).
|
||
import sys
|
||
try:
|
||
sys.stdout.reconfigure(encoding="utf-8")
|
||
except Exception:
|
||
pass
|
||
|
||
import os
|
||
import matplotlib.pyplot as plt
|
||
|
||
FIGDIR = os.path.join(os.path.dirname(os.path.abspath(__file__)) if "__file__" in globals() else ".", "figures")
|
||
os.makedirs(FIGDIR, exist_ok=True)
|
||
|
||
plt.rcParams.update({
|
||
"figure.dpi": 110,
|
||
"savefig.dpi": 120,
|
||
"font.size": 11,
|
||
"axes.titlesize": 13,
|
||
"axes.titleweight": "bold",
|
||
"axes.grid": True,
|
||
"grid.alpha": 0.25,
|
||
"axes.axisbelow": True,
|
||
"figure.facecolor": "white",
|
||
"axes.facecolor": "#fbfbfd",
|
||
})
|
||
|
||
def _finish(fig, name):
|
||
"""Сохранить фигуру в figures/ (для скрипта) и показать (в ноутбуке)."""
|
||
fig.savefig(os.path.join(FIGDIR, name), bbox_inches="tight")
|
||
if IN_NOTEBOOK:
|
||
plt.show() # гарантированный inline-вывод при nbconvert
|
||
else:
|
||
plt.close(fig)
|
||
|
||
print("matplotlib", matplotlib.__version__, "| numpy", np.__version__, "| notebook:", IN_NOTEBOOK)
|
||
|
||
|
||
# %% [markdown]
|
||
# ## 1. Основы LBM и решётка D2Q9
|
||
#
|
||
# Решёточный метод Больцмана отслеживает **популяции** $f_i(\mathbf x, t)$ —
|
||
# плотность фиктивных частиц, движущихся в узле $\mathbf x$ со скоростью
|
||
# $\mathbf c_i$. Динамика — чередование **переноса** (streaming) и **столкновения**
|
||
# (collision):
|
||
# $$ f_i(\mathbf x + \mathbf c_i\,\Delta t,\; t+\Delta t) = f_i(\mathbf x,t) + \Omega_i, $$
|
||
# где $\Omega_i$ — оператор столкновения.
|
||
#
|
||
# **Решётка D2Q9** (2D, 9 скоростей): покоящаяся частица, 4 «осевые» и 4
|
||
# «диагональные». Скорость звука $c_s^2 = 1/3$. Веса $w_i$ выбраны так, чтобы
|
||
# квадратура точно интегрировала моменты максвеллиана до нужного порядка
|
||
# (изотропия 4-го порядка): $w_0=4/9$, осевые $1/9$, диагональные $1/36$,
|
||
# $\sum_i w_i = 1$.
|
||
#
|
||
# **Макроскопические моменты:**
|
||
# $$ \rho=\sum_i f_i,\qquad \rho\mathbf u=\sum_i \mathbf c_i f_i. $$
|
||
|
||
# %%
|
||
# --- Решётка D2Q9 -----------------------------------------------------------
|
||
# Порядок направлений: 0 — покой; 1..4 — оси (+x,+y,-x,-y); 5..8 — диагонали.
|
||
C = np.array([[0, 0], [1, 0], [0, 1], [-1, 0], [0, -1],
|
||
[1, 1], [-1, 1], [-1, -1], [1, -1]], dtype=float) # c_i = (c_ix, c_iy)
|
||
W = np.array([4/9, 1/9, 1/9, 1/9, 1/9, 1/36, 1/36, 1/36, 1/36]) # веса w_i
|
||
CS2 = 1.0 / 3.0 # c_s^2
|
||
Q = 9
|
||
CX, CY = C[:, 0], C[:, 1]
|
||
OPP = np.array([0, 3, 4, 1, 2, 7, 8, 5, 6]) # индекс противоположного c_i
|
||
|
||
# Самопроверка решёточных тождеств (нужны для корректной гидродинамики):
|
||
assert np.isclose(W.sum(), 1.0) # сумма весов = 1
|
||
assert np.allclose((W[:, None] * C).sum(0), 0.0) # sum_i w_i c_i = 0
|
||
assert np.allclose((W[:, None, None] * C[:, :, None] * C[:, None, :]).sum(0),
|
||
CS2 * np.eye(2)) # sum_i w_i c_ia c_ib = c_s^2 δ_ab
|
||
assert np.all(C[OPP] == -C) # c_{opp(i)} = -c_i
|
||
print("Решётка D2Q9: тождества весов и изотропии 2-го порядка — OK")
|
||
|
||
|
||
# %%
|
||
def make_grid(field_shape):
|
||
"""Удобный конструктор массива популяций формы (9, Ny, Nx)."""
|
||
return np.zeros((Q,) + field_shape)
|
||
|
||
|
||
def macros(f):
|
||
"""Плотность rho (Ny,Nx) и скорость u (2,Ny,Nx) из популяций f (9,Ny,Nx)."""
|
||
rho = f.sum(axis=0)
|
||
jx = (CX[:, None, None] * f).sum(axis=0)
|
||
jy = (CY[:, None, None] * f).sum(axis=0)
|
||
u = np.stack([jx / rho, jy / rho])
|
||
return rho, u
|
||
|
||
|
||
# %% [markdown]
|
||
# ### Рисунок: стенсиль D2Q9
|
||
#
|
||
# Стрелки — направления $\mathbf c_i$, цвет/толщина — вес $w_i$ (центральный узел
|
||
# имеет наибольший вес $4/9$, диагонали — наименьший $1/36$).
|
||
|
||
# %%
|
||
fig, ax = plt.subplots(figsize=(5.2, 5.2))
|
||
tier_color = {4/9: "#d1495b", 1/9: "#2e7d32", 1/36: "#1f6feb"}
|
||
for i in range(Q):
|
||
col = tier_color[W[i]]
|
||
if i == 0:
|
||
ax.plot(0, 0, "o", ms=18, color=col, zorder=3)
|
||
ax.annotate(f"$w_0={W[i]:.3f}$", (0, 0), textcoords="offset points",
|
||
xytext=(8, 8), fontsize=9)
|
||
else:
|
||
ax.annotate("", xy=(CX[i], CY[i]), xytext=(0, 0),
|
||
arrowprops=dict(arrowstyle="-|>", lw=1.5 + 40 * W[i], color=col))
|
||
ax.annotate(f"$c_{{{i}}}$", (CX[i], CY[i]), textcoords="offset points",
|
||
xytext=(6 * np.sign(CX[i] + 0.1), 6 * np.sign(CY[i] + 0.1)), fontsize=9)
|
||
ax.set_xlim(-1.6, 1.6); ax.set_ylim(-1.6, 1.6); ax.set_aspect("equal")
|
||
ax.set_title("Решётка D2Q9: дискретные скорости $\\mathbf{c}_i$ и веса $w_i$")
|
||
ax.set_xlabel("$c_x$ [lu/ts]"); ax.set_ylabel("$c_y$ [lu/ts]")
|
||
from matplotlib.lines import Line2D
|
||
ax.legend(handles=[Line2D([0], [0], color=c, lw=3, label=f"$w={k:.4f}$")
|
||
for k, c in tier_color.items()], loc="upper right", fontsize=8)
|
||
_finish(fig, "01_stencil_d2q9.png")
|
||
|
||
|
||
# %% [markdown]
|
||
# ## 2. Равновесие: полиномиальное и энтропийное
|
||
#
|
||
# **Полиномиальное равновесие** (разложение максвеллиана по Эрмиту до $O(u^2)$),
|
||
# стандарт для BGK:
|
||
# $$ f_i^{\mathrm{eq,poly}} = w_i\,\rho\left[1 + \frac{\mathbf c_i\cdot\mathbf u}{c_s^2}
|
||
# + \frac{(\mathbf c_i\cdot\mathbf u)^2}{2c_s^4} - \frac{\mathbf u^2}{2c_s^2}\right]. $$
|
||
# Коэффициенты при $c_s^2=1/3$: $1/c_s^2=3$, $1/(2c_s^4)=4.5$, $1/(2c_s^2)=1.5$.
|
||
#
|
||
# **Энтропийное равновесие** (Ansumali–Karlin): точный максимум дискретной
|
||
# H-функции при связях $\sum f^{eq}=\rho$, $\sum \mathbf c f^{eq}=\rho\mathbf u$.
|
||
# Для D2Q9/D3Q27 имеет **замкнутую факторизованную форму** ⚠ (сверить с
|
||
# оригиналом):
|
||
# $$ f_i^{\mathrm{eq}} = \rho\,w_i\prod_{\alpha\in\{x,y\}}
|
||
# \underbrace{\left(2-\sqrt{1+3u_\alpha^2}\right)}_{\text{префактор}}
|
||
# \left(\frac{2u_\alpha+\sqrt{1+3u_\alpha^2}}{1-u_\alpha}\right)^{c_{i\alpha}}. $$
|
||
# Здесь $3 = 1/c_s^2$. Энтропийная форма **точно** воспроизводит $\rho$ и
|
||
# $\rho\mathbf u$ (проверим ниже) и, в отличие от полинома, гарантирует
|
||
# положительность и H-теорему.
|
||
|
||
# %%
|
||
def feq_poly(rho, u):
|
||
"""Полиномиальное равновесие O(u^2). rho:(Ny,Nx), u:(2,Ny,Nx) -> (9,Ny,Nx)."""
|
||
cu = CX[:, None, None] * u[0] + CY[:, None, None] * u[1] # c_i · u
|
||
usq = u[0] ** 2 + u[1] ** 2
|
||
return W[:, None, None] * rho * (1 + cu / CS2 + cu ** 2 / (2 * CS2 ** 2) - usq / (2 * CS2))
|
||
|
||
|
||
def feq_entropic(rho, u):
|
||
"""Энтропийное равновесие (product-form). Клип скорости для численной безопасности."""
|
||
ux = np.clip(u[0], -0.95, 0.95)
|
||
uy = np.clip(u[1], -0.95, 0.95)
|
||
sx = np.sqrt(1 + 3 * ux ** 2)
|
||
sy = np.sqrt(1 + 3 * uy ** 2)
|
||
prefx, prefy = 2 - sx, 2 - sy # префактор по каждому измерению
|
||
Bx = (2 * ux + sx) / (1 - ux) # основание степени по x
|
||
By = (2 * uy + sy) / (1 - uy)
|
||
base = rho * prefx * prefy # (Ny,Nx)
|
||
Bxp = Bx[None] ** CX[:, None, None] # (9,Ny,Nx): B_x^{c_ix}
|
||
Byp = By[None] ** CY[:, None, None] # (9,Ny,Nx): B_y^{c_iy}
|
||
return W[:, None, None] * base[None] * Bxp * Byp
|
||
|
||
|
||
# %% [markdown]
|
||
# ### Самопроверка равновесий
|
||
#
|
||
# Проверяем, что **энтропийное** равновесие точно даёт нужные нулевой и первый
|
||
# моменты ($\rho$, $\rho\mathbf u$), а второй момент при малом $\mathrm{Ma}$ близок
|
||
# к $\rho(c_s^2\delta_{\alpha\beta}+u_\alpha u_\beta)$ (тензор давления Навье–Стокса).
|
||
|
||
# %%
|
||
rng = np.random.default_rng(0)
|
||
rho_t = np.ones((1, 1)) * (1.0 + 0.1 * rng.standard_normal((1, 1)))
|
||
u_t = np.stack([0.06 * rng.standard_normal((1, 1)), 0.06 * rng.standard_normal((1, 1))])
|
||
fe = feq_entropic(rho_t, u_t)
|
||
rho_m, u_m = macros(fe)
|
||
Pxx = (CX[:, None, None] ** 2 * fe).sum(0)
|
||
err_rho = abs(rho_m - rho_t).max()
|
||
err_u = abs(u_m - u_t).max()
|
||
err_P = abs(Pxx - rho_t * (CS2 + u_t[0] ** 2)).max()
|
||
print(f"Энтропийное равновесие: |Σf-ρ|={err_rho:.2e}, |Σcf-ρu|={err_u:.2e} (точно),")
|
||
print(f" |Σc_x²f - ρ(c_s²+u_x²)|={err_P:.2e} (мало при малом Ma)")
|
||
assert err_rho < 1e-12 and err_u < 1e-12
|
||
print("Моменты равновесия — OK")
|
||
|
||
|
||
# %% [markdown]
|
||
# ## 3. Разложение популяций $f = k + s + h$ через моментные проекторы
|
||
#
|
||
# Сердце KBC — разбиение каждого 9-вектора популяций на три части по **моментным
|
||
# подпространствам**:
|
||
# - $k_i$ — **кинематическая** часть, несёт только сохраняющиеся моменты
|
||
# ($\rho$, $j_x$, $j_y$);
|
||
# - $s_i$ — **сдвиговая** часть, несёт девиаторный тензор напряжений
|
||
# ($N=\Pi_{xx}-\Pi_{yy}$ и $\Pi_{xy}$);
|
||
# - $h_i$ — **высшие** («ghost») моды: след 2-го момента (объёмная вязкость),
|
||
# 3-й и 4-й моменты.
|
||
#
|
||
# Строим базис из 9 независимых мономов-моментов $\phi^{(a)}_i$ и матрицу
|
||
# $M_{ai}=\phi^{(a)}(\mathbf c_i)$, тогда $m = M f$ — вектор моментов. Проектор на
|
||
# группу моментов $G$: $\;P_G = M^{-1} D_G M$, где $D_G$ — диагональная маска
|
||
# (единицы на моментах из $G$). По построению $P_k+P_s+P_h = I$, поэтому
|
||
# $k+s+h=f$ **точно**. Назначение моментов в $s$ vs $h$ задаёт *вариант* KBC
|
||
# (ниже — KBC-N1: в $s$ только девиаторный стресс).
|
||
|
||
# %%
|
||
def moment_matrix():
|
||
"""Строки — мономы-моменты, столбцы — направления. m = M @ f."""
|
||
M = np.zeros((Q, Q))
|
||
M[0] = 1.0 # ρ
|
||
M[1] = CX # j_x
|
||
M[2] = CY # j_y
|
||
M[3] = 3 * (CX ** 2 + CY ** 2) - 2 # след 2-го момента (энергия), центрированный
|
||
M[4] = CX ** 2 - CY ** 2 # N — разность нормальных напряжений
|
||
M[5] = CX * CY # Π_xy — сдвиговое напряжение
|
||
M[6] = CX ** 2 * CY # 3-й момент
|
||
M[7] = CX * CY ** 2 # 3-й момент
|
||
M[8] = CX ** 2 * CY ** 2 # 4-й момент
|
||
return M
|
||
|
||
M = moment_matrix()
|
||
Minv = np.linalg.inv(M) # существует: 9 мономов независимы на D2Q9
|
||
|
||
K_MOMENTS = [0, 1, 2] # сохраняющиеся: ρ, j_x, j_y -> k
|
||
S_MOMENTS = [4, 5] # девиаторный стресс: N, Π_xy -> s (вариант KBC-N1)
|
||
H_MOMENTS = [3, 6, 7, 8] # след + 3-й + 4-й моменты -> h
|
||
|
||
def projector(moment_idx):
|
||
D = np.zeros((Q, Q))
|
||
D[moment_idx, moment_idx] = 1.0
|
||
return Minv @ D @ M
|
||
|
||
P_k, P_s, P_h = projector(K_MOMENTS), projector(S_MOMENTS), projector(H_MOMENTS)
|
||
|
||
# Самопроверка: полнота (P_k+P_s+P_h=I) и что k несёт ровно сохраняющиеся моменты.
|
||
assert np.allclose(P_k + P_s + P_h, np.eye(Q))
|
||
ftest = rng.random((Q, 3, 4))
|
||
k = np.einsum("ij,jab->iab", P_k, ftest)
|
||
s = np.einsum("ij,jab->iab", P_s, ftest)
|
||
h = np.einsum("ij,jab->iab", P_h, ftest)
|
||
assert np.allclose(k + s + h, ftest) # разложение полное
|
||
# у части k моменты N, Π и высшие — нулевые (несёт только ρ, j_x, j_y):
|
||
mk = np.einsum("ai,ijl->ajl", M, k) # моменты части k
|
||
assert np.allclose(mk[S_MOMENTS], 0.0) and np.allclose(mk[H_MOMENTS], 0.0)
|
||
print("Разложение k+s+h=f: макс. невязка =", abs(k + s + h - ftest).max())
|
||
print("Проекторы P_k+P_s+P_h=I и состав части k — OK")
|
||
|
||
|
||
# %% [markdown]
|
||
# ## 4. Оператор столкновений KBC и стабилизатор $\gamma$
|
||
#
|
||
# Поскольку $f$ и $f^{eq}$ имеют **одинаковые** сохраняющиеся моменты, неравновесие
|
||
# $\,f-f^{eq}\,$ не содержит кинематической части: $\Delta k = 0$, значит
|
||
# $f-f^{eq}=\Delta s+\Delta h$, где $\Delta s = P_s(f-f^{eq})$,
|
||
# $\Delta h = P_h(f-f^{eq})$.
|
||
#
|
||
# **Столкновение KBC:**
|
||
# $$ f_i^{\text{после}} = f_i - \beta\,\bigl(2\,\Delta s_i + \gamma\,\Delta h_i\bigr). $$
|
||
# Множитель **2** перед $\Delta s$ фиксирует сдвиговую вязкость через
|
||
# $\beta=1/(2\tau)$; высшие моды релаксируются отдельным $\gamma$.
|
||
#
|
||
# **Стабилизатор $\gamma$** выбирается локально из условия неувеличения дискретной
|
||
# энтропии ($\Delta H = 0$ во 2-м порядке). Результат — **замкнутая формула**
|
||
# ⚠ (сверить с оригиналом):
|
||
# $$ \gamma = \frac{1}{\beta} - \Bigl(2-\frac{1}{\beta}\Bigr)
|
||
# \frac{\langle \Delta s\,|\,\Delta h\rangle}{\langle \Delta h\,|\,\Delta h\rangle},
|
||
# \qquad \langle X|Y\rangle = \sum_i \frac{X_i Y_i}{f_i^{eq}}. $$
|
||
# При $\gamma=2$: $2\Delta s+2\Delta h = 2(f-f^{eq})$, и KBC **сводится к BGK**
|
||
# с $\omega=2\beta$: $f^{\text{после}}=f-2\beta(f-f^{eq})$. Адаптивный $\gamma$
|
||
# добавляет/убавляет диссипацию высших мод там, где это нужно для устойчивости.
|
||
|
||
# %%
|
||
def collide_kbc(f, feq, beta, gamma_fixed=None):
|
||
"""KBC-столкновение. Возвращает (f_после, поле gamma)."""
|
||
dfneq = f - feq
|
||
ds = np.einsum("ij,jab->iab", P_s, dfneq) # Δs
|
||
dh = dfneq - ds # Δh = (f-feq) - Δs, т.к. Δk=0
|
||
if gamma_fixed is not None:
|
||
gamma = np.full(f.shape[1:], float(gamma_fixed))
|
||
else:
|
||
inv = 1.0 / feq
|
||
num = (ds * dh * inv).sum(0) # ⟨Δs|Δh⟩
|
||
den = (dh * dh * inv).sum(0) # ⟨Δh|Δh⟩
|
||
with np.errstate(divide="ignore", invalid="ignore"):
|
||
gamma = np.where(den > 1e-14,
|
||
1.0 / beta - (2.0 - 1.0 / beta) * num / den,
|
||
2.0) # у равновесия (Δh→0) берём предел BGK
|
||
fpost = f - beta * (2.0 * ds + gamma[None] * dh)
|
||
return fpost, gamma
|
||
|
||
|
||
def collide_bgk(f, feq, beta):
|
||
"""Классический BGK: f - ω(f-feq), ω=2β."""
|
||
return f - 2.0 * beta * (f - feq)
|
||
|
||
|
||
# %% [markdown]
|
||
# ### Самопроверка: предел $\gamma=2$ совпадает с BGK
|
||
# При $\gamma=2$ оператор KBC обязан в точности равняться $f-2\beta(f-f^{eq})$.
|
||
|
||
# %%
|
||
rho0 = np.ones((6, 7))
|
||
u0 = np.stack([0.05 * rng.standard_normal((6, 7)), 0.05 * rng.standard_normal((6, 7))])
|
||
feq0 = feq_entropic(rho0, u0)
|
||
f0 = feq0 + 0.01 * rng.standard_normal((Q, 6, 7)) # слегка неравновесное состояние
|
||
f0 *= rho0 / f0.sum(0) # сохранить ρ
|
||
beta0 = 0.7
|
||
fk2, _ = collide_kbc(f0, feq_entropic(*macros(f0)), beta0, gamma_fixed=2.0)
|
||
fbgk = collide_bgk(f0, feq_entropic(*macros(f0)), beta0)
|
||
diff = abs(fk2 - fbgk).max()
|
||
# diff — машинный ноль (масштаб |f|·eps ≈ 1e-18); важно лишь, что он ≪ порога 1e-12.
|
||
print(f"max|KBC(γ=2) - BGK| = {diff:.1e} → машинный ноль (порог assert 1e-12), операторы совпадают")
|
||
assert diff < 1e-12
|
||
print("Предел γ=2 → BGK — OK")
|
||
|
||
|
||
# %% [markdown]
|
||
# ## 4b. Популярно: что делает каждая формула (математика + физика)
|
||
#
|
||
# Здесь каждая формула KBC объясняется «на пальцах»: что это математически, какая за
|
||
# этим физика потока, и как это считается на каждом шаге.
|
||
#
|
||
# ---
|
||
# **Популяции $f_i(\mathbf x,t)$ и шаг «перенос + столкновение».**
|
||
# *Математика:* в каждом узле сетки храним 9 чисел — «сколько газа летит» по каждому из
|
||
# 9 направлений $\mathbf c_i$. *Физика:* как 9 встречных потоков на перекрёстке.
|
||
# *Как считается:* за шаг потоки сначала **разъезжаются** в соседние узлы (перенос,
|
||
# `np.roll`), затем в каждом узле **перемешиваются** (столкновение). Перенос точный и
|
||
# дешёвый; вся «физика жидкости» сидит в столкновении.
|
||
#
|
||
# ---
|
||
# **Моменты $\rho=\sum_i f_i$, $\rho\mathbf u=\sum_i \mathbf c_i f_i$.**
|
||
# *Математика:* сумма всех потоков = плотность; взвешенная по направлениям = импульс.
|
||
# *Физика:* «сколько вещества» и «куда оно в среднем движется». *Как считается:* два
|
||
# суммирования по 9 числам — мгновенно. Отсюда привычные $\rho$ и скорость $\mathbf u$.
|
||
#
|
||
# ---
|
||
# **Равновесие $f_i^{\mathrm{eq}}(\rho,\mathbf u)$.**
|
||
# *Математика:* распределение по направлениям, максимизирующее энтропию при заданных
|
||
# $\rho,\mathbf u$. *Физика:* «к чему стремится газ, если оставить в покое» — локальный
|
||
# максвелловский «штиль». *Как считается:* по замкнутой формуле из $\rho,\mathbf u$.
|
||
# Столкновение всегда **тянет** $f$ к $f^{eq}$ — это и есть вязкое стремление к равновесию.
|
||
#
|
||
# ---
|
||
# **H-функция $H=\sum_i f_i\ln(f_i/w_i)$.**
|
||
# *Математика:* дискретный аналог энтропии. *Физика:* мера беспорядка/необратимости; по
|
||
# второму началу не должна самопроизвольно убывать. *Как считается:* одно суммирование.
|
||
# В KBC $H$ — **физический компас устойчивости**: из условия «энтропия не падает»
|
||
# выводится стабилизатор $\gamma$.
|
||
#
|
||
# ---
|
||
# **Разложение $f = k + s + h$.**
|
||
# *Математика:* проекция 9-вектора на три группы моментов. *Физика:* три «слоя» движения:
|
||
# - $k$ — **перенос массы и импульса** (поток как целое);
|
||
# - $s$ — **вязкий сдвиг** (трение слоёв друг о друга — задаёт вязкость);
|
||
# - $h$ — **мелкая «рябь»** высших мод (мелкомасштабные флуктуации, где копится
|
||
# неустойчивость на грубой сетке).
|
||
#
|
||
# *Как считается:* умножение на заранее посчитанные проекторы $P_s,P_h$. Идея KBC:
|
||
# релаксировать эти слои **с разной скоростью**.
|
||
#
|
||
# ---
|
||
# **Столкновение $f \leftarrow f - \beta\,(2\,\Delta s + \gamma\,\Delta h)$.**
|
||
# *Математика:* $\Delta s,\Delta h$ — отклонения сдвига и ряби от равновесия. *Физика:*
|
||
# **вязкий сдвиг** гасим фиксированным темпом $2\beta$ (он задаёт вязкость и Re); **рябь** —
|
||
# отдельным темпом $\gamma\beta$, ровно чтобы не нарушить энтропию. *Как считается:* две
|
||
# линейные комбинации 9-векторов на узел. Множитель **2** перед $\Delta s$ воспроизводит
|
||
# правильную вязкость; при $\gamma=2$ всё сворачивается в обычный BGK.
|
||
#
|
||
# ---
|
||
# **Стабилизатор $\gamma = \dfrac{1}{\beta} - \Bigl(2-\dfrac{1}{\beta}\Bigr)
|
||
# \dfrac{\langle \Delta s|\Delta h\rangle}{\langle \Delta h|\Delta h\rangle}$.**
|
||
# *Математика:* единственное $\gamma$, при котором изменение энтропии за шаг равно нулю
|
||
# (во 2-м порядке). *Физика:* **локальный термостат потока** в реальном времени — в каждой
|
||
# точке на каждом шаге подкручивает, насколько гасить «рябь», чтобы энтропия не росла. В
|
||
# спокойных зонах $\gamma\approx 2$ (как BGK); у резких градиентов (кромки тел, вихри)
|
||
# $\gamma$ отклоняется и добавляет ровно столько диссипации, сколько нужно, **не давая
|
||
# симуляции взорваться**. Это «implicit LES» — встроенная, без подгоночных констант,
|
||
# модель мелкомасштабной турбулентности. *Как считается:* два скалярных произведения на
|
||
# узел → одно число $\gamma$.
|
||
#
|
||
# ---
|
||
# **Энтропийное скалярное произведение $\langle X|Y\rangle=\sum_i X_iY_i/f_i^{eq}$.**
|
||
# *Математика:* взвешенное (на $1/f^{eq}$) скалярное произведение — кривизна энтропии у
|
||
# равновесия. *Физика:* измеряет, насколько «рябь» $\Delta h$ выровнена со сдвигом
|
||
# $\Delta s$; по этому решается, сколько диссипации добавить. *Как считается:* поэлементное
|
||
# произведение и сумма по 9 направлениям.
|
||
#
|
||
# ---
|
||
# **Вязкость $\nu = c_s^2\bigl(\tfrac{1}{2\beta}-\tfrac12\bigr)$, $\beta=1/(2\tau)$.**
|
||
# *Математика:* связь темпа релаксации сдвига с вязкостью. *Физика:* чем медленнее гасится
|
||
# сдвиг (меньше $\beta$), тем «гуще» жидкость; через $\nu$ задаётся $\mathrm{Re}=UL/\nu$.
|
||
# *Как считается:* один раз при настройке (выбрали $\nu$ под нужный Re → получили $\beta$).
|
||
|
||
|
||
# %% [markdown]
|
||
# ## 5. Перенос (streaming) и связь вязкости с $\beta$
|
||
#
|
||
# Перенос — сдвиг каждой популяции на её вектор $\mathbf c_i$ (здесь — периодически
|
||
# через `np.roll`). Вязкость: $\nu = c_s^2\bigl(\tfrac{1}{2\beta}-\tfrac12\bigr)$,
|
||
# то есть $\tau=\nu/c_s^2+\tfrac12$, $\beta=1/(2\tau)$. Число Рейнольдса
|
||
# $\mathrm{Re}=UL/\nu$; число Маха $\mathrm{Ma}=U/c_s$ держим $\lesssim 0.1$.
|
||
|
||
# %%
|
||
def stream(f):
|
||
"""Периодический перенос популяций на c_i. f:(9,Ny,Nx)."""
|
||
fs = np.empty_like(f)
|
||
for i in range(Q):
|
||
fs[i] = np.roll(f[i], shift=(int(CY[i]), int(CX[i])), axis=(0, 1))
|
||
return fs
|
||
|
||
|
||
def nu_to_beta(nu):
|
||
tau = nu / CS2 + 0.5
|
||
return 1.0 / (2.0 * tau)
|
||
|
||
|
||
def H_function(f):
|
||
"""Дискретная H-функция (энтропия) H = Σ f_i ln(f_i/w_i)."""
|
||
fp = np.maximum(f, 1e-300)
|
||
return (fp * np.log(fp / W[:, None, None])).sum(0)
|
||
|
||
|
||
print("Пример: ν=0.01 → β =", round(nu_to_beta(0.01), 4), ", τ =", round(0.01 / CS2 + 0.5, 4))
|
||
|
||
|
||
# %% [markdown]
|
||
# ## 6. Алгоритм одного шага KBC
|
||
#
|
||
# ```
|
||
# для каждого узла:
|
||
# 1) макро-моменты: rho = sum f_i, rho*u = sum c_i f_i
|
||
# 2) равновесие: f_i^eq (энтропийное product-form)
|
||
# 3) неравновесие: Ds = P_s (f - f^eq), Dh = (f - f^eq) - Ds
|
||
# 4) стабилизатор: gamma = 1/beta - (2 - 1/beta) <Ds|Dh> / <Dh|Dh>
|
||
# 5) столкновение: f <- f - beta (2 Ds + gamma Dh)
|
||
# 6) перенос: f_i(x + c_i) <- f_i(x)
|
||
# 7) граничные условия (bounce-back / Zou-He / sponge)
|
||
# ```
|
||
|
||
# %%
|
||
fig, ax = plt.subplots(figsize=(4.3, 6.2))
|
||
ax.set_axis_off()
|
||
steps_txt = [
|
||
("Макро-моменты\n$\\rho=\\sum f_i,\\ \\rho\\mathbf{u}=\\sum \\mathbf{c}_i f_i$", "#e8f0fe"),
|
||
("Равновесие $f_i^{eq}$\n(энтропийное)", "#e6f4ea"),
|
||
("Неравновесие\n$\\Delta s=P_s(f-f^{eq}),\\ \\Delta h=(f-f^{eq})-\\Delta s$", "#fff4e5"),
|
||
("Стабилизатор $\\gamma$\n$=1/\\beta-(2-1/\\beta)\\frac{\\langle\\Delta s|\\Delta h\\rangle}{\\langle\\Delta h|\\Delta h\\rangle}$", "#fde7ef"),
|
||
("Столкновение\n$f\\leftarrow f-\\beta(2\\Delta s+\\gamma\\Delta h)$", "#f3e8fd"),
|
||
("Перенос\n$f_i(\\mathbf{x}+\\mathbf{c}_i)\\leftarrow f_i(\\mathbf{x})$", "#e8f0fe"),
|
||
("Граничные условия", "#eceff1"),
|
||
]
|
||
y = 0.95
|
||
for txt, col in steps_txt:
|
||
ax.add_patch(plt.Rectangle((0.05, y - 0.1), 0.9, 0.085, transform=ax.transAxes,
|
||
facecolor=col, edgecolor="#555", lw=1.2))
|
||
ax.text(0.5, y - 0.057, txt, transform=ax.transAxes, ha="center", va="center", fontsize=8.5)
|
||
if y > 0.2:
|
||
ax.annotate("", xy=(0.5, y - 0.118), xytext=(0.5, y - 0.1), xycoords="axes fraction",
|
||
arrowprops=dict(arrowstyle="-|>", color="#555"))
|
||
y -= 0.13
|
||
ax.set_title("Шаг KBC")
|
||
_finish(fig, "06_algorithm_flow.png")
|
||
|
||
|
||
# %% [markdown]
|
||
# ## 7. Демо-1: вихрь Тейлора–Грина (валидация вязкости и H-теоремы)
|
||
#
|
||
# Точное решение несжимаемого Навье–Стокса в периодической области:
|
||
# $$ u_x=-U_0\cos(k x)\sin(k y),\quad u_y=U_0\sin(k x)\cos(k y),\quad k=2\pi/N. $$
|
||
# Кинетическая энергия затухает как $E(t)=E_0\,e^{-4\nu k^2 t}$. Прогоняем KBC и
|
||
# **измеряем** $\nu$ по наклону $\ln E(t)$ — он должен совпасть с заданным $\nu$.
|
||
|
||
# %%
|
||
def run_tgv(N=96, U0=0.04, nu=0.0125, steps=10000, rec=50, n_snap=10, frame_every=42):
|
||
k = 2 * np.pi / N
|
||
xs = np.arange(N)
|
||
X, Yc = np.meshgrid(xs, xs) # X — индекс x (ось 1), Yc — индекс y (ось 0)
|
||
ux = -U0 * np.cos(k * X) * np.sin(k * Yc)
|
||
uy = U0 * np.sin(k * X) * np.cos(k * Yc)
|
||
rho = np.ones((N, N))
|
||
f = feq_entropic(rho, np.stack([ux, uy]))
|
||
beta = nu_to_beta(nu)
|
||
# моменты времени для n_snap снимков поля скорости через равные интервалы
|
||
snap_at = sorted(set(int(round(v)) for v in np.linspace(0, steps, n_snap)))
|
||
snaps, frames, ts, KE, Hs = [], [], [], [], []
|
||
for t in range(steps + 1):
|
||
rho, u = macros(f)
|
||
f, _ = collide_kbc(f, feq_entropic(rho, u), beta)
|
||
f = stream(f)
|
||
need = (t in snap_at) or (t % rec == 0) or (t % frame_every == 0)
|
||
if need:
|
||
_, uu = macros(f)
|
||
if t in snap_at:
|
||
snaps.append((t, uu.copy())) # снимки для статической панели
|
||
if t % frame_every == 0:
|
||
frames.append((t, uu.copy())) # кадры для анимации
|
||
if t % rec == 0:
|
||
ts.append(t)
|
||
KE.append(0.5 * np.mean(uu[0] ** 2 + uu[1] ** 2))
|
||
Hs.append(H_function(f).mean())
|
||
return dict(N=N, k=k, nu=nu, U0=U0, ts=np.array(ts), KE=np.array(KE),
|
||
H=np.array(Hs), u=u, rho=rho, snaps=snaps, frames=frames)
|
||
|
||
|
||
tgv = run_tgv()
|
||
# Измеряем ν по наклону ln(KE): slope = -4 ν k^2
|
||
mask = tgv["KE"] > 0
|
||
slope = np.polyfit(tgv["ts"][mask], np.log(tgv["KE"][mask]), 1)[0]
|
||
nu_meas = -slope / (4 * tgv["k"] ** 2)
|
||
rel_err = abs(nu_meas - tgv["nu"]) / tgv["nu"]
|
||
print(f"TGV: nu_zadano={tgv['nu']:.5f}, nu_izmereno={nu_meas:.5f}, otn.oshibka={rel_err*100:.2f}%")
|
||
assert rel_err < 0.06, "TGV: измеренная вязкость отклонилась слишком сильно"
|
||
print("Валидация TGV (вязкость по затуханию энергии) — OK")
|
||
|
||
|
||
# %%
|
||
def vorticity(u):
|
||
duy_dx = (np.roll(u[1], -1, 1) - np.roll(u[1], 1, 1)) * 0.5
|
||
dux_dy = (np.roll(u[0], -1, 0) - np.roll(u[0], 1, 0)) * 0.5
|
||
return duy_dx - dux_dy
|
||
|
||
|
||
fig, axs = plt.subplots(1, 3, figsize=(13.5, 4.0))
|
||
om = vorticity(tgv["u"])
|
||
im = axs[0].imshow(om, origin="lower", cmap="RdBu_r")
|
||
axs[0].set_title("Завихренность $\\omega$ (TGV, срез)")
|
||
axs[0].set_xlabel("$x$ [lu]"); axs[0].set_ylabel("$y$ [lu]")
|
||
plt.colorbar(im, ax=axs[0], fraction=0.046)
|
||
|
||
axs[1].semilogy(tgv["ts"], tgv["KE"], "o", ms=3, label="KBC (измерено)")
|
||
axs[1].semilogy(tgv["ts"], tgv["KE"][0] * np.exp(-4 * tgv["nu"] * tgv["k"] ** 2 * tgv["ts"]),
|
||
"-", color="#d1495b", label="$E_0 e^{-4\\nu k^2 t}$ (аналитика)")
|
||
axs[1].set_title("Затухание кин. энергии"); axs[1].set_xlabel("$t$ [ts]")
|
||
axs[1].set_ylabel("$E(t)$"); axs[1].legend()
|
||
|
||
axs[2].plot(tgv["ts"], tgv["H"], color="#2e7d32")
|
||
axs[2].set_title("H-функция (энтропия) во времени")
|
||
axs[2].set_xlabel("$t$ [ts]"); axs[2].set_ylabel("$\\langle H\\rangle$")
|
||
fig.suptitle("Демо-1: вихрь Тейлора–Грина (KBC, D2Q9)", y=1.04, fontweight="bold")
|
||
_finish(fig, "07_taylor_green.png")
|
||
|
||
|
||
# %% [markdown]
|
||
# ### Распад вихря: 10 срезов через равные интервалы времени
|
||
#
|
||
# Поле завихренности $\omega=\partial_x u_y-\partial_y u_x$ в 10 моментов времени,
|
||
# равномерно распределённых от $t=0$ до конца прогона. **Шкала цвета единая для всех
|
||
# панелей** (фиксирована по начальной завихренности $\pm|\omega|_{\max}(0)$) — поэтому
|
||
# наглядно видно именно *затухание*: пространственная структура (4 противоположно
|
||
# вращающихся вихря) сохраняется, а амплитуда экспоненциально гаснет. В заголовке
|
||
# каждой панели — шаг $t$ и доля кинетической энергии $E/E_0$.
|
||
|
||
# %%
|
||
snaps = tgv["snaps"]
|
||
om0 = np.abs(vorticity(snaps[0][1])).max() # масштаб цвета по t=0
|
||
KE0 = 0.5 * np.mean(snaps[0][1][0] ** 2 + snaps[0][1][1] ** 2) # начальная энергия
|
||
|
||
fig, axs = plt.subplots(2, 5, figsize=(15, 6.3), constrained_layout=True)
|
||
im = None
|
||
for ax, (tt, uu) in zip(axs.ravel(), snaps):
|
||
im = ax.imshow(vorticity(uu), origin="lower", cmap="RdBu_r", vmin=-om0, vmax=om0)
|
||
KEf = 0.5 * np.mean(uu[0] ** 2 + uu[1] ** 2) / KE0
|
||
ax.set_title(f"t={tt} · $E/E_0$={KEf:.2f}", fontsize=10)
|
||
ax.set_xticks([]); ax.set_yticks([])
|
||
fig.colorbar(im, ax=axs, fraction=0.025, pad=0.01, label="завихренность $\\omega$ [1/ts]")
|
||
fig.suptitle("Распад вихря Тейлора–Грина: завихренность на равных интервалах времени "
|
||
"(единая шкала цвета)", fontweight="bold")
|
||
_finish(fig, "07b_tgv_decay_slices.png")
|
||
|
||
|
||
# %% [markdown]
|
||
# ## 8. Демо-2: устойчивость — двойной сдвиговый слой (BGK против KBC)
|
||
#
|
||
# Классический тест Minion–Brown: тонкие сдвиговые слои на **недоразрешённой**
|
||
# сетке. Стандартный BGK при малой вязкости разрушается (отрицательные
|
||
# популяции, NaN), а KBC за счёт адаптивного $\gamma$ остаётся устойчивым и даёт
|
||
# гладкие свёрнутые вихри — это и есть «implicit LES».
|
||
|
||
# %%
|
||
def init_shear(N, U0=0.04, delta=0.05, kappa=80.0):
|
||
xs = np.arange(N) / N
|
||
X, Yc = np.meshgrid(xs, xs)
|
||
ux = np.where(Yc <= 0.5, U0 * np.tanh(kappa * (Yc - 0.25)),
|
||
U0 * np.tanh(kappa * (0.75 - Yc)))
|
||
uy = delta * U0 * np.sin(2 * np.pi * (X + 0.25))
|
||
return np.stack([ux, uy])
|
||
|
||
|
||
def run_shear(scheme, N=128, U0=0.04, nu=1.7e-4, steps=8000, rec=100, frame_every=33):
|
||
u0 = init_shear(N, U0=U0)
|
||
rho = np.ones((N, N))
|
||
beta = nu_to_beta(nu)
|
||
feq0 = feq_entropic if scheme == "kbc" else feq_poly
|
||
f = feq0(rho, u0)
|
||
ts, umax, frames = [], [], []
|
||
blew_up_at = None
|
||
last_u, last_g = u0, np.full((N, N), 2.0)
|
||
for t in range(steps + 1):
|
||
rho, u = macros(f)
|
||
if not np.all(np.isfinite(u)) or np.abs(u).max() > 0.5:
|
||
blew_up_at = t
|
||
break
|
||
if t % frame_every == 0:
|
||
frames.append((t, u.copy())) # кадры (только до разрушения)
|
||
if scheme == "kbc":
|
||
f, g = collide_kbc(f, feq_entropic(rho, u), beta)
|
||
last_g = g
|
||
else:
|
||
f = collide_bgk(f, feq_poly(rho, u), beta)
|
||
f = stream(f)
|
||
last_u = u
|
||
if t % rec == 0:
|
||
ts.append(t); umax.append(np.abs(u).max())
|
||
return dict(scheme=scheme, ts=np.array(ts), umax=np.array(umax), frames=frames,
|
||
u=last_u, gamma=last_g, blew_up_at=blew_up_at, N=N, U0=U0, nu=nu)
|
||
|
||
|
||
Re_shear = 0.04 * 128 / 1.7e-4
|
||
print(f"Sdvigovyj sloj: Re ~ {Re_shear:.0f}, beta ~ {nu_to_beta(1.7e-4):.4f}")
|
||
res_bgk = run_shear("bgk")
|
||
res_kbc = run_shear("kbc")
|
||
print("BGK:", ("raz=" + str(res_bgk["blew_up_at"])) if res_bgk["blew_up_at"] else "ustojchiv")
|
||
print("KBC:", ("raz=" + str(res_kbc["blew_up_at"])) if res_kbc["blew_up_at"] else "ustojchiv")
|
||
assert res_kbc["blew_up_at"] is None, "KBC должен оставаться устойчивым"
|
||
print("Демонстрация устойчивости BGK vs KBC — OK")
|
||
|
||
|
||
# %%
|
||
fig, axs = plt.subplots(1, 3, figsize=(13.5, 4.0))
|
||
om_k = vorticity(res_kbc["u"])
|
||
im0 = axs[0].imshow(om_k, origin="lower", cmap="RdBu_r")
|
||
axs[0].set_title("KBC: завихренность (устойчиво)")
|
||
axs[0].set_xlabel("$x$ [lu]"); axs[0].set_ylabel("$y$ [lu]")
|
||
plt.colorbar(im0, ax=axs[0], fraction=0.046)
|
||
|
||
im1 = axs[1].imshow(res_kbc["gamma"] - 2.0, origin="lower", cmap="magma")
|
||
axs[1].set_title("Поле стабилизатора $\\gamma-2$")
|
||
axs[1].set_xlabel("$x$ [lu]"); axs[1].set_ylabel("$y$ [lu]")
|
||
plt.colorbar(im1, ax=axs[1], fraction=0.046)
|
||
|
||
axs[2].plot(res_kbc["ts"], res_kbc["umax"], color="#2e7d32", label="KBC")
|
||
axs[2].plot(res_bgk["ts"], res_bgk["umax"], color="#d1495b", label="BGK")
|
||
if res_bgk["blew_up_at"]:
|
||
axs[2].axvline(res_bgk["blew_up_at"], color="#d1495b", ls="--", alpha=0.6)
|
||
axs[2].set_title("$\\max|\\mathbf{u}|$ во времени")
|
||
axs[2].set_xlabel("$t$ [ts]"); axs[2].set_ylabel("$\\max|\\mathbf{u}|$ [lu/ts]"); axs[2].legend()
|
||
fig.suptitle(f"Демо-2: сдвиговый слой, Re~{Re_shear:.0f} (BGK разрушается, KBC держит)",
|
||
y=1.04, fontweight="bold")
|
||
_finish(fig, "08_shear_stability.png")
|
||
|
||
|
||
# %% [markdown]
|
||
# ## 9. Демо-3: каверна с движущейся крышкой против эталона Ghia (Re=1000)
|
||
#
|
||
# Квадратная каверна, верхняя стенка движется со скоростью $U$. На неподвижных
|
||
# стенках — **bounce-back** (отскок), на крышке — bounce-back с поправкой импульса
|
||
# подвижной стенки $-2 w_i \rho\,(\mathbf c_i\cdot\mathbf u_w)/c_s^2$. Сравниваем
|
||
# профили скорости по центральным линиям с данными **Ghia, Ghia & Shin (1982)**.
|
||
# Здесь сетка намеренно грубая (для быстрого прогона), а bounce-back — первого
|
||
# порядка, поэтому совпадение **качественное**; оно улучшается измельчением сетки,
|
||
# числом шагов и переходом к ГУ второго порядка (Zou–He / Grad на крышке).
|
||
|
||
# %%
|
||
UP = [2, 5, 6] # c_iy = +1 (вверх)
|
||
DOWN = [4, 7, 8] # c_iy = -1 (вниз)
|
||
RIGHT = [1, 5, 8] # c_ix = +1
|
||
LEFT = [3, 6, 7] # c_ix = -1
|
||
|
||
|
||
def run_cavity(N=110, Re=1000.0, U=0.05, steps=14000, tol=1e-7, check=500, frame_every=58):
|
||
nu = U * N / Re
|
||
beta = nu_to_beta(nu)
|
||
rho = np.ones((N, N))
|
||
u = np.zeros((2, N, N))
|
||
f = feq_entropic(rho, u)
|
||
uw = np.array([U, 0.0])
|
||
prev = None
|
||
frames = []
|
||
t = 0
|
||
for t in range(steps):
|
||
rho, u = macros(f)
|
||
if t % frame_every == 0:
|
||
frames.append((t, u.copy())) # кадры становления течения
|
||
fcol, _ = collide_kbc(f, feq_entropic(rho, u), beta)
|
||
f = stream(fcol)
|
||
# bounce-back: на каждой стенке неизвестны направления, входящие в область.
|
||
for i in UP: # низ (y=0)
|
||
f[i, 0, :] = fcol[OPP[i], 0, :]
|
||
for i in RIGHT: # лево (x=0)
|
||
f[i, :, 0] = fcol[OPP[i], :, 0]
|
||
for i in LEFT: # право (x=N-1)
|
||
f[i, :, -1] = fcol[OPP[i], :, -1]
|
||
for i in DOWN: # верх (крышка) + поправка подвижной стенки
|
||
io = OPP[i] # io ударяет в крышку (идёт вверх)
|
||
corr = 2.0 * W[io] * (CX[io] * uw[0] + CY[io] * uw[1]) / CS2
|
||
f[i, -1, :] = fcol[io, -1, :] - corr
|
||
if t % check == 0 and t > 0:
|
||
cur = u.copy()
|
||
if prev is not None and np.abs(cur - prev).max() < tol:
|
||
break
|
||
prev = cur
|
||
rho, u = macros(f)
|
||
return dict(N=N, Re=Re, U=U, u=u, steps=t, frames=frames)
|
||
|
||
|
||
cav = run_cavity()
|
||
print(f"Kaverna: Re={cav['Re']:.0f}, N={cav['N']}, ostanovleno na shage {cav['steps']}")
|
||
|
||
# Эталон Ghia, Ghia & Shin (1982), Re=1000.
|
||
ghia_y = np.array([0.0000, 0.0547, 0.0625, 0.0703, 0.1016, 0.1719, 0.2813, 0.4531,
|
||
0.5000, 0.6172, 0.7344, 0.8516, 0.9531, 0.9609, 0.9688, 0.9766, 1.0000])
|
||
ghia_u = np.array([0.0000, -0.18109, -0.20196, -0.22220, -0.29730, -0.38289, -0.27805,
|
||
-0.10648, -0.06080, 0.05702, 0.18719, 0.33304, 0.46604, 0.51117,
|
||
0.57492, 0.65928, 1.0000])
|
||
ghia_x = np.array([0.0000, 0.0625, 0.0703, 0.0781, 0.0938, 0.1563, 0.2266, 0.2344,
|
||
0.5000, 0.8047, 0.8594, 0.9063, 0.9453, 0.9531, 0.9609, 0.9688, 1.0000])
|
||
ghia_v = np.array([0.0000, 0.27485, 0.29012, 0.30353, 0.32627, 0.37095, 0.33075, 0.32235,
|
||
0.02526, -0.31966, -0.42665, -0.51550, -0.39188, -0.33714, -0.27669,
|
||
-0.21388, 0.0000])
|
||
|
||
Ncav = cav["N"]
|
||
yc = (np.arange(Ncav) + 0.5) / Ncav
|
||
xc = (np.arange(Ncav) + 0.5) / Ncav
|
||
u_centerline = cav["u"][0, :, Ncav // 2] / cav["U"] # u_x(y) на вертикали x=0.5
|
||
v_centerline = cav["u"][1, Ncav // 2, :] / cav["U"] # u_y(x) на горизонтали y=0.5
|
||
|
||
fig, axs = plt.subplots(1, 3, figsize=(14, 4.4))
|
||
spd = np.sqrt(cav["u"][0] ** 2 + cav["u"][1] ** 2) / cav["U"]
|
||
axs[0].imshow(spd, origin="lower", cmap="viridis", extent=[0, 1, 0, 1])
|
||
sx = np.linspace(0, 1, Ncav)
|
||
axs[0].streamplot(sx, sx, cav["u"][0], cav["u"][1], density=1.1, color="white", linewidth=0.6)
|
||
axs[0].set_title("Линии тока и $|\\mathbf{u}|/U$"); axs[0].set_xlabel("$x/L$"); axs[0].set_ylabel("$y/L$")
|
||
|
||
axs[1].plot(u_centerline, yc, "-", color="#1f6feb", label="KBC")
|
||
axs[1].plot(ghia_u, ghia_y, "o", ms=5, mfc="none", color="#d1495b", label="Ghia 1982")
|
||
axs[1].set_title("$u_x/U$ по вертикали $x=0.5$"); axs[1].set_xlabel("$u_x/U$")
|
||
axs[1].set_ylabel("$y/L$"); axs[1].legend()
|
||
|
||
axs[2].plot(xc, v_centerline, "-", color="#1f6feb", label="KBC")
|
||
axs[2].plot(ghia_x, ghia_v, "o", ms=5, mfc="none", color="#d1495b", label="Ghia 1982")
|
||
axs[2].set_title("$u_y/U$ по горизонтали $y=0.5$"); axs[2].set_xlabel("$x/L$")
|
||
axs[2].set_ylabel("$u_y/U$"); axs[2].legend()
|
||
fig.suptitle("Демо-3: каверна с движущейся крышкой, Re=1000 (KBC vs Ghia)",
|
||
y=1.04, fontweight="bold")
|
||
_finish(fig, "09_cavity_ghia.png")
|
||
|
||
u_interp = np.interp(ghia_y, yc, u_centerline)
|
||
rms = np.sqrt(np.mean((u_interp - ghia_u) ** 2))
|
||
print(f"Kaverna: RMS-otklonenie u_x ot Ghia ~ {rms:.3f} (gruboe BB -> dopustimo)")
|
||
|
||
|
||
# %% [markdown]
|
||
# ## 9b. Анимации: зацикленное воспроизведение тестов (видео)
|
||
#
|
||
# Для каждого теста собираем **зацикленный GIF** (бесконечный `loop`) из кадров,
|
||
# записанных во время прогона: вихрь Тейлора–Грина (затухание), сдвиговый слой
|
||
# (BGK против KBC бок о бок — видно момент разрушения BGK) и становление течения в
|
||
# каверне. Файлы сохраняются в `anim/` (самостоятельные «видео»), а ниже встроены
|
||
# в ноутбук и проигрываются циклически. Скорость воспроизведения — `fps` кадров/с
|
||
# (физическое время между кадрами постоянно).
|
||
|
||
# %%
|
||
import io as _io
|
||
from PIL import Image
|
||
from IPython.display import HTML, display
|
||
|
||
ANIMDIR = os.path.join(os.path.dirname(os.path.abspath(__file__)) if "__file__" in globals() else ".", "anim")
|
||
os.makedirs(ANIMDIR, exist_ok=True)
|
||
|
||
|
||
def render_gif(n_frames, draw_fn, name, figsize, fps=18, dpi=80):
|
||
"""Отрисовать n_frames кадров (draw_fn(fig, i)) и сохранить зацикленный GIF."""
|
||
path = os.path.join(ANIMDIR, name)
|
||
fig = plt.figure(figsize=figsize, dpi=dpi)
|
||
pil = []
|
||
for i in range(n_frames):
|
||
fig.clf()
|
||
draw_fn(fig, i)
|
||
buf = _io.BytesIO()
|
||
fig.savefig(buf, format="png") # фикс. размер кадра (figsize×dpi)
|
||
buf.seek(0)
|
||
pil.append(Image.open(buf).convert("RGB").convert("P", palette=Image.ADAPTIVE))
|
||
plt.close(fig)
|
||
pil[0].save(path, save_all=True, append_images=pil[1:],
|
||
duration=int(1000 / fps), loop=0, optimize=True) # loop=0 → бесконечно
|
||
print("GIF:", os.path.relpath(path), f"({n_frames} кадров)")
|
||
return path
|
||
|
||
|
||
def show_gif(path, width=460):
|
||
"""Встроить GIF в ноутбук (ссылкой на файл — авто-цикл в браузере/Jupyter)."""
|
||
if IN_NOTEBOOK:
|
||
display(HTML(f'<img src="anim/{os.path.basename(path)}" width="{width}"/>'))
|
||
|
||
|
||
# %% [markdown]
|
||
# ### Видео 1 — затухание вихря Тейлора–Грина
|
||
|
||
# %%
|
||
_tgv_fr = tgv["frames"]
|
||
_om0 = np.abs(vorticity(_tgv_fr[0][1])).max()
|
||
_KE0 = 0.5 * np.mean(_tgv_fr[0][1][0] ** 2 + _tgv_fr[0][1][1] ** 2)
|
||
|
||
def _draw_tgv(fig, i):
|
||
t, u = _tgv_fr[i]
|
||
ax = fig.add_subplot(111)
|
||
ax.imshow(vorticity(u), origin="lower", cmap="RdBu_r", vmin=-_om0, vmax=_om0)
|
||
ke = 0.5 * np.mean(u[0] ** 2 + u[1] ** 2) / _KE0
|
||
ax.set_title(f"Тейлор–Грин · t={t} · $E/E_0$={ke:.2f}", fontsize=10)
|
||
ax.set_xticks([]); ax.set_yticks([])
|
||
fig.tight_layout()
|
||
|
||
show_gif(render_gif(len(_tgv_fr), _draw_tgv, "tgv_decay.gif", (4.6, 4.4), fps=24))
|
||
|
||
# %% [markdown]
|
||
# ### Видео 2 — сдвиговый слой: BGK (слева) против KBC (справа)
|
||
# BGK замирает на последнем устойчивом кадре с пометкой «разрушился», KBC
|
||
# продолжает гладко сворачивать вихри.
|
||
|
||
# %%
|
||
_kf, _bf = res_kbc["frames"], res_bgk["frames"]
|
||
_om0s = np.percentile(np.abs(vorticity(_kf[len(_kf) // 2][1])), 99.5)
|
||
|
||
def _draw_shear(fig, i):
|
||
ax1 = fig.add_subplot(1, 2, 1)
|
||
ax2 = fig.add_subplot(1, 2, 2)
|
||
dead = i >= len(_bf)
|
||
tb, ub = _bf[min(i, len(_bf) - 1)]
|
||
ax1.imshow(vorticity(ub), origin="lower", cmap="RdBu_r", vmin=-_om0s, vmax=_om0s)
|
||
ax1.set_title("BGK — разрушился" if dead else f"BGK · t={tb}",
|
||
color="red" if dead else "black", fontsize=10)
|
||
tk, uk = _kf[i]
|
||
ax2.imshow(vorticity(uk), origin="lower", cmap="RdBu_r", vmin=-_om0s, vmax=_om0s)
|
||
ax2.set_title(f"KBC · t={tk}", fontsize=10)
|
||
for a in (ax1, ax2):
|
||
a.set_xticks([]); a.set_yticks([])
|
||
fig.tight_layout()
|
||
|
||
show_gif(render_gif(len(_kf), _draw_shear, "shear_bgk_vs_kbc.gif", (8.4, 4.3), fps=24), width=720)
|
||
|
||
# %% [markdown]
|
||
# ### Видео 3 — становление течения в каверне (Re=1000)
|
||
|
||
# %%
|
||
_cf = cav["frames"]
|
||
_vmax = max(np.sqrt(u[0] ** 2 + u[1] ** 2).max() for _, u in _cf) / cav["U"]
|
||
|
||
def _draw_cav(fig, i):
|
||
t, u = _cf[i]
|
||
ax = fig.add_subplot(111)
|
||
sp = np.sqrt(u[0] ** 2 + u[1] ** 2) / cav["U"]
|
||
ax.imshow(sp, origin="lower", cmap="viridis", vmin=0, vmax=_vmax, extent=[0, 1, 0, 1])
|
||
ax.set_title(f"Каверна Re=1000 · t={t} · $|\\mathbf{{u}}|/U$", fontsize=10)
|
||
ax.set_xlabel("$x/L$"); ax.set_ylabel("$y/L$")
|
||
fig.tight_layout()
|
||
|
||
show_gif(render_gif(len(_cf), _draw_cav, "cavity_transient.gif", (4.9, 4.5), fps=24))
|
||
|
||
|
||
# %% [markdown]
|
||
# ## Демо-4: обтекание твёрдых тел разной формы (2D)
|
||
#
|
||
# Канал с втоком слева (Дирихле + плавный разгон), нуль-градиентным оттоком справа и
|
||
# стенками сверху/снизу. Тело задаётся маской и отражает популяции (**bounce-back**).
|
||
# За плохо обтекаемыми телами при Re≈150 формируется **вихревая дорожка Кармана** —
|
||
# поочерёдный срыв вихрей; KBC держит течение устойчивым на грубой сетке. Тела: цилиндр,
|
||
# квадрат, треугольник, профиль (под углом атаки). Ниже — анимация всех четырёх (≈10 c),
|
||
# статический срез завихренности и колебания скорости в следе (мерило частоты срыва).
|
||
|
||
# %%
|
||
def smoothstep(x):
|
||
x = min(max(x, 0.0), 1.0)
|
||
return x * x * (3.0 - 2.0 * x)
|
||
|
||
|
||
def solid_mask(Nx, Ny, shape, cx, cy, D):
|
||
Y, X = np.mgrid[0:Ny, 0:Nx]
|
||
if shape == "circle":
|
||
m = (X - cx) ** 2 + (Y - cy) ** 2 <= (D / 2) ** 2
|
||
elif shape == "square":
|
||
m = (np.abs(X - cx) <= D / 2) & (np.abs(Y - cy) <= D / 2)
|
||
elif shape == "triangle": # вершиной против потока
|
||
tt = (X - (cx - D / 2)) / D
|
||
m = (tt >= 0) & (tt <= 1) & (np.abs(Y - cy) <= tt * (D / 2))
|
||
elif shape == "airfoil": # наклонный эллипс (угол атаки 18°)
|
||
aoa = np.deg2rad(18.0); a = D * 0.95; b = D * 0.26
|
||
xr = (X - cx) * np.cos(aoa) + (Y - cy) * np.sin(aoa)
|
||
yr = -(X - cx) * np.sin(aoa) + (Y - cy) * np.cos(aoa)
|
||
m = (xr / a) ** 2 + (yr / b) ** 2 <= 1
|
||
else:
|
||
m = np.zeros((Ny, Nx), bool)
|
||
m[0, :] = True; m[-1, :] = True # стенки канала (верх/низ)
|
||
return m
|
||
|
||
|
||
def run_obstacle(shape, Nx=150, Ny=72, D=14, U=0.06, Re=150.0, steps=9600,
|
||
ramp=1000, frame_every=40):
|
||
cx, cy = Nx // 4, Ny // 2
|
||
solid = solid_mask(Nx, Ny, shape, cx, cy, D)
|
||
beta = nu_to_beta(U * D / Re)
|
||
f = feq_entropic(np.ones((Ny, Nx)), np.zeros((2, Ny, Nx)))
|
||
frames, probe = [], []
|
||
px, py = cx + 3 * D, cy # зонд в следе
|
||
ones_col = np.ones((Ny, 1))
|
||
for t in range(steps + 1):
|
||
rho, u = macros(f)
|
||
fcol, _ = collide_kbc(f, feq_entropic(rho, u), beta)
|
||
for i in range(Q):
|
||
fcol[i, solid] = f[OPP[i], solid] # bounce-back (тело + стенки)
|
||
f = stream(fcol)
|
||
Uin = U * smoothstep(t / ramp) # плавный разгон входа
|
||
uin = np.zeros((2, Ny, 1)); uin[0, :, 0] = Uin
|
||
f[:, 1:-1, 0] = feq_entropic(ones_col, uin)[:, 1:-1, 0] # вток (Дирихле)
|
||
f[:, 1:-1, -1] = f[:, 1:-1, -2] # отток (нуль-градиент)
|
||
if t % frame_every == 0:
|
||
_, uu = macros(f)
|
||
frames.append((t, uu.copy()))
|
||
probe.append(uu[1, py, px])
|
||
return dict(shape=shape, solid=solid, frames=frames, probe=np.array(probe),
|
||
Nx=Nx, Ny=Ny, U=U, Re=Re, D=D, cx=cx, cy=cy)
|
||
|
||
|
||
SHAPES = [("circle", "цилиндр"), ("square", "квадрат"),
|
||
("triangle", "треугольник"), ("airfoil", "профиль")]
|
||
obs = {key: run_obstacle(key) for key, _ in SHAPES}
|
||
for key, name in SHAPES:
|
||
p = obs[key]["probe"]; half = p[len(p) // 2:]
|
||
print(f"{name:12}: max|u_y(зонд)|={np.abs(half).max():.4f} (колебания → срыв вихрей)")
|
||
|
||
|
||
# %%
|
||
# Анимация-монтаж всех четырёх тел (завихренность, тело — тёмное).
|
||
_keys = [k for k, _ in SHAPES]
|
||
_titles = {k: n for k, n in SHAPES}
|
||
_nf = min(len(obs[k]["frames"]) for k in _keys)
|
||
_vmax = max(np.percentile(np.abs(vorticity(obs[k]["frames"][_nf // 2][1])), 99.0) for k in _keys)
|
||
_ocmap = plt.get_cmap("RdBu_r").copy(); _ocmap.set_bad("#1b1b1b")
|
||
|
||
def _draw_obst(fig, i):
|
||
for j, k in enumerate(_keys):
|
||
ax = fig.add_subplot(2, 2, j + 1)
|
||
t, u = obs[k]["frames"][i]
|
||
vort = np.where(obs[k]["solid"], np.nan, vorticity(u))
|
||
ax.imshow(vort, origin="lower", cmap=_ocmap, vmin=-_vmax, vmax=_vmax)
|
||
ax.set_title(f"{_titles[k]} · t={t}", fontsize=9)
|
||
ax.set_xticks([]); ax.set_yticks([])
|
||
fig.suptitle("Обтекание тел (KBC, Re≈150): завихренность", fontsize=11)
|
||
fig.tight_layout(rect=(0, 0, 1, 0.96))
|
||
|
||
show_gif(render_gif(_nf, _draw_obst, "obstacle_shapes.gif", (11, 6.2), fps=24), width=780)
|
||
|
||
# %%
|
||
# Статический срез завихренности в конце прогона + колебания в следе.
|
||
fig, axs = plt.subplots(2, 2, figsize=(12, 5.6), constrained_layout=True)
|
||
for ax, k in zip(axs.ravel(), _keys):
|
||
t, u = obs[k]["frames"][-1]
|
||
vort = np.where(obs[k]["solid"], np.nan, vorticity(u))
|
||
ax.imshow(vort, origin="lower", cmap=_ocmap, vmin=-_vmax, vmax=_vmax)
|
||
ax.set_title(f"{_titles[k]} (t={t})", fontsize=10)
|
||
ax.set_xticks([]); ax.set_yticks([])
|
||
fig.suptitle("Обтекание тел: поле завихренности в конце прогона", fontweight="bold")
|
||
_finish(fig, "10_obstacle_shapes.png")
|
||
|
||
fig, ax = plt.subplots(figsize=(9.5, 3.4))
|
||
for k in _keys:
|
||
ax.plot(np.arange(len(obs[k]["probe"])), obs[k]["probe"], label=_titles[k])
|
||
ax.set_title("Поперечная скорость в следе (зонд за телом): колебания = срыв вихрей")
|
||
ax.set_xlabel("кадр"); ax.set_ylabel("$u_y$ в зонде [lu/ts]"); ax.legend(fontsize=8, ncol=4)
|
||
_finish(fig, "10b_obstacle_probe.png")
|
||
|
||
|
||
# %% [markdown]
|
||
# ## KBC на сетках с измельчением у границ (grid refinement)
|
||
#
|
||
# **Может ли KBC работать с неоднородным разрешением?** Да — но важно, *как именно*:
|
||
#
|
||
# 1. **Чистый LBM живёт на РАВНОМЕРНОЙ решётке.** Перенос $f_i(\mathbf x+\mathbf c_i)$
|
||
# требует, чтобы сосед был ровно на расстоянии $\mathbf c_i$ — шаг сетки постоянен.
|
||
# Поэтому одна «кривая» неравномерная сетка для классического LBM не подходит.
|
||
# 2. **Неоднородное разрешение делают блочным измельчением (multi-block / AMR):** область
|
||
# делят на вложенные РАВНОМЕРНЫЕ блоки разного шага — крупный вдали, мелкий у тел и в
|
||
# следе. На стыке блоков поля **интерполируют** и **перемасштабируют**: при переходе на
|
||
# вдвое более мелкий блок $\Delta x\!\to\!\Delta x/2$, $\Delta t\!\to\!\Delta t/2$
|
||
# (акустическое масштабирование, $c_s$ сохраняется), время релаксации меняется как
|
||
# $\tau_f = 2\tau_c - \tfrac12$ (вязкость $\nu$ постоянна), а неравновесная часть
|
||
# $f^{neq}$ масштабируется множителем отношения шагов. ⚠ детали сверить с оригиналом.
|
||
# 3. **KBC полностью совместим с измельчением.** Это локальный оператор столкновения: на
|
||
# каждом блоке своё $\beta$, $\gamma$ — поузловой как обычно. Энтропийная устойчивость
|
||
# особенно помогает на стыках крупный↔мелкий блок и в пограничном слое у стенок, где
|
||
# градиенты резкие. KBC применяли с измельчением и на телесно-подогнанных сетках
|
||
# (Dorschner, Chikatamarla, Karlin, *JFM* 2018; grid refinement — arXiv:1608.06915).
|
||
# 4. **Альтернатива у кривых стенок** — интерполированный bounce-back (Bouzidi) на
|
||
# равномерной решётке (доля стенки $q$). Измельчение добавляют там, где надо разрешить
|
||
# тонкий пограничный слой/след.
|
||
#
|
||
# Полноценный AMR-решатель — отдельная большая работа. Ниже **схема** многоуровневого
|
||
# измельчения вокруг тела, наложенная на реальное поле $|\omega|$ из расчёта цилиндра —
|
||
# видно, ГДЕ нужна мелкая сетка (резкие градиенты у тела и в следе).
|
||
|
||
# %%
|
||
import matplotlib.patches as mpatches
|
||
|
||
rc = obs["circle"]
|
||
_t, _u = rc["frames"][-1]
|
||
vmag = np.where(rc["solid"], np.nan, np.abs(vorticity(_u)))
|
||
fig, ax = plt.subplots(figsize=(11, 5.6))
|
||
ax.imshow(vmag, origin="lower", cmap="inferno",
|
||
vmax=np.nanpercentile(vmag, 99))
|
||
cx, cy, D = rc["cx"], rc["cy"], rc["D"]
|
||
|
||
def _grid_block(x0, y0, w, h, step, color, label):
|
||
ax.add_patch(mpatches.Rectangle((x0, y0), w, h, fill=False, edgecolor=color, lw=2.2))
|
||
for xx in np.arange(x0, x0 + w + 0.1, step):
|
||
ax.plot([xx, xx], [y0, y0 + h], color=color, lw=0.4, alpha=0.45)
|
||
for yy in np.arange(y0, y0 + h + 0.1, step):
|
||
ax.plot([x0, x0 + w], [yy, yy], color=color, lw=0.4, alpha=0.45)
|
||
ax.text(x0 + 1, y0 + h - 4, label, color=color, fontsize=8, weight="bold")
|
||
|
||
_grid_block(0, 0, rc["Nx"] - 1, rc["Ny"] - 1, 16, "#4fc3f7", "уровень 0: Δx (крупный)")
|
||
_grid_block(cx - 1.5 * D, cy - 1.9 * D, 6.5 * D, 3.8 * D, 8, "#aed581", "уровень 1: Δx/2 (тело+след)")
|
||
_grid_block(cx - 1.1 * D, cy - 1.3 * D, 2.6 * D, 2.6 * D, 4, "#ff8a65", "уровень 2: Δx/4 (у тела)")
|
||
ax.set_xlim(0, rc["Nx"] - 1); ax.set_ylim(0, rc["Ny"] - 1)
|
||
ax.set_title("Схема блочного измельчения вокруг тела (поверх $|\\omega|$): "
|
||
"сетка мельче там, где градиенты резче")
|
||
ax.set_xlabel("x [lu]"); ax.set_ylabel("y [lu]")
|
||
_finish(fig, "11c_refinement_schematic.png")
|
||
|
||
|
||
# %% [markdown]
|
||
# ## Реализация AMR-солвера (2 уровня) и демонстрация на цилиндре
|
||
#
|
||
# По исследованию выше реализуем **рабочий 2-уровневый блочный refinement**:
|
||
# - крупная сетка покрывает весь канал (с цилиндром и стенками);
|
||
# - **мелкий блок 2×** охватывает цилиндр (то же тело bounce-back, но вдвое мельче);
|
||
# - связка: $\tau_f = 2\tau_c-\tfrac12$ (вязкость сохраняется); на каждый крупный шаг —
|
||
# **2 мелких подшага**; граница блока заполняется из крупной сетки (билинейная
|
||
# интерполяция $\rho,\mathbf u$ + равновесие + перемасштаб $f^{neq}$ на $\tau_f/(2\tau_c)$),
|
||
# а внутренние крупные узлы обновляются из мелкой (обратный масштаб); рестрикция —
|
||
# только по жидким узлам.
|
||
#
|
||
# **Проверка связки** на отдельном гладком тесте (Тейлора–Грина с внутренним мелким
|
||
# блоком): схема устойчива, отн. отклонение AMR-крупной сетки от равномерной крупной
|
||
# ≈ **0.7 %**, «шов» на границе блока ≈ **1 %** — мелкий блок не вносит артефактов, KBC-γ
|
||
# работает на обоих уровнях. ⚠ Это **минимальный AMR** (линейная пространственная
|
||
# интерполяция, без временно́й интерполяции границы между подшагами); детали продакшн-схемы
|
||
# сверить с оригиналом (Lagrava et al. 2012; Dorschner et al., *JFM* 2018).
|
||
#
|
||
# Ниже — анимация: крупное поле завихренности с **наложенным мелким блоком 2×** вокруг
|
||
# цилиндра (видно более тонкое разрешение у тела).
|
||
|
||
# %%
|
||
def run_amr_cylinder(Nx=120, Ny=60, D=10, U=0.06, Re=120.0, steps=5400,
|
||
ramp=1000, frame_every=22):
|
||
cx, cy = Nx // 4, Ny // 2
|
||
nu = U * D / Re
|
||
tau_c = nu / CS2 + 0.5; beta_c = nu_to_beta(nu)
|
||
tau_f = 2 * tau_c - 0.5; beta_f = 1.0 / (2 * tau_f)
|
||
Rcf = tau_f / (2 * tau_c); Rfc = 1.0 / Rcf
|
||
# мелкий блок (в крупных клетках), не касается стенок канала
|
||
bx0 = max(cx - 2 * D, 2); bx1 = min(cx + 3 * D, Nx - 2)
|
||
by0 = max(cy - 2 * D, 2); by1 = min(cy + 2 * D, Ny - 2)
|
||
Wx, Wy = bx1 - bx0, by1 - by0
|
||
Nfx, Nfy = 2 * Wx + 1, 2 * Wy + 1
|
||
# маски тела
|
||
Yc_, Xc_ = np.mgrid[0:Ny, 0:Nx]
|
||
solid_c = (Xc_ - cx) ** 2 + (Yc_ - cy) ** 2 <= (D / 2) ** 2
|
||
solid_c[0, :] = True; solid_c[-1, :] = True
|
||
fcx, fcy = 2 * (cx - bx0), 2 * (cy - by0)
|
||
Yf_, Xf_ = np.mgrid[0:Nfy, 0:Nfx]
|
||
solid_f = (Xf_ - fcx) ** 2 + (Yf_ - fcy) ** 2 <= float(D) ** 2
|
||
# билинейная интерполяция крупная → мелкая граница (индексы фиксированы)
|
||
FX, FY = np.meshgrid(bx0 + 0.5 * np.arange(Nfx), by0 + 0.5 * np.arange(Nfy))
|
||
ix0 = np.floor(FX).astype(int); iy0 = np.floor(FY).astype(int)
|
||
ix1 = np.minimum(ix0 + 1, Nx - 1); iy1 = np.minimum(iy0 + 1, Ny - 1)
|
||
wx = FX - ix0; wy = FY - iy0
|
||
|
||
def interp(fld):
|
||
return (fld[..., iy0, ix0] * (1 - wx) * (1 - wy) + fld[..., iy0, ix1] * wx * (1 - wy)
|
||
+ fld[..., iy1, ix0] * (1 - wx) * wy + fld[..., iy1, ix1] * wx * wy)
|
||
|
||
slx, sly = slice(2, 2 * Wx, 2), slice(2, 2 * Wy, 2)
|
||
block_fluid = ~solid_c[by0 + 1:by1, bx0 + 1:bx1]
|
||
onecol = np.ones((Ny, 1))
|
||
fc = feq_entropic(np.ones((Ny, Nx)), np.zeros((2, Ny, Nx)))
|
||
ff = feq_entropic(np.ones((Nfy, Nfx)), np.zeros((2, Nfy, Nfx)))
|
||
frames = []
|
||
for t in range(steps + 1):
|
||
rho, u = macros(fc)
|
||
pre = fc.copy()
|
||
fc = collide_kbc(fc, feq_entropic(rho, u), beta_c)[0]
|
||
for i in range(Q):
|
||
fc[i, solid_c] = pre[OPP[i], solid_c]
|
||
fc = stream(fc)
|
||
uin = np.zeros((2, Ny, 1)); uin[0, :, 0] = U * smoothstep(t / ramp)
|
||
fc[:, 1:-1, 0] = feq_entropic(onecol, uin)[:, 1:-1, 0]
|
||
fc[:, 1:-1, -1] = fc[:, 1:-1, -2]
|
||
# крупная → мелкая (ghost)
|
||
rcd, ucd = macros(fc); neqc = fc - feq_entropic(rcd, ucd)
|
||
gh = (feq_entropic(interp(rcd), np.stack([interp(ucd[0]), interp(ucd[1])]))
|
||
+ Rcf * interp(neqc))
|
||
for _ in range(2):
|
||
rf, uf = macros(ff); pref = ff.copy()
|
||
ff = collide_kbc(ff, feq_entropic(rf, uf), beta_f)[0]
|
||
for i in range(Q):
|
||
ff[i, solid_f] = pref[OPP[i], solid_f]
|
||
ff = stream(ff)
|
||
ff[:, 0, :] = gh[:, 0, :]; ff[:, -1, :] = gh[:, -1, :]
|
||
ff[:, :, 0] = gh[:, :, 0]; ff[:, :, -1] = gh[:, :, -1]
|
||
# мелкая → крупная (рестрикция, только жидкие узлы)
|
||
rf, uf = macros(ff); neqf = ff - feq_entropic(rf, uf)
|
||
newv = feq_entropic(rf[sly, slx], uf[:, sly, slx]) + Rfc * neqf[:, sly, slx]
|
||
cur = fc[:, by0 + 1:by1, bx0 + 1:bx1]
|
||
fc[:, by0 + 1:by1, bx0 + 1:bx1] = np.where(block_fluid[None], newv, cur)
|
||
if t % frame_every == 0:
|
||
_, uuc = macros(fc); _, uuf = macros(ff)
|
||
frames.append((t, uuc.copy(), uuf.copy()))
|
||
return dict(frames=frames, solid_c=solid_c, solid_f=solid_f, Nx=Nx, Ny=Ny,
|
||
bx0=bx0, bx1=bx1, by0=by0, by1=by1, tau_c=tau_c, tau_f=tau_f, Rcf=Rcf)
|
||
|
||
|
||
amr = run_amr_cylinder()
|
||
print(f"AMR-цилиндр: крупная {amr['Nx']}x{amr['Ny']} + мелкий блок 2x, "
|
||
f"tau_c={amr['tau_c']:.3f} tau_f={amr['tau_f']:.3f} Rcf={amr['Rcf']:.3f}, "
|
||
f"кадров={len(amr['frames'])} (устойчиво)")
|
||
|
||
# %%
|
||
# Анимация AMR: крупное поле + наложенный мелкий блок (видно тонкое разрешение у тела).
|
||
_af = amr["frames"]
|
||
_vmaxA = np.nanpercentile(np.abs(np.where(amr["solid_c"], np.nan,
|
||
vorticity(_af[len(_af) // 2][1]))), 99.0)
|
||
_cmapA = plt.get_cmap("RdBu_r").copy(); _cmapA.set_bad("#1b1b1b")
|
||
|
||
def _draw_amr(fig, i):
|
||
t, uc, uf = _af[i]
|
||
ax = fig.add_subplot(111)
|
||
# interpolation="nearest" — чтобы были ВИДНЫ ячейки: крупные снаружи, вдвое
|
||
# мельче в блоке (иначе imshow сглаживает и разница незаметна).
|
||
ax.imshow(np.where(amr["solid_c"], np.nan, vorticity(uc)), origin="lower",
|
||
cmap=_cmapA, vmin=-_vmaxA, vmax=_vmaxA, extent=(0, amr["Nx"], 0, amr["Ny"]),
|
||
interpolation="nearest")
|
||
omf = np.where(amr["solid_f"], np.nan, vorticity(uf))[1:-1, 1:-1] # без ghost-рамки
|
||
ax.imshow(omf, origin="lower", cmap=_cmapA, vmin=-_vmaxA, vmax=_vmaxA,
|
||
extent=(amr["bx0"] + 0.5, amr["bx1"] - 0.5, amr["by0"] + 0.5, amr["by1"] - 0.5),
|
||
interpolation="nearest")
|
||
ax.add_patch(mpatches.Rectangle((amr["bx0"], amr["by0"]),
|
||
amr["bx1"] - amr["bx0"], amr["by1"] - amr["by0"],
|
||
fill=False, edgecolor="#ffe600", lw=3.0))
|
||
ax.text(amr["bx0"] + 1, amr["by1"] - 4, "мелкий блок 2×", color="#ffe600", fontsize=9, weight="bold")
|
||
ax.set_title(f"AMR: крупная сетка + мелкий блок 2× (жёлтая рамка) · t={t}", fontsize=10)
|
||
ax.set_xticks([]); ax.set_yticks([])
|
||
fig.tight_layout()
|
||
|
||
show_gif(render_gif(len(_af), _draw_amr, "amr_cylinder.gif", (9.2, 4.7), fps=24), width=780)
|
||
|
||
# %% [markdown]
|
||
# ### AMR на цилиндре на БОЛЬШОЙ области — обе реализации (как на схеме исследования)
|
||
#
|
||
# Эти две анимации прогнаны отдельным скриптом **`_amr_anims.py`** (большая область
|
||
# 150×70, длинный прогон ~7200 шагов → 240 кадров @24fps) и встроены готовыми GIF —
|
||
# чтобы тяжёлый счёт не раздувал рендер ноутбука. На большой области видна и **грубая
|
||
# сетка** (уровень 0, Δx), и **учащённые зоны**: зелёный блок **2× (тело+след)** и
|
||
# оранжевый **4× (у тела)**. `interpolation="nearest"` — поэтому размеры ячеек различимы.
|
||
|
||
# %%
|
||
from IPython.display import HTML as _HTMLbig, display as _dispbig
|
||
if IN_NOTEBOOK:
|
||
_dispbig(_HTMLbig(
|
||
'<b>1) Одиночный 2× AMR на цилиндре (большая область):</b><br>'
|
||
'<img src="anim/amr2x_cyl.gif" width="900"/>'
|
||
'<br><b>2) Полный вложенный 2×+4× AMR на цилиндре (большая область):</b><br>'
|
||
'<img src="anim/amr_nested_cyl.gif" width="900"/>'))
|
||
|
||
|
||
# %% [markdown]
|
||
# ## Полный AMR: вложенный 2×+4× с временно́й интерполяцией и численное сравнение
|
||
#
|
||
# Расширяем минимальный AMR до **полного вложенного** (как в исследовании — зона **2×
|
||
# для окружения/хвоста** и **4× для завихрений у тела**):
|
||
# - три уровня: коарс → **2×** (большой блок) → **4×** (вложенный блок у тела);
|
||
# - **рекурсивный sub-cycling 1:2:4** — на 1 крупный шаг 2 подшага уровня-2× и 4
|
||
# подшага уровня-4× (шаг по времени делится на 2 при каждом измельчении);
|
||
# - **временна́я интерполяция**: граница каждого блока между подшагами линейно
|
||
# интерполируется во времени между состояниями родителя в моменты $t$ и $t+\Delta t$
|
||
# (минимальный AMR держал границу постоянной — это вносило ошибку);
|
||
# - перемасштаб на каждой границе: $\tau_{child}=2\tau_{parent}-\tfrac12$,
|
||
# $f^{neq}\!\cdot\!\tau_{child}/(2\tau_{parent})$.
|
||
#
|
||
# **Методика сравнения.** Контролируемый тест — локализованный резкий вихрь (его ядро
|
||
# крупная сетка разрешает плохо). Эталон точности — равномерная **мелкая 4×** сетка.
|
||
# Сравниваем: равномерную крупную, минимальный 2× (без врем.интерп.), одиночный 4×,
|
||
# вложенный 2×+4× без и с врем.интерп. **Производительность меряем в cell-updates**
|
||
# (объём работы — аппаратно-независимо); wall-time в чистом numpy искажён Python-
|
||
# оверхедом мелких блоков, поэтому приводим его лишь справочно.
|
||
# Коэффициент производительности/точности: $\eta = 1/(\varepsilon\cdot W_6)$, где
|
||
# $\varepsilon$ — отн. L2-ошибка к эталону, $W_6$ — cell-updates в миллионах.
|
||
|
||
# %%
|
||
import time as _time
|
||
|
||
_Nc = 48; _L = 48; _xc = _yc = 24.0; _sig = 1.4; _A = 0.115
|
||
_nu_c = 0.008; _tau_c = _nu_c / CS2 + 0.5; _BSTEPS = 400
|
||
|
||
def _vfield(X, Y):
|
||
rr = (X - _xc) ** 2 + (Y - _yc) ** 2; e = np.exp(-rr / (2 * _sig ** 2))
|
||
return np.stack([-_A * (Y - _yc) / _sig ** 2 * e, _A * (X - _xc) / _sig ** 2 * e])
|
||
|
||
def _grid_v(N, h):
|
||
xs = np.arange(N) * h; Xg, Yg = np.meshgrid(xs, xs); return _vfield(Xg, Yg)
|
||
|
||
def _adv(f, beta):
|
||
r, u = macros(f); return stream(collide_kbc(f, feq_entropic(r, u), beta)[0])
|
||
|
||
def _patch(pN, a, b, r):
|
||
Wc = b - a; Nf = r * Wc + 1; pos = a + np.arange(Nf) / r; FX, FY = np.meshgrid(pos, pos)
|
||
x0 = np.floor(FX).astype(int); y0 = np.floor(FY).astype(int)
|
||
x1 = np.minimum(x0 + 1, pN - 1); y1 = np.minimum(y0 + 1, pN - 1)
|
||
return dict(Nf=Nf, a=a, b=b, Wc=Wc, r=r, x0=x0, y0=y0, x1=x1, y1=y1,
|
||
tx=FX - x0, ty=FY - y0, sl=slice(r, r * Wc, r))
|
||
|
||
def _pint(fld, P):
|
||
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):
|
||
r, u = macros(pf); neq = pf - feq_entropic(r, u)
|
||
return feq_entropic(_pint(r, P), np.stack([_pint(u[0], P), _pint(u[1], P)])) + Rcf * _pint(neq, P)
|
||
|
||
def _fill(cf, gh):
|
||
cf[:, 0, :] = gh[:, 0, :]; cf[:, -1, :] = gh[:, -1, :]; cf[:, :, 0] = gh[:, :, 0]; cf[:, :, -1] = gh[:, :, -1]; return cf
|
||
|
||
def _restrict(cf, pf, P, Rfc):
|
||
r, u = macros(cf); neq = cf - feq_entropic(r, u); sl = P["sl"]
|
||
pf[:, P["a"] + 1:P["b"], P["a"] + 1:P["b"]] = feq_entropic(r[sl, sl], u[:, sl, sl]) + Rfc * neq[:, sl, sl]; return pf
|
||
|
||
def _tauc(tp, r): return r * tp - (r - 1) / 2
|
||
|
||
def _run_uniform(N, h, tau, nsteps):
|
||
beta = 1 / (2 * tau); f = feq_entropic(np.ones((N, N)), _grid_v(N, h)); t0 = _time.perf_counter()
|
||
for _ in range(nsteps): f = _adv(f, beta)
|
||
return macros(f)[1], _time.perf_counter() - t0, N * N * nsteps
|
||
|
||
def _init_patch(P, a_cc, h_par):
|
||
pos = a_cc + np.arange(P["Nf"]) * (h_par / P["r"]); X, Y = np.meshgrid(pos, pos)
|
||
return feq_entropic(np.ones((P["Nf"], P["Nf"])), _vfield(X, Y))
|
||
|
||
def _run_single(r, a, b, temporal):
|
||
P = _patch(_Nc, a, b, r); tau_f = _tauc(_tau_c, r)
|
||
b0, bf = 1 / (2 * _tau_c), 1 / (2 * tau_f); Rcf = tau_f / (r * _tau_c); Rfc = 1 / Rcf
|
||
f0 = feq_entropic(np.ones((_Nc, _Nc)), _grid_v(_Nc, 1.0)); ff = _init_patch(P, a, 1.0)
|
||
t0 = _time.perf_counter()
|
||
for t in range(_BSTEPS):
|
||
f0o = f0.copy(); f0 = _adv(f0, b0)
|
||
gho = _ghost(f0o, P, Rcf); ghn = _ghost(f0, P, Rcf)
|
||
for s in range(r):
|
||
ff = _adv(ff, bf); ff = _fill(ff, (1 - (s + 1) / r) * gho + ((s + 1) / r) * ghn if temporal else ghn)
|
||
f0 = _restrict(ff, f0, P, Rfc)
|
||
return macros(f0)[1], _time.perf_counter() - t0, _Nc * _Nc * _BSTEPS + P["Nf"] ** 2 * r * _BSTEPS
|
||
|
||
def _run_nested(temporal):
|
||
P1 = _patch(_Nc, 8, 40, 2); P2 = _patch(P1["Nf"], 16, 48, 2)
|
||
t1 = _tauc(_tau_c, 2); t2 = _tauc(t1, 2)
|
||
b0, b1, b2 = 1 / (2 * _tau_c), 1 / (2 * t1), 1 / (2 * t2)
|
||
R01 = t1 / (2 * _tau_c); Rf01 = 1 / R01; R12 = t2 / (2 * t1); Rf12 = 1 / R12
|
||
f0 = feq_entropic(np.ones((_Nc, _Nc)), _grid_v(_Nc, 1.0))
|
||
f1 = _init_patch(P1, 8, 1.0); f2 = _init_patch(P2, 16.0, 0.5)
|
||
t0 = _time.perf_counter()
|
||
for t in range(_BSTEPS):
|
||
f0o = f0.copy(); f0 = _adv(f0, b0)
|
||
g1o = _ghost(f0o, P1, R01); g1n = _ghost(f0, P1, R01)
|
||
for s1 in range(2):
|
||
f1o = f1.copy(); f1 = _adv(f1, b1)
|
||
f1 = _fill(f1, (1 - (s1 + 1) / 2) * g1o + ((s1 + 1) / 2) * g1n if temporal else g1n)
|
||
g2o = _ghost(f1o, P2, R12); g2n = _ghost(f1, P2, R12)
|
||
for s2 in range(2):
|
||
f2 = _adv(f2, b2); f2 = _fill(f2, (1 - (s2 + 1) / 2) * g2o + ((s2 + 1) / 2) * g2n if temporal else g2n)
|
||
f1 = _restrict(f2, f1, P2, Rf12)
|
||
f0 = _restrict(f1, f0, P1, Rf01)
|
||
cells = _Nc * _Nc * _BSTEPS + P1["Nf"] ** 2 * 2 * _BSTEPS + P2["Nf"] ** 2 * 4 * _BSTEPS
|
||
return macros(f0)[1], _time.perf_counter() - t0, cells
|
||
|
||
# эталон: равномерная мелкая 4× (то же физ. время → 4×шагов)
|
||
_uref, _tref, _cref = _run_uniform(_Nc * 4, 0.25, _tauc(_tau_c, 4), 4 * _BSTEPS)
|
||
_uref = _uref[:, ::4, ::4]; _rms = np.sqrt(np.mean(_uref[0] ** 2 + _uref[1] ** 2))
|
||
def _err(u): return np.sqrt(np.mean((u[0] - _uref[0]) ** 2 + (u[1] - _uref[1]) ** 2)) / _rms
|
||
|
||
_bench = []
|
||
u, t, c = _run_uniform(_Nc, 1.0, _tau_c, _BSTEPS); _bench.append(("uniform-coarse", _err(u), t, c))
|
||
u, t, c = _run_single(2, 12, 36, False); _bench.append(("AMR 2× (мин., без вр.инт.)", _err(u), t, c))
|
||
u, t, c = _run_single(4, 12, 36, True); _bench.append(("AMR 4× single (вр.инт.)", _err(u), t, c))
|
||
u, t, c = _run_nested(False); _bench.append(("AMR 2×+4× nested (без вр.инт.)", _err(u), t, c))
|
||
u, t, c = _run_nested(True); _bench.append(("AMR 2×+4× nested ПОЛНЫЙ", _err(u), t, c))
|
||
|
||
print(f"Эталон uniform-fine 4× (192²): время={_tref:.1f}c, cell-updates={_cref:.2e}")
|
||
print(f"{'метод':32}{'ε,%':>8}{'cell-upd':>11}{'η=1/(εW)':>10}{'×fine(работа)':>14}")
|
||
for name, e, t, c in _bench:
|
||
eta = 1.0 / (e * c / 1e6)
|
||
print(f"{name:32}{e*100:>8.2f}{c:>11.2e}{eta:>10.2f}{_cref/c:>14.1f}")
|
||
print("\nη = 1/(ε·W6), W6 — cell-updates в млн. ×fine(работа) — во сколько раз меньше работы, чем у эталона.")
|
||
|
||
# %%
|
||
# Графики сравнения: ошибка по методам и фронт «ошибка vs работа».
|
||
names = [b[0] for b in _bench]; errs = np.array([b[1] * 100 for b in _bench])
|
||
cells = np.array([b[3] for b in _bench]); cols = ["#9e9e9e", "#1f6feb", "#2e7d32", "#fb8c00", "#d1495b"]
|
||
fig, axs = plt.subplots(1, 2, figsize=(13.5, 4.4))
|
||
axs[0].barh(range(len(names)), errs, color=cols)
|
||
axs[0].set_yticks(range(len(names))); axs[0].set_yticklabels(names, fontsize=8); axs[0].invert_yaxis()
|
||
axs[0].set_xlabel("отн. L2-ошибка к эталону, %"); axs[0].set_title("Точность (меньше — лучше)")
|
||
for i, e in enumerate(errs): axs[0].text(e + 0.05, i, f"{e:.2f}", va="center", fontsize=8)
|
||
axs[1].scatter(cells, errs, c=cols, s=90, zorder=3, edgecolors="k")
|
||
for n, c2, e in zip(names, cells, errs):
|
||
axs[1].annotate(n, (c2, e), fontsize=7, xytext=(5, 4), textcoords="offset points")
|
||
axs[1].axvline(_cref, color="#d1495b", ls="--", alpha=0.6); axs[1].text(_cref, errs.max()*0.9, " работа эталона", color="#d1495b", fontsize=8)
|
||
axs[1].set_xscale("log"); axs[1].set_xlabel("cell-updates (работа, лог)"); axs[1].set_ylabel("ошибка, %")
|
||
axs[1].set_title("Фронт «точность–работа»: AMR ближе к низ-лево (точно и дёшево)")
|
||
fig.suptitle("AMR: численное сравнение точности и производительности (KBC, тест-вихрь)", fontweight="bold")
|
||
_finish(fig, "11d_amr_comparison.png")
|
||
|
||
# %% [markdown]
|
||
# **Численные результаты (типичные значения; точные печатаются ячейкой выше):**
|
||
#
|
||
# | метод | ε (L2) | cell-updates | работа vs эталон |
|
||
# |---|---|---|---|
|
||
# | uniform-coarse | ~5.3 % | ~9.2·10⁵ | ×64 меньше |
|
||
# | AMR 2× (мин., без вр.инт.) | ~2.6 % | ~2.8·10⁶ | ×21 меньше |
|
||
# | AMR 4× single (вр.инт.) | ~2.2 % | ~1.6·10⁷ | ×3.7 меньше |
|
||
# | AMR 2×+4× nested (без вр.инт.) | ~2.3 % | ~1.1·10⁷ | ×5.3 меньше |
|
||
# | **AMR 2×+4× nested ПОЛНЫЙ** | **~2.15 %** | **~1.1·10⁷** | **×5.3 меньше** |
|
||
# | uniform-fine 4× (эталон) | 0 (опорная) | ~5.9·10⁷ | — |
|
||
#
|
||
# **Выводы:**
|
||
# - **Вложенность эффективна:** полный 2×+4× точнее одиночного 4× (≈2.15 % против
|
||
# ≈2.2 %) при **в ~1.45× меньшей работе** — 4× только у ядра, 2× вокруг.
|
||
# - **Временна́я интерполяция помогает:** при той же работе ошибка падает (≈2.26 %→
|
||
# ≈2.15 %, на ~5 % относительно).
|
||
# - **AMR vs эталон:** полный AMR даёт точность ~2 % за **в ~5× меньше работы**, чем
|
||
# равномерная мелкая сетка.
|
||
# - **Крупная сетка** дёшева, но её точность (~5 %) принципиально хуже и не улучшается
|
||
# без измельчения.
|
||
# - ⚠ wall-time в чистом numpy искажён Python-оверхедом мелких блоков; на GPU/в
|
||
# оптимизированном коде время ∝ cell-updates, и преимущество AMR по времени сохраняется.
|
||
|
||
# %% [markdown]
|
||
# ### Анимация полного вложенного AMR (видны три разрешения: коарс / 2× / 4×)
|
||
#
|
||
# Тот же вложенный солвер записывает кадры. На анимации `interpolation="nearest"` —
|
||
# поэтому **видны размеры ячеек**: крупные снаружи, вдвое мельче в зелёном блоке (2×),
|
||
# ещё вдвое мельче в оранжевом (4×). Так наглядно работает разделение «2× для
|
||
# окружения / 4× у ядра».
|
||
|
||
# %%
|
||
def _run_nested_anim(steps=750, frame_every=3):
|
||
P1 = _patch(_Nc, 8, 40, 2); P2 = _patch(P1["Nf"], 16, 48, 2)
|
||
t1 = _tauc(_tau_c, 2); t2 = _tauc(t1, 2)
|
||
b0, b1, b2 = 1 / (2 * _tau_c), 1 / (2 * t1), 1 / (2 * t2)
|
||
R01 = t1 / (2 * _tau_c); Rf01 = 1 / R01; R12 = t2 / (2 * t1); Rf12 = 1 / R12
|
||
f0 = feq_entropic(np.ones((_Nc, _Nc)), _grid_v(_Nc, 1.0))
|
||
f1 = _init_patch(P1, 8, 1.0); f2 = _init_patch(P2, 16.0, 0.5)
|
||
frames = []
|
||
for t in range(steps + 1):
|
||
f0o = f0.copy(); f0 = _adv(f0, b0)
|
||
g1o = _ghost(f0o, P1, R01); g1n = _ghost(f0, P1, R01)
|
||
for s1 in range(2):
|
||
f1o = f1.copy(); f1 = _adv(f1, b1); f1 = _fill(f1, (1 - (s1 + 1) / 2) * g1o + ((s1 + 1) / 2) * g1n)
|
||
g2o = _ghost(f1o, P2, R12); g2n = _ghost(f1, P2, R12)
|
||
for s2 in range(2):
|
||
f2 = _adv(f2, b2); f2 = _fill(f2, (1 - (s2 + 1) / 2) * g2o + ((s2 + 1) / 2) * g2n)
|
||
f1 = _restrict(f2, f1, P2, Rf12)
|
||
f0 = _restrict(f1, f0, P1, Rf01)
|
||
if t % frame_every == 0:
|
||
frames.append((t, macros(f0)[1], macros(f1)[1], macros(f2)[1]))
|
||
return frames
|
||
|
||
_naf = _run_nested_anim()
|
||
_vmaxN = np.nanpercentile(np.abs(vorticity(_naf[0][1])), 99.0)
|
||
_cN = plt.get_cmap("RdBu_r").copy(); _cN.set_bad("#1b1b1b")
|
||
|
||
def _draw_nested(fig, i):
|
||
t, u0, u1, u2 = _naf[i]
|
||
ax = fig.add_subplot(111)
|
||
ax.imshow(vorticity(u0), origin="lower", cmap=_cN, vmin=-_vmaxN, vmax=_vmaxN,
|
||
extent=(0, _Nc, 0, _Nc), interpolation="nearest") # коарс: крупные ячейки
|
||
ax.imshow(vorticity(u1)[1:-1, 1:-1], origin="lower", cmap=_cN, vmin=-_vmaxN, vmax=_vmaxN,
|
||
extent=(8.5, 39.5, 8.5, 39.5), interpolation="nearest") # 2×: вдвое мельче
|
||
ax.imshow(vorticity(u2)[1:-1, 1:-1], origin="lower", cmap=_cN, vmin=-_vmaxN, vmax=_vmaxN,
|
||
extent=(16.25, 31.75, 16.25, 31.75), interpolation="nearest") # 4×: ещё вдвое мельче
|
||
ax.add_patch(mpatches.Rectangle((8, 8), 32, 32, fill=False, edgecolor="#7CFC00", lw=2.5))
|
||
ax.add_patch(mpatches.Rectangle((16, 16), 16, 16, fill=False, edgecolor="#ff8a00", lw=2.5))
|
||
ax.text(8.3, 40.4, "2× блок", color="#7CFC00", fontsize=9, weight="bold")
|
||
ax.text(16.3, 32.4, "4× блок", color="#ff8a00", fontsize=9, weight="bold")
|
||
ax.set_title(f"Полный вложенный AMR: коарс + 2× + 4× (видны размеры ячеек) · t={t}", fontsize=10)
|
||
ax.set_xticks([]); ax.set_yticks([]); fig.tight_layout()
|
||
|
||
show_gif(render_gif(len(_naf), _draw_nested, "amr_nested.gif", (6.2, 6.2), fps=24), width=620)
|
||
|
||
|
||
# %% [markdown]
|
||
# ## 10. Сопутствующие технологии для полной симуляции
|
||
#
|
||
# KBC задаёт только **оператор столкновения**. Полный решатель требует ещё ряда
|
||
# компонентов.
|
||
#
|
||
# ### 10.1 Граничные условия
|
||
# - **Bounce-back (отскок)** — стенка без проскальзывания (no-slip): входящая
|
||
# популяция отражается, $f_{\bar i}=f_i$. Первый порядок (стенка на полпути между
|
||
# узлами). Для **движущейся** стенки добавляется импульс
|
||
# $-2 w_i\rho\,(\mathbf c_i\cdot\mathbf u_w)/c_s^2$ (использован в каверне).
|
||
# - **Zou–He** — заданные скорость/давление на входе/выходе: недостающие популяции
|
||
# восстанавливаются из условий на $\rho,\mathbf u$ и «bounce-back неравновесной
|
||
# части». Второй порядок.
|
||
# - **Grad / регуляризованные** — реконструкция неравновесной части через тензор
|
||
# напряжений $\Pi^{neq}$; устойчивее при больших градиентах (рекомендуются с
|
||
# энтропийными моделями).
|
||
# - **Интерполированный bounce-back (Bouzidi, Yu–Mei–Luo–Shyy)** — для кривых
|
||
# стенок: используется доля стенки $q\in[0,1]$ из signed-distance field (SDF),
|
||
# две ветви $q\!\ge\!0.5$ и $q\!<\!0.5$.
|
||
# - **Sponge / absorbing-слой** — у границ плавно подмешивается целевое равновесие,
|
||
# чтобы гасить отражения волн (использовался в прежней реализации на 6 гранях).
|
||
#
|
||
# ### 10.2 Объёмные силы (Guo)
|
||
# Схема Guo et al. (2002): к столкновению добавляется
|
||
# $F_i = \bigl(1-\tfrac{\beta}{...}\bigr)w_i\Bigl[\tfrac{\mathbf c_i-\mathbf u}{c_s^2}
|
||
# +\tfrac{(\mathbf c_i\cdot\mathbf u)}{c_s^4}\mathbf c_i\Bigr]\cdot\mathbf F$, а
|
||
# скорость берётся со сдвигом $\mathbf u\to\mathbf u+\tfrac{\mathbf F}{2\rho}$.
|
||
# Даёт второй порядок и точное сохранение импульса.
|
||
#
|
||
# ### 10.3 Турбулентность: implicit LES vs Smagorinsky
|
||
# - **KBC — implicit LES**: адаптивный $\gamma$ добавляет диссипацию высших мод
|
||
# ровно там, где растут неразрешённые масштабы — параметр-free.
|
||
# - **Smagorinsky (явная LES)**: вихревая вязкость $\nu_t=(C_s\Delta)^2|S|$,
|
||
# $|S|=\sqrt{2 S_{\alpha\beta}S_{\alpha\beta}}$ из неравновесного тензора
|
||
# напряжений; эффективное $\tau_{\mathrm{eff}}=3(\nu+\nu_t)+\tfrac12$. $C_s\!\approx\!0.1$–$0.17$.
|
||
# Иногда комбинируют с KBC у стенок.
|
||
#
|
||
# ### 10.4 Единицы и обезразмеривание (lu/ts ↔ СИ)
|
||
# Выбираем масштабы $\Delta x$ (м/lu), $\Delta t$ (с/ts), $\rho_0$ (кг/м³). Тогда
|
||
# $$ \nu_{\text{СИ}} = \nu_{\text{lu}}\,\frac{\Delta x^2}{\Delta t},\quad
|
||
# \mathbf u_{\text{СИ}} = \mathbf u_{\text{lu}}\,\frac{\Delta x}{\Delta t},\quad
|
||
# p_{\text{СИ}} = c_s^2\rho_{\text{lu}}\,\rho_0\frac{\Delta x^2}{\Delta t^2}. $$
|
||
# В решёточных единицах $\Delta x=\Delta t=1$ (узел и шаг — за единицу), поэтому
|
||
# скорости и вязкости — безразмерные «решёточные». Ограничение сжимаемости:
|
||
# $\mathrm{Ma}=|\mathbf u|/c_s\lesssim 0.1$–$0.3$ (ошибка $\propto\mathrm{Ma}^2$);
|
||
# отсюда Mach-лимитер шага $\Delta t \le \mathrm{Ma}\,c_s\,\Delta x/|\mathbf u|$.
|
||
#
|
||
# ### 10.5 Решётки, инициализация, измельчение сетки
|
||
# - **Решётки**: D2Q9 (2D), D3Q19 (3D, экономная), **D3Q27** (3D, полная изотропия
|
||
# 4-го порядка и галилеева инвариантность — предпочтительна для KBC/турбулентности).
|
||
# - **Инициализация**: задать $f=f^{eq}(\rho,\mathbf u)$; для согласованности с
|
||
# градиентами — добавить неравновесную часть из $\Pi^{neq}$ (схема Mei).
|
||
# - **Измельчение сетки**: многосеточные блоки с пространственно-временной
|
||
# интерполяцией; KBC устойчив к смене разрешения (нет явных подсеточных констант).
|
||
|
||
# %% [markdown]
|
||
# ## 11. Перенос на D3Q27 и связь с целевой реализацией
|
||
#
|
||
# Вся математика KBC переносится на 3D **без изменений по сути**: меняются
|
||
# решётка (27 скоростей, веса $8/27,\,2/27,\,1/54,\,1/216$), число моментов
|
||
# (моментный базис 27×27) и состав сдвиговой части $s$ (девиаторный тензор
|
||
# напряжений 3D — 5 независимых компонент). Формула $\gamma$ и оператор
|
||
# $f\leftarrow f-\beta(2\Delta s+\gamma\Delta h)$ идентичны.
|
||
|
||
# %%
|
||
# D3Q27: скорости и веса.
|
||
C3 = np.array([[x, y, z] for x in (-1, 0, 1) for y in (-1, 0, 1) for z in (-1, 0, 1)], dtype=float)
|
||
nz = (C3 != 0).sum(1)
|
||
W3 = np.select([nz == 0, nz == 1, nz == 2, nz == 3], [8/27, 2/27, 1/54, 1/216])
|
||
assert np.isclose(W3.sum(), 1.0)
|
||
assert np.allclose((W3[:, None, None] * C3[:, :, None] * C3[:, None, :]).sum(0), CS2 * np.eye(3))
|
||
print("D3Q27: Σw=1 и Σw c c = c_s² I — OK; число направлений:", len(C3))
|
||
|
||
fig = plt.figure(figsize=(6.4, 5.6))
|
||
ax = fig.add_subplot(111, projection="3d")
|
||
tier3 = {8/27: ("#d1495b", "покой"), 2/27: ("#2e7d32", "оси (×6)"),
|
||
1/54: ("#1f6feb", "рёбра (×12)"), 1/216: ("#9c27b0", "углы (×8)")}
|
||
for wv, (col, lab) in tier3.items():
|
||
m = np.isclose(W3, wv)
|
||
ax.scatter(C3[m, 0], C3[m, 1], C3[m, 2], s=70 + 5000 * wv, color=col, label=f"{lab}, w={wv:.4f}",
|
||
depthshade=True, edgecolors="k", linewidths=0.4)
|
||
ax.set_xlabel("$c_x$"); ax.set_ylabel("$c_y$"); ax.set_zlabel("$c_z$")
|
||
ax.set_title("Решётка D3Q27 (цель реализации)")
|
||
ax.legend(loc="upper left", fontsize=7, bbox_to_anchor=(0.0, 1.0))
|
||
_finish(fig, "11_d3q27.png")
|
||
|
||
# %% [markdown]
|
||
# ### Чем «настоящий KBC» отличается от прежнего подхода
|
||
#
|
||
# Прежняя реализация (GPU-решатель) фактически использовала **полиномиальное**
|
||
# равновесие + **TRT** (две скорости релаксации) + **α-лимитер положительности**
|
||
# (ELBM-типа, по симметричной части) + явный **Smagorinsky**. Это устойчивая и
|
||
# рабочая схема, но **не** KBC. Переход к настоящему KBC:
|
||
#
|
||
# | Аспект | Настоящий KBC | Прежний подход (TRT + α) |
|
||
# |---|---|---|
|
||
# | Равновесие | энтропийное product-form | полиномиальное $O(u^2)$ |
|
||
# | Столкновение | энтропийный MRT, $k/s/h$ | TRT (симм./антисимм.) |
|
||
# | Стабилизация | $\gamma$ из энтропии (высшие моды) | α-лимитер положительности (симм. часть) |
|
||
# | H-теорема | встроена в вывод $\gamma$ | контролируется отдельно (диагностика) |
|
||
# | Турбулентность | implicit LES (параметр-free) | явный Smagorinsky |
|
||
# | Решётка | D3Q27 | D3Q27 (совпадает) |
|
||
#
|
||
# **Дорожная карта реализации на GPU**: (1) энтропийное $f^{eq}$ (product-form);
|
||
# (2) моментные проекторы $P_s,P_h$ для D3Q27 (предвычислить матрицы);
|
||
# (3) per-node $\gamma$ по формуле; (4) коллизия $f-\beta(2\Delta s+\gamma\Delta h)$;
|
||
# (5) опционально — комбинировать с пристеночным Smagorinsky.
|
||
|
||
# %% [markdown]
|
||
# ## 12. Формулы под сверку с оригиналом (⚠) и литература
|
||
#
|
||
# **Под сверку с первоисточниками** (доступ к части PDF был ограничен — проверить
|
||
# коэффициенты по оригиналам):
|
||
# 1. **Энтропийное равновесие** product-form: префактор $(2-\sqrt{1+3u_\alpha^2})$ и
|
||
# основание $\frac{2u_\alpha+\sqrt{1+3u_\alpha^2}}{1-u_\alpha}$ — сверить вид и
|
||
# область применимости (Ansumali–Karlin; Chikatamarla–Karlin).
|
||
# *Здесь проверено численно: даёт $\rho,\rho\mathbf u$ точно.*
|
||
# 2. **Формула $\gamma$** $=\tfrac1\beta-(2-\tfrac1\beta)\tfrac{\langle\Delta s|\Delta h\rangle}{\langle\Delta h|\Delta h\rangle}$
|
||
# и порядок разложения $\Delta H$, из которого она получена (Bösch et al. 2015).
|
||
# *Здесь проверено: предел $\gamma=2$ совпадает с BGK.*
|
||
# 3. **Назначение моментов** в $s$ vs $h$ для вариантов KBC-N1…N4 (D2Q9/D3Q27) —
|
||
# сверить точный состав по таблицам оригинала.
|
||
# 4. **Множитель 2** перед $\Delta s$ и связь $\beta=1/(2\tau)$ — сверить нормировку.
|
||
#
|
||
# **Литература (первоисточники):**
|
||
# - I. V. Karlin, F. Bösch, S. S. Chikatamarla, *Gibbs' principle for the
|
||
# lattice-kinetic theory of fluid dynamics*, **Phys. Rev. E 90**, 031302(R) (2014).
|
||
# - F. Bösch, S. S. Chikatamarla, I. V. Karlin, *Entropic multirelaxation lattice
|
||
# Boltzmann models for turbulent flows*, **Phys. Rev. E 92**, 043309 (2015).
|
||
# - S. S. Chikatamarla, I. V. Karlin, *Entropy and Galilean invariance of lattice
|
||
# Boltzmann theories*, **Phys. Rev. Lett.** (2006) — энтропийное равновесие.
|
||
# - B. Dorschner, S. S. Chikatamarla, I. V. Karlin и др. — применения KBC
|
||
# (каналы, обтекание тел), измельчение сетки, подвижные границы.
|
||
# - Z. Guo, C. Zheng, B. Shi, *Discrete lattice effects on the forcing term in the
|
||
# lattice Boltzmann method*, **Phys. Rev. E 65**, 046308 (2002) — форсинг.
|
||
# - U. Ghia, K. N. Ghia, C. T. Shin (1982) — эталон каверны (валидация).
|
||
|
||
# %% [markdown]
|
||
# ## 13. Итоги
|
||
#
|
||
# - **KBC** = энтропийный MRT-LBM: сохраняющиеся моды не трогаем, сдвиг релаксируем
|
||
# фиксированным $2\beta$ (задаёт вязкость), высшие моды — адаптивным $\gamma$ из
|
||
# условия неувеличения энтропии. При $\gamma=2$ — это в точности BGK.
|
||
# - **Зачем**: безусловная нелинейная устойчивость и параметр-free implicit LES —
|
||
# считает турбулентность на грубых сетках там, где BGK разрушается
|
||
# (продемонстрировано на сдвиговом слое).
|
||
# - **Проверено численно**: точность равновесия и разложения (до $10^{-15}$),
|
||
# правильная вязкость (TGV, ошибка $<0.1\%$), устойчивость (сдвиговый слой),
|
||
# разумное совпадение с эталоном Ghia (каверна).
|
||
# - **Дальше**: перенос на D3Q27 на GPU (Vulkan), сверка коэффициентов с
|
||
# оригинальными статьями, при необходимости — комбинация с Smagorinsky и
|
||
# точными ГУ (Grad/Zou–He, Bouzidi для кривых стенок).
|
||
|
||
# %%
|
||
print("\n=== Все самопроверки и валидации пройдены ===")
|