CC BY
110
УДК 004.021
АЛГОРИТМ РАСЧЕТА ВРЕМЕН ПЕРВЫХ ВСТУПЛЕНИЙ ДЛЯ РЕШЕНИЯ ТРЕХМЕРНЫХ ТОМОГРАФИЧЕСКИХ ЗАДАЧ С ИСПОЛЬЗОВАНИЕМ ВЫЧИСЛИТЕЛЬНЫХ ПЛАТФОРМ С РАЗДЕЛЬНОЙ ПАМЯТЬЮ*
Александр Алексеевич Никитин
Институт нефтегазовой геологии и геофизики им. А. А. Трофимука СО РАН, 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 method, сингулярное разложение, параллельный алгоритм.
PARALLEL IMPLEMENTATION OF SEISMIC TRAVELTIME COMPUTATION FOR 3D TOMOGRAPHY ON MPI PLATFORMS
Alexandr A. Nikitin
Trofimuk Institute of Petroleum Geology and Geophysics SB RAS, 630090, Russia, Novosibirsk, 3 Koptyug Prospect, 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, 3 Koptyug Prospect, 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, 3 Koptyug Prospect, Ph. D., Head of the Laboratory, tel. (383)363-67-14, e-mail: Duchko-vAA@ipgg.sbras.ru
Работа выполнена при частичной поддержке гранта РФФИ (номер проекта 15-05-06752) и грантов компании ВР для молодых ученых.
We present a new parallel implementation of fast-sweeping method of solving eikonal equation for computing seismic traveltimes. The proposed algorithm is the base of developing 3D tomography package for MPI platforms.
Key words: seismic tomography, eikonal equation, fast sweeping method, singular value decomposition, parallel algorithm.
Метод сейсмической томографии получил широкое распространение при обработке региональных и глобальных сейсмологических данных. Повышение эффективности вычислительных алгоритмов является актуальной задачей. Частично это связано с появлением новых алгоритмов, обладающих большей пространственной разрешенностью. Рассматриваемые подходы требуют использования мелких шагов по пространству при аппроксимации исследуемых трехмерных скоростных разрезов. В результате возникает необходимость хранения масштабных трехмерных сеток в оперативной памяти. Другой важной задачей является ускорение вычислений. В частности, в ряде регионов России происходит развертывание региональных систем наблюдения. Резко растет число фиксируемых сейсмических событий и объем данных для обработки. Таким образом, актуальным направлением исследований является поиск походов использования возможностей высокопроизводительных вычислительных платформ с разделенной памятью.
Ключевым этапом разработки параллельных томографических алгоритмов является решение прямой кинематической задачи расчета времен пробега волн и построения томографической матрицы, т. к., собственно, решение системы линейных уравнений является достаточно универсальной процедурой и для ее параллельной реализации существуют стандартные библиотеки. Наиболее универсальным подходом к расчету времен пробега волн является конечно-разностное решение уравнения эйконала и процедура обратного лучевого трассирования на его основе. В силу необходимости учета графа зависимостей значений решения уравнения в узлах расчетной сетки, эффективное распараллеливание вычислений является отдельной научной задачей. Использование разделенной памяти усложняет эту задачу.
В работе предложен алгоритм параллельной реализации конечно-разностного решения уравнения эйконала методом fast-sweeping [1], подходящий для вычислительных систем с распределенной памятью. Результат является новым, т.к. задача реализации рассматриваемого метода на платформах данного типа ранее не рассматривалась. С нашей точки зрения, метод fast-sweeping является в настоящее время единственной схемой, допускающей относительно простое и эффективное распараллеливание на платформах с распределенной памятью. Предложенный параллельный алгоритм расчета пробега был использован авторами при разработке реализации метода сейсмической томографии для решения масштабных трехмерных задач.
Параллельный алгоритм расчета времен пробега
Основу используемой схемы сейсмической томографии составляет расчет времен пробега сейсмических волн с помощью численного решения уравнения эйконала:
|Vu(x)| = f(x), x е Rn (1.1)
В качестве базового численного алгоритма решения уравнения эйконала на регулярных сетках для параллельной реализации в данной работе был выбран Fast Sweeping Method (FSM) [1] - итеративный метод, использующий схемы первого порядка с разностями против потока для дискретизации частных производных и итерации Гаусса-Зейделя с чередующимся направлением обхода. Алгоритм FSM в оригинальной формулировке является последовательным, т.к. каждая итерация цикла обхода сетки с заданным направлением зависит от результата предыдущих итераций данного цикла. Решение в узле с индексами (i, j, k) для заданного направления обхода может быть обновлено после того, как будут обновлены элементы с индексами, меньшими или большими (i, j, k), в зависимости от направления обхода. Таким образом, решение может одновременно обновляться в узлах сетки из множества индексов (i, j, k), задаваемого числом C = dii + djj + dkk, где di,dj,dk е {-1,1} определяются направлением обхода.
В данной работе была разработана следующая новая параллельная реализация алгоритма FSM для работы в системах с распределенной памятью. Для распределения задач по процессам используется декомпозиция трехмерной вычислительной сетки по двум измерениям между процессами и локальная декомпозиция по третьему измерению внутри процессов на подзадачи (рис. 1). При таком распределении данных образуется виртуальная топология «решетка», в которой каждый процесс отвечает за вычисление фиксированного блока сетки {ui,j,k |i e[i1,i2],je[j1,j2],k e[0,K -1]}.
Рис. 1. Схема реализации алгоритма FSM в распределенной памяти:
а) топология коммуникаций между процессами. Указаны ]) координаты процессов в топологии и направление обхода; б) в ячейках указан порядок вычисления блоков подзадач в (^ ]) процессах для заданного направления обхода и направление обмена теневыми гранями
Исходя из зависимостей по данным в алгоритме FSM, между соседними по топологии процессами производится обмен теневыми гранями на границах блоков: в процессе текущей итерации Гаусса-Зейделя по направлению обхода, а до начала следующей итерации завершается обмен гранями против направления текущего обхода. Двумерная декомпозиция между процессами позволяет запускать вычисления в двух соседних процессах по топологии одновременно, а разбиение блока по третьему измерению внутри процесса на подзадачи позволяет избежать пересылок больших объемов данных между процессами и быстрее запускать вычисления в соседних процессах, чем это было бы возможно при использовании только двумерной декомпозиции сетки. Запуск вычисления подзадачи в процессе производится, как только в его память поступят все необходимые данные, то есть будут обновлены теневые грани.
Тестирование алгоритма: эффективность и пример расчета
В ходе тестирования программной реализации приведенного подхода на кластере Новосибирского государственного университета нам удалось добиться эффективности распараллеливания около 70 %. Расчеты проводились на скоростной модели объемом 22 Гб: сетка с размерами 1800 точек в каждом направлении, было задействовано 36 узлов кластера, в каждом по 4 потока. Время расчета составило около 2 минут. Данный результат не является окончательным. Ведутся поиски алгоритмических и технических решений для дальнейшей оптимизации программ. Пример вычислений времен пробега при помощи рассматриваемого алгоритма показан на рис 2. Представлены контуры времен пробега в трехмерной скоростной модели с высококонтрастным включением.
А2В2
у (км)
Рис. 2. Пример расчета времен пробега для сложной трехмерной модели:
а) трехмерная скоростная модель и плоскости сечений; б) двумерные сечения модели по приведенным плоскостям и изолинии времен первых вступлений, рассчитанные предложенным в работе алгоритмом
Томография на платформах с разделенной памятью
В данном разделе мы кратко опишем нашу реализацию метода сейсмической томографии. В настоящее время при численном решении задач сейсмической томографии преобладает матричный подход [3]. При таком подходе уточнение скоростной модели сводится к решению системы линейных алгебраических уравнений (СЛАУ) большого размера:
Ax = b, (1.2)
где b - вектор невязок времен пробега для наблюдаемых и синтетических данных, x - неизвестное приращение медленности, A - томографическая матрица. Эффективные методы построения матрицы (1.2) основаны на использования вычисленных полей времен пробега. В случае классической лучевой сейсмической томографии может быть использована простая и эффективная процедура обратного лучевого трассирования. Она сводится к построению лучей, соединяющих источники с приемниками, путем построения линий тока векторного поля градиента вычисленных времен пробега. Данная процедура была реализована нами для систем с распределенной памятью на основе приведенного выше алгоритма расчета времен пробега сейсмических волн. Томографическая матрица A является разреженной и плохо обусловленной, поэтому для решения системы уравнений требуется регуляризация. Нами был использован метод усеченного сингулярного разложения [4, 5].
Поиск старшей части сингулярного спектра в рассматриваемой ситуации целесообразно проводить при помощи итерационных процедур вычисления
собственных чисел и векторов разреженной матрицы AA типа Арнольди [7]. Для эффективной реализации структур хранения разреженных матриц нами была использована библиотека PETSc [6], а для реализации алгоритма усеченного разложения - библиотека SLEPc [7].
В работе был предложен новый параллельный алгоритм расчета времен первых вступлений сейсмических волн, для которого показана эффективность распараллеливания 70 % при выполнении 36 узлах. Алгоритм является составной частью параллельной реализации метода трехмерной сейсмической томографии, которую планируется применять при обработке сейсмологических данных.
БИБЛИОГРАФИЧЕСКИЙ СПИСОК
1. Zhao H. A fast sweeping method for eikonal equations // Mathematics of computation. -2005. - Т. 74, № 250. - С. 603-627.
2. Detrixhe M., Gibou F., Min Ch. A parallel fast sweeping method for the Eikonal equation // Journal of Computational Physics. - 2013. - Т. 237. - С. 46-55.
3. Natterer F. The mathematics of computerized tomography. Siam. - 1986. - V. 32.
4. Решения уравнений первого рода с компактным оператором в гильбертовых пространствах: существование и устойчивость // Доклады Академии наук. - 1997. - Т. 335, № 3. -С. 308.
5. Кабанник А. В., Орлов Ю. А., Чеверда В. А. Численное решение задачи линейной сейсмической томографии на проходящих волнах: случай неполных данных // Сиб. журн. индустр. матем. - 2004. - Т. 7 (2). - C. 54-67.
6. Balay S. et al. PETSc Users Manual // Argonne National Laboratory, Tech. Rep. ANL-95/11. - Revision 3.5. - 2014.
7. Roman J. E., Campos C., Romero E., Tomas A. SLEPc Users Manual // D. Sistemes Informatics i Computació, Universität Politécnica de Valencia, Tech. Rep. DSIC-II/24/02. - Revision 3.5. - 2014.
© А. А. Никитин, А. С. Сердюков, А. А. Дучков, 2015
