CC BY
17
УДК 550.34.013
DOI: 10.183 03/2618-981X-2018-4-152-161
СРАВНЕНИЕ 3D ВРЕМЕННЫХ И ЧАСТОТНЫХ РЕШАТЕЛЕЙ ДЛЯ ПРОМЫШЛЕННЫХ ЗАДАЧ ПОЛНОГО ОБРАЩЕНИЯ ВОЛНОВЫХ ПОЛЕЙ (FWI)
Сергей Александрович Соловьев
Институт нефтегазовой геологии и геофизики им. А. А. Трофимука СО РАН, 630090, Россия, г. Новосибирск, пр. Академика Коптюга, 3, кандидат физико-математических наук, старший научный сотрудник, тел. (383)330-60-16, e-mail: solovevsa@ipgg.sbras.ru
Виктор Иванович Костин
Институт нефтегазовой геологии и геофизики им. А. А. Трофимука СО РАН, 630090, Россия, г. Новосибирск, пр. Академика Коптюга, 3, кандидат физико-математических наук, старший научный сотрудник, тел. (383)330-60-16, e-mail: KostinVI@ipgg.sbras.ru
Сравнивается производительность двух акустических решателей с использованием реалистичной трехмерной скоростной модели, различного количества кластерных узлов и источников. Для небольшого числа источников эффективнее работает временной решатель (TD), а для большого количества источников - частотной решатель (FD). Аналогично, FD быстрее на малом числе кластерных узлов, а TD - на большом. Таким образом, оба решателя могут быть использованы при разработке инструментов полного обращения волновых полей.
Ключевые слова: полное обращение волновых полей, уравнение Гельмгольца, система линейных алгебраических уравнений (СЛАУ), высокопроизводительные вычисления.
COMPARING 3D TIME- AND FREQUENCY-DOMAIN SOLVERS FOR INDUSTRIAL FULL-WAVEFORM INVERSION (FWI) APPLICATIONS
Sergey A. Solovyev
Trofimuk Institute of Petroleum Geology and Geophysics SB RAS, 3, Prospect Аkademik Koptyug St., Novosibirsk, 630090, Russia, Ph. D., Senior Researcher, phone: (383)330-60-16, e-mail: solovevsa@ipgg.sbras.ru
Victor I. Kostin
Trofimuk Institute of Petroleum Geology and Geophysics SB RAS, 3, Prospect Аkademik Koptyug St., Novosibirsk, 630090, Russia, Ph. D., Senior Researcher, phone: (383)330-60-16, e-mail: KostinVI@ipgg.sbras.ru
We compare the performance of two acoustic solvers by using a realistic 3D velocity model, various number of cluster nodes and shots taken from a 3D seismic acquisition scenario. When the number of shots is small, a time-domain (TD) solver performs best. For large number of shots, the frequency-domain (FD) solver is better. Likewise, FD performs better on smaller number of nodes; TD becomes the better option for a large number of nodes. Therefore, both solvers deserve a place in a modern FWI toolkit.
Key words: full-waveform inversion, Helmholtz equation, System of Linear Algebraic Equations (SLAE), High Performance Computing (HPC).
Введение
Численное моделирование акустических волновых полей является неотъемлемой частью алгоритмов, разработанных для решения задач, возникающих в геофизической разведке. В частности, такие алгоритмы используются в задаче полного обращения волновых полей (от английского Full-Waveform Inversion, FWI), которая обычно выполняется в частотной области [1], [2]. Для макромас-штабного восстановления скоростной модели моделирование выполняется для набора низких частот (до 10 Гц) на каждом шаге итерационного процесса. В таком моделировании поле давления порождается точечным источником, работающим как гармонический осциллятор на заданной частоте. В промышленных масштабах число источников достигает десятки тысяч и более, в то время как приемников обычно еще больше.
Любой подход к моделированию акустического волнового поля в частотной области имеет как преимущества, так и недостатки, которые могут иметь решающее значение в конкретной ситуации. Здесь мы рассмотрим два подхода, сравним их точность и производительность.
Первый основан на моделировании во временной области, с последующим применением преобразования Фурье. Другой подход относится непосредственно к решению уравнения Гельмгольца. В 3D это требует решения системы линейных уравнений с большой разреженной матрицей коэффициентов. Итерационные решатели являются наиболее популярными для решения таких систем
[3]. Сравнительно новые идеи заключаются в сжатии промежуточных данных
[4], [5] в прямых алгоритмах решения систем линейных уравнений с разреженными матрицами коэффициентов. Для сравнения с временным решателем нами и используются эти идеи сжатия.
Во временной области для решения волнового уравнения
д2
- с2 (X,y,z)Ap = f (t)£(x - x,y - ys,z - zs)
dt
используется моделирование поля давления во временной области для исходного источника, определенного в точке (xs, ys, zs), с последующим применением преобразования Фурье
p (ю, x, y, z ) = je~1(0tp (t, x, y, z ) dt
для получения решения в частотной области. Здесь с (x, y, z) - скорость в точке (x, y, z ), f (t) - функция источника, S - дельта функция Дирака, p (t, x, y, z) -исходное давление в точке (x, y, z) во время t. Для численного моделирования
мы используем программный продукт, разработанный консорциумом Seiscope (https://seiscope2.osug.fr). Далее мы будем его называть конечно разностным
временным решателем (TDFD от английского Time-Domain Finite-Difference). При этом для уменьшения влияния границ расчетной области на решения мы используем идеально поглощающий слой (Perfectly Matched Layers, PML) [6]. При этом, после преобразования Фурье получаем одновременно решение для разных частот. Для различных положений источника вычисления выполняются независимо друг от друга. Это приводит к идеальному MPI-параллелизму TDFD решателя относительно числа источников, предполагая, что один кластерный узел решает один источник.
Рассматривая задачу в частотной области, необходимо решить уравнение Гельмгольца:
со2
Au + ~г:-:u = s( х - х, у - ys,z - zs )• (1)
с (х,y, z)
Дополняя это уравнение условиями PML, а также значением нуля на границе, приходим к дифференциально краевой задаче, конечно-разностная аппроксимация которой приводит к системе линейных уравнений
Au = f. (2)
Матрица А этой системы является симметричной, разреженной и комплексной из-за условия PML. Ненулевые элементы матрицы зависят от значения скорости с (х, y, z ) в точках расчетной области и от частоты ю. Вектор правой части f определяется положением источника.
Мультифронтальный прямой подход для решения системы (2) состоит из разложении матрицы с переставленными элементами [7] (P - матрица перестановок) в произведение треугольных L,L и диагональной D матриц
A = P ■ A • Pt = L ■ D • L (3)
с последующим системы:
У = L~lPff; z = D-1 y; u =PL4z.
Факторизация (3) выполняется один раз и может использоваться для всех правых сторон. Это основное преимущество прямого решения краевых задач для уравнения Гельмгольца со многими источниками. Будем называть этот подход конечно-разностным прямым решателем в частотной области (frequency-domain finite-difference, FDFD).
Основная сложность этого подхода - большая заполняемость L-фактора ненулевыми элементами, намного превышающие число ненулевых элементов в исходной матрице A. Для решения этой проблемы мы производим промежуточное сжатие данных [8, 9] в процессе факторизации. Эта техника основана на аппроксимации матрицами малого ранга и использовании иерархического формата данных (Hierarchically Semi Separable, HSS) [4, 5, 10]. Реализацию под распределенные вычислительные системы (кластера) такого алгоритма будем далее называть HSS решатель.
Численные эксперименты
Тесты на производительность и на точность, описанные далее, проведены на суперкомпьютере Shaheen II (https://www.hpc.kaust.edu.sa/content/shaheen-ii) со следующими характеристиками каждого кластерного узла: 2 х Intel® Xeon® CPU E5-2698 v3 @2.3 GHz, 128 GB RAM. Для численного сравнения мы используем скоростную модель Overthrust (OT) размера 9 х 9 х 4.5 km. В этой модели скорость изменяется от 2 300 м/с до 6 000 м/c. Пример сечения модели в плоскости X-Z и Y-Z показан на рис. 1. В табл. 1 приведен список параметров решаемых задач вышеупомянутыми решателями.
X Y
а) б)
Рис. 1. Скоростная модель OT, в (a) сечении X-Z и (b) сечении Y-Z
Таблица 1
Параметры численных экспериментов
V (Hz) h (m) ppw Nx Ny Nz N dt (5)
5 90 5.1 110 110 52 0.6-106 0.00742
7 60 5.4 165 165 78 2.1-106 0.00495
15 30 5.1 330 330 155 17-106 0.00247
Для обеспечения примерно равного количество точек на одну длину волны («ppw» в табл. 1) для частот 5, 7 и 15 Гц, решения вычисляются с использованием разных шагов сетки. Соответствующие количества шагов сетки показаны в столбцах Ых, Ыу и N2 с общими числами точек сетки в столбце N. Столбец dt
содержит значения шага дискретизации времени, используемого для TDFD. Расчеты решателя TDFD выполнялись с использованием временного интервала [0,10 с]. Источники размещались на глубине одного шага сетки в центре области, приемники размещены по всей плоскости Х-У на той же глубине. Реальные части вычисленного волнового поля, взятые на глубинах плоскости приемника, показаны на рис. 2.
х (кт) Х(кт)
Рис. 2. Решение (действительная часть), плоскости поверхностных приемников с использованием решения HSS (слева - для 5 Гц, а справа - для 15 Гц; красным отмечены линии, используемые на рис. 3 и рис. 4)
На рис. 3 и рис. 4 показаны решения TDFD и HSS вдоль выбранных на рис. 2 линий, а также увеличенная разница между двумя решениями. Для решения 5 Гц видна хорошая согласованность между кривыми. При более высокой частоте 15 Гц соглашение несколько хуже (рис. 4), но все же приемлемо. Эта разница может быть вызвана недостаточной дискретизацией и / или различным моделированием высококонтрастных границ в разных решателях.
х 10"'
х(кт)
Рис. 3. Реальная часть решения частотной области для 5 Гц, показанная по красному профилю на рис. 2 (красные и синие линии соответствуют решениям TDFD и HSS соответственно; черная линия представляет собой увеличенную (х5) разницу между ними)
х(кш)
Рис. 4. То же, что и на рис. 3, но для решения 15 Гц
Для сравнения решения полученными HSS и TDFD подходами используем различные метрики: три из них - это относительные разницы двух сеточных решений и и V, рассчитанные с использованием классических норм, а четвертая - на основе корреляционной функции:
Рк (и, ^)
и
Ч
ик
к = 1,2, ^,у(и, V)
(и, v)
т
(4)
Функции Рк (г) и у( г) получены вычислением норм в формуле (4) в шаре
радиуса г с центром в точке источника. Волновые поля ^ V в формуле (4), имеют особенности в местоположениях источника. Поэтому при вычислении мы исключаем шар малого радиуса Г) вокруг источника. Графики функций Рк (г)
и у(г) показаны на рис. 5; Рк обозначает Рк (рнбб,Ртт),к = 1,2,^ и у -у(Рнбб,Ртэеэ). Решения Рнбб и Ртбеб относятся к решателям HSS и TDFD соответственно. Радиус Г) в экспериментах равен 10^. Очевидно, что расхождение между двумя решателями возрастает с увеличением радиуса и частотой. Немного больший скачок между 5 и 15 Гц может предполагать недостаточный размер сетки для модели ОТ с высокими контрастами.
Рис. 5. Разница между решениями, полученными с помощью TDFD и HSS, измеренная как функция от размера области
Сравнение численной производительности
Для сравнения относительной производительности решателей TDFD и HSS используются одни и те же вычислительные ресурсы (число узлов) и число источников. Для конкретной пары значений (), время решения (с)
приведено в табл. 2 в формате дроби, в которой в числителе задано время а в знаменателе - .
Таблица 2
Расчетное время двух решателей, представленных как отношения с
в числителе и в знаменателе
О) ^sllots
tHSS.Cs) 1 128 1 280 12 800
32 330 1,320 13,200 132,000
2.911 2,944 3,237 6,170
64 330 660 6,600 66,000
Nnod.es 2,606 2,629 2,840 4,946
128 330 330 3,300 33,000
2,538 2,558 2,739 4,543
256 330 330 1,650 16,500
2,538 2,548 2,638 3,541
Решатель TDFD идеально масштабируется относительно числа источников N н 0 Время решения для одной правой части на одном узле равно 330 с. Соответственно, 128 источников могут быть посчитаны на 32 узлах за 1320 с, на 128 узлах за 330 с., однако дальнейшее увеличение Nn0ае5 не приводит к уменьшению времени счета. Масштабируемость HSS решателя для одного источника определяется масштабируемостью факторизации и ухудшается с увеличением числа Nn0ае5 . Однако увеличение Nsh0дает дополнительную возможность для улучшения масштабируемости.
Рис. 6 показывает, как решатели FDFD и TDFD дополняют друг друга при сравнении N sh 0^ и Nn0йеs. Синяя линия на рис. 6 определяет «линию равной производительности» двух решателей. Для заданного количества узлов 0 а е!;) эта строка определяет количество источников, которое «достаточно большое», чтобы полностью воспользоваться преимуществами решения HSS и достичь численной производительности TDFD. Например, для задачи со сравнительно небольшим количеством источников, которые должны быть решены на конкретном кластере (точка, близкая к точке С сегмента CD), самым быстрым способом будет использование решателя TDFD. При увеличении количества источников (точка перемещается в D вдоль CD) прямой решатель FDFD становится быстрее. Другой способ взглянуть на это - зафиксировать количество источников, например, вдоль горизонтальной линии АВ. Затем двигаться от точки А к В вдоль горизонтального сегмента АВ. Это означает решение той же задачи
на кластерах при увеличении числа узлов. Для фиксированного количества источников прямой решатель FDFD эффективен на относительно небольших кластерах, но для более крупных кластерах TDFD более эффективен.
С ^nodes
Рис. 6. Относительная производительность решателей TDFD и FDFD, показанных относительно количеством источников (вертикальная ось) и относительно количеству узлов (горизонтальная ось). Синяя линия определяет «линию равной производительности»
Заключение
Проведено сравнение двух акустических решателей в рамках задачи полного обращения волнового поля для разных размеров вычислительных кластеров и разного числа источников. Конечно-разностный временной решатель (TDFD) идеально (линейно) масштабируется при увеличении числа кластерных узлов. С другой стороны, производительность решателя HSS нелинейно зависит от доступных ресурсов, а соотношение выгодности в затратах резко возрастает с увеличением количества источников. Серию численных экспериментов, использующих модель Overtrust и суперкомпьютер Shaheen II, показали существование «линии равной производительности» при сравнении числа источников и доступных кластерных узлов. Решатель TDFD обычно выигрывает в сегменте
сравнительно больших значений отношения Nnodes/ , а HSS выигрывает,
/ N shots
когда отношение меньше.
С годами все большее количество узлов становится доступным (вычислительная мощность становится дешевле), однако этому противодействует постоянное увеличение плотности трасс в сейсмических системах сбора данных, которое приводит к постоянному росту числа источников и приемников на квадратный километр поверхности. Эти две параллельные тенденции, вероятно, будут сохранять потребность обоих решателей в производстве без явного преимущества. В зависимости от доступной вычислительной мощности и количества источников для конкретной FWI задачи один или другой решатель может быть значительно более эффективным. Поэтому оптимальный набор инстру-
ментов FWI должен содержать оба решателя, и выбор наиболее эффективного может быть основан на «линии равной производительности», приведенной в данном исследовании.
Работа выполнена при финансовой поддержке Российского научного фонда (РНФ) грант № 17-17-01128.
БИБЛИОГРАФИЧЕСКИЙ СПИСОК
1. Shin C., Cha Y. H., Waveform inversion in the Laplace domain // Geophysical Journal International. - 2008. - Т. 173. - № 3. - С. 922-931.
2. Virieux J., Operto S., Ben-Hadj-Ali H., Brossier R., Etienne V., Sourber F., Giraud L., Haidar A. Seismic wave modeling for seismic imaging // The Leading Edge. - 2009. - Т. 28. -№ 5. - С. 538-544.
3. Belonosov M., Dmitriev M., Kostin V., Neklyudov D., Tcheverda V. An Iterative Solver for the 3D Helmholtz Equation // Journal of Computational Physics. - 2017. - Vol. 345. - P. 330344. doi: 10.1016/j.jcp.2017.05.026
4. Hackbusch W. A sparse matrix arithmetic based on K-Matrices. Part I: Introduction to K-matrices // Computing. - 1999. - Т. 62. - № 2. - С. 89-108. doi: 10.1007/s006070050015
5. Xia, J. Efficient structured multifrontal factorization for large sparse matrices // SIAM Journal of Scientific Computing, - 2013. - Т. 35. - № 2. - С. A832-A860. doi: 10.1137/1208670.
6. Collino F., Tsogka C. Application of the PML absorbing layer model to the linear elastodynamic problem in anisotropic heterogeneous media: Rapport de recherché. - 1998 - 3471, INRIA, Rocquencourt.
7. Duff I. S., Reid J. K. The multifrontal solution of indefinite sparse symmetric linear systems // ACM Transactions on Mathematical Software (TOMS) - 1983. - Т. 9. - № 3. - С. 302325. doi: 10.1145/356044.356047.
8. Solovyev S., Vishnevsky D., Liu H. 2015, Multifrontal hierarchically semi-separable solver for 3D Helmholtz problem using optimal 27-point finite-difference scheme // 77th EAGE Conference and Exhibition 2015.
9. Kostin V., Solovyev S., Liu H., Bakulin A. HSS cluster-based direct solver for acoustic wave equation // 87th Annual International Meeting, SEG Technical Program Expanded Abstracts. - 2017. P. 4017-4021.
10. Wang S., de Hoop M., Xia J., On 3D modeling of seismic wave propagation via a structured parallel multifrontal direct Helmholtz solver // Geophysical Prospecting. - 2011 -Vol. 59. - P. 857-873, doi: 10.1111/j.1365-2478.2011.00982.x.
REFERENCES
1. Shin C., Cha Y. H., Waveform inversion in the Laplace domain // Geophysical Journal International. - 2008. - T. 173. - № 3. - S. 922-931.
2. Virieux J., Operto S., Ben-Hadj-Ali H., Brossier R., Etienne V., Sourber F., Giraud L., Haidar A. Seismic wave modeling for seismic imaging // The Leading Edge. - 2009. - T. 28. -№ 5. - S. 538-544.
3. Belonosov M., Dmitriev M., Kostin V., Neklyudov D., Tcheverda V. An Iterative Solver for the 3D Helmholtz Equation // Journal of Computational Physics. - 2017. - Vol. 345. - P. 330344. doi: 10.1016/j.jcp.2017.05.026
4. Hackbusch W. A sparse matrix arithmetic based on H-Matrices. Part I: Introduction to H-matrices // Computing. - 1999. - T. 62. - № 2. - S. 89-108. doi: 10.1007/s006070050015
5. Xia, J. Efficient structured multifrontal factorization for large sparse matrices // SIAM Journal of Scientific Computing, - 2013. - T. 35. - № 2. - S. A832-A860. doi: 10.1137/1208670.
6. Collino F., Tsogka C. Application of the PML absorbing layer model to the linear elastodynamic problem in anisotropic heterogeneous media: Rapport de recherché. - 1998 - 3471, INRIA, Rocquencourt.
7. Duff I. S., Reid J. K. The multifrontal solution of indefinite sparse symmetric linear systems // ACM Transactions on Mathematical Software (TOMS) - 1983. - T. 9. - № 3. - S. 302325. doi: 10.1145/356044.356047.
8. Solovyev S., Vishnevsky D., Liu H. 2015, Multifrontal hierarchically semi-separable solver for 3D Helmholtz problem using optimal 27-point finite-difference scheme // 77th EAGE Conference and Exhibition 2015.
9. Kostin V., Solovyev S., Liu H., Bakulin A. HSS cluster-based direct solver for acoustic wave equation // 87th Annual International Meeting, SEG Technical Program Expanded Abstracts. - 2017. P. 4017-4021.
10. Wang S., de Hoop M., Xia J., On 3D modeling of seismic wave propagation via a structured parallel multifrontal direct Helmholtz solver // Geophysical Prospecting. - 2011 -Vol. 59. - P. 857-873, doi: 10.1111/j.1365-2478.2011.00982.x.
© С. А. Соловьев, В. И. Костин, 2018
