添加#pragma omp parallel for后双曲偏微分方程并行程序结果异常求助

你的并行程序结果错误的核心原因是循环顺序导致的数据竞争，以及并行区域的错误选择。

问题根源 #

你当前的循环结构是外层遍历空间维度i，内层遍历时间维度j，并对i循环做并行化：

#pragma omp parallel for num_threads(9)
for(i = 1; i < XMAX - 1; i++){
    for(j = 1; j < TMAX - 1; j++){
        U[i][j + 1] = 2 * U[i][j] - U[i][j - 1] + sigma * (U[i + 1][j] - 2*U[i][j] + U[i - 1][j]) + (TAU*TAU * f(x[i],t[i]));
    }
}

这种情况下，不同线程处理不同的i，但计算U[i][j+1]时依赖U[i-1][j]或U[i+1][j]——这些值由其他线程负责计算，而线程的执行顺序是不确定的，会出现数据未就绪就被读取的竞争问题，导致计算结果错误。

另外注意公式中的t[i]应该是t[j]（时间维度对应j），否则时间项的取值会错误。

正确的并行实现方案 #

波动方程的显式格式要求按时间步推进：先完成所有空间点的第j+1层计算，再进入下一个时间步j。此时每个时间步内的空间点计算相互独立（仅依赖已确定的j和j-1层数据），可以安全并行化。

修改后的代码如下：

// 按时间步j循环，外层串行（时间步必须顺序执行）
for(j = 1; j < TMAX - 1; j++){
    // 并行化空间维度i的循环，每个线程处理部分空间点
    #pragma omp parallel for num_threads(9)
    for(i = 1; i < XMAX - 1; i++){
        // 修正t的索引为j，对应当前时间步
        U[i][j + 1] = 2 * U[i][j] - U[i][j - 1] + sigma * (U[i + 1][j] - 2*U[i][j] + U[i - 1][j]) + (TAU*TAU * f(x[i], t[j]));
    }
}

额外注意事项 #

• 确保初始条件（U[i][0]和U[i][1]）和边界条件（U[0][j]、U[XMAX-1][j]）在并行前已正确初始化，且并行区域内不会被修改。
• 检查sigma的取值是否满足稳定性条件（波动方程显式格式要求sigma = (c*TAU/H)^2 ≤ 1，其中c是波速，H是空间步长），稳定性问题也可能导致结果异常。

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