# Демо-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)")