向量化计算 - 《实验物理的大数据方法》

🌪️ 流体力学计算中向量化优化核心方法与注意事项

🔧 向量化替代循环的核心方法

网格坐标系统构建

# 正确矩阵索引（避免笛卡尔坐标错误）
i_grid, j_grid = np.meshgrid(
    np.arange(height), 
    np.arange(width),
    indexing='ij'  # 关键参数
)

原理：创建(height, width)形状的坐标矩阵，保持与物理场维度一致

物理量批量计算

# 向量化回溯位置计算
phys_x = i_grid + 0.5  # 网格中心物理坐标
phys_y = j_grid + 0.5
back_x = phys_x - velocity_x  # 所有点同时计算
back_y = phys_y - velocity_y

插值参数向量化

# 批量计算双线性插值参数
x0 = np.floor(back_x).astype(int)
y0 = np.floor(back_y).astype(int)
dx = back_x - x0
dy = back_y - y0
weights = np.stack([(1-dx)*(1-dy), dx*(1-dy), (1-dx)*dy, dx*dy], axis=2)

场更新向量操作

# 单指令更新整个场
advected_field = (
    weights[:,:,0] * field[x0, y0] +
    weights[:,:,1] * field[x1, y0] +
    weights[:,:,2] * field[x0, y1] +
    weights[:,:,3] * field[x1, y1]
)

⚠️ 关键注意事项

网格索引陷阱
- indexing='ij'：科学计算标准 (行,列)
- indexing='xy'：笛卡尔标准 (列,行)
  错误案例：
```
# 错误：默认xy导致维度反转
x, y = np.meshgrid(np.arange(width), np.arange(height))
```
边界一致性处理
- 物理坐标转换：网格点位于(i+0.5, j+0.5)
- 安全裁剪：np.clip(pos, 0, size-1.001)避免索引溢出
- 边界条件：速度场用factor=-1，压力场用factor=1

内存管理优化

# 避免中间数组副本
result = np.empty_like(input_array)
np.multiply(a, b, out=result)  # 原地操作

# 大数组分块处理
for i in range(0, height, block_size):
    block = array[i:i+block_size]
    # 处理分块

数值稳定性保障
- CFL条件：max(|u|)Δt/Δx < 1
- 浮点精度：使用np.float64减少累积误差
- 正则化处理：dyes = np.clip(dyes, 0, 1)