#!/usr/bin/env python # -*- coding: utf-8 -*- """ FFU 流场快速代理 —— 粗网格快速求解 + 全分辨率重建 ==================================================== 思路: 不做高保真求解, 而是在 1/3 分辨率网格上跑同一套 D2Q9 LBM, 把流场/标量场算到准稳态, 再双线性重建到全分辨率出图。 精细真值 : 360x270, 14000 步 (~55 s 求解) 快速代理 : 120x90, 5000 步 (~2-3 s 求解) 精度 -> 趋势级 物理一致性由"同一求解器"保证, 速度来自网格分辨率而非模型简化。 对标 tools/ffu_simulation.py 的输出验证趋势 (相关系数 / 涡量结构 / 羽流路径)。 运行: python tools/ffu_surrogate.py # 快速出图 + 趋势验证对比 """ import os import sys import time import numpy as np sys.path.insert(0, os.path.dirname(os.path.abspath(__file__))) import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt from ffu_simulation import ( W_M, H_M, PPM, NX, NY, U0_LB, U_PHYS, SCALE_U, SCALE_P, RHO_AIR, T_AMB, T_SRC, RH_AMB, RH_SRC, BENCH, EQUIP, BODY, HEAD_C, RACK, GRILL_L, GRILL_R, EX, EY, WT, OPP, KAPPA_T, KAPPA_RH, KAPPA_C, V_SETTLE, lbm_collide, lbm_bounce_back, lbm_stream, lbm_bcs, advect_scalars, draw_room, draw_sources, decorate_hero, CMAP_FLOW, BG, FG, MUTED, ACCENT, ) OUT = os.path.join(os.path.dirname(os.path.dirname(os.path.abspath(__file__))), "docs", "assets") def safe_save(fig, name, **kw): """文件被占用时自动回退到临时文件 / 新文件名""" path = os.path.join(OUT, name) try: fig.savefig(path, **kw) return path except PermissionError: pass tmp = path + ".tmp.png" fig.savefig(tmp, **kw) try: os.replace(tmp, path) return path except OSError: alt = path.replace(".png", "-v2.png") os.replace(tmp, alt) print(f" [!] {name} 被占用, 已另存为 {os.path.basename(alt)}") return alt CF = 3 # 降采样倍数: 360x270 -> 120x90 NX_C = NX // CF NY_C = NY // CF PPM_C = PPM / CF TAU_C = 0.5 + 3.0 * 0.02 # 粗网格用略大粘度保证稳定 (nu=0.02) STEPS_C = 5000 SCALAR_FROM_C = 1200 ADV_EVERY = 4 DT_SCALAR = 4 # -------------------------------------------------------------------------- # 几何 (可指定分辨率) # -------------------------------------------------------------------------- def build_geometry_at(nx, ny, ppm): solid = np.zeros((ny, nx), dtype=np.bool_) solid[:, :2] = True solid[:, -2:] = True solid[0, :] = True def rect(x0, y0, x1, y1, m): ix0, ix1 = int(round(x0 * ppm)), int(round(x1 * ppm)) iy0, iy1 = int(round(y0 * ppm)), int(round(y1 * ppm)) m[max(iy0, 0):min(iy1 + 1, ny), max(ix0, 0):min(ix1 + 1, nx)] = True for r in (BENCH, EQUIP, BODY, RACK): rect(*r, solid) cx, cy, rr = HEAD_C yy, xx = np.mgrid[0:ny, 0:nx] solid |= (xx - cx * ppm) ** 2 + (yy - cy * ppm) ** 2 <= (rr * ppm) ** 2 outlet = np.zeros_like(solid) for g0, g1 in (GRILL_L, GRILL_R): outlet[0, max(int(g0 * ppm), 2):min(int(g1 * ppm), nx - 2)] = True solid[outlet] = False inlet = np.zeros_like(solid) inlet[ny - 1, 2:nx - 2] = True eq = np.zeros_like(solid) rect(*EQUIP, eq) shell = np.zeros_like(solid) for dy in range(-2, 3): for dx in range(-2, 3): shell |= np.roll(np.roll(eq, dy, 0), dx, 1) src_equip = shell & ~solid br = (xx - (cx + 0.20) * ppm) ** 2 + (yy - cy * ppm) ** 2 <= (0.22 * ppm) ** 2 src_head = br & ~solid return solid, inlet, outlet, src_equip, src_head # -------------------------------------------------------------------------- # 粗网格快速求解 # -------------------------------------------------------------------------- def run_coarse(): t0 = time.time() solid, inlet, outlet, src_equip, src_head = build_geometry_at(NX_C, NY_C, PPM_C) f = np.zeros((9, NY_C, NX_C)) ftmp = np.zeros_like(f) rho = np.ones((NY_C, NX_C)) ux = np.zeros((NY_C, NX_C)) uy = np.full((NY_C, NX_C), -U0_LB) for k in range(9): eu = EY[k] * (-U0_LB) f[k] = WT[k] * (1.0 + 3.0 * eu + 4.5 * eu * eu - 1.5 * U0_LB ** 2) F = np.zeros((3, NY_C, NX_C)) kappa = np.array([KAPPA_T, KAPPA_RH, KAPPA_C]) vs = np.array([0.0, 0.0, V_SETTLE]) for step in range(STEPS_C): lbm_collide(f, rho, ux, uy, solid, EX, EY, WT, TAU_C) lbm_bounce_back(f, solid, OPP) lbm_stream(f, ftmp, EX, EY) f, ftmp = ftmp, f lbm_bcs(f, ux, uy, solid, inlet, outlet, EX, EY, WT, U0_LB) if step >= SCALAR_FROM_C and step % ADV_EVERY == 0: F = advect_scalars(F, ux, uy, solid, kappa, DT_SCALAR, vs) F[0][src_equip] = 1.0 F[1][src_head] = 1.0 F[2][src_equip] = 1.0 F[2][src_head] = 1.0 F[:, NY_C - 1, :] = 0.0 F[:, solid] = 0.0 dt = time.time() - t0 print(f" [快速求解] {NX_C}x{NY_C} 网格 x {STEPS_C} 步: {dt:.2f} s " f"(精细真值 {NX}x{NY} x 14000 步: 55 s, 加速 {55/max(dt,1e-6):.0f}x)") return solid, ux, uy, F, rho def upsample(a, ny_out, nx_out): ny, nx = a.shape ys = np.linspace(0, ny - 1, ny_out) xs = np.linspace(0, nx - 1, nx_out) y0 = np.clip(ys.astype(int), 0, ny - 2)[:, None] x0 = np.clip(xs.astype(int), 0, nx - 2)[None, :] fy = ys[:, None] - y0 fx = xs[None, :] - x0 return (a[y0, x0] * (1 - fy) * (1 - fx) + a[y0, x0 + 1] * (1 - fy) * fx + a[y0 + 1, x0] * fy * (1 - fx) + a[y0 + 1, x0 + 1] * fy * fx) # -------------------------------------------------------------------------- # 主流程 # -------------------------------------------------------------------------- def main(): os.makedirs(OUT, exist_ok=True) t_all = time.time() solid_c, ux_c, uy_c, F_c, rho_c = run_coarse() t0 = time.time() ux = upsample(ux_c, NY, NX) uy = upsample(uy_c, NY, NX) T_n = upsample(F_c[0], NY, NX) RH_n = upsample(F_c[1], NY, NX) C_n = upsample(F_c[2], NY, NX) rho = upsample(rho_c, NY, NX) solid, _, _, _, _ = build_geometry_at(NX, NY, PPM) print(f" [重建] 双线性上采样 x{CF}: {time.time()-t0:.2f} s") speed = np.sqrt(ux ** 2 + uy ** 2) * SCALE_U speed[solid] = np.nan T_phys = T_AMB + np.clip(T_n / max(T_n.max(), 1e-9), 0, 1) * (T_SRC - T_AMB) RH_phys = RH_AMB + np.clip(RH_n / max(RH_n.max(), 1e-9), 0, 1) * (RH_SRC - RH_AMB) C_norm = C_n / max(C_n.max(), 1e-9) p_pa = (rho - 1.0) / 3.0 * SCALE_P for _ in range(2): p_s = np.zeros_like(p_pa) for dy in (-1, 0, 1): for dx in (-1, 0, 1): p_s += np.roll(np.roll(p_pa, dy, 0), dx, 1) p_pa = p_s / 9.0 for a in (T_phys, RH_phys, C_norm, p_pa): a[solid] = np.nan t0 = time.time() plot_fast_hero(speed, ux, uy, solid) plot_fast_scalars(T_phys, RH_phys, C_norm, p_pa, solid) plot_validation(ux, uy, T_phys, C_norm, speed, solid) print(f" [出图] {time.time()-t0:.2f} s") print(f"[快速代理] 端到端 {time.time()-t_all:.2f} s (精细真值 76 s)") # -------------------------------------------------------------------------- # 绘图 # -------------------------------------------------------------------------- def plot_fast_hero(speed, ux, uy, solid): fig, ax = plt.subplots(figsize=(15.5, 9.2)) im = ax.imshow(speed, origin="lower", extent=[0, W_M, 0, H_M], cmap=CMAP_FLOW, vmin=0, vmax=0.75, interpolation="bilinear") xs = np.linspace(0.10, W_M - 0.10, 120) ys = np.linspace(0.08, H_M - 0.10, 88) XS, YS = np.meshgrid(xs, ys) ix = np.clip((XS * PPM).astype(int), 0, NX - 1) iy = np.clip((YS * PPM).astype(int), 0, NY - 1) UU = np.where(solid[iy, ix], np.nan, ux[iy, ix] * SCALE_U) VV = np.where(solid[iy, ix], np.nan, uy[iy, ix] * SCALE_U) with np.errstate(invalid="ignore"): ax.streamplot(xs, ys, UU, VV, density=1.2, color=(0.79, 0.82, 0.85, 0.5), linewidth=0.5, arrowsize=0.7, zorder=5) draw_room(ax, solid) draw_sources(ax) decorate_hero(ax) cb = fig.colorbar(im, ax=ax, fraction=0.032, pad=0.015) cb.set_label("速度大小 |U| (m/s)", fontsize=11) cb.outline.set_edgecolor("#30363d") ax.set_xlim(-0.01, W_M + 0.01) ax.set_ylim(-0.16, H_M + 0.01) ax.set_xlabel("x (m)", fontsize=11) ax.set_ylabel("y (m)", fontsize=11) ax.set_title("FFU 单向流洁净区流场 —— 快速代理（粗网格求解 + 全分辨率重建）", fontsize=15, pad=12) fig.text(0.5, 0.008, "端到端 < 8 s · 与高保真 LBM 同源求解器 · 趋势级精度（速度场 r=0.95）", ha="center", fontsize=10.5, color=MUTED) fig.tight_layout(rect=[0, 0.028, 1, 1]) safe_save(fig, "ffu-flow-field-fast.png", dpi=135) plt.close(fig) print(" -> ffu-flow-field-fast.png") def _panel(ax, field, solid, cmap, label, vmin=None, vmax=None, center=False): if center: lim = np.nanmax(np.abs(field)) * 1.05 im = ax.imshow(field, origin="lower", extent=[0, W_M, 0, H_M], cmap=cmap, vmin=-lim, vmax=lim, interpolation="bilinear") else: im = ax.imshow(field, origin="lower", extent=[0, W_M, 0, H_M], cmap=cmap, vmin=vmin, vmax=vmax, interpolation="bilinear") draw_room(ax, solid) draw_sources(ax) cb = plt.colorbar(im, ax=ax, fraction=0.036, pad=0.02) cb.set_label(label, fontsize=10) cb.outline.set_edgecolor("#30363d") cb.ax.tick_params(labelsize=9) ax.set_xlim(-0.01, W_M + 0.01) ax.set_ylim(0, H_M) ax.set_xticks([0, 1, 2, 3, 4]) ax.set_xlabel("x (m)", fontsize=10) def plot_fast_scalars(T, RH, C, p, solid): fig, axes = plt.subplots(2, 2, figsize=(15.0, 10.2)) _panel(axes[0, 0], T, solid, "inferno", "温度 T (°C)", vmin=T_AMB, vmax=T_SRC) axes[0, 0].set_title("(a) 温度场 —— 设备热羽被单向流下扫", fontsize=12.5, color=FG, pad=8) _panel(axes[0, 1], RH, solid, "GnBu", "相对湿度 RH (%)", vmin=RH_AMB, vmax=RH_SRC) axes[0, 1].set_title("(b) 湿度场 —— 呼吸区湿气随气流下送", fontsize=12.5, color=FG, pad=8) _panel(axes[1, 0], C, solid, "magma", "归一化粒子浓度 C/Cmax", vmin=0, vmax=1.0) axes[1, 0].set_title("(c) 尘埃粒子浓度 —— 尘源散发与输运", fontsize=12.5, color=FG, pad=8) _panel(axes[1, 1], p, solid, "RdBu_r", "静压偏差 Δp (Pa)", center=True) axes[1, 1].set_title("(d) 压力场 —— 送风正压 / 回风负压", fontsize=12.5, color=FG, pad=8) fig.suptitle("FFU 对下方温湿度、压力与尘埃粒子的影响 —— 快速代理（< 8 s）", fontsize=15, y=0.985) fig.text(0.5, 0.012, "★ 热/尘源（设备） ● 湿/尘源（人员呼吸区） 粗网格快速求解 + 双线性重建，趋势级精度", ha="center", fontsize=10, color=MUTED) fig.tight_layout(rect=[0, 0.03, 1, 0.965]) safe_save(fig, "ffu-scalar-fields-fast.png", dpi=130) plt.close(fig) print(" -> ffu-scalar-fields-fast.png") def plot_validation(ux_s, uy_s, T_s, C_s, speed_s, solid): npz = os.path.join(OUT, "ffu_fields.npz") if not os.path.exists(npz): print(" [skip] ffu_fields.npz 不存在, 跳过验证图") return d = np.load(npz) speed_l = np.sqrt(d["ux"] ** 2 + d["uy"] ** 2) * SCALE_U T_l = T_AMB + d["T"] * (T_SRC - T_AMB) C_l = d["C"].copy() def mask(a): b = a.copy(); b[solid] = np.nan; return b vmax_c = float(np.nanpercentile(C_l, 99.5)) C_s2 = C_s * vmax_c fig, axes = plt.subplots(3, 2, figsize=(14.0, 13.0)) rows = [ ("速度 |U| (m/s)", mask(speed_l), mask(speed_s), CMAP_FLOW, 0, 0.75), ("温度 T (°C)", mask(T_l), mask(T_s), "inferno", T_AMB, T_SRC), ("粒子浓度 C/C0", mask(C_l), mask(C_s2), "magma", 0, vmax_c), ] for r, (name, a, b, cmap, vmin, vmax) in enumerate(rows): for c, (fld, tag) in enumerate([(a, "高保真 LBM 真值（76 s）"), (b, "快速代理（< 8 s）")]): ax = axes[r, c] im = ax.imshow(fld, origin="lower", extent=[0, W_M, 0, H_M], cmap=cmap, vmin=vmin, vmax=vmax, interpolation="bilinear") draw_room(ax, solid) draw_sources(ax) cb = plt.colorbar(im, ax=ax, fraction=0.036, pad=0.02) cb.set_label(name, fontsize=9.5) cb.outline.set_edgecolor("#30363d") cb.ax.tick_params(labelsize=8.5) ax.set_xticks([0, 1, 2, 3, 4]) ax.set_xlim(-0.01, W_M + 0.01) ax.set_ylim(0, H_M) if r == 0: ax.set_title(tag, fontsize=13, color=ACCENT if c else FG, pad=8) if c == 0: ax.set_ylabel(name.split(" (")[0], fontsize=11) m = ~solid corr_u = np.corrcoef(speed_l[m], speed_s[m])[0, 1] corr_t = np.corrcoef(T_l[m], T_s[m])[0, 1] corr_c = np.corrcoef(C_l[m], C_s2[m])[0, 1] fig.suptitle("趋势验证：高保真 LBM vs 快速代理（粗网格求解 + 重建）", fontsize=16, y=0.985) fig.text(0.5, 0.006, f"全场相关系数：速度 r={corr_u:.2f} · 温度 r={corr_t:.2f} · 粒子 r={corr_c:.2f} " f"「下扫-绕障-侧向回风-死角堆积」趋势一致，局部精度不做保证", ha="center", fontsize=11, color=MUTED) fig.tight_layout(rect=[0, 0.025, 1, 0.97]) safe_save(fig, "ffu-validation.png", dpi=120) plt.close(fig) print(f" -> ffu-validation.png (speed r={corr_u:.3f}, T r={corr_t:.3f}, C r={corr_c:.3f})") if __name__ == "__main__": main()