CC BY
16
УДК 004.021:550.34
ОПТИМИЗАЦИЯ ПАРАЛЛЕЛЬНЫХ SWEEPING МЕТОДОВ ЧИСЛЕННОГО РАСЧЕТА ВРЕМЕН ПРОБЕГА СЕЙСМИЧЕСКИХ ВОЛН ДЛЯ ВЫЧИСЛИТЕЛЬНЫХ СИСТЕМ С ОБЩЕЙ ПАМЯТЬЮ
Александр Алексеевич Никитин
Институт нефтегазовой геологии и геофизики им. А. А. Трофимука СО РАН, 630090, Россия, г. Новосибирск, пр. Академика Коптюга, 3, аспирант, тел. (383)330-60-18, e-mail: NikitinAA@ipgg.sbras.ru
Александр Сергеевич Сердюков
Институт нефтегазовой геологии и геофизики им. А. А. Трофимука СО РАН, 630090, Россия, г. Новосибирск, пр. Академика Коптюга, 3, кандидат физико-математических наук, научный сотрудник, тел. (383)335-64-57, e-mail: AleksanderSerdyukov@yandex.ru
Антон Альбертович Дучков
Институт нефтегазовой геологии и геофизики им. А. А. Трофимука СО РАН, 630090, Россия, г. Новосибирск, пр. Академика Коптюга, 3, кандидат физико-математических наук, зав. лабораторией, тел. (383)363-67-14, e-mail: DuchkovAA@ipgg.sbras.ru
Численное решение уравнения эйконала применяется для расчета времен первых вступлений волн в задачах сейсморазведки. В докладе представлен новый способ параллельной реализации алгоритмов Fast Sweeping и Locking Sweeping решения уравнения эйконала для вычислительных систем с общей памятью, позволяющий ускорить вычисления за счет эффективного использования кэша процессора.
Ключевые слова: уравнение эйконала, времена первых вступлений, fast sweeping method, locking sweeping method.
OPTIMIZATION OF PARALLEL SWEEPING METHODS OF NUMERICAL COMPUTATION OF SEISMIC WAVE TRAVEL TIMES FOR SHARED MEMORY COMPUTING SYSTEMS
Alexandr A. Nikitin
Trofimuk Institute of Petroleum Geology and Geophysics SB RAS, 630090, Russia, Novosibirsk, Koptyug Prospect 3, Ph. D., Student, tel. (383)330-60-18, e-mail: NikitinAA@ipgg.sbras.ru
Alexandr S. Serdyukov
Trofimuk Institute of Petroleum Geology and Geophysics SB RAS, 630090, Russia, Novosibirsk, Koptyug Prospect 3, Ph. D., Researcher, tel. (383)335-64-57, e-mail: AleksanderSerdyukov@yandex.ru
Anton A. Duchkov
Trofimuk Institute of Petroleum Geology and Geophysics SB RAS, 630090, Russia, Novosibirsk, Koptyug Prospect 3, Ph. D., Head of the Laboratory, tel. (383)363-67-14, e-mail: DuchkovAA@ipgg.sbras.ru
Numerical solutions to the eikonal equation are used to compute first-arrival travel times of seismic waves. Here we present а new parallel implementation of the Fast Sweeping and the Locking Sweeping methods to speed-up computations due to efficient cache use.
Key words: eikonal equation, first arrival travel times, fast sweeping method, locking sweeping method.
Распространение сейсмических волн описывается системой уравнений эластодинамики. В лучевом приближении возникают независимые системы уравнений для описания распространения продольных и поперечных волн, а их кинематика (расчет времен пробега) описывается двумя уравнениями эйконала для продольных и поперечных волн соответственно, см. разд. 2.4.2 [1].
Решение прямой кинематической задачи является важным элементом во многих алгоритмах обработки сейсмических данных, таких как сейсмическая томография, миграция и др. Нелинейность уравнения эйконала приводит к возникновению неустойчивостей в ходе расчетов, что требует регуляризации численной схемы, т. е. использования так называемых «вязких» решений [2]. Эти решения позволяют найти только времена первых вступлений волн, что оказывается достаточным для многих приложений. Разработаны численные методы решения уравнения эйконала, например, такие как Fast Sweeping [3] или Fast Marching [4]. В случае больших скоростных моделей возникает необходимость в параллельных алгоритмах для быстрого решения прямой кинематической задачи, что и является целью данной работы.
Уравнение эйконала является нелинейным дифференциальным уравнением, которое имеет несколько форм, см. разд. 3.1.1 [1]. Рассмотрим уравнение эйконала в следующей форме:
где - неизвестная функция, описывающая время пробега волны в точку x, /(х) - заданная положительная функция медленности (величина, обратная к скорости распространения волны) в точке х, О. - расчетная область пространства Яп, Г - подобласть в О (точка или область вокруг сейсмического источника) с границей г?Г с заданным фиксированным значением времени пробега
Fast Sweeping Method (FSM) [3] является итерационным алгоритмом решения уравнения эйконала, подходящим для эффективной параллельной реализации. В начале вычислительного процесса в области источника решение задается аналитическим способом, а в расчетной области полагается равным бесконечности. Для дискретизации уравнения эйконала используется противопоточная конечно-разностная аппроксимация, которая в двумерном случае имеет вид:
| Vi(x)| = /(х), хеО\Гс Rn,
с заданными краевыми условиями:
t(x) = g(x), х е аг,
g(x).
_t. . Vnew-1 ■ ■ )+Il= if. h1
new
hmin = niin(//_|/,// + |/ ), tjmin = rnin(/,. /-b/,. / + l )>
где ^ - шаг сетки дискретизации функций t(x) и f(x). В трехмерном случае про-тивопоточная аппроксимация имеет аналогичный вид. Приведенное выражение
представляет собой квадратное алгебраическое уравнение относительно t^/ -
нового значения времени пробега в соответствующем узле декартовой сетки. Использование данной аппроксимации обеспечивает правильный порядок пересчета решения в рамках конечно-разностного шаблона.
Для поддержания правильного порядка пересчета во всей расчетной области в FSM на каждой итерации производится смена направлений обхода, соответствующих координатным осям с положительными или отрицательными направлениями вдоль каждой из координат. Это позволяет уточнять времена прихода для областей, где направление движения волн меняется. Так, в двумерном случае получаются четыре направления обхода сетки: 1) i=1:I, j=1:J, 2) i=I:1, j=1:J, 3) i=I:1, j=J:1, 4) i=1:I, j=J:1, в трехмерном - аналогично восемь направлений. При обходе в каждом направлении обновляются только те узлы сетки, в которых новое вычисленное значение времени меньше предыдущего.
Locking Sweeping Method (LSM) [5] является модификацией FSM, направленной на устранение необходимости расчета нового значения в узлах сетки, которые заведомо не будут обновлены на текущей итерации. Для этого в данном методе используются так называемые «замки», или флаги состояния. На шаге инициализации закрыты замки всех узлов, кроме источника и его соседних узлов по шаблону. Если замок находится в состоянии «закрыт» для текущего узла, то значение в нем не вычисляется и алгоритм переходит к следующему узлу. Иначе вычисляется новое значение, и, если оно меньше предыдущего, открываются замки у всех соседей по шаблону, у которых текущее время больше нового значения. Затем в текущем узле замок закрывается. В остальном данный метод полностью аналогичен FSM и требует столько же итераций для сходимости.
Узлы сетки в обоих методах можно разбить на множества «независимых» узлов, значения в которых можно вычислять параллельно. Примером такого разбиения являются диагонали в двумерном случае и диагональные плоскости в трехмерном, которые используются в параллельной реализации FSM [6], обозначаемой далее DFSM. Вычисления узлов с диагональных плоскостей (рис. 1, слева) распределяются по потокам исполнения. Между вычислениями различных плоскостей осуществляется глобальная синхронизация потоков. Как показало наше тестирование (см. ниже), данный алгоритм обладает низкой эффективностью по сравнению с последовательным алгоритмом FSM вследствие низкой локальности обращений к памяти внутри потоков, что не позволяет эффективно использовать кэш-память процессора.
Для решения этой проблемы нами были предложены новые оптимизированные под эффективную работу с кэшем параллельные реализации FSM и LSM для вычислительных систем с общей памятью, обозначаемые далее как BFSM и BLSM соответственно. Предлагаемые алгоритмы основаны на разбие-
нии сетки на блоки (рис. 1, в центре), вычисления внутри которых производятся последовательно, что обеспечивает хорошую локальность обращений к памяти. Для повышения эффективности параллельного исполнения вместо глобальной синхронизации задаются индивидуальные зависимости между всеми блоками. Блок может быть взят потоком на вычисление, как только будут вычислены его соседние блоки, предшествующие ему по текущему направлению обхода. Например, при прямом направлении обхода (по направлению координатных осей) блок (i, j, k) может быть посчитан после блоков (i - 1, j, k), (i, j - 1, k) и (i, j, k -1) (рис. 1, справа). В нашей параллельной реализации для обеспечения синхронизации вычисления блоков мы используем конструкцию task depend стандарта OpenMP 4.0, позволяющую указывать непосредственные зависимости между задачами.
3D декомпозиция
Рис. 1. Слева пример диагональной плоскости в DFSM. В центре трехмерная декомпозиция в BFSM/BLSM. Справа пример зависимостей между блоками в BFSM/BLSM для заданного направления обхода
На рис. 2. представлены результаты тестирования производительности рассмотренных реализаций FSM и LSM на кластере НГУ, на узлах HP BL2x220c G7 с двумя 6-ядерными процессорами Intel Xeon X5670. Результаты приведены для однородной скоростной модели на сетке 600x600x600 с источником в центре области.
Рис. 2. Результаты тестирования производительности. Пунктирной и штрихпунктирной линиями показано время работы последовательных реализаций FSM и LSM
Как видно из результатов, предлагаемые параллельные реализации позволяют достичь высокой эффективности по сравнению с последовательными методами и существующей параллельной реализацией FSM [6]. На 12 ядрах эффективность BFSM составила 90 %, а BLSM - 76 %. Более низкая эффективность BLSM может быть вызвана неравномерностью объема вычислений в блоках в алгоритме LSM. Тем не менее эта реализация оказалась самой быстрой из рассмотренных. В обоих реализациях при тестировании был выбран размер блока 1х10х600.
В дальнейшем планируется использование разработанных алгоритмов для вычисления времен первых вступлений при решении томографических задач в разведочных и сейсмологических приложениях.
Исследование выполнено при частичной финансовой поддержке стипендии Президента РФ для молодых ученых и аспирантов № СП-2899.2015.5 и гранта РФФИ № 15-05-06752 А.
БИБЛИОГРАФИЧЕСКИЙ СПИСОК
1. Cerveny V. Seismic ray theory. - Cambridge University Press, 2001. - 713 p.
2. Crandall M.G., Lions P.L. Viscosity solutions of Hamilton-Jacobi equations // Transactions of the American Mathematical Society. - 1983. - Vol. 277. - N 1. - P. 1-42.
3. Zhao H. A fast sweeping method for eikonal equations // Mathematics of computation. -2005. - Vol. 74. - N 250. - P. 603-627.
4. Sethian J.A. Level set methods and fast marching methods. - Cambridge university press, 1999. - 404 p.
5. Bak S., McLaughlin J., Renzi D. Some improvements for the fast sweeping method // SIAM Journal on Scientific Computing. - 2010. - Vol. 32. - N 5. - P. 2853-2874.
6. Detrixhe M., Gibou F., Min C. A parallel fast sweeping method for the Eikonal equation // Journal of Computational Physics. - 2013. -Vol. 237. - P. 46-55.
© А. А. Никитин, А. С. Сердюков, А. А. Дучков, 2016
