fortran - 与串行和附加并行运行相比，并行仿真在一些时间步后给出不同的结果

Question

我正在尝试使用 OpenMP 并行运行涡流模拟代码。这些类似于粒子模拟，其中在每个时间步长，必须根据其速度计算下一个时间步长的涡流位置，该速度由当前时间步长的所有其他涡流的位置确定。一旦离开域，涡流就会被删除。我比较并行版本代码和串行版本代码在每个时间步的涡流数量，并多次运行每个版本。

对于串行版本，涡流计数在每个时间步都完全匹配。对于并行情况，在几十个时间步长内，所有运行都与串行情况匹配，之后，每个并行运行都显示出差异，但与串行情况保持在 7-10% 的误差范围内(如中所示下面的结果链接)。我知道这可能是由于并行情况下由于不同线程之间的分布导致计算步骤顺序不同而导致的舍入误差，但误差真的应该高达 10% 吗？

我只在并行 do 构造中使用了 reduction 子句。整个代码中唯一的并行区域是函数 vblob() 内，它位于模块内，我从主代码中调用它。 vblob()内的所有函数调用都是ixi()，fxi()在这个模块之外。

function vblob(blobs,xj,gj)
    complex(8), intent(in) :: blobs(:,:), xj
    complex(8) :: delxi, delxic, di, gvic, xi
    real(8), intent(in) :: gj
    real(8) :: vblob(2)
    integer :: p

    gvic = 0.0; delxi = 0.0; delxic = 0.0; di = 0.0; xi = 0.0
    !$omp parallel do private(xi,delxic,delxi,di) shared(xj) reduction(+:gvic)
    do p = 1, size(blobs,1)
      xi = ixi(blobs(p,1))
      delxic = xj-conjg(xi)
      delxi = xj-xi
      di = del*fxi(xi)
      gvic = gvic + real(blobs(p,2))*1/delxic
      if (abs(delxi) .gt. 1.E-4) then
        gvic = gvic +  (-1)*real(blobs(p,2))*1/delxi
      end if
    end do
    !$omp end parallel do
    gvic = j*gvic*fxi(xj)/(2*pi)
    vblob(1) = real(gvic)
    vblob(2) = -imag(gvic)

  end function vblob

如果我构建并行代码的方式是错误的，那么错误应该从前几个时间步开始就出现了，对吧？

(在这个result中可以看出，'blobs'和'sheets'只是涡旋元素的类型，蓝线是总元素。P和S分别代表Parallel和Serial， R代表runs，实心的plot marker是串行代码，空心的是并行代码的三轮)

编辑:当我将变量的数值精度更改为 real(4) 时，结果中的分歧发生在比上面的 real(8) 情况更早的时间步长。所以这显然是一个舍入错误问题。

TLDR:我想与任何其他在一定时间步长范围内看到此类结果的人一起澄清这一点，其中并行代码在前几个时间步长匹配然后发散？

Answer 1

您的代码基本上总结了 gvic 中的很多术语。浮点运算不是关联的，即 (a+b)+c 由于四舍五入而与 a+(b+c) 不同。此外，根据条款上的值和符号，每次操作可能会严重损失精度。参见 here对于真正必读的主题。

当顺序循环计算时(假设没有聪明的编译器优化):

gvic = (...((((g_1 + g_2) + g_3) + g_4) + g_5) + ...)

其中 g_i 是通过迭代 i 添加到 gvic 的值，并行版本计算:

gvic = t_0 + t_1 + t_2 + ... t_(#threads-1)

其中t_i 是线程i 中gvic 的累积私有(private)值(OpenMP 中的线程即使在Fortran 中也是0 编号)。未指定不同 t_i 的缩减顺序。 OpenMP 实现可以自由选择它认为合适的任何内容。即使所有的 t_i 都按顺序求和，结果仍然与顺序循环计算的结果不同。不稳定的数值算法在并行化时特别容易产生不同的结果。

这是您很难完全避免的事情，而是学会控制或简单地接受其后果。在许多情况下，问题的数值解无论如何都是一个近似值。您应该关注守恒或统计属性。例如，遍历分子动力学模拟可能会并行产生完全不同的相轨迹，但总能量或热力学平均值等值将非常接近(除非存在一些严重的算法错误或非常糟糕的数值不稳定性)。

旁注 - 您现在进入这个领域实际上很幸运，因为大多数 CPU 使用标准的 32 位和 64 位浮点运算。多年前，当 x87 流行时，浮点运算以 80 位内部精度完成，最终结果将取决于值离开并重新进入 FPU 寄存器的次数。