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>
266 lines
14 KiB
Python
266 lines
14 KiB
Python
# Демо-5: блочное измельчение (AMR) на компонентах solver_2x_sdf (amr.patch/ghost/fill/restrict —
|
||
# ТЕ ЖЕ функции, что в финальной модели). Три части:
|
||
# 1) минимальный 2-уровневый AMR на цилиндре (без временно́й интерполяции) → anim/amr_cylinder.gif
|
||
# 2) бенчмарк 5 вариантов против эталона uniform-fine 4× (тест-вихрь) → figures/11d_amr_comparison.png
|
||
# 3) анимация полного вложенного 2×+4× с временно́й интерполяцией → anim/amr_nested.gif
|
||
# Числа → out/amr_bench.txt|npz. Запуск: python amr_demo_gpu.py
|
||
import os
|
||
import time as _time
|
||
import numpy as np
|
||
import _common as cm
|
||
import amr as Amr # из solver_2x_sdf (sys.path настроен в _common)
|
||
|
||
cp = cm.cp; xp = cm.xp; E = cm.E; Coll = cm.Coll; Strm = cm.Strm
|
||
plt = cm.plt; vorticity = cm.vorticity; OPP = cm.OPP; Q = cm.Q; CS2 = cm.CS2
|
||
import matplotlib.patches as mpatches
|
||
|
||
def _sync():
|
||
cp.cuda.Device().synchronize()
|
||
|
||
# ===== часть 1: минимальный 2-уровневый AMR на цилиндре =========================================
|
||
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 = cm.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)
|
||
P = Amr.patch(Nx, Ny, bx0, bx1, by0, by1, 2)
|
||
Nfx, Nfy = P["Nfx"], P["Nfy"]
|
||
# маски тела (host → device); стенки канала входят в solid_c
|
||
Yc_, Xc_ = np.mgrid[0:Ny, 0:Nx]
|
||
solid_c_np = (Xc_ - cx)**2 + (Yc_ - cy)**2 <= (D/2)**2
|
||
solid_c_np[0, :] = True; solid_c_np[-1, :] = True
|
||
fcx, fcy = 2*(cx - bx0), 2*(cy - by0)
|
||
Yf_, Xf_ = np.mgrid[0:Nfy, 0:Nfx]
|
||
solid_f_np = (Xf_ - fcx)**2 + (Yf_ - fcy)**2 <= float(D)**2
|
||
solid_c = cp.asarray(solid_c_np); solid_f = cp.asarray(solid_f_np)
|
||
block_fluid = cp.asarray(~solid_c_np[by0+1:by1, bx0+1:bx1])
|
||
onecol = cp.ones((Ny, 1), cm.DTYPE)
|
||
fc = E.feq(cp.ones((Ny, Nx), cm.DTYPE), cp.zeros((2, Ny, Nx), cm.DTYPE))
|
||
ff = E.feq(cp.ones((Nfy, Nfx), cm.DTYPE), cp.zeros((2, Nfy, Nfx), cm.DTYPE))
|
||
frames = []
|
||
prog = cm.Viz.make_progress("amr-cyl")
|
||
for t in range(steps + 1):
|
||
rho, u = E.macros(fc)
|
||
pre = fc.copy()
|
||
fc = Coll.kbc_collide(fc, E.feq(rho, u), beta_c)
|
||
for i in range(Q):
|
||
fc[i, solid_c] = pre[OPP[i], solid_c]
|
||
fc = Strm.stream(fc)
|
||
uin = cp.zeros((2, Ny, 1), cm.DTYPE); uin[0, :, 0] = U*cm.smoothstep(t/ramp)
|
||
fc[:, 1:-1, 0] = E.feq(onecol, uin)[:, 1:-1, 0]
|
||
fc[:, 1:-1, -1] = fc[:, 1:-1, -2]
|
||
# МИНИМАЛЬНЫЙ AMR: один ghost на оба подшага (без временно́й интерполяции)
|
||
gh = Amr.ghost(fc, P, Rcf)
|
||
for _ in range(2):
|
||
rf, uf = E.macros(ff); pref = ff.copy()
|
||
ff = Coll.kbc_collide(ff, E.feq(rf, uf), beta_f)
|
||
for i in range(Q):
|
||
ff[i, solid_f] = pref[OPP[i], solid_f]
|
||
ff = Strm.stream(ff)
|
||
ff = Amr.fill(ff, gh)
|
||
fc = Amr.restrict(ff, fc, P, Rfc, block_fluid)
|
||
if t % frame_every == 0:
|
||
frames.append((t, cm.to_cpu(E.macros(fc)[1]), cm.to_cpu(E.macros(ff)[1])))
|
||
prog(t, steps)
|
||
return dict(frames=frames, solid_c=solid_c_np, solid_f=solid_f_np, 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'])} (устойчиво)")
|
||
|
||
_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)
|
||
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]
|
||
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()
|
||
|
||
cm.render_gif(len(_af), _draw_amr, "amr_cylinder.gif", (9.2, 4.7), fps=24)
|
||
|
||
# ===== часть 2: бенчмарк точность/работа (тест-вихрь, без тел) ==================================
|
||
_Nc = 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 = cp.exp(-rr/(2*_sig**2))
|
||
return cm.B.pack(((-_A*(Y - _yc)/_sig**2*e).astype(cm.DTYPE),
|
||
(_A*(X - _xc)/_sig**2*e).astype(cm.DTYPE)))
|
||
|
||
def _grid_v(N, h):
|
||
xs = cp.arange(N, dtype=cm.DTYPE)*h
|
||
Xg, Yg = cp.meshgrid(xs, xs)
|
||
return _vfield(Xg, Yg)
|
||
|
||
def _adv(f, beta):
|
||
r, u = E.macros(f)
|
||
return Strm.stream(Coll.kbc_collide(f, E.feq(r, u), beta))
|
||
|
||
def _sq_patch(pN, a, b, r):
|
||
return Amr.patch(pN, pN, a, b, a, b, r)
|
||
|
||
def _allfluid(P):
|
||
return cp.ones((P["by"]-P["ay"]-1, P["bx"]-P["ax"]-1), bool)
|
||
|
||
def _tauc(tp, r): return r*tp - (r - 1)/2
|
||
|
||
def _init_patch(P, a_phys, h_par):
|
||
pos = a_phys + cp.arange(P["Nfx"], dtype=cm.DTYPE)*(h_par/P["r"])
|
||
X, Y = cp.meshgrid(pos, pos)
|
||
return E.feq(cp.ones((P["Nfy"], P["Nfx"]), cm.DTYPE), _vfield(X, Y))
|
||
|
||
def _run_uniform(N, h, tau, nsteps):
|
||
beta = 1/(2*tau); f = E.feq(cp.ones((N, N), cm.DTYPE), _grid_v(N, h))
|
||
_sync(); t0 = _time.perf_counter()
|
||
for _ in range(nsteps):
|
||
f = _adv(f, beta)
|
||
_sync()
|
||
return E.macros(f)[1], _time.perf_counter() - t0, N*N*nsteps
|
||
|
||
def _run_single(r, a, b, temporal):
|
||
P = _sq_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
|
||
fl = _allfluid(P)
|
||
f0 = E.feq(cp.ones((_Nc, _Nc), cm.DTYPE), _grid_v(_Nc, 1.0)); ff = _init_patch(P, a, 1.0)
|
||
_sync(); t0 = _time.perf_counter()
|
||
for t in range(_BSTEPS):
|
||
f0o = f0.copy(); f0 = _adv(f0, b0)
|
||
gho = Amr.ghost(f0o, P, Rcf); ghn = Amr.ghost(f0, P, Rcf)
|
||
for s in range(r):
|
||
ff = _adv(ff, bf)
|
||
w = (s + 1)/r
|
||
ff = Amr.fill(ff, (1 - w)*gho + w*ghn if temporal else ghn)
|
||
f0 = Amr.restrict(ff, f0, P, Rfc, fl)
|
||
_sync()
|
||
return E.macros(f0)[1], _time.perf_counter() - t0, _Nc*_Nc*_BSTEPS + P["Nfx"]**2*r*_BSTEPS
|
||
|
||
def _run_nested(temporal, steps=_BSTEPS, frame_every=0):
|
||
P1 = _sq_patch(_Nc, 8, 40, 2); P2 = _sq_patch(P1["Nfx"], 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
|
||
fl1 = _allfluid(P1); fl2 = _allfluid(P2)
|
||
f0 = E.feq(cp.ones((_Nc, _Nc), cm.DTYPE), _grid_v(_Nc, 1.0))
|
||
f1 = _init_patch(P1, 8, 1.0); f2 = _init_patch(P2, 16.0, 0.5)
|
||
frames = []
|
||
_sync(); t0 = _time.perf_counter()
|
||
for t in range(steps + (1 if frame_every else 0)):
|
||
f0o = f0.copy(); f0 = _adv(f0, b0)
|
||
g1o = Amr.ghost(f0o, P1, R01); g1n = Amr.ghost(f0, P1, R01)
|
||
for s1 in range(2):
|
||
f1o = f1.copy(); f1 = _adv(f1, b1)
|
||
w1 = (s1 + 1)/2
|
||
f1 = Amr.fill(f1, (1 - w1)*g1o + w1*g1n if temporal else g1n)
|
||
g2o = Amr.ghost(f1o, P2, R12); g2n = Amr.ghost(f1, P2, R12)
|
||
for s2 in range(2):
|
||
f2 = _adv(f2, b2)
|
||
w2 = (s2 + 1)/2
|
||
f2 = Amr.fill(f2, (1 - w2)*g2o + w2*g2n if temporal else g2n)
|
||
f1 = Amr.restrict(f2, f1, P2, Rf12, fl2)
|
||
f0 = Amr.restrict(f1, f0, P1, Rf01, fl1)
|
||
if frame_every and t % frame_every == 0:
|
||
frames.append((t, cm.to_cpu(E.macros(f0)[1]), cm.to_cpu(E.macros(f1)[1]),
|
||
cm.to_cpu(E.macros(f2)[1])))
|
||
_sync()
|
||
cells = _Nc*_Nc*steps + P1["Nfx"]**2*2*steps + P2["Nfx"]**2*4*steps
|
||
if frame_every:
|
||
return frames
|
||
return E.macros(f0)[1], _time.perf_counter() - t0, cells
|
||
|
||
print("\nБенчмарк AMR (эталон uniform-fine 4×, физ. время одинаково)...")
|
||
_uref, _tref, _cref = _run_uniform(_Nc*4, 0.25, _tauc(_tau_c, 4), 4*_BSTEPS)
|
||
_uref = _uref[:, ::4, ::4]
|
||
_rms = float(cp.sqrt((_uref[0]**2 + _uref[1]**2).mean()))
|
||
|
||
def _err(u):
|
||
return float(cp.sqrt(((u[0] - _uref[0])**2 + (u[1] - _uref[1])**2).mean()))/_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 2x (мин., без вр.инт.)", _err(u), t, c))
|
||
u, t, c = _run_single(4, 12, 36, True); _bench.append(("AMR 4x single (вр.инт.)", _err(u), t, c))
|
||
u, t, c = _run_nested(False); _bench.append(("AMR 2x+4x nested (без вр.инт.)", _err(u), t, c))
|
||
u, t, c = _run_nested(True); _bench.append(("AMR 2x+4x nested ПОЛНЫЙ", _err(u), t, c))
|
||
|
||
out_lines = [f"reference uniform-fine 4x (192^2): time={_tref:.2f}s cell-updates={_cref:.3e}"]
|
||
print(f"Эталон uniform-fine 4× (192²): время={_tref:.1f}с, 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}")
|
||
out_lines.append(f"{name}: eps_pct={e*100:.3f} cells={c:.3e} eta={eta:.2f} work_vs_fine={_cref/c:.1f}")
|
||
|
||
np.savez(os.path.join(cm.OUTDIR, "amr_bench.npz"),
|
||
names=np.array([b[0] for b in _bench]),
|
||
eps=np.array([b[1] for b in _bench]),
|
||
wall=np.array([b[2] for b in _bench]),
|
||
cells=np.array([b[3] for b in _bench]),
|
||
cref=_cref, tref=_tref)
|
||
|
||
# фигура: ошибка по методам + фронт «точность–работа»
|
||
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")
|
||
cm.finish(fig, "11d_amr_comparison.png")
|
||
|
||
# ===== часть 3: анимация полного вложенного 2×+4× ===============================================
|
||
_naf = _run_nested(True, steps=750, frame_every=3)
|
||
_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")
|
||
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")
|
||
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()
|
||
|
||
cm.render_gif(len(_naf), _draw_nested, "amr_nested.gif", (6.2, 6.2), fps=24)
|
||
|
||
out_lines.append("PASS AMR: вложенный 2x+4x точнее и дешевле; временная интерполяция снижает ошибку")
|
||
cm.write_out("amr_bench.txt", out_lines)
|
||
print("DONE (amr demo)")
|