Frank14f/CelerisLab
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
main
CelerisLab/docs/audit/kernel_layer_audit.md

11 KiB
Raw Permalink Blame History

Sub-agent A: CUDA Kernel & Fix Verification

1. 旧修正确认

结论: [已确认] Curved Bouzidi 时序 — aux_kernels.cu:26-65 + stepper.py:71 + one_step_double.cu:20

理由: stepper.py:71 中 _launch_curved() 在 OneStep 之前调用。aux_kernels.cu 的 CurvedBoundaryKernel 写入 f.ddf_gpu（当前 buffer），其后 OneStep 从同一 buffer 做 stream_pull_load — 同一时间步内完成"前写入→后拉取"。one_step_double.cu:20 中的 apply_boundary_pull 对 is_curved(fl) 直接 return，确保流体节点不会二次覆盖。调用顺序确认无误。

结论: [已确认] q>=0.5 分支时间层 — curved_boundary.cuh:73-75

理由: q >= 0.5 分支读取 load_ddf(fi, index_f(k_f, dir_opp))，其中 fi 是当前步骤的 DDF buffer。旧代码读的是前一时间步的缓冲区（fi_in），现在读的是当前的 fi，即上一时间步碰撞后的 post-collision 数据。Bouzidi 两分支算法都要求同一时间层的 post-collision 数据 [Bou01]，当前写法满足此要求。

结论: [已确认] Moving wall 修正按 q 分支 — curved_boundary.cuh:30-43

理由: bouzidi_linear_moving_correction 函数三路分支：

  • fallback_class != BOUZIDI → 2 * alpha_ci_dot_uw（半格加移动修正）
  • q < 0.5 → 2 * alpha_ci_dot_uw（Bou01 公式）
  • q >= 0.5 → alpha_ci_dot_uw / q（Bou01 公式） 各项系数与文献一致，alpha_ci_dot_uw = 3 * w_i * (c_i · u_w)。

结论: [已确认] Forcing 预因子统一 — collision_srt.cuh:15, collision_trt.cuh:35, collision_mrt.cuh:41/116

理由: 三个文件的类内路径都使用了 c_tau = 1.0f - 0.5f * omega 并将 c_tau * Fin[i] 增量加到碰撞输出中。collide_dispatch()（helpers.cuh:33-37/46-50）在 d_params.fx/fy/fz非零时调用 compute_guo_forcing() 生成 Fin，然后传入对应碰撞函数。三者完全一致。

结论: [已确认] TRT outlet NEQ 全分布重构 — pressure_neq.cuh:45-51(D2Q9):89-95(D3Q19)

理由: #if COLLISION_MODEL == 0 || COLLISION_MODEL == 1 把 TRT 纳入全分布 damped NEQ 分支（与 SRT 同级），对所有 NQ 方向做 f[i] = feq_tar[i] + beta * fneq，其中 fneq = f_neb[i] - feq_neb[i]，beta = OUTLET_SRT_NEQ_DAMP。已不再是"少量未知方向"路径。

结论: [已确认] add_vortex() 除以 rho — initializers.py:55-58

理由: ux_old = sum(f[i] * cx[i]) / rho_safe，uy_old = sum(f[i] * cy[i]) / rho_safe。旧 bug 是动量除以 rho 这一步缺失，现在存在 rho_safe 分母。

结论: [已确认] Init flag overlay 顺序 — simulation.py:139-159 + init_flow.cu:48-58 + body/manager.py:115-127

理由: 顺序为 build_channel_flags()（干净通道）→ build_flags()（叠加物体）→ upload_flags()（上传GPU）→ stepper.initialize()（运行 init kernel 保持 obstacle flag）。init_flow.cu:finalize_domain_flag 对 is_obstacle(fl) 直接返回原 flag。sync_to_gpu(rebuild_flags=False) 在初始化时不重复构建 flags。_rest_nonfluid() 的二次重置已移除。

2. 待验证项复查

结论: [已确认] MRT D2Q9 方向序与符号 — collision_mrt.cuh:46-54,79-92

理由: 经子代 agent 验证：

  • m[3] = f1−f2 + f5−f6 + f7−f8 = ρ·ux（与 macro.cuh:36 的 compute_rho_u() 完全一致）
  • m[5] = f3−f4 + f5−f6 − f7 + f8 = ρ·uy（一致）
  • m[7] = f1+f2−f3−f4 = ρ·(ux²−uy²)（一致）
  • 逆变换系数 g[0] += (dm0 − dm1 + dm2)/9 等均符合 M·M⁻¹ = I
  • meq[4] = −ρ·ux、meq[6] = −ρ·uy 符号正确（正交基结构）
  • 无符号错误。MRT D2Q9 的 paired 方向排列下的矩变换自洽。

