并行仿真与串行并行运行和附加并行运行相比，经过一定的时间步骤，得到了不同的结果。

Question

我试图运行一个代码的涡旋模拟并行使用OpenMP。这类似于粒子模拟，在每一时间步骤中，涡旋在下一时间步长的位置必须根据它的速度来计算，而速度是由其他旋涡在当前时间步长上的位置决定的。旋涡一旦离开域就会被删除。我比较了并行版本的代码和串行版本的代码在每一时间步骤上的涡数，并多次运行每个版本。

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

我只在并行do结构中使用了减缩条款。整个代码中唯一的并行区域是函数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

如果我构造并行代码的方式是错误的，那么错误应该出现在最初的几个时间步骤本身，对吗？

(在这个https://i.stack.imgur.com/Kg6DQ.png中可以看到，'blobs‘和'sheets’只是涡旋元素的类型，蓝线是总元素。P和S分别代表并行和串行，R代表运行。THe实心绘图标记是串行代码，空心标记是并行代码的三次运行)

编辑:当我将变量的数值精度改为实(4)时，结果的差异发生在比实际(8)更早的时间步骤上。因此，这显然是一个全面的错误问题。

TLDR:我想向那些在一系列时间步骤中看到这样一个结果的其他人澄清这一点，在这些步骤中，并行代码与前几个时间步骤匹配，然后就会发散？

Answer 1

从本质上说，您的代码在gvic中概括了许多术语。浮点运算是不相联的，即由于四舍五入，(a+b)+c和a+(b+c)不一样.此外，根据数值和符号的条件，可能会有严重的精度损失在每一个操作。有关这一主题的强制性阅读，请参见这里。

当顺序循环计算时(没有给出聪明的编译器优化)：

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的累积私有值( OpenMP中的线程甚至在Fortran中为0编号)。减少不同t_is的顺序未指定。OpenMP实现可以自由选择它认为好的任何东西。即使所有的t_i按顺序求和，结果仍然与顺序循环计算的结果不同。不稳定的数值算法在并行化时特别容易产生不同的结果。

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

顺便提醒一下，当大多数CPU使用标准的32位和64位浮点算法时，您现在很幸运地进入了这个字段。几年前，当x87还是一件事的时候，浮点操作是以80位的内部精度完成的，最终结果将取决于一个值离开和重新进入FPU寄存器的次数。