使用OMP并行化洪水模拟器循环结果不一致问题求助

在自主开发的洪水模拟器中，尝试用OpenMP并行化核心循环以提升性能，但并行后结果与串行执行不一致，以下是串行代码、并行尝试代码及问题分析和解决建议。

串行核心代码 #

#define accessMat( arr, exp1, exp2 ) arr[ (int)(exp1) * columns + (int)(exp2) ] 
#define PRECISION       1000000
#define FIXED(a)    ( (int)((a) * PRECISION) )

max_spillage_iter = 0.0;
for ( row_pos=0; row_pos<rows; row_pos++ ) {
    for ( col_pos=0; col_pos<columns; col_pos++ ) {
        // If the cell has spillage
        if( accessMat( spillage_flag, row_pos, col_pos ) == 1 ) {

            // Eliminate the spillage from the origin cell
            accessMat( water_level, row_pos, col_pos ) -= FIXED( accessMat( spillage_level, row_pos, col_pos ) / SPILLAGE_FACTOR );

            // Compute termination condition: Maximum cell spillage during the iteration
            if ( accessMat( spillage_level, row_pos, col_pos ) / SPILLAGE_FACTOR > max_spillage_iter ) {
                max_spillage_iter = accessMat( spillage_level, row_pos, col_pos ) / SPILLAGE_FACTOR;
            }
            // Statistics: Record maximum cell spillage during the scenario and its time
            if ( accessMat( spillage_level, row_pos, col_pos ) / SPILLAGE_FACTOR > max_spillage_scenario ) {
                max_spillage_scenario = accessMat( spillage_level, row_pos, col_pos ) / SPILLAGE_FACTOR;
                max_spillage_minute = minute;
            }
        }

        // Accumulate spillage from neighbors
        for ( cell_pos= 0; cell_pos < CONTIGUOUS_CELLS; cell_pos++ ) {
            int depths = CONTIGUOUS_CELLS;
            accessMat( water_level, row_pos, col_pos ) += FIXED( accessMat3D( spillage_from_neigh, row_pos, col_pos, cell_pos ) / SPILLAGE_FACTOR );
        }
    }
}

并行尝试代码 #

#pragma omp parallel for collapse(2),private(row_pos, col_pos, cell_pos),shared(water_level, minute),reduction(max:local_max_spillage_iter, local_max_spillage_scenario)
for ( row_pos=0; row_pos<rows; row_pos++ ) {
    for ( col_pos=0; col_pos<columns; col_pos++ ) {
        // If the cell has spillage
        if( accessMat( spillage_flag, row_pos, col_pos ) == 1 ) {

            // Eliminate the spillage from the origin cell
            accessMat( water_level, row_pos, col_pos ) -= FIXED( accessMat( spillage_level, row_pos, col_pos ) / SPILLAGE_FACTOR );

            // Compute termination condition: Maximum cell spillage during the iteration
            if ( accessMat( spillage_level, row_pos, col_pos ) / SPILLAGE_FACTOR > max_spillage_iter ) {
                max_spillage_iter = accessMat( spillage_level, row_pos, col_pos ) / SPILLAGE_FACTOR;
            }
            // Statistics: Record maximum cell spillage during the scenario and its time
            if ( accessMat( spillage_level, row_pos, col_pos ) / SPILLAGE_FACTOR > max_spillage_scenario ) {
                max_spillage_scenario = accessMat( spillage_level, row_pos, col_pos ) / SPILLAGE_FACTOR;
                max_spillage_minute = minute;
            }
        }

        // Accumulate spillage from neighbors
        for ( cell_pos= 0; cell_pos < CONTIGUOUS_CELLS; cell_pos++ ) {
            int depths = CONTIGUOUS_CELLS;
            #pragma omp atomic
            accessMat( water_level, row_pos, col_pos ) += FIXED( accessMat3D( spillage_from_neigh, row_pos, col_pos, cell_pos ) / SPILLAGE_FACTOR );
        }
    }
}

问题分析与解决建议 #

1. 全局最大值变量的竞争条件 #

并行代码中，max_spillage_iter和max_spillage_scenario是全局共享变量，直接在多线程中修改会导致竞争条件：多个线程同时读写这些变量时，会出现值覆盖或读取不一致的情况。虽然声明了reduction(max:local_max_spillage_iter, local_max_spillage_scenario)，但代码中并未使用这些局部变量，而是直接操作全局变量，导致reduction机制失效。

修正方案：
每个线程使用局部变量存储当前线程内的最大值，最后通过reduction合并到全局变量：

#pragma omp parallel for collapse(2) \
private(row_pos, col_pos, cell_pos, local_spill) \
shared(water_level, minute, spillage_flag, spillage_level, spillage_from_neigh) \
reduction(max:max_spillage_iter, max_spillage_scenario)
for ( row_pos=0; row_pos<rows; row_pos++ ) {
    for ( col_pos=0; col_pos<columns; col_pos++ ) {
        if( accessMat( spillage_flag, row_pos, col_pos ) == 1 ) {
            local_spill = accessMat( spillage_level, row_pos, col_pos ) / SPILLAGE_FACTOR;
            // 消除自身溢水
            accessMat( water_level, row_pos, col_pos ) -= FIXED(local_spill);
            
            // 更新线程局部的迭代最大值（reduction会自动合并到全局max_spillage_iter）
            if (local_spill > max_spillage_iter) {
                max_spillage_iter = local_spill;
            }
            
            // 场景最大值需要额外处理，因为关联了时间戳
            if (local_spill > max_spillage_scenario) {
                #pragma omp critical
                {
                    // 临界区内再次检查，避免线程间竞争覆盖
                    if (local_spill > max_spillage_scenario) {
                        max_spillage_scenario = local_spill;
                        max_spillage_minute = minute;
                    }
                }
            }
        }

        // 累加邻居溢水：每个cell由唯一线程处理，无需atomic
        for ( cell_pos= 0; cell_pos < CONTIGUOUS_CELLS; cell_pos++ ) {
            accessMat( water_level, row_pos, col_pos ) += FIXED( accessMat3D( spillage_from_neigh, row_pos, col_pos, cell_pos ) / SPILLAGE_FACTOR );
        }
    }
}

2. 不必要的atomic操作导致的异常 #

并行代码中对water_level的累加操作添加了#pragma omp atomic，但在collapse(2)的并行循环中，每个(row_pos, col_pos)单元格由唯一线程负责处理，不存在多线程同时修改同一个单元格的情况，atomic操作完全多余，反而会引入不必要的同步开销，甚至可能因为原子操作的内存模型导致结果偏差。

修正方案：直接移除#pragma omp atomic，保留原有的累加逻辑即可。

3. 场景最大值时间戳的同步问题 #

max_spillage_minute与max_spillage_scenario强关联，当更新全局场景最大值时，必须保证时间戳与最大值的一致性。多线程环境下，若多个线程同时尝试更新，会导致时间戳与最大值不匹配的问题。

修正方案：使用#pragma omp critical临界区包裹场景最大值和时间戳的更新逻辑，确保同一时间只有一个线程能修改这两个变量，并且在临界区内再次检查最大值，避免线程竞争导致的错误覆盖。

4. 变量作用域与冗余定义 #

并行代码中int depths = CONTIGUOUS_CELLS;是冗余定义，未被使用，建议删除以优化代码。

关键优化总结 #

• 利用OpenMP的reduction机制处理迭代最大值，避免全局变量竞争；
• 用临界区保护场景最大值和时间戳的更新，保证数据一致性；
• 移除不必要的atomic操作，减少同步开销；
• 提取重复计算的local_spill变量，避免重复计算，提升性能。

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