结论: [已确认] Esopull 邻壁行处理 — one_step_esopull.cu vs one_step_double.cu

理由: 逐行对比 one_step_esopull.cu:76-151 与 one_step_double.cu:65-144。Double-buffer 在 line 88-94 有显式 is_fluid(fl) && (y==1 || y==NY-2) 块调用 apply_wall_bb_y_pull，而 Esopull 没有对应分支。但这不构成 bug — Esopull 交替读写模式（load_f_esopull / store_f_esopull，esopull_single_buffer.cuh:33-77）天然避免了从 y=0 壁面节点直接拉取垃圾数据：

  • 偶数步：f[3] 读取本地 fi[n, 4]（前一步该节点的 -y 方向），f[4] 读取 fi[j[3], 3]（y=2 的 +y 方向）
  • 奇数步同理交错 两个时间步的综合效果等价于半格 BB 的反射，无需显式修正。壁面节点 (y=0) 经 apply_boundary_esopull → bounce_back_swap 处理，不会积累垃圾。

结论: [无法确认: 需运行验证] Force 提取与符号约定 — curved_boundary.cuh:90-97 + obs.cuh + manager.py:328-338

理由: 链路追踪如下：

  1. curved_boundary.cuh:91-97: fx = c_x * (f_toward + f_reflected) 累加到 obs[obs_force_index(body_id, 0)]
  2. obs.cuh:11: obs_force_index(id_obj, d) = OBS_FORCE0_FLOATS + id_obj * DIM + d（起点=0）
  3. manager.py:read_force(): 返回 self.obs_pinned[i0:i0 + d]，无符号反转
  • 存储的是流体动量交换量 (fluid momentum exchange)，不是物体受力的直接值（牛顿第三定律要求 F_body = -F_fluid）。如果外部 Cd/Cl 归一化时把 obs 直读值当做物体受力，符号会反转。需用已知算例（如静止圆柱的阻力系数）确认当前实践是否在外部做了隐含的取反。

3. 架构缺陷

结论: [未改] Curved boundary 仍假设圆形/球形 — curved_boundary.cuh + body/objects.py:110-291

理由: Cylinder.get_curved_list() 对每个 cut link 调用 find_circle_intersection / find_sphere_ray_segment，完全依赖于圆/球几何参数（center + radius）。没有通用多边形/三角网格接口。旧审计标记的 [待重构] 未处理。

结论: [未改] 3D 旋转为 z 轴占位 — aux_kernels.cu:58-63

理由: CurvedBoundaryKernel 的 D3Q19 分支中 Ww = 0.0f，且 Uw = -omega * ry; Vw = omega * rx 只包含 xy 平面转动。任何具有 z 分量的刚体旋转都不会产生正确的壁面速度。

结论: [已备注] Bouzidi-TRT 不相容注释 — curved_boundary.cuh:17-21

理由: 注释明确说明 "plain linear Bouzidi interpolation ... is not a TRT-parametrized curved-wall family"，位置精准且语意清晰。

4. 新发现

4a. _no_force 碰撞变体为死代码 — collision_srt.cuh:25, collision_trt.cuh:57, collision_mrt.cuh:95

问题: collide_srt_no_force、collide_trt_no_force、collide_mrt_no_force 三个函数在编译后的任何内核路径中均不被调用。collide_dispatch (helpers.cuh:15-69) 始终通过带 Fin 数组的普通变体执行碰撞，当外力为零时通过 zero_forcing(Fin) 将 Fin 全清零。死代码约 40 行。 影响: 无功能影响，增加维护成本。zero_forcing 仅用于栈初始化，并非无用调用。

4b. BC_MOVING / BC_PERIODIC 有定义无处理 — flags.cuh + 全内核搜

问题: FLAG_BC_MOVING (0x0060) 和 FLAG_BC_PERIODIC (0x0050) 在 flags.cuh 有定义，Python 侧 descriptors.py 也有对应常量。但 CUDA 内核中没有针对 is_moving() 或 is_periodic() 的分支处理：

  • one_step_double.cu:apply_boundary_pull 只有 is_curved/is_inlet/is_outlet/BBS 分支
  • 任何单元格若被标记 BC_MOVING（非 curved 路径）将落入 bounce_back_swap() 分支，得到 std half-way BB 而非移动壁面速度修正 影响: 只有当 obstacle 使用 BC_CURVED 配合 cl_body_id 进入 curved boundary kernel 才能获得正确移动壁面速度。BC_MOVING 直接标记在通道壁面或其他固体节点上无效果。

