Параллельное моделирование дает разные результаты после некоторых временных шагов по сравнению с последовательным и дополнительным параллельным запуском.

Question

Я пытаюсь запустить код для моделирования вихрей параллельно с использованием OpenMP. Это похоже на моделирование частиц, где на каждом временном шаге положение вихря на следующем временном шаге должно вычисляться на основе его скорости, которая определяется положениями всех других вихрей на текущем временном шаге. Вихри удаляются, когда они покидают область. Я сравниваю количество вихрей на каждом временном шаге для параллельной версии кода с последовательной версией кода и запускаю каждую версию несколько раз.

Для последовательных версий счетчики вихрей точно совпадают на каждом временном шаге. Для параллельного случая все прогоны совпадают с последовательным случаем для нескольких десятков временных шагов, после чего каждый параллельный прогон показывает разницу, но остается в пределах ошибки 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

Если способ, которым я построил параллельный код, неверен, то ошибки должны появиться с первых нескольких сами временные шаги, верно?

(Как видно из этого результата , «капли» и «листы» - это просто типы вихревых элементов, синяя линия - это итоговые элементы. P и S обозначают параллельный и последовательный соответственно, а R обозначает прогоны. solid маркеры на графике представляют собой серийный код, а полые - три прогона параллельного кода)

РЕДАКТИРОВАТЬ: когда я вместо этого меняю числовую точность своих переменных на реальную (4), дивергенция c в результатах происходит на более раннем временном шаге, чем реальный (8) случай выше. Итак, это явно проблема с ошибкой округления.

TL; DR: Я хочу прояснить это с любым, кто видел такой результат в диапазоне временных шагов, где параллельный код совпадает для первых нескольких временных шагов и потом расходится?

Answer 1

Ваш код, по сути, суммирует множество терминов в gvic. Арифметика с плавающей запятой c не ассоциативна, то есть (a+b)+c не совпадает с a+(b+c) из-за округления. Кроме того, в зависимости от значений и знаков в терминах может быть серьезная потеря точности в каждой операции. См. здесь для действительно обязательного чтения по теме.

В то время как последовательный l oop вычисляет (без умной оптимизации компилятора):

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

, где g_i - значение, добавленное к gvic с помощью итерации i, параллельная версия вычисляет:

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

, где t_i - накопленное частное значение gvic в thread i (потоки в OpenMP нумеруются даже в Fortran). Порядок уменьшения различных t_i s не указан. Реализация OpenMP может свободно выбирать то, что считает приемлемым. Даже если все t_i суммировать по порядку, результат все равно будет отличаться от результата, вычисленного последовательным l oop. Нестабильные численные алгоритмы исключительно склонны давать разные результаты при распараллеливании.

Это то, чего вы вряд ли сможете полностью избежать, вместо этого научитесь контролировать или просто жить с его последствиями. Во многих случаях численное решение проблемы в любом случае является приближением. Вам следует сосредоточиться на сохраняемых или статистических свойствах. Например, моделирование молекулярной динамики ergodi c может параллельно генерировать совершенно другую фазовую траекторию, но такие значения, как полная энергия или средние значения термодинамики c, будут довольно близкими (если только не будет какой-либо серьезной ошибки алгоритма c или очень плохая числовая нестабильность).

Примечание - вам действительно повезло ввести это поле сейчас, когда большинство процессоров используют стандартные 32- и 64-битные арифметические операции с плавающей запятой c. Годы a go, когда был x87, операции с плавающей запятой выполнялись с 80-битной внутренней точностью, и конечный результат зависел от того, сколько раз значение покидает и повторно входит в регистры FPU.