#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""空床流场诊断: 检查顶部是否有非物理进流 / 几何开口, 并采样 5 度倾斜效果。"""
import os
import sys

import numpy as np

sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
import fbed_simulation as fs

D = np.load(os.path.join(fs.OUT, "fbed_fields.npz"))
UX, UY, UZ, SOLID = D["ux"], D["uy"], D["uz"], D["solid"]
UX, UY, UZ = UX * fs.SCALE_U, UY * fs.SCALE_U, UZ * fs.SCALE_U
fluid = ~SOLID
NZ, NY, NX = SOLID.shape

print("=== 出口层 (z 顶部) 速度 ===")
for iz in (NZ - 1, NZ - 2):
    m = fluid[iz]
    u = UZ[iz][m]
    print(f" iz={iz} z={fs.ZC[iz]:.3f}  n={m.sum()}  uz[min,mean,max]="
          f"[{u.min():.2f}, {u.mean():.2f}, {u.max():.2f}] m/s  "
          f"inflow(<0)占比={100*np.mean(u<0):.1f}%")

print("\n=== 沸腾段顶部区域 (z>1.35) 是否有非物理入流 ===")
band = fluid & (fs.ZC[:, None, None] > 1.35)
uzb = UZ[band]
print(f" n={uzb.size}  uz mean={uzb.mean():.2f}  inflow<0 占比={100*np.mean(uzb<0):.2f}%")

print("\n=== 采样疑似'进风口'位置 (yz 切面左上亮斑) ===")
for (x, y, z) in [(0.15, -0.45, 1.2), (0.15, -0.35, 1.1), (2.25, -0.45, 1.2),
                  (0.15, -0.50, 0.9), (1.05, -0.45, 1.2)]:
    ix, iy, iz = fs._xi(x), fs._yi(y), fs._zi(z)
    if fluid[iz, iy, ix]:
        print(f" (x={x},y={y},z={z}): |U|={np.hypot(np.hypot(UX[iz,iy,ix],UY[iz,iy,ix]),UZ[iz,iy,ix]):.2f} "
              f"ux={UX[iz,iy,ix]:+.2f} uy={UY[iz,iy,ix]:+.2f} uz={UZ[iz,iy,ix]:+.2f}")
    else:
        print(f" (x={x},y={y},z={z}): SOLID")

print("\n=== 5 度倾斜效果: 布风板正上方薄层水平速度 (y 平均) ===")
for z in (0.04, 0.06, 0.10, 0.15):
    iz = fs._zi(z)
    m = fluid[iz] & (np.abs(fs.YC)[:, None] < fs.BED_W / 2)
    ux_m = UX[iz][m].mean()
    uz_m = UZ[iz][m].mean()
    print(f" z={z:.2f}  <ux>={ux_m:+.3f}  <uz>={uz_m:+.3f} m/s  "
          f"偏角={np.degrees(np.arctan2(ux_m, uz_m)):+.2f} deg (理论 +5 deg)")

print("\n=== 布风板各段流量 (按 x 分 4 段) ===")
izp = fs._zi(0.0)
cell_a = fs.DX ** 2
for k in range(4):
    x0, x1 = 0.6 * k, 0.6 * (k + 1)
    ix0, ix1 = fs._xi(x0), fs._xi(x1)
    sub = UZ[izp, :, ix0:ix1] * (fluid[izp, :, ix0:ix1])
    q = sub.sum() * cell_a * 3600.0
    print(f" 第{k+1}段 x[{x0:.1f},{x1:.1f}]: Q={q:.0f} m3/h")