气动涡格法求解器C移植的Off-by-one错误排查求助

问题背景 #

将已验证正确的Python气动涡格法求解器移植到C语言后，程序仅在特殊翼型（如NACA 9999）下出现异常：常规翼型（如NACA 2412）的计算结果与Python完全一致，但特殊翼型的升阻比（L/D）错误。

定位过程 #

1. 排查发现问题源于build_system函数生成的A矩阵：常规翼型下矩阵与Python版本一致，特殊翼型下最后一行的最后一个元素存在Off-by-one错误，其余元素及b向量均正常，该错误最终导致求解结果异常。
2. 进一步追踪到错误出现在s.A[NUM_PANELS][NUM_PANELS]的累加逻辑中，gv_l.b的值被取反。
3. 根源锁定为transform_to_local函数的输出：C代码中的l_l.b（对应Python的y_local）符号与Python版本相反。输入参数仅存在1e-19级别的浮点数差异（属于正常精度范畴），dx、dy及r.a（对应Python的x_local）均一致，仅r.b符号相反。
4. 已确认翼型点、面板参数及面板顺序完全一致。

核心原因分析 #

• 排除double与numpyfloat64的精度差异：两者均为64位双精度浮点数，精度特性完全一致，1e-19的差异是不同语言运算顺序导致的正常舍入误差，不会直接引发符号翻转。
• 问题本质是局部坐标系变换的符号逻辑不一致：
  • 涡格法中局部坐标系通常以面板切线为x轴，法线为y轴，法线方向的定义（如指向翼型外侧/内侧、向上/向下）直接决定y_local的符号。
  • 极端厚弦比的翼型（如NACA 9999）面板法线计算更容易触发临界情况：接近0的数值因微小舍入误差，触发代码中依赖符号的分支逻辑，导致符号翻转。
• 可能的代码细节问题：
  • Python与C代码中向量叉乘的顺序不同，导致法线方向计算结果符号相反；
  • 局部坐标系变换公式中对y轴方向的定义不一致；
  • 代码中存在直接判断数值符号的逻辑（如if (val > 0)），未考虑微小浮点数误差的影响。

解决方案建议 #

1. 对齐局部坐标系定义
   • 对比Python与C代码中面板法线方向的计算逻辑：检查向量叉乘的顺序（例如(x2-x1, y2-y1)叉乘(0,1)还是(1,0)），确保两者的法线方向定义完全一致（如均指向翼型外侧）。
2. 修复transform_to_local的符号错误
   • 若确认C代码的r.b符号与Python相反，直接在函数中对结果取反，或调整坐标系变换的公式/矩阵。例如：

     // 假设需要翻转y_local符号
     r.b = -r.b;
     

3. 添加临界值处理逻辑
   • 针对接近0的浮点数，避免直接依赖符号判断，引入epsilon（如1e-12）过滤微小误差：

     const double EPS = 1e-12;
     if (fabs(r.b) < EPS) {
         r.b = 0.0;
     }
     // 再根据正确的坐标系定义确认符号
     

4. 针对性验证极端场景
   • 手动计算NACA 9999最后一个面板的transform_to_local输出，对比Python与C的结果，明确符号差异的触发条件；
   • 单独验证build_system中对gv_l.b的使用逻辑，确保符号处理与Python完全匹配。

验证步骤 #

1. 修改后重新生成A矩阵，确认最后一行最后一个元素与Python版本完全一致；
2. 计算NACA 9999的升阻比，验证结果与Python匹配；
3. 测试其他极端翼型（如超高厚弦比、不对称翼型），确保问题彻底解决。

内容的提问来源于stack exchange，提问作者Pazzel