4c. Outlet NEQ 的 OUTLET_MODE 嵌套逻辑可读性隐患 — pressure_neq.cuh:25-62

问题:

#if OUTLET_MODE == 1        → 纯拷贝未知方向
#else                       → 普通路径
  #if COLLISION_MODEL==0||1 → 全分布 NEQ（含 TRT）
  #elif OUTLET_MODE == 2    → 混合模式
  #else                     → 少量未知方向重构（默认 OUTLET_MODE=0 路径）
  #endif
#endif

当 OUTLET_MODE=0、COLLISION_MODEL=2（MRT）时，代码进入最内层 #else 分支（少量未知方向重构），而非全分布 NEQ。这与注释宣称的"SRT 和 TRT 使用全分布 NEQ"一致（MRT 未承诺），但 MRT outlet 路径与 SRT/TRT 行为不同，可能导致 MRT 结果系统性偏差。 影响: MRT outlet 行为与 SRT/TRT 不一致，应至少加注释说明此差异，或行为对齐。

4d. collide_inlet_ghost 对 y=0/NY-1 的过滤 — one_step_double.cu:102-103, one_step_esopull.cu:106-107

问题: collide_inlet_ghost = is_inlet(fl) && interior_y && inlet_scheme_uses_post_collision_ghost()，其中 interior_y = (y>0 && y<NY-1)。y=0 和 y=NY-1 的 inlet 节点不会触发 ghost 碰撞。但 inlet_scheme_uses_post_collision_ghost() 只在 INLET_SCHEME==0（Zou-He）时返回 true，而 x=0 inlet 节点原本处于 y=0 或 y=NY-1 时，已被 build_channel_flags 覆盖为 SOLID|BC_WALL（角点优先），因此这些节点本身不是 inlet，过滤是安全的。 影响: 当前无实际影响（角点 wall 覆盖 inlet），但代码隐含的逻辑依赖比较脆弱。如果未来修改建标记序，可能在 y=0/ymax inlet 节点产生未碰撞 ghost 状态。

4e. 传感器归一化在 manager 层正确实现 — aux_kernels.cu:67-99 + manager.py:348-365

结论: SensorKernel 做逐格点求和，ObjectManager.read_sensor(normalize=True) 除以 sensor_cell_counts[body_id]。已实现，正确。

4f. inlet_target_u 对 y=0 和 y=NY-1 的保护 — inlet/common.cuh:34

影响: y_clamped = fminf(NY-2, fmaxf(1.0, y)) 使 y=0 节点获得 y=1 的 inlet 速度。对正确定义的 SOLID|BC_WALL 角点无实际影响，属于保护性逻辑。

4g. compute_omega_minus 在 ω⁺ = 2.0 时分母为零 — collision_trt.cuh:24-26

问题: 1.0f / omega_plus - 0.5f 当 omega_plus = 2.0 时为 0.5 - 0.5 = 0，导致除零。实际路径中 omega_col 在 collide_dispatch 中被钳位（helpers.cuh:56），假设 OMEGA_COLLISION_MAX < 2.0 则为安全。 影响: 低风险（依赖外部宏正确配置）。

4h. west_velocity_rho_closure 在 u_target ≥ 1.0 时除零 — inlet/common.cuh:47-48

问题: rho = sum(...) / (1.0f - ux_target)，当 ux_target >= 1.0 时除数为 0 或负数。物理上入口马赫数小于 1 可避免，但无运行时保护。 影响: 低风险（物理约束），但崩溃行为不如显式断言清晰。

总结

类别 数量 关键项
旧修正确认 7/7 全部确认正确
待验复查 3 2 确认, 1 需运行时验证
架构缺陷 3 圆形硬编码、3D占位、TRT注释齐全
新发现 8 4a死代码、4b未处理BC、4c MRT outlet差异、4d脆弱逻辑依赖、4e已实现、4f无害保护、4g除零边界、4h无保护输入

最高优先级: 4b (BC_MOVING 无处理路径)、4c (MRT outlet 路径与 SRT/TRT 不一致)、force 符号约定需运行验证。